Plotting facing "arc" or "fan" style trees with phytools

I’m just playing around with some stuff on a Saturday, but here’s a quick tutorial on plotting facing "arc" or "fan" style trees with opposite orientation. (No reason. It just looks cool.)

For the example, I’ll use a phylogeny & dataset for liolaemid lizards from Esquerre et al. (2019).

library(phytools)
data("liolaemid.tree")
data("liolaemid.data")

Let’s generate something to graph on the tree: in this case, we could say ancestral state estimates under two different models.

For my first reconstruction I’ll use a standard extended Mk ("ARD") model.

For my second reconstruction I’ll use the \(\Gamma\) variable rate model that I’ve been posting about for the past couple of days.

## get the character
parity_mode<-setNames(liolaemid.data$parity_mode,
  rownames(liolaemid.data))
## fit Mk model
parity_mk<-fitMk(liolaemid.tree,parity_mode,
  model="ARD",pi="fitzjohn",rand_start=TRUE)
parity_mk

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           O        V
## O -0.036703 0.036703
## V  0.000000 0.000000
## 
## Fitted (or set) value of pi:
## O V 
## 1 0 
## due to treating the root prior as (a) nuisance.
## 
## Log-likelihood: -64.270457 
## 
## Optimization method used was "nlminb"
## 
## R thinks it has found the ML solution.

## fit Gamma model
parity_gamma<-fitgammaMk(liolaemid.tree,parity_mode,
  model="ARD",pi="fitzjohn",min.alpha=0.01,
  rand_start=TRUE)
parity_gamma

## Object of class "fitgammaMk".
## 
## Fitted (or set) value of Q:
##           O        V
## O -0.056514 0.056514
## V  0.000000 0.000000
## 
## Fitted (or set) value of alpha rate heterogeneity
## (with 4 rate categories): 0.086767
## 
## Fitted (or set) value of pi:
## O V 
## 1 0 
## due to treating the root prior as (a) nuisance.
## 
## Log-likelihood: -63.435413 
## 
## Optimization method used was "nlminb"
## 
## R thinks it has found the ML solution.

If we compare the two models, we’ll find that there’s actually not all that much support for the \(\Gamma\) model.

anova(parity_mk,parity_gamma)

##                 log(L) d.f.      AIC    weight
## parity_mk    -64.27046    2 132.5409 0.5411457
## parity_gamma -63.43541    3 132.8708 0.4588543

Let’s run ancr to get our marginal states under each of the two models.

parity_mk.anc<-ancr(parity_mk)
parity_mk.anc

## Marginal ancestral state estimates:
##     O V
## 258 1 0
## 259 1 0
## 260 1 0
## 261 1 0
## 262 1 0
## 263 1 0
## ...
## 
## Log-likelihood = -64.270457

parity_gamma.anc<-ancr(parity_gamma)
parity_gamma.anc

## Marginal ancestral state estimates:
##     O V
## 258 1 0
## 259 1 0
## 260 1 0
## 261 1 0
## 262 1 0
## 263 1 0
## ...
## 
## Log-likelihood = -63.435413

Lastly, let’s graph our facing trees.

par(mfrow=c(2,1))
cols<-setNames(viridisLite::viridis(n=2),
  levels(parity_mode))
cex<-apply(parity_mk.anc$ace,1,
  function(x) if(any(x>0.95)) 0.2 else 0.5)
plot(parity_mk.anc,
  args.plotTree=list(type="arc",arc_height=0.5,
    mar=c(0.1,0.1,2.1,0.1)),
  args.nodelabels=list(piecol=cols,cex=cex),
  legend=FALSE)
mtext("a) ancestral states under standard Mk model",line=0,adj=0.05)
xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim
ylim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$y.lim
cex<-apply(parity_gamma.anc$ace,1,
  function(x) if(any(x>0.95)) 0.2 else 0.5)
plot(parity_gamma.anc,
  args.plotTree=list(type="arc",arc_height=0.5,
    xlim=xlim,ylim=-ylim[2:1],
    tips=1:Ntip(liolaemid.tree)+Ntip(liolaemid.tree),
    part=1,mar=c(0.1,0.1,2.1,0.1)),
  args.nodelabels=list(piecol=cols,cex=cex),
  legend=FALSE)
mtext(expression(paste("b) ancestral states with ",Gamma,
  "-distributed edge rates")),line=0,adj=0.05)
legend("bottomright",c("oviparity","viviparity"),
  pch=16,col=cols,cex=1,pt.cex=2,bty="n")

The plot is cool – but there’s very little disagreement between the two reconstructions, probably in large part because the best-fitting model has an irreversible (0 -> V) Q matrix in both cases.Perhaps I should’ve chosen a better example!