R/sim.data.FCGDINA.R

In ForceChoice: Forced-Choice Modeling Based on Item Response Theory and Cognitive Diagnostic Models

#' Simulate Forced-Choice GDINA Data
#'
#' @description
#' Generates forced-choice response data from a cognitive diagnostic item
#' model. The item-level model may be DINA, DINO, ACDM, or GDINA. The resulting
#' statement endorsement probabilities are transformed into forced-choice
#' response probabilities for \code{"RANK"}, \code{"MOLE"}, or \code{"PICK"}
#' blocks using the same sequential Luce--Plackett mechanism used by
#' \code{\link{sim.data.FCMIRT}} and \code{\link{sim.data.FCGGUM}}.
#'
#' @section Data-Generating Process:
#'
#' The simulation proceeds in four steps:
#' \enumerate{
#'   \item Generate or validate the Q-matrix and forced-choice block
#'         structure.
#'   \item Generate latent attribute profiles \eqn{\boldsymbol{\alpha}_n}.
#'   \item Generate item-level CDM endorsement probabilities using guessing
#'         and slipping-style bounds \code{P0} and \code{P1}; then recover the
#'         corresponding CDM delta parameters from the design matrix.
#'   \item Convert item endorsement probabilities into block-level
#'         forced-choice pattern probabilities and sample observed responses.
#' }
#'
#' Attribute profiles can be generated from a higher-order latent trait model,
#' a thresholded multivariate normal model, a uniform distribution over all
#' attribute patterns, or supplied directly through \code{control$alpha}.
#'
#' @param N.person Number of persons.
#' @param N.block Number of forced-choice blocks.
#' @param I.block Number of items per block when \code{control$block.items} is
#'   not supplied. Must be at least 2.
#' @param D Number of latent attributes. If \code{control$Q.matrix} is supplied,
#'   \code{D} is reset to \code{ncol(control$Q.matrix)}.
#' @param model CDM item model: \code{"DINA"}, \code{"DINO"}, \code{"ACDM"}, or
#'   \code{"GDINA"}.
#' @param fc.type Forced-choice response type: \code{"RANK"}, \code{"MOLE"}, or
#'   \code{"PICK"}. A scalar value is recycled to all blocks; a vector of
#'   length \eqn{B} may be supplied for mixed block formats.
#' @param control Optional named list controlling the data-generating process.
#'   Supported entries are described below.
#'
#' @section Control List:
#'
#' \describe{
#'   \item{\code{Q.matrix}}{Optional \eqn{I \times D} binary Q-matrix. If
#'         omitted, a Q-matrix is generated by \code{sim.data.Q.CDM()}.}
#'   \item{\code{block.items}}{Optional list of item indices per block. If
#'         omitted, blocks are formed sequentially using \code{N.block} and
#'         \code{I.block}.}
#'   \item{\code{single}}{Logical passed to \code{sim.data.Q.CDM()} when
#'         \code{Q.matrix} is generated internally. Default \code{TRUE}.}
#'   \item{\code{alpha}}{Optional \eqn{N \times D} binary matrix of true
#'         attribute profiles. If supplied, latent attributes are not
#'         generated.}
#'   \item{\code{distribution}}{Latent attribute distribution when
#'         \code{alpha} is not supplied. Options are \code{"horder"}
#'         (default; alias \code{"higher.order"}), \code{"mvnorm"}, and
#'         \code{"uniform"}.}
#'   \item{\code{delta1}, \code{delta0}, \code{theta}}{Higher-order
#'         discrimination, threshold, and person trait values used when
#'         \code{distribution = "horder"}. Defaults are generated from
#'         log-normal, normal, and standard normal distributions.}
#'   \item{\code{Corr}, \code{threshold}}{Correlation matrix and thresholds
#'         used when \code{distribution = "mvnorm"}.}
#'   \item{\code{gs}}{Optional \eqn{I \times 2} matrix with columns
#'         \code{P0} and \code{P1}; \code{P0} is the lower endorsement
#'         probability and \code{P1} is the upper endorsement probability for
#'         each statement. Defaults are sampled from \code{U(0, .2)} and
#'         \code{U(.8, 1)}.}
#'   \item{\code{mono.constraint}}{Logical; whether to enforce monotonic item
#'         probabilities for GDINA item generation. Default \code{TRUE}.}
#' }
#'
#' @return An object of class \code{"data.FCGDINA"}, a list containing:
#' \describe{
#'   \item{\code{data}}{\eqn{N \times B} character matrix of forced-choice
#'         responses, such as \code{"1>3>2"} for ranking blocks.}
#'   \item{\code{response}}{\eqn{N \times B} integer matrix of 1-based pattern
#'         indices matching \code{patterns}.}
#'   \item{\code{alpha}}{\eqn{N \times D} binary matrix of true latent
#'         attribute profiles.}
#'   \item{\code{gs}}{\eqn{I \times 2} matrix of item probability bounds
#'         \code{P0} and \code{P1}.}
#'   \item{\code{delta.list}}{List of true CDM delta vectors, one per
#'         statement.}
#'   \item{\code{design.matrix.list}}{List of CDM design matrices used to map
#'         delta parameters to item endorsement probabilities.}
#'   \item{\code{Q.matrix}}{\eqn{I \times D} Q-matrix.}
#'   \item{\code{block.items}}{List of statement indices per forced-choice
#'         block.}
#'   \item{\code{model}}{CDM item model used for simulation.}
#'   \item{\code{patterns}}{List of observed response-pattern matrices per
#'         block.}
#'   \item{\code{patterns.total}}{List of full-ranking pattern matrices per
#'         block.}
#'   \item{\code{prob.item}}{\eqn{N \times I} matrix of item-level endorsement
#'         probabilities for sampled persons.}
#'   \item{\code{prob.states}}{\eqn{I \times 2^D} matrix of item endorsement
#'         probabilities for every attribute pattern.}
#'   \item{\code{prob}}{\eqn{N \times \sum_b P_b} matrix of forced-choice
#'         response-pattern probabilities.}
#'   \item{\code{N.person}, \code{N.block}, \code{I.block}, \code{D},
#'         \code{I.states}, \code{fc.type}}{Final data-generating dimensions
#'         and response type.}
#'   \item{\code{call}, \code{arguments}}{Matched call and effective
#'         simulation arguments.}
#' }
#'
#' @references
#' de la Torre, J. (2011). The generalized DINA model framework.
#' \emph{Psychometrika}, 76(2), 179--199.
#' \doi{10.1007/s11336-011-9207-7}
#'
#' @seealso \code{\link{fit.FCGDINA}}, \code{\link{model.FCGDINA}},
#'   \code{\link{sim.data.FCMIRT}}, \code{\link{sim.data.FCGGUM}},
#'   \code{\link{sim.data.FCDCM}}
#'
#' @examples
#' set.seed(123)
#' sim <- sim.data.FCGDINA(N.person = 20, N.block = 3, I.block = 2,
#'                         D = 2, model = "DINA", fc.type = "RANK")
#'
#' str(sim$data)
#' head(sim$response)
#' head(sim$alpha)
#' sim$Q.matrix
#'
#' # Mixed response formats across blocks.
#' sim.mix <- sim.data.FCGDINA(
#'   N.person = 12, N.block = 3, I.block = 3, D = 2,
#'   model = "GDINA", fc.type = c("RANK", "MOLE", "PICK")
#' )
#' vapply(sim.mix$patterns, nrow, integer(1))
#'
#' @export
sim.data.FCGDINA <- function(N.person = 1000, N.block = 10, I.block = 2,
                             D = 3, model = "GDINA", fc.type = "RANK",
                             control = NULL) {

  call <- match.call()

  if (is.null(control)) {
    control <- list()
  }

  # ---- Model validation ----
  model <- toupper(model[1L])
  model.valid <- c("DINA", "DINO", "ACDM", "GDINA")
  if (!model %in% model.valid) {
    stop("'model' must be one of: ", paste(model.valid, collapse = ", "), call. = FALSE)
  }

  # ---- Scalar integer validation ----
  N.person <- check_integer_scalar(N.person, "N.person", 1L)
  D        <- check_integer_scalar(D,        "D",        1L)
  I.block  <- check_integer_scalar(I.block,  "I.block",  2L)
  N.block  <- check_integer_scalar(N.block,  "N.block",  1L)

  # ---- Q-matrix and block configuration ----
  Q.matrix    <- control$Q.matrix
  block.items <- control$block.items

  res <- resolve_fc_blocks(
    Q.matrix       = Q.matrix,
    block.items    = block.items,
    I.block        = I.block,
    N.block        = N.block,
    D              = D,
    q.valid.values = c(0, 1),
    q.row.check    = function(r) any(r > 0)
  )
  I.states   <- res$I.states
  N.block    <- res$N.block
  I.block    <- res$I.block
  block.items <- res$block.items
  Q.matrix   <- res$Q.matrix
  D          <- res$D

  # ---- fc.type validation ----
  fc.type <- normalize_fc_type(fc.type, N.block)

  # ---- Generate Q-matrix if not supplied ----
  if (is.null(Q.matrix)) {
    single <- get_ctrl("single", TRUE, control)
    Q.matrix <- sim.data.Q.CDM(D, I.states, single = single)
  }

  # ---- Alpha generation ----
  alpha <- control$alpha
  if (!is.null(alpha)) {
    alpha <- as.matrix(alpha)
    if (nrow(alpha) != N.person || ncol(alpha) != D ||
        anyNA(alpha) || !all(alpha %in% c(0, 1))) {
      stop("'control$alpha' must be a ", N.person, " x ", D, " binary matrix.", call. = FALSE)
    }
    storage.mode(alpha) <- "integer"
  } else {
    distribution <- tolower(get_ctrl("distribution", "horder", control))
    if (distribution == "higher.order") {
      distribution <- "horder"
    }
    if (!distribution %in% c("horder", "mvnorm", "uniform")) {
      stop("'control$distribution' must be one of: horder, mvnorm, uniform.", call. = FALSE)
    }

    if (distribution == "horder") {
      delta1 <- as.numeric(get_ctrl("delta1", rlnorm(D, meanlog = 0.25, sdlog = 0.25), control))
      delta0 <- as.numeric(get_ctrl("delta0", rnorm(D), control))
      theta  <- as.numeric(get_ctrl("theta",  rnorm(N.person), control))

      if (length(delta1) != D) {
        stop("'control$delta1' must have length ", D, ".", call. = FALSE)
      }
      if (length(delta0) != D) {
        stop("'control$delta0' must have length ", D, ".", call. = FALSE)
      }
      if (length(theta) != N.person) {
        stop("'control$theta' must have length ", N.person, ".", call. = FALSE)
      }

      prob.alpha <- distribution_higher_order(
        theta = theta, delta1 = delta1, delta0 = delta0
      )
      alpha <- (prob.alpha > matrix(runif(N.person * D), N.person, D)) * 1L
    } else if (distribution == "mvnorm") {
      Corr <- validate_corr_matrix(control$Corr, D)
      threshold <- as.numeric(get_ctrl("threshold", seq_len(D) / (D + 1), control))
      if (length(threshold) != D) {
        stop("'control$threshold' must have length ", D, ".", call. = FALSE)
      }
      alpha <- distribution_multi_normal(N.person, Corr, threshold)
    } else {
      patterns.all <- attributepattern(D)
      alpha <- distribution_uniform(N.person, patterns.all)
    }
  }

  # ---- Guessing/slipping parameters ----
  gs <- control$gs
  if (is.null(gs)) {
    gs <- cbind(
      P0 = runif(I.states, 0.00, 0.20),
      P1 = runif(I.states, 0.80, 1.00)
    )
  } else {
    gs <- as.matrix(gs)
    if (nrow(gs) != I.states || ncol(gs) != 2L ||
        anyNA(gs) || !is.numeric(gs)) {
      stop("'control$gs' must be an I.states x 2 numeric matrix.", call. = FALSE)
    }
    colnames(gs) <- c("P0", "P1")
  }

  if (any(gs[, 1L] > gs[, 2L])) {
    stop("Each row of 'gs' must satisfy P0 <= P1.", call. = FALSE)
  }

  # ---- Item-level endorsement probabilities and delta parameters ----
  patterns.full <- attributepattern(D)
  pattern.names <- pattern_key(patterns.full)

  prob.states <- matrix(
    NA_real_,
    nrow = I.states,
    ncol = nrow(patterns.full),
    dimnames = list(paste0("item", seq_len(I.states)), pattern.names)
  )

  delta.list        <- vector("list", I.states)
  design.matrix.list <- vector("list", I.states)

  for (i in seq_len(I.states)) {
    att <- which(Q.matrix[i, ] > 0)
    patterns.i <- if (length(att) > 0) attributepattern(length(att)) else matrix(0L, 1, 0)

    # Generate valid CDM item response probabilities
    mono <- get_ctrl("mono.constraint", TRUE, control)
    p.i <- make_item_prob_cdm(
      patterns        = patterns.i,
      model           = model,
      p0              = gs[i, 1L],
      p1              = gs[i, 2L],
      mono.constraint = mono
    )

    # Build design matrix and back-solve for delta (true model parameters)
    X.i <- get_design_matrix_cdm(patterns.i, model)
    delta.i <- qr.solve(X.i, p.i)

    # Map reduced-pattern probabilities back to full pattern space
    if (length(att) > 0) {
      local.key <- pattern_key(patterns.full[, att, drop = FALSE])
    } else {
      local.key <- rep("", nrow(patterns.full))
    }
    prob.states[i, ] <- p.i[local.key]

    delta.list[[i]]        <- delta.i
    design.matrix.list[[i]] <- X.i
  }

  alpha.names <- pattern_key(alpha)
  alpha.index <- match(alpha.names, pattern.names)

  prob.item <- t(prob.states[, alpha.index, drop = FALSE])
  colnames(prob.item) <- paste0("item", seq_len(I.states))

  # ---- FC permutation patterns ----
  pat <- generate_fc_permutation_patterns(block.items, fc.type)
  patterns.total <- pat$patterns.total
  patterns       <- pat$patterns

  # ---- FC pattern probabilities ----
  prob <- model.FCGDINA(
    alpha = alpha,
    delta.list = delta.list,
    design.matrix.list = design.matrix.list,
    Q.matrix = Q.matrix,
    patterns.total = patterns.total,
    patterns = patterns,
    prob.item = prob.item
  )

  # ---- Response sampling ----
  samples <- sample_fc_responses(prob, patterns, N.person, N.block)
  response <- samples$response
  data     <- samples$data

  # ---- Row/column names ----
  rownames(alpha) <- paste0("person", seq_len(N.person))
  colnames(alpha) <- paste0("alpha", seq_len(D))
  rownames(response) <- rownames(data) <- rownames(alpha)
  colnames(response) <- colnames(data) <- paste0("block.", seq_len(N.block))

  rownames(Q.matrix) <- paste0("item", seq_len(I.states))
  colnames(Q.matrix) <- paste0("dim", seq_len(D))

  rownames(gs) <- paste0("item", seq_len(I.states))

  rownames(prob.item) <- rownames(alpha)
  rownames(prob.states) <- paste0("item", seq_len(I.states))

  # ---- Output assembly ----
  data.obj <- list(
    data              = data,
    response          = response,
    alpha             = alpha,
    gs                = gs,
    delta.list        = delta.list,
    design.matrix.list = design.matrix.list,
    Q.matrix          = Q.matrix,
    block.items       = block.items,
    model             = model,
    patterns          = patterns,
    patterns.total    = patterns.total,
    prob.item         = prob.item,
    prob.states       = prob.states,
    prob              = prob,
    N.person          = N.person,
    N.block           = N.block,
    I.block           = I.block,
    D                 = D,
    I.states          = I.states,
    fc.type           = fc.type,
    call              = call,
    arguments         = list(
      N.person = N.person,
      N.block  = N.block,
      I.block  = I.block,
      D        = D,
      model    = model,
      fc.type  = fc.type,
      control  = control
    )
  )

  class(data.obj) <- "data.FCGDINA"

  return(data.obj)
}


