This reference covers simulate_irt(), which generates item-response data under a known Item Response Theory model with ground-truth item and person parameters. It supports the 1PL (Rasch), 2PL, and 3PL dichotomous models and the graded response model (GRM) for ordered polytomous items. Each example below is small, fast, and seeded so the results are reproducible.

simulate_irt()

Generate item-response data under a known IRT model. The returned saqr_sim object carries the simulated $data (a person-by-item response grid) and the true item/person parameters in $params, so parameter recovery can be checked directly.

Signature

args(simulate_irt)
#> function (n_persons = 500, n_items = 20, model = c("2PL", "1PL", 
#>     "3PL", "GRM"), n_categories = 2, a = NULL, b = NULL, c = NULL, 
#>     theta = NULL, seed = NULL) 
#> NULL

Example — 1PL (Rasch)

In the 1PL model every item shares the same discrimination (fixed at 1); items differ only in difficulty.

fit_1pl <- simulate_irt(n_persons = 300, n_items = 10, model = "1PL", seed = 1)
fit_1pl
#> saqr_sim [irt]  300 x 10  (seed=1)
#>   params: model, a, b, c, theta, thresholds, n_categories 
#>   cols:   item1, item2, item3, item4, item5, item6, ..., item10
head(fit_1pl$data)
str(fit_1pl$params, max.level = 1)
#> List of 7
#>  $ model       : chr "1PL"
#>  $ a           : Named num [1:10] 1 1 1 1 1 1 1 1 1 1
#>   ..- attr(*, "names")= chr [1:10] "item1" "item2" "item3" "item4" ...
#>  $ b           : Named num [1:10] 1.0436 0.0991 -0.4541 -0.6558 -0.0359 ...
#>   ..- attr(*, "names")= chr [1:10] "item1" "item2" "item3" "item4" ...
#>  $ c           : Named num [1:10] 0 0 0 0 0 0 0 0 0 0
#>   ..- attr(*, "names")= chr [1:10] "item1" "item2" "item3" "item4" ...
#>  $ theta       : num [1:300] -0.626 0.184 -0.836 1.595 0.33 ...
#>  $ thresholds  : NULL
#>  $ n_categories: int 2

Example — 2PL

The 2PL model adds an item-specific discrimination parameter a.

fit_2pl <- simulate_irt(n_persons = 400, n_items = 15, model = "2PL", seed = 2)
head(fit_2pl$data)
str(fit_2pl$params, max.level = 1)
#> List of 7
#>  $ model       : chr "2PL"
#>  $ a           : Named num [1:15] 0.938 0.442 0.738 0.78 1.293 ...
#>   ..- attr(*, "names")= chr [1:15] "item1" "item2" "item3" "item4" ...
#>  $ b           : Named num [1:15] -1.313 1.108 -1.008 1.576 0.291 ...
#>   ..- attr(*, "names")= chr [1:15] "item1" "item2" "item3" "item4" ...
#>  $ c           : Named num [1:15] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..- attr(*, "names")= chr [1:15] "item1" "item2" "item3" "item4" ...
#>  $ theta       : num [1:400] -0.8969 0.1848 1.5878 -1.1304 -0.0803 ...
#>  $ thresholds  : NULL
#>  $ n_categories: int 2

Example — 3PL

The 3PL model adds a lower-asymptote (guessing) parameter c (auto-set to 0.2).

fit_3pl <- simulate_irt(n_persons = 400, n_items = 12, model = "3PL", seed = 3)
head(fit_3pl$data)
fit_3pl$params$c
#>  item1  item2  item3  item4  item5  item6  item7  item8  item9 item10 item11 
#>    0.2    0.2    0.2    0.2    0.2    0.2    0.2    0.2    0.2    0.2    0.2 
#> item12 
#>    0.2

Example — GRM (graded response, 4 ordered categories)

For polytomous items, responses take values 0 ... n_categories - 1 and the true category thresholds are returned in $params$thresholds.

fit_grm <- simulate_irt(n_persons = 400, n_items = 8, model = "GRM",
                        n_categories = 4, seed = 4)
head(fit_grm$data)
str(fit_grm$params, max.level = 1)
#> List of 7
#>  $ model       : chr "GRM"
#>  $ a           : Named num [1:8] 0.701 1.147 0.931 1.276 0.803 ...
#>   ..- attr(*, "names")= chr [1:8] "item1" "item2" "item3" "item4" ...
#>  $ b           : Named num [1:8] 0.77 1.588 0.351 -0.534 -0.285 ...
#>   ..- attr(*, "names")= chr [1:8] "item1" "item2" "item3" "item4" ...
#>  $ c           : Named num [1:8] 0 0 0 0 0 0 0 0
#>   ..- attr(*, "names")= chr [1:8] "item1" "item2" "item3" "item4" ...
#>  $ theta       : num [1:400] 0.217 -0.542 0.891 0.596 1.636 ...
#>  $ thresholds  : num [1:8, 1:3] -0.73 0.088 -1.149 -2.034 -1.785 ...
#>  $ n_categories: int 4
fit_grm$params$thresholds
#>             [,1]        [,2]      [,3]
#> [1,] -0.73035779  0.76964221 2.2696422
#> [2,]  0.08798766  1.58798766 3.0879877
#> [3,] -1.14865030  0.35134970 1.8513497
#> [4,] -2.03433804 -0.53433804 0.9656620
#> [5,] -1.78487864 -0.28487864 1.2151214
#> [6,] -1.55272416 -0.05272416 1.4472758
#> [7,] -1.43555660  0.06444340 1.5644434
#> [8,] -2.07866676 -0.57866676 0.9213332

Parameter recovery

For a 2PL fit, item proportion-correct should track true item difficulty b negatively (easier items are answered correctly more often), and person total score should track true ability theta positively.

prop_correct <- colMeans(as.data.frame(fit_2pl))
cor(prop_correct, fit_2pl$params$b)        # strongly negative
#> [1] -0.9722128

total_score <- rowSums(as.data.frame(fit_2pl))
cor(total_score, fit_2pl$params$theta)     # strongly positive
#> [1] 0.8258373

A plot of item proportion-correct against true difficulty b makes the recovery visible: the points fall along a clear downward trend.

plot(fit_2pl$params$b, prop_correct,
     pch = 19, col = "#2c7fb8",
     xlab = "True item difficulty  b",
     ylab = "Observed proportion correct",
     main = sprintf("2PL recovery: cor = %.2f",
                    cor(prop_correct, fit_2pl$params$b)))
abline(lm(prop_correct ~ fit_2pl$params$b), col = "#d95f0e", lwd = 2)

Optional: estimating with mirt

When the mirt package is installed, the simulated data can be fit to confirm the generating parameters are recoverable by an external estimator. This chunk is skipped if mirt is not available.

mod <- mirt::mirt(as.data.frame(fit_2pl), 1, itemtype = "2PL", verbose = FALSE)
est_b <- mirt::coef(mod, simplify = TRUE, IRTpars = TRUE)$items[, "b"]
cor(est_b, fit_2pl$params$b)     # estimated vs. true difficulty
#> [1] 0.9882692