Pages

Saturday, November 17, 2012

Testing for Pagel's λ < 1.0

A recent commenter on the phytools blog asks:I’m working on a comparative phylogenetic project and want to show that my data has phylogenetic signal. Using Pagel’s lambda I get lambda – 0.854. Is this value significant enough in that it’s safe to model the variable as Brownian motion?

I interpret this to mean "how can we test the alternative hypothesis that Pagel's λ is less then 1.0; against the null that λ=1.0." This is easy if we accept the contention that our likelihood ratio should be distributed as a Χ2 with 1 degree of freedom (for one parameter, λ, estimated in the more complex model), an asymptotic property of likelihoods. (There may be some very good reasons - which we'll set aside for the moment - to question this assertion, for instance see a recent article by Carl Boettiger et al.)

This can be done most naturally using the fitContinuous function in 'geiger' - but we can also do it using entirely phytools functions, and since I developed and maintain phytools (not geiger), I'll demo that option.

First, let's simulate a tree & data for the demo, using the λ=0.854 obtained by the commenter. We will use geiger::lambdaTree for this part:

> library(phytools); library(geiger)
> # simulate tree
> tree
> # simulate data with lambda=0.854
> x
OK, now let's fit both models, compute the likelihood ratio, and then compare it to a Χ2 distribution with one degree of freedom. We need to use a little bit of a work-around to fit the BM model using brownie.lite:

An additional word of caution, however. Sampling error in the estimation of species means can also have the effect of depressing phylogenetic signal - so if the species means are estimated with a limited number of samples per species, it might be wise to conduct this analysis including estimated within species sampling variances. These can be supplied using the argument se in both phylosig and brownie.lite.

Yes - "se" is the standard error of the mean for each species. This is normally calculated for each species as the within species SD divided by the square root of the sample size.

In the case above, we would not reject the null - but the LR was not 15 (as you say) in this particular stochastic simulation. (The LR for a P of 0.135 would be around 2.24, or so. Since this is a simulation - individual results will vary.) The rejection level for alpha=0.05 with df=1 is about 3.84.

Hi Liam, to follow up on this interesting thread, what would be your suggestion to test the other end of the spectrum, that lambda is (or not) significantly different from 0? Would a comparison with a white noise model do the trick? Or using lambdaTree() in geiger to transform the phylogeny to a star tree?I enjoy the blog!CheersAlejandro

About this blog

This web-log chronicles the development of new tools for phylogenetic analyses in the phytools R package. Unless you a reading a very recent page of the blog, I recommend that you install the latest CRAN version of phytools (or latest beta release) before attempting to replicate any of the analyses of this site. That is because the linked functions may be archived, and very likely have been replaced by newer versions.