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!