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?

Ancestral states under the discretized diffusion approximation using bounded_bm in phytools

As many followers & readers of this blog will undoubtedly know, I’ve been thinking & writing endlessly about what we refer to as the “discretized diffusion approximation” for computing the likelihood in continuous & discrete character models (first instroduced to phylogenetic comparative biology by Boucher & Démery 2016) for well over two years. We recently published a pre-print that I encourage anyone interested in this topic to check out!

In this post, I just wanted to briefly explore ancestral state estimation using this approach which I’ve implemented in the phytools ancr generic method.

First, let’s just see how it works in general. To do that, I’m going to simulate a tree & some data via unbounded Brownian motion under which our MLEs for ancestral states would normally be very easy to obtain (e.g., using phytools::fastAnc).

## load phytools
library(phytools)

## simulate a tree
phy<-pbtree(n=26,scale=1,tip.label=LETTERS)
phy

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##   A, B, C, D, E, F, ...
## 
## Rooted; includes branch length(s).

## simulate some trait data
x<-fastBM(phy)
x

##           A           B           C           D           E           F           G 
## -0.52238323 -1.08476410 -0.73893741 -0.75336367 -0.94397490 -0.67163255 -1.21387274 
##           H           I           J           K           L           M           N 
## -0.65653756 -1.16240551 -0.87838336 -0.93716934 -0.80657269 -0.96539003 -1.09064863 
##           O           P           Q           R           S           T           U 
## -1.80940444 -0.99874761 -2.07552774 -1.83396425 -0.77599997 -2.04509690 -0.67626465 
##           V           W           X           Y           Z 
## -0.09662133  0.57076545  0.47261265  0.26072599 -0.03490109

Next, let’s fit an unbounded Brownian model using the discrete approximation. This can be done in phytools with the function bounded_bm just by allowing the limits or bounds of the trait to fall well outside the observed trait distribution. (This is technically a “bounded” model, but since almost no probability density hits the bounds it is effectively unbounded BM.)

## fit unbounded model
bm<-bounded_bm(phy,x,root="mle",levs=400)
bm

## Object of class "bounded_bm" based on
##    a discretization with k = 400 levels.
## 
## Unwrapped (i.e., bounded) model
## 
## Set or estimated bounds: [ -3.398674 , 1.893912 ]
## 
## Fitted model parameters:
##   sigsq: 0.639077 
##      x0: -0.69284 
## 
## Log-likelihood: -17.871792 
## 
## R thinks it has found the ML solution.

If we compare this to the fit obtained using the continuous likelihood function (rather than our discretized approximation) we should find that the parameter estimates & log-likelihood (very, very nearly) match. This similar will only increase with levs (or equivalently, as \(\delta\), the bin width of our discretization, decreases towards zero).

geiger::fitContinuous(phy,x,model="BM")

## GEIGER-fitted comparative model of continuous data
##  fitted 'BM' model parameters:
## 	sigsq = 0.642594
## 	z0 = -0.692552
## 
##  model summary:
## 	log-likelihood = -17.940158
## 	AIC = 39.880316
## 	AICc = 40.402055
## 	free parameters = 2
## 
## Convergence diagnostics:
## 	optimization iterations = 100
## 	failed iterations = 0
## 	number of iterations with same best fit = 100
## 	frequency of best fit = 1.000
## 
##  object summary:
## 	'lik' -- likelihood function
## 	'bnd' -- bounds for likelihood search
## 	'res' -- optimization iteration summary
## 	'opt' -- maximum likelihood parameter estimates

This is old news.

What you might not have realized is that our model object from bounded_bm can be passed to the generic method ancr and we ought to get back out MLE ancestral states – with a small catch that I’ll get to!

This is what I mean….

## ancestral state estimates from the discrete approximation
anc.d<-ancr(bm)
print(anc.d,printlen=5)

