A hasty and incomplete intro to comparative methods

comparative-intro.R

# PICS and PGLS - oh my!
# Will Pearse - 2018-04-03

# Headers
library(ape)
library(caper)
library(geiger)
library(phytools)
library(mvMORPH)

#############################
# 'The Comparative Method' ##
#############################
# Simulate a phylogeny under a Yule model
tree <- sim.bdtree()

# Simulate a trait under a Brownian motion model
trait <- sim.char(tree, .5, model="BM", root=2)[,,1]
#...picking parameters to make the plotting easier

# Plot that trait
plot(tree, show.tip.label=FALSE)
tiplabels(pch=20, cex=abs(trait), col=ifelse(trait>0, "black","red"))

# Now let's have a response trait that is a function of that trait
non.phylo.response <- trait * 2 + rnorm(length(trait))
plot(non.phylo.response ~ trait)

# Do an (OLD) PIC-based analysis
pic.resp <- pic(non.phylo.response, tree)
pic.trait <- pic(trait, tree)
plot(pic.resp ~ pic.trait)
summary(lm(pic.resp ~ pic.trait))
abline(lm(pic.resp ~ pic.trait), col="red", lwd=3)
#...great, a positive correlation, but these parameters aren't perfect...

# Do a (NEW) PGLS-based analysis
data <- data.frame(species=names(trait), resp=non.phylo.response, trait=trait)
c.data <- comparative.data(tree, data, species)
summary(pgls(resp ~ trait, data=c.data, lambda="ML"))
#...two things to notice: (1) the coefficients are correct
#...                      (2) the lambda isn't significant. WHY?

# Let's try it again, but now with phylogenetically-correlated errors...
phylo.response <- trait * 2 + sim.char(tree, .25, model="BM", root=0)[,,1]
pic.phy.response <- pic(phylo.response, tree)
plot(pic.phy.response ~ pic.trait) #...hmmm, this looks different...
summary(lm(pic.phy.response ~ pic.trait))
#...hmmm, these parameters look... better...
data$phy.resp <- phylo.response
c.data <- comparative.data(tree, data, species)
summary(pgls(phy.resp ~ trait, data=c.data, lambda="ML"))
#...still fine. But our lambda value is different...
#...(if you get fitting errors,) SCALE YOUR DATA

# What have we learned?
# (1) Use PGLS, not PICs
# (2) The Comparative Method controls for residual phylogenetic error
#     (Read Uyeda et al. 2018 bioRxiv)

# A brief digression on Pagel's "other" transformations
# - the true tree is in red
tree <- sim.bdtree(n=30)
par(mfrow=c(3,3))
# Lambda - Brownian motion
plot(rescale(tree,"lambda",0.01), show.tip.label=FALSE, main=expression(paste("Pagel's ",lambda," = 0.01")))
plot(rescale(tree,"lambda",0.5), show.tip.label=FALSE, main=expression(paste("Pagel's ",lambda," = 0.5")))
plot(rescale(tree,"lambda",1), show.tip.label=FALSE, main=expression(paste("Pagel's ",lambda," = 1")),edge.col="red")
# Delta - accelerated/OU evolution
plot(rescale(tree,"delta",0.01), show.tip.label=FALSE, main=expression(paste("Pagel's ",delta," = 0.01")))
plot(rescale(tree,"delta",1), show.tip.label=FALSE, main=expression(paste("Pagel's ",delta," = 1")),edge.col="red")
plot(rescale(tree,"delta",2), show.tip.label=FALSE, main=expression(paste("Pagel's ",delta," = 2")))
# Kappa - punctuated evolution and/or ecological speciation
plot(rescale(tree,"kappa",0), show.tip.label=FALSE, main=expression(paste("Pagel's ",kappa," = 0")))
plot(rescale(tree,"kappa",1), show.tip.label=FALSE, main=expression(paste("Pagel's ",kappa," = 1")),edge.col="red")
plot(rescale(tree,"kappa",2), show.tip.label=FALSE, main=expression(paste("Pagel's ",kappa," = 2")))