R/exactRLRT.R

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

#' Restricted Likelihood Ratio Tests for additive and linear mixed models
#' 
#' This function provides an (exact) restricted likelihood ratio test based on
#' simulated values from the finite sample distribution for testing whether the
#' variance of a random effect is 0 in a linear mixed model with known
#' correlation structure of the tested random effect and i.i.d. errors.
#' 
#' Testing in models with only a single variance component require only the
#' first argument \code{m}. For testing in models with multiple variance
#' components, the fitted model \code{m} must contain \bold{only} the random
#' effect set to zero under the null hypothesis, while \code{mA} and \code{m0}
#' are the models under the alternative and the null, respectively. For models
#' with a single variance component, the simulated distribution is exact if the
#' number of parameters (fixed and random) is smaller than the number of
#' observations. Extensive simulation studies (see second reference below)
#' confirm that the application of the test to models with multiple variance
#' components is safe and the simulated distribution is correct as long as the
#' number of parameters (fixed and random) is smaller than the number of
#' observations and the nuisance variance components are not superfluous or
#' very small. We use the finite sample distribution of the restricted
#' likelihood ratio test statistic as derived by Crainiceanu & Ruppert (2004).
#' 
#' No simulation is performed if the observed test statistic is 0. (i.e., if the
#' fit of the model fitted under the alternative is indistinguishable from the
#' model fit under H0), since the  p-value is always 1 in this case.
#' 
#' @param m The fitted model under the alternative or, for testing in models
#' with multiple variance components, the reduced model containing only the
#' random effect to be tested (see Details), an \code{lme}, \code{lmerMod} or
#' \code{spm} object
#' @param mA The full model under the alternative for testing in models with
#' multiple variance components
#' @param m0 The model under the null for testing in models with multiple
#' variance components
#' @param seed input for \code{set.seed}
#' @param nsim Number of values to simulate
#' @param log.grid.hi Lower value of the grid on the log scale. See
#' \code{\link{exactRLRT}}.
#' @param log.grid.lo Lower value of the grid on the log scale. See
#' \code{\link{exactRLRT}}.
#' @param gridlength Length of the grid. See \code{\link{exactLRT}}.
#' @param parallel The type of parallel operation to be used (if any). If
#' missing, the default is "no parallelization").
#' @param ncpus integer: number of processes to be used in parallel operation:
#' typically one would chose this to the number of available CPUs. Defaults to
#' 1, i.e., no parallelization.
#' @param cl An optional parallel or snow cluster for use if parallel = "snow".
#' If not supplied, a cluster on the local machine is created for the duration
#' of the call.
#' @return A list of class \code{htest} containing the following components:
#' @return A list of class \code{htest} containing the following components:
#' \itemize{
#' \item \code{statistic} the observed likelihood ratio
#' \item \code{p} p-value for the observed test statistic
#' \item \code{method} a character string indicating what type of test was
#' performed and how many values were simulated to determine the critical value
#' \item \code{sample} the samples from the null distribution returned by
#' \code{\link{RLRTSim}}
#' }
#' @author Fabian Scheipl, bug fixes by Andrzej Galecki, updates for
#' \pkg{lme4}-compatibility by Ben Bolker
#' @seealso \code{\link{RLRTSim}} for the underlying simulation algorithm;
#' \code{\link{exactLRT}} for likelihood based tests
#' @references Crainiceanu, C. and Ruppert, D. (2004) Likelihood ratio tests in
#' linear mixed models with one variance component, \emph{Journal of the Royal
#' Statistical Society: Series B},\bold{66},165--185.
#' 
#' Greven, S., Crainiceanu, C., Kuechenhoff, H., and Peters, A. (2008)
#' Restricted Likelihood Ratio Testing for Zero Variance Components in Linear
#' Mixed Models, \emph{Journal of Computational and Graphical Statistics},
#' \bold{17} (4): 870--891.
#' 
#' Scheipl, F., Greven, S. and Kuechenhoff, H. (2008) Size and power of tests
#' for a zero random effect variance or polynomial regression in additive and
#' linear mixed models.  \emph{Computational Statistics & Data Analysis},
#' \bold{52}(7):3283--3299.
#' @keywords htest
#' @examples
#' 
#' data(sleepstudy, package = "lme4")
#' mA <- lme4::lmer(Reaction ~ I(Days-4.5) + (1|Subject) + (0 + I(Days-4.5)|Subject), 
#'   data = sleepstudy)
#' m0 <- update(mA, . ~ . - (0 + I(Days-4.5)|Subject))
#' m.slope  <- update(mA, . ~ . - (1|Subject))
#' #test for subject specific slopes:
#' exactRLRT(m.slope, mA, m0)
#' 
#' library(mgcv)
#' data(trees)
#' #test quadratic trend vs. smooth alternative
#' m.q<-gamm(I(log(Volume)) ~ Height + s(Girth, m = 3), data = trees, 
#'   method = "REML")$lme
#' exactRLRT(m.q)
#' #test linear trend vs. smooth alternative
#' m.l<-gamm(I(log(Volume)) ~ Height + s(Girth, m = 2), data = trees, 
#'   method = "REML")$lme
#' exactRLRT(m.l)
#' 
#' @export exactRLRT
#' @importFrom stats anova cov2cor logLik quantile 
#' @importFrom utils packageVersion
'exactRLRT' <- function(m, mA = NULL, m0 = NULL, seed = NA, 
  nsim = 10000, log.grid.hi = 8, log.grid.lo = -10, gridlength = 200,
  parallel = c("no", "multicore", "snow"), 
  ncpus = 1L, cl = NULL) {
  if (inherits(m, "spm")) {
    m <- m$fit
    class(m) <- "lme"
  }
  if (any(class(m) %in% c("amer", "mer")))
    stop("Models fit with package <amer> or versions of <lme4> below 1.0 are no longer supported.")
  c.m <- class(m)
  if (!any(c.m %in% c("lme", "lmerMod", "merModLmerTest", "lmerModLmerTest"))) 
    stop("Invalid <m> specified. \n")
  if (any(c.m %in% c("merModLmerTest", "lmerModLmerTest"))) 
    c.m <- "lmerMod"
  if ("REML" != switch(c.m, 
    lme = m$method, 
    lmerMod = ifelse(lme4::isREML(m), "REML", "ML"))){
    message("Using restricted likelihood evaluated at ML estimators.")
    message("Refit with method=\"REML\" for exact results.")	
  }
  
  d <- switch(c.m, lme = extract.lmeDesign(m), 
    lmerMod = extract.lmerModDesign(m))
  X <- d$X
  qrX <- qr(X)
  Z <- d$Z
  y <- d$y
  Vr <- d$Vr
  if (all(Vr == 0)) {
    # this only happens if the estimate of the tested variance component is 0. 
    # since we still want chol(cov2cor(Vr)) to work, this does the trick.
    diag(Vr) <- 1
  }
  K <- ncol(Z)
  n <- nrow(X)
  p <- ncol(X)
  if (is.null(mA) && is.null(m0)) {
    if (length(d$lambda) != 1 || d$k != 1) 
      stop("multiple random effects in model - 
                 exactRLRT needs <m> with only a single random effect.")
    #2*restricted ProfileLogLik under H0: lambda=0
    res <- qr.resid(qrX, y)
    R <- qr.R(qrX)
    detXtX <- det(t(R) %*% R)
    reml.H0 <- -((n - p) * log(2 * pi) + (n - p) * log(sum(res^2)) + 
        log(detXtX) + (n - p) - (n - p) * log(n - p))
    #observed value of the test-statistic
    reml.obs <- 2 * logLik(m, REML = TRUE)[1]
    rlrt.obs <- max(0, reml.obs - reml.H0)
    lambda <- d$lambda
  } else {
    nonidentfixmsg <- 
      "Fixed effects structures of <mA> and <m0> not identical.
        REML-based inference not appropriate."
    if (c.m == "lme") {
      if (any(mA$fixDF$terms != m0$fixDF$terms)) 
        stop(nonidentfixmsg)
    } else {
      if (c.m == "mer") {
        if (any(mA@X != m0@X)) 
          stop(nonidentfixmsg)
      } else {
        if (c.m == "lmerMod") {
          if (any(lme4::getME(mA,"X") != lme4::getME(m0,"X")))
            stop(nonidentfixmsg)
        }
      }     
    }
    lmer_nm <- if (utils::packageVersion("lme4")<="1.1.21") "Df" else "npar"
    ## bug fix submitted by Andrzej Galecki 3/10/2009
    DFx <- switch(c.m, lme = anova(mA,m0)$df, 
                  lmerMod = anova(mA, m0, refit = FALSE)[[lmer_nm]]) 
    if (abs(diff(DFx)) > 1) {
      stop("Random effects not independent - covariance(s) set to 0 under H0.\n
                 exactRLRT can only test a single variance.\n")
    }
    rlrt.obs <- max(0, 2 * (logLik(mA, REML = TRUE)[1] - 
        logLik(m0, REML = TRUE)[1]))
  }
  p <- if (rlrt.obs != 0) {
    sample <- RLRTSim(X, Z, qrX = qrX, sqrt.Sigma = chol(cov2cor(Vr)), 
      lambda0 = 0, seed = seed, nsim = nsim, 
      log.grid.hi = log.grid.hi, 
      log.grid.lo = log.grid.lo, gridlength = gridlength, 
      parallel = match.arg(parallel), 
      ncpus = ncpus, cl = cl)
    if (quantile(sample, 0.9) == 0) {
      warning("Null distribution has mass ", mean(sample == 
          0), " at zero.\n")
    }
    mean(rlrt.obs < sample)
  } else {
    message("Observed RLRT statistic is 0, no simulation performed.")
    nsim <- 0
    sample <- NULL
    1
  }  
  RVAL <- list(statistic = c(RLRT = rlrt.obs), p.value = p, 
    method = paste("simulated finite sample distribution of RLRT.\n
                                (p-value based on", 
      nsim, "simulated values)"), sample = sample)
  class(RVAL) <- "htest"
  return(RVAL)
}