## Ancestral character estimates from "bounded_bm" object:
##        27        28        29        30        31     
##  -0.69284 -0.909587 -0.809253 -0.858961 -0.854084 ....
## 
## Lower & upper 95% CIs:
##        lower     upper
## 27  -0.69284  -0.69284
## 28 -1.394107  -0.42821
## 29 -1.208867 -0.414979
## 30 -1.169172 -0.547293
## 31 -1.155941 -0.547293
##         ....      ....

## ancestral states using continuous likelihood function
anc.c<-fastAnc(phy,x,CI=TRUE)
print(anc.c,printlen=5)

## Ancestral character estimates using fastAnc:
##         27       28        29        30       31     
##  -0.692552 -0.90973 -0.810438 -0.860237 -0.85374 ....
## 
## Lower & upper 95% CIs:
##        lower     upper
## 27 -1.378038 -0.007065
## 28 -1.438079  -0.38138
## 29 -1.220804 -0.400071
## 30 -1.174338 -0.546135
## 31 -1.167611 -0.539869
##         ....      ....

Let’s plot one set of point estimates against the other.

par(mar=c(5.1,4.1,1.1,1.1))
plot(anc.c$ace,anc.d$ace,bty="n",pch=21,
  bg=make.transparent("blue",0.5),
  xlab="estimates from fastAnc",
  ylab="estimates from bounded_bm > ancr",
  las=1,cex.axis=0.8)
grid()
abline(a=0,b=1,lty="dashed")

I don’t need any fancy statistics to see that I’m getting the right estimates here!

There is one teeny, tiny problem though.

It becomes evident when we closely inspect our model object from ancr as follows.

print(anc.d,prinlen=5)

## Ancestral character estimates from "bounded_bm" object:
##        27        28        29        30        31        32     
##  -0.69284 -0.909587 -0.809253 -0.858961 -0.854084 -0.884419 ....
## 
## Lower & upper 95% CIs:
##        lower     upper
## 27  -0.69284  -0.69284
## 28 -1.394107  -0.42821
## 29 -1.208867 -0.414979
## 30 -1.169172 -0.547293
## 31 -1.155941 -0.547293
## 32 -1.169172 -0.600219
##         ....      ....

Here we see that the upper & lower bounds of the confidence intervals on the global root state are both equal to the point estimate which is clearly wrong. Indeed, this is typically the most uncertain internal node state, rather than the least!

This happens because to compute the likelihood under our discretized diffusion approximation and have it match the continuous density, when we get to the root state we have to fix it to its MLE.

This has the effect of meaning that when we go ahead & undertake ASR using this fitted model, instead of just conditioning on the MLE of \(\sigma^2\), we’re also conditioning on the MLE of \(x_0\), the root state. This is not what we want to be doing at all!

Well, in this case it doesn’t really matter because we can solve the continuous density anyway; however, for other models where that’s intractable, what I might tentatively recommend instead is changing the root prior in estimation if we intend to use our fitted model in ASR.

Here’s what that would look like for our twenty-six taxon phylogeny of earlier.

bm_flatroot<-bounded_bm(phy,x,levs=400,root="flat")
bm_flatroot

## Object of class "bounded_bm" based on
##    a discretization with k = 400 levels.
## 
## Unwrapped (i.e., bounded) model
## 
## Set or estimated bounds: [ -3.398674 , 1.893912 ]
## 
## Fitted model parameters:
##   sigsq: 0.664673 
##      x0: -0.692566 
## 
## Log-likelihood: -19.682352 
## 
## R thinks it has found the ML solution.

anc.d_flatroot<-ancr(bm_flatroot)
print(anc.d_flatroot,printlen=5)

## Ancestral character estimates from "bounded_bm" object:
##         27        28        29        30        31     
##  -0.692566 -0.909521 -0.809241 -0.858956 -0.854082 ....
## 
## Lower & upper 95% CIs:
##        lower     upper
## 27 -1.380876 -0.004803
## 28 -1.433802 -0.388516
## 29 -1.222098 -0.401747
## 30 -1.169172 -0.547293
## 31 -1.169172 -0.547293
##         ....      ....

