Update to fitPagel binary trait evolution model plotting method to facilitate its use custom multi-panel figures

During my workshop & seminar visit to Lund University in Sweden this week, a phytools user asked me about combining the subplots of a graphed fitted Pagel ’94 (correlational) binary trait evolution model with a graphed tree.

One limitation of the existing implementation is that it used par(mfrow=c(1,2)) to create subplots, and thus was not very flexible when wanting to generate a multi-panel figure differing from the default.

Luckily, this was easy enough to fix – and can already be obtained by installing phytools from it’s GitHub page, e.g., using the devtools or remotes packages.

Having done that, let’s try it out.

To start, I’ll load phytools and a demo dataset from Benun Sutton and Wilson (2019).

library(phytools)
packageVersion("phytools")

## [1] '1.5.13'

data(bonyfish.tree)
data(bonyfish.data)

Next, I’ll pull out my two binary traits for using in phytools::fitPagel.

## extract discrete characters
spawning_mode<-setNames(bonyfish.data$spawning_mode,
    rownames(bonyfish.data))
paternal_care<-setNames(bonyfish.data$paternal_care,
    rownames(bonyfish.data))

Now I’ll fit a total of four models – an independent model, and three dependent models: one in which parent care depends on spawning mode (but not the reverse), a second in which spawning mode depends on parental care (but not the reverse), and a third in which the two traits evolve interdependently.

bonyfish.depx<-fitPagel(bonyfish.tree,paternal_care,
	spawning_mode,dep.var="x",pi="fitzjohn")
bonyfish.depy<-fitPagel(bonyfish.tree,paternal_care,
	spawning_mode,dep.var="y",pi="fitzjohn")
bonyfish.full<-fitPagel(bonyfish.tree,paternal_care,
	spawning_mode,pi="fitzjohn")
anova(bonyfish.depx,bonyfish.depy,bonyfish.full)

##                  log(L) d.f.      AIC     weight
## independent   -60.78126    4 129.5625 0.03856507
## bonyfish.depx -57.53577    6 127.0715 0.13399940
## bonyfish.depy -56.05585    6 124.1117 0.58860816
## bonyfish.full -54.95787    8 125.9157 0.23882738

This shows us that there is little support for the independent model – but the three dependent models have substantial model weights (although the latter two dependent models seem to be quite a bit better than the first).

Finally, let’s create a multipanel figure in which we show: the phylogeny and trait data (using phytools::plotTree.datamatrix); and then three figure panels showing each one of the fitted dependent models!

To subdivide our plot into multiple panels, we’ll use the base R graphics function layout.

layout(matrix(c(1,1,1,2,3,4),3,2))
object<-plotTree.datamatrix(bonyfish.tree,bonyfish.data,
	fsize=0.5,yexp=1,header=FALSE,
	palettes=c("YlOrRd","PuBuGn"))

## Loading required package: RColorBrewer

leg<-legend(x="topleft",names(object$colors$spawning_mode),
	cex=0.7,pch=22,pt.bg=object$colors$spawning_mode,
	pt.cex=1.5,bty="n",title="spawning mode")
## add a second legend for trait 2
leg<-legend(x=leg$rect$left+4.7,y=leg$rect$top-leg$rect$h,
	names(object$colors$paternal_care),cex=0.7,
	pch=22,pt.bg=object$colors$paternal_care,pt.cex=1.5,
	bty="n",title="paternal care")
plot(bonyfish.depx,show="dependent",main="")
mtext(paste("a) parental care depends on spawning mode",
	"\n  (but not the reverse), AIC =",
	round(bonyfish.depx$dependent.AIC,2)),
	adj=0,cex=0.7,line=-0.25)
plot(bonyfish.depy,show="dependent",main="")
mtext(paste("b) spawning mode depends on parental care",
	"\n  (but not the reverse), AIC =",
	round(bonyfish.depy$dependent.AIC,2)),
	adj=0,cex=0.7,line=-0.25)
plot(bonyfish.full,show="dependent",main="")
mtext(paste("c) interdependent trait evolution model",
	"\n  for spawning mode & parental care, AIC =",
	round(bonyfish.full$dependent.AIC,2)),
	adj=0,cex=0.7,line=-0.25)

Neat!