Today I tweeted about simulating & visualizing what I’ll call a “semi-threshold model.” The threshold model is one in which the observed trait is discretely valued but underlain by an unobserved continuous character called “liability.” Every time liability crosses some fixed threshold value, the discrete character changes state. To learn more about the threshold model in phylogenetic comparative biology, check out Felsenstein (2005, 2012), Revell (2014), or Chapter 7 of Revell & Harmon (2022).
A semi-threshold model, on the other hand, imagines that on some part of its range liability can be observed as a measurable quantitative trait. On another part of its range, by contrast, liability is not observable but nonetheless continues to evolve. A simple example of a trait that might evolve like this could be “horn length” which can increase or decrease, but when it shrinks to zero the horn is absent and can no longer be measured. The ability to produce a horn, however, might continue to evolve even if the horn is no longer present and measurable.
Here is my tweet.
Working on a "semi-threshold" model for #phytools. This is a trait evolution scenario where on one side of a threshold, our trait is present & measurable -- while on the other side it's absent, but its liability continues to evolve. Thought I'd make a go of visualizing. 👇 pic.twitter.com/MSYCURd4qO
— Liam Revell (@phytools_liam) December 10, 2024
Just for posteriority (and because I’ll undoubtedly want to re-create a similar figure in the future), here’s my code.
## load packages
library(phytools)
## set number of taxa
N<-8
## simulate a tree
tree<-pbtree(n=N,scale=100,tip.label=LETTERS[1:N])
## add a lot of non-branching nodes along each edge
tt<-map.to.singleton(make.era.map(tree,seq(0,100,by=1)))
## simulate liability (x)
x<-fastBM(tt,internal=TRUE,a=5)
## map the presence/absence of x as >0 or <0
maps<-list()
for(i in 1:nrow(tt$edge)){
if(x[tt$edge[i,1]]<0 && x[tt$edge[i,2]]<0)
maps[[i]]<-setNames(tt$edge.length[i],"0")
else if(x[tt$edge[i,1]]>0 && x[tt$edge[i,2]]>0)
maps[[i]]<-setNames(tt$edge.length[i],"1")
else {
d<-x[tt$edge[i,2]]-x[tt$edge[i,1]]
if(d>0)
maps[[i]]<-setNames(abs(c(-x[tt$edge[i,1]]/d,
x[tt$edge[i,2]]/d))*tt$edge.length[i],
0:1)
else
maps[[i]]<-setNames(abs(c(x[tt$edge[i,1]]/d,
x[tt$edge[i,2]]/d))*tt$edge.length[i],
1:0)
}
}
tt$maps<-maps
class(tt)<-c("simmap","phylo")
## create phenogram plot
par(mar=c(5.1,5.1,1.1,1.1),mfrow=c(2,1))
tmp<-phenogram(tt,x,spread.cost=c(1,0),
lty=setNames(c("solid","solid"),0:1),
colors=setNames(c("grey","blue"),0:1),
xlim=c(-5,120),ylim=rep(max(abs(x)),2)*c(-1,1),
ylab="liability / phenotype",las=1,cex.axis=0.9)
lines(c(0,100),c(0,0))
## label the threshold
text(0,-0.06*max(abs(x)),"threshold",cex=0.8,pos=4,font=3)
## label the trait ranges
text(x=-7,y=mean(c(0,max(abs(x)))),"trait present",srt=90,
cex=0.8,font=3)
lines(c(-3,-5,-5,-3),c(0.5,0.5,max(abs(x)),max(abs(x))))
text(x=-7,y=mean(c(0,-max(abs(x)))),"trait absent",srt=90,
cex=0.8,font=3)
lines(c(-3,-5,-5,-3),c(-0.5,-0.5,-max(abs(x)),-max(abs(x))))
## create contMap plot
obj<-contMap(tt,x[1:Ntip(tt)],anc.states=x[1:tt$Nnode+Ntip(tt)],
lims=rep(max(abs(x)),2)*c(-1,1),plot=FALSE,method="user")
## re-color so that any value of trait < 0 is grey
obj<-setMap(obj,c(rep("grey",500),
colorRampPalette(c("white","blue"))(500)))
## plot contMap
plot(obj,mar=c(1.1,5.1,0.1,1.1),legend=FALSE,
xlim=c(-5,120),ftype="off",ylim=c(-0.2*N,N))
segments(x0=rep(100,Ntip(tt))+1,x1=tmp[,"x"],
y0=1:Ntip(obj$tree),y1=1:Ntip(obj$tree),lty="dotted")
text(obj$tree$tip.label,x=tmp[,"x"],y=1:Ntip(obj$tree),
pos=4)
## add legend
add.color.bar(60,cols=colorRampPalette(c("white","blue"))(100),
prompt=FALSE,lims=c(0,max(abs(x))),x=0,y=-0.1*N,
title="quantitative trait value (if present)",subtitle="")
lines(x=c(0,10),y=rep(-0.2*N,2),lwd=6,lend=4)
lines(x=c(0,10),y=rep(-0.2*N,2),lwd=4,lend=4,col="grey")
text(x=11,y=-0.2*N,"trait absent",pos=4)
That’s it.
