Over the past few days I’ve posted about, and gave an R function for, a (so far as I know) novel method to modeling discrete-state-dependent rate heterogeneity in a continuous character on a phylogeny based on the discrete approximation of Boucher & DÃ©mery (2016).

I’ve written more about Boucher & DÃ©mery (2016) in multiple earlier posts to this blog, which I would encourage readers to check out (1, 2, 3, 4, 5, 6, 7, 8).

The new function (`fitmultiBM`

) is now in *phytools* on GitHub. Please use with considerable caution. I’m still working out the bugs!

The idea of this method & function is to jointly model the evolution of a discrete character (evolving via a standard M*k* stochastic process) and a continuous trait (whose rate of evolution \(\sigma^2\) is, in turn, influenced by the discrete character).

The standard way to investigate this type of hypothesis in phylogenetic comparative biology has been to *first* map the discrete character history on the tree (using, for example, stochastic character mapping), and *then* fit a regime-dependent model following O’Meara et al. (2006).

The disadvantage of this two-step process is that by first mapping regimes onto the tree in a way that considers only the discrete trait, we’re ignoring any information that our continuous character might contain about the discrete character history. Averaging across a set of such maps can in turn result in a biased underestimate of the difference in rate between regimes (e.g., see Revell 2013).

The discrete approximation allows us to simultaneously model our discrete and continuous characters. I predict that this will provide both more accurate estimates of the state-dependent rate variation in our continuous character, as well as a more accurate reconstruction of the evolutionary process for the discrete trait. I hope to explore this in future posts.

Before getting to that, though, I thought it might be neat to demo the method for a couple of empirical cases. These have been selected *entirely* due to the convenience of the data as both involve datasets that Luke Harmon & I used for our 2022 book with Princeton University Press.

We can begin, of course, by loading *phytools*.

```
library(phytools)
```

```
packageVersion("phytools")
```

```
## [1] '2.2.12'
```

(To follow along, you’ll need this version of *phytools* or later.)

Our first example involves a phylogeny and dataset for elopomorph eels from Collar et al. (2014). I’m going to start off by reading these two files directly from the book website.

```
eel_tree<-read.tree(
file="http://www.phytools.org/Rbook/8/elopomorph.tre")
eel_data<-read.csv(
file="http://www.phytools.org/Rbook/8/elopomorph.csv",
row.names=1,stringsAsFactors=TRUE)
```

Our phylogeny contains 61 tips and 60 nodes.

```
eel_tree
```

```
##
## Phylogenetic tree with 61 tips and 60 internal nodes.
##
## Tip labels:
## Moringua_edwardsi, Kaupichthys_nuchalis, Gorgasia_taiwanensis, Heteroconger_hassi, Venefica_proboscidea, Anguilla_rostrata, ...
##
## Rooted; includes branch lengths.
```

Our data set has one discrete character (feeding mode: i.e., whether each eel species is a suction feeder or a biter) and a second continuous trait, overall body size. Let’s pull these out. Please keep in mind that I have no *a priori* reason to imagine that feeding mode influences the rate of body size evolution in eels, but these are the data we have!

```
head(eel_data)
```

```
## feed_mode Max_TL_cm
## Albula_vulpes suction 104
## Anguilla_anguilla suction 50
## Anguilla_bicolor suction 120
## Anguilla_japonica suction 150
## Anguilla_rostrata suction 152
## Ariosoma_anago suction 60
```

```
feed_mode<-setNames(eel_data$feed_mode,
rownames(eel_data))
lnTL<-setNames(log(eel_data$Max_TL_cm),
rownames(eel_data))
```

Let’s graph our data.

```
par(mfrow=c(1,2))
FMODE<-to.matrix(feed_mode[eel_tree$tip.label],
levels(feed_mode))
plotTree(eel_tree,mar=c(5.1,4.1,1.1,1.1),
fsize=0.4,ftype="i")
tiplabels(pie=FMODE,piecol=hcl.colors(2),cex=0.4)
par(mar=c(5.1,4.1,1.1,1.1))
phenogram(eel_tree,lnTL,spread.cost=c(1,0),ftype="i",
fsize=0.4,ylab="log(max length)",quiet=TRUE,las=1,
cex.axis=0.7)
tiplabels(pie=FMODE,piecol=hcl.colors(2),cex=0.4)
legend("topleft",levels(feed_mode),pch=21,
pt.bg=hcl.colors(2),bty="n",cex=0.8,pt.cex=1.5)
```

We’re now ready to fit our discrete-state-dependent multi-rate Brownian motion evolution model using the discrete approximation of Boucher & DÃ©mery (2016).

As I do that, one small feature that I’ll demo here is a hidden option called `plot_model`

. When we set `plot_model=TRUE`

