R/compare.R

In INLAvaan: Approximate Bayesian Latent Variable Analysis

#' Compare Bayesian Models Fitted with INLAvaan
#'
#' Compare two or more Bayesian SEM fitted with INLAvaan, reporting model-fit
#' statistics and (optionally) fit indices side by side.
#'
#' @details
#' All models appear in the comparison table. When incremental fit indices
#' (BCFI, BTLI, BNFI) are requested via `fit.measures`, they are scaled
#' against the independence (null) model, fitted once on the data of the
#' first model and shared by every model in the table (see
#' [bfit_indices()]).
#'
#' The default table always includes:
#'
#'   - **npar**: Number of free parameters.
#'   - **Marg.Loglik**: Approximated marginal log-likelihood.
#'   - **logBF**: Natural-log Bayes factor relative to the best model.
#'   - **DIC** / **pD**: Deviance Information Criterion and effective number
#'     of parameters (when the fit computed the DIC, i.e. `test` included
#'     `"dic"` during fitting; the default `"standard"` does).
#'
#' Set `fit.measures` to a character vector of measure names (anything
#' returned by [fitMeasures()][lavaan::fitMeasures]) to append extra columns.
#' Use `fit.measures = "all"` to include every available measure.
#'
#' Set `loo = TRUE` to compare models by leave-one-out cross-validation
#' (see [loo()]). This appends **ELPD** / **SE** (the Taylor expected log
#' predictive density and its standard error), **p_loo**, and,
#' against the best-ELPD model, the difference **elpd_diff** with its
#' *paired* standard error **se_diff** computed from the pointwise
#' contributions (the appropriate uncertainty for nested or same-data
#' comparisons). Every model is scored at one common Taylor order, the
#' lowest any of them can supply: if some unit of some model has no
#' second-order term, all models are compared at first order, since
#' otherwise a change of estimator between models would read as a
#' difference between the models themselves. The order used is stated when
#' the table is printed. The table is then sorted by descending ELPD. All models
#' must be fitted to the same data with matching units; units are paired
#' by id rather than by row order, so fits that stack groups differently
#' -- a pooled fit against a multigroup fit, or multigroup fits with
#' different group orderings -- still pair up unit by unit. For
#' missing-data (FIML) fits, "the same data" also means the same observed
#' entries: each unit is scored on the entries it has, so comparisons
#' require identical missingness patterns across models. All models
#' must also share the
#' score flavour (see [loo()]): mixing fits with modelled covariates
#' (`fixed.x = FALSE`, joint scores) and fixed covariates
#' (`fixed.x = TRUE`, conditional scores) is refused. Joint scores
#' additionally require identical variable sets across models, while
#' conditional scores require only matching outcome variables -- covariate
#' sets may differ, which is the covariate-selection setting. Stored LOO
#' results (`test` including `"loo"` or `"full"`, or [add_loo()]) are
#' reused.
#'
#' `anova()` is disabled for `INLAvaan` fits -- there is no direct Bayesian
#' analogue of the classical likelihood-ratio test -- and points here instead.
#'
#' @param x,y,... Two or more [INLAvaan] (or `inlavaan_internal`) objects
#'   fitted to the same data.
#' @param fit.measures Character vector of additional fit-measure names to
#'   include (e.g. `"BRMSEA"`, `"BCFI"`). Use `"all"` to include every
#'   measure returned by [fitMeasures()][lavaan::fitMeasures]. The default
#'   (`NULL`) shows only the core comparison statistics.
#' @param loo Logical; if `TRUE`, compare models by leave-one-out
#'   cross-validation with paired standard errors (see Details). Defaults to
#'   `FALSE`.
#' @param object An [INLAvaan] object (the `anova()` method, which is disabled
#'   and redirects to `compare()`).
#'
#' @return A data frame of class `compare.inlavaan_internal` containing model
#'   fit statistics, sorted by descending marginal log-likelihood (or by
#'   descending ELPD when `loo = TRUE`).
#'
#' @references <https://lavaan.ugent.be/tutorial/groups.html>
#'
#' @seealso [fitmeasures()], [bfit_indices()]
#'
#' @example inst/examples/ex-model_comparison.R
#' @export
setGeneric("compare", function(x, y, ..., fit.measures = NULL, loo = FALSE) {
  standardGeneric("compare")
})

