1. Overview
Once item and person parameters have been estimated, a psychometric
study usually continues with a battery of follow-up analyses: evaluating
measurement precision, building score-conversion tables, checking person
and item fit, screening for local dependence and differential item
functioning, summarizing reliability, and sometimes probing for latent
subpopulations. This vignette walks through the post-estimation tools in
tirt. Every function accepts the item-parameter data
frame produced by binary_irt(),
polytomous_irt(), or mixed_irt(), so the
workflow is uniform across models.
We start from a small simulated 2PL data set and a quick calibration.
set.seed(2025)
sim <- sim_irt(
n_people = 600,
item_structure = list(list(model = "2PL", n_items = 12))
)
#> ----------------------------------------------------------------
#> Starting Simulation for N = 600 examinees...
#> >> Ability (Theta): Generated from N(mean=0.00, sd=1.00).
#> >> Block 1: 12 items using 2PL
#> - Discrimination (a): Default values used (Fixed at 1).
#> - Difficulty (b): Default values used (Fixed at 0).
#> ----------------------------------------------------------------
#> Constructing final data frames...
#> Simulation Complete.
#> Summary: 12 items, 600 examinees.
#> ----------------------------------------------------------------
fit <- binary_irt(sim$resp, model = "2PL", method = "EM",
control = list(max_iter = 20, verbose = FALSE))
head(fit$item_params)
#> item discrimination discrimination_se difficulty difficulty_se number
#> 1 item_1 1.010 0.108 -0.057 0.089 600
#> 2 item_2 0.946 0.105 -0.001 0.094 600
#> 3 item_3 0.974 0.106 0.073 0.092 600
#> 4 item_4 0.866 0.102 -0.099 0.102 600
#> 5 item_5 0.906 0.103 -0.070 0.098 600
#> 6 item_6 0.970 0.106 0.016 0.092 600
#> pvalue
#> 1 0.512
#> 2 0.500
#> 3 0.485
#> 4 0.518
#> 5 0.513
#> 6 0.4972. Item and Test Information
item_info() returns the Fisher information of each item
across a grid of ability values, and test_info() sums it
into the test information function together with the conditional
standard error of measurement (SEM).
theta_grid <- seq(-3, 3, by = 0.5)
# Item information (items in rows, theta in columns)
info <- item_info(fit$item_params, theta = theta_grid)
round(info[1:3, ], 3)
#> -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
#> item_1 0.047 0.074 0.110 0.156 0.205 0.243 0.255 0.236 0.194 0.145 0.101 0.067
#> item_2 0.047 0.070 0.102 0.140 0.180 0.212 0.224 0.212 0.180 0.140 0.102 0.070
#> item_3 0.043 0.066 0.098 0.139 0.183 0.220 0.237 0.227 0.195 0.151 0.109 0.075
#> 3
#> item_1 0.043
#> item_2 0.047
#> item_3 0.049
# Test information function and conditional SEM
tif <- test_info(fit$item_params, theta = theta_grid)
tif
#> theta test_info sem reliability
#> 1 -3.0 0.5479329 1.3509408 0.3539772
#> 2 -2.5 0.8379575 1.0924184 0.4559178
#> 3 -2.0 1.2398965 0.8980640 0.5535508
#> 4 -1.5 1.7458850 0.7568193 0.6358187
#> 5 -1.0 2.2902138 0.6607878 0.6960684
#> 6 -0.5 2.7337536 0.6048119 0.7321730
#> 7 0.0 2.9109843 0.5861113 0.7443099
#> 8 0.5 2.7414786 0.6039592 0.7327260
#> 9 1.0 2.2991850 0.6594973 0.6968948
#> 10 1.5 1.7505310 0.7558143 0.6364338
#> 11 2.0 1.2393432 0.8982645 0.5534405
#> 12 2.5 0.8341527 1.0949069 0.4547891
#> 13 3.0 0.5430402 1.3570131 0.3519287
# Where does the test measure most precisely?
tif$theta[which.max(tif$test_info)]
#> [1] 03. Summed-Score to Theta Conversion Table
Operational programs report scale scores from raw (summed) scores.
score_table() produces the conversion using the
Lord-Wingersky recursion, supporting expected a posteriori
("EAP"), weighted likelihood ("WLE"), and
maximum likelihood ("MLE") scoring.
# EAP conversion table (0 to 12 correct)
score_table(fit$item_params, method = "EAP")
#> summed_score theta se
#> 1 0 -1.845 0.643
#> 2 1 -1.468 0.605
#> 3 2 -1.131 0.572
#> 4 3 -0.826 0.548
#> 5 4 -0.541 0.532
#> 6 5 -0.269 0.522
#> 7 6 -0.004 0.519
#> 8 7 0.262 0.522
#> 9 8 0.534 0.531
#> 10 9 0.819 0.547
#> 11 10 1.124 0.571
#> 12 11 1.461 0.604
#> 13 12 1.840 0.643
# Maximum-likelihood conversion
score_table(fit$item_params, method = "MLE")
#> summed_score theta se
#> 1 0 -4.000 NA
#> 2 1 -2.471 1.096
#> 3 2 -1.654 0.810
#> 4 3 -1.129 0.679
#> 5 4 -0.713 0.625
#> 6 5 -0.348 0.595
#> 7 6 -0.004 0.588
#> 8 7 0.339 0.595
#> 9 8 0.703 0.624
#> 10 9 1.117 0.678
#> 11 10 1.639 0.783
#> 12 11 2.449 1.054
#> 13 12 4.000 NA4. Person Fit
person_fit() computes the lz standardized
log-likelihood index to identify examinees whose response patterns are
unlikely under the model.
pf <- person_fit(sim$resp, fit$item_params, fit$person_params)
head(pf)
#> person n_items theta loglik lz flag
#> 1 1 12 0.668 -7.585 0.102 FALSE
#> 2 2 12 0.325 -8.121 0.069 FALSE
#> 3 3 12 0.044 -8.410 -0.763 FALSE
#> 4 4 12 0.399 -8.259 -0.259 FALSE
#> 5 5 12 0.406 -8.186 -0.159 FALSE
#> 6 6 12 -1.591 -5.518 0.012 FALSE
# Number of examinees flagged as potentially misfitting
sum(pf$flag, na.rm = TRUE)
#> [1] 05. Item Fit
item_fit() returns the infit and outfit mean-square
statistics and their standardized versions, for dichotomous and
polytomous items alike.
item_fit(sim$resp, fit$item_params)
#> item n outfit outfit_z infit infit_z
#> 1 item_1 600 0.861 -3.893 0.892 -3.746
#> 2 item_2 600 0.875 -3.783 0.903 -3.518
#> 3 item_3 600 0.863 -3.989 0.901 -3.517
#> 4 item_4 600 0.894 -3.476 0.916 -3.257
#> 5 item_5 600 0.887 -3.550 0.909 -3.431
#> 6 item_6 600 0.871 -3.779 0.899 -3.621
#> 7 item_7 600 0.871 -3.849 0.902 -3.530
#> 8 item_8 600 0.848 -3.991 0.880 -3.998
#> 9 item_9 600 0.891 -3.482 0.911 -3.392
#> 10 item_10 600 0.834 -4.241 0.876 -4.036
#> 11 item_11 600 0.829 -4.278 0.875 -4.012
#> 12 item_12 600 0.869 -3.843 0.900 -3.5616. Local Dependence (Yen’s Q3)
ld_stats() computes Yen’s Q3 residual correlations for
every item pair. Large positive values signal that a pair depends on
something beyond the common trait, which is the situation a testlet
model is designed to absorb.
q3 <- ld_stats(sim$resp, fit$item_params)
round(q3[1:5, 1:5], 3)
#> item_1 item_2 item_3 item_4 item_5
#> item_1 1.000 -0.039 -0.004 -0.067 -0.040
#> item_2 -0.039 1.000 -0.098 0.001 -0.155
#> item_3 -0.004 -0.098 1.000 -0.056 -0.105
#> item_4 -0.067 0.001 -0.056 1.000 -0.018
#> item_5 -0.040 -0.155 -0.105 -0.018 1.000
# Largest absolute residual correlation
attr(q3, "max_abs_q3")
#> [1] 0.18412577. Differential Item Functioning
dif() screens dichotomous items for DIF using the
Mantel-Haenszel procedure (with the ETS delta effect size and A/B/C
flags) and logistic regression (which separates uniform and non-uniform
DIF). Here we plant DIF in item 3.
resp_dif <- sim$resp
grp <- rep(c("Reference", "Focal"), each = 300)
flip <- grp == "Focal" & resp_dif[[3]] == 1
resp_dif[[3]][flip] <- rbinom(sum(flip), 1, 0.55)
dif(resp_dif, group = grp)[, c("item", "MH_delta", "ETS_class", "LR_p")]
#> item MH_delta ETS_class LR_p
#> 1 item_1 0.175 A 0.6064
#> 2 item_2 -0.235 A 0.1189
#> 3 item_3 1.910 C 0.0000
#> 4 item_4 -0.587 A 0.2569
#> 5 item_5 0.599 A 0.3319
#> 6 item_6 0.155 A 0.3628
#> 7 item_7 0.434 A 0.4699
#> 8 item_8 -0.781 A 0.1236
#> 9 item_9 -0.820 A 0.2681
#> 10 item_10 -0.144 A 0.9380
#> 11 item_11 -0.262 A 0.0020
#> 12 item_12 -0.475 A 0.73138. Reliability
reliability() reports the empirical (marginal)
reliability from person estimates and their standard errors, the
model-based marginal reliability from the test information function, and
classical Cronbach’s alpha from the raw responses.
reliability(person_params = fit$person_params,
data = sim$resp,
item_params = fit$item_params)
#> Index Value
#> 1 empirical_reliability 0.7256
#> 2 marginal_reliability 0.6886
#> 3 cronbach_alpha 0.7070
#> 4 n_items 12.0000
#> 5 n_persons 600.00009. Test Characteristic Curve
tcc() returns the expected score of each item and the
test characteristic curve, which maps ability onto the number-correct
metric.
curves <- tcc(fit$item_params, theta = seq(-3, 3, by = 1))
curves$test_curve
#> theta expected_score
#> 1 -3 0.6204589
#> 2 -2 1.5004155
#> 3 -1 3.2870410
#> 4 0 6.0113887
#> 5 1 8.7381246
#> 6 2 10.5220784
#> 7 3 11.393812710. Multidimensional and Mixed Simulation
Two new simulators extend the data-generation toolkit.
sim_mirt() produces compensatory multidimensional data
(aligned with mirt_binary()), and sim_tirt()
produces forms that mix independent items with testlets (aligned with
irt_trt()).
# Two correlated dimensions, simple structure
mdat <- sim_mirt(
n_people = 400,
dimension = 2,
Sigma = matrix(c(1, 0.4, 0.4, 1), 2, 2),
item_structure = list(
list(model = "M2PL", n_items = 6, dims = 1),
list(model = "M2PL", n_items = 6, dims = 2)
)
)
#> ----------------------------------------------------------------
#> Starting Multidimensional Simulation (N = 400, D = 2)...
#> >> Ability (Theta): Generated from MVN with mean 0.00 and supplied Sigma.
#> >> Block 1: 6 items (M2PL), loading on dim(s) 1
#> >> Block 2: 6 items (M2PL), loading on dim(s) 2
#> ----------------------------------------------------------------
#> Constructing final data frames...
#> Simulation Complete.
#> Summary: 12 items, 400 examinees, 2 dimensions.
#> ----------------------------------------------------------------
head(mdat$true_params)
#> item_id block model categories a_Dim1 a_Dim2 d guessing
#> 1 item_1 1 M2PL 2 0.8192317 0 0.8491260 0
#> 2 item_2 1 M2PL 2 1.7109341 0 0.5145288 0
#> 3 item_3 1 M2PL 2 1.7031725 0 -1.0493497 0
#> 4 item_4 1 M2PL 2 0.9961762 0 1.4578654 0
#> 5 item_5 1 M2PL 2 1.5794204 0 0.5032779 0
#> 6 item_6 1 M2PL 2 1.2982607 0 1.3152575 0# Independent items plus two testlets
tdat <- sim_tirt(
n_people = 400,
item_structure = list(
list(model = "2PL", n_items = 6),
list(model = "2PLT", n_items = 4, testlet_id = "P1", testlet_var = 0.6),
list(model = "GPCT", n_items = 3, categories = 3, testlet_id = "P2")
)
)
#> ================================================================
#> STARTING MIXED IRT/TRT SIMULATION (N = 400)
#> ================================================================
#> >> Ability (Theta): Generated from N(mean=0.00, sd=1.00).
#> >> Testlet Effects (Gamma):
#> - Testlet 'P1': Generated Gamma ~ N(0, 0.60) [User Var].
#> - Testlet 'P2': Generated Gamma ~ N(0, 0.50) [Default].
#> >> Block 1: 6 items (Model: 2PL, Independent)
#> >> Block 2: 4 items (Model: 2PLT, Testlet: P1)
#> >> Block 3: 3 items (Model: GPCT, Testlet: P2)
#> ================================================================
#> Constructing final data frames...
#> Simulation Complete.
#> Summary: 13 items (2 testlets), 400 examinees.
#> ================================================================
tdat$true_item_params[, c("item_id", "model", "testlet")]
#> item_id model testlet
#> 1 item_1 2PL <NA>
#> 2 item_2 2PL <NA>
#> 3 item_3 2PL <NA>
#> 4 item_4 2PL <NA>
#> 5 item_5 2PL <NA>
#> 6 item_6 2PL <NA>
#> 7 item_7 2PLT P1
#> 8 item_8 2PLT P1
#> 9 item_9 2PLT P1
#> 10 item_10 2PLT P1
#> 11 item_11 GPCT P2
#> 12 item_12 GPCT P2
#> 13 item_13 GPCT P211. Mixture (Latent-Class) IRT
mixture_irt() fits a mixture Rasch or 2PL model, in
which the population is a blend of latent classes with class-specific
item parameters.
set.seed(11)
N <- 300; J <- 8
b1 <- seq(-1.5, 1.5, length.out = J); b2 <- rev(b1)
theta <- rnorm(N); cls <- rep(1:2, each = N / 2)
rmat <- matrix(0, N, J)
for (i in 1:N) {
b <- if (cls[i] == 1) b1 else b2
rmat[i, ] <- rbinom(J, 1, 1 / (1 + exp(-(theta[i] - b))))
}
mdf <- as.data.frame(rmat); names(mdf) <- paste0("I", 1:J)
mix <- mixture_irt(mdf, n_class = 2, model = "Rasch",
control = list(max_iter = 40, verbose = FALSE))
mix$class_params
#> class proportion
#> 1 1 0.4873
#> 2 2 0.5127
mix$model_fit
#> Index Value
#> 1 LogLikelihood -1533.969
#> 2 AIC 3101.938
#> 3 BIC 3164.902
#> 4 n_class 2.000
#> 5 entropy 0.61312. Final Comment
Together with the estimation, calibration, and equating functions
described in the companion vignette, these tools cover the routine
post-estimation workflow of an operational testing program. For details
on any function, use the help system, for example
?score_table.