Bootstrap pvalues: strange adjustment -- pbkrtest PBmodcomp LRT

Question

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).

Answer 1

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).