#' @name compare
#' @rdname compare
#' @aliases compare,INLAvaan-method
#' @export
setMethod(
  "compare",
  "INLAvaan",
  function(x, y, ..., fit.measures = NULL, loo = FALSE) {
    mc <- match.call()
    dots <- list(...)

    # x = baseline, y + unnamed ... = models to compare
    model_objs <- c(list(x, y), dots)
    model_exprs <- c(
      list(mc$x, mc$y),
      as.list(mc)[-1][
        !names(as.list(mc)[-1]) %in%
          c("x", "y", "fit.measures", "loo")
      ]
    )

    modnames <- vapply(model_exprs, deparse, character(1))

    compare_impl(
      models = model_objs,
      modnames = modnames,
      fit.measures = fit.measures,
      loo = loo
    )
  }
)

#' @exportS3Method compare inlavaan_internal
compare.inlavaan_internal <- function(
  x,
  y,
  ...,
  fit.measures = NULL,
  loo = FALSE
) {
  mc <- match.call()
  dots <- list(...)

  model_objs <- c(list(x, y), dots)
  model_exprs <- c(
    list(mc$x, mc$y),
    as.list(mc)[-1][
      !names(as.list(mc)[-1]) %in%
        c("x", "y", "fit.measures", "loo")
    ]
  )

  modnames <- vapply(model_exprs, deparse, character(1))

  compare_impl(
    models = model_objs,
    modnames = modnames,
    fit.measures = fit.measures,
    loo = loo
  )
}

# ---- Internal workhorse ------------------------------------------------------

