R/CIfunction_paramMLvarCompRob.R

In confintROB: Confidence Intervals for Robust and Classical Linear Mixed Model Estimators

#Title
# @title Parametric Confidence Interval for varComprob objects
#
# @param model an object of class varComprob
# @param nsim number of bootstrap samples, positive integer
# @param max.tries number of times to try to produce a valid model fit
#   before giving up
# @param .export passed on to \code{\link[foreach]{foreach}}
#
# @return Returns a parametric Confidence Interval
#' @importFrom stats rnorm
#' @importFrom methods is
# @export
parametric.varComprob <-
  function(model,
           nsim,
           max.tries = 100,
           .export,
           data,
           varComprob.random,
           ...) {
    sample <- createParamSampleFunction(model, data)
    fit <-
      createFitFunction.varComprob(model, data, varComprob.random)
    bootstrap(model, nsim, max.tries, .export, sample, fit)
  }

#Title
# @title Generation of samples by parametric bootstrap for varComprob objects
#
# @param model an object of class varComprob
# @param data the original dataset
#
# @return Returns a sample generated by parametric bootstrap
#' @importFrom stats rnorm model.response
#' @importFrom methods is
#' @importFrom mvtnorm rmvnorm
createParamSampleFunction <- function(model, data) {
  if (is.null(model[["model"]])) {
    stop("model$model is null, please run varComprob with argument 'model=TRUE'")
  }
  
  X <- model.matrix(model, data)
  ntot <- NROW(X)
  nsubj <- dim(model$X)[1]
  nobs <- dim(model$X)[2]
  pred_fixed <- drop(X %*% model$beta)
  
  # build covariance matrix for observations of a single subject
  sigK <- 0
  for (i in seq_along(model$K)) {
    pred_i <- model$eta[i] * as.matrix(model$K[i])[[1]]
    sigK <- sigK + pred_i
  }
  I_err <- model$eta0 * diag(nsubj)
  SIGMA <- sigK + I_err
  
  # compute reverse permutation used to bring
  # errors generated ordered by observations
  # into the order as they are in data
  groups <- model[["model"]][["(groups)"]]
  revperm <- order(order(groups[, 1], groups[, 2]))
  
  sample <- function(n) {
    E_ordered <-
      rmvnorm(
        n = n * nobs,
        mean = rep(0, nsubj),
        sigma = SIGMA,
        method = "eigen"
      )
    e_ordered_by_sample <-
      split(E_ordered, rep(1:n, times = nsubj, each = nobs))
    yboot_by_sample <-
      lapply(e_ordered_by_sample, function(e)
        pred_fixed + e[revperm])
    return(yboot_by_sample)
  }
  return(sample)
}