fit.MIRT: Fit the Multidimensional Item Response Theory (MIRT) Model

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

Fit the Multidimensional Item Response Theory (MIRT) Model

Description

Fits a multidimensional extension of the 1PL, 2PL, 3PL, or 4PL item response model to binary response data. Two estimation backends are provided: full Bayesian inference via Hamiltonian Monte Carlo (Stan) and a fast stochastic-EM algorithm (iStEM) that scales to large data sets.

Usage

fit.MIRT(
  data,
  model = "2PL",
  D = NULL,
  Q.matrix = NULL,
  method = c("iStEM", "stan"),
  control.model = NULL,
  control.method = NULL
)

Arguments

data	An N \times I matrix of binary responses coded as 0 (incorrect) and 1 (correct). Rows index persons, columns index items. Missing values are not allowed.
model	Character string specifying the model type. Accepts both internal codes ("m1pl", "m2pl", "m3pl", "m4pl") and user-friendly aliases ("Rasch", "1PL", "2PL", "3PL", "4PL"). Case-insensitive. Default is "2PL".
D	Integer; number of latent dimensions (D \ge 1). If NULL (default), D is inferred from Q.matrix; if both are NULL, an error is raised. When D = 1, the model reduces to unidimensional IRT.
Q.matrix	An optional I \times D binary matrix. Entry q_{id} = 1 indicates that the d-th dimension loads on item i. For 2PL–4PL models, q_{id} = 1 frees a_{id} and q_{id} = 0 fixes a_{id} = 0. For the current M1PL backend, slopes are fixed to 1 and the Q-matrix does not zero out slope entries. If NULL (default), a triangular identification pattern is used: all items load on all dimensions except the last D items, which follow a lower-triangular structure for rotational invariance resolution.
method	Estimation method: "iStEM" (fast stochastic EM for large data; default) or "stan" (full Bayesian HMC).
control.model	A named list of model-level hyperparameters. Supported entries: a.mu, a.sigma Prior location and scale for \log a_{id} (log-normal). Defaults: 0.25, 0.25. b.mu, b.sigma Prior mean and SD for b_i (normal). Defaults: 0, 1. c.mu, c.sigma Uniform support bounds [c_{\min}, c_{\max}] for c_i. Defaults: 0, 0.35. d.mu, d.sigma Uniform support bounds [d_{\min}, d_{\max}] for d_i. Defaults: 0.65, 1. theta.mu Prior mean vector for \boldsymbol{\theta}_j. Default: rep(0, D). L Theta grid size per dimension for marginal log-likelihood computation and iStEM block Gibbs sampling. Default adapts to D (e.g., 61 for D = 1, 31 for D = 2, 15 for D = 3). theta.lower, theta.upper Bounds for the theta grid used by marginal log-likelihood computation and iStEM block Gibbs sampling. Defaults: -6, 6. use.prior Logical (iStEM only). If FALSE, the item-prior penalty term is omitted from the item update, recovering an approximate maximum-likelihood item step. Default is TRUE.
control.method	A named list of method-specific tuning parameters. Common entries (used by both Stan and iStEM): cores Number of CPU cores for parallel chains (Stan) or ignored (iStEM). Default: the number of chains. vis Logical; if TRUE (default), prints progress information to the console. seed Random seed for reproducibility. Default: a random integer. Stan-specific entries: chains Number of MCMC chains (default: 2). iter Total iterations per chain (default: 5000). warmup Warmup/burn-in iterations per chain (default: iter / 2). thin Thinning interval (default: 1). init Initial values: "random" (default) for uniform(-2, 2) initialization, or a list of initial values per chain. algorithm MCMC algorithm: "HMC" (default), "HMC", or "Fixed_param". adapt_delta Target average acceptance probability (NUTS; default: 0.95). max_treedepth Maximum tree depth (NUTS; default: 10). stepsize Initial step size for the leapfrog integrator. If not set, Stan determines an appropriate value automatically during warmup. Only applicable when algorithm = "HMC" or "NUTS". int_time Total integration time for each HMC leapfrog trajectory. The number of steps is int_time / stepsize. Only applicable when algorithm = "HMC". metric The mass matrix for HMC sampling. Can be a unit vector ("unit_e"), a diagonal matrix ("diag_e"), or a dense matrix ("dense_e"). Stan defaults to "diag_e". adapt_engaged Logical; if TRUE (default), the warmup adaptation is enabled. Set to FALSE to disable step-size and mass-matrix adaptation during warmup. adapt_init_buffer Number of initial warmup iterations used for pure exploration before adaptation begins (default: 25). adapt_term_buffer Number of final warmup iterations where adaptation is frozen so the sampler can converge to the stationary distribution (default: 50). adapt_window Number of warmup iterations between adaptation updates. Larger values reduce the frequency of adaptation (default: 25). iStEM-specific entries: M Number of burn-in batches retained for convergence diagnosis (default: 10). Must be at least 2. Larger values improve the Geweke diagnostic stability but increase computation during burn-in. B Batch size: number of stochastic EM iterations per batch (default: 20). Each iteration samples all persons' \boldsymbol{\theta}_j blocks and updates all items. burnin.maxitr Maximum number of burn-in batches (default: 100). If convergence criteria are not met before this limit, the algorithm proceeds with a warning. Increase this value if convergence warnings appear consistently. maxitr Maximum number of total batches including post-burn-in iterations (default: 2000). The algorithm stops when either the MC error criterion is met or this limit is reached. eps1 Geweke z-score convergence threshold for the burn-in phase (default: 1.5). The burn-in phase ends when \sum z^2 / K < \epsilon_1 or the MC error criterion is satisfied. Lower values enforce stricter convergence but prolong burn-in. eps2 Monte Carlo error tolerance for the final chain (default: 0.4). The algorithm continues sampling until \max(d_k \cdot N) < \epsilon_2, where d_k is the batch-means variance for each parameter. Lower values yield more precise estimates. frac1, frac2 Fractions used in the Geweke convergence diagnostic (defaults: 0.1, 0.5). frac1 defines the proportion of the chain used for the early segment; frac2 defines the proportion for the late segment. These follow standard practice from the coda package. corr.optim.maxit Maximum L-BFGS-B iterations for the constrained unit-diagonal correlation update (default: 50). optim.maxit Maximum L-BFGS-B iterations per item during the item-parameter update step (default: 50). Increase if item-level optimization warnings appear. fix.corr Logical; if TRUE, the inter-trait correlation matrix is fixed to the identity during estimation (default: FALSE). Setting to TRUE enforces orthogonal latent dimensions. estimate.se Logical; if TRUE (default), standard errors are computed as the batch-means-based standard deviation of the final Monte Carlo chain. Set to FALSE to skip SE computation (slightly faster). a.lower, a.upper Bounds on each a_{id} (defaults: 1e-4, 6). For 2PL–4PL models, discrimination parameters are constrained to this interval during L-BFGS-B optimization. b.lower, b.upper Bounds on b_i (defaults: -6, 6). Item difficulty/intercept parameters are constrained to this interval during optimization. c.lower, c.upper Bounds on c_i for 3PL/4PL models. Defaults are inherited from the prior support bounds (control.model$c.mu = 0 and control.model$c.sigma = 0.35, respectively). Override these to impose tighter or wider bounds on the lower-asymptote parameters during optimization. d.lower, d.upper Bounds on d_i for 4PL models. Defaults are inherited from the prior support bounds (control.model$d.mu = 0.65 and control.model$d.sigma = 1, respectively). Override these to constrain the upper-asymptote parameters during optimization.

