lmeresampler (version 0.2.4)

bootstrap.merMod: Bootstrap Nested Linear Mixed-Effects Models

Description

Perform various bootstrap process for nested linear mixed effects (LMEs) models including: parametric, residual, cases, wild, and REB bootstraps.

Usage

# S3 method for merMod
bootstrap(
  model,
  .f = extract_parameters,
  type,
  B,
  resample,
  reb_type,
  hccme,
  aux.dist,
  orig_data = NULL,
  .refit = TRUE,
  rbootnoise = 0
)
# S3 method for lme
bootstrap(
  model,
  .f = extract_parameters,
  type,
  B,
  resample,
  reb_type,
  hccme,
  aux.dist,
  orig_data = NULL,
  .refit = TRUE,
  rbootnoise = 0
)
bootstrap(
  model,
  .f,
  type,
  B,
  resample = NULL,
  reb_type = NULL,
  hccme = NULL,
  aux.dist = NULL,
  orig_data = NULL,
  .refit = TRUE,
  rbootnoise = 0
)

Value

The returned value is an object of class "lmeresamp". This is a list with the following elements:

• observed: the estimated values for the model parameters

• model: the fitted model object

• .f: the function call

• replicates: a \(B \times p\) data frame of bootstrap values for each of the p model parameters,

• stats: a tibble containing the observed, rep.mean (bootstrap mean), se (bootstrap standard error), and bias values for each model parameter,

• B: the number of bootstrap resamples performed

• data: the data with which the model was fit

• seed: a vector of randomly generated seeds that are used by the bootstrap

• type: the type of bootstrap executed

• call: the call to bootstrap() that the user

• message: a list of length B giving any messages generated during refitting. An entry will be NULL if no message was generated.

• warning: a list of length B giving any warnings generated during refitting. An entry will be NULL if no message was generated.

• error: a list of length B giving any errors generated during refitting. An entry will be NULL if no message was generated.

Arguments

model

The model object you wish to bootstrap.

.f

A function returning the statistic(s) of interest.

type

A character string indicating the type of bootstrap that is being requested. Possible values are "parametric", "residual", "case", "wild", or "reb" (random effect block bootstrap).

B

The number of bootstrap resamples.

resample

A logical vector specifying whether each level of the model should be resampled in the cases bootstrap. The levels should be specified from the highest level (largest cluster) of the hierarchy to the lowest (observation-level); for example for students within a school, specify the school level first, then the student level.

reb_type

Specification of what random effect block bootstrap version to implement. Possible values are 0, 1 or 2.

hccme

either "hc2" or "hc3", indicating which heteroscedasticity consistent covariance matrix estimator to use.

aux.dist

one of "mammen", "rademacher", "norm", "webb", or "gamma" indicating which auxiliary distribution to draw the errors from

orig_data

the original data frame. This should be specified if variables are transformed within the formula for glmer() or lmer() and the case bootstrap is used.

.refit

a logical value indicating whether the model should be refit to the bootstrap resample, or if the simulated bootstrap resample should be returned. Defaults to TRUE.

rbootnoise

a numeric value between 0-1 indicating the strength of technical 2-level noise added in relation to the 1-level variation (in standard deviations) during residual bootstrapping. Minuscule noise, such as rbootnoise = 0.0001, can be used to avoid errors with singular matrices when exactly the same values are replicated during the bootstrapping, or when the model being processed fails to return any 2-level variation. Currently applicable only with lme4::lmer models. The feature has been tested with 2-level random-intercept models with predictors. Defaults to 0 (i.e. the feature is not used by default).

Details

All of the below methods have been implemented for nested linear mixed-effects models fit by lmer (i.e., an lmerMod object) and lme (i.e., an lmerMod object). Details of the bootstrap procedures can be found in the help file for that specific function.

References

