Post-hoc tests for a generalized phylogenetic ANOVA or ANCOVA

An R phylogenetics user recently asked me the following aboug phylogenetic generalized least squares ANOVA using nlme::gls.

”[I fit the model] pgls.lambda <-gls(Y~X1+X2+X3,correlation=corPagel(1,tree)) which I think is fine. Do you know a way to apply a post-hoc test to such a model? I’ve tried multcomp::glht(pgls.lambda,linfct=mcp("tension"="Tukey")) but it doesn’t work.”

Actually, I had not done this before – although in principle it should be relatively straightforward to compute t or z values from the estimated coefficients and their standard errors, and then compare them in the normal way while controlling for experiment-wise error.

On the other hand, the internet had another solution, which involved the simple act of creating several new methods for the "gls" object class such that multcomp::glht worked.

Here are the methods.

model.matrix.gls<-function(object,...){
	model.matrix(terms(object),data=getData(object),...)
}
model.frame.gls<-function(object,...){
	model.frame(formula(object),data=getData(object),...)
}
terms.gls<-function(object,...){
	terms(model.frame(object),...)
}

To see how this works, I will reiterate the phylogenetic ANCOVA of Chapter 3 from my recent book with Luke Harmon. For a simpler ANOVA model, this method should work just as well.

## load packages
library(phytools)
library(nlme)
library(multcomp)
## read data from file
primate.data<-read.csv("http://www.phytools.org/Rbook/3/primateEyes.csv",
    row.names=1,stringsAsFactors=TRUE)
primate.tree<-read.tree("http://www.phytools.org/Rbook/3/primateEyes.phy")
## fit generalized
spp<-rownames(primate.data)
corBM<-corBrownian(phy=primate.tree,form=~spp)
primate.ancova<-gls(log(Orbit_area)~log(Skull_length)+
    Activity_pattern,data=primate.data,
    correlation=corBM)
primate.ancova

## Generalized least squares fit by REML
##   Model: log(Orbit_area) ~ log(Skull_length) + Activity_pattern 
##   Data: primate.data 
##   Log-restricted-likelihood: 48.22693
## 
## Coefficients:
##               (Intercept)         log(Skull_length)   Activity_patternDiurnal 
##                -0.7310348                 1.4717247                -0.1060187 
## Activity_patternNocturnal 
##                 0.2935941 
## 
## Correlation Structure: corBrownian
##  Formula: ~spp 
##  Parameter estimate(s):
## numeric(0)
## Degrees of freedom: 90 total; 86 residual
## Residual standard error: 0.2933082

Now let’s do our post-hoc tests using glht.

post.hoc<-glht(primate.ancova,linfct=mcp(Activity_pattern="Tukey"))
summary(post.hoc)

## 
## 	 Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: gls(model = log(Orbit_area) ~ log(Skull_length) + Activity_pattern, 
##     data = primate.data, correlation = corBM)
## 
## Linear Hypotheses:
##                             Estimate Std. Error z value Pr(>|z|)    
## Diurnal - Cathemeral == 0   -0.10602    0.07773  -1.364    0.347    
## Nocturnal - Cathemeral == 0  0.29359    0.12139   2.419    0.039 *  
## Nocturnal - Diurnal == 0     0.39961    0.09706   4.117   <0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

In another post I’ll look at whether this post hoc test has appropriate experiment-wise type I error.