Comparative methods in R - Ilhabela

Continuous & discrete character evolution using the threshold model

Jul 5, 2015
Ancestral character estimation under the threshold model

In this tutorial, we're going to cover two different analyses:

  1. We're going to estimate the correlation between two binary characters, as well as between a binary & a continuous character, assuming a threshold model of character evolution for our discrete character using both threshBayes in the phytools package; and
  2. We're going to reconstruct the ancestral states of a multistate, naturally ordered character on the tree assuming the threshold model.

If you want to following along, but do not have a fast enough computer to run the MCMC, you can download and load the .Rdata file thresh.Rdata and then run everything except the MCMC.

Correlation between discrete & continuous characters under the threshold model

First, we're going to use Bayesian MCMC method to sample the evolutionary parameters of our model from their joint posterior probability distribution. These parameters include evolutionary rates (σ2s) for each continuous character; and the correlation between our discrete and continuous character. Since this is a binary model with a single threshold - we don't have to worry about the relative positions of the thresholds.

We'll start by simulating a tree & two correlated Brownian variables on the tree.

For illustrative purposes, let's start by using threshBayes on two continuous characters & show that the result lines up closely with our MLE of the correlation.

## first, load phytools
library(phytools)

## Loading required package: ape
## Loading required package: maps

tree<-pbtree(n=50,scale=1)
r<-0.75 # simulate using a high correlation
V<-matrix(c(1,r,r,1),2,2)
X<-sim.corrs(tree,V)
## set some parameters for the MCMC
sample<-500
ngen<-200000 ## in 'real' studies this should be larger
burnin<-0.2*ngen
## we can start by running our Bayesian MCMC on the original data
mcmcCC<-threshBayes(tree,X,types=c("cont","cont"),
    ngen=ngen,control=list(sample=sample))

## gen 1000
## gen 2000
## gen 3000
## gen 4000
## gen 5000
## gen 6000
## gen 7000
## gen 8000
## gen 9000
## gen 10000
## gen 11000
## gen 12000
## gen 13000
## gen 14000
## gen 15000
## gen 16000
## gen 17000
## gen 18000
## gen 19000
## gen 20000
## gen 21000
## gen 22000
## gen 23000
## gen 24000
## gen 25000
## gen 26000
## gen 27000
## gen 28000
## gen 29000
## gen 30000
## gen 31000
## gen 32000
## gen 33000
## gen 34000
## gen 35000
## gen 36000
## gen 37000
## gen 38000
## gen 39000
## gen 40000
## gen 41000
## gen 42000
## gen 43000
## gen 44000
## gen 45000
## gen 46000
## gen 47000
## gen 48000
## gen 49000
## gen 50000
## gen 51000
## gen 52000
## gen 53000
## gen 54000
## gen 55000
## gen 56000
## gen 57000
## gen 58000
## gen 59000
## gen 60000
## gen 61000
## gen 62000
## gen 63000
## gen 64000
## gen 65000
## gen 66000
## gen 67000
## gen 68000
## gen 69000
## gen 70000
## gen 71000
## gen 72000
## gen 73000
## gen 74000
## gen 75000
## gen 76000
## gen 77000
## gen 78000
## gen 79000
## gen 80000
## gen 81000
## gen 82000
## gen 83000
## gen 84000
## gen 85000
## gen 86000
## gen 87000
## gen 88000
## gen 89000
## gen 90000
## gen 91000
## gen 92000
## gen 93000
## gen 94000
## gen 95000
## gen 96000
## gen 97000
## gen 98000
## gen 99000
## gen 100000
## gen 101000
## gen 102000
## gen 103000
## gen 104000
## gen 105000
## gen 106000
## gen 107000
## gen 108000
## gen 109000
## gen 110000
## gen 111000
## gen 112000
## gen 113000
## gen 114000
## gen 115000
## gen 116000
## gen 117000
## gen 118000
## gen 119000
## gen 120000
## gen 121000
## gen 122000
## gen 123000
## gen 124000
## gen 125000
## gen 126000
## gen 127000
## gen 128000
## gen 129000
## gen 130000
## gen 131000
## gen 132000
## gen 133000
## gen 134000
## gen 135000
## gen 136000
## gen 137000
## gen 138000
## gen 139000
## gen 140000
## gen 141000
## gen 142000
## gen 143000
## gen 144000
## gen 145000
## gen 146000
## gen 147000
## gen 148000
## gen 149000
## gen 150000
## gen 151000
## gen 152000
## gen 153000
## gen 154000
## gen 155000
## gen 156000
## gen 157000
## gen 158000
## gen 159000
## gen 160000
## gen 161000
## gen 162000
## gen 163000
## gen 164000
## gen 165000
## gen 166000
## gen 167000
## gen 168000
## gen 169000
## gen 170000
## gen 171000
## gen 172000
## gen 173000
## gen 174000
## gen 175000
## gen 176000
## gen 177000
## gen 178000
## gen 179000
## gen 180000
## gen 181000
## gen 182000
## gen 183000
## gen 184000
## gen 185000
## gen 186000
## gen 187000
## gen 188000
## gen 189000
## gen 190000
## gen 191000
## gen 192000
## gen 193000
## gen 194000
## gen 195000
## gen 196000
## gen 197000
## gen 198000
## gen 199000
## gen 200000

