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.

Comparing a discrete character dependent multi-regime OU joint model to OUwie

In some recent posts to this blog (e.g., 1, 2, 3, 4) I have described a new discrete character dependent multi-optimum model in phytools that uses the finite-space or discrete diffusion approximation of Boucher & Démery (2016; also see our in review bioRxiv pre-print).

Even though they implement different models, I thought it might be interesting to compare this new method, fitmultiOU, to a fixed-regime multi-optimum OU model fit using the popular OUwie package by Jeremy Beaulieu and Brian O’Meara.

Since our new method integrates over uncertainty in the regime history by jointly optimizing an Mk transition process for the regimes along with the multi-optimum stochastic process of our continuous trait, my thesis is that the parameter estimates we obtain should be highly similar between the two models if our discrete trait changes infrequently.

To explore this, we can start by simulating a tree, a discrete character history, and some data using phytools as follows.

## load phytools
library(phytools)

## simulate a tree
N<-100 ## number of taxa
phy<-pbtree(n=N,scale=10)
phy

## 
## Phylogenetic tree with 100 tips and 99 internal nodes.
## 
## Tip labels:
##   t2, t3, t89, t90, t35, t53, ...
## 
## Rooted; includes branch length(s).

## set the transition matrix of our
## discrete trait
q<-0.05
k<-2
Q<-matrix(q,k,k,
  dimnames=list(letters[1:k],letters[1:k]))
diag(Q)<-0
diag(Q)<--rowSums(Q)
Q

##       a     b
## a -0.05  0.05
## b  0.05 -0.05

## simulate trait history
sim_tree<-sim.history(phy,Q,anc=sample(letters[1:k],1))

## Done simulation(s).

cols<-setNames(hcl.colors(n=k),letters[1:k])
plot(sim_tree,cols,ftype="off",lwd=2,
  direction="upwards")
par(lend=1)
legend("bottomleft",letters[1:k],lwd=4,
  col=hcl.colors(n=k),cex=0.7,bty="n")

This is great because our tree has very few changes in the character. My hypothesis is that this will make parameter estimates of the OU process very similar between a fixed regime model, as in OUwie, and our new fitmultiOU method.

Now let’s set the generating conditions of our discrete character dependentn multi-\(\theta\) Ornstein-Uhlenbeck process.

The parameters \(\alpha\) and \(\sigma^2\) will be the same across all regimes (though we still have to specify \(k = 2\) of each for our generator), while \(\theta\) will vary according to the state of our discrete trait.

## set alpha
alpha<-setNames(rep(0.3,k),letters[1:k])
alpha

##   a   b 
## 0.3 0.3

## set sigma-squared
sig2<-setNames(rep(0.1,k),letters[1:k])
sig2

##   a   b 
## 0.1 0.1

# set theta
theta<-setNames(c(-0.5,2),letters[1:k])
theta

##    a    b 
## -0.5  2.0

At this point we can simulate our continuous trait using phytools::multiOU.

## simulate continuous trait
x<-multiOU(sim_tree,alpha,sig2,theta,a0=0)
head(x)

##         t2         t3        t89        t90        t35        t53 
## -0.3416784 -0.6347880 -0.3128115 -0.5001694 -0.7747887 -0.5183583

Our discrete character is already simulated, but we need to pull of the tip values from the sim_tree "simmap" object to use them in our analysis.

## pull off discrete trait
y<-as.factor(getStates(sim_tree,"tips"))
head(y)

##  t2  t3 t89 t90 t35 t53 
##   a   a   a   a   a   a 
## Levels: a b

For completeness, let’s start by fitting our null model using fitmultiOU.

We can then confirm that this fitted model matches (to a reasonable degree – they will only converge exactly as levs goes towards \(\infty\)) what we’d obtain using geiger::fitContinuous and phytools::fitMk under the same model assumptions. I'm going to set levs = 200 for this analysis, but this actually takes a very long time to run (much more than twice as long as levs = 100).

## fit null model
fit_null<-fitmultiOU(phy,x,y,model="ER",levs=200,
  parallel=TRUE,ncores=10,root="mle",trace=1,
  null_model=TRUE)

## iter	theta	alpha	sigsq	q[1]	log(L)
## 0	1.9456	0.0721	1.1329	0.0914	-175.6004 
## 100	0.6244	0.0116	0.0897	0.0271	-98.3007 
## 200	0.6123	0.0053	0.0892	0.0283	-98.2222 
## 271	0.6338	0.0004	0.0860	0.0292	-98.1838 
## Done optimizing.

fit_null

## Object of class "fitmultiOU" based on
##     a discretization with k = 200 levels.
## 
## Fitted multi-theta OU model parameters:
##  levels: [ a, b ]
##   theta: [ 0.6338 ]
##   alpha: 4e-04 
##   sigsq: 0.086 
## 
## Estimated Q matrix:
##             a           b
## a -0.02921177  0.02921177
## b  0.02921177 -0.02921177
## 
## Log-likelihood: -98.1838 
## 
## R thinks optimization may not have converged.

## fit single regime OU model using fitContinuous
ou_fit<-geiger::fitContinuous(phy,x,model="OU")
ou_fit

Now, again, let's compare to geiger::fitContinuous and phytools::fitMk.

## GEIGER-fitted comparative model of continuous data
##  fitted 'OU' model parameters:
## 	alpha = 0.000000
## 	sigsq = 0.085106
## 	z0 = 0.625022
## 
##  model summary:
## 	log-likelihood = -68.250476
## 	AIC = 142.500953
## 	AICc = 142.750953
## 	free parameters = 3
## 
## Convergence diagnostics:
## 	optimization iterations = 100
## 	failed iterations = 0
## 	number of iterations with same best fit = 50
## 	frequency of best fit = 0.500
## 
##  object summary:
## 	'lik' -- likelihood function
## 	'bnd' -- bounds for likelihood search
## 	'res' -- optimization iteration summary
## 	'opt' -- maximum likelihood parameter estimates

## fit Mk model using fitMk
mk_fit<-fitMk(phy,y,model="ER",pi="equal")
mk_fit

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           a         b
## a -0.028064  0.028064
## b  0.028064 -0.028064
## 
## Fitted (or set) value of pi:
##   a   b 
## 0.5 0.5 
## due to treating the root prior as (a) flat.
## 
## Log-likelihood: -29.168076 
## 
## Optimization method used was "nlminb"
## 
## R thinks it has found the ML solution.

