Introduction to irtQ
Overview
The irtQ package fits unidimensional item response
theory (IRT) models to test data that may contain both dichotomous and
polytomous items. The dichotomous models are the one-, two-, and
three-parameter logistic models (1PLM, 2PLM, and 3PLM), and the
polytomous models are the graded response model (GRM) and the
generalized partial credit model (GPCM). A single test can mix these
models.
The package covers the steps of a typical IRT analysis: estimating
item parameters by marginal maximum likelihood with the EM algorithm
(est_irt(), est_mg()), calibrating pretest
items with fixed item parameters or fixed abilities
(est_irt(), est_item()), estimating examinees’
abilities (est_score()), evaluating model-data fit
(irtfit(), sx2_fit()), and detecting
differential item functioning (rdif(),
crdif(), grdif(), catsib()). It
also provides functions for item and test information
(info()), item response curves (traceline()),
classification accuracy and consistency (cac_lee(),
cac_rud()), multistage test panels
(panel_info(), reval_mst(),
run_mst()), and classical test theory
(ctt()).
This vignette walks through a short workflow: read item parameters,
simulate responses, calibrate the items, estimate abilities, check fit,
and compute test information. Every example uses the scaling constant
D = 1.
Simulating item responses
simdat() generates item responses from item metadata and
a vector of abilities. Here 500 examinees are drawn from a standard
normal distribution.
set.seed(2026)
theta <- rnorm(500, mean = 0, sd = 1)
resp <- simdat(x = x, theta = theta, D = 1)
dim(resp)
#> [1] 500 25
Estimating item parameters
est_irt() calibrates the items. The model and the number
of score categories are given for all items at once. A Beta(5, 16) prior
on the guessing parameters, which is the default, stabilizes the 3PLM
estimates; it is written out here for clarity.
fit <- est_irt(
data = resp, D = 1, model = "3PLM", cats = 2, item.id = x$id,
use.gprior = TRUE, gprior = list(dist = "beta", params = c(5, 16)),
verbose = FALSE
)
fit
#>
#> Call:
#> est_irt(data = resp, D = 1, model = "3PLM", cats = 2, item.id = x$id,
#> use.gprior = TRUE, gprior = list(dist = "beta", params = c(5,
#> 16)), verbose = FALSE)
#>
#> Item parameter estimation using MMLE-EM.
#> 29 E-step cycles were completed using 49 quadrature points.
#> First-order test: Convergence criteria are satisfied.
#> Second-order test: Solution is a possible local maximum.
#> Computation of variance-covariance matrix:
#> Variance-covariance matrix of item parameter estimates is obtainable.
#>
#> Log-likelihood: -7122.816
Printing the fitted object reports the number of EM cycles, the
number of quadrature points, the first- and second-order convergence
checks, and the log-likelihood. A fuller estimation summary is available
from summary(fit).
getirt() extracts parts of the fitted object, such as
the item parameter estimates.
par_est <- getirt(fit, what = "par.est")
head(par_est)
#> id cats model par.1 par.2 par.3
#> 1 CMC1 2 3PLM 0.548348 1.52622329 0.21322043
#> 2 CMC2 2 3PLM 1.834283 -1.18009847 0.18392144
#> 3 CMC3 2 3PLM 1.194246 0.80370044 0.21057443
#> 4 CMC4 2 3PLM 1.104963 -0.40524060 0.25026086
#> 5 CMC5 2 3PLM 1.184828 -0.02633171 0.14063801
#> 6 CMC6 2 3PLM 1.580076 0.53583179 0.08719593
Estimating abilities
est_score() estimates the abilities of the examinees. A
fitted est_irt object can be passed directly, and the
response data stored in the object are used. The code below uses maximum
likelihood estimation (ML) and compares the estimates with the abilities
used to simulate the data. The range argument keeps the ML
estimates within [-4, 4], which matters for examinees who answer all
items correctly or all incorrectly.
scores <- est_score(fit, method = "ML", range = c(-4, 4))
head(scores$est.theta)
#> [1] 0.08603911 -1.30225167 0.08313340 0.08013606 -0.28721183 -4.00000000
cor(scores$est.theta, theta)
#> [1] 0.8777982
Model-data fit
sx2_fit() computes the S-X2 item fit statistic, and
irtfit() computes chi-square and likelihood-ratio
statistics, infit and outfit, and the contingency tables behind the
residual plots.
sx2 <- sx2_fit(fit)
head(sx2$fit_stat)
#> id chisq df crit.val p
#> 1 CMC1 11.342 16 26.296 0.788
#> 2 CMC2 6.096 11 19.675 0.867
#> 3 CMC3 19.788 16 26.296 0.230
#> 4 CMC4 15.647 15 24.996 0.406
#> 5 CMC5 19.984 15 24.996 0.173
#> 6 CMC6 13.013 15 24.996 0.601
irtfit() groups the examinees by the ML estimates from
the previous section, and the plot shows the raw and standardized
residuals of the first item.
fit_irt <- irtfit(
x = fit, score = scores$est.theta, group.method = "equal.width",
n.width = 10, loc.theta = "middle"
)
plot(fit_irt, item.loc = 1, type = "both", ci.method = "wald",
show.table = FALSE, ylim.sr.adjust = TRUE)
