Stochastic character mapping under the threshold model using fitThresh in phytools

Unfortunately, I’m not usually so responsive to email inquiries these days (sorry!), but… just yesterday a phytools user contacted me with the subject line “getting a simmap style tree from a threshold model result” and the following email text:

“I have a weird question regarding your threshold model functions. I have a continuous trait I am making discrete bins out of to test whether multivariate rates of shape evolution vary with the levels of this (now discrete) trait using mvMORPH. fitThresh makes a lot of sense for reconstructing this character, and results in ancestral values that look different/more sensible than make.simmap (based on a reconstruction of the underlying continuous character). But I need the threshold model results in some sort of format that will work with mvBM (e.g. "simmap")…. Is there a not too complicated way to convert the ‘full’ Q matrix into transitions between the original discrete character, or do you have a better recommendation for making fitThresh play nice with the trait fitting functions in phytools/mvMORPH/OUwie/etc? If there is no good answer right now no worries, just couldn’t sus it out on my own.”

Indeed, this is something that I’d thought of already – that is, that the discrete diffusion / finite-space approximation that we use in fitThresh and other new methods of phytools would enable new sorts of stochastic character mapping, including of the threshold model and, indeed, of actual continuous traits.

Before I go on to show how to do this, I will note that this is an approximation of the corresponding stochastic character map for the threshold process, only inasmuch as technically the threshold is crossed \(\infty\) times by the underlying fractal BM process when it moves from one side of the threshold to the other. (If you didn’t already know this, just trust me. It’s true.)

So, how do we do this?

Well, so far I have not yet automated this in software, but it’s pretty straightforward if we go through it step by step.

I’m going to start with the simplest example, which is the threshold model with two states, let’s call them "a" and "b".

## load phytools
library(phytools)

First, let’s simulate some data under this process.

## simulate tree
phy<-pbtree(n=80,scale=1)
phy

## 
## Phylogenetic tree with 80 tips and 79 internal nodes.
## 
## Tip labels:
##   t28, t31, t32, t23, t24, t19, ...
## 
## Rooted; includes branch length(s).

## simulate liabilities
liability<-fastBM(phy,a=0.5)
head(liability,10)

##        t28        t31        t32        t23        t24        t19        t25        t26 
##  0.9810162  1.1667724  0.9394853  0.5019343  0.2609373 -0.1475927 -0.4670823  0.2442964 
##        t27        t33 
## -0.7697432 -1.0546414

## "threshold" liabilities
thresh<-function(x) if(x<0) "a" else "b"
x<-sapply(liability,thresh)
head(x,20)

## t28 t31 t32 t23 t24 t19 t25 t26 t27 t33 t35 t37 t38 t79 t80  t9 t29 t30 t46 t47 
## "b" "b" "b" "b" "b" "a" "a" "b" "a" "a" "a" "a" "a" "a" "a" "b" "a" "a" "a" "a"

Now let’s proceed to fit the threshold model to these data using fitThresh as follows.

thresh_fit<-fitThresh(phy,x)
thresh_fit

## Object of class "fitThresh".
## 
##     Set value of sigsq (of the liability) = 1.0
## 
## 	  Set or estimated threshold(s) =  [ 0 ]*
## 
##     Log-likelihood: -33.901524 
## 
## (*lowermost threshold is fixed)

(Totally just for fun, let’s fit & compare a symmetric Mk model to see if they’re distinguishable based on these data.)

mk_model<-fitMk(phy,x,model="ER")
anova(mk_model,thresh_fit)

##               log(L) d.f.      AIC    weight
## mk_model   -36.04421    1 74.08842 0.1050167
## thresh_fit -33.90152    1 69.80305 0.8949833

Now, to undertake stochastic mapping using our threshold model, we need to pull out the implicit Mk model that’s hidden within our "fitThresh" object.

mkm<-thresh_fit$mk_fit

(I’m not going to print it because the Q matrix is \(200 \times 200\), but readers following along should feel free to go for it.)

Next, because this hidden object is not quite the same as a standard "fitMk" object, let’s add the element root.prior as follows.

mkm$root.prior<-"fitzjohn"

Now, believe it or not, we’re totally ready for stochastic mapping.

When we run the following code we’re probably going to see the message Warning in rstate(p/sum(p)) : Some probabilities (slightly?) < 0. Setting p < 0 to zero. This is just because our matrix is so big & some of the elements are close to zero, so we can safely ignore it.

I’m only going to do one stochastic character map here. Typically we should do many of these, of course. This would be accomplished by modifying nsim and then iterating all subsequent steps over each stochastic map tree.

smp<-simmap(mkm,nsim=1,pi="mle")

## Warning in rstate(p/sum(p)): Some probabilities (slightly?) < 0. Setting p < 0 to zero.
## Warning in rstate(p/sum(p)): Some probabilities (slightly?) < 0. Setting p < 0 to zero.
## Warning in rstate(p/sum(p)): Some probabilities (slightly?) < 0. Setting p < 0 to zero.
## ...
## Warning in rstate(p/sum(p)): Some probabilities (slightly?) < 0. Setting p < 0 to zero.

