#' @rdname bootstrap
#' @export
bootstrap.lme <- function(model, .f = extract_parameters, type, B, resample, reb_type, hccme, aux.dist, orig_data = NULL, .refit = TRUE, rbootnoise = 0){
if (rbootnoise != 0) {
stop("'rbootnoise' implemented currently only for LME4 models. Do not define or use default 0.")
}
switch(type,
parametric = parametric_bootstrap.lme(model, .f, B, .refit),
residual = resid_bootstrap.lme(model, .f, B, .refit, rbootnoise),
case = case_bootstrap.lme(model, .f, B, resample, orig_data, .refit),
reb = reb_bootstrap.lme(model, .f, B, reb_type, .refit),
wild = wild_bootstrap.lme(model, .f, B, hccme, aux.dist, .refit))
}
#' @rdname parametric_bootstrap
#' @export
#' @importFrom nlmeU simulateY
parametric_bootstrap.lme <- function(model, .f, B, .refit = TRUE){
# Match function
if (.refit) {
.f <- match.fun(.f)
t0 <- .f(model)
}
# Extract fixed effects
# model.fixef <- nlme::fixef(model)
ystar <- nlmeU::simulateY(model, nsim = B)
row.names(ystar) <- 1:model$dims$N
ystar <- data.frame(ystar)
if(!.refit) return(ystar)
refits <- refit_lme(ystar = ystar, model = model, .f = .f)
.bootstrap.completion(model, tstar = refits$tstar, B, .f, type = "parametric", warnings = refits$warnings)
}
#' @rdname case_bootstrap
#' @export
case_bootstrap.lme <- function(model, .f, B, resample, orig_data = NULL, .refit = TRUE){
data <- model$data
# data$.id <- seq_len(nrow(data))
clusters <- c(names(model$groups), ".id")
if(length(clusters) != length(resample))
stop("'resample' is not the same length as the number of grouping variables.
Please specify whether to resample the data at each level of grouping.")
refits <- purrr::map(integer(B), function(x)
.resample_refit.cases(model = model, .f = .f, dat = data, cluster = clusters, resample = resample, .refit = .refit))
tstar <- purrr::map(refits, ~.x)
if(!.refit) return(tstar)
.bootstrap.completion(model, tstar = tstar, B, .f, type = "case",
warnings = collect_warnings(.resample_refit.cases))
}
#' @rdname reb_bootstrap
#' @inheritParams bootstrap
#' @export
reb_bootstrap.lme <- function(model, .f, B, reb_type, .refit = TRUE){
if(missing(reb_type)){
reb_type <- 0
warning("'reb_type' unspecified, performing REB 0 bootstrap")
}
if(ncol(model$groups) > 1) {
stop("The REB bootstrap has not been adapted for 3+ level models.")
}
if(!.refit & reb_type == 2) {
stop(".refit == FALSE is not available with reb_type = 2.")
}
.f <- match.fun(.f)
# Set up for bootstrapping
setup <- .setup.lme(model, type = "reb", reb_type = reb_type)
# Generate bootstrap responses
y.star <- replicate(
n = B,
.resample.reb.lme(
Xbeta = setup$Xbeta,
Zlist = setup$Zlist,
Uhat = setup$b,
estar.vec = as.numeric(setup$e),
flist = setup$flist,
levs = setup$levs
)
)
ystar <- as.data.frame(y.star)
if(!.refit) return(ystar)
# Extract bootstrap statistics
if(reb_type == 2) .f <- extract_parameters.lme
refits <- refit_lme(ystar = ystar, model = model, .f = .f)
tstar <- refits$tstar
# Extract original statistics
t0 <- .f(model)
# Postprocessing for REB/2
if(reb_type == 2)
tstar <- .postprocess.reb2(t0, tstar, nbeta = length(nlme::fixef(model)), B = B)
# Format for return
.bootstrap.completion(model, tstar = tstar, B, .f,
type = paste("reb", reb_type, sep = ""),
warnings = refits$warnings)
}
#' @rdname resid_bootstrap
#' @inheritParams bootstrap
#' @export
resid_bootstrap.lme <- function(model, .f, B, .refit = TRUE, rbootnoise){
.f <- match.fun(.f)
setup <- .setup.lme(model, type = "residual")
ystar <- as.data.frame(
replicate(
n = B,
.resample.cgr.lme(
b = setup$b,
e = setup$e,
Xbeta = setup$Xbeta,
Zlist = setup$Zlist,
vclist = setup$vclist,
sig0 = setup$sig0
)
)
)
if(!.refit) return(ystar)
refits <- refit_lme(ystar = ystar, model = model, .f = .f)
.bootstrap.completion(model, tstar = refits$tstar, B, .f, type = "residual", warnings = refits$warnings)
}
#' @rdname wild_bootstrap
#' @inheritParams bootstrap
#' @export
wild_bootstrap.lme <- function(model, .f, B, hccme = c("hc2", "hc3"),
aux.dist = c("mammen", "rademacher",
"norm", "webb", "gamma"),
.refit = TRUE){
.f <- match.fun(.f)
hccme <- match.arg(hccme)
aux.dist <- match.arg(aux.dist)
setup <- .setup.lme(model, type = "wild")
ystar <- as.data.frame(
replicate(
n = B,
.resample.wild(
Xbeta = setup$Xbeta,
mresid = setup$mresid,
.hatvalues = setup$.hatvalues,
hccme = hccme,
aux.dist = aux.dist,
n.lev = setup$n.lev,
flist = setup$flist
)
)
)
if(!.refit) return(ystar)
refits <- refit_lme(ystar = ystar, model = model, .f = .f)
.bootstrap.completion(model, tstar = refits$tstar, B, .f, type = "wild", warnings = refits$warnings)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.