You can compare the model parameter estimates, but let’s also assure ourselves that the likelihood seems to be converging on the same value as follows.

## compute null log(L) from fitContinuous & fitMk results
null_logL<-logLik(ou_fit)+logLik(mk_fit)
attr(null_logL,"df")<-4 ## fix d.f.
null_logL

## [1] -97.41855
## attr(,"df")
## [1] 4

## compare to fitmultiOU
logLik(fit_null)

## [1] -98.18378
## attr(,"df")
## [1] 4

This is pretty close. Again, we would expect these two values to get even closer for higher levs, but (as currently implemented) this already takes a really long time to run!

OK. Now let’s bring OUwie into the picture. We can start by loading the package, which I recently updated from CRAN.

## load OUwie
library(OUwie)

Now for OUwie we need to put our data in a special data frame format as follows.

## compile our data for OUwie
ouwie.data<-data.frame(
  Genus_species=names(x),
  Reg=y,
  X=x)
head(ouwie.data)

##     Genus_species Reg          X
## t2             t2   a -0.3416784
## t3             t3   a -0.6347880
## t89           t89   a -0.3128115
## t90           t90   a -0.5001694
## t35           t35   a -0.7747887
## t53           t53   a -0.5183583

Why don’t we start by simply re-fitting our null OU model in OUwie?

This should give us a result that quite closely matches what we obtained using geiger::fitContinuous.

I’ll still give it our known discrete character data & history, but I set model = "OU1" to specify that I want a single \(\theta\) model only.

## fit OUwie null model
fitOU.smp<-OUwie(sim_tree,ouwie.data,model="OU1",
  simmap.tree=TRUE,root.station=FALSE)

## Warning: An algorithm was not specified. Defaulting to computing the determinant and inversion of the vcv.

## Initializing... 
## Finished. Begin thorough search... 
## Finished. Summarizing results.

fitOU.smp

## 
## Fit
##        lnL     AIC    AICc      BIC model ntax
##  -68.25048 142.501 142.751 150.3165   OU1  100
## 
## 
## Rates
##        alpha     sigma.sq 
## 1.524044e-08 8.510672e-02 
## 
## Optima
##                  1
## estimate 0.6250224
## se       0.3348563
## 
## 
## Half life (another way of reporting alpha)
##    alpha 
## 45480781 
## 
## Arrived at a reliable solution

Hopefully, we see that this fitted model pretty closely matches what we got using fitContinuous earlier.

Next, I’m going to go ahead & fit our discrete character dependent multi- \(\theta\) OU model using phytools::fitmultiOU.

This is the model that I’ve been blogging about recently - but, to remind the reader, this is a joint discrete & continuous trait model, not the fixed regime model of OUwie – but I’m hypothesizing that our continuous trait model parameter estimates should pretty closely match what we’d get from OUwie using the true history (or, for that matter, a stochastic character history), just because our discrete character changes so infrequently on the tree.

To start with, I'm going to try to get reasonable starting values for my fitmultiOU parameters, as follows.

## identify sensible starting parameter values
init<-setNames(
  c(
    mean(x[y==levels(y)[1]]),
    mean(x[y==levels(y)[2]]),
    log(2)/max(nodeHeights(phy)),
    var(x)/max(nodeHeights(phy)),
    fitMk(phy,y,model="ER")$rates),
  c("theta[a]","theta[b]",
    "alpha","sigsq","q[1]"))
init

##    theta[a]    theta[b]       alpha       sigsq        q[1] 
## -0.27841916  1.61842864  0.06931472  0.12828760  0.02806402

Then I can go ahead and fit the joint model.

## fit discrete trait dependent model
fit_mou<-fitmultiOU(phy,x,y,model="ER",levs=200,
  parallel=TRUE,ncores=10,root="mle",trace=1,
  maxit=2000,init=init)

## iter	the[a]	the[b]	alpha	sigsq	q[1]	log(L)
## 0	-0.2784	1.6184	0.0693	0.1283	0.0281	-92.1417 
## 100	-0.4558	1.8575	0.3260	0.1157	0.0268	-72.4242 
## 200	-0.4520	1.8537	0.3254	0.1143	0.0267	-72.4172 
## 300	-0.4531	1.8535	0.3252	0.1145	0.0269	-72.4170 
## 400	-0.4520	1.8523	0.3246	0.1142	0.0271	-72.4163 
## 404	-0.4519	1.8523	0.3246	0.1141	0.0271	-72.4163 
## Done optimizing.

fit_mou

## Object of class "fitmultiOU" based on
##     a discretization with k = 200 levels.
## 
## Fitted multi-theta OU model parameters:
##  levels: [ a, b ]
##   theta: [ -0.4519, 1.8523 ]
##   alpha: 0.3246 
##   sigsq: 0.1141 
## 
## Estimated Q matrix:
##             a           b
## a -0.02713072  0.02713072
## b  0.02713072 -0.02713072
## 
## Log-likelihood: -72.4163 
## 
## R thinks it has found the ML solution.

Finally, we can fit our fixed regime model using the OUwie package. Here, don’t forget that though I will be using the true discrete character history, this is almost never known in practice! (Indeed, if we genuinely knew the discrete character history, I would always recommend fitting a fixed regime model, not a joint model.)

## fit multi-regime OU model using OUwie
fitOUM.smp<-OUwie(sim_tree,ouwie.data,model="OUM",
  simmap.tree=TRUE,root.station=FALSE)

## Warning: An algorithm was not specified. Defaulting to computing the determinant and inversion of the vcv.

## Initializing... 
## Finished. Begin thorough search... 
## Finished. Summarizing results.

fitOUM.smp

## 
## Fit
##        lnL      AIC     AICc      BIC model ntax
##  -44.47594 96.95187 97.37292 107.3726   OUM  100
## 
## 
## Rates
##                  b         a
## alpha    0.3007500 0.3007500
## sigma.sq 0.1198701 0.1198701
## 
## Optima
##                   b          a
## estimate 1.75307314 -0.5184579
## se       0.09854245  0.1132707
## 
## 
## Half life (another way of reporting alpha)
##        b        a 
## 2.304729 2.304729 
## 
## Arrived at a reliable solution

Cool.

Now, if we look closely at our results, we should see that our parameter estimates do in fact match up fairly well, keeping in mind, of course, that these are not the same models – the OUwie model is based on fixed regimes, while our fitmultiOU function implements a joint discrete trait & continuous character evolutionary model.

