A couple of days ago a colleague asked me about running multiple optimization iterations of the new *phytools* function `fitmultiBM`

, as well as some other stuff pertaining to how it can be used.

I thought I would create a quick demo of this and some other things. I’m still working on some of the details of this method and it’s implementation. Nonetheless, let’s go….

Readers following on should install *phytools* 2.3-2 from GitHub as I intend to use some features of `fitmultiBM`

absent from the current CRAN version.

```
library(phytools)
```

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

```
## [1] '2.3.2'
```

For this example, I’m going to use data for the totally real fish group called ‘grunters’ (Terapontidae). The phylogeny and trait data for this demo come from Davis et al. (2016) and can be downloaded from the website for our book.

```
grunter_tree<-read.tree(
file="http://www.phytools.org/Rbook/11/terapontidae.phy")
grunter_tree
```

```
##
## Phylogenetic tree with 38 tips and 37 internal nodes.
##
## Tip labels:
## Terapon_jarbua, Pelsartia_humeralis, Mesopristes_cancellatus, Mesopristes_argenteus, Rhyncopelates_oxyrhynchus, Terapon_puta, ...
##
## Rooted; includes branch lengths.
```

```
grunter_data<-read.csv(
file="http://www.phytools.org/Rbook/11/terapontidae.csv",
row.names=1)
head(grunter_data)
```

```
## Diet Animal Plant PC1 PC2
## Amniataba_affinis 2 0.2818512 -0.2818512 0.2913522 -0.21526981
## Amniataba_caudavittatus 1 2.5866893 -2.5866893 -0.6336029 0.15468468
## Amniataba_percoides 2 -0.2330491 0.2330491 0.1912102 0.07985197
## Bidyanus_bidyanus 3 -2.1972246 2.1972246 0.4823617 -0.78346644
## Bidyanus_welchi 2 0.8472979 -0.8472979 0.2891397 -0.77832681
## Hannia_greenwayi_1 2 1.1970524 -1.1970524 -0.1860913 -0.29919632
```

In our dataset I can see that our discrete character, `Diet`

, is coded as a numeric variable. We could leave it like this, but I know from the study that state `1`

is carnivore, `2`

is omnivore, and `3`

is herbivore, so let’s re-code our character to match these states.

```
diet_cat<-setNames(grunter_data$Diet,rownames(grunter_data))
diet<-setNames(
as.factor(c("carnivore","omnivore","herbivore")[diet_cat]),
names(diet_cat))
```

Here’s our discrete trait mapped onto the tips of the tree.

```
plotTree(grunter_tree,direction="upwards",fsize=0.8,
ftype="i",offset=3)
tiplabels(pie=to.matrix(diet[grunter_tree$tip.label],
levels(diet)),cex=0.5,offset=1.5,
piecol=hcl.colors(n=3,palette="plasma")[c(2,1,3)])
legend("bottomright",levels(diet),pch=21,
pt.bg=hcl.colors(n=3,palette="plasma")[c(2,1,3)],
cex=0.8,pt.cex=2,inset=0.03,bty="n",horiz=TRUE)
```

The continuous character we’re going to focus on for this analysis is `PC1`

in the data frame. This is a principle component from an analysis of dental morphology. (For more details, see Davis et al. 2016)

```
pc1<-setNames(grunter_data$PC1,rownames(grunter_data))
head(pc1)
```

```
## Amniataba_affinis Amniataba_caudavittatus Amniataba_percoides Bidyanus_bidyanus
## 0.2913522 -0.6336029 0.1912102 0.4823617
## Bidyanus_welchi Hannia_greenwayi_1
## 0.2891397 -0.1860913
```

Let’s proceed to fit our multi-rate model.

In this case, I’m going to begin by assuming that evolution is *ordered* such that carnivore \(\leftrightarrow\) omnivore \(\leftrightarrow\) herbivore, but in which the back-and-forth transition rates between adjacent levels are allowed to be asymmetric and different between each pair of states. (This is the thing that we couldn’t have done with the CRAN version of *phytools*!)

```
## design matrix for discrete character model
model<-matrix(c(0,0,1,0,0,2,3,4,0),3,3,
dimnames=list(levels(diet),levels(diet)))
model
```

```
## carnivore herbivore omnivore
## carnivore 0 0 3
## herbivore 0 0 4
## omnivore 1 2 0
```

```
pc1_multirate<-fitmultiBM(grunter_tree,pc1,diet,
model=model,parallel=TRUE,levs=100,plot_model=TRUE)
```

