Ancestral character estimation under the threshold model

Since posting about it on the anoleannals.org blog, I have been working on a first pass at ancestral character estimation under the threshold model.

Under the threshold model, the phenotypic trait of interest has a discretely valued presentation - for instance, color or habitat*; but a continuous, unmeasurable, underlying liability. (*Note, of course, that a characteristic like habitat might be judged by measuring a number of quantitative features about the environment - this is not a serious problem.) When the liability exceeds some threshold value, the discretely valued state of the observable character trait changes. So, for instance, the figure below (which I already used in a previous post) shows evolution of liability on a phylogenetic tree. When the liability crosses a threshold, the observable trait (here shown as a set of colors) changes state.

For more details on the threshold model, which derives originally from quantitative genetics, I would encourage you to check out Felsenstein (2012) and/or my prior posts on the subject (1, 2, 3).

The model provides several advantages. One is that it may be inherently more realistic than the discrete Markov model that we typically use for discretely-valued character traits. This model is the same as the "nucleotide model" that we use for DNA sequence evolution in phylogeny inference, and thus may not be well suited to complex traits typical of organismal phenotypes.

OK, so in principle, applying the threshold model to ancestral character estimation shouldn't be that different from applying it other problems. If we are to use Bayesian MCMC (as I have done in this case), we should just sample the liabilities for the tips of the tree and internal nodes; compute their probability under our model for evolution of the liabilities; and then multiply this probability by 1.0 if the liabilities predict our observed phenotypes, and 0.0 otherwise. Since it is not too interesting to reconstruct a trait that only has a binary presentation, I have also added MCMC sampling of the "positions" of the thresholds. (Our liabilities are unobserved, and thus scaleless, so these positions can only be interpreted in relative terms.)

Unfortunately, I have to run - so I will post more about this later, but below I've plotted the results from a single MCMC run of this method. The colored circles at the tips are the tip states for the discrete character, and the pie graphs at internal nodes represent the posterior probability of the internal nodes between in each of the four observed states, under the threshold model.

> source("ancThresh.R")
> tree<-pbtree(n=100)
> l<-fastBM(tree,internal=TRUE) # simulate liability
> # discretize!
> x<-rep("blue",length(l))
> x[l>0]<-"green"
> x[l>1]<-"red"
> x[l>2]<-"yellow"
> names(x)<-names(l)
> # run the MCMC
> result<-ancThresh(tree,x[1:length(tree$tip)],ngen=100000, control=list(sample=1000))
**** NOTE: no sequence provided, using alphabetical or numerical order
[1] "gen 1000"
[1] "gen 2000"
...

Pretty cool! I will post a lot more on this very soon.