The lme4 and nlme packages have made
fitting nested linear mixed-effects models quite easy. Using the
functionality of these packages we can easily use maximum likelihood or
restricted maximum likelihood to fit our model and conduct inference
using our parametric toolkit. In practice, the assumptions of our model
are often violated to such a degree that leads to biased estimators and
incorrect standard errors. In these situations, resampling methods such
as the bootstrap can be used to obtain consistent estimators of the bias
and standard errors for inference. lmeresampler provides an
easy way to bootstrap nested linear-mixed effects models using the
parametric, residual, cases, wild, or random effects block (REB)
bootstraps fit using either lme4 or nlme. The
output from lmeresampler is an lmeresamp
object. Moreover, the lmeresampler package has its own
print(), summary(), confint(),
and plot() functions. All of these functions follow the
syntax of generic print(), summary(),
confint(), and plot() functions; thus, the
commands may be used with lmeresamp objects.
You can install the latest released version from CRAN with:
install.packages("lmeresampler")or the latest development version from GitHub:
if(!require(devtools)) install.packages("devtools")
devtools::install_github("aloy/lmeresampler")Below is a short example of how to execute a parametric bootstrap for
models fit either by lme4 or by nlme. For a
complete guide outlining the purpose of lmeresampler, the
structure of lmeresamp objects, each of the bootstrap
types, how to implement parallelization with bootstrap()
calls, directions for future developers, and examples, please see the
package vignette, available on the lmeresampler
website.
First, fit the necessary model(s) with the LME package of your
choice. Here we show the same model fit first with lme4,
then with nlme:
library(lmeresampler)
library(lme4)
vcmodA <- lmer(mathAge11 ~ mathAge8 + gender + class + (1 | school), data = jsp728)
library(nlme)
vcmodB <- lme(mathAge11 ~ mathAge8 + gender + class, random = ~1|school, data = jsp728)Now, to perform a parametric bootstrap, first specify the model to be
used, then the function to return the parameters of interest, then
type = "parametric, and finally, the number of bootstrap
resamples to be computed:
# let's set .f = fixef to specify that we want only the fixed effects bootstrapped
# lme4
lmer_par_boot <- bootstrap(vcmodA, .f = fixef, type = "parametric", B = 100)
# nlme
lme_par_boot <- bootstrap(vcmodB, .f = fixef, type = "parametric", B = 100)The four other bootstrap types (residual, cases, wild, and REB) may be executed in a similar way, with some minor changes.
A full list of references may be found in the package vignette.