New example for multi-rate Brownian motion evolution method

I recently described a new penalized likelihood method to model variable-rate Brownian evolution on the tree.

To find out more about this method, you can either read a bioRxiv pre-print that I wrote about it, or check out my Evolution 2021 meeting presentation on YouTube.

I just added a new example to the documentation of the method in phytools. This example uses relative gape width evolution in the buccal morphology on a phylogeny of 28 species of centrarchids or sunfish. The example came out really nicely, so I thought I'd also show it on here.

First, let's load the data and tree:

data(sunfish.tree)
sunfish.tree

## 
## Phylogenetic tree with 28 tips and 27 internal nodes.
## 
## Tip labels:
##  Acantharchus_pomotis, Lepomis_gibbosus, Lepomis_microlophus, Lepomis_punctatus, Lepomis_miniatus, Lepomis_auritus, ...
## 
## The tree includes a mapped, 2-state discrete character
## with states:
##  non, pisc
## 
## Rooted; includes branch lengths.

data(sunfish.data)
head(sunfish.data)

##                      feeding.mode gape.width buccal.length
## Acantharchus_pomotis         pisc      0.114        -0.009
## Lepomis_gibbosus              non     -0.133        -0.009
## Lepomis_microlophus           non     -0.151         0.012
## Lepomis_punctatus             non     -0.103        -0.019
## Lepomis_miniatus              non     -0.134         0.001
## Lepomis_auritus               non     -0.222        -0.039

Next, let's pull out the character trait of interest – gape.width.

gw<-setNames(sunfish.data$gape.width,
    rownames(sunfish.data))
gw

##    Acantharchus_pomotis        Lepomis_gibbosus     Lepomis_microlophus       Lepomis_punctatus 
##                   0.114                  -0.133                  -0.151                  -0.103 
##        Lepomis_miniatus         Lepomis_auritus      Lepomis_marginatus       Lepomis_megalotis 
##                  -0.134                  -0.222                  -0.187                  -0.073 
##         Lepomis_humilis     Lepomis_macrochirus         Lepomis_gulosus     Lepomis_symmetricus 
##                   0.024                  -0.191                   0.131                   0.013 
##       Lepomis_cyanellus      Micropterus_coosae      Micropterus_notius     Micropterus_treculi 
##                  -0.002                   0.045                   0.097                   0.056 
##   Micropterus_salmoides  Micropterus_floridanus Micropterus_punctulatus    Micropterus_dolomieu 
##                   0.056                   0.096                   0.092                   0.035 
## Centrarchus_macropterus      Enneacantus_obesus       Pomoxis_annularis  Pomoxis_nigromaculatus 
##                  -0.007                   0.016                  -0.004                   0.105 
##  Archolites_interruptus    Ambloplites_ariommus   Ambloplites_rupestris   Ambloplites_cavifrons 
##                  -0.024                   0.135                   0.086                   0.130

Now I'll run our penalized likelihood method. I'm going to set the smoothing parameter, λ, to 1.0, but then I'll also try a couple of other values and see the effect.

fitBM<-multirateBM(sunfish.tree,gw,
    lambda=1)

## Beginning optimization....
## Optimization iteration 1. Using "L-BFGS-B" optimization method.
## Best (penalized) log-likelihood so far: 84.3134 
## Optimization iteration 2. Using "Nelder-Mead" optimization method.
## Best (penalized) log-likelihood so far: 84.3134 
## Optimization iteration 3. Using "BFGS" optimization method.
## Best (penalized) log-likelihood so far: 84.3135 
## Done optimization.

plot(fitBM,ftype="i",fsize=0.8,lwd=6,
    outline=TRUE)

This shows that the rate of evolution of relative gape width is much higher in the genus Lepomis compared to other parts of the tree. Let's try again with two other values of λ and see if this pattern is robust.

fitLambda0.1<-multirateBM(sunfish.tree,gw,
    lambda=0.1)

## Beginning optimization....
## Optimization iteration 1. Using "L-BFGS-B" optimization method.
## Best (penalized) log-likelihood so far: 39.261 
## Done optimization.

fitLambda10<-multirateBM(sunfish.tree,gw,
    lambda=10)

## Beginning optimization....
## Optimization iteration 1. Using "L-BFGS-B" optimization method.
## Best (penalized) log-likelihood so far: 564.272 
## Optimization iteration 2. Using "Nelder-Mead" optimization method.
## Best (penalized) log-likelihood so far: 564.272 
## Optimization iteration 3. Using "BFGS" optimization method.
## Best (penalized) log-likelihood so far: 564.273 
## Done optimization.

par(mfrow=c(1,2))
plot(fitLambda0.1,ftype="i",fsize=0.8,
    mar=c(1.1,1.1,2.1,1.1))
mtext(expression(paste("a) ",lambda," = 0.1")),
    line=0,adj=0)
plot(fitLambda10,ftype="i",fsize=0.8,
    mar=c(1.1,1.1,2.1,1.1))
mtext(expression(paste("b) ",lambda," = 10")),
    line=0,adj=0)

This is neat because the pattern shows up regardless of which value of λ we use.

Cool.