Monday, April 5, 2021

Updated plotting method for variable-rate penalized likelihood model in phytools

Inspired by comments from my friend & collaborator Lars Schmitz, I just pushed some updates to the variable rate Brownian motion penalized likelihood method that I blogged about here last week.

The idea of this approach, remember, is to fit a variable-rate Brownian motion model to our tree and data, in which the rate parameter itself (σ2) evolves from branch to branch in the tree (by a Brownian process on the log scale). Because this model would be non-identifiable if we allowed σ2 for the evolutionary rate to assume any value, we use a term (the probability of the rates under our model), and a user-specified penalty coefficient (λ), to determine how much σ2 should be permitted to change through time and among lineages. Larger values of the penalty term λ result in less changes; whereas lower values of λ permit σ2 to evolve from edge to edge on the phylogeny with relatively little cost.

In this update the main thing that I did was add much more user control of the color schema used for plotting our fitted object from the multirateBM function, as well as more options for the style of the plotted tree.

Here's a quick example using log-transformed home range size for mammals. These data come from Garland et al. (1992).

## Loading required package: ape
## Loading required package: maps
## Beginning optimization....
## Optimization iteration 1. Using "Nelder-Mead" optimization method.
## Optimization iteration 2. Using "BFGS" optimization method.
## Done optimization.
## Multi-rate Brownian model using multirateBM.
## Fitted rates:
##  U._maritimus U._arctos U._americanus N._narica P._lotor M._mephitis M._meles
##      0.036776  0.026262      0.029172  0.000897 0.000818    0.021497 0.023634
##  C._lupus     
##  0.546833 ....
## lambda penalty term: 1
## log-likelihood:  -79.564297 
## AIC:  355.128593 
## R thinks it has found a solution.

Now let's make our plot using the default settings. When I do it, I'll also return the plotted object to the user as follows:


plot of chunk unnamed-chunk-2

## Object of class "multirateBM_plot" containing:
## (1) A phylogenetic tree with 49 tips and 48 internal nodes.
## (2) A mapped set of Brownian evolution rates.

Now I'm going to re-plot this object in a different style:


plot of chunk unnamed-chunk-3

That's kind of neat, right.

Now let's change the color too. We'll do this using the phytools generic function setMap as follows:


plot of chunk unnamed-chunk-4

We can also use a totally arbitrary gradient if we want. Here's a grey-scale example:


plot of chunk unnamed-chunk-5

That's pretty cool.

For those interested in such things, some of the things that you'll notice in the source code of this update are: (1) a new object class ("multirateBM_plot") containing the plot object, rather than the fitted model; and (2) use of to call the plot methods (mostly, phytools::plotSimmap) internally. is a very handy function when we want to allow users to have a lot of flexibility in passing their own argument values to internally used functions, but some values for these arguments would break our function. Check it out!

No comments:

Post a Comment

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