R/fit.MGPCM.R

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

#' Fit the Multidimensional Generalized Partial Credit Model (MGPCM)
#'
#' @description
#' Fits a multidimensional generalization of the partial credit model to
#' polytomous (multi-category) response data. Two estimation backends are
#' provided: full Bayesian inference via Hamiltonian Monte Carlo (Stan) and
#' a fast stochastic-EM algorithm (iStEM).
#'
#' @section Model Specification:
#'
#' Let \eqn{Y_{ij} \in \{0, 1, \dots, K_i - 1\}} denote the categorical
#' response of person \eqn{j = 1, \dots, N} to item \eqn{i = 1, \dots, I},
#' where \eqn{K_i \ge 2} is the number of response categories for item
#' \eqn{i}. Let \eqn{\boldsymbol{\theta}_j} be the \eqn{D}-dimensional
#' latent trait vector.
#'
#' The category response probability implemented in both the C++ and Stan
#' backends is a softmax over category scores:
#' \deqn{
#'   P(Y_{ij} = k \mid \boldsymbol{\theta}_j) =
#'   \frac{
#'     \exp\{k\,\eta_{ij} + d_{ik}\}
#'   }{
#'     \sum_{r=0}^{K_i-1} \exp\{r\,\eta_{ij} + d_{ir}\}
#'   },
#'   \qquad k = 0, 1, \dots, K_i - 1,
#' }
#' where
#' \deqn{
#'   \eta_{ij} = \sum_{d=1}^{D} a_{id}\,\theta_{jd}.
#' }
#' The discrimination parameters satisfy \eqn{a_{id} > 0} when
#' \eqn{q_{id} = 1} and \eqn{a_{id} = 0} when \eqn{q_{id} = 0}. The
#' \eqn{d_{ik}} values are category intercepts, not cumulative step
#' difficulties; \eqn{d_{i0} = 0} is fixed for identification and
#' \eqn{d_{i1}, \dots, d_{i,K_i-1}} are free.
#'
#' This formulation nests the standard (unidimensional) generalized partial
#' credit model (Muraki, 1992) when \eqn{D = 1} and all \eqn{q_{i1} = 1}.
#'
#' \strong{Prior distributions:}
#' \describe{
#'   \item{\eqn{a_{id}} (free)}{Log-normal: \eqn{\log a_{id} \sim N(\mu_a, \sigma_a^2)}}
#'   \item{\eqn{d_{ik}} (free)}{Normal:
#'         \eqn{d_{ik} \sim N(\mu_d, \sigma_d^2)} for \eqn{k \ge 1};
#'         \eqn{d_{i0}=0} is fixed. The same normal penalty is used by
#'         Stan and by the iStEM item update when \code{use.prior = TRUE}.}
#'   \item{\eqn{\boldsymbol{\theta}_j}}{Multivariate normal:
#'         \eqn{\boldsymbol{\theta}_j \sim N_D(\mathbf{0}, \boldsymbol{\Sigma})}}
#' }
#'
#' @section Estimation Methods:
#'
#' \describe{
#'   \item{\strong{Stan} (\code{method = "stan"}):}{
#'     Full Bayesian inference via HMC. The joint posterior is sampled
#'     with multiple chains, providing posterior means, standard deviations,
#'     and \eqn{\hat{R}} diagnostics.}
#'   \item{\strong{iStEM} (\code{method = "iStEM"}):}{
#'     Improved Stochastic EM alternating finite-grid block Gibbs
#'     person-sampling and item-parameter optimization (L-BFGS-B).
#'     Convergence monitored via Geweke diagnostics and batch-means MC error.}
#' }
#'
#' @param data An \eqn{N \times I} matrix of integer responses coded
#'   \eqn{0, 1, \dots, K_i - 1}. Rows index persons, columns index items.
#' @param D Integer; number of latent dimensions (\eqn{D \ge 1}).
#'   Default is \code{2}.
#' @param Q.matrix An optional \eqn{I \times D} binary matrix. Entry
#'   \eqn{q_{id} = 1} frees \eqn{a_{id}}; \eqn{q_{id} = 0} fixes it to 0.
#'   Default uses a triangular identification structure.
#' @param length.poly Optional integer scalar or length-\eqn{I} vector giving
#'   the number of categories per item. If \code{NULL}, counts are inferred
#'   from observed maxima; supply this argument when a valid category is
#'   unobserved in the sample.
#' @param method Estimation method: \code{"iStEM"} (default) or
#'   \code{"stan"}.
#' @param control.model A named list of model-level hyperparameters.
#'   Supported entries:
#'   \describe{
#'     \item{\code{a.mu}, \code{a.sigma}}{Prior location and scale for
#'           \eqn{\log a_{id}} (log-normal). Defaults: 0.25, 0.25.
#'           Controls the prior mean and spread of the discrimination
#'           parameters across dimensions.}
#'     \item{\code{d.mu}, \code{d.sigma}}{Prior mean and SD for free
#'           category intercepts \eqn{d_{ik}} (\eqn{k \ge 1}; normal).
#'           Defaults: 0, 1. The first category intercept
#'           \eqn{d_{i0}=0} is fixed for identification.}
#'     \item{\code{theta.mu}}{Prior mean vector for
#'           \eqn{\boldsymbol{\theta}_j}. Default: \code{rep(0, D)}.
#'           A vector of length \eqn{D} specifying the prior mean for
#'           each latent dimension.}
#'     \item{\code{L}}{Theta grid size per dimension for marginal
#'           log-likelihood computation and iStEM block Gibbs sampling.
#'           Default adapts to \eqn{D}
#'           (e.g., 61 for D = 1, 31 for D = 2, 15 for D = 3).
#'           Larger grids increase numerical precision at the cost of
#'           exponential growth in computation (\eqn{L^D} total nodes).}
#'     \item{\code{theta.lower}, \code{theta.upper}}{Bounds for the
#'           theta grid used by marginal log-likelihood computation and
#'           iStEM block Gibbs sampling. Defaults: -6, 6.}
#'   }
#' @param control.method A named list of method-specific tuning parameters.
#'   Common entries (used by both Stan and iStEM):
#'   \describe{
#'     \item{\code{cores}}{Number of CPU cores for parallel chains
#'           (Stan) or ignored (iStEM). Default: the number of
#'           \code{chains}.}
#'     \item{\code{vis}}{Logical; if \code{TRUE} (default), prints progress
#'           information to the console.}
#'     \item{\code{seed}}{Random seed for reproducibility.
#'           Default: a random integer.}
#'   }
#'   Stan-specific entries:
#'   \describe{
#'     \item{\code{chains}}{Number of MCMC chains (default: 2).}
#'     \item{\code{iter}}{Total iterations per chain (default: 5000).}
#'     \item{\code{warmup}}{Warmup/burn-in iterations per chain
#'           (default: \code{iter / 2}).}
#'     \item{\code{thin}}{Thinning interval (default: 1).}
#'     \item{\code{init}}{Initial values: \code{"random"} (default)
#'           for uniform(-2, 2) initialization, or a list of initial
#'           values per chain.}
#'     \item{\code{algorithm}}{MCMC algorithm: \code{"HMC"} (default),
#'           \code{"HMC"}, or \code{"Fixed_param"}.}
#'     \item{\code{adapt_delta}}{Target average acceptance probability
#'           (NUTS; default: 0.95). Values closer to 1 reduce step size
#'           and improve sampling for difficult posteriors.}
#'     \item{\code{max_treedepth}}{Maximum tree depth (NUTS; default: 10).
#'           Increase if "max treedepth exceeded" warnings appear.}
#'     \item{\code{stepsize}}{Initial step size for the leapfrog
#'           integrator (auto-tuned by Stan if not set).}
#'     \item{\code{int_time}}{Total integration time for HMC trajectories
#'           (only when \code{algorithm = "HMC"}).}
#'     \item{\code{metric}}{Mass matrix type: \code{"unit_e"},
#'           \code{"diag_e"} (default), or \code{"dense_e"}.}
#'     \item{\code{adapt_engaged}}{Logical; if \code{TRUE} (default),
#'           warmup adaptation is enabled.}
#'     \item{\code{adapt_init_buffer}, \code{adapt_term_buffer},
#'           \code{adapt_window}}{Warmup adaptation scheduling parameters
#'           (defaults: 25, 50, 25).}
#'   }
#'   iStEM-specific entries:
#'   \describe{
#'     \item{\code{M}}{Number of burn-in batches retained for Geweke
#'           convergence diagnosis (default: 10; must be \eqn{\ge 2}).}
#'     \item{\code{B}}{Batch size: MCMC iterations per batch (default: 20).}
#'     \item{\code{burnin.maxitr}}{Maximum burn-in batches (default: 100).}
#'     \item{\code{maxitr}}{Maximum total batches (default: 2000).}
#'     \item{\code{eps1}}{Geweke z-score convergence threshold
#'           (default: 1.5).}
#'     \item{\code{eps2}}{Monte Carlo error tolerance (default: 0.4).}
#'     \item{\code{frac1}, \code{frac2}}{Fractions for the Geweke
#'           diagnostic (defaults: 0.1, 0.5).}
#'     \item{\code{corr.optim.maxit}}{Maximum L-BFGS-B iterations for the
#'           constrained unit-diagonal correlation update (default: 50).}
#'     \item{\code{optim.maxit}}{Maximum L-BFGS-B iterations per item
#'           (default: 50).}
#'     \item{\code{fix.corr}}{Logical; fix correlations to identity
#'           (default: \code{FALSE}).}
#'     \item{\code{estimate.se}}{Logical; compute standard errors from
#'           final MC chain (default: \code{TRUE}).}
#'     \item{\code{a.lower}, \code{a.upper}}{Bounds on \eqn{a_{id}}
#'           (defaults: 1e-4, 6).}
#'     \item{\code{d.lower}, \code{d.upper}}{Bounds on free category
#'           intercepts \eqn{d_{ik}} for \eqn{k \ge 1}. Defaults are -8
#'           and 8 in the iStEM backend.}
#'   }
#'
#' @return An object of class \code{"MGPCM"} with components:
#' \describe{
#'   \item{\code{npar}}{Number of free parameters.}
#'   \item{\code{method}}{\code{"stan"} or \code{"iStEM"}.}
#'   \item{\code{theta}}{List with \code{est}, \code{se}, \code{Rhat}
#'         (\eqn{N \times D}).}
#'   \item{\code{par}}{List with \code{est}, \code{se}, \code{Rhat},
#'         \code{free} (\eqn{I \times (D + K_{max})}). Columns are
#'         \code{a1..aD, d0, d1, ..., d_{K_{max}-1}}.}
#'   \item{\code{Corr}}{List with \code{est}, \code{se}, \code{Rhat}
#'         (\eqn{D \times D}).}
#'   \item{\code{length.poly}}{Integer vector of per-item category counts.}
#'   \item{\code{logLik}}{Marginal log-likelihood (class \code{"logLik"}).}
#'   \item{\code{call}, \code{arguments}}{Call and argument records.}
#' }
#'
#' @references
#' Muraki, E. (1992). A generalized partial credit model: Application of an
#'   EM algorithm. \emph{Applied Psychological Measurement}, 16(2), 159--176.
#'   \doi{10.1177/014662169201600206}
#'
#' Yao, L., & Schwarz, R. D. (2006). A multidimensional partial credit model
#'   for polytomous data. \emph{Applied Psychological Measurement}, 30(4),
#'   295--318.
#'
#' @seealso
#' \code{\link{sim.data.MGPCM}}, \code{\link{get.fit.index.MGPCM}},
#' \code{\link{logLik.MGPCM}}, \code{\link{rotate}}
#'
#' @examples
#' sim <- sim.data.MGPCM(N = 20, I = 6, D = 2, length.poly = 4)
#' fit <- fit.MGPCM(sim$response, D = 2, method = "iStEM",
#'                  control.method = list(
#'                    vis = FALSE, seed = 123,
#'                    M = 2, B = 2, burnin.maxitr = 2,
#'                    maxitr = 3, eps1 = 10, eps2 = 10,
#'                    estimate.se = FALSE))
#' head(fit$theta$est)
#' fit$par$est[1:5, ]
#' gof <- get.fit.index(fit)
#' summary(gof)
#'
#' @export
fit.MGPCM <- function(data, D = NULL, Q.matrix = NULL, length.poly = NULL,
                      method = c("iStEM", "stan"),
                      control.model = NULL,
                      control.method = NULL) {

  call <- match.call()
  method <- match.arg(method)
  data <- istem_prepare_response(data, binary = FALSE)
  length.poly <- istem_prepare_length_poly(data, length.poly)
  D <- resolve_D(D, Q.matrix = Q.matrix, I = ncol(data))
  Q.matrix <- istem_prepare_01_q(
    I = ncol(data), D = D, Q.matrix = Q.matrix, triangular = TRUE
  )
  control.model  <- fc_as_control_list(control.model, "control.model")
  control.method <- fc_as_control_list(control.method, "control.method")

  if (method == "stan") {
    return(fit.MGPCM.stan(
      response = data,
      D = D,
      Q.matrix = Q.matrix,
      length.poly = length.poly,
      control.model = control.model,
      control.method = control.method,
      .call = call
    ))
  }

  fit.MGPCM.iStEM(
    response = data,
    D = D,
    Q.matrix = Q.matrix,
    length.poly = length.poly,
    control.model = control.model,
    control.method = control.method,
    .call = call
  )
}