Null distribution by simulation for brownie.lite()

I just posted a new version of brownie.lite() on my R-phylogenetics page (direct link to code here). Just to remind you, this function implements the non-censored evolutionary rates test of O'Meara et al. (2006). In the method we fit two different models using the technique of maximum likelihood - one in which one evolutionary rate prevails, and another in which the evolutionary rate shifts one or multiple times between two or more evolutionary rates on different branches of the tree. We can then compare the fit of these models by, for instance, comparing 2 × their log-likelihood ratio to a Χ²-distribution.

We could also compare the likelihood ratio to likelihood ratios obtained by simulating under the null hypothesis. (In fact, I believe that the authors suggested this in their original paper.) I have now implemented this test in brownie.lite(). I am also presently conducting a mini-experiment to compare the type I error rates of each method. Since the likelihood ratio is only asymptotically Χ²-distributed for large n, we might expect that the simulation test would be more conservative than the Χ² test. So far I have not found a lot of evidence to suggest that it is for n=10 and n=20.

I'll also confess to a small fudge in the implementation of this method. Occasionally (i.e., about 1% of the time or so), brownie.lite() fails to converge during optimization. Rather than fix this, to avoid this problem when trying to compute the P-value by simulation, I just discarded any simulated datasets that were troublesome to optimize. I don't see how this would result in any bias, but I don't know. In the near future, I will fix this problem in brownie.lite() (and hopefully clean up the function in general, as it is pretty messy).

If you'd like to participate in my mini experiment, I encourage you to download the new version of brownie.lite() and try and run the following code (hopefully increasing nrep and varying nsp. Please post the results.

require(phytools)
require(geiger)
source("brownie.lite.R") # load new version of brownie.lite()
typeI.x2<-0; typeI.sim<-0
nrep<-100; nsp=20
for(i in 1:nrep){
  # simulate tree
  tree<-drop.tip(birthdeath.tree(b=1,d=0,taxa.stop=(nsp+1)), as.character(nsp+1))
  # rescale tree
  tree$edge.length<-tree$edge.length/mean(diag(vcv(tree)))
  # simulate rate regime
  x<-sim.char(tree,model.matrix=list(matrix(c(-1,1,1,-1),2,2)), model="discrete")[,,1]
  while(length(unique(x))==1)
    x<-sim.char(tree,model.matrix=list(matrix(c(-1,1,1,-1), 2,2)),model="discrete")[,,1]
  mtree<-make.simmap(tree,x)
  # simulate data under the null
  y<-fastBM(tree)
  # fit the model & compute X^2 P-value
  result.x2<-brownie.lite(mtree,y,test="chisq")
  # fit the model & compute simulation P-value
  result.sim<-brownie.lite(mtree,y,test="simulation")
  # count the frequence of P<0.05
  typeI.x2<-typeI.x2+(result.x2$P.chisq<0.05)/nrep
  typeI.sim<-typeI.sim+(result.sim$P.sim<0.05)/nrep
  # something to look at while this is running
  # (should converge to 0.05)
  print(paste("X^2 type I error:",typeI.x2*nrep/i, "; Simulation type I error:",typeI.sim*nrep/i))
}