New generic plot method for "ancr" object class with marginal ancestral state reconstructions

I’ve recently added an S3 plot method for the "ancr" object class produced by phytools generic marginal ancestral reconstruction function.

The idea is to very simply automate a model-fitting → comparison → ancestral state reconstruction → visualization workflow.

To illustrate it, I’m going to use a phylogeny and dataset of diet evolution in grunters from Davis et al. (2016). (Grunters are perciform fish. Don’t ask me for any more details!)

## load phytools
library(phytools)

packageVersion("phytools")

## [1] '1.8.4'

## load tree and data
terapontid.tre<-read.tree(file="http://www.phytools.org/Rbook/11/terapontidae.phy")
terapontid.dat<-read.csv(file="http://www.phytools.org/Rbook/11/terapontidae.csv",
  row.names=1)
diet<-setNames(c("carnivory","omnivory","herbivory")[terapontid.dat$Diet],
  rownames(terapontid.dat))

OK, now I’m going to fit four models: an equal-rates model (ER); an all-rates-different model (ARD); and an ordered model.

Instead of using fitMk, I’ll going to use phytools::fitHRM which fits a hidden-rates model, but can also be used to fit various standard Mk models (if ncat=1).

fitER<-fitHRM(terapontid.tre,diet,model="ER",ncat=1,pi="fitzjohn",
  parallel=TRUE,logscale=TRUE,ncores=10)

## 
## This is the design matrix of the fitted model.
## Does it make sense?
## 
##           carnivory herbivory omnivory
## carnivory         0         1        1
## herbivory         1         0        1
## omnivory          1         1        0
## 
## Opened cluster with 10 cores.
## Running optimization iterations in parallel.
## Please wait....

fitARD<-fitHRM(terapontid.tre,diet,model="ARD",ncat=1,pi="fitzjohn",
  parallel=TRUE,logscale=TRUE,ncores=10)

## 
## This is the design matrix of the fitted model.
## Does it make sense?
## 
##           carnivory herbivory omnivory
## carnivory         0         1        2
## herbivory         3         0        4
## omnivory          5         6        0
## 
## Opened cluster with 10 cores.
## Running optimization iterations in parallel.
## Please wait....

fitOrdered<-fitHRM(terapontid.tre,diet,ncat=1,pi="fitzjohn",
  umbral=TRUE,parallel=TRUE,logscale=TRUE,
  ordered=TRUE,order=c("carnivory","omnivory","herbivory"),
  ncores=10)

## 
## This is the design matrix of the fitted model.
## Does it make sense?
## 
##           carnivory herbivory omnivory
## carnivory         0         0        1
## herbivory         0         0        2
## omnivory          3         4        0
## 
## Opened cluster with 10 cores.
## Running optimization iterations in parallel.
## Please wait....

Now let’s plot our three fitted models.

par(mfrow=c(3,1))
plot(fitER,width=TRUE,color=TRUE,xlim=c(-1.8,1),ylim=c(-1,1),
  spacer=0.2,mar=c(0,0,1.1,0))
mtext("a) ER model",line=0,adj=0.1)
plot(as.Qmatrix(fitOrdered),width=TRUE,color=TRUE,
  xlim=c(-1.8,1),ylim=c(-1,1),spacer=0.2,mar=c(0,0,1.1,0))
mtext("b) Ordered model",line=0,adj=0.1)
plot(fitARD,width=TRUE,color=TRUE,xlim=c(-1.8,1),ylim=c(-1,1),
  spacer=0.2,mar=c(0,0,1.1,0))
mtext("b) ARD model",line=0,adj=0.1)

We can compare them using an anova call.

aov_grunters<-anova(fitER,fitOrdered,fitARD)

##               log(L) d.f.      AIC     weight
## object     -38.03035    1 78.06069 0.04666782
## fitOrdered -32.16432    4 72.32864 0.81981518
## fitARD     -31.97917    6 75.95834 0.13351699

ancr computes model-averaged marginal ancestral states.

grunters_ancr<-ancr(aov_grunters)
grunters_ancr

## Marginal ancestral state estimates:
##    carnivory herbivory omnivory
## 39  0.018704  0.966901 0.014396
## 40  0.019044  0.965501 0.015455
## 41  0.020037  0.950306 0.029658
## 42  0.022603  0.933903 0.043494
## 43  0.028993  0.460921 0.510086
## 44  0.311920  0.175454 0.512625
## ...
## 
## Log-likelihood = -33.035963

Finally, let’s plot our "ancr" object.

plot(grunters_ancr)

Not bad, right?