R: bootstrapped mixed model binary logistic regression

Question

I need to bootstrap my mixed model binary logistic regression. The model itself works fine (and is approved and corrected by an expert friend), but the bootstrapped version is buggy. The bootstrapped version was previously approved by another expert friend (in CrossValidated but later mods removed my post saying it does not belong on CrossValidated). But the same code happened to work for a simple fixed-effects multiple logistic regression (although in that case too there were lots of warnings similar to the warnings here [except this single warning which is for the lmer() function: "In mer_finalize(ans) : false convergence (8)").

Could you please let me know where the error resides and how to debug it?

Many thanks.

My code is (I temporarily kept the replicate numbers too low to debug the code):

library(boot)
library(lme4)

mixedGLM <- function(formula, data, indices) {
        d <- data[indices, ]
        (fit <- lmer(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2 
                     + (1 | PatientID) + (0 + Trt | PatientID)
                     , family=binomial(logit), d))
        return(coef(fit))
      }

results <- boot(data=MixedModelData4 , statistic = mixedGLM, R= 2, formula= DV~Demo1 +Demo2 +Demo3 +Demo4 +Trt)

. . . My errors are:

Error in t.star[r, ] <- res[[r]] : 
  incorrect number of subscripts on matrix
In addition: Warning messages:
1: In mer_finalize(ans) : false convergence (8)
2: glm.fit: algorithm did not converge 
3: glm.fit: fitted probabilities numerically 0 or 1 occurred 
4: glm.fit: fitted probabilities numerically 0 or 1 occurred 
5: In mer_finalize(ans) : false convergence (8) 

. . . Also could you please tell me how to make the boot() function give P values too??! It just gives beta and SE and bias and CI, but I need the P values too.

Many thanks.

--------------------------------------------------- Developing Story -----------------------------------------------------

Ok I gladly ran the nice code of Henrik. But the code did not quite finish running. First it gave this error:

Fitting 17 lmer() models:
[...
Error: pwrssUpdate did not converge in 30 iterations
In addition: Warning message:
In mixed(DV ~ (Demo1 + Demo2 + Demo3 + Demo4 + Trt)^2 + (1 | PatientID) +  :
  Due to missing values, reduced number of observations to 90
> (results2 <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2
+ results3 <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2 

Then I removed the first parentheses block and revised the syntax to this one:

results3 <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2 
                 + (0 + Trt | PatientID),
                 family=binomial(logit), data = MixedModelData4,
                 method = "PB", args.test = list(nsim = 2))

This time the test passed the first step (fitting the models) but failed at obtaining P values, again giving the same errors and warnings:

Fitting 17 lmer() models:
[.................]
Obtaining 16 p-values:
[....
Error: pwrssUpdate did not converge in 30 iterations
In addition: Warning messages:
1: In mixed(DV ~ (Demo1 + Demo2 + Demo3 + Demo4 + Trt)^2 + (0 + Trt |  :
  Due to missing values, reduced number of observations to 90
2: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations
3: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations
4: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations
5: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations
6: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations

I have no idea how to debug it, or if the problem is my dataset? I should add that my dataset is fully mean-centered (all variables). The DV is only negated (since mean centering disallowed R to work and negating would do the same for a binary outcome).

---------------------------------------------------------- Update -------------------------------------------------------------

I changed the PB value of METHOD to LRT (as Henrik recommended) and the process of fitting the models finished but the process of obtaining the P values didn't start:

> results4 <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2 
+                   + (0 + Trt | PatientID),
+                   family=binomial(logit), data = MixedModelData4,
+                   method = "LRT", args.test = list(nsim = 2))
Fitting 17 lmer() models:
[.................]
Warning message:
In mixed(DV ~ (Demo1 + Demo2 + Demo3 + Demo4 + Trt)^2 + (0 + Trt |  :
  Due to missing values, reduced number of observations to 90

It turned out the P values are not obtained by bootstrapping when LRT is being used. Therefore, the results were already ready (although non-bootstrapped).

Answer 1

If you want p-values from a GLMM with a parametric bootstrap you can use function mixed from package afex which obtains them via pbkrtest::PBmodcomp:

library(afex)
results <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2 
                     + (1 | PatientID) + (0 + Trt | PatientID),
                     family=binomial(logit), data = d,
                     method = "PB", args.test = list(nsim = 1000))

You could even first define a local cluster (i.e., use multiple cores):

cl <- makeCluster(rep("localhost", 4))
results <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2 
                     + (1 | PatientID) + (0 + Trt | PatientID),
                     family=binomial(logit), data = d,
                     method = "PB", args.test = list(nsim = 1000, cl = cl))

It is probably the best to install the development versions of all three packages (as the current version of pbkrtest is designed for lme4 1.0 which is not yet on cran):

• afex from r-forge: install.packages("afex", repos="http://R-Forge.R-project.org")
• pbkrtest from: http://people.math.aau.dk/~sorenh/software/pbkrtest/devel/
• lme4 from github: https://github.com/lme4/lme4