New phytools version on CRAN; bug fix in make.simmap(...,Q="mcmc")

A new version of phytools is now on CRAN, although it will probably take a few days for binaries to be built and then percolate through all the mirror repositories.

Unfortunately, as soon it this version was accepted, I discovered a small bug with the function make.simmap. make.simmap is a popular phytools function that implements the method of stochastic character mapping. The bug is present in only the full hierarchical Bayesian method (Q="mcmc"), and was introduced in the latest version of phytools because I know use the function fitMk to compute the likelihood of the Mk model internally.

Here is how the bug manifests:

library(phytools)

## Loading required package: ape
## Loading required package: maps

packageVersion("phytools")

## [1] '0.5.0'

## load data
data(anoletree)
## pull out the data for tips so we can test make.simmap
x<-getStates(anoletree,"tips")
## default method (empirical Bayesian method)
eb.trees<-make.simmap(anoletree,x,nsim=100,model="ER")

## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##             CG          GB          TC          TG          Tr          Tw
## CG -0.13884868  0.02314145  0.02314145  0.02314145  0.02314145  0.02314145
## GB  0.02314145 -0.13884868  0.02314145  0.02314145  0.02314145  0.02314145
## TC  0.02314145  0.02314145 -0.13884868  0.02314145  0.02314145  0.02314145
## TG  0.02314145  0.02314145  0.02314145 -0.13884868  0.02314145  0.02314145
## Tr  0.02314145  0.02314145  0.02314145  0.02314145 -0.13884868  0.02314145
## Tw  0.02314145  0.02314145  0.02314145  0.02314145  0.02314145 -0.13884868

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

##        CG        GB        TC        TG        Tr        Tw 
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

## Done.

## try full Bayesian method (should fail)
fb.trees<-make.simmap(anoletree,x,nsim=100,model="ER",Q="mcmc")

## Running MCMC burn-in. Please wait....

## Error in if (p.odds >= runif(n = 1)) {: argument is of length zero

I have fixed this bug and the fixed version can be already be installed from GitHub using devtools. Here I will demo it by loading the function from source:

source("../phytools/R/make.simmap.R")
fb.trees<-make.simmap(anoletree,x,nsim=100,model="ER",Q="mcmc")

## Running MCMC burn-in. Please wait....
## Running 10000 generations of MCMC, sampling every 100 generations.
## Please wait....
## 
## make.simmap is simulating with a sample of Q from
## the posterior distribution
## 
## Mean Q from the posterior is
## Q =
##             CG          GB          TC          TG          Tr          Tw
## CG -0.11975495  0.02395099  0.02395099  0.02395099  0.02395099  0.02395099
## GB  0.02395099 -0.11975495  0.02395099  0.02395099  0.02395099  0.02395099
## TC  0.02395099  0.02395099 -0.11975495  0.02395099  0.02395099  0.02395099
## TG  0.02395099  0.02395099  0.02395099 -0.11975495  0.02395099  0.02395099
## Tr  0.02395099  0.02395099  0.02395099  0.02395099 -0.11975495  0.02395099
## Tw  0.02395099  0.02395099  0.02395099  0.02395099  0.02395099 -0.11975495
## and (mean) root node prior probabilities
## pi =
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

## Done.

obj<-summary(fb.trees)
plot(obj,fsize=0.6,ftype="i",ylim=c(-2,Ntip(anoletree)))
colors<-setNames(palette()[1:length(unique(x))],sort(unique(x)))
add.simmap.legend(colors=colors,x=0,y=-2,prompt=FALSE,vertical=FALSE)

This fix can be obtained most easily by installing phytools from GitHub using devtools. In a new R session, run:

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

That's it.