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?