I would expect this similarity to hold whenever our discrete character history is pretty unambiguous, usually because the character changes quite infrequently on the tree, but to diminish as the rate or number of changes of our discrete trait increases.

That's all folks!

Visualizing a multi-regime OU process on the tree using phytools

In a recent tweet about phytools’ new discrete character dependent multi-optimum OU model, I included an illustration of this process that was not derived from the original blog post.

More on a discrete character dependent multi-optimum OU model in #phytools: https://t.co/LpLHskJoo6. pic.twitter.com/JrmAlmQpbh

— Liam Revell (@phytools_liam) July 13, 2026

Just in case it might become useful later, I decided to post the code for that cool illustration here. In contrast to my prior use, here I decided to set different values for \(\alpha\) & \(\sigma^2\), as well as for \(\theta\), between the three regimes.

## load packages
library(phytools)

In this first chunk I’m going to generate the data and objects I need for the plot.

## simulate a tree
tree<-pbtree(n=100,scale=10)
## add singleton (unbranching) nodes
tt<-map.to.singleton(make.era.map(tree,
  limits=seq(0,10,length.out=1001)))
## set the transition matrix for the discrete
## trait
q<-0.2
Q<-matrix(c(
  -2*q,q,q,
  q,-2*q,q,
  q,q,-2*q),3,3,
  dimnames=list(letters[1:3],letters[1:3]))
## generate a character history of the discrete
## trait
s.tt<-sim.history(tt,Q,anc="a")

## Done simulation(s).

## set the parameters of the multi-regime OU
## process
theta<-setNames(c(-0.5,0.9,2.2),letters[1:3])
sigsq<-setNames(c(0.12,0.05,0.08),letters[1:3])
alpha<-setNames(c(0.7,0.4,0.8),letters[1:3])
## generate data for tips (and nodes) under the
## multi-regime OU process
x<-multiOU(s.tt,alpha=alpha,sig2=sigsq,
  theta=theta,a0=theta["a"],internal=TRUE)

Now for our plot.

## set colors
cols<-setNames(hcl.colors(n=3),letters[1:3])
## set plot layout
layout(matrix(c(1,2),2,1),heights=c(0.4,0.6))
## plot discrete character history on the tree
plot(s.tt,cols,ftype="off",mar=c(1.1,4.1,2.1,2.1),
  xlim=c(0,11))
mtext("a)",adj=0,line=0)
## set margins for second subplot
par(mar=c(5.1,4.1,2.1,2.1))
## plot continuous character history
phenogram(s.tt,x,ftype="off",
  colors=make.transparent(cols,0.5),
  xlim=c(0,11),cex.axis=0.8,las=1)
mtext("b)",adj=0,line=1.5)
## add the predicted equilibrium distribution
## at the tips
for(i in 1:length(theta)){
  stat_theta<-dnorm(seq(min(x),max(x),
    length.out=200),mean=theta[i],
    sd=sqrt(sigsq[i]/(2*alpha[i])))
  stat.norm<-stat_theta/max(stat_theta)
  polygon(x=c(0,stat.norm,0)+10,
    y=c(min(x),seq(min(x),max(x),
      length.out=200),max(x)),border=FALSE,
    col=make.transparent(cols[i],0.5))
}

Nothing to it!

More on a discrete character dependent multi-optimum OU model in phytools

Good morning blog readers.

Recently I’ve posted about a new discrete character dependent multi- \(\theta\) (i.e., multi-optimum) OU model in phytools (e.g., 1, 2, 3).

Since this is very new it should only be used with the utmost caution. Nonetheless, I thought I’d post up a quick demo of how it works, and how to do a model comparison to a simpler model of joint discrete & continuous trait evolution but without dependence. (We also might consider an alternative null model with hidden characters, but this will have to be covered in a future post!)

This will only work on recent (at the time of writing) versions of phytools, so we can start by loading the package & checking which version we have.

library(phytools)

## Loading required package: ape

## Loading required package: maps

packageVersion("phytools")

## [1] '2.6.3'

To fit this model I’m going to need some data. With none readily at hand, I’m going to use phytools to simulate some.

We can start with a tree:

N<-200 ## number of taxa
phy<-pbtree(n=N,scale=10)
phy

## 
## Phylogenetic tree with 200 tips and 199 internal nodes.
## 
## Tip labels:
##   t6, t7, t44, t57, t136, t177, ...
## 
## Rooted; includes branch length(s).

Next, we want a generating discrete character history for our multi-regime OU process. Note that though we use this regime history for simulation, in an empirical case it would’ve been unknown, so shall naturally be set aside when we move forward to estimation.

## set the transition matrix of our
## discrete trait
q<-0.2
Q<-matrix(c(
  -2*q,q,q,
  q,-2*q,q,
  q,q,-2*q),3,3,
  dimnames=list(letters[1:3],letters[1:3]))
Q

##      a    b    c
## a -0.4  0.2  0.2
## b  0.2 -0.4  0.2
## c  0.2  0.2 -0.4

k<-nrow(Q) ## trait levels
k

## [1] 3

sim_tree<-sim.history(phy,Q,anc="a")

## Done simulation(s).

sim_tree

## 
## Phylogenetic tree with 200 tips and 199 internal nodes.
## 
## Tip labels:
## 	t6, t7, t44, t57, t136, t177, ...
## 
## The tree includes a mapped, 3-state discrete character
## with states:
## 	a, b, c
## 
## Rooted; includes branch lengths.

Let’s plot our generating tree as follows.

cols<-setNames(hcl.colors(n=3),letters[1:k])
plot(sim_tree,cols,ftype="off",lwd=1,
  direction="upwards")
par(lend=1)
legend("bottomleft",letters[1:3],lwd=3,
  col=hcl.colors(n=k),cex=0.8,bty="n")

Next, we can set the generating conditions for our continous trait simulation. Our model allows for multiple \(\theta\) by discrete character state, but assumes constant \(\alpha\) and \(\sigma^2\) across the \(k\) levels of our discrete trait, so let’s simulate that.

alpha<-setNames(rep(0.3,k),letters[1:k])
alpha

##   a   b   c 
## 0.3 0.3 0.3

sig2<-setNames(rep(0.1,k),letters[1:k])
sig2

##   a   b   c 
## 0.1 0.1 0.1

theta<-setNames(c(-0.5,1,2),letters[1:k])
theta

##    a    b    c 
## -0.5  1.0  2.0

Now we’re nearly ready to simulate our continuous trait. To do that, I’ll use phytools::multiOU as I have in prior posts.

