R/model_comparison.R

In dyadMLM: Tools for Dyadic Multilevel Models

#' Compare nested glmmTMB models fitted to equivalent data
#'
#' Performs a likelihood-ratio test for two nested `glmmTMB` models. The models
#' may use ordinary data frames or [dyadMLM_data][prepare_dyad_data()]
#' objects, and their calls do not need to refer to the same R object. Models may
#' be supplied in either order. The model with fewer estimated parameters is
#' shown first in the result.
#'
#' @param model1,model2 Two fitted `glmmTMB` models to compare.
#' @param alpha Significance level used for the printed conclusion. Must be one
#'   number between 0 and 1.
#'
#' @details
#' Both model calls must use named data-frame objects that remain available when
#' the models are compared. The checks assume these objects have not been
#' modified since fitting. All ordinary data columns must be identical,
#' including their types and attributes. For `dyadMLM_data`, package-generated
#' columns may differ, but the original columns must be identical. Ordinary and
#' prepared data may be compared with each other. Dyad metadata are checked when
#' both models use `dyadMLM_data`. The function also checks fitted rows,
#' outcomes, weights and offsets, model family and link, maximum-likelihood
#' estimation, and model convergence. Each model must use the same untransformed
#' response column.
#'
#' Because the test requires maximum-likelihood fits, models fitted with
#' `priors` in `glmmTMB::glmmTMB()` are rejected. Their estimates are maximum a
#' posteriori, not maximum likelihood.
#'
#' These checks establish that the models use equivalent observations. They
#' cannot establish that one model is mathematically nested within the other.
#' The caller remains responsible for supplying genuinely nested models. The
#' usual chi-squared reference distribution may also be inappropriate when
#' tested variance parameters are on the boundary.
#' `alpha` affects only the printed conclusion, not the test results.
#'
#' @return An `anova`-style data frame containing model degrees of freedom,
#'   information criteria, log-likelihoods, the likelihood-ratio statistic, and
#'   its chi-squared p-value. When printed, a short conclusion interprets the
#'   test at the selected significance level.
#'
#' @examples
#' if (requireNamespace("glmmTMB", quietly = TRUE)) {
#'   restricted_data <- prepare_dyad_data(
#'     dyads_cross,
#'     dyad = coupleID,
#'     member = personID,
#'     role = gender,
#'     # All three observed compositions in `dyads_cross` are detected and retained
#'     # by default. This example focuses on `female-male` dyads, so we restrict the
#'     # analysis here.
#'     keep_compositions = "female-male",
#'     include_arbitrary_member_contrast = TRUE,
#'     seed = 123
#'   )
#'   full_data <- restricted_data
#'
#'   # Compare an exchangeable model with a distinguishable model.
#'   restricted_model <- glmmTMB::glmmTMB(
#'     closeness ~ 1 +
#'       us(1 | coupleID) +
#'       us(0 + .member_contrast_arbitrary | coupleID),
#'     dispformula = ~ 0,
#'     data = restricted_data
#'   )
#'
#'   full_model <- glmmTMB::glmmTMB(
#'     closeness ~ 0 + gender +
#'       us(0 + gender | coupleID),
#'     dispformula = ~ 0,
#'     data = full_data
#'   )
#'
#'   # Test equal means and residual variances jointly.
#'   compare_nested_models(restricted_model, full_model)
#' }
#'
#' @export
compare_nested_models <- function(model1, model2, alpha = 0.05) {
  if (!is.numeric(alpha) || length(alpha) != 1L ||
      !is.finite(alpha) || alpha <= 0 || alpha >= 1) {
    stop("`alpha` must be one finite number between 0 and 1.", call. = FALSE)
  }

  validate_comparison_model(model1, "model1")
  validate_comparison_model(model2, "model2")

  # Save the calling environment so we can find the original data objects.
  caller_env <- parent.frame()
  model1_data <- comparison_model_data(model1, caller_env, "model1")
  model2_data <- comparison_model_data(model2, caller_env, "model2")
  validate_comparison_data(
    model1,
    model2,
    model1_data,
    model2_data
  )

  # Use the argument expressions as row labels in the result.
  labels <- c(
    deparse1(substitute(model1)),
    deparse1(substitute(model2))
  )
  return(likelihood_ratio_comparison(model1, model2, labels, alpha))
}

