lme4 (version 1.0-6)

bootMer: Model-based (Semi-)Parametric Bootstrap for Mixed Models

Description

Perform model-based (Semi-)parametric bootstrap for mixed models.

Usage

bootMer(x, FUN, nsim = 1, seed = NULL, use.u = FALSE,
    type = c("parametric", "semiparametric"),
    verbose = FALSE, .progress = "none", PBargs = list())

Arguments

x

fitted *lmer() model, see lmer, glmer, etc.

FUN

a function(x), computating the statistic of interest, which must be a numeric vector, possibly named.

nsim

number of simulations, positive integer; the bootstrap $B$ (or $R$).

seed

optional argument to set.seed.

use.u

logical, indicating whether the spherized random effects should be simulated / bootstrapped as well. If TRUE, they are not changed, and all inference is conditional on these values. If FALSE, new normal deviates are draw

type

character string specifying the type of bootstrap, "parametric" or "semiparametric"; partial matching is allowed. Note that the semiparametric bootstrap is currently an experimental feature, and therefore may not be stabl

verbose

logical indicating if progress should print output

.progress

character string - type of progress bar to display. Default is "none"; the function will look for a relevant *ProgressBar function, so "txt" will work in general; "tk" is available if the

PBargs

a list of additional arguments to the progress bar function (the package authors like list(style=3)).

Value

• an object of S3 class "boot", compatible with boot package's boot() result.

Details

The semi-parametric variant is not yet implemented, and we only provide a method for lmer and glmer results.

The working name for bootMer() was simulestimate(), as it is an extension of simulate, but we want to emphasize its potential for valid inference.

• Ifuse.uisFALSEandtypeis"parametric", each simulation generates new values of both thespherizedrandom effects$u$and the i.i.d. errors$\epsilon$, usingrnorm()with parameters corresponding to the fitted modelx.
• Ifuse.uisTRUEandtype=="parametric", only the i.i.d. errors (or, for GLMMs, response values drawn from the appropriate distributions) are resampled, with the values of$u$staying fixed at their estimated values.
• Ifuse.uisTRUEandtype=="semiparametric", the i.i.d. errors are sampled from the distribution of (response) residuals. (For GLMMs, the resulting sample will no longer have the same properties as the original sample, and the method may not make sense; a warning is generated.) Note that the semiparametric bootstrap is currently an experimental feature, and therefore may not be stable.
• The case whereuse.uisFALSEandtype=="semiparametric"is not implemented; Morris (2002) suggests that resampling from the estimated values of$u$is not good practice.

References

Davison, A.C. and Hinkley, D.V. (1997) Bootstrap Methods and Their Application. Cambridge University Press.

Morris, J. S. (2002). The BLUPs Are Not best When It Comes to Bootstrapping. Statistics & Probability Letters 56(4): 425--430. doi:10.1016/S0167-7152(02)00041-X.

See Also

• For inference, including confidence intervals,profile-methods.
• boot(), and thenboot.cifrom packageboot.

Examples

fm01ML <- lmer(Yield ~ 1|Batch, Dyestuff, REML = FALSE)
## see ?"profile-methods"
mySumm <- function(.) { s <- sigma(.)
    c(beta =getME(., "beta"), sigma = s, sig01 = unname(s * getME(., "theta"))) }
(t0 <- mySumm(fm01ML)) # just three parameters
## alternatively:
mySumm2 <- function(.) {
c(beta=fixef(.),sigma=sigma(.),sig01=unlist(VarCorr(.)))
}

set.seed(101)
## 3.8s (on a 5600 MIPS 64bit fast(year 2009) desktop "AMD Phenom(tm) II X4 925"):
system.time( boo01 <- bootMer(fm01ML, mySumm, nsim = 100) )

## to "look" at it
require("boot") ## a recommended package, i.e. *must* be there
boo01
## note large estimated bias for sig01
## (~30\% low, decreases _slightly_ for nsim = 1000)

## extract the bootstrapped values as a data frame ...
head(as.data.frame(boo01))

## ------ Bootstrap-based confidence intervals ------------

## intercept
(bCI.1 <- boot.ci(boo01, index=1, type=c("norm", "basic", "perc")))# beta

## Residual standard deviation - original scale:
(bCI.2  <- boot.ci(boo01, index=2, type=c("norm", "basic", "perc")))
## Residual SD - transform to log scale:
(bCI.2l <- boot.ci(boo01, index=2, type=c("norm", "basic", "perc"),
                   h = log, hdot = function(.) 1/., hinv = exp))

## Among-batch variance:
(bCI.3 <- boot.ci(boo01, index=3, type=c("norm", "basic", "perc")))# sig01

## Graphical examination:
plot(boo01,index=3)

## Check stored values from a longer (1000-replicate) run:
load(system.file("testdata","boo01L.RData",package="lme4"))
plot(boo01L,index=3)