Bayesian Integrative models of Trait Evolution (BITE)

The BITE package implements the JIVE model and other Bayesian models aimed at understanding traits evolution.

The JIVE model implements a joint estimation of the evolution of intra- and inter-specific variance for continuous traits in a phylogenetic framework. JIVE is described in Kostikova et al. (2016) and is being further developed by T. Gaboriau.

Install BITE package:

library(devtools)
install_github("theogab/bite")

BITE requires the following libraries: ape, phytools, coda, sm, vioplot

Load the library and set the working directory, where the output files will be saved:

library(bite)
setwd("my_path")

A JIVE analysis requires: 1) a phylogenetic tree (ultrametric) 2) a table listing all occurrences with species name and individual trait value 3) [optional] a discrete trait that can be on the tree to define evolutionary regimes

Load example files, a phylogeny and a table with multiple measurements for each species:

data(Anolis_tree)
data(Anolis_traits)

Setting up a model

JIVE can use different models for the evolution of species mean and intra-specific trait variance. The data and the models are stored in an object that is then used to run the analysis. For example:

jive_obj <- make_jive(Anolis_tree, Anolis_traits, model.mean = "BM", model.var = "OU")

specifies a Brownian model of evolution (BM) for the mean and an Ornstein-Uhlenbeck model (OU) for the variance. The white noise model (WN) is also available.

You can plot the input data using: plot_jive(jive_obj)

Launch analysis

The command mcmc_bite launches the Bayesian analysis of trait evolution. The output of the analysis is saved to a test file (file name provided with the flag log.file).

mcmc_bite(jive_obj, log.file = "output_mBM_vOU.log", ngen = 10000, sampling.freq = 100)

The parameters ngen and sampling.freq define the number of MCMC iterations and the sampling frequency, respectively. Note that 10,000 iterations are likely to be too few for the algorithm to find a reliable solution.

Plot output

Import the results in R

logfile <- "output_mBM_vOU.log"
res <- read.csv(logfile, header = T, sep = "\t")

Plot the rate of evolution of species means (BM model):

plot_mcmc_bite(res, var = "mean.bm.sig.")

Plot the parameters of evolution of species variance (OU model). Rate parameter:

plot_mcmc_bite(res, var = "var.ou.sig.")

Optimal log variance:

plot_mcmc_bite(res, var = "var.ou.the.")

Stationary variance (sv = alpha / 2*sig):

plot_mcmc_bite(res, var = "var.ou.sv.")

---

Model testing using JIVE

We can test different models of evolution for trait means and variances and compute their statistical support using marginal likelihoods.

Model 1

Define the JIVE model and output file:

jive_obj <- make_jive(Anolis_tree, Anolis_traits, model.mean = "BM", model.var = "OU")
logfile <- "output_mBM_vOU_TI.log"

When running the analysis the command ncat defines the number of discrete categories for integrating the marginal likelihood:

mcmc_bite(jive_obj, log.file = logfile, sampling.freq = 1000, ngen = 1e+05, 
    ncat = 10, burnin = 1000)

The flag burnin = 100 specifies that the first 1000 iterations should be discarded as burnin.

Read output and compute marginal likelihood:

res <- read.csv(logfile, header = T, sep = "\t")
mlik_mBM_vOU <- marginal_lik(res)
mlik_mBM_vOU

Model 2

Define a new model with BM evolution for both mean and variance:

jive_obj <- make_jive(Anolis_tree, Anolis_traits, model.mean = "BM", model.var = "BM")
logfile <- "output_mBM_vBM_TI.log"
mcmc_bite(jive_obj, log.file = logfile, sampling.freq = 1000, print.freq = 1000, 
    ngen = 1e+05, ncat = 10)

Load results and compute marginal likelihood:

res <- read.csv(logfile, header = T, sep = "\t")
mlik_mBM_vBM <- marginal_lik(res, burnin = 10)
mlik_mBM_vBM

Which model is best supported?

---

Trait-dependent JIVE models

The JIVE model can be adapted to infer different model parameters in different parts of the tree following the evolution of a discrete trait. For example, Kostikova et al. used such a model to test whether the variance in ecological niche varied between (2016) annual and perennial plant lineages.

Load the Anolis data and read the table with a discrete trait assigned to each species, and turn it into a named vector

data(Anolis_tree)
data(Anolis_traits)
tr <- read.table("Anolis_discrete_trait.txt", row.names = 1, h = T)
trait <- unlist(tr)
names(trait) <- row.names(tr)

Load the library phytools and generate a stochastic map of the trait

library(phytools)
mapped_tree = make.simmap(Anolis_tree, trait)
plotSimmap(mapped_tree)

Define a JIVE model using BM for trait mean with independent rates of evolution for each state of the discrete trait. For trait variance we use an OU model with independent optima for each state of the discrete trait.

my.jive <- make_jive(mapped_tree, Anolis_traits, model.mean = c("BM", "sigma"), 
    model.var = c("OU", "theta"))
mcmc_bite(my.jive, log.file = "Anolis_trait_JIVE.log", sampling.freq = 100, 
    ngen = 10000)

Are OU optima significantly different?

post <- read.csv("Anolis_trait_JIVE.log", header = T, sep = "\t")
delta_optima = post$var.ou.the.0 - post$var.ou.the.1
hist(delta_optima)

Calculate the posterior probability of optimum.0 > optimum.1

length(delta_optima[delta_optima > 0])/length(delta_optima)