Value

An object of class "MIRT" containing the following components:

npar

Integer; number of free parameters (= freely estimated item parameters + free correlation elements).

method

Character; "stan" or "iStEM".

theta

List with matrices est, se, Rhat (N \times D); person parameter estimates.

par

List with matrices est, se, Rhat, free (I \times (D+3)); item parameter arrays. Columns are a1..aD, b, c, d.

Corr

List with matrices est, se, Rhat (D \times D); inter-trait correlation matrix.

Q.matrix

The I \times D Q-matrix used.

stan.obj

The stanfit object (Stan only; NULL for iStEM).

MCMC.obj

The list returned by rstan::extract() (Stan only; NULL for iStEM).

logLik

The marginal log-likelihood computed via Gauss–Hermite-type quadrature (class "logLik").

call

The matched call.

arguments

List of arguments used in fitting.

iStEM

List of iStEM diagnostic quantities (iStEM only), including burn-in size, convergence flags, and theta grid length.

Model Specification

Let Y_{ij} \in \{0, 1\} denote the binary response of person j = 1, \dots, N to item i = 1, \dots, I, and let \boldsymbol{\theta}_j = (\theta_{j1}, \dots, \theta_{jD})' denote the D-dimensional latent trait vector. The item response function (IRF) for the four model variants is:

M1PL (Rasch / one-parameter logistic):

P(Y_{ij} = 1 \mid \boldsymbol{\theta}_j) = \frac{1}{1 + \exp\bigl[-\bigl(\sum_{d=1}^{D} a_{id}\,\theta_{jd} - b_i\bigr)\bigr]}, \qquad a_{id} = 1

The M1PL implementation fixes the slope to 1 and is restricted to D = 1; fixed unit slopes do not identify separate dimensions.

M2PL (two-parameter logistic):

P(Y_{ij} = 1 \mid \boldsymbol{\theta}_j) = \frac{1}{1 + \exp\bigl[-\bigl(\sum_{d=1}^{D} a_{id}\,\theta_{jd} - b_i\bigr)\bigr]}, \qquad a_{id} > 0 \text{ if } q_{id} = 1,\; a_{id} = 0 \text{ otherwise}

The discrimination (slope) parameters a_{id} are freely estimated subject to the Q-matrix pattern and a log-normal prior.