Once again, this is a big object because it contains all 200 of the finely discretized levels of our diffusion approximation – effectively, a stochastic map of liabilities.

To turn this into a stochastic character map of our original discrete trait we need to simply merge liabilities on each side of the single threshold as follows.

ss<-colnames(mkm$data)
merged.smp<-mergeMappedStates(smp,ss[which(as.numeric(ss)<0)],"a")
merged.smp<-mergeMappedStates(merged.smp,ss[which(as.numeric(ss)>0)],"b")
merged.smp

## 
## Phylogenetic tree with 80 tips and 79 internal nodes.
## 
## Tip labels:
## 	t28, t31, t32, t23, t24, t19, ...
## 
## The tree includes a mapped, 2-state discrete character
## with states:
## 	a, b
## 
## Rooted; includes branch lengths.

Great. Let’s plot it.

plot(merged.smp,direction="upwards",ftype="off",
  lwd=3,colors=setNames(hcl.colors(n=2),c("a","b")))
legend("bottomleft",c("a","b"),col=hcl.colors(n=2),
  lwd=2,bty="n")

We can see that this is clearly capturing a form of the threshold model history of our threshold character by seeing just how different it looks compared to running the same analysis with a standard Mk model!

mk.smp<-simmap(mk_model,nsim=1)
plot(mk.smp,direction="upwards",ftype="off",
  lwd=3,colors=setNames(hcl.colors(n=2),c("a","b")))
legend("bottomleft",c("a","b"),col=hcl.colors(n=2),
  lwd=2,bty="n")

We see, for instance, that threshold crossing usually involves many switches back & forth under the threshold model, and none in the Mk model, which makes perfect sense!

Things get only a little more complicated when I want to fit the multi-state model.

Let’s see if I can run through that.

First, I’ll reuse my original tree & liabilities to simulate some tip data as follows:

y<-threshState(liability,setNames(c(0,1,2),letters[1:3]))
head(y)

## t28 t31 t32 t23 t24 t19 
## "b" "c" "b" "b" "b" "a"

Next, fit the multi-state threshold model as follows.

ms_thresh<-fitThresh(phy,y,sequence=letters[1:3])
ms_thresh

## Object of class "fitThresh".
## 
##     Set value of sigsq (of the liability) = 1.0
## 
## 	  Set or estimated threshold(s) =  [ -0.7568, 0.341537 ]*
## 
##     Log-likelihood: -51.847661 
## 
## (*lowermost threshold is fixed)

Note that the lower estimated threshold is not a separately estimable parameter (otherwise the model would be non-identifiable). It’s normally set to zero, but fitThresh does something different that I’ll probably fix in the future. (It centers the whole liability distribution on zero instead. Since the liabilities are unitless and scaleless this doesn’t matter at all for the model fit or relative positions of the thresholds, but our results are a bit more annoying to interpret.)

Pull out our hidden Mk object & do stochastic mapping, as before.

mkn<-ms_thresh$mk_fit
mkn$root.prior<-"fitzjohn"
smp<-simmap(mkn,nsim=1)

## Warning in rstate(p/sum(p)): Some probabilities (slightly?) < 0. Setting p < 0 to zero.
## Warning in rstate(p/sum(p)): Some probabilities (slightly?) < 0. Setting p < 0 to zero.
## Warning in rstate(p/sum(p)): Some probabilities (slightly?) < 0. Setting p < 0 to zero.
## ...
## Warning in rstate(p/sum(p)): Some probabilities (slightly?) < 0. Setting p < 0 to zero.

Now, repeat the mergeMappedStates step, but this time using the estimated thresholds.

ss<-mkn$states
merged.smp.2<-mergeMappedStates(smp,
  ss[which(as.numeric(ss)<ms_thresh$threshold[1])],"a")
merged.smp.2<-mergeMappedStates(merged.smp.2,
  ss[which(as.numeric(ss)>=ms_thresh$threshold[1]&
      as.numeric(ss)<ms_thresh$threshold[2])],"b")
merged.smp.2<-mergeMappedStates(merged.smp.2,
  ss[which(as.numeric(ss)>=ms_thresh$threshold[2])],"c")
merged.smp.2

## 
## Phylogenetic tree with 80 tips and 79 internal nodes.
## 
## Tip labels:
## 	t28, t31, t32, t23, t24, t19, ...
## 
## The tree includes a mapped, 3-state discrete character
## with states:
## 	a, b, c
## 
## Rooted; includes branch lengths.

Let’s plot it!

plot(merged.smp.2,direction="upwards",ftype="off",
  lwd=3,colors=setNames(hcl.colors(n=3),c("a","b","c")))
legend("bottomleft",c("a","b","c"),col=hcl.colors(n=3),
  lwd=3,bty="n")

Wow. That’s beautiful!

Very cool.