Chuan Tang - 1 year ago 68
R Question

# R: obtain coefficients&CI from bootstrapping mixed-effect model results

The working data looks like:

``````set.seed(1234)
df <- data.frame(y = rnorm(1:30),
fac1 = as.factor(sample(c("A","B","C","D","E"),30, replace = T)),
fac2 = as.factor(sample(c("NY","NC","CA"),30,replace = T)),
x = rnorm(1:30))
``````

The
`lme`
model is fitted as:

``````library(lme4)
mixed <- lmer(y ~ x + (1|fac1) + (1|fac2), data = df)
``````

I used
`bootMer`
to run the parametric bootstrapping and I can successfully obtain the coefficients (intercept) and SEs for fixed&random effects:

``````mixed_boot_sum <- function(data){s <- sigma(data)
c(beta = getME(data, "fixef"), theta = getME(data, "theta"), sigma = s)}

mixed_boot <- bootMer(mixed, FUN = mixed_boot_sum, nsim = 100, type = "parametric", use.u = FALSE)
``````

My first question is how to obtain the coefficients(slope) of each individual levels of the two random effects from the bootstrapping results
`mixed_boot`
?

I have no problem extracting the coefficients(slope) from
`mixed`
model by using
`augment`
function from
`broom`
package, see below:

``````library(broom)
mixed.coef <- augment(mixed, df)
``````

However, it seems like
`broom`
can't deal with
`boot`
class object. I can't use above functions directly on
`mixed_boot`
.

I also tried to modify the
`mixed_boot_sum`
`mmList`
( I thought this would be what I am looking for), but R complains as：

``````Error in bootMer(mixed, FUN = mixed_boot_sum, nsim = 100, type = "parametric",  :
bootMer currently only handles functions that return numeric vectors
``````

Furthermore, is it possible to obtain CI of both fixed&random effects by specifying
`FUN`
as well?

Now, I am very confused about the correct specifications for the
`FUN`
in order to achieve my needs. Any help regarding to my question would be greatly appreciated!

My first question is how to obtain the coefficients(slope) of each individual levels of the two random effects from the bootstrapping results mixed_boot ?

I'm not sure what you mean by "coefficients(slope) of each individual level". `broom::augment(mixed, df)` gives the predictions (residuals, etc.) for every observation. If you want the predicted coefficients at each level I would try

``````mixed_boot_coefs <- function(fit){
unlist(coef(fit))
}
``````

which for the original model gives

``````mixed_boot_coefs(mixed)
## fac1.(Intercept)1 fac1.(Intercept)2 fac1.(Intercept)3 fac1.(Intercept)4
##        -0.4973925        -0.1210432        -0.3260958         0.2645979
## fac1.(Intercept)5           fac1.x1           fac1.x2           fac1.x3
##        -0.6288728         0.2187408         0.2187408         0.2187408
##           fac1.x4           fac1.x5 fac2.(Intercept)1 fac2.(Intercept)2
##         0.2187408         0.2187408        -0.2617613        -0.2617613
##  ...
``````

If you want the resulting object to be more clearly named you can use:

``````flatten <- function(cc) setNames(unlist(cc),
outer(rownames(cc),colnames(cc),
function(x,y) paste0(y,x)))

mixed_boot_coefs <- function(fit){
unlist(lapply(coef(fit),flatten))
}
``````

When run through `bootMer`/`confint`/`boot::boot.ci` these functions will give confidence intervals for each of these values (note that all of the slopes `facW.xZ` are identical across groups because the model assumes random variation in the intercept only). In other words, whatever information you know how to extract from a fitted model (conditional modes/BLUPs [`ranef`], predicted intercepts and slopes for each level of the grouping variable [`coef`], parameter estimates [`fixef`, `getME`], random-effects variances [`VarCorr`], predictions under specific conditions [`predict`] ...) can be used in `bootMer`'s `FUN` argument, as long as you can flatten its structure into a simple numeric vector.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download