Let’s once again compare these states (that is, the point estimates at nodes) to the continuous likelihood function values.

par(mar=c(5.1,4.1,1.1,1.1))
plot(anc.c$ace,anc.d_flatroot$ace,bty="n",pch=21,
  bg=make.transparent("blue",0.5),
  xlab="estimates from fastAnc",
  ylab="estimates from bounded_bm > ancr",
  las=1,cex.axis=0.8)
grid()
abline(a=0,b=1,lty="dashed")

Once again, it’s evident that we have the right ML states.

Now, for fun, let’s compare the lower & upper confidence bounds as well….

par(mfrow=c(1,2),mar=c(5.1,4.1,2.1,1.1))
plot(anc.c$CI95[,1],anc.d_flatroot$CI95[,1],bty="n",pch=21,
  bg=make.transparent("blue",0.5),
  xlab="lower CI bound from fastAnc",
  ylab="lower CI bound from bounded_bm > ancr",
  las=1,cex.axis=0.8)
grid()
abline(a=0,b=1,lty="dashed")
mtext("a) lower CI bound",adj=0)
plot(anc.c$CI95[,2],anc.d_flatroot$CI95[,2],bty="n",pch=21,
  bg=make.transparent("blue",0.5),
  xlab="upper CI bound from fastAnc",
  ylab="upper CI bound from bounded_bm > ancr",
  las=1,cex.axis=0.8)
grid()
abline(a=0,b=1,lty="dashed")
mtext("b) upper CI bound",adj=0)

Dang. I didn’t expect that to work this well!

Nice.

Of course, this actually gets interesting when we start to introduce models (such as bounded BM) without tractable solutions for the likelihood. Perhaps I’ll demo that in another post!

Generic plot method for the "fastAnc" object class in phytools

Sometimes I forget about all the hidden features of phytools.

One that I’d nearly forgotten about today is an S3 generic plot method that I wrote for the "fastAnc" continuous character ancestral state reconstruction method object class that I wrote just a few months ago.

Here’s a very simple demo.

## load phytools
library(phytools)

## load a tree & some data
data(primate.data)
data(primate.tree)

head(primate.data)

##                                  Group Skull_length Optic_foramen_area Orbit_area
## Allenopithecus_nigroviridis Anthropoid         98.5                7.0      298.7
## Alouatta_palliata           Anthropoid        109.8                5.3      382.3
## Alouatta_seniculus          Anthropoid        108.0                8.0      359.4
## Aotus_trivirgatus           Anthropoid         60.5                3.1      297.4
## Arctocebus_aureus            Prosimian         49.5                1.2      134.8
## Arctocebus_calabarensis      Prosimian         53.8                1.6      156.6
##                             Activity_pattern Activity_pattern_code
## Allenopithecus_nigroviridis          Diurnal                     0
## Alouatta_palliata                    Diurnal                     0
## Alouatta_seniculus                   Diurnal                     0
## Aotus_trivirgatus                  Nocturnal                     2
## Arctocebus_aureus                  Nocturnal                     2
## Arctocebus_calabarensis            Nocturnal                     2

Let’s analyze primate skull length, on a log-scale.

lnSkullLength<-setNames(log(primate.data$Skull_length),
  rownames(primate.data))

Here’s our ancestral state reconstruction under BM.

primate_anc<-fastAnc(primate.tree,lnSkullLength,CI=TRUE)
print(primate_anc,printlen=6)

## Ancestral character estimates using fastAnc:
##        91       92       93       94       95       96     
##  4.175934 4.186044 4.431545 4.724982 4.681729 4.665113 ....
## 
## Lower & upper 95% CIs:
##       lower    upper
## 91 3.814174 4.537694
## 92 3.828175 4.543914
## 93 4.116413 4.746677
## 94 4.477337 4.972627
## 95 4.483873 4.879584
## 96 4.494131 4.836094
##        ....     ....