x<-multiOU(sim_tree,alpha,sig2,theta,a0=0)
head(x)

##        t6        t7       t44       t57      t136      t177 
## 1.9418267 0.5788125 1.2039455 0.7955580 0.6064443 0.9103976

Though we’ve simulated our discrete character history already, for our analysis we’ll use just the tip states, so let’s pull those into a factor vector using phytools::getStates.

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

##   t6   t7  t44  t57 t136 t177 
##    a    b    c    a    c    c 
## Levels: a b c

Awesome. Now let’s first fit our null model using fitmultiOU.

fit_null<-fitmultiOU(phy,x,y,model="ER",levs=100,
  parallel=TRUE,ncores=10,root="mle",trace=1,
  null_model=TRUE)

## iter	theta	alpha	sigsq	q[1]	log(L)
## 0	1.0601	0.2025	0.2682	0.0177	-409.2704 
## 100	0.4702	0.0182	0.0846	0.2622	-287.2005 
## 200	0.5735	0.0150	0.0829	0.2569	-287.0941 
## 279	0.5743	0.0141	0.0829	0.2555	-287.0493 
## Done optimizing.

fit_null

## Object of class "fitmultiOU" based on
##     a discretization with k = 100 levels.
## 
## Fitted multi-theta OU model parameters:
##  levels: [ a, b, c ]
##   theta: [ 0.5743 ]
##   alpha: 0.0141 
##   sigsq: 0.0829 
## 
## Estimated Q matrix:
##            a          b          c
## a -0.5109335  0.2554668  0.2554668
## b  0.2554668 -0.5109335  0.2554668
## c  0.2554668  0.2554668 -0.5109335
## 
## Log-likelihood: -287.0493 
## 
## R thinks it has found the ML solution.

Let’s confirm that our parameter estimates and log-likelihood match what we would’ve obtained using a geiger::fitContinuous and phytools::fitMk. This isn’t hard.

ou_fit<-geiger::fitContinuous(phy,x,model="OU")
ou_fit

## GEIGER-fitted comparative model of continuous data
##  fitted 'OU' model parameters:
## 	alpha = 0.012919
## 	sigsq = 0.081833
## 	z0 = 0.655514
## 
##  model summary:
## 	log-likelihood = -89.743851
## 	AIC = 185.487702
## 	AICc = 185.610151
## 	free parameters = 3
## 
## Convergence diagnostics:
## 	optimization iterations = 100
## 	failed iterations = 0
## 	number of iterations with same best fit = 50
## 	frequency of best fit = 0.500
## 
##  object summary:
## 	'lik' -- likelihood function
## 	'bnd' -- bounds for likelihood search
## 	'res' -- optimization iteration summary
## 	'opt' -- maximum likelihood parameter estimates

mk_fit<-fitMk(phy,y,model="ER",pi="equal")
mk_fit

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           a         b         c
## a -0.511178  0.255589  0.255589
## b  0.255589 -0.511178  0.255589
## c  0.255589  0.255589 -0.511178
## 
## Fitted (or set) value of pi:
##        a        b        c 
## 0.333333 0.333333 0.333333 
## due to treating the root prior as (a) flat.
## 
## Log-likelihood: -195.760975 
## 
## Optimization method used was "nlminb"
## 
## R thinks it has found the ML solution.

null_logL<-logLik(ou_fit)+logLik(mk_fit)
null_logL

## [1] -285.5048
## attr(,"df")
## [1] 3

This should be very close to the values we obtained in fit_null. In fact, the two values are a bit farther apart than I'm comfortable with, but would undoubtedly converge if we were to increase `levs`. We should do this with some caution, though, because it very substantially is going to increase our run time.

Finally, we can fit our discrete character dependent multi- \(\theta\) OU model!

fit_mou<-fitmultiOU(phy,x,y,model="ER",levs=20,
  parallel=TRUE,ncores=10,root="mle",trace=1,
  maxit=2000)

## iter	the[a]	the[b]	the[c]	alpha	sigsq	q[1]	log(L)
## 0	-0.3575	0.3680	1.8215	0.2403	0.0115	0.2066	-309.6850 
## 100	-0.3023	1.0148	1.6594	0.2779	0.0615	0.1947	-284.8126 
## 200	-0.2191	1.1288	1.9057	0.2837	0.0611	0.1861	-284.1910 
## 300	-0.2140	1.1372	1.9328	0.2794	0.0554	0.1914	-284.0423 
## 400	-0.2355	1.1808	1.9946	0.2713	0.0499	0.1976	-283.9893 
## 500	-0.2229	1.1901	2.0192	0.2716	0.0514	0.2043	-283.9538 
## 600	-0.2202	1.2210	1.9799	0.2748	0.0519	0.2054	-283.9183 
## 700	-0.2172	1.2083	1.9723	0.2747	0.0528	0.2056	-283.9137 
## 800	-0.2164	1.2088	1.9738	0.2751	0.0528	0.2053	-283.9135 
## 900	-0.2156	1.2094	1.9724	0.2753	0.0527	0.2052	-283.9134 
## 913	-0.2157	1.2097	1.9725	0.2752	0.0527	0.2053	-283.9134 
## Done optimizing.

Just because I know that optimization of this model is difficult, I'm a bit suspicious we may not have converged on the true MLE. Let’s try again, but with “sensible” starting values for all our different model parameters.

init<-setNames(
  c(
    mean(x[y==levels(y)[1]]),
    mean(x[y==levels(y)[2]]),
    mean(x[y==levels(y)[3]]),
    log(2)/max(nodeHeights(phy)),
    var(x)/max(nodeHeights(phy)),
    fitMk(phy,y,model="ER")$rates),
  c("theta[a]","theta[b]","theta[c]",
    "alpha","sigsq","q[1]"))
init

##   theta[a]   theta[b]   theta[c]      alpha      sigsq       q[1] 
## 0.23239962 0.64856597 0.93721136 0.06931472 0.04357107 0.25558924

fit_mou<-fitmultiOU(phy,x,y,model="ER",levs=100,
  parallel=TRUE,ncores=10,root="mle",trace=1,
  maxit=2000,init=init)