, the function graphs the structure of the discretized model we’re fitting. I built this mostly for the purposes of debugging, but it’s generally useful for those readers of this blog interested in how these models are structured. The way to interpret the visualization is as a very large \(k \cdot n \times k \cdot n\) in size (for \(k\) levels of our discrete trait and \(n\) bins for our continuous character) model design matrix, in which each different color is a different estimated parameter in the model – either a different rate of evolution for our continuous character (on the sub- and superdiagonals) or a rate of transition between discrete states (elsewhere).

```
eel_sdmodel<-fitmultiBM(eel_tree,lnTL,feed_mode,
model="ARD",parallel=TRUE,levs=100,
plot_model=TRUE)
```

Let’s print our results.

```
eel_sdmodel
```

```
## Object of class "fitmultiBM" based on
## a discretization with k = 100 levels.
##
## Fitted multi-rate BM model parameters:
## levels: [ bite, suction ]
## sigsq: [ 0.0098, 0.0214 ]
## x0: 4.1992
##
## Estimated Q matrix:
## bite suction
## bite -0.013215826 0.013215826
## suction 0.008103207 -0.008103207
##
## Log-likelihood: -100.7688
##
## R thinks it has found the ML solution.
```

We can see from the fitted model that the ML estimated \(\sigma^2\) rate of log(maximum body size) evolution is about twice as high in suction feeding compared to biting eels. What’s more difficult to say, however, is whether this model has significantly better explanatory power than one in which there was no difference in the rate of evolution of size between our two groups.

We can (now) fit this model in `fitmultiBM`

by setting the hidden option `null_model=TRUE`

. Let’s try it.

```
eel_nullmodel<-fitmultiBM(eel_tree,lnTL,feed_mode,
model="ARD",parallel=TRUE,levs=100,
null_model=TRUE)
```

```
eel_nullmodel
```

```
## Object of class "fitmultiBM" based on
## a discretization with k = 100 levels.
##
## Fitted multi-rate BM model parameters:
## levels: [ bite, suction ]
## sigsq: [ 0.0142, 0.0142 ]
## x0: 4.1992
##
## Estimated Q matrix:
## bite suction
## bite -0.01488082 0.01488082
## suction 0.01357827 -0.01357827
##
## Log-likelihood: -102.2585
##
## R thinks it has found the ML solution.
```

Note, simply as an aside, that both our parameter estimates and log-likelihood under the null model should very closely match those estimated independent (e.g., using `geiger::fitContinuous`

and `fitMk`

), as follows.

```
geiger::fitContinuous(eel_tree,lnTL)$opt$sigsq
```

```
## [1] 0.01444782
```

```
phytools::fitMk(eel_tree,feed_mode,
model="ARD",pi="mle")$rates
```

```
## [1] 0.01357832 0.01488083
```

(One issue I haven’t worked out is apparent above. To get the correct log-likelihood for the continuous data I have to set the root to its MLE, but for the moment I’m doing that for both the discrete *and* continuous trait at once, which I’m not sure is desirable. That’s why I needed to set `pi="mle"`

in `fitMk`

to get the the “correct” parameter estimates. I should be able to figure out how to fix this!)

```
logLik(geiger::fitContinuous(eel_tree,lnTL))+
logLik(fitMk(eel_tree,feed_mode,model="ARD",
pi="mle"))
```

```
## [1] -102.2191
## attr(,"df")
## [1] 2
```

To compare the two models, let’s use a generic `anova`

call as follows.

```
anova(eel_nullmodel,eel_sdmodel)
```

```
## log(L) d.f. AIC weight
## eel_nullmodel -102.2585 4 212.5170 0.3799782
## eel_sdmodel -100.7688 5 211.5377 0.6200218
```

Since these two models are nested, we could also undertake a likelihood-ratio test. Let’s do that using the CRAN package *lmtest*.

```
lmtest::lrtest(eel_nullmodel,eel_sdmodel)
```

