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!