Minor improvements to fitPagel for fitting a correlated binary trait evolution model in phytools

Earlier today, a colleague contacted me about interpreting a result from phytools::fitPagel (to fit the Pagel ’94 correlated binary trait evolution model).

I suggested she check for convergence to the ML solution for each model, but then was dismayed to realize that I hadn’t programmed fitPagel in such a way where it was at all straightforward to do so! Indeed, even though fitPagel uses fitMk (or, optionally, geiger::fitDiscrete) internally, it strips out only the parts of the "fitMk" object I deemed important at the time (mainly, the Q matrix and log-likelihood).

Since this function gets quite a lot of use, I thought it should minimally return the fitted "fitMk" class model objects and a message indicating whether or not R thought converge had been achieved – so now it does precisely that. To get this feature, users must simply update phytools from GitHub.

To demonstrate, I’m going to use a simple dataset – the same one that I used in my book with Luke Harmon – for spawning mode and paternal care of bony fish based on Benun Sutton & Wilson (2019).

We can start by loading phytools & checking the package version number.

## load package
library(phytools)

## check version number
packageVersion("phytools")

## [1] '2.4.6'

Next, let’s load our data as follows.

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

Now we’re ready to fit our models using fitPagel. The function fits both an independent model as well as the trait-dependent Pagel ’94 model.

## fit correlational model
bonyfish.pagel<-fitPagel(bonyfish.tree,paternal_care,
  spawning_mode,rand_start=TRUE)
bonyfish.pagel

## 
## Pagel's binary character correlation test:
## 
## Assumes "ARD" substitution model for both characters
## 
## Independent model rate matrix:
##            male|group male|pair none|group none|pair
## male|group  -0.003030  0.003030   0.000000  0.000000
## male|pair    0.001935 -0.001935   0.000000  0.000000
## none|group   0.002229  0.000000  -0.005259  0.003030
## none|pair    0.000000  0.002229   0.001935 -0.004164
## 
## Dependent (x & y) model rate matrix:
##            male|group male|pair none|group none|pair
## male|group   -0.02674  0.006940   0.019800  0.000000
## male|pair     0.00000  0.000000   0.000000  0.000000
## none|group    0.00000  0.000000  -0.004150  0.004150
## none|pair     0.00000  0.003123   0.002899 -0.006022
## 
## Model fit:
##             log-likelihood      AIC
## independent      -62.00739 132.0148
## dependent        -55.50975 127.0195
## 
## Hypothesis test result:
##   likelihood-ratio:  12.9953 
##   p-value:  0.0112989 
## 
## Model fitting method used was fitMk 
## 
## R thinks one or both likelihood optimizations did not converge.

According to the model print out, it appears that one or the other of our independent or dependent model optimizations may not have converged to the ML solution.

To find out which, we need to inspect the fitted models which are now (but were not before) returned as the model object element mk_fits.

bonyfish.pagel$mk_fits

## $independent
## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##            male|group male|pair none|group none|pair
## male|group  -0.003030  0.003030   0.000000  0.000000
## male|pair    0.001935 -0.001935   0.000000  0.000000
## none|group   0.002229  0.000000  -0.005259  0.003030
## none|pair    0.000000  0.002229   0.001935 -0.004164
## 
## Fitted (or set) value of pi:
## male|group  male|pair none|group  none|pair 
##       0.25       0.25       0.25       0.25 
## due to treating the root prior as (a) flat.
## 
## Log-likelihood: -62.007389 
## 
## Optimization method used was "nlminb"
## 
## R thinks it has found the ML solution.
## 
## 
## $dependent
## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##            male|group male|pair none|group none|pair
## male|group   -0.02674  0.006940   0.019800  0.000000
## male|pair     0.00000  0.000000   0.000000  0.000000
## none|group    0.00000  0.000000  -0.004150  0.004150
## none|pair     0.00000  0.003123   0.002899 -0.006022
## 
## Fitted (or set) value of pi:
## male|group  male|pair none|group  none|pair 
##       0.25       0.25       0.25       0.25 
## due to treating the root prior as (a) flat.
## 
## Log-likelihood: -55.509755 
## 
## Optimization method used was "nlminb"
## 
## R thinks optimization may not have converged.

This reveals that R thinks that the dependent model did not converge.

We can fit both models using different random starting values again as follows.

bonyfish.pagel<-fitPagel(bonyfish.tree,paternal_care,
  spawning_mode,rand_start=TRUE)
bonyfish.pagel

## 
## Pagel's binary character correlation test:
## 
## Assumes "ARD" substitution model for both characters
## 
## Independent model rate matrix:
##            male|group male|pair none|group none|pair
## male|group  -0.003030  0.003030   0.000000  0.000000
## male|pair    0.001935 -0.001935   0.000000  0.000000
## none|group   0.002229  0.000000  -0.005259  0.003030
## none|pair    0.000000  0.002229   0.001935 -0.004164
## 
## Dependent (x & y) model rate matrix:
##            male|group male|pair none|group none|pair
## male|group  -0.287481  0.105817   0.181664  0.000000
## male|pair    0.000000  0.000000   0.000000  0.000000
## none|group   0.000000  0.000000  -0.004022  0.004022
## none|pair    0.000000  0.003114   0.002871 -0.005985
## 
## Model fit:
##             log-likelihood      AIC
## independent      -62.00739 132.0148
## dependent        -55.35818 126.7164
## 
## Hypothesis test result:
##   likelihood-ratio:  13.2984 
##   p-value:  0.00990611 
## 
## Model fitting method used was fitMk 
## 
## R thinks both model optimizations converged.

This tells us that R thinks the model has converged.

Now, if we could not get fitPagel to optimize, perhaps we’d like to re-run our model in fitMk directly. By exporting the fitted model objects, we’ve made that much easier. That’s because our object now contains both the data (as state combinations, rather than the original traits) and a model design matrix, as well as our tree.

Let’s see.

refit_dependent<-fitMk(
  tree=bonyfish.pagel$tree,
  x=bonyfish.pagel$mk_fits[[2]]$data,
  model=bonyfish.pagel$mk_fits[[2]]$index.matrix)
refit_dependent

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##            male|group male|pair none|group none|pair
## male|group  -0.288362  0.106083   0.182279  0.000000
## male|pair    0.000000  0.000000   0.000000  0.000000
## none|group   0.000000  0.000000  -0.004022  0.004022
## none|pair    0.000000  0.003114   0.002871 -0.005985
## 
## Fitted (or set) value of pi:
## male|group  male|pair none|group  none|pair 
##       0.25       0.25       0.25       0.25 
## due to treating the root prior as (a) flat.
## 
## Log-likelihood: -55.358184 
## 
## Optimization method used was "nlminb"
## 
## R thinks optimization may not have converged.

That worked. We could’ve likewise run multiple optimization iterations, used different routines, or assumed a different root prior.

Lastly, let’s plot our fitted model.

## plot fitted models
plot(bonyfish.pagel,lwd.by.rate=TRUE)

A general clade box function for phytools

A few days ago, I posted some code illustrating the geometry of drawing a “box” around a clade on an arc- or fan-style plotted tree in R.

Here’s my tweet about the post which includes a nice illustration of its application:

On the #phytools blog, clade boxes for circular "arc" or "fan" style trees in R: https://t.co/2WqGrfaibB. pic.twitter.com/EJvp41LTxI

— Liam Revell (@phytools_liam) January 20, 2025

I asked a colleague if he thought this might be worth adding to phytools and his response was “yes, definitely useful.” The working function for that is below.

cladebox<-function(node,col="#0000FF40",...){
  pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
  phy<-list(edge=pp$edge,Nnode=pp$Nnode,
    tip.label=1:pp$Ntip)
  class(phy)<-"phylo"
  dd<-getDescendants(phy,node)
  dd<-dd[which(dd<=Ntip(phy))]
  if(hasArg(tip_buffer)) tip_buffer<-list(...)$tip_buffer
  else tip_buffer<-NULL
  if(pp$type%in%c("fan","arc")){
    if(node==(Ntip(phy)+1)) stem<-Inf
    else {
      parent<-phy$edge[which(phy$edge[,2]==node),1]
      stem<-sqrt(pp$xx[node]^2+pp$yy[node]^2)-
        sqrt(pp$xx[parent]^2+pp$yy[parent]^2)
    }
    if(hasArg(npoints)) npoints<-list(...)$npoints
    else npoints<-30
    if(hasArg(part)) part<-list(...)$part
    else {
      if(all(pp$xx>=0)&&all(pp$yy>=0)) part<-0.25
      else if(all(pp$yy>=0)) part<-0.5
      else part<-1
    }
    theta<-atan(pp$yy[dd]/pp$xx[dd])
    theta[pp$xx[dd]<0]<-pi+theta[pp$xx[dd]<0]
    theta[pp$xx[dd]>0&pp$yy[dd]<0]<-
      2*pi+theta[pp$xx[dd]>0&pp$yy[dd]<0]
    max_theta<-max(theta)+0.4*2*pi/Ntip(phy)*
      max(c(part,0.5))
    min_theta<-min(theta)-0.4*2*pi/Ntip(phy)*
      max(c(part,0.5))
    max_h<-max(sqrt(pp$xx^2+pp$yy^2))
    ff<-if(part>0.5) 1.6 else 
      if(0.25<part&&part<=0.5) 1.2 else 0.8
    if(is.null(tip_buffer)&&pp$show.tip.label){
      tip_buffer<-ff*(pp$x.lim[2]-max_h)
    } else tip_buffer<-0.02*max_h
    h1<-sqrt((pp$xx[node]^2+pp$yy[node]^2))-
      min(c(0.5*stem,0.02*max_h))
    h2<-max(sqrt(pp$xx[dd]^2+pp$yy[dd]^2))+tip_buffer
    tt<-seq(min_theta,max_theta,length.out=npoints)
    out_arc.x<-h2*cos(tt)
    out_arc.y<-h2*sin(tt)
    in_arc.x<-h1*cos(tt)
    in_arc.y<-h1*sin(tt)
    xx<-c(out_arc.x,in_arc.x[npoints:1])
    yy<-c(out_arc.y,in_arc.y[npoints:1])
    polygon(xx,yy,lwd=5,col=col,border=FALSE)
  } else if(pp$type%in%c("phylogram","cladogram")){
    if(pp$direction=="rightwards"){
      if(node==(Ntip(phy)+1)) stem<-0
      else {
        parent<-phy$edge[which(phy$edge[,2]==node),1]
        stem<-pp$xx[node]-pp$xx[parent]
      }
      max_h<-max(pp$xx)
      tip_buffer<-
        if(is.null(tip_buffer)&&pp$show.tip.label) 
          pp$x.lim[2]-max_h
        else 0.02*max_h
      h1<-pp$xx[node]-min(c(0.5*stem,0.02*max_h))
      h2<-max(pp$xx[dd])+tip_buffer
      xx<-c(h1,h2,h2,h1)
      yy<-c(rep(min(pp$yy[dd])-0.4,2),
        rep(max(pp$yy[dd])+0.4,2))
      polygon(xx,yy,border=FALSE,col=col)
    } else if(pp$direction=="leftwards"){
      if(node==(Ntip(phy)+1)) stem<-0
      else {
        parent<-phy$edge[which(phy$edge[,2]==node),1]
        stem<-pp$xx[parent]-pp$xx[node]
      }
      max_h<-max(pp$xx)
      tip_buffer<-
        if(is.null(tip_buffer)&&pp$show.tip.label) 
          pp$x.lim[1]
      else -0.02*max_h
      h1<-pp$xx[node]+min(c(0.02*max_h,0.5*stem))
      h2<-max(pp$xx[dd])+tip_buffer
      xx<-c(h1,h2,h2,h1)
      yy<-c(rep(min(pp$yy[dd])-0.4,2),
        rep(max(pp$yy[dd])+0.4,2))
      polygon(xx,yy,border=FALSE,col=col)
    } else if(pp$direction=="upwards"){
      if(node==(Ntip(phy)+1)) stem<-0
      else {
        parent<-phy$edge[which(phy$edge[,2]==node),1]
        stem<-pp$yy[node]-pp$yy[parent]
      }
      max_h<-max(pp$yy)
      tip_buffer<-
        if(is.null(tip_buffer)&&pp$show.tip.label) 
          pp$y.lim[2]-max_h
      else 0.02*max_h
      xx<-c(rep(min(pp$xx[dd])-0.4,2),
        rep(max(pp$xx[dd])+0.4,2))
      h1<-pp$yy[node]-min(c(0.02*max_h,0.5*stem))
      h2<-max(pp$yy[dd])+tip_buffer
      yy<-c(h1,h2,h2,h1)
      polygon(xx,yy,border=FALSE,col=col)
    } else if(pp$direction=="downwards"){
      if(node==(Ntip(phy)+1)) stem<-0
      else {
        parent<-phy$edge[which(phy$edge[,2]==node),1]
        stem<-pp$yy[parent]-pp$yy[node]
      }
      max_h<-max(pp$yy)
      tip_buffer<-
        if(is.null(tip_buffer)&&pp$show.tip.label) 
          pp$y.lim[1]
      else -0.02*max_h
      xx<-c(rep(min(pp$xx[dd])-0.4,2),
        rep(max(pp$xx[dd])+0.4,2))
      h1<-pp$yy[node]+min(c(0.02*max_h,0.5*stem))
      h2<-max(pp$yy[dd])+tip_buffer
      yy<-c(h1,h2,h2,h1)
      polygon(xx,yy,border=FALSE,col=col)
    }
  }
  invisible(list(x=xx,y=yy))
}

