R/method-fitmeasures.R

In INLAvaan: Approximate Bayesian Latent Variable Analysis

# --- Bayesian fit index helpers -----------------------------------------------

# Saturated log-likelihood (constant for ML; sum across groups). lavaan
# stores it in the fit's h1 slot for every model type it supports, including
# two-level and missing-data fits, so that is read first. The single-level
# formulas below remain as a fallback for objects without an h1 slot. Under
# FIML they use the pattern-based formula to stay on the same scale as
# inlav_model_loglik (which delegates to lavaan___lav_model_loglik).
compute_loglik_sat <- function(object, lavsamplestats, lavdata) {
  h1 <- tryCatch(object@h1$logl$loglik, error = function(e) NULL)
  if (is.numeric(h1) && length(h1) == 1L && is.finite(h1)) {
    return(h1)
  }
  # nocov start
  ngroups <- lavdata@ngroups
  logl_sat <- 0
  for (g in seq_len(ngroups)) {
    if (lavsamplestats@missing.flag) {
      # nocov start
      logl_sat <- logl_sat +
        lavaan___lav_mvn_mi_loglik_samp(
          yp = lavsamplestats@missing[[g]],
          mu = lavsamplestats@mean[[g]],
          sigma_1 = lavsamplestats@cov[[g]],
          x_idx = lavsamplestats@x.idx[[g]],
          x_mean = lavsamplestats@mean.x[[g]],
          x_cov = lavsamplestats@cov.x[[g]]
        )
    } else {
      # nocov end
      logl_sat <- logl_sat +
        lavaan___lav_mvn_loglik_samp(
          sample_mean = lavsamplestats@mean[[g]],
          sample_cov = lavsamplestats@cov[[g]],
          sample_nobs = lavsamplestats@nobs[[g]],
          mu = lavsamplestats@mean[[g]],
          sigma_1 = lavsamplestats@cov[[g]],
          x_idx = lavsamplestats@x.idx[[g]],
          x_mean = lavsamplestats@mean.x[[g]],
          x_cov = lavsamplestats@cov.x[[g]]
        )
    }
  }
  logl_sat
  # nocov end
}

# Per-sample deviance chi-square:  chisq_s = 2 * (loglik_sat - loglik(x_s))
# This equals N * F_ML(x_s).
compute_chisq_dev <- function(
  object,
  x_samp,
  lavmodel,
  lavsamplestats,
  lavdata,
  lavoptions,
  lavcache
) {
  loglik_sat <- compute_loglik_sat(object, lavsamplestats, lavdata)
  vapply(
    seq_len(nrow(x_samp)),
    function(i) {
      ll_i <- inlav_model_loglik(
        x_samp[i, ],
        lavmodel,
        lavsamplestats,
        lavdata,
        lavoptions,
        lavcache
      )
      2 * (loglik_sat - ll_i)
    },
    numeric(1)
  )
}

# Absolute fit indices (vectorised over posterior samples) ---------------------
compute_BRMSEA <- function(nonc, df, N, Ngr) {
  sqrt(nonc / (df * N)) * sqrt(Ngr)
}
compute_BGammaHat <- function(nonc, nvar_total, N) {
  nvar_total / (nvar_total + 2 * nonc / N)
}
compute_adjBGammaHat <- function(BGammaHat, p, df) {
  1 - (p / df) * (1 - BGammaHat)
}
compute_BMc <- function(nonc, N) exp(-0.5 * nonc / N)

# Incremental fit indices (vectorised) ----------------------------------------
compute_BCFI <- function(nonc, nonc_null) 1 - nonc / nonc_null
compute_BTLI <- function(adj_dev, df, adj_dev_null, df_null) {
  tli_null <- adj_dev_null / df_null
  denom <- tli_null - 1
  out <- (tli_null - adj_dev / df) / denom
  # A baseline whose own ratio is 1 leaves nothing to scale by
  out[abs(denom) < 1e-8] <- NA_real_
  out
}
compute_BNFI <- function(adj_dev, adj_dev_null) {
  (adj_dev_null - adj_dev) / adj_dev_null
}