compare_impl <- function(
  models,
  modnames,
  fit.measures = NULL,
  loo = FALSE
) {
  # Normalise to internal objects, keeping originals for fitMeasures()
  originals <- models
  internals <- lapply(models, function(m) {
    if (inherits(m, "INLAvaan")) {
      m@external$inlavaan_internal
    } else if (inherits(m, "inlavaan_internal")) {
      m # nocov
    } else {
      cli_abort(
        # nocov
        "Each model must be an {.cls INLAvaan} or {.cls inlavaan_internal} object."
      )
    }
  })

  nmod <- length(internals)
  npar <- vapply(internals, function(m) length(m$theta_star), integer(1))
  marg_ll <- vapply(internals, function(m) m$mloglik, numeric(1))
  DIC_vec <- vapply(internals, function(m) m$DIC$dic %||% NA_real_, numeric(1))
  pD_vec <- vapply(internals, function(m) m$DIC$pD %||% NA_real_, numeric(1))

  # Marginal likelihoods, Bayes factors, and DIC are only comparable
  # between fits with the same mean treatment: without a mean structure
  # the saturated means carry an improper flat prior whose arbitrary
  # normalisation cancels within that flavour but is orphaned against a
  # fit with modelled (proper-prior) means. LOO comparisons are unaffected
  # (leave-one-out conditionals are proper under both treatments).
  ms_vec <- vapply(
    internals,
    function(m) isTRUE(m$lavmodel@meanstructure),
    logical(1)
  )
  if (length(unique(ms_vec)) > 1L) {
    hint <- if (isTRUE(loo)) {
      c(
        "i" = "Interpret only the ELPD columns; leave-one-out conditionals
         are proper under both treatments."
      )
    } else {
      c(
        "i" = "Use {.code compare(..., loo = TRUE)} for a comparison that
         is valid across mean treatments."
      )
    }
    cli_warn(c(
      "Comparing fits with and without a mean structure: marginal
       log-likelihoods, Bayes factors, and DIC are not comparable across
       the two mean treatments (the flat-prior normalisation of the
       saturated means does not cancel).",
      hint
    ))
  }

  best_ll <- max(marg_ll)
  logBF <- marg_ll - best_ll

  out <- data.frame(
    Model = modnames,
    npar = npar,
    Marg.Loglik = marg_ll,
    logBF = round(logBF, 3),
    stringsAsFactors = FALSE
  )

  # Append DIC/pD if any model has them
  if (!all(is.na(DIC_vec))) {
    out$DIC <- round(DIC_vec, 3)
    out$pD <- round(pD_vec, 3)
  }

  # Append extra fit measures if requested
  if (!is.null(fit.measures)) {
    has_inlavaan <- vapply(originals, function(m) is(m, "INLAvaan"), logical(1))
    if (!all(has_inlavaan)) {
      # nocov start
      cli_warn(
        "Fit measures require {.cls INLAvaan} objects; skipping for {.cls inlavaan_internal} models."
      )
    } else {
      # nocov end
      # One independence baseline, fitted on the first model's data, serves
      # every model's incremental indices. FALSE skips the refit when none
      # is requested.
      need_incr <- identical(fit.measures, "all") ||
        any(c("BCFI", "BTLI", "BNFI") %in% fit.measures)
      baseline_obj <- if (need_incr) {
        tryCatch(
          fit_independence_baseline(originals[[1]]),
          error = function(e) NULL
        )
      } else {
        FALSE
      }

      fm_list <- lapply(originals, function(m) {
        tryCatch(
          fitMeasures(
            m,
            fit.measures = fit.measures,
            baseline.model = baseline_obj
          ),
          error = function(e) NULL
        )
      })
      # Union of all measure names
      all_names <- unique(unlist(lapply(fm_list, names)))
      for (nm in all_names) {
        out[[nm]] <- vapply(
          fm_list,
          function(fm) {
            if (is.null(fm) || is.na(fm[nm])) NA_real_ else round(fm[nm], 4)
          },
          numeric(1)
        )
      }
    }
  }

  # LOO comparison with paired standard errors
  if (isTRUE(loo)) {
    loo_list <- lapply(internals, function(m) {
      if (!is.null(m$loo)) m$loo else loo.inlavaan_internal(m)
    })

    # Paired differences are only meaningful for matching units on the
    # same data. Units are matched by id (case number for LOSO, cluster
    # position for LOCO), not by row order, so fits that stack groups
    # differently -- e.g. a pooled fit against a multigroup fit -- still
    # pair up unit by unit.
    pu1 <- loo_list[[1L]]$per_unit
    align <- lapply(loo_list, function(l) {
      if (
        !identical(l$type, loo_list[[1L]]$type) ||
          nrow(l$per_unit) != nrow(pu1)
      ) {
        return(NULL)
      }
      idx <- match(pu1$unit, l$per_unit$unit)
      if (anyNA(idx) || !identical(l$per_unit$nobs[idx], pu1$nobs)) {
        return(NULL)
      }
      idx
    })
    if (any(vapply(align, is.null, logical(1)))) {
      cli_abort(
        "LOO comparison requires models fitted to the same data, with
         matching units."
      )
    }

    # Joint and conditional scores live on different scales and must not
    # meet in one comparison
    flavs <- vapply(
      loo_list,
      function(l) l$flavour %||% "joint",
      character(1)
    )
    if (length(unique(flavs)) > 1L) {
      cli_abort(c(
        "LOO comparison cannot mix joint and conditional scores.",
        "i" = "Fit all models with the same {.code fixed.x} setting."
      ))
    }

    # Joint scores are densities over all modelled variables, so the
    # variable sets must match; conditional scores are densities over the
    # outcomes only, so covariate sets may differ but the outcomes must match
    score_vars <- lapply(internals, function(m) {
      ov <- sort(unique(unlist(m$lavdata@ov.names)))
      if (flavs[1L] == "conditional") {
        setdiff(ov, unlist(m$lavdata@ov.names.x))
      } else {
        ov
      }
    })
    same_vars <- vapply(
      score_vars,
      identical,
      logical(1),
      y = score_vars[[1L]]
    )
    if (!all(same_vars)) {
      if (flavs[1L] == "conditional") {
        cli_abort(
          "Conditional LOO comparison requires models for the same set of
           outcome variables (the covariate sets may differ)."
        )
      }
      cli_abort(c(
        "Joint LOO comparison requires models for the same set of observed
         variables.",
        "i" = "To drop a covariate's effect, keep the variable in the model
         without the path (e.g. {.code fb ~ w1} plus {.code w2 ~~ w1}), or
         compare {.code fixed.x = TRUE} fits, which are scored
         conditionally."
      ))
    }

    # One approximation order across the whole comparison. A model whose
    # units all admit a second-order term would be scored at second order on
    # its own, but pitting that against a model forced to first order would
    # read the change of estimator as a difference between the models, so the
    # comparison drops to the lowest order any of them can supply.
    order_2 <- all(vapply(loo_list, function(l) isTRUE(l$use_second), TRUE))
    n_forced <- sum(!vapply(loo_list, function(l) isTRUE(l$use_second), TRUE))

    elpd <- vapply(
      loo_list,
      function(l) if (order_2) l$elpd_2 else l$elpd_1,
      numeric(1)
    )
    best <- which.max(elpd)
    # The same pointwise contributions the reported ELPDs are summed from,
    # aligned to the first model's unit order for pairing, so that se_diff is
    # the standard error of the elpd_diff actually reported
    pw <- lapply(seq_along(loo_list), function(k) {
      v <- loo_headline_pointwise(loo_list[[k]]$per_unit, order_2)
      v[align[[k]]]
    })
    n_units <- nrow(pu1)

    out$ELPD <- round(elpd, 3)
    out$SE <- round(
      vapply(
        loo_list,
        function(l) if (order_2) l$se_2 else l$se_1,
        numeric(1)
      ),
      3
    )
    out$p_loo <- round(
      vapply(
        loo_list,
        function(l) if (order_2) l$p_loo_2 else l$p_loo_1,
        numeric(1)
      ),
      3
    )
    out$elpd_diff <- round(elpd - elpd[best], 3)
    out$se_diff <- round(
      vapply(
        seq_along(loo_list),
        function(k) {
          if (k == best) {
            return(0)
          }
          sqrt(n_units * var(pw[[k]] - pw[[best]]))
        },
        numeric(1)
      ),
      3
    )
  }

  if (isTRUE(loo)) {
    out <- out[order(-out$ELPD), ]
  } else if (is.null(fit.measures)) {
    out <- out[order(-out$Marg.Loglik), ]
  }

  rownames(out) <- NULL
  attr(out, "fit_measures_used") <- !is.null(fit.measures)
  attr(out, "loo_used") <- isTRUE(loo)
  if (isTRUE(loo)) {
    attr(out, "loo_order") <- if (order_2) 2L else 1L
    attr(out, "loo_n_forced") <- n_forced
    attr(out, "loo_n_models") <- length(loo_list)
  }
  class(out) <- c("compare.inlavaan_internal", class(out))
  out
}

#' @exportS3Method print compare.inlavaan_internal
print.compare.inlavaan_internal <- function(x, ...) {
  cat("Bayesian Model Comparison (INLAvaan)\n")
  if (isTRUE(attr(x, "loo_used"))) {
    ord <- attr(x, "loo_order")
    cat(
      "Models ordered by ELPD (Taylor LOO, ",
      if (identical(ord, 2L)) "second" else "first",
      "-order)\n",
      sep = ""
    )
  } else {
    cat("Models ordered by marginal log-likelihood\n")
  }
  if (isTRUE(attr(x, "loo_used"))) {
    cat("elpd_diff/se_diff are paired differences vs the best model\n")
    n_forced <- attr(x, "loo_n_forced")
    if (!is.null(n_forced) && n_forced > 0L) {
      cat(
        "Scored at first order throughout: ",
        n_forced,
        " of ",
        attr(x, "loo_n_models"),
        " models have units with no second-order term\n",
        sep = ""
      )
    }
  }
  cat("\n")
  print.data.frame(x, row.names = FALSE)
  invisible(x)
}