Thursday, August 18, 2011

New phylogenetic ANOVA function with posthoc tests

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.

9 comments:

  1. 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?

    ReplyDelete
  2. Hi 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

    ReplyDelete
  3. Hi 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.

    ReplyDelete
  4. Hi Lian,
    I have a simple question: what is the name of your posthoc test? Is Tukey the default?
    Thanks!

    ReplyDelete
    Replies
    1. Hi Sanuel.
      phylANOVA 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

      Delete
  5. Hey Liam,

    Is 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

    ReplyDelete
    Replies
    1. Hi Rick.

      You 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

      Delete
  6. Hi Liam,

    Is 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)

    ReplyDelete