Carpenter, J. R., Goldstein, H. and Rasbash, J. (2003) A novel bootstrap procedure for assessing the relationship between class size and achievement. Journal of the Royal Statistical Society. Series C (Applied Statistics), 52, 431--443.

Chambers, R. and Chandra, H. (2013) A random effect block bootstrap for clustered data. Journal of Computational and Graphical Statistics, 22, 452--470.

Morris, J. S. (2002) The BLUPs are not "best" when it comes to bootstrapping. Statistics and Probability Letters, 56, 425--430.

Van der Leeden, R., Meijer, E. and Busing F. M. (2008) Resampling multilevel models. In J. de Leeuw and E. Meijer, editors, Handbook of Multilevel Analysis, pages 401--433. New York: Springer.

Bates, D., Maechler, M., Bolker, W., Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67, 1--48. doi:10.18637/jss.v067.i01.

Modugno, L., & Giannerini, S. (2015). The Wild Bootstrap for Multilevel Models. Communications in Statistics -- Theory and Methods, 44(22), 4812--4825.

See Also

• parametric_bootstrap, resid_bootstrap, case_bootstrap, reb_bootstrap, wild_bootstrap for more details on a specific bootstrap.

• bootMer in the lme4 package for an implementation of (semi-)parametric bootstrap for mixed models.

Examples

library(lme4) 
vcmodA <- lmer(mathAge11 ~ mathAge8 + gender + class + (1 | school), data = jsp728)

## you can write your own function to return stats, or use something like 'fixef'
mySumm <- function(.) { 
  s <- getME(., "sigma")
    c(beta = getME(., "beta"), sigma = s, sig01 = unname(s * getME(., "theta"))) 
}

## running a parametric bootstrap 
set.seed(1234)
boo1 <- bootstrap(model = vcmodA, .f = mySumm, type = "parametric", B = 20)

## to print results in a formatted way
print(boo1)

if (FALSE) {
## running a cases bootstrap - only resampling the schools
boo2 <- bootstrap(model = vcmodA, .f = mySumm, type = "case", B = 100, resample = c(TRUE, FALSE))

## running a cases bootstrap - resampling the schools and students within the school
boo3 <- bootstrap(model = vcmodA, .f = mySumm, type = "case", B = 100, resample = c(TRUE, TRUE))

## running a residual bootstrap
boo4 <- bootstrap(model = vcmodA, .f = mySumm, type = "residual", B = 100)

## running an REB0 bootstrap
boo5 <- bootstrap(model = vcmodA, .f = mySumm, type = "reb", B = 100, reb_typ = 0)

## Running the Wild bootstrap
boo6 <- bootstrap(model = vcmodA, .f = mySumm, type = "wild", B= 100, 
                  hccme = "hc2", aux.dist = "mammen")

## Running a bootstrap in parallel via foreach
library(foreach)
library(doParallel)
set.seed(1234)
numCores <- 2

cl <- makeCluster(numCores, type = "PSOCK") # make a socket cluster
doParallel::registerDoParallel(cl)          # how the CPU knows to run in parallel

b_parallel <- foreach(B = rep(250, 2), 
                       .combine = combine_lmeresamp,
                       .packages = c("lmeresampler", "lme4")) %dopar% {
                         bootstrap(vcmodA, .f = fixef, type = "parametric", B = B)
                       }
stopCluster(cl)

## Running a bootstrap in parallel via parLapply
cl <- makeCluster(numCores, type = "PSOCK") # make a socket cluster
doParallel::registerDoParallel(cl)          # how the CPU knows to run in parallel

boot_mod <- function(...) {
  library(lme4)
  library(lmeresampler)
  vcmodA <- lmer(mathAge11 ~ mathAge8 + gender + class + (1 | school), data = jsp728)
  bootstrap(vcmodA, .f = fixef, type = "parametric", B = 250)
  
}

result <- parLapply(cl, seq_len(2), boot_mod)
b_parallel2 <- do.call("combine_lmeresamp", result)

stopCluster(cl)
}