R: Bootstrapped binary mixed-model logistic regres

2020-05-29 05:57发布

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

2条回答
Juvenile、少年°
2楼-- · 2020-05-29 06:58

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
## 
查看更多
祖国的老花朵
3楼-- · 2020-05-29 07:02

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')])))
查看更多
登录 后发表回答