Bug fix for color bar in "contMap" & "densityMap" plots

It was recently reported to me that the color gradient legend labels in the plotting method for "contMap" are messed up. This was evidently due to a bug I introduced the option to plot the color legend vertically.

I just pushed a bug fix (which is basically just a deletion of the introduced bug.

First, here is the bug:

library(phytools)
obj<-contMap(tree,x,plot=FALSE)
plot(obj,legend=2.1)

Note the rightmost legend label is 'floating' too high about the legend line.

Now the fix. Since the function add.color.bar is called internally, R won't use it unless we load the source code of the other functions as well:

source("https://raw.githubusercontent.com/liamrevell/phytools/master/R/plotBranchbyTrait.R")
source("https://raw.githubusercontent.com/liamrevell/phytools/master/R/densityMap.R")
source("https://raw.githubusercontent.com/liamrevell/phytools/master/R/contMap.R")
plot(obj,legend=2.1)

This can also be installed from GitHub of course:

library(devtools)
install_github("liamrevell/phytools")

On the accuracy of reconstructed tip states using make.simmap

Earlier (1, 2) I blogged about reconstructing tip states when they are unknown.

What one needs to do to show that this method is working as designed, for discrete traits, at least, is not to merely measure whether the most probable state is correct - but to compare the posterior probabilities to the frequency of times said state is in fact correct. In other words, we expect a tip state that is estimated to have a posterior probability of being in state a of 0.6, for instance, to in fact be in state a 60% of the time. Another way of saying this is that if we take all tip states with PP(a) = x, x×100% of the time these tips should have state a. This was pointed out to me as a good way to measure the performance of a Bayesian estimation method by Brian O'Meara a few yeara ago, and, of course, he was correct.

Since posterior probabilities can take any value on [0,1], in practice this means we need to bin our posterior probabilities.

Here, that's exactly what I have done. First, I simulated trees & discrete character data 200 times; next I arbitrarily set 40% (40/100 taxa) of tip values to be unknown; then I sampled their true values using stochastic mapping; and, finally, I binned the posterior probabilities and asked what with what frequency, say, tips with posterior probabilities between PP(a)=0.05 & PP(a)=0.1 were in fact in state a, and so on. Here I elected to fix the transition matrix, Q, to its true value for stochastic mapping - essentially to avoid confounding error in the estimation of tip states given Q with estimation of Q itself.

Here is what I did & found piece by piece.

Load packages & simulate trees. Set Q:

library(phytools)
## simulate trees & set Q
trees<-pbtree(n=200,scale=1,nsim=200)
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-letters[1:2]
Q

##    a  b
## a -1  1
## b  1 -1

Create vectors to store our results. I have to count all instances in which the posterior probabilities were between x & y (denoms), as well as the frequencies for each of those intervals in which each state was correct.

nums<-denoms<-
    setNames(rep(0,20),paste(seq(0,0.95,0.05),"<->",
    seq(0.05,1,0.05),sep=""))

Now let's run our simulations:

for(i in 1:length(trees)){
    missing<-sample(trees[[i]]$tip.label,40)
    x<-sim.history(trees[[i]],Q,message=FALSE)$states
    X<-to.matrix(x,letters[1:2])
    X[missing,]<-rep(0.5,2)
    m<-make.simmap(trees[[i]],X,nsim=100,Q=Q,message=FALSE)
    tips<-summary(m)$tips[,letters[1:2]]
    n<-round(20*tips[missing,1]*to.matrix(x[missing],
        letters[1:2])[,1])
    d<-round(20*tips[missing,1])
    for(j in 1:length(n)){
        nums[n[j]]<-nums[n[j]]+1
        denoms[d[j]]<-denoms[d[j]]+1
    }
    if(i%%20==0&&i!=0) cat(".\n") else cat(".")
}

## ....................
## ....................
## ....................
## ....................
## ....................
## ....................
## ....................
## ....................
## ....................
## ....................

cat("Done simulations.\n")

## Done simulations.

Finally, we have our results!

print(nums/denoms)

##   0<->0.05 0.05<->0.1 0.1<->0.15 0.15<->0.2 0.2<->0.25 0.25<->0.3 
## 0.05668449 0.07736390 0.15746421 0.19537275 0.22800000 0.27450980 
## 0.3<->0.35 0.35<->0.4 0.4<->0.45 0.45<->0.5 0.5<->0.55 0.55<->0.6 
## 0.31952663 0.37956204 0.43801653 0.48739496 0.42500000 0.58955224 
## 0.6<->0.65 0.65<->0.7 0.7<->0.75 0.75<->0.8 0.8<->0.85 0.85<->0.9 
## 0.69677419 0.63068182 0.71984436 0.78932584 0.86315789 0.89079563 
## 0.9<->0.95   0.95<->1 
## 0.94976077 0.98012232

plot(seq(0.025,0.975,0.05),nums/denoms,xlim=c(0,1),ylim=c(0,1),
    xlab="Posterior probability (state 'a')",
    ylab="Observed frequency in state 'a'",pch=21,bg="grey",
    cex=1.2)
lines(c(0,1),c(0,1),lty="dashed",lwd=2,
    col=make.transparent("blue",0.4))

If our results are on a 1:1 line, as I hope will be the case, then most likely the method is working as designed.

Stochastic mapping discrete character histories with missing data at some tips

Stochastic mapping in phytools permits some tips to have unknown states. This is done by using a prior probability distribution on the tips that is flat across all possible states. This leaves the posterior probabilities at internal nodes largely unaffected compared to just dropping said tips from the tree. (Although I do not know for sure that they could not theoretically be affected, it is hard to see how.)

Here's an example using simulated data.

Load libraries:

library(phytools)
library(phangorn)

Simulate data including missing tip values for some taxa:

## simulate tree & data
N<-60 ## number of tips
n.miss<-20 ## number with missing values
tree<-pbtree(n=N,scale=1)
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-c("a","b")
x<-sim.history(tree,Q)$states

## Done simulation(s).

X<-to.matrix(x,c("a","b"))
## missing
missing<-sample(tree$tip.label,n.miss)
tips<-sapply(missing,function(x,y) which(y==x),
    y=tree$tip.label)
X[missing,]<-rep(0.5,2)

Now perform stochastic mapping:

## stochastic mapping
## with missing
m1<-make.simmap(tree,X,nsim=500)

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##            a          b
## a -0.8958661  0.8958661
## b  0.8958661 -0.8958661
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##   a   b 
## 0.5 0.5

## Done.

plot(sum1<-summary(m1),colors=setNames(c("blue","red"),c("a","b")),
    ftype="off")

## without missing
m2<-make.simmap(drop.tip(tree,missing),
    X[setdiff(tree$tip.label,missing),],
    nsim=500)

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##           a         b
## a -0.895866  0.895866
## b  0.895866 -0.895866
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##   a   b 
## 0.5 0.5

## Done.

plot(sum2<-summary(m2),colors=setNames(c("blue","red"),c("a","b")),
    ftype="off")

We can also visualize the concordance of internal node reconstructions using cophylo:

obj<-cophylo(tree,drop.tip(tree,missing),rotate=FALSE)
plot(obj,link.type="curved",fsize=c(0.6,0.7))
tiplabels.cophylo(pie=sum1$tips,
    piecol=setNames(c("blue","red"),c("a","b")),
    which="left",cex=0.3)
nodelabels.cophylo(pie=sum1$ace,
    piecol=setNames(c("blue","red"),c("a","b")),
    which="left",cex=0.5)
tiplabels.cophylo(pie=sum2$tips,
    piecol=setNames(c("blue","red"),c("a","b")),
    which="right",cex=0.3)
nodelabels.cophylo(pie=sum2$ace,
    piecol=setNames(c("blue","red"),c("a","b")),
    which="right",cex=0.5)

Finally, we can compare the nodes by mapping one to the other across the two analyses, and then comparing the posterior probabilities directly:

## compare nodes
M<-matchNodes(tree,drop.tip(tree,missing),"distances")
plot(sum1$ace[as.character(M[!is.na(M[,2]),1]),],
    sum2$ace,xlab="with missing",ylab="prune missing")
lines(c(0,1),c(0,1),lty="dashed",col="grey",lwd=2)
lines(c(0,1),c(0,1),lty="dashed",col="grey",lwd=2)

Note that for the case of computing the likelihood of a tree given a data pattern at the tip, we would set the probabilities of each state equal to 1.0. Here I'm computing the likelihood of a particular character pattern at the nodes & missing tips (& a model), given the data for the tips so I have set the prior probabilities for the missing nodes to be 1/(number of states). I'm not sure this makes a difference (other than numerically affecting the computed likelihood); however I suppose it could.