# Bootstrap pvalues: strange adjustment -- pbkrtest PBmodcomp LRT

As far as I can tell, PBmodcomp is doing some straight up bootstrapping, similar to this Faraway example for lmer models mmod and rmod (maximal and reduced models):

lrstat <- numeric(1000)
for(i in 1:1000){
rmath <- unlist(simulate(rmod))
bmod <- refit(mmod, rmath)
smod <- refit(rmod, rmath)
lrstat[i] <- 2*(logLik(bmod)-logLik(smod))
}
pvalue <- mean(lrstat > olrt)


That Faraway pvalue is a bit off, see http://www.ncbi.nlm.nih.gov/pmc/articles/PMC379178/

The PBmodcomp documentation is in line with the above link as far as calculating the pvalue http://www.jstatsoft.org/v59/i09/paper.

However, in the code for PBmodcomp something else is going on as far as calculating the pvalue:

refpos <- ref[ref>0]
nsim <- length(ref)
npos <- length(refpos)
n.extreme <- sum(tobs < refpos)
p.PB  <- (1+n.extreme) / (1+npos)


Instead of using all iterations in the denominator, the code is only using the iterations that resulted in a positive value.

This biases the pvalue up. In some tests I'm running, I'm getting .12 instead of .08, which is quite a difference.

I'm wondering if someone can explain this different scaling factor and provide a reference. I've seen many examples of bootstrapping like the Faraway above, but never this.

My guess is that the justification may be that the bootstrapped test statistic here should always be positive because it should be asymptotically chi square distributed. Even so, using that to correct the p-value in such a strong way makes me uncomfortable without more justification (especially since it departs from what appears to be standard practice).