R/generate_cfc_lavaan_syntax.R

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

#' Generate lavaan Syntax for Compositional Forced-Choice (CFC) Models
#'
#' @description Generates the syntax required to fit a continuous Compositional
#' Forced-Choice model (log-ratio differences of utilities) in lavaan,
#' strictly adhering to Brown's (2016) formulation.
#'
#' @param n_blocks Integer. Number of blocks.
#' @param block_size Integer. Number of items per block.
#' @param key_matrix A data.frame indicating item trait and item sign.
#' @param trait_col Character string. The name of the column in \code{key_matrix} indicating the trait
#'        measured by the item.
#' @param key_col Character string. The name of the column in \code{key_matrix} indicating the keying
#'        direction of the item (e.g., positive/negative).
#' @param cor_matrix Optional matrix. Starting values for trait correlations.
#' @param force_positive_variances Logical. If TRUE, prevents Heywood cases by forcing
#'        the utility uniquenesses to be > 0.001. Defaults to TRUE.
#'
#' @return A character string containing the lavaan model syntax.
#' @export
generate_cfc_lavaan_syntax <- function(n_blocks,
                                       block_size,
                                       key_matrix,
                                       trait_col,
                                       key_col,
                                       cor_matrix = NULL,
                                       force_positive_variances = TRUE) {

  n_items <- n_blocks * block_size
  n_obs <- n_blocks * (block_size - 1) # N-1 observed outcomes per block

  trait_names <- unique(as.character(key_matrix[[trait_col]]))
  n_traits <- length(trait_names)

  if (is.null(cor_matrix)) {
    cor_matrix <- diag(n_traits)
    rownames(cor_matrix) <- trait_names; colnames(cor_matrix) <- trait_names
  }

  # Build Observed Variable Names (e.g., y1y3, y2y3)
  obs_names <- character(n_obs)
  obs_idx <- 1
  for (b in 1:n_blocks) {
    ref_item <- b * block_size
    for (i in 1:(block_size - 1)) {
      item_i <- (b - 1) * block_size + i
      obs_names[obs_idx] <- paste0("y", item_i, "y", ref_item)
      obs_idx <- obs_idx + 1
    }
  }

  syn <- c("# --- COMPOSITIONAL FORCED-CHOICE MODEL (Brown, 2016) ---")

  # 1. First-Order Utilities (Defining the contrasts)
  syn <- c(syn, "\n# 1. First-Order Utilities (y* = t_i - t_k)")
  for (b in 1:n_blocks) {
    ref_item <- b * block_size
    block_obs_names <- character(block_size - 1)

    # Non-reference items load 1 onto their respective log-ratio
    for (i in 1:(block_size - 1)) {
      item_i <- (b - 1) * block_size + i
      y_name <- paste0("y", item_i, "y", ref_item)
      block_obs_names[i] <- y_name
      syn <- c(syn, paste0("t", item_i, " =~ 1 * ", y_name))
    }

    # Reference item loads -1 onto all log-ratios in its block
    terms_ref <- paste0("-1 * ", block_obs_names, collapse = " + ")
    syn <- c(syn, paste0("t", ref_item, " =~ ", terms_ref))
  }

  # 2. Second-Order Structure (Traits load on Utilities)
  syn <- c(syn, "\n# 2. Second-Order Structure (Traits load onto Utilities)")
  for (t_idx in 1:n_traits) {
    t_name <- trait_names[t_idx]
    items_for_t <- which(as.character(key_matrix[[trait_col]]) == t_name)

    terms <- sapply(items_for_t, function(i) {
      sign_val <- key_matrix[[key_col]][i]
      paste0("NA*t", i, " + start(", sign_val, ")*t", i)
    })
    syn <- c(syn, paste0(t_name, " =~ ", paste(terms, collapse = " + ")))
  }

  # 3. Trait Variances & Covariances
  syn <- c(syn, "\n# 3. Trait Variances (Fixed to 1) and Covariances")
  for (t_name in trait_names) {
    syn <- c(syn, paste0(t_name, " ~~ 1 * ", t_name))
  }
  if (n_traits > 1) {
    for (i in 1:(n_traits - 1)) {
      for (j in (i + 1):n_traits) {
        t1 <- trait_names[i]; t2 <- trait_names[j]
        syn <- c(syn, paste0(t1, " ~~ start(", cor_matrix[t1, t2], ") * ", t2))
      }
    }
  }

  # 4. Utility Residual Variances (Freely Estimated - See Footnote 2!)
  syn <- c(syn, "\n# 4. Utility Residual Variances (Freely Estimated)")
  for (i in 1:n_items) {
    if (force_positive_variances) {
      syn <- c(syn, paste0("t", i, " ~~ var_t", i, " * t", i))
    } else {
      syn <- c(syn, paste0("t", i, " ~~ t", i))
    }
  }

  # 5. Observed Variable Residual Variances (Fixed to 0) & Intercepts
  syn <- c(syn, "\n# 5. Observed Variable Residual Variances (Fixed to 0) & Intercepts")
  for (y_name in obs_names) {
    syn <- c(syn, paste0(y_name, " ~~ 0 * ", y_name)) # No pair-specific noise
    syn <- c(syn, paste0(y_name, " ~ 1"))            # Estimate intercepts (Gamma)
  }

  # 6. Prevent Heywood Cases (Applied to Utilities, not the observed variables!)
  if (force_positive_variances) {
    syn <- c(syn, "\n# 6. Prevent Heywood Cases (Utility Variances > 0)")
    for (i in 1:n_items) {
      syn <- c(syn, paste0("var_t", i, " > 0.0001"))
    }
  }

  return(paste(syn, collapse = "\n"))
}