Fitting a discrete or hidden character multi-rate continuous trait evolution model using the discrete approximation in phytools

For a few months, phytools has had a function called fitmultiBM which uses the discrete approximation of Boucher & Démery (2016) to jointly model the evolution of a discrete and continuous trait in which the state of the former affects the rate of evolution of the latter.

This is conceptually similar to O’Meara et al. (2006), but in which we jointly model the evolution of our discrete and continuous traits – rather than first mapping the discrete character onto the tree (e.g., via stochastic mapping), and then fitting a continuous trait model in which \(\sigma^2\) depends on the mapped levels.

I’ve written a number of prior blog posts about this model (e.g., 1,2, 3, 4, and others), and I gave a talk about it (and other related analyses) at the Evolution meetings in Montreal last summer which can be viewed on YouTube here.

I genuinely think this is one of the most important new methods in phytools this year, so I wanted to revisit it to show just how powerful fitmultiBM can be both in fitting a discrete character dependent variable rate continuous trait evolution model, and also in identifying shifts in continuous character rate even where no specific discrete character dependence was postulated.

Note that I won’t focus on measuring parameter estimation via multiple simulations, which I did in a prior post here – but I’d encourage readers to reiterate the simulations & analyses of this post with different conditions.

To start with, I’ll begin by generating data under our model.

Remember, this is an evolutionary scenario in which a discrete trait (evolving via a simple Mk process of any flavor) influences the rate of evolution of a separate continuous character, evolving by variable-rate Brownian motion.

To create my discrete character history, I can use phytools::sim.history. Then, later, to simulate my continuous trait, I’ll use phytools::sim.rates.

Remember that when we go to fit our model, however, we only know the discrete character states at the tips of the tree – not the full history we used for simulation!

## load phytools
library(phytools)
## simulate a tree
tree<-pbtree(n=250,scale=10)
## discrete character transition matrix
Q<-0.1*matrix(c(
  -2,1,1,
  1,-2,1,
  1,1,-2),3,3,
  dimnames=list(c("a","b","c"),
    c("a","b","c")))
Q

##      a    b    c
## a -0.2  0.1  0.1
## b  0.1 -0.2  0.1
## c  0.1  0.1 -0.2

From the value of Q that I intend to use for simulation, the reader should be able to see that my character’s evolution is unordered and equal-rate. I could’ve simulated and fit any arbitrary flavor of Mk model, but for computational reasons I won’t get into (that, surprisingly, don’t have much to do with the number of parameters to be estimated) asymmetric models like "ARD" are much slower to fit to data than are symmetric models.

## simulate discrete character history
tt<-sim.history(tree,Q,anc="a")

## Done simulation(s).

tt

## 
## Phylogenetic tree with 250 tips and 249 internal nodes.
## 
## Tip labels:
## 	t16, t21, t123, t124, t38, t49, ...
## 
## The tree includes a mapped, 3-state discrete character
## with states:
## 	a, b, c
## 
## Rooted; includes branch lengths.

Next, I’ll set the three values of \(\sigma^2\) that correspond to the three levels of my discrete character. Here, I’ll set these at \(\sigma_{a}^2 = 0.5\), \(\sigma_{b}^2 = 5\), and \(\sigma_{c}^2 = 20\). (No particular reason – except that I’m intentionally simulating a large difference in rate!)

## set values of sigma-squared
sig2<-setNames(c(0.5,5,20),c("a","b","c"))
sig2

##    a    b    c 
##  0.5  5.0 20.0

To simulate rate-varying Brownian motion trait evolution using phytools we can use the function sim.rates as follows.

## simulate rate-varying continuous trait
x<-sim.rates(tt,sig2)

To have a visual of what we’ve done, here’s a tree (on the left) showing the true history of my discrete character (on the left) and a reconstruction of the continuous trait projected into trait space on the right.

layout(matrix(c(1,2),1,2),widths=c(0.45,0.55))
cols<-setNames(hcl.colors(n=3),c("a","b","c"))
plot(tt,cols,ftype="off",mar=c(4.1,1.1,1.1,1.1))
par(mar=c(4.1,5.1,1.1,1.1))
phenogram(tt,x,colors=cols,ftype="off",las=1,
  cex.axis=0.8)
legend("bottomleft",c(
  expression(paste(sigma^2,"[a] = 0.5")),
  expression(paste(sigma^2,"[b] = 5.0")),
  expression(paste(sigma^2,"[c] = 20"))),
  col=cols,lwd=3,bty="n",cex=0.8)

Before we analyze our data, we need to peel off just the tip values of our discrete trait. Remember, that’s all we’re giving fitmultiBM (along with our continuous trait in x) to fit our model!

y<-as.factor(getStates(tt,"tips"))

Finally, let’s fit our discrete character dependent rate-varying continuous character evolution model using fitmultiBM as follows.

sd_fit<-fitmultiBM(tree,x,y,model="ER",
  parallel=TRUE,levs=100,root="nuisance")

sd_fit

