R/emmeans_mixed_model.R

In VetResearchLMM: Linear Mixed Models - An Introduction with Applications in Veterinary Research

#' Estimated Marginal Means for a Fitted Mixed Model
#'
#' @description
#' Compute estimated marginal means or pairwise comparisons for a fitted linear
#' mixed model using the optional \pkg{emmeans} package.
#'
#' @param model A fitted model object, typically from \code{lme4::lmer()},
#'   \code{lmerTest::lmer()}, or \code{nlme::lme()}.
#' @param specs Specifications for the marginal means, passed to
#'   \code{emmeans::emmeans()}; for example, \code{~ dose} or
#'   \code{~ dose | Drug}.
#' @param pairwise Logical. If \code{FALSE}, return the estimated marginal
#'   means. If \code{TRUE}, return contrasts computed from those marginal means.
#' @param method Contrast method passed to \code{emmeans::contrast()} when
#'   \code{pairwise = TRUE}. The default is \code{"pairwise"}.
#' @param adjust Multiplicity adjustment passed to \code{emmeans::contrast()}
#'   when \code{pairwise = TRUE}. The default is \code{"tukey"}. Use
#'   \code{NULL} to leave the adjustment unspecified.
#' @param ... Additional arguments passed to \code{emmeans::emmeans()}.
#'
#' @return
#' An \code{emmGrid} object from \pkg{emmeans}. With \code{pairwise = FALSE},
#' this contains estimated marginal means. With \code{pairwise = TRUE}, this
#' contains the requested contrasts.
#'
#' @details
#' Estimated marginal means, also called least-squares means, summarize model
#' predictions for factor levels after accounting for the fitted model
#' structure. They are useful after mixed model fitting because fixed-effect
#' coefficients are often expressed relative to contrast coding, while marginal
#' means and their contrasts are closer to the scientific comparisons shown in
#' the book examples.
#'
#' This helper complements \code{\link{report_mixed_model}()}. Use
#' \code{\link{report_mixed_model}()} for narrative model interpretation and
#' \code{emmeans_mixed_model()} for post hoc inference, estimated marginal
#' means, and pairwise comparisons.
#'
#' The helper keeps \pkg{emmeans} optional. It does not refit the model or
#' change the estimates; it delegates marginal-mean calculations to
#' \code{emmeans::emmeans()} and, when requested, contrasts to
#' \code{emmeans::contrast()}.
#'
#' @examples
#' if (requireNamespace("lme4", quietly = TRUE) &&
#'     requireNamespace("emmeans", quietly = TRUE)) {
#'   data(ex125, package = "VetResearchLMM")
#'   fit <- lme4::lmer(
#'     Pcv ~ dose * Drug + (1 | Region / Drug),
#'     data = ex125,
#'     REML = TRUE
#'   )
#'
#'   emmeans_mixed_model(fit, ~ dose | Drug, lmer.df = "asymptotic")
#'   emmeans_mixed_model(
#'     fit,
#'     ~ dose | Drug,
#'     pairwise = TRUE,
#'     lmer.df = "asymptotic"
#'   )
#' }
#'
#' @references
#' Lenth, R. V. (2024). \pkg{emmeans}: Estimated Marginal Means, aka
#' Least-Squares Means. R package. See \code{utils::citation("emmeans")}.
#'
#' Duchateau, L., Janssen, P., and Rowlands, G. J. (1998).
#' \emph{Linear Mixed Models: An Introduction with Applications in Veterinary
#' Research}. International Livestock Research Institute.
#'
#' @seealso
#' \code{\link{report_mixed_model}()}, \code{emmeans::emmeans()},
#' \code{emmeans::contrast()}.
#'
#' @export
emmeans_mixed_model <- function(model,
                                specs,
                                pairwise = FALSE,
                                method = "pairwise",
                                adjust = "tukey",
                                ...) {
  if (missing(model) || is.null(model)) {
    stop("`model` must be a fitted model object.", call. = FALSE)
  }

  if (missing(specs) || is.null(specs)) {
    stop("`specs` must be supplied.", call. = FALSE)
  }

  if (!is.logical(pairwise) || length(pairwise) != 1L || is.na(pairwise)) {
    stop("`pairwise` must be a single TRUE or FALSE value.", call. = FALSE)
  }

  if (!requireNamespace("emmeans", quietly = TRUE)) {
    stop(
      "The optional package `emmeans` is required. Install it with install.packages(\"emmeans\").",
      call. = FALSE
    )
  }

  emmeans_fun <- getExportedValue("emmeans", "emmeans")
  emm <- emmeans_fun(model, specs = specs, ...)

  if (!pairwise) {
    return(emm)
  }

  contrast_fun <- getExportedValue("emmeans", "contrast")
  contrast_args <- list(object = emm, method = method)
  if (!is.null(adjust)) {
    contrast_args$adjust <- adjust
  }

  do.call(contrast_fun, contrast_args)
}