## iter	the[a]	the[b]	the[c]	alpha	sigsq	q[1]	log(L)
## 0	0.2324	0.6486	0.9372	0.0693	0.0436	0.2556	-315.4681 
## 100	-0.5480	0.6211	2.3112	0.1711	0.0816	0.1765	-270.5138 
## 200	-0.9475	0.9025	2.0608	0.1907	0.0725	0.2221	-268.0952 
## 300	-0.8971	0.7990	2.0222	0.1960	0.0707	0.2176	-267.9648 
## 400	-0.8941	0.8285	2.0684	0.1935	0.0701	0.2202	-267.9573 
## 500	-0.8831	0.8318	2.0694	0.1936	0.0707	0.2166	-267.9496 
## 600	-0.8828	0.8724	2.0711	0.1903	0.0708	0.2171	-267.9394 
## 700	-0.8728	0.8737	2.0625	0.1894	0.0710	0.2152	-267.9342 
## 800	-0.8513	0.8355	1.9719	0.2003	0.0724	0.2116	-267.9037 
## 900	-0.8589	0.8295	1.9641	0.2014	0.0720	0.2139	-267.9002 
## 1000	-0.8613	0.8330	1.9623	0.2014	0.0719	0.2149	-267.8995 
## 1100	-0.8608	0.8220	1.9628	0.2017	0.0713	0.2154	-267.8972 
## 1200	-0.8600	0.8276	1.9672	0.2008	0.0712	0.2156	-267.8957 
## 1300	-0.8269	1.0429	2.0894	0.1860	0.0700	0.2134	-267.8096 
## 1400	-0.8193	1.0484	2.1319	0.1836	0.0691	0.2161	-267.7998 
## 1500	-0.8143	1.0503	2.1171	0.1844	0.0696	0.2132	-267.7953 
## 1600	-0.7951	1.0509	2.1326	0.1842	0.0690	0.2107	-267.7829 
## 1700	-0.7939	1.0513	2.1340	0.1851	0.0689	0.2095	-267.7812 
## 1723	-0.7939	1.0513	2.1335	0.1852	0.0689	0.2094	-267.7812 
## Done optimizing.

fit_mou

## Object of class "fitmultiOU" based on
##     a discretization with k = 100 levels.
## 
## Fitted multi-theta OU model parameters:
##  levels: [ a, b, c ]
##   theta: [ -0.7939, 1.0513, 2.1335 ]
##   alpha: 0.1852 
##   sigsq: 0.0689 
## 
## Estimated Q matrix:
##            a          b          c
## a -0.4188483  0.2094242  0.2094242
## b  0.2094242 -0.4188483  0.2094242
## c  0.2094242  0.2094242 -0.4188483
## 
## Log-likelihood: -267.7812 
## 
## R thinks it has found the ML solution.

Here it seems that we've definitely done much better.

Remember, the generating values of \(\theta\) were as follows:

theta

##    a    b    c 
## -0.5  1.0  2.0

Our estimates are remarkably close to the generating conditions of our simulation!

Naturally, we might be interested to know whether our discrete character dependent model better explains our continuous trait data than the independent null model. This comparison is very easy.

anova(fit_null,fit_mou)

##             log(L) d.f.      AIC       weight
## fit_null -287.0493    4 582.0985 3.166496e-08
## fit_mou  -267.7812    6 547.5624 1.000000e+00

This tells us that basically all of the weight of evidence falls on our discrete character dependent multi-regime OU model compared to the null model.

That’s all there is to it!

Preliminary, but very cool.

Null model for discrete character dependent multi-θ Ornstein-Uhlenbeck model in phytools

A few weeks ago I posted about a discrete character dependent multi-\(\theta\) Ornstein-Uhlenbeck model using the discrete diffusion approximation of Boucher & Démery (2016) and our bioRxiv pre-print. (As in my prior post, I recommend that those interested in the technical details of this general approach check out our pre-print.)

Since that time, I’ve been corresponding with a colleague who has been trying to use the prototype fitmultiOU function. Consequently, I decided to write today’s post demonstrating how we fit the “null model” of joint continuous trait OU & discrete character Mk evolution but without discrete character dependence using phytools. I’m going to illustrate (I hope) that we get the same (to a reasonable extent of numerical precision) parameter estimates & likelihood as we would obtain from geiger::fitContinuous under a single \(\theta\) OU model and phytools::fitMk. Note that in a prior GitHub update of phytools I was inadvertently trying to separately estimate \(\theta\) and \(x_0\), the root state – but these are not identifiable under a single regime OU model, so I have set them equal to one another in phytools \(\geq\) 2.6-3.

OK. Let’s get started.

## load phytools & check package version
library(phytools)

## should be phytools >= 2.6-3
packageVersion("phytools")

## [1] '2.6.3'

For this demo I’m going to simulate a pretty small tree. This is because I’m going to set levs = 200 for estimation and this will take a while. In practice, we probably need larger trees to fit a multi-regime OU model.

## simulate tree
N<-60
phy<-pbtree(n=N,scale=10)
phy

## 
## Phylogenetic tree with 60 tips and 59 internal nodes.
## 
## Tip labels:
##   t9, t10, t53, t54, t19, t49, ...
## 
## Rooted; includes branch length(s).

Next, let’s specify a generating transition matrix of our discrete trait, Q.

## set generating Q matrix for discrete character
q<-0.2
Q<-matrix(c(
  -q,q,
  q,-q),2,2,
  dimnames=list(letters[1:2],letters[1:2]))
Q

##      a    b
## a -0.2  0.2
## b  0.2 -0.2

## get number of levels of discrete trait for simulation
k<-nrow(Q)
k

## [1] 2

We can go ahead & simulate a discrete character history of our trait. Since I’m actually going to simulate under the null, I could’ve also used sim.Mk here.

## simulate TRUE discrete character history
sim_tree<-sim.history(phy,Q,anc="a")

## Done simulation(s).

sim_tree

## 
## Phylogenetic tree with 60 tips and 59 internal nodes.
## 
## Tip labels:
## 	t9, t10, t53, t54, t19, t49, ...
## 
## The tree includes a mapped, 2-state discrete character
## with states:
## 	a, b
## 
## Rooted; includes branch lengths.

## visualize generating discrete character history
cols<-setNames(hcl.colors(n=k),letters[1:k])
plot(sim_tree,cols,ftype="off",lwd=2,
  direction="upwards")
par(lend=1)
legend("bottomleft",letters[1:k],lwd=3,
  col=hcl.colors(n=k),cex=0.8,bty="n")

So far so good.

Next, I’m going to specify my simulation conditions of the continuous trait. Once again, since I’m actually simulating under the null model it isn’t necessary to use phytools::multiOU here, but I will anyway – just with the same values of \(\alpha\), \(\sigma^2\), and \(\theta\) for each of my discrete character levels.

