R: Bootstrapped binary mixed-model logistic regression using bootMer() of the new lme4 package

Question

I want to use the new bootMer() feature of the new lme4 package (the developer version currently). I am new to R and don't know which function should I write for its FUN argument. It says it needs a numerical vector, but I have no idea what that function will perform. So I have a mixed-model formula which is cast to the bootMer(), and have a number of replicates. So I don't know what that external function does? Is it supposed to be a template for bootstrapping methods? Aren't bootstrapping methods already implemented in he bootMer? So why they need an external "statistic of interest"? And which statistic of interest should I use?

Is the following syntax proper to work on? R keeps on error generating that the FUN must be a numerical vector. I don't know how to separate the estimates from the "fit" and even should I do that in the first place? I can just say I am lost with that "FUN" argument. Also I don't know should I pass the mixed-model glmer() formula using the variable "Mixed5" or should I pass some pointers and references? I see in the examples that X (the first argument of bootMer() is a *lmer() object. I wanted to write *Mixed5 but it rendered an error.

Many thanks.

My code is:

library(lme4)
library(boot)

(mixed5 <- glmer(DV ~ (Demo1 +Demo2 +Demo3 +Demo4 +Trt)^2 
                 + (1 | PatientID) + (0 + Trt | PatientID)
                 , family=binomial(logit), MixedModelData4))


FUN <- function(formula) {
  fit <- glmer(DV ~ (Demo1 +Demo2 +Demo3 +Demo4 +Trt)^2 
               + (1 | PatientID) + (0 + Trt | PatientID)
               , family=binomial(logit), MixedModelData4)
  return(coef(fit))
}

result <- bootMer(mixed5, FUN, nsim = 3, seed = NULL, use.u = FALSE,
        type = c("parametric"),
        verbose = T, .progress = "none", PBargs = list())

result
FUN
fit

And the error:

Error in bootMer(mixed5, FUN, nsim = 3, seed = NULL, use.u = FALSE, type = c("parametric"),  : 
  bootMer currently only handles functions that return numeric vectors

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

I edited the code like what Ben instructed. The code ran very good but the SEs and Biases were all zero. Also do you know how to extract P values from this output (strange to me)? Should I use mixed() of afex package?

My revised code:

library(lme4)
library(boot)

(mixed5 <- glmer(DV ~ (Demo1 +Demo2 +Demo3 +Demo4 +Trt)^2 
                + (0 + Trt | PatientID)
                 , family=binomial(logit), MixedModelData4))


FUN <- function(fit) {
  fit <- glmer(DV ~ (Demo1 +Demo2 +Demo3 +Demo4 +Trt)^2 
               + (1 | PatientID) + (0 + Trt | PatientID)
               , family=binomial(logit), MixedModelData4)
  return(fixef(fit))
}

result <- bootMer(mixed5, FUN, nsim = 3)

result

-------------------------------------------------------- Update 2 -----------------------------------------------------

I also tried the following but the code generated warnings and didn't give any result.

(mixed5 <- glmer(DV ~ Demo1 +Demo2 +Demo3 +Demo4 +Trt 
                 + (1 | PatientID) + (0 + Trt | PatientID)
                 , family=binomial(logit), MixedModelData4))

FUN <- function(mixed5) {
  return(fixef(mixed5))}

result <- bootMer(mixed5, FUN, nsim = 2)

Warning message:

In bootMer(mixed5, FUN, nsim = 2) : some bootstrap runs failed (2/2)
> result

Call:
bootMer(x = mixed5, FUN = FUN, nsim = 2)

Bootstrap Statistics :
WARNING: All values of t1* are NA
WARNING: All values of t2* are NA
WARNING: All values of t3* are NA
WARNING: All values of t4* are NA
WARNING: All values of t5* are NA
WARNING: All values of t6* are NA

-------------------------------------------------------- Update 3 -----------------------------------------------------

This code as well generated warnings:

FUN <- function(fit) {
  return(fixef(fit))}

result <- bootMer(mixed5, FUN, nsim = 2)

The warnings and results:

Warning message:
In bootMer(mixed5, FUN, nsim = 2) : some bootstrap runs failed (2/2)
> result

Call:
bootMer(x = mixed5, FUN = FUN, nsim = 2)

Bootstrap Statistics :
WARNING: All values of t1* are NA
WARNING: All values of t2* are NA
WARNING: All values of t3* are NA
WARNING: All values of t4* are NA
WARNING: All values of t5* are NA
WARNING: All values of t6* are NA

Esterify

Answer 1

There are basically two (simple) confusions here.

• The first is between coef() (which returns a list of matrices) and fixef() (which returns a vector of the fixed-effect coefficients): I assume that fixef() is what you wanted, although you might want something like c(fixef(mixed),unlist(VarCorr(mixed))).
• the second is that FUN should take a fitted model object as input ...

For example:

library(lme4)
library(boot)

mixed <- glmer(incidence/size ~ period + (1|herd),
               weights=size, data=cbpp, family=binomial)

FUN <- function(fit) {
    return(fixef(fit))
}

result <- bootMer(mixed, FUN, nsim = 3)

result

## Call:
## bootMer(x = mixed, FUN = FUN, nsim = 3)
## Bootstrap Statistics :
##      original      bias    std. error
## t1* -1.398343 -0.20084060  0.09157886
## t2* -0.991925  0.02597136  0.18432336
## t3* -1.128216 -0.03456143  0.05967291
## t4* -1.579745 -0.08249495  0.38272580
## 

Klimesh

I revised the code accordingly (the update under my question) and it ran smoothly and errorless. However I have 2 question: all the biases and SEs are zero. 2: what are these lines? (I mean t1, t2...)? Raw bootstrapped resamples? Is there any way to extract estimates and P values from this output? – Esterify

(1) you should not re-run your model within FUN -- just use the fixef() call by itself. (2) t1, t2, ... in the output are the values of the fixed effect estimates. (3) You should definitely look at ?confint.merMod ... – Klimesh

Thank you so much for all these gems. I also used this: "FUN <- function(fit) { return(fixef(fit))" [see my update 2] but no estimates were generated. But I will fix it THANKS to you :) – Esterify

oh I see sorrycsorry! I have given it a fit object (instead of "fit" object) – Esterify

Uh I passed "fit" to it, but other warnings happened and the results were artifact warnings. (detailed in update 3) – Esterify

In the last example, you only did two bootstrap simulations, and both of them failed. It is not terribly unusual for some bootstrap realizations to fail in re-fitting, but if the model & data are reasonable it shouldn't be all of them. Can you re-try with a more reasonable number (say, 20? for real inference you should use 100-200 or more, but 20 should be enough for the next stage of testing). (If you have multiple cores you might also want to turn on some of the parallel functionality ...) – Klimesh

Thank you very mych dear Ben. I set the nsim to 50 and later at 1000. In both cases, the same warning happened within 10 seconds. "Warning message: In bootMer(mixed5, FUN, nsim = 1000) : some bootstrap runs failed (1000/1000)" – Esterify

OK, then; I don't see anything obviously wrong, which means the next step is a reproducible example. Can you either (1, preferably) post a reproducible example <tinyurl.com/reproducible-000> or (2, if necessary) send me your data? – Klimesh

(I had set the number of simulation too low to make debugging as short as possible but I would use the numbers you suggested in my experiment) :) Thank you again. But I don't know what causes the problem. Do I need to install the latest Developing version of BOOT package too? perhaps the CRAN version is not compatible with the new lme4. – Esterify

