Computing the posterior probabilities for tip nodes from ancThresh

The phytools function ancThresh, which conducts ancestral state reconstruction of a discrete character under the threshold model using Bayesian MCMC (Revell In press), also allows uncertain tip nodes in the tree. Uncertainty in the tip states are treated as a set of prior probability distributions on the states of terminal species in the tree. In this case, we can also compute posterior probabilities that each tip is in each state. These are automatically plotted by ancThresh if ancThresh(...,control=list(tipcol="estimated"). The burn-in that is used to compute these posterior probabilities is the user supplied burn-in in ancThresh(...,control=list(burnin)).

Unfortunately, and unlike the posterior probabilities for internal nodes, these probabilities are not returned to the user nor are they especially easy to compute. This is because the list component $mcmc only contains the implied states at internal nodes; and thus to get the posterior probabilities for tips, we need to use $liab and $par (which contains the sampled positions of the thresholds for each generation of the MCMC) to get the states.

Here is some code that can be used to compute these probabilities. In the example, mcmc is our object returned by ancThresh:

burnin<-200000 ## for instance
## find the row with the first post-burnin sample
ii<-which(mcmc$par[,"gen"]==burnin)
## number of tips
n<-length(tree$tip.label)
## states for the discrete character
states<-colnames(mcmc$par)[2:(ncol(mcmc$par)-1)]
## create and populate our matrix of posterior probabilities
PP<-matrix(NA,n,length(states),dimnames=
  list(colnames(mcmc$liab[1:n]),states))
for(i in 1:n){
  x<-vector(length=nrow(mcmc$liab)-ii+1)
  for(j in ii:nrow(mcmc$liab))
    x[j-ii+1]<-threshState(mcmc$liab[j,i],
      mcmc$par[j,2:(ncol(mcmc$par)-1)])
  PP[i,]<-summary(factor(x,levels=states))/length(x)
}

That's it.