## set simulation conditions for continuous trait
alpha<-setNames(rep(0.3,k),letters[1:k])
alpha

##   a   b 
## 0.3 0.3

sig2<-setNames(rep(0.1,k),letters[1:k])
sig2

##   a   b 
## 0.1 0.1

theta<-setNames(c(-0.5,-0.5),letters[1:k])
theta

##    a    b 
## -0.5 -0.5

(I’m including the next step only for people who might like to adapt this code to simulate different levels of \(\theta\) for the two different discrete character states.)

## get root state from discrete character history
root_state<-getStates(sim_tree,"nodes")[1]
root_state

##  61 
## "a"

## generating continuous character using multiOU
X<-multiOU(sim_tree,alpha,sig2,theta,
  a0=theta[root_state])
head(X)

##         t9        t10        t53        t54        t19        t49 
## -0.5425203 -0.2690402 -0.4526387 -0.2366645 -0.1525741 -0.9399830

We’ve simulated our discrete and continuous traits; however, we still need to pull our discrete trait off the tree using phytools::getStates.

## pull discrete trait off sim_tree using getStates
Y<-as.factor(getStates(sim_tree,"tips"))
head(Y)

##  t9 t10 t53 t54 t19 t49 
##   a   a   a   a   b   b 
## Levels: a b

OK. Now, to start let’s quickly fit our continuous OU model using geiger::fitContinuous, our discrete model using phytools::fitMk, and then add the log-likelihoods.

## get null log(L) using fitContinuous & fitMk
ou_fit<-geiger::fitContinuous(phy,X,model="OU")
ou_fit

## GEIGER-fitted comparative model of continuous data
##  fitted 'OU' model parameters:
## 	alpha = 0.549895
## 	sigsq = 0.152174
## 	z0 = -0.592858
## 
##  model summary:
## 	log-likelihood = -20.294564
## 	AIC = 46.589128
## 	AICc = 47.017699
## 	free parameters = 3
## 
## Convergence diagnostics:
## 	optimization iterations = 100
## 	failed iterations = 0
## 	number of iterations with same best fit = 43
## 	frequency of best fit = 0.430
## 
##  object summary:
## 	'lik' -- likelihood function
## 	'bnd' -- bounds for likelihood search
## 	'res' -- optimization iteration summary
## 	'opt' -- maximum likelihood parameter estimates

mk_fit<-fitMk(phy,Y,model="ER",pi="equal")
mk_fit

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           a         b
## a -0.188436  0.188436
## b  0.188436 -0.188436
## 
## Fitted (or set) value of pi:
##   a   b 
## 0.5 0.5 
## due to treating the root prior as (a) flat.
## 
## Log-likelihood: -34.29465 
## 
## Optimization method used was "nlminb"
## 
## R thinks it has found the ML solution.

(In an earlier version of this post I had set pi="mle", but then realized that pi="equal" matched our joint model, so I re-ran it.)

null_logL<-logLik(ou_fit)+logLik(mk_fit)
null_logL

## [1] -54.58921
## attr(,"df")
## [1] 3

Having done this, I’m ready to fit this same null model using fitmultiOU and null_model=TRUE. This is a joint model, but in which the discrete character has no effect on our continuous character’s evolutionary mode. Warning: this takes a while!!

## now fit the null model using fitmultiOU
fit_null<-fitmultiOU(phy,X,Y,model="ER",levs=200,
  parallel=TRUE,ncores=10,root="mle",trace=1,
  null_model=TRUE)

## iter	theta	alpha	sigsq	q[1]	log(L)
## 0	0.0991	0.0802	0.0040	0.0662	-430.5097 
## 100	-0.6013	0.3368	0.0960	0.1504	-55.8024 
## 200	-0.6002	0.5436	0.1250	0.2330	-55.1476 
## 300	-0.5919	0.5487	0.1494	0.1888	-54.5434 
## 331	-0.5912	0.5499	0.1496	0.1883	-54.5433 
## Done optimizing.

Here’s our fitted joint model.

fit_null

## Object of class "fitmultiOU" based on
##     a discretization with k = 200 levels.
## 
## Fitted multi-theta OU model parameters:
##  levels: [ a, b ]
##   theta: [ -0.5912 ]
##   alpha: 0.5499 
##   sigsq: 0.1496 
## 
## Estimated Q matrix:
##            a          b
## a -0.1883095  0.1883095
## b  0.1883095 -0.1883095
## 
## Log-likelihood: -54.5433 
## 
## R thinks it has found the ML solution.

Let’s do a quick comparison of parameter estimates. (Even though this is a pretty small tree, so we don’t expect our parameter estimates to match the generating values too closely, I’ll throw those in as well.)

## compare parameter estimates
obj<-cbind(
  setNames(c(alpha[1],sig2[1],theta[1],q),
    c("alpha","sigsq","theta","q")),
  unlist(c(ou_fit$opt[c("alpha","sigsq","z0")],
    q=mk_fit$rates)),
  unlist(list(alpha=fit_null$alpha,
    sigsq=fit_null$sigsq,
    z0=fit_null$theta,
    q=fit_null$rates)))
colnames(obj)<-c("generating","estimated","fitmultiOU")
obj

##       generating  estimated fitmultiOU
## alpha        0.3  0.5498946  0.5498956
## sigsq        0.1  0.1521741  0.1495582
## theta       -0.5 -0.5928578 -0.5912444
## q            0.2  0.1884358  0.1883095

Finally, we can compare log-likelihoods. Once again these should be similar & converge in the limit as we increase levs.

## compare log likelihoods
null_logL ## ignore the df

## [1] -54.58921
## attr(,"df")
## [1] 3

logLik(fit_null)

## [1] -54.5433
## attr(,"df")
## [1] 4

This is pretty cool. We expect both the parameter estimates and the log-likelihoods to increase in similarity with levs and some micro-experimentation (so far) confirms this.

More on this model soon to come!

A discrete-character-dependent multi-θ Ornstein-Uhlenbeck model in phytools

I’m still working out the bugs, but I just added a new function, fitmultiOU, to phytools 2.6-1 on GitHub for fitting a discrete-character-dependent multi-regime Ornstein-Uhlenbeck using the discrete diffusion approximation of Boucher & Démery (2016) – that we also describe in great detail in a new bioRxiv pre-print. I encourage anyone interested in this powerful tool to go & check out our pre-print!