```
## Likelihood ratio test
##
## Model 1: eel_nullmodel
## Model 2: eel_sdmodel
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 -102.26
## 2 5 -100.77 1 2.9793 0.08434 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

Even though AIC *slightly* (though non-significantly, according to our likelihood ratio test) favors the multi-rate model, this is pretty weak evidence to support feeding mode dependent evolution of body size in eels!

There are a couple of other models that we might like to consider.

One possibility is that the rate of evolution of body size *is* heterogeneous, so our null model is *indeed* wrong – but this rate heterogeneity has nothing at all to do with feeding mode.

To examine this possibility we might fit a hidden-rate model in which our discrete character evolves but only the hidden character influences body size evolution. I’m going to do that by keeping `null_model=TRUE`

but updating `ncat=2`

. This is a “null model” in the sense that, as we’ll see in a sec, under this model the rate of evolution of log(maximum length) varies throughout the tree, but not as a function of our discrete character: feeding mode.

```
eel_nullmodel.hrm<-fitmultiBM(eel_tree,lnTL,feed_mode,
model="ARD",ncat=2,parallel=TRUE,levs=100,
null_model=TRUE)
```

```
eel_nullmodel.hrm
```

```
## Object of class "fitmultiBM" based on
## a discretization with k = 100 levels.
##
## Fitted multi-rate BM model parameters:
## levels: [ bite, bite*, suction, suction* ]
## sigsq: [ 0.0304, 0.0026, 0.0304, 0.0026 ]
## x0: 4.368
##
## Estimated Q matrix:
## bite bite* suction suction*
## bite -0.04059342 0.02571261 0.01488080 0.00000000
## bite* 0.01968139 -0.03456219 0.00000000 0.01488080
## suction 0.01357809 0.00000000 -0.03929071 0.02571261
## suction* 0.00000000 0.01357809 0.01968139 -0.03325948
##
## Log-likelihood: -96.8294
##
## R thinks it has found the ML solution.
```

Lastly, we can also consider a hidden-rate model in which both feeding mode *and* our hidden character influence the rate of evolution of eel body size.

```
eel_sdmodel.hrm<-fitmultiBM(eel_tree,lnTL,feed_mode,
model="ARD",ncat=2,parallel=TRUE,levs=100)
```

```
eel_sdmodel.hrm
```

```
## Object of class "fitmultiBM" based on
## a discretization with k = 100 levels.
##
## Fitted multi-rate BM model parameters:
## levels: [ bite, bite*, suction, suction* ]
## sigsq: [ 0.0231, 0.0011, 0.0399, 0.0044 ]
## x0: 4.3117
##
## Estimated Q matrix:
## bite bite* suction suction*
## bite -0.03540772 0.02098248 0.01442524 0.00000000
## bite* 0.01817748 -0.03260272 0.00000000 0.01442524
## suction 0.01189854 0.00000000 -0.03288102 0.02098248
## suction* 0.00000000 0.01189854 0.01817748 -0.03007601
##
## Log-likelihood: -95.6653
##
## R thinks it has found the ML solution.
```

As a word of warning, as we progress up through these models they become both more parameter complex *and* also involve larger & larger model matrices for the discrete approximation. That means to reiterate the analyses above (depending on the speed of your comput
er) you may need to start it running, go have dinner, and then come back to it!

Here’s a comparison of all four models, generally as they increase in model complexity.

```
anova(eel_nullmodel,eel_sdmodel,eel_nullmodel.hrm,
eel_sdmodel.hrm)
```

```
## log(L) d.f. AIC weight
## eel_nullmodel -102.25848 4 212.5170 0.05291132
## eel_sdmodel -100.76884 5 211.5377 0.08633699
## eel_nullmodel.hrm -96.82939 7 207.6588 0.60046674
## eel_sdmodel.hrm -95.66532 9 209.3306 0.26028495
```

This shows us that the hidden-state heterogeneous rate model *without* state-dependent body size evolution is the best supported, accounting for its parameter complexity. This would indicate that though there may be rate heterogeneity of body size, it is not well-explained by our hypothesized discrete character.

For our second example, I’m going to use diel activity pattern and skull length in primates with data from Kirk & Kay (2004).

These data and tree are even easy to get our hands on than for our previous example as they now come packaged with *phytools*!

```
data(primate.tree)
primate.tree
```

```
##
## Phylogenetic tree with 90 tips and 89 internal nodes.
##
## Tip labels:
## Allenopithecus_nigroviridis, Cercopithecus_mitis, Cercopithecus_cephus, Cercopithecus_petaurista, Erythrocebus_patas, Chlorocebus_aethiops, ...
##
## Rooted; includes branch lengths.
```

```
data(primate.data)
lnSkull<-setNames(log(primate.data$Skull_length),
rownames(primate.data))
activity<-setNames(primate.data$Activity_pattern,
rownames(primate.data))
```

Once again, let’s graph our data.

```
par(mfrow=c(1,2))
ACTIVITY<-to.matrix(activity[primate.tree$tip.label],
levels(activity))
plotTree(primate.tree,mar=c(5.1,4.1,1.1,1.1),
fsize=0.3,ftype="i")
cols<-setNames(hcl.colors(n=3)[c(2,3,1)],
colnames(ACTIVITY))
tiplabels(pie=ACTIVITY,piecol=cols,cex=0.4)
par(mar=c(5.1,4.1,1.1,1.1))
phenogram(primate.tree,lnSkull,spread.cost=c(1,0),ftype="i",
fsize=0.3,ylab="log(skull length)",quiet=TRUE,las=1,
cex.axis=0.7)
tiplabels(pie=ACTIVITY,piecol=cols,cex=0.4)
legend("topleft",colnames(ACTIVITY),pch=21,
pt.bg=cols,bty="n",cex=0.8,pt.cex=1.5)
```