R/boot.R

In jrklasen/relMM: Mixed Models incopareting relationship matrices

#' Parametric bootstrapping
#'
#' Parametric bootstrapping for mixed model with relationship matrix
#'
#' @param model a \code{\link[pedigreemm]{pedigreemm}} model object.
#' @param func a function which extracts the parameter of interest from the 
#' model. This function should have one argument - the model object.
#' @param R the number of bootstrapping runs.
#' @param cores the number of cores for parallel computing at Unix-like 
#' machines.
#'
#' @importFrom boot boot
#' @importFrom stats update simulate
#' @export
mmBoot <- function(model, func, R = 1000, cores = 1) {
  out <- boot(list(data = model@frame, model = model, func = func),
              statistic = function(x) {
                fit <- update(x$model, data = x$data)
                x$func(fit)
              }, 
              ran.gen = function(x, ...) {
                x$data[, 1] <- unlist(simulate(x$model, use.u = TRUE))
                x
              },
              sim = "parametric", R = R,
              parallel = "multicore", ncpus = cores)
  out
}


#' Beta Approximation Confidence Intervals
#'
#' Using the bootstrapping mean and variance to calculate beta distributed two-sided 
#' confidence intervals.
#'
#' @param object a bootstrap object returned from a call to boot.
#' @param conf	a scalar or vector containing the confidence level of the 
#' required interval.
#' @param index the index of the statistic of interest within the output of a 
#' call to boot.out$statistic.
#'
#' @importFrom stats qbeta var
#' @export
beta.ci <- function(object, conf = 0.95, index = 1) {
  stopifnot(class(object) == "boot")
  p <- (1 + conf) / 2
  x <- object$t[, index[1]]
  mu0 <- object$t0[index[1]]
  # var (including mean bias) 
  var0 <- var(x) + (mean(x) - mu0)^2 * length(x) / (length(x) - 1)
  # from mean and var to alpha and beta
  a <- ((1 - mu0) / var0 - 1 / mu0) * mu0^2
  b <- a * (1 / mu0 - 1)
  out <- c(qbeta(1 - p, a, b), mu0, qbeta(p, a, b))
  names(out) <- c(paste("lower", conf, "limit"), "mean", paste("upper", conf, "limit"))
  out
}