#' Compute FCGDINA forced-choice pattern probabilities
#'
#' @description
#' Computes block-level forced-choice response probabilities from known
#' attribute profiles and item-level CDM delta parameters. This is the
#' deterministic probability engine used by \code{\link{sim.data.FCGDINA}} for
#' response generation and is useful for checking a fitted or simulated
#' FCGDINA parameterization.
#'
#' @details
#' Item endorsement probabilities are first computed by
#' \code{\link{compute_fcgdina_item_prob}}. For person \eqn{n} and statement
#' \eqn{i}, the probability is
#' \deqn{
#'   p_{ni} = \mathbf{x}_i(\boldsymbol{\alpha}_n)' \boldsymbol{\delta}_i,
#' }
#' where \eqn{\mathbf{x}_i(\boldsymbol{\alpha}_n)} is the relevant row of the
#' item CDM design matrix. The item probabilities are then converted to
#' block-level forced-choice pattern probabilities by
#' \code{\link{forced_choice_from_agree}}.
#'
#' @param alpha \eqn{N \times D} binary matrix of attribute profiles.
#' @param delta.list List of per-item delta parameter vectors. The length must
#'   equal \code{nrow(Q.matrix)}.
#' @param design.matrix.list List of per-item CDM design matrices generated by
#'   \code{\link{get_design_matrix_cdm}}.
#' @param Q.matrix \eqn{I \times D} binary Q-matrix.
#' @param patterns.total List of full-ranking pattern matrices.
#' @param patterns List of observed pattern matrices.
#' @param prob.item Optional precomputed item endorsement-probability matrix.
#'   Supplying it avoids recomputation in simulation workflows.
#'
#' @return An \eqn{N \times \sum_b P_b} matrix of forced-choice response
#'   probabilities, where \eqn{P_b} is the number of observable patterns in
#'   block \eqn{b}.
#'
#' @seealso \code{\link{sim.data.FCGDINA}},
#'   \code{\link{compute_fcgdina_item_prob}},
#'   \code{\link{forced_choice_from_agree}}
#'
#' @examples
#' sim <- sim.data.FCGDINA(N.person = 10, N.block = 2, I.block = 2,
#'                         D = 2, model = "DINA")
#' prob <- model.FCGDINA(
#'   alpha = sim$alpha,
#'   delta.list = sim$delta.list,
#'   design.matrix.list = sim$design.matrix.list,
#'   Q.matrix = sim$Q.matrix,
#'   patterns.total = sim$patterns.total,
#'   patterns = sim$patterns
#' )
#' stopifnot(all(abs(rowSums(prob) - length(sim$block.items)) < 1e-8))
#'
#' @export
#' @keywords internal
model.FCGDINA <- function(alpha, delta.list, design.matrix.list,
                          Q.matrix, patterns.total, patterns,
                          prob.item = NULL) {
  if (is.null(prob.item)) {
    prob.item <- compute_fcgdina_item_prob(
      alpha = alpha, delta.list = delta.list,
      design.matrix.list = design.matrix.list,
      Q.matrix = Q.matrix
    )
  }
  forced_choice_from_agree(prob.item, patterns.total, patterns)
}