## we can compute our MLE of the correlation
V.mle<-phyl.vcv(X,vcv(tree),lambda=1)$R
V.mle[1,2]/sqrt(V.mle[1,1]*V.mle[2,2])

## [1] 0.6792364

## here is our "post burnin" mean from the posterior sample
mean(mcmcCC$par[(burnin/sample+1):nrow(mcmcCC$par),"r"])

## [1] 0.6722659

## plot our likelihood profile
plot(mcmcCC$par[,"logL"],type="l",xlab="generation",ylab="logL")

plot of chunk unnamed-chunk-1 

## plot our posterior density for the correlation
plot(density(mcmcCC$par[(burnin/sample+1):nrow(mcmcCC$par),
    "r"],bw=0.1),xlab="r",main="posterior density for r")
lines(c(r,r),c(0,1000),lty="dashed")

plot of chunk unnamed-chunk-2

OK, this was just a proof-of-concept exercise. Now let's try the same thing after having converted one of our trait values into a threshold character.

## convert X[,2] to a binary character
X[,2]<-as.numeric(sapply(X[,2],threshState,
    thresholds=setNames(c(0,Inf),0:1)))
mcmcCD<-threshBayes(tree,X,types=c("cont","disc"),
    ngen=ngen,control=list(sample=sample))

