Trait challenge solution

First, read in the tree and trait data. We will also check to make sure the names match.

library(geiger)

## Loading required package: ape

gruntTree<-read.tree("~/Desktop/grunts/grunts.phy")
gruntData<-read.csv("~/Desktop/grunts/grunts.csv", row.names=1)
name.check(gruntTree, gruntData)

## [1] "OK"

Now we will try to do question 1.

1. (P)GLS. Is there a correlation between traits 1 and 2 (ignoring the phylogeny - that is, using a standard correlation). Then use PGLS to determine if there is an evolutionary correlation between traits 1 and 2 in the dataset. Before figuring out your final results, use AIC to determine if it is a good idea to use an OU or Pagel’s lambda for the PGLS.

library(nlme)

plot(trait2~trait1, data=gruntData)

pglsResult<-gls(trait1~trait2, cor=corBrownian(phy=gruntTree), data=gruntData, method="ML")
pglsResultOU<-gls(trait1~trait2, cor=corMartins(1, phy=gruntTree, fixed=F), data=gruntData, method="ML")
pglsResultPagel<-gls(trait1~trait2, cor=corPagel(1, phy=gruntTree, fixed=F), data=gruntData, method="ML")

anova(pglsResult, pglsResultOU, pglsResultPagel)

##                 Model df      AIC      BIC     logLik   Test  L.Ratio
## pglsResult          1  3 59.81018 65.54625 -26.905088
## pglsResultOU        2  4 30.42696 38.07506 -11.213482 1 vs 2 31.38321
## pglsResultPagel     3  4 22.83320 30.48129  -7.416599
##                 p-value
## pglsResult
## pglsResultOU     <.0001
## pglsResultPagel

# the main one is the lowest AIC, choose that

anova(pglsResultPagel)

## Denom. DF: 48
##             numDF  F-value p-value
## (Intercept)     1 0.004714  0.9455
## trait2          1 6.053490  0.0175

# yes there is a relationship.

Now try problem 2:

2. Models for continuous character evolution. Fit three models (BM, OU, and EB) to trait1, trait2, and trait3 in your datamatrix. Summarize your results in a table.

results<-matrix(nrow=3, ncol=3)

for(i in 3:5) {
  xx <- gruntData[,i]
  names(xx) <- rownames(gruntData)
  bm<-fitContinuous(gruntTree, xx)
  ou<-fitContinuous(gruntTree, xx, model="OU")
  eb<-fitContinuous(gruntTree, xx, model="EB")

  aicc<-c(bm$opt$aicc, ou$opt$aicc, eb$opt$aicc)
  aiccD <- aicc - min(aicc)
  aw <- exp(-0.5 * aiccD)
  aiccW <- aw/sum(aw)

  results[i-2,]<-aiccW

}

## Loading required package: parallel

## Warning in fitContinuous(gruntTree, xx, model = "OU"): Parameter estimates appear at bounds:
##  alpha

## Warning in fitContinuous(gruntTree, xx, model = "OU"): Parameter estimates appear at bounds:
##  alpha

## Warning in fitContinuous(gruntTree, xx, model = "OU"): Parameter estimates appear at bounds:
##  alpha

rownames(results)<-colnames(gruntData)[3:5]
colnames(results)<-c("bm", "ou", "eb")

# Table of AICc weights
results

##               bm        ou        eb
## trait1 0.4583683 0.3940379 0.1475937
## trait2 0.5292929 0.3002759 0.1704312
## trait3 0.5230986 0.3084647 0.1684367

Last, number 3:

3. Ancestral character states. Reconstruct ancestral character states for the column habitat_names (which has reef versus non-reef) and trait1.

library(phytools)

## Loading required package: maps

## Warning in FUN(X[[35L]], ...): failed to assign RegisteredNativeSymbol for
## getData to getData since getData is already defined in the 'phangorn'
## namespace

habitat<-gruntData[,2]
names(habitat)<-rownames(gruntData)
fitER<-ace(habitat,gruntTree,model="ER",type="discrete")

cols<-setNames(palette()[1:length(unique(habitat))],sort(unique(habitat)))

plotTree(gruntTree,type="fan",fsize=0.8,ftype="i")

## setEnv=TRUE for this type is experimental. please be patient with bugs

nodelabels(node=1:gruntTree$Nnode+Ntip(gruntTree),
    pie=fitER$lik.anc,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(habitat,sort(unique(habitat))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
    y=-max(nodeHeights(gruntTree)),fsize=0.8)

t1 <- gruntData[,"trait1"]
names(t1) <- rownames(gruntData)
fit <- fastAnc(gruntTree, t1, vars=T, CI=T)
obj<-contMap(gruntTree,t1,plot=FALSE)
plot(obj,type="fan",legend=0.7*max(nodeHeights(gruntTree)))