#' Compute FCGDINA item-level endorsement probabilities
#'
#' @description
#' Computes statement endorsement probabilities for supplied attribute
#' profiles under an FCGDINA item-level CDM parameterization.
#'
#' @param alpha \eqn{N \times D} binary matrix of attribute profiles.
#' @param delta.list List of per-item delta parameter vectors.
#' @param design.matrix.list List of per-item design matrices.
#' @param Q.matrix \eqn{I \times D} binary Q-matrix.
#' @param gs Optional \eqn{I \times 2} matrix of lower and upper probability
#'   bounds. If supplied, computed probabilities are clamped to
#'   \code{gs[i, 1]} and \code{gs[i, 2]} for item \eqn{i}.
#' @param eps Numeric probability clamp used when \code{gs = NULL}. Default
#'   \code{1e-12}.
#' @param link Item probability link. \code{"identity"} uses
#'   \eqn{p = X\delta}; \code{"logit"} uses
#'   \eqn{p = \mathrm{logit}^{-1}(X\delta)}. The latter is used by the Stan
#'   FCGDINA backend.
#'
#' @return An \eqn{N \times I} matrix of item endorsement probabilities.
#'
#' @details
#' The function maps each full attribute pattern in \code{alpha} to the
#' reduced pattern required by each item-specific Q-vector, evaluates the
#' corresponding design row, and multiplies by that item's delta vector.
#' This function does not sample responses.
#'
#' @seealso \code{\link{model.FCGDINA}}, \code{\link{sim.data.FCGDINA}}
#'
#' @keywords internal
compute_fcgdina_item_prob <- function(alpha, delta.list, design.matrix.list,
                                       Q.matrix, gs = NULL, eps = 1e-12,
                                       link = c("identity", "logit")) {
  link <- match.arg(link)
  alpha <- as.matrix(alpha)
  N <- nrow(alpha)
  D <- ncol(alpha)
  I <- nrow(Q.matrix)

  patterns.full <- attributepattern(D)
  pattern.names <- pattern_key(patterns.full)
  alpha.names <- pattern_key(alpha)
  alpha.idx <- match(alpha.names, pattern.names)

  prob.item <- matrix(NA_real_, N, I)

  for (i in seq_len(I)) {
    att <- which(Q.matrix[i, ] > 0)
    X <- design.matrix.list[[i]]
    delta <- delta.list[[i]]

    if (length(att) > 0) {
      reduced.key <- pattern_key(patterns.full[, att, drop = FALSE])
      reduced.idx <- match(reduced.key[alpha.idx],
                           pattern_key(attributepattern(length(att))))
    } else {
      reduced.idx <- rep(1L, N)
    }

    eta.i <- as.vector(X %*% delta)
    p.i <- if (link == "logit") stats::plogis(eta.i) else eta.i

    if (!is.null(gs)) {
      p.i <- pmin(pmax(p.i, gs[i, 1L]), gs[i, 2L])
    } else {
      p.i <- pmin(pmax(p.i, eps), 1 - eps)
    }

    prob.item[, i] <- p.i[reduced.idx]
  }

  prob.item
}