# Reconstruct lavoptions suitable for inlav_model_loglik() from the INLAvaan
# object (whose @Options$estimator was changed to "Bayes").
reconstruct_lavoptions <- function(object) {
  opts <- object@Options
  opts$estimator <- object@external$inlavaan_internal$lavmodel@estimator
  opts
}

# Compute rescaled chi-square, adjusted deviance, df, and noncentrality for
# a single model under the chosen rescaling method.
#
# Returns list(chisq, adj_dev, df, nonc, N_adj, pD) where chisq/adj_dev/nonc
# are per-sample vectors and df/N_adj/pD are scalars.
compute_rescaled_quantities <- function(
  object,
  x_samp,
  lavmodel,
  lavsamplestats,
  lavdata,
  lavoptions,
  lavcache,
  p,
  rescale
) {
  N <- lavsamplestats@ntotal
  Ngr <- lavdata@ngroups
  npar <- object@Fit@npar

  chisq_dev <- compute_chisq_dev(
    object,
    x_samp,
    lavmodel,
    lavsamplestats,
    lavdata,
    lavoptions,
    lavcache
  )

  if (rescale == "devM") {
    pD <- object@external$inlavaan_internal$DIC$pD
    if (is.null(pD) || pD <= 0 || pD >= p) {
      pD <- npar
    } # nocov
    adj_dev <- chisq_dev - pD # obs - reps
    df <- p - pD
    N_adj <- N
  } else {
    # MCMC: use classical chi-square = (N-1)/N * deviance
    pD <- npar
    chisq_dev <- (N - 1) / N * chisq_dev
    adj_dev <- chisq_dev # reps = 0
    df <- p - npar
    N_adj <- N - Ngr # EQS-style: Min1 = TRUE
  }

  nonc <- pmax(adj_dev - df, 0)
  list(
    chisq = chisq_dev,
    adj_dev = adj_dev,
    df = df,
    nonc = nonc,
    N_adj = N_adj,
    pD = pD
  )
}

# ---------------------------------------------------------------------------
# Independence baseline for the incremental indices
# ---------------------------------------------------------------------------

# Keys of the free parameters of a parameter table, for structural equality
free_param_keys <- function(pt) {
  i <- pt$free > 0
  grp <- pt$group %||% rep(1L, length(pt$lhs))
  lvl <- pt$level %||% rep(1L, length(pt$lhs))
  sort(paste(pt$lhs[i], pt$op[i], pt$rhs[i], grp[i], lvl[i]))
}

# TRUE when every free parameter is a variance or an intercept, which is
# what the independence model consists of
is_independence_partable <- function(pt) {
  i <- pt$free > 0
  all((pt$op[i] == "~~" & pt$lhs[i] == pt$rhs[i]) | pt$op[i] == "~1")
}

# Fit the independence (null) model that the incremental indices BCFI, BTLI
# and BNFI are scaled against, on the same data and likelihood options as
# `object`. lavaan writes that model's parameter table, and the fitted
# object's data and sample-statistics slots are reused directly, so no data
# frame is re-read. The indices need only this model's posterior draws and
# its pD, so the marginals are Gaussian and the VB shift is skipped, which
# makes the fit take a fraction of a second even for many items.
fit_independence_baseline <- function(object, nsamp = NULL) {
  int <- object@external$inlavaan_internal
  pt0 <- lavaan::lav_partable_independence(object)
  opt <- object@Options
  inlavaan(
    model = pt0,
    data = NULL,
    slotData = object@Data,
    slotSampleStats = object@SampleStats,
    missing = opt$missing %||% "default",
    meanstructure = opt$meanstructure %||% "default",
    fixed.x = opt$fixed.x %||% "default",
    conditional.x = opt$conditional.x %||% "default",
    likelihood = opt$likelihood %||% "default",
    estimator = int$lavmodel@estimator,
    marginal_method = "marggaus",
    vb_correction = FALSE,
    test = "dic",
    nsamp = nsamp %||% int$nsamp %||% 1000L,
    verbose = FALSE
  )
}

