R/bootstrap_lmer.R

In LiamDBailey/MyFuncs:

#' Bootstrap confidence intervals for (g)lmer models.
#'
#' @param model (g)lmer model for confidence interval estimation.
#' @param n Number of bootstrap iterations
#' @param conf Size of confidence interval. 0.95 by default.
#' @param ... Other arguments passed to bootMer in lme4.
#'
#' @return A dataframe with confidence interval information.
#' @export
#' @import lme4
#' @import boot
#' @importFrom stats sigma
#'
#' @examples
#' #Fit model
#' library(lme4)
#' library(boot)
#' data(sleepstudy)
#' data(cbpp)
#'
#' fm1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
#'
#' #Calculate 95% CIs with parametric boostrap and 100 iterations
#' bootstrap_lmer(fm1, n = 100, conf = 0.95, type = "parametric")
#'
#' #Do the same with glmer
#' gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
#' data = cbpp, family = binomial)
#'
#' #Calculate 95% CIs with parametric boostrap and 100 iterations
#' bootstrap_lmer(gm1, n = 100, conf = 0.95, type = "parametric")
#'

bootstrap_lmer <- function(model, n, conf = 0.95, ...){

  #Build a formula to extract parameter estimates AND random effects estimates.
  Extract <- function(.) {
    c(beta = fixef(.), sigma=sigma(.), sig01=sqrt(unlist(VarCorr(.))))
  }

  #Bootstrap with n iterations and extract parameter estimates and random effects estimates.
  #Bootstrapping options are passed with ...
  bootstrap_test <- lme4::bootMer(model, Extract, nsim = n, ...)

  #Determine confidence intervals
  #Create loop to extract variable estimates and CIs.
  output <- as.data.frame(summary(model)$coefficient)
  output$low <- NA
  output$hi <- NA

  for(i in 1:nrow(output)){

    CI <- boot::boot.ci(bootstrap_test, conf = conf, index = i, type=c("norm"))$normal

    output$low[i] <- CI[2]

    output$hi[i]  <- CI[3]

  }

  return(output)

}