mixed parametric bootstrapping with optimx::nlminb

[Peter Vavra]

I’m having trouble getting a glmm fitted using optimix::nlminb (the default optimizer doesn’t converge) and calculating the p-values using parametric bootstrapping, using the mixed() function. Interestingly, using likelihood ratio tests works just fine with the “same” code.

this one works without any errors:

m6.LRT <- mixed( decision ~ expectation.mean * expectation.var * c.offer + (1 +  c_c.offer  | subject.id) + (0 + c_expectation.mean1 | subject.id),
                   family="binomial",
                   data=data.merged,
                   method="LRT",
                   cl=my.cluster,
                   control=glmerControl(optimizer = "optimx",calc.derivs = FALSE, 
                                        optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE))

and this one (I only changed method="LRT" to method="PB") does not work:

m6.PB <- mixed( decision ~ expectation.mean * expectation.var * c.offer + (1 +  c_c.offer  | subject.id) + (0 + c_expectation.mean1 | subject.id),
                   family="binomial",
                   data=data.merged,
                   method="PB",
                   cl=my.cluster,
                   control=glmerControl(optimizer = "optimx",calc.derivs = FALSE, 
                                        optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE))

Instead, I get the following error message:
Error in optwrap(object@optinfo$optimizer, ff, x0, lower = lower, control = control$optCtrl, : must specify 'method' explicitly for optimx

Any ideas what I can do to make it work?

[henrik]

Does it help to load optimx at the clusters before running the code?
clusterEvalQ(my.cluster, library(optimx))

If not, please let me know and I will investigate this further.

[Peter Vavra]

It’s not an issue with the multi-core/cluster environment.

I’ve tried the suggested code, and also simply commented out the cl=my.cluster line – still no success.

However, running it on a single core gave a clue which might be helpful for debugging this: In my case, 8 models need to be fit. That step seems to work, as all eight “progress dots” (the [........]) are there and the error message pops up when the function proceeds to calculate the p-values, i.e. just after

Obtaining 7 p-values: 
[

[henrik]

The problem appears to be in the function calculating the parametric bootstrap pbkrtest::PBmodcomp:

require(lme4)
require(optimx)
require(pbkrtest)
## Linear mixed effects model:
data(beets, package="pbkrtest")
sug   <- lmer(sugpct ~ block + sow + harvest + (1|block:harvest), data=beets, REML=FALSE, control =lmerControl(optimizer = "optimx",optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
sug.h <- update(sug, .~. -harvest)
PBmodcomp(sug, sug.h, nsim=50)
# Error in optwrap(object@optinfo$optimizer, ff, x0, lower = lower, control = control$optCtrl,  : 
#   must specify 'method' explicitly for optimx

I will contact the author and let you know what comes out of it. I can do little at the time.

[henrik]

After contacting pbkrtest author Søren Højsgaard he told me that it was actually a bug in lme4::refit. So I looked into it and wrote a bug fix that was just merged into lme4: https://github.com/lme4/lme4/pull/421

So either you wait until the new version of lme4 appears on CRAN (and I have no idea when this might be) or you install the development version: library("devtools"); install_github("lme4/lme4",dependencies=TRUE)
Note that installing the development version might require other dependencies depending on your OS (e.g., RTools on Windows or XCode on Mac).

[Author]

Posts