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).
 A: There is a now a note in the package documentation that answers this:
Note
It can happen that some values of the LRT statistic in the reference distribution are negative. When
this happens one will see that the number of used samples (those where the LRT is positive) are
reported (this number is smaller than the requested number of samples).
In theory one can not have a negative value of the LRT statistic but in practice on can: We speculate
that the reason is as follows:  We simulate data under the small model and fit both the small and
the large model to the simulated data.  Therefore the large model represents - by definition - an
overfit;  the model has superfluous parameters in it.  Therefore the fit of the two models will for
some simulated datasets be very similar resulting in similar values of the log-likelihood.  There is
no guarantee that the the log-likelihood for the large model in practice always will be larger than
for the small (convergence problems and other numerical issues can play a role here).
To look further into the problem, one can use the
PBrefdist()
function for simulating the refer-
ence distribution (this reference distribution can be provided as input to
PBmodcomp()
). Inspection
sometimes reveals that while many values are negative, they are numerically very small. In this case
one may try to replace the negative values by a small positive value and then invoke
PBmodcomp()
to get some idea about how strong influence there is on the resulting p-values.  (The p-values get
smaller this way compared to the case when only the originally positive values are used).
