R/sim.data.FCMIRT.R

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

#' Simulate FCMIRT Data
#'
#' @description
#' Simulates forced-choice ranking data using a MIRT item response model
#' (\code{"m1pl"} through \code{"m4pl"}) as the item-level endorsement model.
#' Forced-choice responses are generated by applying the sequential
#' Luce/Plackett choice rule to the logits of the item-level endorsement
#' probabilities, with RANK, MOLE, and PICK patterns handled by
#' \code{\link{model.FCMIRT}}.
#' For identification, item difficulty/intercept parameters \code{b} are
#' centered within each forced-choice block so that their block sum is zero.
#'
#' @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.
#' @param D Number of latent dimensions. If \code{control$Q.matrix} is supplied,
#'   \code{D} is reset to \code{ncol(control$Q.matrix)}.
#' @param model MIRT model type: \code{"m1pl"}, \code{"m2pl"}, \code{"m3pl"}, or
#'   \code{"m4pl"}.
#' @param fc.type Forced-choice response type: \code{"RANK"}, \code{"MOLE"}, or
#'   \code{"PICK"}. A scalar value is recycled to all blocks.
#' @param control Optional list. Supported entries are \code{Q.matrix},
#'   \code{block.items}, and \code{Corr}.
#'
#' @return An object of class \code{"data.FCMIRT"}, a list containing
#'   \code{data}, \code{response}, \code{theta}, \code{par}, \code{Q.matrix},
#'   \code{block.items}, \code{Corr}, \code{model}, \code{patterns},
#'   \code{patterns.total}, \code{probability}, \code{prob}, and
#'   data-generating arguments.
#'
#' @seealso \code{\link{fit.FCMIRT}}, \code{\link{model.FCMIRT}}
#'
#' @examples
#' set.seed(123)
#' sim <- sim.data.FCMIRT(N.person = 20, N.block = 3, I.block = 2,
#'                        D = 2, model = "m2pl", fc.type = "RANK")
#' str(sim$data)
#' head(sim$response)
#' dim(sim$prob)
#'
#' @export
sim.data.FCMIRT <- function(N.person = 1000, N.block = 10, I.block = 2,
                            D = 3, model = "m2pl", fc.type = "RANK",
                            control = NULL) {

  call <- match.call()

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

  # ---- Model validation ----
  model <- tolower(model)
  model.valid <- c("m1pl", "m2pl", "m3pl", "m4pl")
  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)

  # ---- Corr and theta ----
  Corr  <- validate_corr_matrix(control$Corr, D)
  theta <- generate_theta_mvn(N.person, D, Corr)

  # ---- Auto-generate Q-matrix if not supplied ----
  if (is.null(Q.matrix)) {
    Q.matrix <- sim.Q.matrix.FC(I.states = I.states, D = D,
                                block.items = block.items,
                                allow.negative = FALSE)
  }

  # ---- MIRT item parameters ----
  if (model == "m1pl") {
    a <- Q.matrix
  } else {
    a <- matrix(rlnorm(I.states * D, meanlog = 0.25, sdlog = 0.25),
                nrow = I.states, ncol = D) * Q.matrix
  }

  b <- rnorm(I.states, 0, 1)
  for (items in block.items) {
    b[items] <- b[items] - mean(b[items])
  }
  c <- runif(I.states, 0.00, 0.35)
  d <- runif(I.states, 0.65, 1.00)

  if (model %in% c("m1pl", "m2pl")) {
    c <- rep(0.0, I.states)
    d <- rep(1.0, I.states)
  } else if (model == "m3pl") {
    d <- rep(1.0, I.states)
  }

  par <- cbind(a, b, c, d)

  # ---- Item-level endorsement probabilities ----
  probability <- model.MIRT(theta = theta, par = par)

  # ---- 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.FCMIRT(theta = theta, par = par,
                       patterns.total = patterns.total,
                       patterns = patterns)

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

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

  rownames(par) <- paste0("item", seq_len(I.states))
  colnames(par) <- c(paste0("a", seq_len(D)), "b", "c", "d")

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

  rownames(Corr) <- colnames(Corr) <- paste0("dim", seq_len(D))

  # ---- Output assembly ----
  data.obj <- list(
    data           = data,
    response       = response,
    theta          = theta,
    par            = par,
    Q.matrix       = Q.matrix,
    block.items    = block.items,
    Corr           = Corr,
    model          = model,
    patterns       = patterns,
    patterns.total = patterns.total,
    probability    = probability,
    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.FCMIRT"

  return(data.obj)
}


#' Compute FCMIRT Forced-Choice Pattern Probabilities
#'
#' @param theta N x D matrix of latent trait values.
#' @param par I x (D + 3) MIRT item parameter matrix, where b is the
#'   item difficulty and the item-level linear predictor is a' theta - b.
#' @param patterns.total List of full-ranking pattern matrices.
#' @param patterns List of observed pattern matrices.
#'
#' @return N x sum(pattern counts) probability matrix. For each block,
#'   item endorsement probabilities are first computed by \code{model.MIRT},
#'   transformed to logits, and then combined with a sequential
#'   Luce/Plackett rule. MOLE and PICK probabilities are obtained by summing
#'   compatible full-ranking probabilities and normalizing over the observed
#'   block patterns.
#' @keywords internal
model.FCMIRT <- function(theta, par, patterns.total, patterns) {
  cpp_model_FCMIRT(theta, par, patterns.total, patterns)
}