Example 4: Alternative Evolutionary Models

Eric W. Goolsby edited this page May 11, 2016 · 4 revisions

Alternative evolutionary models: Ornstein-Uhlenbeck, Early-Burst

Again using our simulated phylogeny and dataset with missing observations from the example in Example 1: Getting Started), first we fit a model of evolution using the `phylopars` function.

``````library(Rphylopars)
library(phytools) # for simulating pure-birth phylogenies
set.seed(21) # Set the seed for reproducible results
trait_cov <- matrix(c(4,2,2.2,2,3,1.5,2.2,1.5,2.5),nrow = 3,ncol = 3)
tree <- pbtree(n = 20)
p_BM <- phylopars(trait_data = sim_data\$trait_data,tree = sim_data\$tree)
``````

Next we fit two alternative (non-Brownian) evolutionary models: Ornstein-Uhlenbeck (Hansen 1997) and Early-Burst (Harmon et al. 2010). As with Pagel's lambda, we fit a single tree transformation parameter jointly to all traits, rather than fitting trait-specific parameters. This means that the same alpha (for OU) or rate (for EB) is shared across all traits. For the OU model, this is like a multivariate OU where the alpha matrix is diagonal and with identical values across the diagonal.

``````p_OU <- phylopars(trait_data = sim_data\$trait_data,tree = sim_data\$tree,model = "OU")
p_OU # Estimated trait covariance, OU alpha, and stationary covariance
``````

The phylogenetic trait variance-covariance is the \$\boldsymbol{\Sigma}\$ matrix of the OU process, and the stationary covariance is that same matrix divided by `2*alpha`.

``````p_EB <- phylopars(trait_data = sim_data\$trait_data,tree = sim_data\$tree,model = "EB")
p_EB # Estimated trait covariance and EB rate parameter
``````

Now we may perform model selection via information-based metrics such as AIC and BIC.

``````AIC(p_BM)
AIC(p_OU)
AIC(p_EB)

BIC(p_BM)
BIC(p_OU)
BIC(p_EB)
``````

Because a lower score indicates a better fit, Brownian motion is the best-supported model by both AIC and BIC.

Now we simulate a new dataset under an OU model and perform model selection among BM, EB, and OU models.

``````set.seed(2) # Set the seed for reproducible results
trait_cov <- matrix(c(4,2,2.2,2,3,1.5,2.2,1.5,2.5),nrow = 3,ncol = 3)
tree <- pbtree(n = 30)
sim_data_OU <- simtraits(v = trait_cov,tree = tree,
nmissing=10,model = "OU",parameters = list(alpha = 1.5))

p_BM <- phylopars(trait_data = sim_data_OU\$trait_data,tree = sim_data_OU\$tree,model = "BM")
p_EB <- phylopars(trait_data = sim_data_OU\$trait_data,tree = sim_data_OU\$tree,model = "EB")
p_OU <- phylopars(trait_data = sim_data_OU\$trait_data,tree = sim_data_OU\$tree,model = "OU")

BIC(p_BM)
BIC(p_OU)
BIC(p_EB)
``````

Here, the OU model has a lower score, so OU is selected as the best-supported evolutionary model.

Multivariate Ornstein-Uhlenbeck

We may also fit the multivariate Ornstein-Uhlenbeck model (`model="mvOU"`), which fits an alpha matrix (i.e., the adaptation rate parameter) to the data. (In the previous section, this matrix was constrained to be proportional to the identity matrix.) The alpha matrix can be a full matrix, in which case it influences both phylogenetic correlation between species and correlated evolution between traits. Alternatively, the alpha matrix can be assumed to be diagonal (by setting `full_alpha=FALSE` with `model="mvOU"`), which represents adaptation acting on individual traits only, and in which case alpha influences phylogenetic correlation between species but not correlations between traits.

``````p_mvOU <- phylopars(trait_data = sim_data_OU\$trait_data,tree = sim_data_OU\$tree,model = "mvOU")
p_mvOU # Multivariate OU model (with alpha as full matrix)

p_mvOU_diag <- phylopars(trait_data = sim_data_OU\$trait_data,tree = sim_data_OU\$tree,
model = "mvOU",full_alpha = FALSE)
p_mvOU_diag # Multivariate OU model (with diagonal alpha)
``````

Here, the stationary covariance matrix is a multivariate generalization of \$\boldsymbol{\Sigma}\$ divided by `2*alpha` (see Supplement 1 of Clavel et. al 2015).

Again, we may perform model selection using AIC or BIC.

``````BIC(p_BM)
BIC(p_OU)
BIC(p_EB)
BIC(p_mvOU)
BIC(p_mvOU_diag)
``````

Consistent with simulated conditions, the jointly fit OU model `p_OU` is best-supported (with a shared alpha scalar for all traits, that is, a diagonal alpha matrix where all diagonal entries are equal).

To simulate evolution under a multivariate Ornstein-Uhlenbeck model, refer to the `mvSIM` function in the `mvMORPH` package.

References

Clone this wiki locally
You can’t perform that action at this time.
Press h to open a hovercard with more details.