# Model and data checks -----------------------------------------------------

validate_comparison_model <- function(model, name) {
  if (!inherits(model, "glmmTMB")) {
    stop(
      sprintf("`%s` must be a fitted `glmmTMB` model.", name),
      call. = FALSE
    )
  }
  if (isTRUE(model$modelInfo$REML)) {
    stop(
      sprintf(
        "The `%s` model was fitted with REML. Refit both models with `REML = FALSE` before comparing them.",
        name
      ),
      call. = FALSE
    )
  }
  # An empty prior specification is stored as a data frame without rows.
  if (NROW(model$modelInfo$priors) > 0L) {
    stop(
      sprintf(
        "The `%s` model was fitted with priors. Refit both models without `priors` before comparing them.",
        name
      ),
      call. = FALSE
    )
  }
  if (!is.null(model$fit$convergence) &&
      !isTRUE(model$fit$convergence == 0)) {
    stop(
      sprintf(
        "The `%s` model did not converge. Diagnose and refit it before comparing models.",
        name
      ),
      call. = FALSE
    )
  }
  # pdHess is NULL when the Hessian was not calculated, e.g. with se = FALSE.
  if (!is.null(model$sdr$pdHess) && !isTRUE(model$sdr$pdHess)) {
    stop(
      sprintf(
        paste0(
          "The `%s` model has a non-positive-definite Hessian matrix. ",
          "Diagnose and refit it before comparing models."
        ),
        name
      ),
      call. = FALSE
    )
  }

  invisible(model)
}

comparison_model_data <- function(model, caller_env, argument) {
  data_call <- model$call$data
  if (!is.symbol(data_call)) {
    stop(
      sprintf(
        "The `%s` model must have been fitted with a named data frame in `data`.",
        argument
      ),
      call. = FALSE
    )
  }

  data_name <- as.character(data_call)
  model_formula <- stats::formula(model, component = "cond")

  # A model fitted inside a helper may keep its data with its formula.
  # Search there first, then where compare_nested_models() was called.
  environments <- list(environment(model_formula), caller_env)

  for (search_env in environments) {
    if (!is.environment(search_env)) {
      next
    }
    if (exists(data_name, envir = search_env, inherits = TRUE)) {
      candidate <- get(data_name, envir = search_env, inherits = TRUE)
      if (is.data.frame(candidate)) {
        return(candidate)
      }
    }
  }

  stop(
    sprintf(
      paste0(
        "Could not recover the named data frame used by `%s`. ",
        "Keep it available when comparing models."
      ),
      argument
    ),
    call. = FALSE
  )
}

