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:
## t8, t10, t11, t6, t1, t34, ...
##
## 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:
## t8, t10, t11, t6, t1, t34, ...
##
## 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)
## t8 t10 t11 t6 t1 t34
## -0.4708026 -0.4431907 -0.8007267 -0.6463636 -0.1806929 -0.8320607
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)
## t8 t10 t11 t6 t1 t34
## a a a a a 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.351120
## sigsq = 0.083454
## z0 = -0.737954
##
## model summary:
## log-likelihood = -6.492717
## AIC = 18.985434
## AICc = 19.414005
## free parameters = 3
##
## Convergence diagnostics:
## optimization iterations = 100
## failed iterations = 0
## number of iterations with same best fit = 56
## frequency of best fit = 0.560
##
## 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="mle")
mk_fit
## Object of class "fitMk".
##
## Fitted (or set) value of Q:
## a b
## a -0.076687 0.076687
## b 0.076687 -0.076687
##
## Fitted (or set) value of pi:
## a b
## 1 0
## due to treating the root prior as (a) it's MLE.
##
## Log-likelihood: -23.216988
##
## Optimization method used was "nlminb"
##
## R thinks it has found the ML solution.
null_logL<-logLik(ou_fit)+logLik(mk_fit)
null_logL
## [1] -29.7097
## 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.5320 0.0823 0.0279 0.0918 -43.0590
## 100 -0.7229 0.2822 0.0732 0.0646 -30.6777
## 200 -0.7393 0.3519 0.0827 0.0822 -30.3418
## 249 -0.7396 0.3518 0.0827 0.0820 -30.3418
## 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.7396 ]
## alpha: 0.3518
## sigsq: 0.0827
##
## Estimated Q matrix:
## a b
## a -0.08204587 0.08204587
## b 0.08204587 -0.08204587
##
## Log-likelihood: -30.3418
##
## 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.35111999 0.35181166
## sigsq 0.1 0.08345362 0.08270091
## theta -0.5 -0.73795424 -0.73960081
## q 0.2 0.07668651 0.08204587
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] -29.7097
## attr(,"df")
## [1] 3
logLik(fit_null)
## [1] -30.34178
## 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. I'm still trying to figure out a couple of things about how to deal with the root.
More on this model soon to come!
