The ‘gls.ancova’ function is an adaptation of the standard ANCOVA procedure so as to allow testing differences in intercepts among groups while accounting for phylogenetic relatedness among species. This adaptation consists of using the phylogenetic tree to describe the covariance among data points, thus rendering the procedure to be a phylogenetic ANCOVA (‘pANCOVA’).
This function can also be used to test the assumption of homogeneity of slopes that underpins ANCOVA.
The implementation of gls.ancova uses standard generalized least-squares procedures only. The advantage of this approach is that it is unbiased for sample size, is expected to appropriately penalize for increased parameterization using degrees of freedom, and can be generalized towards more complex models by including additional indicator variables.
When generalizing the procedure towards more complex models, this function can be used to test whether including a separate intercept and/or slope for a particular group provides a better fit to the data than a less parameterized model. In such a case, this procedure can thus test whether individual species, or a group of species, deviate significantly from allometric predictions.
Reference: Smaers JB & Rohlf FJ. Testing species’ deviations from allometric predictions using the phylogenetic regression. Evolution. 70 (5): 1145-1149.
Below two case studies that exemplify the use of pANCOVA.
As a first step, we’ll load the data needed to perform the analysis.
require(phytools); require(phytools); require(geiger); require(nlme); require(evomap)
data<-primateData; tree<-primateTree
Y<-"Brain"; X<-"Body"
data<-data[,c(which(colnames(data)==Y),which(colnames(data)==X)),drop=F]
data<-na.omit(data); data<-log(data)
tree<-treedata(tree,data,sort=T,warnings=T)$phy #match the tree to the data
data<-as.data.frame(treedata(tree,data,sort=T,warnings=T)$data) #match the data to the tree
colnames(data)<-c("Dependent","Independent") #Note that naming the variables 'Dependent' and 'Independent' serves only to standardize the procedure below across different data sets.
Case study 1: Platyrrhines versus catarrhines
Research question: Do platyrrhines (New World monkeys) indicate a different brain to body allometric scaling pattern than catarrhines (Old World monkeys and apes)?
Set groups to be compared:
Catarrhines<-getTips(tree,findMRCA(tree,c("Homo_sapiens","Allenopithecus_nigroviridis")))
Platyrrhines<-getTips(tree,findMRCA(tree,c("Callicebus_moloch","Alouatta_palliata")))
Plot the two groups being compared:
plot(data$Dependent~data$Independent,col="white",pch=19,xlab="", ylab="",asp=1,cex.lab=2)
pGLS.plotGrade("Dependent","Independent",data,tree,model="BM",group=Catarrhines,col="green",lwd=5,cex=1,pch=19)
pGLS.plotGrade("Dependent","Independent",data,tree,model="BM",group=Platyrrhines,col="grey",lwd=5,cex=1,pch=19)
Prepare indicator variables by allocating the groups that are compared to a different factor element:
#For differences in slope:
grpS<-rep("A",length(rownames(data)))
grpS[Catarrhines]<-"B"
grpS<-as.factor(grpS)
names(grpS)<-rownames(data)
#For differences in intercept:
grpI<-rep("A",length(rownames(data)))
grpI[Catarrhines]<-"B"
grpI<-as.factor(grpI)
names(grpI)<-rownames(data)
pANCOVA tests a model with fewer parameters against a model with more parameters. The baseline test comprises testing against a model with one slope and one intercept:
Model<-model.matrix(as.formula(Dependent~Independent),data)
Three types of models can be conceived that have more parameters than the baseline model (those that vary slope but not intercept; vary intercept but not slope; vary both):
#(1) Differences in slopes, holding intercept constant: Model_S<-model.matrix(as.formula(Dependent~grpS:Independent),data) #(2) Differences in intercept, holding slopes constant: Model_I<-model.matrix(as.formula(Dependent~grpI + Independent),data) #(3) Differences in slopes and differences in intercept: Model_SI<-model.matrix(as.formula(Dependent~grpI + grpS:Independent),data)
The following tests each of the above three models against the baseline model of one slope and one intercept:
#(1) Differences in slopes, holding intercept constant: gls.ancova(Dependent~Independent,vcv(tree),Model,Model_S) #(2) Differences in intercept, holding slopes constant: gls.ancova(Dependent~Independent,vcv(tree),Model,Model_I) #(3) Differences in slopes and differences in intercept: gls.ancova(Dependent~Independent,vcv(tree),Model,Model_SI)
Additionally, one can also test against models other than the baseline one slope – one intercept model, as long as the reduced model (listed first in the gls.ancova function) is the one with fewer parameters. For example, one could test the model that varies intercept only against the model that varies both intercept and slope :
gls.ancova(Dependent~Independent,vcv(tree),Model_I,Model_SI)
Here none of these results are significant. Therefore, we accept the model with fewer parameters. In this case, this is the one slope – one intercept model
Conclusion: Platyrrhines and catarrhines do not differ in their allometric patterning of brain size and body size.
Case study 2: Humans versus non-human primates
Research question: Do humans indicate a different brain to body scaling pattern than non-human primates?
Note that this case study allocates a single species to a particular group. In doing so, this case study tests whether a single species is a significant outlier in the allometric scaling relationship between brain size and body size. Such a test can only involve testing for a difference in intercept (i.e. whether the group ‘human’ has a different mean after removing variation in the covariates).
Set groups to be compared:
Human<-which(tree$tip.label=="Homo_sapiens")
nonHumanPrimates<-setdiff(getTips(tree,findMRCA(tree,c("Homo_sapiens","Callicebus_moloch"))),Human)
Plot the two groups being compared:
Color=rep("white",length(tree$tip.label)) Color[Human]<-"green"
plot(data$Dependent~data$Independent,col=Color,pch=19,xlab="", ylab="",asp=1,cex=2)
pGLS.plotGrade("Dependent","Independent",data,tree,model="BM",group=nonHumanPrimates,col="grey",lwd=5,cex=1,pch=19)
Set group allocation for the intercepts:
grpI<-rep("A",length(rownames(data)))
grpI[Human]<-"B"
grpI<-as.factor(grpI)
names(grpI)<-rownames(data)
Test against the baseline one-slope one-intercept model:
Model<-model.matrix(as.formula(Dependent~Independent),data) Model_I<-model.matrix(as.formula(Dependent~grpI + Independent),data) gls.ancova(Dependent~Independent,vcv(tree),Model,Model_I)
Conclusion: Humans are significant outliers in the brain to body scaling relationship.
The above procedures can be extrapolated to testing any combination of models with any combination of parameters (though note that the reduced model should always have fewer parameters). Models of evolution can be incorporated by adjusting the variance-covariance matrix of the tree according to relevant model parameters (e.g. the lambda parameter).
smaerslab.com 