```
## This is the design of your custom discrete-trait model:
## carnivore herbivore omnivore
## carnivore 0 0 3
## herbivore 0 0 4
## omnivore 1 2 0
```

(My apologies for the random palette. To my colorblind eye, the `'carnivore'`

and `'herbivore'`

continuous character rate shades look virtually identical, but I assure you they are *different* in the fitted model! My non-colorblind wife tells me they *do* look different to her eye.)

```
pc1_multirate
```

```
## Object of class "fitmultiBM" based on
## a discretization with k = 100 levels.
##
## Fitted multi-rate BM model parameters:
## levels: [ carnivore, herbivore, omnivore ]
## sigsq: [ 0.0035, 0.0293, 0.1598 ]
## x0: -0.406
##
## Estimated Q matrix:
## carnivore herbivore omnivore
## carnivore -0.03639316 0.00000000 0.03639316
## herbivore 0.00000000 -0.05789947 0.05789947
## omnivore 0.07757236 0.05445299 -0.13202536
##
## Log-likelihood: -61.0404
##
## R thinks it has found the ML solution.
```

According to this fitted model, the condition `omnivore`

has the highest \(\sigma^2\) rate of trait evolution for `pc1`

, `carnivore`

is the lowest, and `herbivore`

is in between.

Now, `fitmultiBM`

can be somewhat difficult to optimize – but does not run multiple optimization iterations by default. It *does*, on the other hand, use random initial values for optimization (unlike `fitMk`

), so undertaking multiple optimization is as simple as writing a `for`

loop! (One might think we could use the *foreach* package to parallelize multiple optimization iterations across cores as I have shown elsewhere on this blog; however, we *cannot* because we are parallelizing computation of the likelihood in each optimization iteration already!)

```
fits<-list()
niter<-10
for(i in 1:niter){
capture.output(
fits[[i]]<-fitmultiBM(grunter_tree,pc1,diet,
model=model,parallel=TRUE,levs=100)
)
}
```

(Calling `capture.output`

simply has the effect of suppressing a print-out of our design matrix for every iteration of optimization! If we don’t care about that, we can simply remove it.)

Now let’s get the log-likelihoods for each fitted model.

```
lnL<-sapply(fits,logLik)
lnL
```

```
## [1] -61.04035 -61.04035 -61.04035 -61.04035 -61.04035 -61.04035 -61.04035 -61.04035 -61.04035
## [10] -61.04035
```

Depending on our fortune, we may find that these are all the same – or that some are higher than others.

We can pull out the best-fitting model from the set as follows.

```
best_multirate<-fits[[which.max(lnL)[1]]]
```

(The simple reason I use `which.max(...)[1]`

instead of `which.max(...)`

is in case we’ve converged on the same solution with *precisely* the same log-likelihood in more than one optimization iteration. Indeed, we should *hope* that this is the case!)

```
best_multirate
```

```
## Object of class "fitmultiBM" based on
## a discretization with k = 100 levels.
##
## Fitted multi-rate BM model parameters:
## levels: [ carnivore, herbivore, omnivore ]
## sigsq: [ 0.0035, 0.0293, 0.1598 ]
## x0: -0.406
##
## Estimated Q matrix:
## carnivore herbivore omnivore
## carnivore -0.03639381 0.00000000 0.03639381
## herbivore 0.00000000 -0.05789971 0.05789971
## omnivore 0.07757370 0.05445375 -0.13202745
##
## Log-likelihood: -61.0404
##
## R thinks it has found the ML solution.
```

For this experiment, most likely we’ve found that this fitted model is the same as the one we obtained earlier. Just in case it’s *better* let’s check & replace it if that’s the case.

```
pc1_multirate<-
if(logLik(best_multirate)>logLik(pc1_multirate)) best_multirate else
pc1_multirate
```

Quite often, and for more complex models and larger trees, it *will* take multiple optimization iterations to converge on the correct ML solution!

Lastly, at least for now, let’s compare this fitted model to our standard null model in which the rate of evolution of `pc1`

is *the same* for each of the three levels of our discrete trait.

```
pc1_nullmodel<-fitmultiBM(grunter_tree,pc1,diet,
model=model,parallel=TRUE,levs=100,plot_model=TRUE,
null_model=TRUE)
```

```
## This is the design of your custom discrete-trait model:
## carnivore herbivore omnivore
## carnivore 0 0 3
## herbivore 0 0 4
## omnivore 1 2 0
```

