Plotting a fitted Mk model with different arrow widths for different estimated transition rates

Yesterday I posted about graphing a fitted Mk model showing relative rate differences for the different transition types using different colors on a blue (slow) to red (fast) color gradient.

But what if we could also show rate differences using different arrow widths?

In fact, this is already implemented for the "fitPagel" object class. I thought it could be very useful for phytools users to add it for the standard Mk model as well.

Note that it doesn't matter whether we fit our model using phytools::fitMk or geiger::fitDiscrete (or another function, for that matter). As long as we can convert our fitted model to a "Qmatrix" object, we're good to go!

Here's an example, using the same Anolis dataset as I employed yesterday.

We can start by loading our data and fitting the model using fitMk.

## load phytools
library(phytools)
packageVersion("phytools")

## [1] '1.1.16'

## load data
anole.tree<-read.tree(file="http://www.phytools.org/Rbook/1/Anolis.tre")
anole.ecomorphs<-read.csv(file="http://www.phytools.org/Rbook/1/ecomorph.csv",
    row.names=1)
ecomorph.tree<-drop.tip(anole.tree,
    geiger::name.check(anole.tree,anole.ecomorphs)$tree_not_data)
ecomorph<-setNames(anole.ecomorphs$ecomorph,row.names(anole.ecomorphs))
## fit ARD model
fitARD<-fitMk(ecomorph.tree,ecomorph,model="ARD",pi="fitzjohn")
fitARD

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##          CG       GB        TC        TG       Tr        Tw
## CG 0.000000 0.000000  0.000000  0.000000 0.000000  0.000000
## GB 0.000000 0.000000  0.000000  0.000000 0.000000  0.000000
## TC 0.278023 0.000000 -0.495184  0.000000 0.217161  0.000000
## TG 0.078225 0.332274  0.090038 -0.991693 0.072734  0.418422
## Tr 0.000000 0.000000  0.000000  0.000000 0.000000  0.000000
## Tw 0.112733 0.334595  0.642925  0.000000 0.000000 -1.090252
## 
## Fitted (or set) value of pi:
## CG GB TC TG Tr Tw 
##  0  0  0  1  0  0 
## due to treating the root prior as (a) nuisance.
## 
## Log-likelihood: -68.226274 
## 
## Optimization method used was "nlminb"

Next, let's plot our fitted model. To turn on variable arrow widths, we'll set the argument width=TRUE.

plot(fitARD,width=TRUE)
title("Fitted transition rates\nunder ARD model",
    font.main=3,line=-1)

That's not bad.

Now let's try both color=TRUE and width=TRUE.

plot(fitARD,width=TRUE,color=TRUE)
title("Fitted transition rates\nunder ARD model",
    font.main=3,line=-1)

Neat.

Lastly, since the rates along edges make the plot a bit crowded, and are superflous to the information of the color legend anyway, let's eliminate them.

This can be done using the argument text=FALSE. We'll also set show.zeros=FALSE to plot only arrows with non-zero* transition rates. (*To a certain tolerance level that can be controlled by the user.)

plot(fitARD,width=TRUE,color=TRUE,text=FALSE,
    show.zeros=FALSE)
title("Fitted transition rates\nunder ARD model",
    font.main=3,line=-1)

Awesome.

I purposefully allow the arrow widths to be larger for this plot since we don't need to worry about overwriting the rates.

Lastly, let's try a different example using a dataset from my book with Luke Harmon.

In this case, our data consist of a phylogeny & diet trait dataset for terapontid fishes (grunters) from a publication by Davis et al (2016).

terapontid.tree<-read.tree(file="http://www.phytools.org/Rbook/11/terapontidae.phy")
terapontid.data<-read.csv(file="http://www.phytools.org/Rbook/11/terapontidae.csv",
    row.names=1)
diet<-setNames(c("piscivorous","omnivorous","herbivorous")[terapontid.data$Diet],
    rownames(terapontid.data))
fitDiet<-fitMk(terapontid.tree,diet,model="ARD",pi="fitzjohn")
fitDiet

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##             herbivorous omnivorous piscivorous
## herbivorous   -0.066364   0.066364    0.000000
## omnivorous     0.000000  -0.092228    0.092228
## piscivorous    0.004065   0.000000   -0.004065
## 
## Fitted (or set) value of pi:
## herbivorous  omnivorous piscivorous 
##    0.903924    0.048445    0.047631 
## due to treating the root prior as (a) nuisance.
## 
## Log-likelihood: -31.979172 
## 
## Optimization method used was "nlminb"

plot(fitDiet,width=TRUE,color=TRUE,text=FALSE,cex.traits=0.8,spacer=0.15)
title(main="Estimated diet transition rates\nfor grunter fishes",
    font.main=3,line=-1)

Without imposing any particular structure, the fitted model clearly indicates a directional evolutionary process from herbivorous → omnivorous → piscivorous, which kind of makes sense, I guess. (I'm not an ichthyologist!)

That's all!