validate_comparison_data <- function(model1, model2,
                                     model1_data, model2_data) {
  if (nrow(model1_data) != nrow(model2_data)) {
    stop(
      "The two data objects have different numbers of rows (",
      nrow(model1_data), " and ", nrow(model2_data),
      "). Fit both models from the same analysis data.",
      call. = FALSE
    )
  }

  both_prepared <- inherits(model1_data, "dyadMLM_data") &&
    inherits(model2_data, "dyadMLM_data")

  if (both_prepared) {
    # Both preparations must agree on the columns that define the dyadic data.
    model1_meta <- attr(model1_data, "dyadMLM")
    model2_meta <- attr(model2_data, "dyadMLM")
    for (field in c("dyad", "member", "role", "time")) {
      if (!identical(model1_meta[[field]], model2_meta[[field]])) {
        stop(
          sprintf(
            "The two prepared datasets use different `%s` variables.",
            field
          ),
          call. = FALSE
        )
      }
    }
  }

  # Generated dyadMLM columns may differ across model parameterizations. Use
  # their recorded names rather than treating a prefix as column ownership.
  model1_columns <- names(model1_data)
  model2_columns <- names(model2_data)
  if (inherits(model1_data, "dyadMLM_data")) {
    generated_columns <- dyad_generated_columns(
      attr(model1_data, "dyadMLM")
    )$column
    model1_columns <- model1_columns[!model1_columns %in% generated_columns]
  }
  if (inherits(model2_data, "dyadMLM_data")) {
    generated_columns <- dyad_generated_columns(
      attr(model2_data, "dyadMLM")
    )$column
    model2_columns <- model2_columns[!model2_columns %in% generated_columns]
  }

  if (!setequal(model1_columns, model2_columns)) {
    stop(
      "The two data objects do not contain the same original columns. ",
      "Prepare both models from data with the same original columns.",
      call. = FALSE
    )
  }

  for (column in model1_columns) {
    if (!identical(model1_data[[column]], model2_data[[column]])) {
      stop(
        sprintf(
          paste0(
            "Column `%s` differs between the two data objects. ",
            "Fit both models from the same unchanged analysis data."
          ),
          column
        ),
        call. = FALSE
      )
    }
  }

  # Inspect the left side directly so transformations such as log(y) fail.
  model1_formula <- stats::formula(model1, component = "cond")
  model2_formula <- stats::formula(model2, component = "cond")
  model1_response <- model1_formula[[2L]]
  model2_response <- model2_formula[[2L]]

  if (!is.symbol(model1_response) || !is.symbol(model2_response)) {
    stop(
      "Only an untransformed outcome stored in the data is supported. ",
      "Create the transformed outcome as a new data column and refit both models with that column on the left-hand side.",
      call. = FALSE
    )
  }

  model1_response <- as.character(model1_response)
  model2_response <- as.character(model2_response)
  if (!identical(model1_response, model2_response)) {
    stop(
      "The two models use different outcome variables (`",
      model1_response, "` and `", model2_response,
      "`). Use the same outcome in both models.",
      call. = FALSE
    )
  }
  if (!model1_response %in% model1_columns) {
    stop(
      "Only an untransformed outcome stored in the data is supported. ",
      "Create the transformed outcome as a new data column and refit both models with that column on the left-hand side.",
      call. = FALSE
    )
  }

  # Row names show which observations glmmTMB kept after handling missing data.
  model1_rows <- row.names(model1$frame)
  model2_rows <- row.names(model2$frame)
  if (!identical(model1_rows, model2_rows)) {
    stop(
      "The two models were fitted to different observation rows. ",
      "Check missing-data patterns and refit both models on the same rows.",
      call. = FALSE
    )
  }

  # Compare what the models actually fitted, not only the current data objects.
  model1_outcome <- unname(stats::model.response(model1$frame))
  model2_outcome <- unname(stats::model.response(model2$frame))
  if (!identical(model1_outcome, model2_outcome)) {
    stop(
      "The fitted outcome values differ between the two models. ",
      "Refit both models from the same unchanged analysis data.",
      call. = FALSE
    )
  }

  # No weights means all observations have weight one.
  model1_weights <- stats::model.weights(model1$frame)
  model2_weights <- stats::model.weights(model2$frame)
  if (is.null(model1_weights)) model1_weights <- rep(1, length(model1_rows))
  if (is.null(model2_weights)) model2_weights <- rep(1, length(model2_rows))
  if (!identical(as.numeric(model1_weights), as.numeric(model2_weights))) {
    stop(
      "The two models use different observation weights. ",
      "Use identical weights when comparing their likelihoods.",
      call. = FALSE
    )
  }

  # glmmTMB stores separate offsets for its three model components.
  offset_names <- c("offset", "zioffset", "dispoffset")
  model1_offsets <- model1$obj$env$data[offset_names]
  model2_offsets <- model2$obj$env$data[offset_names]
  if (!identical(model1_offsets, model2_offsets)) {
    stop(
      "The two models use different offsets. ",
      "Use identical offsets when comparing their likelihoods.",
      call. = FALSE
    )
  }

  model1_family <- model1$modelInfo$family
  model2_family <- model2$modelInfo$family
  same_family <- identical(model1_family$family, model2_family$family)
  same_link <- identical(model1_family$link, model2_family$link)
  if (!same_family || !same_link) {
    stop(
      "The two models must use the same family and link. Found `",
      model1_family$family, "` with `", model1_family$link,
      "` and `", model2_family$family, "` with `", model2_family$link,
      "`.",
      call. = FALSE
    )
  }

  invisible(TRUE)
}

