A couple of days ago on this blog I posted an entry on comparing a discrete-character-dependent multi-regime joint OU model to OUwie when the regime history is known.
For those not following closely at home, the joint discrete & continuous trait model that I’m talking about uses a finite-space / discretized diffusion approximation as documented in prior posts to this blog (e.g., 1, 2, 3, 4, 5), following Boucher & Démery (2016; also see our in review bioRxiv pre-print).
A more typical workflow with OUwie (and other similar regime-based methods), is to generate a set of stochastic character maps of the discrete character regime following Huelsenbeck et al. (2003), fit a multi-regime OU model to each tree in the set, and then average the results across maps. Back in 2013, I pointed out that (for rate-heterogeneous models) this would tend to result in an underestimate of the difference in rates, \(\sigma^2\), between regimes. The same seems very likely to be true for multi-regime OU models, but (to my knowledge) this has not been studied explicitly. We might hope or expect that this bias would go away if we can jointly model the discrete & continuous characters.
The magnitude of this bias is also likely to vary as a function of the transition rate for the discrete trait. In other words, the higher the transition rates in Q of our discrete trait, the more uncertain our discrete character history and the more it will vary among stochastic character maps (and the more each map will differ from the unknown true history of our trait).
On a Sunday morning, I thought it would be interesting to engineer a simple comparison between this traditional OUwie workflow and fitmultiOU both for circumstances in which the transition rate of the discrete character is low & high.
To begin with, let’s load the packages we intend to use.
## load packages
library(phytools)
library(OUwie)
Now, I’m going to simulate a tree & a discrete character history on that tree with a low rate (keeping in mind, of course, that there is no special numeric value of \(q\) that is “low” or “high” – it all depends on the total depth of our tree).
phy<-pbtree(n=250,scale=10)
phy
##
## Phylogenetic tree with 250 tips and 249 internal nodes.
##
## Tip labels:
## t16, t21, t123, t124, t38, t49, ...
##
## Rooted; includes branch length(s).
q<-0.05
k<-2
low.Q<-matrix(q,k,k,dimnames=list(letters[1:k],letters[1:k]))
diag(low.Q)<-0
diag(low.Q)<--rowSums(low.Q)
low.Q
## a b
## a -0.05 0.05
## b 0.05 -0.05
low.sim_tree<-sim.history(phy,low.Q,
anc=y0<-sample(letters[1:k],1))
## Done simulation(s).
low.sim_tree
##
## Phylogenetic tree with 250 tips and 249 internal nodes.
##
## Tip labels:
## t16, t21, t123, t124, t38, t49, ...
##
## The tree includes a mapped, 2-state discrete character
## with states:
## a, b
##
## Rooted; includes branch lengths.
## plot this tree (just for fun)
cols<-setNames(hcl.colors(n=k),letters[1:k])
plot(low.sim_tree,cols,ftype="off",lwd=1,
direction="upwards",ylim=c(-1,10))
par(lend=1)
legend("bottomleft",letters[1:k],lwd=4,
col=hcl.colors(n=k),cex=0.7,bty="n")
Great. Clearly we’ve simulated a very low rate of discrete character evolution on this tree!
Now let’s simulate a multi-regime OU process conditioning on this true discrete character history. To do this we’ll use multiOU in phytools. Even though we’re using the same value of \(\alpha\) and \(\sigma^2\) for both of our discrete character regimes we need to set values of these two parameters for each of the two character levels.
## set alpha
alpha<-setNames(rep(0.6,k),letters[1:k])
alpha
## a b
## 0.6 0.6
## set sigma-squared
sig2<-setNames(rep(0.3,k),letters[1:k])
sig2
## a b
## 0.3 0.3
Now let’s set \(\theta_a\) and \(\theta_b\). I’m going to use values of \(0\) and \(5\) for the two different discrete character states.
# set theta
theta<-setNames(c(0,5),letters[1:k])
theta
## a b
## 0 5
Now, I’ll simulate continuous trait data under the multi-regime model using our true discrete character history with phytools::multiOU. I’m going to assume the starting value is at \(\theta\) for the root regime. (This is where the a0=theta[y0] comes from!)
## simulate continuous trait
low.x<-multiOU(low.sim_tree,alpha,sig2,theta,a0=theta[y0])
head(low.x)
## t16 t21 t123 t124 t38 t49
## 0.6548297 0.7918578 -0.1488759 0.4502061 0.3959087 0.2302126
Terrific. So far so good.
Now we’re about to imagine that our discrete character history is unobserved, so let’s pull of just the tip values of our discrete trait from low.sim_tree as follows.
low.y<-as.factor(getStates(low.sim_tree,type="tips"))
head(low.y)
## t16 t21 t123 t124 t38 t49
## a a a a a a
## Levels: a b
Next, let’s generate a set of stochastic character histories for our discrete state using phytools::simmap.
low.mk_fit<-fitMk(phy,low.y,model="ER",
pi="fitzjohn")
low.mk_fit
## Object of class "fitMk".
##
## Fitted (or set) value of Q:
## a b
## a -0.036172 0.036172
## b 0.036172 -0.036172
##
## Fitted (or set) value of pi:
## a b
## 0.32705 0.67295
## due to treating the root prior as (a) nuisance.
##
## Log-likelihood: -76.278772
##
## Optimization method used was "nlminb"
##
## R thinks it has found the ML solution.
low.smps<-simmap(low.mk_fit)
low.smps
## 100 phylogenetic trees with mapped discrete characters
Let’s plot a few of these.
par(mfrow=c(5,5))
nulo<-sapply(low.smps[1:25],plot,lwd=1,col=cols,ftype="off",
direction="upwards")
Something that readers will probably notice immediately is that all of these discrete character histories fairly closely resemble each other & the true history of our trait.
Next, let’s fit OUwie models to each of them. That’ll be fun. To speed it along, I’m going to parallelize across stochastic map trees using the foreach package.
## make our OUwie data frame
low.ouwie_data<-data.frame(
Genus_species=names(low.x),
Reg=low.y,
X=low.x)
head(low.ouwie_data)
## Genus_species Reg X
## t16 t16 a 0.6548297
## t21 t21 a 0.7918578
## t123 t123 a -0.1488759
## t124 t124 a 0.4502061
## t38 t38 a 0.3959087
## t49 t49 a 0.2302126
## load foreach and doParallel
library(foreach)
library(doParallel)
## set up our cluster
ncores<-detectCores()-2
ncores
## [1] 14
mc<-makeCluster(ncores,type="PSOCK")
mc
## socket cluster with 14 nodes on host 'localhost'
registerDoParallel(cl=mc)
## optimize across all 100 stochastic maps
low.ouwie_fits<-foreach(i=1:length(low.smps))%dopar%{
OUwie::OUwie(low.smps[[i]],low.ouwie_data,model="OUM",
simmap.tree=TRUE,root.station=FALSE)
}
stopCluster(mc)
The object we’ve created is a long list of fitted “OUwie” models. We can pull out just the specific results that we’re interested in to make it easier to see the general pattern.
foo<-function(x) setNames(
c(
x$theta[,1],
x$solution[,1],
x$loglik
),
c(
"the[a]","the[b]","alpha","sigsq","log(L)"
)
)
low.ouwie_results<-t(sapply(low.ouwie_fits,foo))
rownames(low.ouwie_results)<-1:nrow(low.ouwie_results)
options(scipen=5) ## just for printing
round(low.ouwie_results[1:25,],digits=2)
## the[a] the[b] alpha sigsq log(L)
## 1 0.08 4.84 0.66 0.42 -182.89
## 2 0.11 4.95 0.49 0.36 -194.30
## 3 0.19 5.07 0.41 0.41 -224.99
## 4 0.20 4.95 0.63 0.48 -206.92
## 5 0.00 4.75 0.67 0.52 -209.01
## 6 0.11 4.95 0.44 0.40 -216.47
## 7 0.08 4.87 0.47 0.41 -214.46
## 8 0.05 4.86 0.48 0.36 -192.84
## 9 0.01 4.91 0.40 0.41 -227.67
## 10 0.03 4.82 0.62 0.41 -186.94
## 11 0.07 4.89 0.50 0.36 -193.14
## 12 -0.10 4.86 0.33 0.38 -232.95
## 13 0.03 4.86 0.59 0.37 -178.45
## 14 -0.03 4.78 0.43 0.38 -213.53
## 15 0.03 4.88 0.48 0.31 -175.21
## 16 0.10 4.97 0.50 0.33 -180.42
## 17 0.17 4.95 0.49 0.48 -228.08
## 18 0.12 4.92 0.60 0.48 -210.78
## 19 0.04 4.82 0.48 0.37 -200.46
## 20 -0.07 4.78 0.42 0.34 -199.99
## 21 0.21 5.29 0.34 0.41 -239.15
## 22 0.05 4.87 0.57 0.40 -193.95
## 23 0.14 5.06 0.47 0.39 -207.78
## 24 0.22 5.21 0.37 0.39 -227.25
## 25 0.10 4.90 0.61 0.43 -196.33
options(scipen=0)
This is pretty cool because it shows that all of our fitted models are pretty close to the generating conditions of \(\mathbf{\theta} = [0, 5]\).
colMeans(low.ouwie_results)
## the[a] the[b] alpha sigsq log(L)
## 0.06588211 4.91247554 0.49325269 0.40862951 -208.88237689
We can also compare this to phytools::fitmultiOU, which I predict will get a similar model.
init<-setNames(
c(
mean(low.x[low.y==levels(low.y)[1]]),
mean(low.x[low.y==levels(low.y)[2]]),
log(2)/max(nodeHeights(phy)),
var(low.x)/max(nodeHeights(phy)),
fitMk(phy,low.y,model="ER")$rates),
c("theta[a]","theta[b]",
"alpha","sigsq","q[1]"))
init
## theta[a] theta[b] alpha sigsq q[1]
## 0.14104895 4.40690078 0.06931472 0.52619893 0.03687111
low.fit_mou<-fitmultiOU(phy,low.x,low.y,model="ER",
levs=100,parallel=TRUE,ncores=ncores,root="mle",
trace=1,maxit=2000,init=init)
## iter the[a] the[b] alpha sigsq q[1] log(L)
## 0 0.1410 4.4069 0.0693 0.5262 0.0369 -378.6136
## 100 -0.3327 4.7524 0.3626 0.1910 0.0282 -258.4337
## 200 0.0864 4.8329 0.5574 0.2065 0.0202 -247.0498
## 300 0.1137 4.9956 0.5351 0.2351 0.0203 -245.0731
## 400 0.0581 4.8858 0.6185 0.2743 0.0387 -242.3865
## 500 0.0999 4.9000 0.6364 0.2739 0.0330 -241.9284
## 600 0.1005 4.9045 0.6323 0.2720 0.0331 -241.9258
## 666 0.1009 4.9043 0.6320 0.2719 0.0331 -241.9258
## Done optimizing.
low.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.1009, 4.9043 ]
## alpha: 0.632
## sigsq: 0.2719
##
## Estimated Q matrix:
## a b
## a -0.03312624 0.03312624
## b 0.03312624 -0.03312624
##
## Log-likelihood: -241.9258
##
## R thinks it has found the ML solution.
Even though both analyses gave as parameter estimates quite close to the generating values, joint estimation seems to be even more accurate... particularly with regard to \(\alpha\) and \(\sigma^2\). That's cool.
Now let’s consider the case of a high rate of transition for our discrete character.
q<-0.8
k<-2
high.Q<-matrix(q,k,k,dimnames=list(letters[1:k],
letters[1:k]))
diag(high.Q)<-0
diag(high.Q)<--rowSums(high.Q)
high.Q
## a b
## a -0.8 0.8
## b 0.8 -0.8
Now the discrete character is expected to change a lot in the true history.
high.sim_tree<-sim.history(phy,high.Q,
anc=y0)
## Done simulation(s).
high.sim_tree
##
## Phylogenetic tree with 250 tips and 249 internal nodes.
##
## Tip labels:
## t16, t21, t123, t124, t38, t49, ...
##
## The tree includes a mapped, 2-state discrete character
## with states:
## a, b
##
## Rooted; includes branch lengths.
## plot this tree (just for fun)
cols<-setNames(hcl.colors(n=k),letters[1:k])
plot(high.sim_tree,cols,ftype="off",lwd=1,
direction="upwards",ylim=c(-1,10))
par(lend=1)
legend("bottomleft",letters[1:k],lwd=4,
col=hcl.colors(n=k),cex=0.7,bty="n")
(Let’s generate continuous trait data under the same \(\alpha\), \(\sigma^2\) and \(\mathbf{\theta}}=[0,5]\) as before, but with our new, high Q history.)
## simulate continuous trait
high.x<-multiOU(high.sim_tree,alpha,sig2,theta,
a0=theta[y0])
head(high.x)
## t16 t21 t123 t124 t38 t49
## 1.142639 1.634035 1.209871 1.537554 2.871243 -1.072081
Of course, now we might expect our stochastic character histories to differ more from each other and from the true history.
high.y<-as.factor(getStates(high.sim_tree,type="tips"))
head(high.y)
## t16 t21 t123 t124 t38 t49
## a a a a a a
## Levels: a b
high.mk_fit<-fitMk(phy,high.y,model="ER",
pi="fitzjohn")
high.mk_fit
## Object of class "fitMk".
##
## Fitted (or set) value of Q:
## a b
## a -0.458785 0.458785
## b 0.458785 -0.458785
##
## Fitted (or set) value of pi:
## a b
## 0.499556 0.500444
## due to treating the root prior as (a) nuisance.
##
## Log-likelihood: -159.945824
##
## Optimization method used was "nlminb"
##
## R thinks it has found the ML solution.
high.smps<-simmap(high.mk_fit)
high.smps
## 100 phylogenetic trees with mapped discrete characters
Much as we did earlier, let’s plot a few of these just to see what we’ve got!
par(mfrow=c(5,5))
nulo<-sapply(high.smps[1:25],plot,lwd=1,col=cols,
ftype="off",direction="upwards")
It should be pretty evident, I think, that (even though they seem to have arisen under the same process) in this case the specific details of our discrete character histories vary quite widely one from the other, as well as from our generating history.
I expect that the consequence of this will be that each one of our stochastic character histories is likely to add error (and perhaps bias) to the estimating of the continuous trait multi-regime OU process. Let’s see if that’s true.
## make our OUwie data frame
high.ouwie_data<-data.frame(
Genus_species=names(high.x),
Reg=high.y,
X=high.x)
head(high.ouwie_data)
## Genus_species Reg X
## t16 t16 a 1.142639
## t21 t21 a 1.634035
## t123 t123 a 1.209871
## t124 t124 a 1.537554
## t38 t38 a 2.871243
## t49 t49 a -1.072081
Once again, we’ll parallelize across maps using foreach::foreach.
mc<-makeCluster(ncores,type="PSOCK")
registerDoParallel(cl=mc)
## optimize across all 100 stochastic maps
high.ouwie_fits<-foreach(i=1:length(high.smps))%dopar%{
OUwie::OUwie(high.smps[[i]],high.ouwie_data,model="OUM",
simmap.tree=TRUE,root.station=FALSE)
}
stopCluster(mc)
Let’s summarize our results.
foo<-function(x) setNames(
c(
x$theta[,1],
x$solution[,1],
x$loglik
),
c(
"the[a]","the[b]","alpha","sigsq","log(L)"
)
)
high.ouwie_results<-t(sapply(high.ouwie_fits,foo))
rownames(high.ouwie_results)<-1:nrow(high.ouwie_results)
(Here’s just the first 25, again.)
options(scipen=5) ## just for printing
round(high.ouwie_results[1:25,],digits=2)
## the[a] the[b] alpha sigsq log(L)
## 1 0.67 3.66 0.20 0.73 -350.92
## 2 0.64 3.55 0.19 0.74 -352.68
## 3 0.55 3.73 0.19 0.71 -349.64
## 4 1.11 3.40 0.19 0.75 -355.78
## 5 0.59 3.63 0.19 0.72 -350.10
## 6 1.24 3.54 0.22 0.80 -354.80
## 7 1.51 3.03 0.19 0.78 -360.86
## 8 1.02 3.40 0.21 0.78 -354.09
## 9 1.37 4.09 0.21 0.76 -350.52
## 10 0.94 3.56 0.20 0.75 -353.15
## 11 0.94 4.29 0.21 0.71 -342.30
## 12 1.14 3.25 0.22 0.81 -356.49
## 13 1.66 3.62 0.17 0.74 -359.53
## 14 1.46 3.64 0.20 0.77 -356.67
## 15 1.02 3.13 0.16 0.72 -358.30
## 16 0.31 3.84 0.22 0.72 -341.27
## 17 1.54 3.59 0.19 0.77 -358.10
## 18 1.43 4.06 0.18 0.73 -355.63
## 19 1.53 3.73 0.21 0.78 -355.44
## 20 1.35 4.16 0.18 0.72 -352.88
## 21 1.53 3.08 0.18 0.77 -360.42
## 22 0.42 3.49 0.19 0.73 -350.55
## 23 1.42 3.85 0.21 0.79 -354.72
## 24 0.23 3.70 0.17 0.68 -350.52
## 25 1.02 3.36 0.20 0.77 -355.54
options(scipen=0)
Now let’s get an average across trees.
colMeans(high.ouwie_results)
## the[a] the[b] alpha sigsq log(L)
## 1.0776599 3.7045410 0.2024869 0.7556554 -352.7203710
As we predicted, the estimates of \(\theta\)), in particular, are biased towards each other -- but we also see that \(\alpha\) is biased downwards and \(\sigma^2\) upwards.
Finally, let’s compare this to phytools::fitmultiOU.
init<-setNames(
c(
mean(high.x[high.y==levels(high.y)[1]]),
mean(high.x[high.y==levels(high.y)[2]]),
log(2)/max(nodeHeights(phy)),
var(high.x)/max(nodeHeights(phy)),
fitMk(phy,high.y,model="ER")$rates),
c("theta[a]","theta[b]",
"alpha","sigsq","q[1]"))
high.fit_mou<-fitmultiOU(phy,high.x,high.y,model="ER",
levs=100,parallel=TRUE,ncores=ncores,root="mle",
trace=1,maxit=2000,init=init)
## iter the[a] the[b] alpha sigsq q[1] log(L)
## 0 1.9220 3.2096 0.0693 0.1930 0.4588 -622.5731
## 100 1.1219 4.2495 0.4753 0.6045 0.5548 -490.4077
## 200 -0.6367 4.8417 0.4181 0.2777 0.6418 -473.6047
## 300 -0.4362 4.8024 0.4288 0.2951 0.6848 -473.2549
## 400 -0.5098 5.0866 0.4541 0.1942 0.7713 -472.2389
## 500 -0.4734 5.0244 0.4641 0.2102 0.7708 -472.1732
## 600 -0.4698 5.0304 0.4678 0.2093 0.7854 -472.1688
## 700 -0.4622 5.0253 0.4686 0.2081 0.7829 -472.1673
## 732 -0.4625 5.0253 0.4682 0.2082 0.7824 -472.1673
## Done optimizing.
high.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.4625, 5.0253 ]
## alpha: 0.4682
## sigsq: 0.2082
##
## Estimated Q matrix:
## a b
## a -0.7823842 0.7823842
## b 0.7823842 -0.7823842
##
## Log-likelihood: -472.1673
##
## R thinks it has found the ML solution.
This is pretty astonishing & seems to make quite a convincing case for joint estimation.
Admittedly, this is a very small experiment, of course – however, it does seem to confirm the thesis that if our discrete character changes infrequently, then the traditional two-step process of generating stochastic character maps & then fitting a fixed regime model to each map should work just fine. On the other hand, it also generally affirms Revell (2013, which perhaps ought to be a bit better cited) that as the rate of transition in our character goes up, the two-step process becomes biased and tends to underestimate the difference between the regimes. (This also may affect the measurement of \(\alpha\), which is underestimated; and \(\sigma^2\), which is overestimated.)
This bias goes away when we jointly model \(x\) and \(y\) using our new finite-space approximation, which is very cool IMO.
OK, more on this later!
