Friday, November 11, 2011

New version of phyl.RMA() with fixed lambda

I just posted a new version of my R function for phylogenetic reduced major axis regression (RMA), called phyl.RMA(). Since the RMA slope is merely computed as the variance ratio, this function just computes phylogenetic variances and then takes their ratio. The function also (for method="lambda") first performs joint estimation of Pagel's λ (e.g., see Freckleton et al. 2002). Direct link to the code is here.

The new feature of this version of the function (as requested by a user) is merely that the value of λ can be fixed by the user rather than estimated. This is done by calling phyl.RMA(...,fixed=TRUE,lambda=XX). So, for instance, try the following:

> require(geiger)
> require(phytools)
> source("phyl.RMA.R")
> tree<-drop.tip(birthdeath.tree(b=1,d=0,taxa.stop=201),"201")
> X<-fastBM(lambdaTree(tree,0.7),nsim=2)
> phyl.RMA(X[,1],X[,2],tree,method="lambda")
$RMA.beta
[1] -0.0866432 1.0590274
$V
x y
x 0.97406967 -0.03095092
y -0.03095092 1.09245711
$lambda
[1] 0.734092
$logL
[1] -770.6879
$resid
[,1]
21 -5.863653843
101 -5.841695560
...


If we wanted to know, say, if the estimated value of λ (here ~0.73) was significantly different than 0.0 (or any other value), we can fix λ at 0.0 and recalculate:

> phyl.RMA(X[,1],X[,2],tree,method="lambda",fixed=T,lambda=0)
$RMA.beta
[1] -0.7000868 0.9720490

$V
x y
x 0.75527418 -0.09025328
y -0.09025328 0.71364296

$lambda
[1] 0

$logL
[1] -845.7699

$resid
[,1]
21 -5.01116336
101 -4.87370275
...


Cool!

Unfortunately, it has just been reported to me that we have had an escaped anole from the common garden Anolis carolinensis rearing experiment that is presently taking place in my lab. Time to go lizard hunting!

1 comment:

  1. Turns out that the tales of escape had been greatly exaggerated - the missing lizard was found exactly where it should've been, in his/her cage.

    ReplyDelete

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.