fitmultiOU builds off of the post I created a couple of weeks ago on this blog showing how to calculate the likelihood under a simple, single- \(\theta\) OU model on a tree using the discrete approximation. fitmultiOU just takes this a step further in exactly the same way as fitmultiBM and fitmultiTrend.

As previously noted, I’m definitely still working out the kinks. Indeed, I wouldn’t trust it farther than I can throw it. Hopefully that will change soon. Nonetheless, preliminary tests are promising.

To see this, let’s start by simulating a discrete character on a tree of 1000 taxa that assumes two different levels, a and b, as follows. (I’m choosing to simulate on a relatively large tree on purpose just to give us the best chance of recovering our generating parameter values if our method works.)

library(phytools)

## Loading required package: ape

## Loading required package: maps

N<-1000
tree<-pbtree(n=N,scale=10)
tree

## 
## Phylogenetic tree with 1000 tips and 999 internal nodes.
## 
## Tip labels:
##   t358, t660, t787, t788, t659, t210, ...
## 
## Rooted; includes branch length(s).

q<-0.2
Q<-matrix(c(-q,q,q,-q),2,2,
  dimnames=list(letters[1:2],letters[1:2]))
Q

##      a    b
## a -0.2  0.2
## b  0.2 -0.2

sim_tree<-sim.history(tree,Q,anc="a")

## Done simulation(s).

cols<-setNames(hcl.colors(n=2),letters[1:2])
plot(sim_tree,cols,ftype="off",lwd=1,
  direction="upwards")
par(lend=1)
legend("bottomleft",letters[1:2],lwd=2,
  col=hcl.colors(n=2),cex=0.6,bty="n")

To be clear, we've generated this trait history for simulation -- but for estimation, we'll use only the tip states & the tree.

Next, let’s choose our different simulation parameters. So far, I’m only allowing \(\theta\), the position of th “optimum” in an OU model, to vary as a function of the discrete trait. Nonetheless, our simulator, multiOU, requires that we supply vectors for \(\alpha\), \(\sigma^2\), and \(\theta\), even if we intend to vary only \(\theta\) by discrete character regime.

I’m going to set \(\alpha = 0.3\), which corresponds to a phylogenetic “half-life” of \(\log(2) / \alpha \approx 2.3\) (our total tree height is 10), \(\sigma^2 = 0.1\), and \(\theta = [-0.5, 1.5]\) for the two different regimes.

alpha<-setNames(rep(0.3,2),letters[1:2])
alpha

##   a   b 
## 0.3 0.3

sig2<-setNames(rep(0.1,2),letters[1:2])
sig2

##   a   b 
## 0.1 0.1

theta<-setNames(c(-0.5,1.5),letters[1:2])
theta

##    a    b 
## -0.5  1.5

Now we’re ready for our simulation.

multiOU allows us to separately specify the root state and the root regime, but because I’m not sure these are separately estimable, I’m going to set them to be equal to one another. (That is, I’ll assume that the initial condition of the simulation for the continuous trait is on the “optimum” as prescribed by the discrete character condition at the root.)

## get the root state
root_state<-getStates(sim_tree,"nodes")[1]
root_state

## 1001 
##  "a"

Now I can simulate my continuous trait, x.

x<-multiOU(sim_tree,alpha,sig2,theta,a0=theta[root_state])
head(x)

##     t358     t660     t787     t788     t659     t210 
## 1.336066 1.620262 1.661155 1.800279 1.728212 1.764376

Our analysis is going to use as input the tip states for our discrete trait, \(y\) – not, importantly, the regimes we generated using sim.history, as mentioned previously.

Since only the object sim_tree has this information, we should pull it off using phytools::getStates.

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

## t358 t660 t787 t788 t659 t210 
##    b    b    b    b    b    b 
## Levels: a b

As a final step, I’m going to try to identify reasonable starting values for our optimization process. This will just make finding the MLEs a bit more efficient.

init<-setNames(
  c(
    mean(x[y==levels(y)[1]]),mean(x[y==levels(y)[2]]),
    log(2)/max(nodeHeights(tree)),
    var(x)/max(nodeHeights(tree)),
    fitMk(tree,y)$rates),
  c("theta[a]","theta[b]","alpha","sigsq","q[1]"))
init

##    theta[a]    theta[b]       alpha       sigsq        q[1] 
## -0.01785477  0.80715486  0.06931472  0.05112673  0.23197867

Now let’s try to optimize our model using fitmultiOU. I’ll turn on trace to get periodic progress reports about our optimization process. As a word of warning, this takes a while. (The following analysis took several hours to run on my laptop.)

fit_mou<-fitmultiOU(tree,x,y,model="ER",levs=100,
  parallel=TRUE,root="mle",init=init,trace=1)

## iter	the[a]	the[b]	alpha	sigsq	q[1]	log(L)
## 0	-0.0179	0.8072	0.0693	0.0511	0.2320	-932.8691 
## 100	-0.8977	1.1648	0.3001	0.0832	0.2364	-724.9863 
## 200	-0.6454	1.2625	0.3260	0.0865	0.2239	-724.2291 
## 300	-0.4044	1.3876	0.3480	0.0904	0.2229	-723.3257 
## 400	-0.4038	1.4186	0.3340	0.0856	0.2210	-722.8809 
## 500	-0.3719	1.4540	0.3327	0.0834	0.2263	-722.6922 
## 600	-0.3316	1.6840	0.3190	0.0774	0.2268	-721.7071 
## 700	-0.3316	1.6656	0.3145	0.0814	0.2383	-720.9715 
## 800	-0.3259	1.6652	0.3145	0.0817	0.2379	-720.9686 
## 874	-0.3266	1.6667	0.3141	0.0816	0.2378	-720.9684 
## Done optimizing.

Here's our fitted model.

fit_mou

## Object of class "fitmultiOU" based on
##     a discretization with k = 100 levels.
## 
## Fitted multi-theta OU model parameters:
##  levels: [ a, b ]
##   theta: [ -0.3266, 1.6667 ]
##   alpha: 0.3141 
##   sigsq: 0.0816 
## 
## Estimated Q matrix:
##           a         b
## a -0.237812  0.237812
## b  0.237812 -0.237812
## 
## Log-likelihood: -720.9684 
## 
## R thinks it has found the ML solution.

Well, that’s pretty uncanny as a first try. Again, fitmultiOU is a prototype and is shared only with a warning to use at one’s own risk. Nonetheless, I think it’s right, so that’s pretty cool.

For when I tweet about this, here’s a visualization of multi-regime OU (in general, not this specific case). The code for this visualization comes from a different recent post to this blog.