# Resolve the `baseline.model` argument of bfit_indices(): NULL fits the
# independence model (none for a model that already is one), FALSE skips
# the incremental indices, and a supplied fit is checked and used.
resolve_baseline_model <- function(object, baseline.model, nsamp = NULL) {
  if (isFALSE(baseline.model)) {
    return(NULL)
  }
  if (is.null(baseline.model)) {
    if (is_independence_partable(object@ParTable)) {
      return(NULL)
    }
    return(tryCatch(
      fit_independence_baseline(object, nsamp),
      error = function(e) {
        cli_warn(c(
          "Could not fit the independence baseline, so BCFI, BTLI and BNFI
           are not reported.",
          "x" = conditionMessage(e)
        ))
        NULL
      }
    ))
  }
  if (!is(baseline.model, "INLAvaan")) {
    cli_abort(
      "{.arg baseline.model} must be an {.cls INLAvaan} object, or
       {.code FALSE} to skip the incremental indices."
    )
  }
  if (
    identical(
      free_param_keys(object@ParTable),
      free_param_keys(baseline.model@ParTable)
    )
  ) {
    cli_warn(
      "{.arg baseline.model} has the same free parameters as {.arg object},
       so BCFI, BTLI and BNFI are zero by construction."
    )
  }
  baseline.model
}

# ---------------------------------------------------------------------------
# bfit_indices: compute per-sample Bayesian fit index vectors and return
# an S3 object of class "bfit_indices" with a summary() and print() method.
# ---------------------------------------------------------------------------

