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.
