I just posted a new function to conduct the phylogenetic ANOVA of Garland et al. (1993). Luke Mahler has been doing phylogenetic ANOVA, but wanted to add posthoc comparison of the means between groups. I had previously programmed this in C - but I thought it would be easy enough to do in R as well. Direct link to the code is here. Note that the function, phylANOVA(), borrows a little bit from Luke Harmon's "geiger" function phy.anova(); and the function tTests() (which phylANOVA() calls internally) borrows code from the R "stats" function pairwise.t.test().
The phylogenetic ANOVA is pretty straightforward. First, we fit our ANOVA model (using lm()), and then we compute the F-statistic (using anova()). ANOVA, like most statistical tests, assumes independence - so the P-value returned by anova() will be incorrect. Garland et al.'s simple solution was to obtain the null distribution by simulation, which we can do by first calling X<-fastBM() a single time, and then computing the value of F for each column of X.
The posthoc tests are only a little bit more complicated. Here, we just compute a matrix of t-values for the real data and for each simulated data vector. Note that, for now, the method computes a single pooled standard-deviation which is used for all groups. For each pairwise comparison, we then count the number of simulated t-values with absolute values (for a two-sided test) that are larger than our corresponding empirical t. At the end, though, we have also conducted a large number of tests and that needs to be taken into consideration to control the experiment-wise error rate.
Fortunately in R this is straightforward to accomplish using the p.adjust() function in the "stats" package. p.adjust() implements a whole slew of methods for multiple test correction, and all of them can be called from phylANOVA(...,posthoc=TRUE,p.adj). For instance, standard Bonferroni correction is called by phylANOVA(...,p.adj="bonferroni") and sequential-Bonferroni (the default, also called the Holm-Bonferroni method) is called by phylANOVA(...,p.adj="holm").
Luke reports that the function seems to work properly.
Hi Liam, thanks for making this available; it looks very useful. Does it assume that the trait is evolving under BM evolution? I'd like to use it on some data that seem to best fit OU (having run fitContinuous on it; branch lengths of tree scaled using the method of Brusatte et al. 2008 using Graeme Lloyd's script). Is this a problem for phylANOVA?
ReplyDeleteHi Roger. I think an implicit assumption is BM evolution of the residual error controlling for the group effect. (A strict BM assumption wouldn't make any sense because under BM there would be no tendency for species in the same group to have similar phenotypes.) I hope that is helpful, but I would encourage you to check out the primary literature on this method for more specifics. Best of luck, Liam
ReplyDeleteHi Liam, thanks for the great code! I thought you might be interested to know that the post-hoc tests fail when one of the groups is represented by a single case (it returns NA for all pairs, even ones with multiple cases per group). Not sure if there is an easy way around this.
ReplyDeleteIs there a solution to this yet?
DeleteHello, i've seen several people with this problem but haven't seen a solution yet, is there one? Thank you
DeleteHi I am actually having this same issue right now and would love to hear if there is a solution!
DeleteStill having the same problem here
DeleteHi Lian,
ReplyDeleteI have a simple question: what is the name of your posthoc test? Is Tukey the default?
Thanks!
Hi Sanuel.
DeletephylANOVA computes t-values for all pairwise comparisons (although it could also just use the raw difference in the mean - it wouldn't matter), then obtains p-values for each comparison via simulation. To correct for multiple tests it uses p.adjust and the default method ("holm") is the default, although this can be changed.
- Liam
Hi Everyone,
DeleteAlso still looking for a solution for this.
Thank you!
Hey Liam,
ReplyDeleteIs it possible to do a phylogentic 2-way ANOVA? I can't seem to get either aov.phylo or phylANOVA to work.
Thanks,
Rick Simpson
Hi Rick.
DeleteYou might want to do this using gls in the nlme package, i.e.:
fit<-gls(y~x1+x2+x1*x2,correlation=corBrownian(1,tree))
for instance. Search for gls on my blog or the web for more information.
- Liam
Liam, I searched your blog for notes on how to conduct a post hoc test for a 2-way ANOVA and couldn't find anything. I can get the gls code to run but still need a way to do contrasts. Thanks!
DeleteI would search for how to do a post-hoc test in nlme or with the function gls. This is not a phytools function, so I don't know - but I suspect that methods have been written for precisely this purpose. - Liam
DeleteHi Liam,
ReplyDeleteIs there a way to perform a phylogenetic ANCOVA?
The following sample function doesn't work since the formula is required to be in the form data~group
aov.phylo(data~groups + covariate, tree, nsim=50)
Dear Liam,
ReplyDeleteFirst of all, thanks a lot for this blog. I am very new in R, and it has helped me a lot.
I wonder if there is a way to include more than one individual per species/tip and run a two-way ANOVA (with post hoc tests) using phylANOVA. Specifically, I’m trying to test for interspecific differences in a trait with species and (e.g.) reproductive status as factors.
Thanks,
Henrique