(The difference between this model structure and our multi-rate model is that the super- and sub-diagonals no longer differ by discrete state. Do you see it?)
```
pc1_nullmodel
```

```
## Object of class "fitmultiBM" based on
## a discretization with k = 100 levels.
##
## Fitted multi-rate BM model parameters:
## levels: [ carnivore, herbivore, omnivore ]
## sigsq: [ 0.0408, 0.0408, 0.0408 ]
## x0: -0.344
##
## Estimated Q matrix:
## carnivore herbivore omnivore
## carnivore -0.02005748 0.000000000 0.02005748
## herbivore 0.00000000 -0.067402676 0.06740268
## omnivore 0.10521478 0.008055272 -0.11327005
##
## Log-likelihood: -67.9894
##
## R thinks it has found the ML solution.
```

We can compare these two models, one against the other as follows.

```
anova(pc1_nullmodel,pc1_multirate)
```

```
## log(L) d.f. AIC weight
## pc1_nullmodel -67.98944 6 147.9789 0.007039976
## pc1_multirate -61.04035 8 138.0807 0.992960024
```

This tells us that there’s a lot more support in these data for a scenario in which \(\sigma^2\) of our continuous trait varies as a function of our discrete trait than for one in which the rate is constant throughout the tree.

Note that this null model is *the same* and should (asympotically, as we increase `levs`

) have the *same likelihood*, as a standard M*k* and constant rate Brownian model for our discrete and continuous traits, respectively.

Let’s see that this is the case.

```
mk_fit<-fitMk(grunter_tree,diet,model=model,
pi="mle")
mk_fit
```

```
## Object of class "fitMk".
##
## Fitted (or set) value of Q:
## carnivore herbivore omnivore
## carnivore -0.020063 0.000000 0.020063
## herbivore 0.000000 -0.067403 0.067403
## omnivore 0.105219 0.008056 -0.113275
##
## Fitted (or set) value of pi:
## carnivore herbivore omnivore
## 0 1 0
## due to treating the root prior as (a) it's MLE.
##
## Log-likelihood: -32.16097
##
## Optimization method used was "nlminb"
##
## R thinks it has found the ML solution.
```

```
bm_fit<-geiger::fitContinuous(grunter_tree,pc1)
bm_fit
```

```
## GEIGER-fitted comparative model of continuous data
## fitted 'BM' model parameters:
## sigsq = 0.044206
## z0 = -0.360878
##
## model summary:
## log-likelihood = -36.372527
## AIC = 76.745055
## AICc = 77.087912
## 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
```

We should see that the parameter estimates are almost exactly the same as in our `pc1_nullmodel`

object. Now let's check the likelihood.

```
## total likelihood
logLik(mk_fit)+logLik(bm_fit)
```

```
## [1] -68.5335
## attr(,"df")
## [1] 4
```

To the extent that this log-likelihood is *not* the same as we obtained for our null model using `fitmultiBM`

, this difference should asymptotically disappear as we increase `levs`

for the discrete approximation. This can get a bit computationally heavy, but let’s just try `levs=300`

and see what happens.

```
fitmultiBM(grunter_tree,pc1,diet,
model=model,parallel=TRUE,levs=300,
null_model=TRUE)
```

```
## This is the design of your custom discrete-trait model:
## carnivore herbivore omnivore
## carnivore 0 0 3
## herbivore 0 0 4
## omnivore 1 2 0
```

```
## Object of class "fitmultiBM" based on
## a discretization with k = 300 levels.
##
## Fitted multi-rate BM model parameters:
## levels: [ carnivore, herbivore, omnivore ]
## sigsq: [ 0.0436, 0.0436, 0.0436 ]
## x0: -0.3647
##
## Estimated Q matrix:
## carnivore herbivore omnivore
## carnivore -0.02006333 0.000000000 0.02006333
## herbivore 0.00000000 -0.067402967 0.06740297
## omnivore 0.10521778 0.008056971 -0.11327475
##
## Log-likelihood: -68.3704
##
## R thinks it has found the ML solution.
```

Now they're even closer, and this trend will continue as we increase `levs`

. As a valuable sidenote, even with `levs=300`

I left this all day (no exaggeration) to run!

Anyway, this is neat. I’ll leave it here for now, but will definitely be circling back to this method and function soon! Follow this space if you’re interested.

## No comments:

## Post a Comment

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.