I'll be presenting on ancestral character estimation under the threshold model at the upcoming symposium entitled "Mathematics for an evolving biodiversity" in Montreal, Quebec in a couple of weeks. Last night, I thought it might be neat to take my existing function, bmPlot (which does discrete-time simulation & visualization under the threshold model), and modify it so that we could also visualize evolution of a threshold character up the branches of the tree.
The way I did this was by tracking the evolution of liability, translating liability into the threshold trait, and then aggregating adjacent generations with the same state into a stochastic-map style tree. Code for the new version of bmPlot is here.
One of the predictions of the threshold model is that the character should change back & forth across the threshold many times very rapidly when near the threshold. (Actually, infinitely many times in theory, under continuous time; remember, ours will be a discrete-time simulation.) That means that (compared to a 'nucleotide-model' type simulation) there some periods when the character changes very rapidly many times; and other long periods when the character doesn't change at all. Let's see how that bears out in simulation.
simulating with both taxa-stop (n) and time-stop (t) is
performed via rejection sampling & may be slow
51 trees rejected before finding a tree
> X<-bmPlot(tree,type="threshold",thresholds=c(0,0.5,2), anc=1,sig2=1/1000,ngen=max(nodeHeights(tree)), return.tree=TRUE)
So our prediction appears to be exactly borne out - large parts of the tree, including very long branches, have few to no changes; however in areas when a lineage is close to one of the thresholds, it tends to change back & forth rapidly. Cool.