Package {bootmlm}

bootmlm: A package for bootstrap resampling with multilevel data.

Description

Currently, the bootMer() function in the lme4 package only implements the parametric bootstrap and a limited version of semiparametric (or residual) bootstrap. The bootmlm package provides the function bootstrap_mer, which performs various parametric, residual, and case bootstrap resampling for fitted model objects with the lme4 package.

Limitations

• Currently only support multilevel models (a.k.a. linear mixed-effects models) fitted by lmer() with the lme4 package. Support for categorical outcome (i.e., generalized linear mixed-effects models fitted by glmer()) and for models fitted with the nlme package will be added in the future.

• Random effect block bootstrap (type = 'reb') and case bootstrap (type = 'case') only support two-level models.

• Bias-corrected and accelerated bootstrap (using empinf_mer()) only supports two-level models.

Author(s)

Maintainer: Mark Lai marklhc@gmail.com [copyright holder]

Authors:

• Mark Lai marklhc@gmail.com [copyright holder]

See Also

Useful links:

• https://github.com/marklhc/bootmlm

• https://marklhc.github.io/bootmlm/

• Report bugs at https://github.com/marklhc/bootmlm/issues

Run Various Bootstrap for Mixed Models.

Description

Run multilevel parametric, residual, and case bootstrap with different options

Usage

bootstrap_mer(
  x,
  FUN,
  nsim = 1,
  type = c("parametric", "residual", "residual_cgr", "residual_trans", "reb", "case"),
  corrected_trans = FALSE,
  lv1_resample = FALSE,
  reb_scale = FALSE,
  .progress = FALSE,
  verbose = FALSE,
  ...
)

Arguments

Details

bootstrap_mer performs different bootstrapping methods to fitted model objects using the lme4 package. Currently, only models fitted using lmer is supported.

Value

An object of S3 class "boot", compatible with boot package's boot(). It contains the following components:

See the documentation in for link[boot]{boot}() for the other components.

References

Carpenter, J. R., Goldstein, H., & 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. https://doi.org/10.1111/1467-9876.00415

Chambers, R., & Chandra, H. (2013). A random effect block bootstrap for clustered data. Journal of Computational and Graphical Statistics, 22(2), 452–470. https://doi.org/10.1080/10618600.2012.681216

Davison, A. C. and Hinkley, D. V. (1997). Bootstrap methods and their application. Cambridge, UK: Cambridge University Press.

Morris, J. S. (2002). The BLUPs are not "best" when it comes to bootstrapping. Statistics & Probability Letters, 56(4), 425–430. https://doi.org/10.1016/S0167-7152(02)00041-X

Van der Leeden, R., Meijer, E., & Busing, F. M. T. A. (2008). Resampling multilevel models. In J. de Leeuw & E. Meijer (Eds.), Handbook of multilevel Analysis (pp. 401–433). New York, NY: Springer.

See Also

• boot for single-level bootstrapping,

• bootMer for parametric and semi-parametric bootstrap implemented in lme4, and

• boot.ci for getting bootstrap confidence intervals and plot.boot for plotting the bootstrap distribution.

Examples

library(lme4)
fm01ML <- lmer(Yield ~ (1 | Batch), Dyestuff, REML = FALSE)
mySumm <- function(x) {
  c(getME(x, "beta"), sigma(x))
}
# Covariance preserving residual bootstrap
boo01 <- bootstrap_mer(fm01ML, mySumm, type = "residual", nsim = 100)
# Plot bootstrap distribution of fixed effect
library(boot)
plot(boo01, index = 1)
# Get confidence interval
boot.ci(boo01, index = 2, type = c("norm", "basic", "perc"))
# BCa using influence values computed from `empinf_mer`
boot.ci(boo01, index = 2, type = "bca", L = empinf_mer(fm01ML, mySumm, 2))

Bootstrap confidence intervals for Two-Level Mixed Models

Description

This is a wrapper for getting CIs for multiple parameters after running bootstrap_mer, to be consistent with a similar method for the bootMer class.

Usage

## S3 method for class 'boot'
confint(
  object,
  parm,
  level = 0.95,
  type = c("norm", "basic", "perc", "bca"),
  L = NULL,
  ...
)

Arguments

Value

A matrix (or vector) with columns giving lower and upper confidence limits for each parameter.

Examples