## gen 1000
## gen 2000
## gen 3000
## gen 4000
## gen 5000
## gen 6000
## gen 7000
## gen 8000
## gen 9000
## gen 10000
## gen 11000
## gen 12000
## gen 13000
## gen 14000
## gen 15000
## gen 16000
## gen 17000
## gen 18000
## gen 19000
## gen 20000
## gen 21000
## gen 22000
## gen 23000
## gen 24000
## gen 25000
## gen 26000
## gen 27000
## gen 28000
## gen 29000
## gen 30000
## gen 31000
## gen 32000
## gen 33000
## gen 34000
## gen 35000
## gen 36000
## gen 37000
## gen 38000
## gen 39000
## gen 40000
## gen 41000
## gen 42000
## gen 43000
## gen 44000
## gen 45000
## gen 46000
## gen 47000
## gen 48000
## gen 49000
## gen 50000
## gen 51000
## gen 52000
## gen 53000
## gen 54000
## gen 55000
## gen 56000
## gen 57000
## gen 58000
## gen 59000
## gen 60000
## gen 61000
## gen 62000
## gen 63000
## gen 64000
## gen 65000
## gen 66000
## gen 67000
## gen 68000
## gen 69000
## gen 70000
## gen 71000
## gen 72000
## gen 73000
## gen 74000
## gen 75000
## gen 76000
## gen 77000
## gen 78000
## gen 79000
## gen 80000
## gen 81000
## gen 82000
## gen 83000
## gen 84000
## gen 85000
## gen 86000
## gen 87000
## gen 88000
## gen 89000
## gen 90000
## gen 91000
## gen 92000
## gen 93000
## gen 94000
## gen 95000
## gen 96000
## gen 97000
## gen 98000
## gen 99000
## gen 100000
## gen 101000
## gen 102000
## gen 103000
## gen 104000
## gen 105000
## gen 106000
## gen 107000
## gen 108000
## gen 109000
## gen 110000
## gen 111000
## gen 112000
## gen 113000
## gen 114000
## gen 115000
## gen 116000
## gen 117000
## gen 118000
## gen 119000
## gen 120000
## gen 121000
## gen 122000
## gen 123000
## gen 124000
## gen 125000
## gen 126000
## gen 127000
## gen 128000
## gen 129000
## gen 130000
## gen 131000
## gen 132000
## gen 133000
## gen 134000
## gen 135000
## gen 136000
## gen 137000
## gen 138000
## gen 139000
## gen 140000
## gen 141000
## gen 142000
## gen 143000
## gen 144000
## gen 145000
## gen 146000
## gen 147000
## gen 148000
## gen 149000
## gen 150000
## gen 151000
## gen 152000
## gen 153000
## gen 154000
## gen 155000
## gen 156000
## gen 157000
## gen 158000
## gen 159000
## gen 160000
## gen 161000
## gen 162000
## gen 163000
## gen 164000
## gen 165000
## gen 166000
## gen 167000
## gen 168000
## gen 169000
## gen 170000
## gen 171000
## gen 172000
## gen 173000
## gen 174000
## gen 175000
## gen 176000
## gen 177000
## gen 178000
## gen 179000
## gen 180000
## gen 181000
## gen 182000
## gen 183000
## gen 184000
## gen 185000
## gen 186000
## gen 187000
## gen 188000
## gen 189000
## gen 190000
## gen 191000
## gen 192000
## gen 193000
## gen 194000
## gen 195000
## gen 196000
## gen 197000
## gen 198000
## gen 199000
## gen 200000

## what is our estimate of the correlation?
mean(mcmcCD$par[(burnin/sample+1):nrow(mcmcCD$par),"r"])

## [1] 0.5516791

## plot our likelihood profile
plot(mcmcCD$par[,"logL"],type="l",xlab="generation",
    ylab="logL")

plot of chunk unnamed-chunk-3 

## plot our estimated posterior density for the correlation
plot(density(mcmcCD$par[(burnin/sample+1):nrow(mcmcCD$par),
    "r"],bw=0.1),xlab="r",main="posterior density for r")
lines(c(r,r),c(0,1000),lty="dashed")

plot of chunk unnamed-chunk-4

Finally, let's analyze two binary characters under the same model. Remember, each time we're losing information about the underlying correlation between liabilities. Nonetheless….

## convert X[,2] to a binary character
X[,1]<-as.numeric(sapply(X[,1],threshState,
    thresholds=setNames(c(0,Inf),0:1)))
mcmcDD<-threshBayes(tree,X,types=c("disc","disc"),
    ngen=ngen,control=list(sample=sample))

## gen 1000
## gen 2000
## gen 3000
## gen 4000
## gen 5000
## gen 6000
## gen 7000
## gen 8000
## gen 9000
## gen 10000
## gen 11000
## gen 12000
## gen 13000
## gen 14000
## gen 15000
## gen 16000
## gen 17000
## gen 18000
## gen 19000
## gen 20000
## gen 21000
## gen 22000
## gen 23000
## gen 24000
## gen 25000
## gen 26000
## gen 27000
## gen 28000
## gen 29000
## gen 30000
## gen 31000
## gen 32000
## gen 33000
## gen 34000
## gen 35000
## gen 36000
## gen 37000
## gen 38000
## gen 39000
## gen 40000
## gen 41000
## gen 42000
## gen 43000
## gen 44000
## gen 45000
## gen 46000
## gen 47000
## gen 48000
## gen 49000
## gen 50000
## gen 51000
## gen 52000
## gen 53000
## gen 54000
## gen 55000
## gen 56000
## gen 57000
## gen 58000
## gen 59000
## gen 60000
## gen 61000
## gen 62000
## gen 63000
## gen 64000
## gen 65000
## gen 66000
## gen 67000
## gen 68000
## gen 69000
## gen 70000
## gen 71000
## gen 72000
## gen 73000
## gen 74000
## gen 75000
## gen 76000
## gen 77000
## gen 78000
## gen 79000
## gen 80000
## gen 81000
## gen 82000
## gen 83000
## gen 84000
## gen 85000
## gen 86000
## gen 87000
## gen 88000
## gen 89000
## gen 90000
## gen 91000
## gen 92000
## gen 93000
## gen 94000
## gen 95000
## gen 96000