## Object of class "fitmultiBM" based on
##     a discretization with k = 100 levels.
## 
## Fitted multi-rate BM model parameters:
##  levels: [ a, b, c ]
##   sigsq: [ 0.426, 4.9639, 16.6551 ]
##      x0: -0.0172 
## 
## Estimated Q matrix:
##             a           b           c
## a -0.18884059  0.09442029  0.09442029
## b  0.09442029 -0.18884059  0.09442029
## c  0.09442029  0.09442029 -0.18884059
## 
## Log-likelihood: -812.2065 
## 
## R thinks it has found the ML solution.

What’s very cool to see from this fitted model is just how close our parameter estimates for both the \(\sigma^2\) Brownian rates of continuous trait evolution, and the value of Q, are to their generating values.

But what I posit may be even cooler still is to find out what happens when we fit a hidden-rates model in which we postulate three unobserved rate categories of our continuous character, with untracked evolution between the three levels.

This can be done by simply dropping the discrete character out of our model and adding the argument ncat (for number of hidden rate categories).

Since we know the true number of rate categories is three, we can set ncat=3. (In an empirical case we might explore multiple values of ncat and compare them.)

hrm_fit<-fitmultiBM(tree,x,ncat=3,
  model.hrm="ER",parallel=TRUE,levs=100,
  root="nuisance")

hrm_fit

## Object of class "fitmultiBM" based on
##     a discretization with k = 100 levels.
## 
## Fitted multi-rate BM model parameters:
##  levels: [ 0, 0*, 0** ]
##   sigsq: [ 16.6899, 6.7361, 0.5816 ]
##      x0: -0.2012 
## 
## Estimated Q matrix:
##               0          0*         0**
## 0   -0.14282147  0.07141073  0.07141073
## 0*   0.07141073 -0.14282147  0.07141073
## 0**  0.07141073  0.07141073 -0.14282147
## 
## Log-likelihood: -632.116 
## 
## R thinks it has found the ML solution.

This is very cool.

As I’ve found in prior posts, we see parameter estimates for the \(\sigma^2\) rates that quite closely match the generating heterogeneity. Interestingly, we also find transition rates for our unobserved discrete trait that reasonably closely approximate the true underlying transition rates in our Brownian rate categories. Neat!

Lastly, let’s do ancestral state estimation and see how closely our estimated node and tip rate categories match their generating values! (Take note that this code only works because I know my tip data are in the order of the tip labels of the tree. Otherwise, we should sort them!)

anc_hrm<-ancr(hrm_fit,tips=TRUE)

par(mfrow=c(1,2))
nn<-getStates(tt,"nodes")
plotTree(tree,ftype="off",lwd=1,xlim=c(-2,10.25),
  mar=c(0.1,1.1,2.1,0.1))
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
cols<-setNames(hcl.colors(n=3),c("a","b","c"))
points(pp$xx[1:tree$Nnode+Ntip(tree)],
  pp$yy[1:tree$Nnode+Ntip(tree)],pch=16,
  col=cols[nn],cex=1)
for(i in 1:Ntip(tree)){
  polygon(x=pp$xx[i]+c(0,0.5,0.5,0)+0.1,
    y=pp$yy[i]+c(-0.5,-0.5,0.5,0.5),
    border="transparent",col=cols[y[i]])
}
labs<-c(
  bquote(paste(sigma[a]^2,' = ',.(round(sig2[1],2)))),
  bquote(paste(sigma[b]^2,' = ',.(round(sig2[2],2)))),
  bquote(paste(sigma[c]^2,' = ',.(round(sig2[3],2)))))
legend("topleft",legend=labs,pch=16,col=cols,bty="n",
  cex=0.8,pt.cex=1.2,y.intersp=1.5)
mtext("a) true rate states",adj=0,line=1)
plotTree(tree,ftype="off",lwd=1,xlim=c(10.25,-2),
  mar=c(0.1,0.1,2.1,1.1))
cols<-setNames(rainbow(n=3),c("0","0*","0**"))
par(fg="transparent")
nodelabels(pie=anc_hrm$ace$discrete[251:499,],piecol=cols,
  cex=0.5)
par(fg="black")
labs<-c(
  bquote(paste(sigma^2,' = ',.(round(hrm_fit$sigsq[1],2)))),
  bquote(paste(sigma^2,' = ',.(round(hrm_fit$sigsq[2],2)))),
  bquote(paste(sigma^2,' = ',.(round(hrm_fit$sigsq[3],2)))))
legend("topright",legend=labs,pch=16,col=cols,bty="n",
  cex=0.8,pt.cex=1.2,y.intersp=1.5)
for(i in 1:Ntip(tree)){
  tmp<-c(0,cumsum(anc_hrm$ace$discrete[i,]))
  for(j in 1:ncol(anc_hrm$ace$discrete)){
    polygon(x=pp$xx[i]+0.5*c(tmp[j],tmp[j+1],
      tmp[j+1],tmp[j])+0.1,
      y=pp$yy[i]+c(-0.5,-0.5,0.5,0.5),
      border="transparent",col=cols[j])
  }
}
mtext("b) estimated hidden rates",adj=0,
  line=1)

OK. That’s not bad.