Skip to content

Instantly share code, notes, and snippets.

@willpearse
Created April 13, 2018 19:50
Show Gist options
  • Save willpearse/5c0b0e88c8d81a8e1c1c7a75e7f58992 to your computer and use it in GitHub Desktop.
Save willpearse/5c0b0e88c8d81a8e1c1c7a75e7f58992 to your computer and use it in GitHub Desktop.
A hasty and incomplete intro to comparative methods
# 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")))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment