Updated default orientation of plotted Q matrix, for discrete character evolution (Mk models) -- and a new option to add them to an existing plot!

I recently learned on Twitter of a new paper by Kate Thomas et al. (2022) modeling the evolution of pupil shape and orientation across amphibians.

Our paper led by @katiethomas10 & @RaynaCBell investigating evolution of the diverse pupil shapes of amphibians is out in @BiolJLinnSoc. We found pupil shape to be evolutionarily labile & that circular pupils are linked to fossorial and aquatic lifestyles🐸https://t.co/cbOXbIfUgA pic.twitter.com/bi81CyrA6R

— Ryan Schott (@rkschott) September 7, 2022

One thing that jumped out at me straight away, and reminded me of something that I've been meaning to do for a long time, was the orientation of their graphed Q matrix (and its superposition on top of the plotted tree).

Even though I have recently added a bunch of new cool functionality to the plot.Qmatrix method in phytools (e.g., 1, 2, 3), the default orientation of the plotted Q matrix (that is, how it is rotated around the origin) is fixed.

Since the default orientation is to start at 0 and then evenly space the trait space in an even-angle fashion, especially for low dimensional (e.g., 3 × 3) Q matrices, this can create a weird visual impression that the plotted matrix is somehow 'off-kilter'.

For example, take a look at this plotted 3 × 3 matrix from an earlier post on this blog.

Here's what the update looks like – using the argument rotate=0 to show the “old” plot orientation as well as the new default.

library(phytools)

## Loading required package: ape

## Loading required package: maps

packageVersion("phytools")

## [1] '1.2.3'

Q2<-as.Qmatrix(matrix(c(
    -0.9,0.9,
    1.2,-1.2),2,2,byrow=TRUE,
    dimnames=list(0:1,0:1)))
Q3<-as.Qmatrix(matrix(c(
    -0.9,0.9,0.0,
    1.2,-2.1,0.9,
    0.0,1.2,-1.2),3,3,byrow=TRUE,
    dimnames=list(0:2,0:2)))
Q4<-as.Qmatrix(matrix(c(
    -0.9,0.9,0.0,0.0,
    1.2,-2.1,0.9,0.0,
    0.0,1.2,-2.1,0.9,
    0.0,0.0,1.2,-1.2),4,4,byrow=TRUE,
    dimnames=list(0:3,0:3)))
Q5<-as.Qmatrix(matrix(c(
    -0.9,0.9,0.0,0.0,0.0,
    1.2,-2.1,0.9,0.0,0.0,
    0.0,1.2,-2.1,0.9,0.0,
    0.0,0.0,1.2,-2.1,0.9,
    0.0,0.0,0.0,1.2,-1.2),5,5,byrow=TRUE,
    dimnames=list(0:4,0:4)))
Q6<-as.Qmatrix(matrix(c(
    -0.9,0.9,0.0,0.0,0.0,0.0,
    1.2,-2.1,0.9,0.0,0.0,0.0,
    0.0,1.2,-2.1,0.9,0.0,0.0,
    0.0,0.0,1.2,-2.1,0.9,0.0,
    0.0,0.0,0.0,1.2,-2.1,0.9,
    0.0,0.0,0.0,0.0,1.2,-1.2),6,6,byrow=TRUE,
    dimnames=list(0:5,0:5)))
Q7<-as.Qmatrix(matrix(c(
    -0.9,0.9,0.0,0.0,0.0,0.0,0.0,
    1.2,-2.1,0.9,0.0,0.0,0.0,0.0,
    0.0,1.2,-2.1,0.9,0.0,0.0,0.0,
    0.0,0.0,1.2,-2.1,0.9,0.0,0.0,
    0.0,0.0,0.0,1.2,-2.1,0.9,0.0,
    0.0,0.0,0.0,0.0,1.2,-2.1,0.9,
    0.0,0.0,0.0,0.0,0.0,1.2,-1.2),7,7,byrow=TRUE,
    dimnames=list(0:6,0:6)))
par(mfrow=c(3,4))
plot(Q2,rotate=0)
mtext("a) 2-state \"old\"",line=0,adj=0)
plot(Q2)
mtext("b) 2-state \"new\"",line=0,adj=0)
plot(Q3,rotate=0)
mtext("c) 3-state \"old\"",line=0,adj=0)
plot(Q3)
mtext("d) 3-state \"new\"",line=0,adj=0)
plot(Q4,rotate=0)
mtext("e) 4-state \"old\"",line=0,adj=0)
plot(Q4)
mtext("f) 4-state \"new\"",line=0,adj=0)
plot(Q5,rotate=0)
mtext("g) 5-state \"old\"",line=0,adj=0)
plot(Q5)
mtext("h) 5-state \"new\"",line=0,adj=0)
plot(Q6,rotate=0)
mtext("i) 6-state \"old\"",line=0,adj=0)
plot(Q6)
mtext("j) 6-state \"new\"",line=0,adj=0)
plot(Q7,rotate=0)
mtext("k) 7-state \"old\"",line=0,adj=0)
plot(Q7)
mtext("l) 7-state \"new\"",line=0,adj=0)

Cool.

In addition, the tweet also shows a visual in which the Q matrix is graphed on top of the plotted tree.

I decided to add this functionality to plot.Qmatrix too, using the logical argument add.

In this example, I'm going to use a phylogeny & squamate digit number example subsampled from a dataset by Brandley et al. (2008).

In this case, rather than paste my fitted Q matrix right on top of the plotted tree, I'll graph my tree through only 75% of the space, and then add the fitted Mk model to the remaining quadrant.

The code below also includes a hacky method to visualize a discrete trait at the tips using thicken lines.

The arguments I used to adjust the position of the plotted Q matrix where xlim and ylim.

For instance, by setting the first value of ylim to be a relatively small (negative) number, and the latter value to be large, we get the Q matrix into the lower half (and so on). The matrix will always plot between -1 and 1 on both axes.

## read data & tree
sqData<-read.csv("http://www.phytools.org/Rbook/6/squamate-data.csv",row.names=1)
sqTree<-read.nexus("http://www.phytools.org/Rbook/6/squamate.tre")
chk<-geiger::name.check(sqTree,sqData)
sqTree.pruned<-drop.tip(sqTree,chk$tree_not_data)
sqData.pruned<-sqData[!(rownames(sqData)%in%
    chk$data_not_tree),,drop=FALSE]
toes<-setNames(as.factor(sqData.pruned[,"rear.toes"]),
    rownames(sqData.pruned))
## fit directional (loss only) trait evolution model
directional.model<-matrix(c(
    0,0,0,0,0,0,
    1,0,0,0,0,0,
    0,2,0,0,0,0,
    0,0,3,0,0,0,
    0,0,0,4,0,0,
    0,0,0,0,5,0),6,6,byrow=TRUE,
dimnames=list(0:5,0:5))
fitDirectional<-fitMk(sqTree.pruned,toes,model=directional.model,
    pi="fitzjohn")
fitDirectional

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##          0         1         2         3         4         5
## 0 0.000000  0.000000  0.000000  0.000000  0.000000  0.000000
## 1 0.024241 -0.024241  0.000000  0.000000  0.000000  0.000000
## 2 0.000000  0.103141 -0.103141  0.000000  0.000000  0.000000
## 3 0.000000  0.000000  0.109233 -0.109233  0.000000  0.000000
## 4 0.000000  0.000000  0.000000  0.083290 -0.083290  0.000000
## 5 0.000000  0.000000  0.000000  0.000000  0.004369 -0.004369
## 
## Fitted (or set) value of pi:
## 0 1 2 3 4 5 
## 0 0 0 0 0 1 
## due to treating the root prior as (a) nuisance.
## 
## Log-likelihood: -118.387049 
## 
## Optimization method used was "nlminb"

## create plot
plotTree(sqTree.pruned,type="fan",ftype="off",
    color="grey",lwd=1,part=0.75)
obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
n<-Ntip(sqTree.pruned)
cols<-setNames(grey.colors(n=length(levels(toes)),0.9,0.1),
    levels(toes))
par(lend=3)
for(i in 1:Ntip(sqTree.pruned)){
    cc<-if(obj$xx[i]>0) 10 else -10
    th<-atan(obj$yy[i]/obj$xx[i])
    segments(obj$xx[i],obj$yy[i],
        obj$xx[i]+cc*cos(th),
        obj$yy[i]+cc*sin(th),
        lwd=20,
        col=cols[toes[sqTree.pruned$tip.label[i]]])
}
legend("topleft",levels(toes),pch=15,col=cols,pt.cex=1.5,
    cex=0.8,bty="n",title="number of digits")
plot(fitDirectional,width=TRUE,add=TRUE,show.zeros=FALSE,
    text=FALSE,color=TRUE,
    xlim=c(-4,0.5),ylim=c(-1,3.5),signif=1)
text(0,1.1,"fitted Mk model\nfor digit number evolution")

OK. It works!