#' Bayesian Fit Indices
#'
#' Compute posterior distributions of Bayesian fit indices for an INLAvaan
#' model, analogous to [blavaan::blavFitIndices()].
#'
#' @param object An object of class [INLAvaan].
#' @param baseline.model The baseline (null) model that the incremental fit
#'   indices (BCFI, BTLI, BNFI) are scaled against. `NULL` (default) fits the
#'   independence model on the same data and options automatically, as
#'   lavaan does: every observed variable keeps its variance (and intercept)
#'   and nothing correlates. That fit uses Gaussian marginals and no VB
#'   shift, since only its posterior draws and pD are needed, and takes a
#'   fraction of a second. Supply an [INLAvaan] object to use another
#'   baseline, or `FALSE` to skip the incremental indices.
#' @param rescale Character string controlling how the Bayesian chi-square
#'   is rescaled. `"devM"` (default) subtracts pD from the deviance at each
#'   sample. `"MCMC"` uses the classical chi-square and classical df at each
#'   sample.
#' @param nsamp Number of posterior samples to draw. Defaults to the value
#'   used when fitting the model.
#' @param samp_copula Logical. When `TRUE` (default), posterior samples are
#'   drawn using the copula method with the fitted marginals. When `FALSE`,
#'   samples are drawn from the Gaussian (Laplace) approximation.
#' @param x An object of class `bfit_indices` (for `print`).
#' @param ... Additional arguments passed to methods.
#'
#' @returns An S3 object of class `"bfit_indices"` containing:
#' \describe{
#'   \item{`indices`}{Named list of numeric vectors (one per posterior sample)
#'     for each computed fit index.}
#'   \item{`details`}{List with `chisq` (per-sample deviance), `df`, `pD`,
#'     `rescale`, and `nsamp`.}
#' }
#' Use [summary()] to obtain a table of posterior summaries (Mean, SD,
#' quantiles, Mode) for each index.
#'
#' @seealso [lavaan::fitMeasures()], [blavaan::blavFitIndices()],
#'   [fitmeasures()], [compare()]
#'
#' @examples
#' \donttest{
#' HS.model <- "
#'   visual  =~ x1 + x2 + x3
#'   textual =~ x4 + x5 + x6
#'   speed   =~ x7 + x8 + x9
#' "
#' utils::data("HolzingerSwineford1939", package = "lavaan")
#' fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
#'             verbose = FALSE)
#'
#' # Absolute fit indices
#' bf <- bfit_indices(fit)
#' bf
#' summary(bf)
#' }
#'
#' @export
bfit_indices <- function(
  object,
  baseline.model = NULL,
  rescale = c("devM", "MCMC"),
  nsamp = NULL,
  samp_copula = TRUE
) {
  rescale <- match.arg(rescale)
  if (!is(object, "INLAvaan")) {
    cli_abort("{.arg object} must be an {.cls INLAvaan} object.")
  }

  int <- object@external$inlavaan_internal
  lavmodel <- int$lavmodel
  lavsamplestats <- int$lavsamplestats
  lavdata <- int$lavdata

  nsamp <- nsamp %||% int$nsamp %||% 500L
  method <- if (isTRUE(samp_copula)) int$marginal_method else "sampling"
  samp <- sample_params(
    theta_star = int$theta_star,
    Sigma_theta = int$Sigma_theta,
    method = method,
    approx_data = int$approx_data,
    pt = int$partable,
    lavmodel = lavmodel,
    nsamp = nsamp,
    R_star = int$R_star
  )
  x_samp <- samp$x_samp

  if (lavmodel@estimator != "ML") {
    # nocov
    cli_abort("Bayesian fit indices are only supported for the ML estimator.")
  }
  if (rescale == "devM" && !has_test(int, "dic")) {
    cli_abort(
      "DIC not available. Refit with {.arg test} including {.val dic}
       (part of the default {.val standard}), or use
       {.code rescale = \"MCMC\"}."
    )
  }

  lavoptions <- reconstruct_lavoptions(object)
  lavcache <- object@Cache
  N <- lavsamplestats@ntotal
  Ngr <- lavdata@ngroups
  nvar <- lavmodel@nvar
  # Number of sample moments, counted as lavaan counts them for the model's
  # degrees of freedom: per group and level, without the moments of fixed
  # exogenous covariates
  p <- lavaan::lav_partable_ndat(object@ParTable)

  rq <- compute_rescaled_quantities(
    object,
    x_samp,
    lavmodel,
    lavsamplestats,
    lavdata,
    lavoptions,
    lavcache,
    p,
    rescale
  )

  indices <- list()

  if (rq$df > 0) {
    indices$BRMSEA <- compute_BRMSEA(rq$nonc, rq$df, rq$N_adj, Ngr)
    bgh <- compute_BGammaHat(rq$nonc, sum(nvar), rq$N_adj)
    indices$BGammaHat <- bgh
    indices$adjBGammaHat <- compute_adjBGammaHat(bgh, p, rq$df)
    indices$BMc <- compute_BMc(rq$nonc, rq$N_adj)
  } # else: df == 0, no absolute fit indices (nocov – saturated model)

  # Incremental indices, scaled against the independence model unless the
  # caller supplies a baseline or asks to skip them
  baseline.model <- resolve_baseline_model(object, baseline.model, nsamp)
  if (!is.null(baseline.model)) {
    bint <- baseline.model@external$inlavaan_internal

    bmethod <- if (isTRUE(samp_copula)) bint$marginal_method else "sampling"
    samp_null <- sample_params(
      theta_star = bint$theta_star,
      Sigma_theta = bint$Sigma_theta,
      method = bmethod,
      approx_data = bint$approx_data,
      pt = bint$partable,
      lavmodel = bint$lavmodel,
      nsamp = nsamp,
      R_star = bint$R_star
    )
    x_samp_null <- samp_null$x_samp
    n_use <- min(nrow(x_samp), nrow(x_samp_null))

    rq_null <- compute_rescaled_quantities(
      baseline.model,
      x_samp_null[seq_len(n_use), , drop = FALSE],
      bint$lavmodel,
      bint$lavsamplestats,
      bint$lavdata,
      reconstruct_lavoptions(baseline.model),
      baseline.model@Cache,
      p,
      rescale
    )

    adj_dev_use <- rq$adj_dev[seq_len(n_use)]
    nonc_use <- rq$nonc[seq_len(n_use)]

    indices$BCFI <- compute_BCFI(nonc_use, rq_null$nonc)
    indices$BTLI <- compute_BTLI(
      adj_dev_use,
      rq$df,
      rq_null$adj_dev,
      rq_null$df
    )
    indices$BNFI <- compute_BNFI(adj_dev_use, rq_null$adj_dev)
  }

  structure(
    list(
      indices = indices,
      details = list(
        chisq = rq$chisq,
        df = rq$df,
        pD = rq$pD,
        rescale = rescale,
        nsamp = nrow(x_samp)
      )
    ),
    class = "bfit_indices"
  )
}

#' @rdname bfit_indices
#' @exportS3Method summary bfit_indices
#' @usage \method{summary}{bfit_indices}(object, ...)
summary.bfit_indices <- function(object, ...) {
  summ_one <- function(x) {
    x <- x[is.finite(x)]
    if (length(x) < 3) {
      # nocov
      return(rep(NA_real_, 8))
    }
    dens <- stats::density(x)
    c(
      Mean = mean(x),
      SD = stats::sd(x),
      `2.5%` = stats::quantile(x, 0.025, names = FALSE),
      `25%` = stats::quantile(x, 0.250, names = FALSE),
      `50%` = stats::quantile(x, 0.500, names = FALSE),
      `75%` = stats::quantile(x, 0.750, names = FALSE),
      `97.5%` = stats::quantile(x, 0.975, names = FALSE),
      Mode = dens$x[which.max(dens$y)]
    )
  }
  tab <- as.data.frame(do.call(rbind, lapply(object$indices, summ_one)))
  attr(tab, "header") <- paste0(
    "Posterior summary of ",
    object$details$rescale,
    "-based Bayesian fit indices (nsamp = ",
    object$details$nsamp,
    "):"
  )
  class(tab) <- c("lavaan.data.frame", "data.frame")
  tab
}

#' @rdname bfit_indices
#' @exportS3Method print bfit_indices
print.bfit_indices <- function(x, ...) {
  tab <- summary(x)
  hdr <- attr(tab, "header")
  if (!is.null(hdr)) {
    cat(hdr, "\n\n")
  }
  eap <- vapply(x$indices, mean, numeric(1), na.rm = TRUE)
  print(round(eap, 3))
  invisible(x)
}

# --- Main fitMeasures function ------------------------------------------------

inlav_fit_measures <- function(
  object,
  fit.measures = "all",
  baseline.model = NULL,
  ...
) {
  dots <- list(...)
  rescale <- match.arg(dots$rescale %||% "devM", c("devM", "MCMC"))

  # Has the model converged?
  if (object@Fit@npar > 0L & !object@optim$converged) {
    # nocov
    cli_alert_warning("Optimiser did not converge.") # nocov
  } # nocov

  out <- vector("numeric")
  out["npar"] <- object@Fit@npar
  out["margloglik"] <- object@external$inlavaan_internal$mloglik

  int <- object@external$inlavaan_internal
  if (has_test(int, "ppp")) {
    out["ppp"] <- int$ppp
  }
  if (has_test(int, "dic")) {
    out["dic"] <- int$DIC$dic
    out["p_dic"] <- int$DIC$pD
  }

  # Validate baseline.model early (before tryCatch)
  if (
    !is.null(baseline.model) &&
      !isFALSE(baseline.model) &&
      !is(baseline.model, "INLAvaan")
  ) {
    cli_abort(
      "{.arg baseline.model} must be an {.cls INLAvaan} object, or
       {.code FALSE} to skip the incremental indices."
    )
  }

  # The independence baseline is fitted only when an incremental index is
  # actually wanted, so a request for absolute indices alone stays free
  incr_measures <- c("BCFI", "BTLI", "BNFI")
  need_incr <- identical(fit.measures, "all") ||
    any(incr_measures %in% fit.measures)
  if (is.null(baseline.model) && !need_incr) {
    baseline.model <- FALSE
  }

  # Bayesian fit indices (BRMSEA, BGammaHat, etc.)
  bfi <- tryCatch(
    bfit_indices(object, baseline.model, rescale),
    error = function(e) NULL
  )
  if (!is.null(bfi)) {
    for (nm in names(bfi$indices)) {
      out[nm] <- mean(bfi$indices[[nm]], na.rm = TRUE)
    }
  }

  # LOO measures: free when stored with the fit (test = "loo"/"waic"/"full",
  # or add_loo()); otherwise computed on demand, and only when requested by
  # name -- the LOO
  # computation is fresh work and cannot be cached into `object` (S4 copy
  # semantics), so it never silently inflates a bare fitMeasures() call
  loo_measures <- c("elpd_loo", "se_loo", "p_loo", "looic")
  res_loo <- object@external$inlavaan_internal$loo
  if (
    is.null(res_loo) &&
      !identical(fit.measures, "all") &&
      any(loo_measures %in% fit.measures)
  ) {
    res_loo <- tryCatch(loo(object), error = function(e) NULL)
  }
  if (!is.null(res_loo)) {
    out["elpd_loo"] <- res_loo$estimates["elpd_loo", "Estimate"]
    out["p_loo"] <- res_loo$estimates["p_loo", "Estimate"]
    out["looic"] <- res_loo$estimates["looic", "Estimate"]
    out["se_loo"] <- res_loo$estimates["looic", "SE"]
    # The order used is a property of the numbers, not of when they were
    # computed. A stored LOO warned at fit time, so without this a reloaded
    # fit would hand back first-order measures without a word.
    n_loo <- res_loo$n_units
    if (isTRUE(res_loo$second_order) && res_loo$n_ok < n_loo) {
      n_bad <- n_loo - res_loo$n_ok
      cli_warn(c(
        "{n_bad} of {n_loo} units {qty(n_bad)}{?has/have} no second-order
         term ({.field k_max} at or above 1).",
        "i" = "{.field elpd_loo}, {.field p_loo} and {.field looic} are
               reported at first order. See {.fun loo}."
      ))
    }
  }

  # WAIC: free when stored with the fit; otherwise computed on demand (a
  # full Taylor pass) when requested by name
  waic_measures <- c("elpd_waic", "se_waic", "p_waic", "waic")
  res_waic <- object@external$inlavaan_internal$waic
  if (
    is.null(res_waic) &&
      !identical(fit.measures, "all") &&
      any(waic_measures %in% fit.measures)
  ) {
    res_waic <- tryCatch(waic(object), error = function(e) NULL)
  }
  if (!is.null(res_waic)) {
    out["elpd_waic"] <- res_waic$estimates["elpd_waic", "Estimate"]
    out["p_waic"] <- res_waic$estimates["p_waic", "Estimate"]
    out["waic"] <- res_waic$estimates["waic", "Estimate"]
    out["se_waic"] <- res_waic$estimates["waic", "SE"]
  }

  # Filter if specific measures requested
  if (!identical(fit.measures, "all")) {
    idx <- which(names(out) %in% fit.measures)
    if (length(idx) == 0L) {
      cli_abort("No matching fit measures found.")
    }
    out <- out[idx]
  }

  class(out) <- c("fitmeasures.inlavaan_internal", "numeric")
  out
}

#' @exportS3Method print fitmeasures.inlavaan_internal
print.fitmeasures.inlavaan_internal <- function(x, ...) {
  nm <- names(x)

  # Apply conditional formatting
  formatted_values <- sapply(seq_along(x), function(i) {
    val <- x[i]
    name <- nm[i]

    if (name == "npar") {
      return(as.character(as.integer(round(val))))
    } else if (startsWith(name, "grad_")) {
      # nocov start
      # Use formatC to force scientific and maintain 3 significant digits
      return(formatC(val, digits = 2, format = "e"))
    } else {
      # nocov end
      # Round to 3 decimal places
      return(formatC(val, digits = 3, format = "f", drop0trailing = FALSE))
    }
  })

  # Set names back onto the formatted character vector
  names(formatted_values) <- nm

  # Print using the standard named vector style without quotes
  print(formatted_values, quote = FALSE, right = TRUE)

  invisible(x)
}

#' Fit Measures for a Latent Variable Model estimated using INLA
#'
#' @param object An object of class [INLAvaan].
#' @param fit_measures If `"all"`, all fit measures available will be returned. If
#'   only a single or a few fit measures are specified by name, only those are
#'   computed and returned. The LOO measures `"elpd_loo"`, `"se_loo"`,
#'   `"p_loo"` and `"looic"` (see [loo()]), and the WAIC measures
#'   `"elpd_waic"`, `"se_waic"`, `"p_waic"` and `"waic"` (see [waic()]), are
#'   included in `"all"` only when stored with the fit (`test` including
#'   `"loo"`, `"waic"` or `"full"` in [inlavaan()], or [add_loo()]);
#'   otherwise they are computed on demand when requested by name, and
#'   recomputed on every call -- store the result with `fit <- add_loo(fit)`
#'   (or call [loo()]/[waic()] directly) for repeated access. INLAvaan's
#'   stable spelling `fit.measures` is also accepted.
#' @param baseline_model The baseline (null) model for the incremental fit
#'   indices (BCFI, BTLI, BNFI). `NULL` (default) fits the independence model
#'   automatically, as lavaan does. Supply an [INLAvaan] object to use
#'   another baseline, or `FALSE` to skip the incremental indices. INLAvaan's
#'   stable spelling `baseline.model` is also accepted; see [bfit_indices()].
#' @param h1_model Ignored (included for compatibility with the lavaan
#'   generic).
#' @param fm_args Ignored (included for compatibility with the lavaan
#'   generic).
#' @param output Ignored (included for compatibility with the lavaan
#'   generic).
#' @param level Ignored (included for compatibility with the lavaan
#'   generic).
#' @param ... Additional arguments. Currently supports:
#' \describe{
#'   \item{`rescale`}{Character string controlling how the Bayesian chi-square
#'     is computed, following `blavaan::blavFitIndices()`. Options are `"devM"`
#'     (default) which uses the deviance rescaled by `pD` from DIC, or `"MCMC"`
#'     which uses the classical chi-square (`(N-1) * F_ML`) and classical
#'     degrees of freedom (`p - npar`) at each posterior sample.}
#' }
#'
#' @returns A named numeric vector of fit measures.
#'
#' @seealso [bfit_indices()], [compare()], [diagnostics()]
#'
#' @examples
#' \donttest{
#' HS.model <- "
#'   visual  =~ x1 + x2 + x3
#'   textual =~ x4 + x5 + x6
#'   speed   =~ x7 + x8 + x9
#' "
#' utils::data("HolzingerSwineford1939", package = "lavaan")
#' fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
#'             verbose = FALSE)
#'
#' # All available fit measures
#' fitMeasures(fit)
#'
#' # Specific measures
#' fitMeasures(fit, c("npar", "dic", "p_dic", "ppp"))
#' }
#'
#' @usage
#' \S4method{fitMeasures}{INLAvaan}(object, fit_measures = "all",
#'   baseline_model = NULL, h1_model = NULL, fm_args,
#'   output = "vector", level = NULL, ...)
#'
#' @importMethodsFrom lavaan fitmeasures fitMeasures
#' @name fitmeasures
#' @rdname fitmeasures
#' @aliases fitMeasures,INLAvaan-method
#' @rawNamespace exportMethods(fitMeasures)
NULL

#' @usage
#' \S4method{fitmeasures}{INLAvaan}(object, fit_measures = "all",
#'   baseline_model = NULL, h1_model = NULL, fm_args,
#'   output = "vector", level = NULL, ...)
#'
#' @name fitmeasures
#' @rdname fitmeasures
#' @aliases fitmeasures,INLAvaan-method
#' @rawNamespace exportMethods(fitmeasures)
NULL

# With lavaan >= 0.7 required, the generics' argument spellings are fixed
# (fit_measures, baseline_model), so the methods can be registered the
# ordinary way at build time; the load-time re-registration this replaced
# existed only to serve both lavaan generations at once.
#
# INLAvaan's own documented/stable parameter names remain fit.measures and
# baseline.model (used internally, e.g. by compare()). Those don't match the
# generic's formals, so a caller using them falls through into "..."
# unmatched rather than being silently dropped; recover them from there.
inlav_fitmeasures_method <- function(
  object,
  fit_measures = "all",
  baseline_model = NULL,
  h1_model = NULL,
  fm_args = list(
    standard.test = "default",
    scaled.test = "default",
    rmsea.ci.level = 0.90,
    rmsea.close.h0 = 0.05,
    rmsea.notclose.h0 = 0.08,
    robust = TRUE,
    cat.nonpd = "na"
  ),
  output = "vector",
  level = NULL,
  ...
) {
  dots <- list(...)
  if ("fit.measures" %in% names(dots)) {
    fit_measures <- dots[["fit.measures"]]
    dots[["fit.measures"]] <- NULL
  }
  if ("baseline.model" %in% names(dots)) {
    baseline_model <- dots[["baseline.model"]]
    dots[["baseline.model"]] <- NULL
  }
  do.call(
    inlav_fit_measures,
    c(
      list(
        object,
        fit.measures = fit_measures,
        baseline.model = baseline_model
      ),
      dots
    )
  )
}

setMethod("fitMeasures", "INLAvaan", inlav_fitmeasures_method)
setMethod("fitmeasures", "INLAvaan", inlav_fitmeasures_method)