R/fit.TIRT.R

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

#' Fit the Thurstonian Item Response Theory (TIRT) Model
#'
#' @description
#' Fits the Thurstonian IRT model for forced-choice data. Unlike the
#' FCMIRT/FCGGUM models that operate on the item-level endorsement
#' probabilities, the TIRT model works directly at the level of paired
#' comparisons: each pairwise comparison is treated as a binary outcome
#' governed by a latent difference process with a probit link.
#'
#' @section Model Specification:
#'
#' \strong{Pairwise utility difference.}
#' For a pair of statements \eqn{(i, k)} in block \eqn{b}, the latent
#' comparative judgment is:
#' \deqn{
#'   y_{j,ik}^* =
#'   -\gamma_{ik} +
#'   \lambda_i \sum_{d=1}^{D} q_{id}\,\theta_{jd} -
#'   \lambda_k \sum_{d=1}^{D} q_{kd}\,\theta_{jd} +
#'   \varepsilon_{j,i} - \varepsilon_{j,k},
#' }
#' where:
#' \itemize{
#'   \item \eqn{\gamma_{ik}} is the pairwise intercept (unfolding threshold
#'         parameter). Skew-symmetry requires \eqn{\gamma_{ik} = -\gamma_{ki}}.
#'   \item \eqn{\lambda_i > 0} is the loading magnitude of statement
#'         \eqn{i}; the sign of the statement direction is carried by
#'         \eqn{q_{id}}.
#'   \item \eqn{q_{id} \in \{-1, 0, 1\}} is the Q-matrix entry; each row
#'         has exactly one non-zero entry.
#'   \item \eqn{\varepsilon_{j,i} \sim N(0, \psi_i^2)} is the uniqueness
#'         (error) of statement \eqn{i}.
#' }
#'
#' \strong{Binary pairwise response.} The observed response is
#' \deqn{
#'   Y_{j,ik} = \begin{cases}
#'     1, & \text{if } y_{j,ik}^* \ge 0 \quad (\text{prefer } i \text{ over } k), \\
#'     0, & \text{otherwise}.
#'   \end{cases}
#' }
#'
#' Under the normality of errors, the probability is:
#' \deqn{
#'   P(Y_{j,ik} = 1 \mid \boldsymbol{\theta}_j) =
#'   \Phi\!\left(
#'     \frac{
#'       -\gamma_{ik} + \lambda_i\,\mathbf{q}_i'\boldsymbol{\theta}_j -
#'       \lambda_k\,\mathbf{q}_k'\boldsymbol{\theta}_j
#'     }{
#'       \sqrt{\psi_i^2 + \psi_k^2}
#'     }
#'   \right),
#' }
#' where \eqn{\Phi(\cdot)} is the standard normal CDF.
#'
#' \strong{Identification constraints.} Several standard constraints are
#' applied:
#' \itemize{
#'   \item \strong{Case A} (\eqn{I_{\text{block}} = 2, D > 2}): All
#'         \eqn{\psi_i^2 = 0.5} are fixed.
#'   \item \strong{Case B} (\eqn{D = 2, I_{\text{block}} = 2}): The first
#'         block's \eqn{\lambda_i} are fixed to 0.80, all
#'         \eqn{\psi_i^2 = 0.5} are fixed.
#'   \item \strong{General case}: The last item in each block has
#'         \eqn{\psi_i^2 = 1.0} fixed. Pairwise \eqn{\gamma_{ik}}
#'         parameters for comparisons that are never observed under the
#'         supplied RANK/MOLE/PICK design are fixed to \code{gamma.mu}.
#' }
#'
#' \strong{Correlation structure.} The inter-trait correlation matrix
#' \eqn{\boldsymbol{\Sigma}} is estimated as in the MIRT model. In the Stan
#' file, \eqn{\boldsymbol{\theta}_j} is centered at zero; the iStEM backend
#' uses \code{control.model$theta.mu} as the mean in the normal kernel.
#'
#' @section Estimation Methods:
#'
#' \describe{
#'   \item{\strong{Stan} (\code{method = "stan"}):}{
#'     Full Bayesian inference via HMC. The probit likelihood is evaluated
#'     over all pairwise comparisons.}
#'   \item{\strong{iStEM} (\code{method = "iStEM"}):}{
#'     Improved Stochastic EM. Person sampling uses finite-grid block Gibbs;
#'     structural parameters (\eqn{\lambda}, \eqn{\psi^2}, \eqn{\gamma})
#'     updated via optimization or closed-form steps.}
#' }
#'
#' @param data An \eqn{N \times B} forced-choice data matrix with ranking
#'   strings (e.g., \code{"2>1>3"}).
#' @param Q.matrix An \eqn{I \times D} matrix with entries \eqn{-1, 0, 1}.
#'   Each row must have exactly one non-zero entry that assigns the statement
#'   to a specific dimension with a sign indicating the scoring direction.
#' @param block.items A list of length \eqn{B} giving item indices in each
#'   block. Parsed from \code{data} if \code{NULL}.
#' @param fc.type Forced-choice response type: \code{"RANK"} (default),
#'   \code{"MOLE"}, or \code{"PICK"}.
#' @param method Estimation method: \code{"iStEM"} (default) or
#'   \code{"stan"}.
#' @param control.model A named list of model-level hyperparameters for the
#'   Thurstonian IRT model. Supported entries:
#'   \describe{
#'     \item{\code{lambda.alpha}, \code{lambda.beta}}{Lower and upper bounds
#'           for free loading magnitudes \eqn{\lambda_i}. Defaults: 0.40,
#'           0.90. Stan uses
#'           \eqn{\lambda_i \sim U(\lambda_{\min}, \lambda_{\max})}, where
#'           these two controls supply \eqn{\lambda_{\min}} and
#'           \eqn{\lambda_{\max}};
#'           iStEM uses the same values as optimization bounds unless
#'           \code{lambda.lower} or \code{lambda.upper} is supplied in
#'           \code{control.method}.}
#'     \item{\code{psi.mu}, \code{psi.sigma}}{Prior mean and SD for
#'           uniqueness standard deviations \eqn{\psi_i} (normal, truncated
#'           to positive values). Defaults: 0.80, 0.30. The prior
#'           \eqn{\psi_i \sim N_+(0.80, 0.30^2)} reflects the expectation
#'           that uniqueness SDs are moderate in size. Note that the model
#'           parameterizes the uniqueness variance \eqn{\psi_i^2}, so
#'           the prior on the SD scale is transformed accordingly.}
#'     \item{\code{gamma.mu}, \code{gamma.sigma}}{Prior mean and SD for
#'           pairwise intercepts \eqn{\gamma_{ik}} (normal). Defaults: 0,
#'           0.8. The skew-symmetry constraint \eqn{\gamma_{ik} = -\gamma_{ki}}
#'           is automatically enforced. Unobserved pairwise comparisons
#'           (in MOLE/PICK designs) have their \eqn{\gamma} parameters fixed
#'           to \code{gamma.mu}.}
#'     \item{\code{theta.mu}}{Mean vector for the iStEM normal kernel for
#'           \eqn{\boldsymbol{\theta}_j}. Default: \code{rep(0, D)}. The
#'           current Stan model accepts this value in its data list but
#'           centers \eqn{\boldsymbol{\theta}_j} at zero.}
#'     \item{\code{L}}{Theta grid size per dimension for marginal
#'           log-likelihood computation and iStEM block Gibbs sampling.
#'           Default adapts to \eqn{D}.}
#'     \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). The TIRT model's pairwise probit
#'           likelihood typically mixes well; however, for high-dimensional
#'           models (\eqn{D > 5}) or models with many fixed parameters,
#'           consider increasing to 0.9.}
#'     \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 for
#'           burn-in (default: 1.5).}
#'     \item{\code{eps2}}{Monte Carlo error tolerance for the final
#'           chain (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
#'           parameter block during structural parameter updates
#'           (default: 50). Applies to the optimization of \eqn{\lambda},
#'           \eqn{\psi^2}, and \eqn{\gamma} parameters.}
#'     \item{\code{fix.corr}}{Logical; if \code{TRUE}, the inter-trait
#'           correlation matrix is fixed to the identity (default:
#'           \code{FALSE}). Useful when orthogonal dimensions are
#'           theoretically expected.}
#'     \item{\code{estimate.se}}{Logical; if \code{TRUE} (default),
#'           standard errors are computed from the final Monte Carlo
#'           chain.}
#'     \item{\code{lambda.lower}, \code{lambda.upper}}{Bounds on free
#'           loading magnitudes during iStEM optimization. Defaults are
#'           \code{lambda.alpha} and \code{lambda.beta} from
#'           \code{control.model}.}
#'     \item{\code{psi.lower}, \code{psi.upper}}{Bounds on uniqueness
#'           standard deviations \eqn{\psi_i} during iStEM optimization
#'           (defaults: 1e-4 and 5). The reported \code{par} matrix stores
#'           \eqn{\psi_i^2}.}
#'     \item{\code{gamma.lower}, \code{gamma.upper}}{Bounds on pairwise
#'           intercepts \eqn{\gamma_{ik}} during iStEM optimization
#'           (defaults: -6 and 6).}
#'   }
#'
#' @return An object of class \code{"TIRT"} containing:
#' \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 + 1)}). Columns include the
#'         loading \eqn{\lambda_i} on the design dimension and
#'         \eqn{\psi_i^2}.}
#'   \item{\code{Corr}}{List with \code{est}, \code{se}, \code{Rhat}
#'         (\eqn{D \times D}).}
#'   \item{\code{gamma.matrix}}{List with \code{est}, \code{se}, \code{Rhat}
#'         (\eqn{I \times I}); pairwise intercept matrix.}
#'   \item{\code{pairs.matrix}, \code{pairs.value}, \code{response},
#'         \code{response.total}}{Pairwise and response data structures.}
#'   \item{\code{logLik}}{Marginal log-likelihood (class \code{"logLik"}).}
#' }
#'
#' @references
#' Brown, A., & Maydeu-Olivares, A. (2011). Item response modeling of
#'   forced-choice questionnaires. \emph{Educational and Psychological
#'   Measurement}, 71(3), 460--502. \doi{10.1177/0013164410375112}
#'
#' Maydeu-Olivares, A., & Brown, A. (2010). Item response modeling of
#'   paired comparison and ranking data. \emph{Multivariate Behavioral
#'   Research}, 45(6), 935--974.
#'
#' Thurstone, L. L. (1927). A law of comparative judgment.
#'   \emph{Psychological Review}, 34(4), 273--286.
#'
#' @seealso
#' \code{\link{sim.data.TIRT}}, \code{\link{get.fit.index.TIRT}},
#' \code{\link{logLik.TIRT}}
#'
#' @examples
#' sim <- sim.data.TIRT(N.person = 20, N.block = 3, I.block = 2,
#'                      D = 2, fc.type = "RANK")
#' fit <- fit.TIRT(sim$data, Q.matrix = sim$Q.matrix,
#'                 block.items = sim$block.items,
#'                 fc.type = "RANK", 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)
#' cor(fit$theta$est, sim$theta)
#' gof <- get.fit.index(fit)
#' summary(gof)
#'
#' @export
fit.TIRT <- function(data, Q.matrix, block.items = NULL,
                     fc.type = NULL,
                     method = c("iStEM", "stan"),
                     control.model = NULL,
                     control.method = NULL) {

  call <- match.call()
  method <- match.arg(method)

  block.items <- resolve_block_items(block.items, data)
  fc.type     <- resolve_fc_type(fc.type, data, block.items)

  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.TIRT.stan(
      data = data,
      Q.matrix = Q.matrix,
      block.items = block.items,
      fc.type = fc.type,
      control.model = control.model,
      control.method = control.method,
      .call = call
    ))
  }

  fit.TIRT.iStEM(
    data = data,
    Q.matrix = Q.matrix,
    block.items = block.items,
    fc.type = fc.type,
    control.model = control.model,
    control.method = control.method,
    .call = call
  )
}