Tuesday, October 16, 2012

Estimating Pagel's λ for a set of traits in a matrix or data frame

Earlier, I posted about estimating Pagel's λ jointly for a set of traits. A phytools user recently submitted the following comment:

One question I had was how to calculate lambda on multiple traits? I have code that gets me close in that it outputs the lambda values but not the log values or p-value.

This could be done, of course, using a for loop, something like:
> require(phytools); require(geiger)
> # simulate tree & data
> tree<-pbtree(n=100)
> X<-fastBM(lambdaTree(tree,0.7),nsim=10)
> lambda<-matrix(NA,ncol(X),4,dimnames=list(NULL, c("lambda","logL","logL0","P")))
> if(is.data.frame(X)) X<-as.matrix(X)
> for(i in 1:nrow(lambda))
lambda[i,]<-unlist(phylosig(tree,X[,i],method="lambda", test=T))
> lambda
        lambda      logL     logL0            P
[1,] 0.8176128 -176.9442 -199.7361 1.462441e-11
[2,] 0.5321105 -183.5490 -192.2505 3.023510e-05
[3,] 0.7868250 -177.3286 -200.6997 8.096301e-12
[4,] 0.6772734 -182.8100 -198.7038 1.719778e-08
[5,] 0.4893630 -183.6516 -193.9559 5.634189e-06
[6,] 0.7477579 -167.3094 -194.6396 1.432188e-13
[7,] 0.7160706 -182.1853 -207.2614 1.422751e-12
[8,] 0.6029530 -184.5244 -200.5472 1.505981e-08
[9,] 0.6654206 -181.8013 -201.3363 4.088944e-10
[10,] 0.7553931 -182.0996 -200.9483 8.261227e-10

Alternatively, however, we could also do the following using sapply:
> X<-X[tree$tip.label,]
> X<-as.data.frame(X)
> lambda<-t(sapply(X,phylosig,tree=tree,method="lambda", test=T))
[1] "x has no names; assuming x is in the same order as tree$tip.label"
...
> lambda
   lambda    logL      logL0     P          
V1  0.8176128 -176.9442 -199.7361 1.462441e-11
V2  0.5321105 -183.549  -192.2505 3.02351e-05
V3  0.786825  -177.3286 -200.6997 8.096301e-12
V4  0.6772734 -182.81   -198.7038 1.719778e-08
V5  0.489363  -183.6516 -193.9559 5.634189e-06
V6  0.7477579 -167.3094 -194.6396 1.432188e-13
V7  0.7160706 -182.1853 -207.2614 1.422751e-12
V8  0.602953  -184.5244 -200.5472 1.505981e-08
V9  0.6654206 -181.8013 -201.3363 4.088944e-10
V10 0.7553931 -182.0996 -200.9483 8.261227e-10

That's it.

No comments:

Post a Comment