I’m not going to peel through all the layers of the onion, but some of the major ideas include that the method will plot an appropriate clade box depending on our plot style, and leave a buffer around the box whose width depends on the amount of space our plotting label left in our device for the tip labels. It even adjusts the “rootward” buffer depending on how closely spaced the MRCA of our focal is to its parent.

Let’s test it out.

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

## 
## Phylogenetic tree with 90 tips and 89 internal nodes.
## 
## Tip labels:
##   Xenomystus_nigri, Chirocentrus_dorab, Talismania_bifurcata, Alepocephalus_tenebrosus, Misgurnus_bipartitus, Opsariichthys_bidens, ...
## 
## Rooted; includes branch length(s).

plotTree(bonyfish.tree,direction="upwards",lwd=1,
  fsize=0.6,ftype="i")
nodelabels(cex=0.5,bg="white")

Unfortunately, I don’t have any interesting clades to label in this tree, so what I’m going to do is draw boxes around the clades descended from the nodes numbered 96, 112, 126, 139, and 168 of this plot as follows.

plotTree(bonyfish.tree,direction="upwards",lwd=1,
  fsize=0.6,ftype="i")
cladebox(node=96)
cladebox(node=112)
cladebox(node=139)
cladebox(node=168)

We can also draw our clade boxes first (as I illustrated last time), and then overlay our tree. Here’s what that might look like for the same tree & set of clade boxes.

plotTree(bonyfish.tree,direction="upwards",lwd=1,
  fsize=0.6,ftype="i",plot=FALSE)
cols<-RColorBrewer::brewer.pal(n=4,"Pastel1")
cladebox(node=96,col=cols[1])
cladebox(node=112,col=cols[2])
cladebox(node=139,col=cols[3])
cladebox(node=168,col=cols[4])
plotTree(bonyfish.tree,direction="upwards",lwd=1,
  fsize=0.6,ftype="i",add=TRUE)

Of course, the point of my prior post on this topic was how to solve the much more complicated task of overlaying clade boxes on top of a circular fan or arc-style tree.

Naturally, our function does that too.

par(fg="transparent")
plotTree(bonyfish.tree,direction="upwards",lwd=1,
  fsize=0.7,ftype="i",type="arc",arc_height=0.25,
  color="transparent")
par(fg="black")
cols<-RColorBrewer::brewer.pal(n=4,"Pastel1")
cladebox(node=96,col=cols[1])
cladebox(node=112,col=cols[2])
cladebox(node=139,col=cols[3])
cladebox(node=168,col=cols[4])
plotTree(bonyfish.tree,direction="upwards",lwd=1,
  fsize=0.7,ftype="i",type="arc",arc_height=0.25,
  add=TRUE)

The function also invisibly returns the coordinates of the polygons, so we can even overlay them in a custom way if we want.

par(fg="transparent")
plotTree(bonyfish.tree,direction="upwards",
  fsize=0.7,ftype="i",type="arc",arc_height=0.25,
  lwd=3,color="grey")
par(fg="black")
cols<-RColorBrewer::brewer.pal(n=4,"Pastel1")
cb_96<-cladebox(node=96,col="transparent")
polygon(cb_96,col=make.transparent(cols[1],0.4),
  density=20,border=cols[1],lwd=2)
cb_112<-cladebox(node=112,col="transparent")
polygon(cb_112,col=make.transparent(cols[2],0.4),
  density=20,border=cols[2],lwd=2)
cb_139<-cladebox(node=139,col="transparent")
polygon(cb_139,col=make.transparent(cols[3],0.4),
  density=20,border=cols[3],lwd=2)
cb_168<-cladebox(node=168,col="transparent")
polygon(cb_168,col=make.transparent(cols[4],0.4),
  density=20,border=cols[4],lwd=2)
plotTree(bonyfish.tree,direction="upwards",lwd=1,
  fsize=0.7,ftype="i",type="arc",arc_height=0.25,
  add=TRUE)

Overall, I think this is kind of cool.