phy<-pbtree(n=40,scale=10)
tt<-map.to.singleton(make.era.map(phy,
  limits=seq(0,10,length.out=1001)))
Q<-matrix(c(-0.2,0.2,0.2,-0.2),2,2,
  dimnames=list(letters[1:2],letters[1:2]))
s.tt<-sim.history(tt,Q,anc="a")

## Done simulation(s).

cols<-setNames(hcl.colors(n=2),letters[1:2])
theta<-setNames(c(-0.5,1.5),letters[1:2])
sigsq<-setNames(rep(0.1,2),letters[1:2])
alpha<-setNames(rep(0.3,2),letters[1:2])
xx<-multiOU(s.tt,alpha=alpha,sig2=rep(sigsq,2),
  theta=theta,a0=theta["a"],internal=TRUE)
layout(matrix(c(1,2),2,1),heights=c(0.4,0.6))
plot(s.tt,cols,ftype="off",mar=c(1.1,4.1,2.1,2.1),
  xlim=c(0,11))
mtext("a)",adj=0,line=0)
par(mar=c(5.1,4.1,2.1,2.1))
phenogram(s.tt,xx,ftype="off",
  colors=make.transparent(cols,0.5),
  xlim=c(0,11),cex.axis=0.8,las=1)
mtext("b)",adj=0,line=1.5)
for(i in 1:length(theta)){
  stat_theta<-dnorm(seq(min(xx),max(xx),length.out=200),
    mean=theta[i],sd=sqrt(sigsq[i]/(2*alpha[i])))
  stat_theta<-c(0,stat_theta,0)
  stat.norm<-stat_theta/max(stat_theta)*1
  polygon(x=stat.norm+10,
    y=c(min(xx),seq(min(xx),max(xx),
    length.out=200),max(xx)),border=FALSE,
    col=make.transparent(cols[i],0.5))
}

Visualizing a discrete-character-dependent multi-regime OU process on a tree with phytools

Working on something else, I was recently reminded of the fact that phytools includes a function to simulate a multi- \(\theta\) (i.e., multi-“regime”) Ornstein-Uhlenbeck process called multiOU.

multiOU does a continuous-time simulation, but can (optionally) output node states, so I thought it would be fun to show how the function can be used (in combination with lots of non-branching nodes in our tree) to graphically illustrate multi-optimum OU.

This doesn’t serve any special purpose outside the visual, but perhaps someone currently preparing to present at Evolution 2026 in Cleveland in a couple of weeks will find this useful for their slideshow!

## load phytools
library(phytools)

We can start by simulating a tree. Here I’m going to scale the total length of this tree to 100 time units.

tree<-pbtree(n=40,scale=100)

Next I’d like to add a bunch of non-branching nodes along the length of my tree from the root to every tip. I’m going to evenly space these & the way I’m going to do it is by wrapping a function call of phytools::make.era.map with another map.to.singleton. This first generates “era” mapped regimes on the tree in "simmap" format and then converts these to nodes.

tt<-map.to.singleton(make.era.map(tree,
  limits=seq(0,100,length.out=1001)))

We can see that this worked as follows.

plotTree(tt,ftype="off",lwd=1)
get("last_plot.phylo",envir=.PlotPhyloEnv)->pp
points(pp$xx,pp$yy,pch=16,col=make.transparent("blue",0.2),
  cex=0.5)

(I hope it’s evident that the blur of blue dots along each edge of the tree & non-branching nodes!)

Next we’re going to simulate our regimes on the tree as a continuous-time Markov chain using phytools::sim.history. I’m just going to simulate two regimes, a and b, with a switching rate between regimes of \(q=0.02\). Keep in mind, we should do this on our tree with unbranching nodes, tt.

## our Q matrix
Q<-matrix(c(-0.02,0.02,0.02,-0.02),2,2,
  dimnames=list(letters[1:2],letters[1:2]))
## our regime tree
s.tt<-sim.history(tt,Q,anc="a")

## Done simulation(s).

To get a sense of the history we’ve simulated, let’s plot it. This can be done in the standard way because phytools::plot.simmap has no problem with unbranching nodes.

cols<-setNames(hcl.colors(n=2),letters[1:2])
plot(s.tt,cols,ftype="off")

So far, so good.

Now, we’ll set our simulation conditions for the multi-regime OU. The function allows us to simulate different \(\alpha\), \(\sigma^2\), and \(\theta\) for each of our regimes, but here we’ll hold the former two parameters, \(\alpha\) & \(\sigma^2\), constant between regimes and set \(\theta = [-1,1]\) as follows.

theta<-setNames(c(-1,1),letters[1:2])
sigsq<-setNames(rep(0.4,2),letters[1:2])
alpha<-setNames(rep(0.7,2),letters[1:2])
x<-multiOU(s.tt,alpha=alpha,sig2=rep(sigsq,2),
  theta=theta,a0=theta["a"],internal=TRUE)

Note that the vector, x, that I simulated contains not only the states for the tips but also the conditions of all the thousands of internal nodes of the tree – the overwhelming majority of which are unbranching.

head(x)

##          t4         t25         t26          t8          t9         t18 
## -1.47750858  0.01917008  1.05348083  1.71793727  0.71903373  1.17666809

length(x)

## [1] 14574

Now, for our plot!

For fun, I’m also going to add the stationary distribution at the tips based on our generating values of \(\theta=[-1,1]\), \(\alpha=0.7\), and \(\sigma^2=0.4\).

layout(matrix(c(1,2),2,1),heights=c(0.4,0.6))
plot(s.tt,cols,ftype="off",mar=c(1.1,4.1,2.1,2.1),
  xlim=c(0,110))
mtext("a)",adj=0,line=0)
par(mar=c(5.1,4.1,2.1,2.1))
phenogram(s.tt,x,ftype="off",
  colors=make.transparent(cols,0.5),
  xlim=c(0,110),cex.axis=0.8,las=1)
mtext("b)",adj=0,line=1.5)
for(i in 1:length(theta)){
  stat_theta<-dnorm(seq(min(x),max(x),length.out=200),
    mean=theta[i],sd=sqrt(sigsq[i]/(2*alpha[i])))
  stat.norm<-stat_theta/max(stat_theta)*10
  polygon(x=stat.norm+100,y=seq(min(x),max(x),
    length.out=200),border=FALSE,
    col=make.transparent(cols[i],0.5))
}

That’s pretty much the idea. Cool, right?