M3PL (three-parameter logistic):

P(Y_{ij} = 1 \mid \boldsymbol{\theta}_j) = c_i + (1 - c_i) \times \frac{1}{1 + \exp\bigl[-\bigl(\sum_{d=1}^{D} a_{id}\,\theta_{jd} - b_i\bigr)\bigr]}

where c_i \in [0, 1) is the lower-asymptote (pseudo-guessing) parameter.

M4PL (four-parameter logistic):

P(Y_{ij} = 1 \mid \boldsymbol{\theta}_j) = c_i + (d_i - c_i) \times \frac{1}{1 + \exp\bigl[-\bigl(\sum_{d=1}^{D} a_{id}\,\theta_{jd} - b_i\bigr)\bigr]}

where c_i \in [0, 1) is the lower asymptote and d_i \in (0, 1] is the upper asymptote (1 - slippage).

In the unidimensional case (D = 1), these reduce to the standard 1PL–4PL models. For D > 1, the Q-matrix governs which dimensions load on each item, enabling both exploratory (default triangular identification) and confirmatory (user-specified Q-matrix) structures.

Prior distributions (Bayesian / MAP):

a_{id} (free)

Log-normal: \log a_{id} \sim N(\mu_a, \sigma_a^2). Defaults: \mu_a = 0.25, \sigma_a = 0.25.

b_i

Normal: b_i \sim N(\mu_b, \sigma_b^2). Defaults: \mu_b = 0, \sigma_b = 1.

c_i

Uniform: c_i \sim U(c_{\min}, c_{\max}). Defaults: c_{\min} = 0, c_{\max} = 0.35.

d_i

Uniform: d_i \sim U(d_{\min}, d_{\max}). Defaults: d_{\min} = 0.65, d_{\max} = 1.

\boldsymbol{\theta}_j

Multivariate normal: \boldsymbol{\theta}_j \sim N_D(\boldsymbol{\mu}_\theta, \boldsymbol{\Sigma}), with \boldsymbol{\mu}_\theta = \mathbf{0} and \boldsymbol{\Sigma} a correlation matrix.

Estimation Methods

Stan (method = "stan"):

Full Bayesian inference via Hamiltonian Monte Carlo (NUTS/HMC). The joint posterior of all parameters is explored using multiple Markov chains. Output includes posterior means, standard deviations, and \hat{R} convergence diagnostics for all parameters. Marginal log-likelihood is approximated via Gauss–Hermite-type quadrature over the latent space.

iStEM (method = "iStEM"):

An improved Stochastic EM (iStEM) algorithm that alternates between (a) sampling \boldsymbol{\theta}_j from its finite-grid full conditional as one D-dimensional block using the person's complete response likelihood and (b) maximizing the complete-data posterior for item parameters via L-BFGS-B. The inter-trait correlation matrix \boldsymbol{\Sigma} is estimated by constrained unit-diagonal normal-likelihood optimization. Convergence is monitored using Geweke (1992) convergence diagnostics and batch-means Monte Carlo error estimates. Standard errors are obtained from the final Monte Carlo chain.

References

Reckase, M. D. (2009). Multidimensional Item Response Theory. Springer. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/978-0-387-89976-3")}

Birnbaum, A. (1968). Some latent trait models and their use in inferring an examinee's ability. In F. M. Lord & M. R. Novick, Statistical theories of mental test scores (pp. 397–479). Addison-Wesley.

Barton, M. A., & Lord, F. M. (1981). An upper asymptote for the three-parameter logistic item-response model. ETS Research Report Series, 1981(1), i–21.

Geweke, J. (1992). Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. In J. M. Bernardo et al. (Eds.), Bayesian Statistics 4 (pp. 169–193). Oxford University Press.

See Also

good.of.fit for goodness-of-fit evaluation, get.fit.index.MIRT for MIRT-specific fit indices, sim.data.MIRT for simulating data from this model, rotate for post-hoc rotation of MIRT solutions, logLik.MIRT for marginal log-likelihood extraction.

Examples

# Simulate data from a 2-dimensional 2PL model
sim <- sim.data.MIRT(N = 20, I = 6, D = 2, model = "m2pl")

# Fit via iStEM (fast, large-data friendly)
fit_istem <- fit.MIRT(sim$response, model = "m2pl", 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))


# stan code, long time
# Fit via Stan (full Bayesian inference, computationally heavier)
fit_stan <- fit.MIRT(sim$response, model = "m2pl", D = 2,
                     method = "stan",
                     control.method = list(
                       chains = 1, iter = 200, warmup = 100,
                       cores = 1, seed = 123))


# Extract results
print(fit_istem$par$est)       # item parameter estimates
print(fit_istem$theta$est)     # person trait estimates
print(fit_istem$Corr$est)      # factor correlation matrix

# Compute goodness-of-fit indices
gof <- get.fit.index(fit_istem)
summary(gof)

ForceChoice documentation built on Sept. 13, 2026, 1:06 a.m.