I can send you my data but it is confidential i cannot post it online :) – Esterify

let us continue this discussion in chat – Klimesh

Answer 2

This might be the same problem, that I reported as an issue here. At least it leads to the same, unhelpful error message and took me a while too.

That would mean you have missings in your data, which lmer ignores but which kill bootMer.

Try:

(mixed5 <- glmer(DV ~ (Demo1 +Demo2 +Demo3 +Demo4 +Trt)^2 
                 + (1 | PatientID) + (0 + Trt | PatientID)
                 , family=binomial(logit), na.omit(MixedModelData4[,c('DV','Demo1','Demo2','Demo3','Trt','PatientId')])))

Livia

Thanks a lot. However, I don't have any missing data. But I should say my data might be pathological. – Esterify

@Esterify You probably would've caught that. But now at least this simple answer is up here, for people who google the error message. I don't really get what you mean by pathological data? Do you get the same problems with a simple model with the same data? – Livia

I appreciate your act Ruben. Also I suggest you to follow Ben Bolker's instruction to me, and upload your data, or a similar data to an online repo, so that the package developers can work on the bug. (See our conversation with Ben Bolker, and read his instructions).... Pathological data means data suffering from complete / partial separation, or serious multicollinearity, or both. – Esterify

@vic Ah ok. If you follow the link to the GitHub issue, you'll see that I already created a reproducible example with sample data - it really takes just a few missings in otherwise ok data. – Livia