library(lme4)
fm01ML <- lmer(Yield ~ (1 | Batch), Dyestuff, REML = FALSE)
mySumm <- function(x) {
  c(getME(x, "beta"), sigma(x))
}

# residual bootstrap
boo_resid <- bootstrap_mer(fm01ML, mySumm, type = "residual", nsim = 100)
confint(boo_resid, type = "bca", L = empinf_merm(fm01ML, mySumm))

Deviance Function for Multilevel Models

Description

Creates a deviance function for a fitted model object, using \theta and sigma as the parameters.

Usage

devfun_mer(x)

devfun_mer2(x)

Arguments

Details

The built-in function(s) in lme4 for generating deviance function are not exported and rely on C++ code. This function is mainly used to obtain Hessian and asymptotic covariance matrix of the random effects.

Value

A deviance function with one argument th_sig that takes input of a numeric vector corresponding to the some estimated values of \theta and \sigma. For devfun_mer2, a function with one argument theta that profiles out \sigma and only takes input of elements for \theta.

References

Bates, D., Maechler, M., Bolker, B. M., & Walker, S. C. Fitting linear mixed-effects models using lme4. Retrieved from https://www.jstatsoft.org/article/view/v067i01

Implementation in lme4pureR: https://github.com/lme4/lme4pureR/blob/master/R/pls.R

Examples

library(lme4)
fm01ML <- lmer(Yield ~ (1 | Batch), Dyestuff, REML = FALSE)
dd <- devfun_mer(fm01ML)
# Asymptotic variance-covariance matrix of (theta, sigma):
2 * solve(numDeriv::hessian(dd, c(fm01ML@theta, sigma(fm01ML))))

Empirical Influence Values for Two-Level Mixed Models

Description

This function calculates the empirical influence values for a statistic in a given fitted model object using the delete-m_j jackknife.

Usage

empinf_mer(x, FUN, index = 1)

empinf_merm(x, FUN)

Arguments

Details

empinf_mer computes non-parametric influence function of models fitted using lmer by deleting one cluster at a time. See van der Leeden, Meijer, and Busing (2008, pp. 420–422) for more information. Whereas empinf_mer computes influence values for a specified position (as specified with the index argument) of the output of FUN, empinf_merm computes influence values for every element in FUN(x).

Value

A numeric vector with length equals to number of clusters of x containing the weighted influence value of each cluster.

References

Van der Leeden, R., Meijer, E., & Busing, F. M. T. A. (2008). Resampling multilevel models. In J. de Leeuw & E. Meijer (Eds.), Handbook of multilevel Analysis (pp. 401–433). New York, NY: Springer.

Examples

library(lme4)
fm01ML <- lmer(Yield ~ (1 | Batch), Dyestuff, REML = FALSE)
# Define function for intraclass correlation
icc <- function(x) 1 / (1 + 1 / getME(x, "theta")^2)
empinf_mer(fm01ML, icc)
empinf_mer(fm01ML, fixef)

Synthetic Pupil Popularity Dataset

Description

A synthetic two-level dataset with pupils nested within schools, generated to mimic the structure and parameters of the popular dataset from Hox (2010). It can be used to demonstrate multilevel bootstrap methods without depending on external data sources.

Usage

pop_syn

Format

A data frame with 2000 rows and 5 variables:

pupil

Pupil identification number within school (integer).

school

School identification number, 1–100 (integer).

popular

Pupil popularity score on a 0–10 scale (numeric).

sex

Pupil sex: 0 = boy, 1 = girl (integer).

texp

Teacher experience in years (numeric).

Details

The dataset was generated by fitting a two-level random intercept model to the original popular data and simulating new data from the estimated parameters:

• Fixed effects: intercept = 3.56, sex = 0.84, texp = 0.093

• School-level random intercept SD: 0.69

• Residual SD: 0.68

Continuous predictions were rounded to the nearest integer and clamped to the [0, 10] range to match the original Likert-type scale.

Source

Simulated from parameters estimated from the popular dataset in: Hox, J. J. (2010). Multilevel analysis: Techniques and applications (2nd ed.). Routledge. Data available at https://stats.oarc.ucla.edu/stat/stata/examples/mlm_ma_hox/popular.dta.

Profile Likelihood Confidence Interval for Intraclass Correlation

Description

Compute confidence intervals for the intraclass correlation of a model fit of class merMod-class.