# Likelihood-ratio test and output -----------------------------------------

likelihood_ratio_comparison <- function(model1, model2, labels, alpha) {
  model1_log_lik <- stats::logLik(model1)
  model2_log_lik <- stats::logLik(model2)
  model1_df <- attr(model1_log_lik, "df")
  model2_df <- attr(model2_log_lik, "df")

  valid_df <- length(model1_df) == 1L && length(model2_df) == 1L &&
    is.finite(model1_df) && is.finite(model2_df)
  if (!valid_df) {
    stop(
      "Could not determine the models' numbers of parameters.",
      call. = FALSE
    )
  }
  if (model1_df == model2_df) {
    stop(
      "The two models must have different numbers of estimated parameters.",
      call. = FALSE
    )
  }

  if (model1_df < model2_df) {
    restricted <- model1
    full <- model2
    restricted_log_lik <- model1_log_lik
    full_log_lik <- model2_log_lik
    restricted_df <- model1_df
    full_df <- model2_df
  } else {
    restricted <- model2
    full <- model1
    restricted_log_lik <- model2_log_lik
    full_log_lik <- model1_log_lik
    restricted_df <- model2_df
    full_df <- model1_df
    labels <- rev(labels)
  }

  restricted_log_lik <- as.numeric(restricted_log_lik)
  full_log_lik <- as.numeric(full_log_lik)
  statistic <- 2 * (full_log_lik - restricted_log_lik)

  # Allow only a negligible negative value caused by numerical rounding.
  if (!is.finite(statistic) || statistic < -sqrt(.Machine$double.eps)) {
    stop(
      paste0(
        "The larger model has a lower log-likelihood than the smaller model. ",
        "Check convergence and nesting."
      ),
      call. = FALSE
    )
  }

  statistic <- max(0, statistic)
  df_difference <- full_df - restricted_df
  p_value <- stats::pchisq(statistic, df_difference, lower.tail = FALSE)

  formatted_p <- format.pval(p_value, digits = 3, eps = 0.001)
  if (startsWith(formatted_p, "<")) {
    p_text <- paste("p <", substring(formatted_p, 2L))
  } else {
    p_text <- paste("p =", formatted_p)
  }

  level_text <- paste0(format(100 * alpha, trim = TRUE), "%")
  if (p_value < alpha) {
    conclusion <- sprintf(
      paste0(
        "Conclusion (%s level): The likelihood-ratio test provides evidence ",
        "that `%s` fits better than `%s` (%s)."
      ),
      level_text, labels[2], labels[1], p_text
    )
  } else {
    conclusion <- sprintf(
      paste0(
        "Conclusion (%s level): The likelihood-ratio test finds no clear ",
        "improvement from `%s` to `%s` (%s). This does not establish equal fit."
      ),
      level_text, labels[1], labels[2], p_text
    )
  }

  # Build the table here because anova.glmmTMB() rejects different data names.
  out <- data.frame(
    Df = c(restricted_df, full_df),
    AIC = c(stats::AIC(restricted), stats::AIC(full)),
    BIC = c(stats::BIC(restricted), stats::BIC(full)),
    logLik = c(restricted_log_lik, full_log_lik),
    deviance = -2 * c(restricted_log_lik, full_log_lik),
    Chisq = c(NA_real_, statistic),
    `Chi Df` = c(NA_real_, df_difference),
    `Pr(>Chisq)` = c(NA_real_, p_value),
    check.names = FALSE,
    row.names = labels
  )

  attr(out, "heading") <- c(
    "Likelihood-ratio test for nested models fitted to equivalent data",
    "Assumes mathematical nesting and an appropriate chi-squared reference distribution.",
    ""
  )
  attr(out, "conclusion") <- conclusion
  attr(out, "alpha") <- alpha
  class(out) <- c("dyadMLM_model_comparison", "anova", "data.frame")
  out
}

#' @export
#' @noRd
print.dyadMLM_model_comparison <- function(x, ...) {
  # Print the usual ANOVA table first, then the plain-language conclusion.
  NextMethod()
  cat("\n", attr(x, "conclusion"), "\n", sep = "")
  invisible(x)
}