There's a fairly subtle distinction here between fitting reduced models, which is what
PBmodcomp() does and what you want to do to get p-values etc., and fitting models with varying values of a focal parameter, which is what
confint() and friends do.
- Fitting reduced and full models (that differ in their fixed effects) via REML and comparing them does not work; that is, if we set one of the parameters to zero and refit, we get a structurally different model that doesn't correspond to the full model we started with
- Taking a REML model and re-evaluating/refitting it if we change the values of the focal (and other) parameters but don't set them specifically to zero is fine; this doesn't change the structure of the model. (Although I see that we do convert the REML fit to ML when profiling over the fixed effects in
lme4 ... hmmm ...)
In particular, if you're looking to compute confidence intervals for a
[g]lmer fit by parametric bootstrapping, you should use
confint(., method="boot") rather than
pbkrtest::PBmodcomp. This doesn't actually do anything with a reduced model; rather, it simulates values from the full model and re-fits.