New version of evolvcv.lite for multistate mapped trait

The existing phytools function evolvcv.lite fits the model of evol.vcv and Revell & Collar (2009), but with various constraints on the evolutionary rate matrices. Specifically, for two quantitative traits and a binary mapped character, evolvcv.lite fits the following four models: the same set of rates and single correlation for both states of the binary character (model 1 - this is the one rate matrix model in evol.vcv); different rates but the same correlation (model 2); different correlations but the same rates (model 3); and, finally, different correlations and rates (model 4 - this is the two rate matrix model from evol.vcv). The reason that I have not so far implemented these partially constrained models for more than two quantitative traits is because doing so constitutes a very challenging problem computationally. However there is no reason why the function could not be expanded to a mapped multistate character (rather than simply a binary character, as in the present implementation). In fact, a phytools user and blog reader recently requested this addition and I have finally gotten around to implementing and debugging this.

This version also fixes a bug in the previous edition which miscalculated the number of parameters (by -1 in all cases), and thus the AIC scores as well.

The new version of this function is here. I have also built a new version of phytools (here) which can be downloaded and installed from source.

One neat trick that I learned in the process of programming this update was how to add corresponding elements of list of matrices. Say, for instance, that X is a list of matrices of identical dimensions. It turns out that sumX<-Reduce("+",X) gives the sum of the matrices in X. Who knew!

Let's try out the new version of evolvcv.lite with an example.

> # install new version of phytools
> install.packages("phytools_0.1-82.tar.gz",type="source", repos=NULL)
Installing package(s) into ‘C:/R/win-library/2.14’
...
* DONE (phytools)
> require(phytools)
Loading required package: phytools
...
> # simulate tree & discrete character history
> tree<-sim.history(pbtree(n=100,scale=1), Q=matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3))
> # simulate uncorrelated traits with three different rates
> # depending on the state
> X<-cbind(sim.rates(tree,c(1,2,3)),sim.rates(tree,c(1,2,3)))
> # fit models
> R<-evolvcv.lite(tree,X)
> # in the following I have hidden some results for brevity
> # only the full output for the best fitting model is shown
> R
$model1
$model1$description
[1] "common rates, common correlation"
...
$model1$AIC
[1] 379.6677

$model2
$model2$description
[1] "different rates, common correlation"
$model2$R
$model2$R$`1`
         [,1]     [,2]
[1,] 1.218973 0.096771
[2,] 0.096771 1.244826
$model2$R$`2`
         [,1]     [,2]
[1,] 1.6109460 0.1451599
[2,] 0.1451599 2.1194570
$model2$R$`3`
         [,1]     [,2]
[1,] 1.6255930 0.2015308
[2,] 0.2015308 4.0483984
$model2$logLik
[1] -181.8455
$model2$convergence
[1] 0
$model2$k
[1] 9
$model2$AIC
[1] 381.6911

$model3
$model3$description
[1] "common rates, different correlation"
...
$model3$AIC
[1] 382.892

$model4
$model4$description
[1] "no common structure"
...
$model4$AIC
[1] 385.1049

Cool. The best fitting model is our generating model in this case: common correlation but different rates depending on the state of our multistate character.