Monday, May 21, 2012

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.

1 comment:

  1. Le magasin a résolu des problèmes pour moi, le service est très patient et c’est un magasin parfait.parfait gucci homme chaussures Quand j'ai reçu cette montre, j'ai pensé que c'était une bonne réplique.parfait gucci tracksuits J'aime cette boutique en ligne, cette montre est simple et généreuse. Très satisfait de cet achat.

    ReplyDelete