Now let’s use the plot method, as follows. Note that most of the arguments of both phytools::plotTree and phytools::add.color.bar can be passed internally.

h<-max(nodeHeights(primate.tree))
plot(primate_anc,direction="upwards",ftype="i",
  fsize=0.6,title="log(skull length)",
  offset=1,legend="bottomleft",
  ylim=c(-0.1*h,1.35*h))

That’s pretty useful, I think!

ansi_phylo ANSI text phylogeny simulator in R

While I was taking a break from fiddling around with the endeavor of implementing a discrete-character-dependent multi-$\theta$ OU model (if such a thing is even possible) using the discretized diffusion approximation, I stumbled on an old phytools called ansi_phylo for ANSI text style phylogenetic tree plotting in R.

It creates a pretty fun visual, so I thought it might be worth a re-post. Here’s a quick demo:

## load packages
library(phytools)
library(phangorn)

## load data
data(Laurasiatherian)
Laurasiatherian

## 47 sequences with 3179 character and 1605 different site patterns.
## The states are a c g t

## estimate tree using ML in phangorn
mammal_mle<-pml_bb(Laurasiatherian,model="K80")

## optimize edge weights:  -54315.47 --> -54130.45 
## optimize rate matrix:  -54130.45 --> -51361.62 
## ...
## optimize topology:  -52200.08 --> -51341.22  NNI moves:  13 
## Ratchet iteration  100 , best pscore so far: -51341.22 
## optimize rate matrix:  -51341.22 --> -51341.22 
## optimize edge weights:  -51341.22 --> -51341.22 
## optimize topology:  -51341.22 --> -51341.22  NNI moves:  0

mammal_mle

## model: K80 
## loglikelihood: -51341.22 
## unconstrained loglikelihood: -17300.92 
## 
## Rate matrix:
##          a        c        g        t
## a 0.000000 1.000000 5.025414 1.000000
## c 1.000000 0.000000 1.000000 5.025414
## g 5.025414 1.000000 0.000000 1.000000
## t 1.000000 5.025414 1.000000 0.000000
## 
## Base frequencies:  
##    a    c    g    t 
## 0.25 0.25 0.25 0.25

## get ML tree and root using Platypus as outgroup
mammal_mle.phy<-root(mammal_mle$tree,outgroup="Platypus",
  resolve.root=TRUE)

ansi_phylo(mammal_mle.phy,horizontal="=")

That’s it!

Here, let’s try another one.

data(yeast)
yeast

## 8 sequences with 127026 character and 8899 different site patterns.
## The states are a c g t

yeast_mle<-pml_bb(yeast,model="K80")

## optimize edge weights:  -739078.1 --> -734615.7 
## optimize rate matrix:  -734615.7 --> -715533.2 
## ...
## optimize topology:  -718794.8 --> -715391.4  NNI moves:  1 
## Ratchet iteration  100 , best pscore so far: -715391.4 
## optimize rate matrix:  -715391.4 --> -715387.4 
## optimize edge weights:  -715387.4 --> -715387.3 
## optimize topology:  -715387.3 --> -715387.3  NNI moves:  0 
## optimize rate matrix:  -715387.3 --> -715387.3 
## optimize edge weights:  -715387.3 --> -715387.3

yeast_mle

## model: K80 
## loglikelihood: -715387.3 
## unconstrained loglikelihood: -659104.8 
## 
## Rate matrix:
##          a        c        g        t
## a 0.000000 1.000000 3.291576 1.000000
## c 1.000000 0.000000 1.000000 3.291576
## g 3.291576 1.000000 0.000000 1.000000
## t 1.000000 3.291576 1.000000 0.000000
## 
## Base frequencies:  
##    a    c    g    t 
## 0.25 0.25 0.25 0.25

rooted.yeast_tree<-midpoint(yeast_mle$tree)

ansi_phylo(rooted.yeast_tree,vertical="*")