Usage

prof_ci_icc(x, level = 0.95, eps_max = 1e+06)

Arguments

Details

This function uses the uniroot function and determine the lower and upper limit by evaluating the profile deviance (for ML) or the -2 profile REML criterion. It works by obtaining the interval for \theta = \tau / \sigma and transforming the two limits.

The resulting interval is bounded by zero for its lower limit.

Value

A named numeric vector with two values showing the lower and upper limit of the intraclass correlation.

Examples

library(lme4)
fm01ML <- lmer(Yield ~ (1 | Batch), Dyestuff, REML = FALSE)
# Profile likelihood 95% CI for the ICC
prof_ci_icc(fm01ML)
# 90% CI
prof_ci_icc(fm01ML, level = 0.90)

Score Functions and Case-wise Derivatives

Description

Score Functions and Case-wise Derivatives

Usage

scores_mer(x, level = 2)

Arguments

Value

A numeric matrix of score contributions. If level = 2, rows correspond to level-2 clusters (aggregated); if level = 1, rows correspond to level-1 observations. Columns correspond to model parameters in order: fixed effects, variance components, residual variance.

Examples

library(lme4)
fm01ML <- lmer(Yield ~ (1 | Batch), Dyestuff, REML = FALSE)
# Cluster-level (level-2) scores
scores_mer(fm01ML, level = 2)
# Observation-level (level-1) scores
scores_mer(fm01ML, level = 1)

Get the square root of a matrix using eigenvalue decomposition, and solve for a linear system

Description

Solve for x in the linear system Ax = b, where A is a symmetric matrix square root of M with eigenvalue decomposition such that AA = M.

Usage

solve_eigen_sqrt(M, b)

Arguments

Value

A numeric vector or matrix x satisfying Ax = b, where A is the symmetric square root of M.

Asymptotic Covariance Matrix for Cholesky Factor of Random Effects

Description

Asymptotic Covariance Matrix for Cholesky Factor of Random Effects

Usage

vcov_theta(x)

Arguments

Value

A symmetric matrix giving the asymptotic covariance matrix of the Cholesky factor \theta of the random-effects covariance.

See Also

vcov_vc

Examples

library(lme4)
fm01ML <- lmer(Yield ~ (1 | Batch), Dyestuff, REML = FALSE)
# Asymptotic covariance of the Cholesky factor theta
vcov_theta(fm01ML)

Asymptotic Covariance Matrix for Random Effects

Description

Return the asymptotic covariance matrix of random effect standard deviations (or variances) for a fitted model object, using the Hessian evaluated at the (restricted) maximum likelihood estimates.

Usage

vcov_vc(x, sd_cor = TRUE, print_names = TRUE)

Arguments

Details

Although it's easy to obtain the Hessian for \theta, the relative Cholesky factor, in lme4, there is no easy way to obtain the Hessian for the variance components. This function uses devfun_mer() to obtain the Hessian (H) of variance components (or standard deviations, SD), and then obtain the asymptotic covariance matrix as -2 H^{-1}.

Value

A (q + 1) * (q + 1) symmetric matrix of the covariance matrix of (\tau, \sigma) (if sd_cor = TRUE) or (\tau^2, \sigma^2) (if sd_cor = FALSE), where q is the the number of estimated random-effects components (excluding \sigma). For example, for a model with random slope, \tau = (intercept SD, intercept-slope correlation, slope SD).

See Also

vcov.merMod for covariance matrix of fixed effects, confint.merMod for confidence intervals of all parameter estimates, and devfun_mer for the underlying function to produce the deviance function.

Examples

library(lme4)
data(Orthodont, package = "nlme")
fm1 <- lmer(distance ~ age + (age | Subject), data = Orthodont)
vc <- VarCorr(fm1)
# Standard deviation only
print(vc, comp = c("Std.Dev"))
# Asymptotic variance-covariance matrix of (tau, sigma):
vcov_vc(fm1, sd_cor = TRUE)


# Compare with (parametric) bootstrap results :
get_sdcor <- function(x) {
  as.data.frame(lme4::VarCorr(x), order = "lower.tri")[ , "sdcor"]
}
boo <- bootstrap_mer(fm1, get_sdcor, type = "parametric", nsim = 200L)
# There might be failures in some resamples
cov(boo$t, use = "complete.obs")