## gen 97000
## gen 98000
## gen 99000
## gen 100000
## gen 101000
## gen 102000
## gen 103000
## gen 104000
## gen 105000
## gen 106000
## gen 107000

## gen 108000
## gen 109000
## gen 110000
## gen 111000
## gen 112000
## gen 113000
## gen 114000
## gen 115000
## gen 116000
## gen 117000
## gen 118000
## gen 119000
## gen 120000
## gen 121000
## gen 122000
## gen 123000
## gen 124000
## gen 125000
## gen 126000
## gen 127000
## gen 128000
## gen 129000
## gen 130000
## gen 131000
## gen 132000
## gen 133000
## gen 134000
## gen 135000
## gen 136000
## gen 137000
## gen 138000
## gen 139000
## gen 140000
## gen 141000
## gen 142000
## gen 143000
## gen 144000
## gen 145000
## gen 146000
## gen 147000
## gen 148000
## gen 149000
## gen 150000
## gen 151000
## gen 152000
## gen 153000
## gen 154000
## gen 155000
## gen 156000
## gen 157000
## gen 158000
## gen 159000
## gen 160000
## gen 161000
## gen 162000
## gen 163000
## gen 164000
## gen 165000
## gen 166000
## gen 167000
## gen 168000
## gen 169000
## gen 170000
## gen 171000
## gen 172000
## gen 173000
## gen 174000
## gen 175000
## gen 176000
## gen 177000
## gen 178000
## gen 179000
## gen 180000
## gen 181000
## gen 182000
## gen 183000
## gen 184000
## gen 185000
## gen 186000
## gen 187000
## gen 188000
## gen 189000
## gen 190000
## gen 191000
## gen 192000
## gen 193000
## gen 194000
## gen 195000
## gen 196000
## gen 197000
## gen 198000
## gen 199000
## gen 200000

## what is our estimate of the correlation?
mean(mcmcDD$par[(burnin/sample+1):nrow(mcmcDD$par),"r"])

## [1] 0.4689516

## plot our likelihood profile
plot(mcmcDD$par[,"logL"],type="l",xlab="generation",
    ylab="logL")

plot of chunk unnamed-chunk-5 

## plot our posterior sample for the correlation
plot(density(mcmcDD$par[(burnin/sample+1):nrow(mcmcDD$par),
    "r"],bw=0.1),xlab="r",main="posterior density for r")
lines(c(r,r),c(0,1000),lty="dashed")

plot of chunk unnamed-chunk-6

So we can see that we can see that threshBayes does a reasonable job at estimating the correlation between binary & continuous characters under the threshold model.

Finally, in 'real' empirical studies, it is always wise to ensure that you have a 'good' sample from the posterior distribution - for example, by computing the effective sample size of your posterior sample.

## let's use coda to compute effective sample sizes & HPDs
require(coda)

## Loading required package: coda

## first, analysis A
rA<-mcmcCC$par[(burnin/sample+1):nrow(mcmcCC$par),"r"]
class(rA)<-"mcmc"
effectiveSize(rA)

## var1 
##  321

HPDinterval(rA)

##          lower     upper
## var1 0.5230285 0.8187862
## attr(,"Probability")
## [1] 0.9501558

## second, analysis B
rB<-mcmcCD$par[(burnin/sample+1):nrow(mcmcCD$par),"r"]
class(rB)<-"mcmc"
effectiveSize(rB)

##     var1 
## 36.28832

HPDinterval(rB)

##          lower     upper
## var1 0.1365123 0.8757108
## attr(,"Probability")
## [1] 0.9501558

## finally, analysis C
rC<-mcmcDD$par[(burnin/sample+1):nrow(mcmcDD$par),"r"]
class(rC)<-"mcmc"
effectiveSize(rC)

