R/score_cfc_lavaan.R

In autoFC: Automatic Toolkit for Construction, Optimization, Scoring and Simulation of Forced-Choice Tests

#' Full-Information Closed-Form Trait Scoring for Compositional FC Models
#'
#' @description Calculates exact, full-information MAP trait scores and Standard Errors
#' for a fitted compositional forced-choice model, preserving all within-block local dependencies.
#'
#' @param fit A fitted lavaan object from a compositional model, generated from \code{cfa()} function
#' with the syntax built from \code{generate_cfc_lavaan_syntax()}.
#' @param data The log-ratio dataset used to fit the model.
#' @param trait_names Optional character vector of the traits to score.
#'
#' @return A data frame containing the exact trait scores, followed by all standard errors (SE) of the traits.
#' @export
score_cfc_lavaan <- function(fit, data, trait_names = NULL) {

  # 1. Extract Parameter Matrices
  est <- lavaan::lavInspect(fit, "est")

  y_names <- lavaan::lavNames(fit, "ov")
  all_lv <- lavaan::lavNames(fit, "lv")
  utils_names <- grep("^t[0-9]+$", all_lv, value = TRUE)

  if (is.null(trait_names)) {
    trait_names <- setdiff(all_lv, utils_names)
  }
  D <- length(trait_names)
  N <- nrow(data)

  # Extract loadings and build marginal loading matrix: Lambda_tilde = L_yt * B_teta
  L_yt <- est$lambda[y_names, utils_names, drop = FALSE]
  L_teta <- est$beta[utils_names, trait_names, drop = FALSE]
  Lambda_tilde <- L_yt %*% L_teta

  # Extract utility residual covariance (Psi_t) and observed intercepts (nu_y)
  Psi_t <- est$psi[utils_names, utils_names, drop = FALSE]
  nu_y <- est$nu[y_names, 1, drop = FALSE]

  # Extract Trait Prior Covariance
  Omega <- est$psi[trait_names, trait_names, drop = FALSE]
  Omega_inv <- solve(Omega)

  data_mat <- as.matrix(data[, y_names, drop = FALSE])

  # Pre-allocate output matrix
  clean_scores <- matrix(NA, nrow = N, ncol = D * 2)
  colnames(clean_scores) <- c(trait_names, paste0(trait_names, "_SE"))

  # 2. Closed-Form Full-Information Calculations
  for (i in 1:N) {
    y_vec <- data_mat[i, ]
    valid <- !is.na(y_vec)
    if (sum(valid) == 0) next

    y_v <- y_vec[valid]
    L_v <- Lambda_tilde[valid, , drop = FALSE]
    L_yt_v <- L_yt[valid, , drop = FALSE]
    nu_v <- nu_y[valid, 1, drop = FALSE]

    # --- Exact Local Dependency Math ---
    # Implied covariance of observed variable residuals: Theta_implied = L_yt * Psi_t * L_yt'
    Theta_implied <- L_yt_v %*% Psi_t %*% t(L_yt_v)

    # We add a tiny ridge to the diagonal to ensure perfect numerical stability during inversion
    diag(Theta_implied) <- diag(Theta_implied) + 1e-9
    Theta_implied_inv <- solve(Theta_implied)

    # Posterior Covariance matrix: Sigma = (Prior_inv + L' * Theta_implied_inv * L)^-1
    Sigma_post <- solve(Omega_inv + t(L_v) %*% Theta_implied_inv %*% L_v)

    # Posterior Mean (MAP scores): Mu = Sigma * L' * Theta_implied_inv * (y - intercept)
    mu_post <- Sigma_post %*% t(L_v) %*% Theta_implied_inv %*% (y_v - nu_v)

    # Standard Errors
    se_post <- sqrt(diag(Sigma_post))

    clean_scores[i, 1:D] <- as.numeric(mu_post)
    clean_scores[i, (D + 1):(D * 2)] <- se_post
  }

  return(as.data.frame(clean_scores))
}