R/extract.lmeDesign.R

In RLRsim: Exact (Restricted) Likelihood Ratio Tests for Mixed and Additive Models

#' Extract the Design of a linear mixed model
#'
#' These functions extract various elements of the design of a fitted
#' \code{lme}-, \code{mer} or \code{lmerMod}-Object.  They are called by
#' \code{exactRLRT} and \code{exactLRT}.
#'
#'
#' @aliases extract.lmerModDesign extract.lmeDesign
#' @param m a fitted \code{lme}- or \code{merMod}-Object
#' @return a a list with components
#' \itemize{
#' \item \code{Vr} estimated covariance of the random effects divided by the
#' estimated variance of the residuals
#' \item \code{X} design of the fixed effects
#' \item \code{Z} design of the random effects
#' \item \code{sigmasq} variance of the residuals
#' \item \code{lambda} ratios of the variances of the random effects and the
#' variance of the residuals
#' \item \code{y} response variable
#' }
#' @author Fabian Scheipl, \code{extract.lmerModDesign} by Ben Bolker.
#' Many thanks to Andrzej Galecki and Tomasz Burzykowski for bug fixes.
#' @keywords utilities
#' @examples
#'
#' library(nlme)
#' design <- extract.lmeDesign(lme(distance ~ age + Sex, data = Orthodont,
#'                              random = ~ 1))
#' str(design)
#'
#' @export extract.lmeDesign
#' @importFrom stats complete.cases formula model.frame model.matrix
#' @importFrom nlme getGroups
#' @importFrom mgcv tensor.prod.model.matrix
extract.lmeDesign <- function(m) {
  start.level = 1

  data <- if (any(!complete.cases(m$data))) {
    warning("Removing incomplete cases from supplied data.")
    m$data[complete.cases(m$data), ]
  } else m$data
  grps <- getGroups(m)
  n <- length(grps)
  X <- list()
  grp.dims <- m$dims$ncol
  Zt <- model.matrix(m$modelStruct$reStruct, data)
  cov <- as.matrix(m$modelStruct$reStruct)
  i.col <- 1
  n.levels <- length(m$groups)
  Z <- matrix(0, n, 0)
  if (start.level <= n.levels) {
    for (i in 1:(n.levels - start.level + 1)) {
      if (length(levels(m$groups[[n.levels - i + 1]])) != 1) {
        grps <- m$groups[[n.levels - i + 1]]
        X[[1]] <- model.matrix(
          ~ grps - 1,
          contrasts.arg = list(grps = "contr.treatment")
        )
      } else {
        X[[1]] <- matrix(1, n, 1)
      }
      X[[2]] <- as.matrix(Zt[, i.col:(i.col + grp.dims[i] - 1)])
      i.col <- i.col + grp.dims[i]
      Z <- cbind(tensor.prod.model.matrix(X), Z)
    }
    Vr <- matrix(0, ncol(Z), ncol(Z))
    start <- 1
    for (i in 1:(n.levels - start.level + 1)) {
      k <- n.levels - i + 1
      for (j in 1:m$dims$ngrps[i]) {
        stop <- start + ncol(cov[[k]]) - 1
        Vr[ncol(Z) + 1 - (stop:start), ncol(Z) + 1 - (stop:start)] <- cov[[k]]
        start <- stop + 1
      }
    }
  }
  X <- if (inherits(m$call$fixed, "name") && !is.null(m$data$X)) {
    m$data$X
  } else {
    model.matrix(formula(eval(m$call$fixed)), data)
  }
  y <- as.vector(
    matrix(m$residuals, ncol = NCOL(m$residuals))[, NCOL(m$residuals)] +
      matrix(m$fitted, ncol = NCOL(m$fitted))[, NCOL(m$fitted)]
  )
  return(list(
    Vr = Vr, #Cov(RanEf)/Var(Error)
    X = X,
    Z = Z,
    sigmasq = m$sigma^2,
    lambda = unique(diag(Vr)),
    y = y,
    k = n.levels
  ))
}