##    var1 
## 34.5236

HPDinterval(rC)

##           lower     upper
## var1 0.08149168 0.8877881
## attr(,"Probability")
## [1] 0.9501558

In addition to this, we can also use Rthreshml in the Rphylip package to call threshml from the PHYLIP package internally. This is much faster & should produce similar etimates of the evolutionary covariances under most circumstances. It can also compute the covariances between multiple traits.

Ancestral character estimation under the threshold model

We can also do ancestral character reconstruction for discrete characters using the threshold model.

## simulate tree & data
n<-50
ngen<-100000
sample<-500
tree<-pbtree(n=n,scale=10)
x<-fastBM(tree,sig2=1,a=0.5,internal=TRUE)
th<-sapply(x,threshState,thresholds=setNames(c(0,0.5,2,Inf),
    letters[1:4]))
mcmc<-ancThresh(tree,th[1:n],ngen=ngen,sequence=letters[1:4],
    control=list(sample=sample,plot=FALSE))

## MCMC starting....
## gen 1000
## gen 2000
## gen 3000
## gen 4000
## gen 5000
## gen 6000
## gen 7000
## gen 8000
## gen 9000
## gen 10000
## gen 11000
## gen 12000
## gen 13000
## gen 14000
## gen 15000
## gen 16000
## gen 17000
## gen 18000
## gen 19000
## gen 20000
## gen 21000
## gen 22000
## gen 23000
## gen 24000
## gen 25000
## gen 26000
## gen 27000
## gen 28000
## gen 29000
## gen 30000
## gen 31000
## gen 32000
## gen 33000
## gen 34000
## gen 35000
## gen 36000
## gen 37000
## gen 38000
## gen 39000
## gen 40000
## gen 41000
## gen 42000
## gen 43000
## gen 44000
## gen 45000
## gen 46000
## gen 47000
## gen 48000
## gen 49000
## gen 50000
## gen 51000
## gen 52000
## gen 53000
## gen 54000
## gen 55000
## gen 56000
## gen 57000
## gen 58000
## gen 59000
## gen 60000
## gen 61000
## gen 62000
## gen 63000
## gen 64000
## gen 65000
## gen 66000
## gen 67000
## gen 68000
## gen 69000
## gen 70000
## gen 71000
## gen 72000
## gen 73000
## gen 74000
## gen 75000
## gen 76000
## gen 77000
## gen 78000
## gen 79000
## gen 80000
## gen 81000
## gen 82000
## gen 83000
## gen 84000
## gen 85000
## gen 86000
## gen 87000
## gen 88000
## gen 89000
## gen 90000
## gen 91000
## gen 92000
## gen 93000
## gen 94000
## gen 95000
## gen 96000
## gen 97000
## gen 98000
## gen 99000
## gen 100000

plotTree(tree,ftype="off")
tiplabels(pie=to.matrix(th[1:n],letters[1:4]),
    piecol=palette()[1:4],cex=0.3)
nodelabels(pie=mcmc$ace,piecol=palette()[1:4],cex=0.8)

plot of chunk unnamed-chunk-8

Theoretically, we can even use this method to identify the relative positions of thresholds:

colMeans(mcmc$par[(0.2*ngen/sample):(ngen/sample)+1,
    letters[1:4]])

##         a         b         c         d 
## 0.0000000 0.3308686 1.8797306       Inf

par(mfrow=c(2,1)); par(mar=c(4.1,5.1,2.1,2.1))
plot(density(mcmc$par[(0.2*ngen/sample):(ngen/sample)+1,
    letters[2]],bw=0.5),xlab="threshold a->b",main="")
lines(c(0.5,0.5),c(0,1000),lty="dashed")
plot(density(mcmc$par[(0.2*ngen/sample):(ngen/sample)+1,
    letters[3]],bw=0.5),xlab="threshold b->c",main="")
lines(c(2,2),c(0,1000),lty="dashed")

plot of chunk unnamed-chunk-9

Save the workspace:

save.image(file="thresh.Rdata")

Written by Liam J. Revell. Last updated July 5, 2015