A reusable Bayesian IRT workflow for instrument validation
Chi Zhang & Maria Pampaka
2026-03-21
Setup and Package/dataset Loading
# Load packages
source(here::here("scripts", "00_packages.R"))Model Formulas: Three Competing Dimensional Structures
In brms, IRT models are specified as generalized linear mixed models (GLMMs). Item parameters are modelled as random intercepts grouped by item, and person parameters as random effects grouped by person id. The three structural hypotheses translate directly into different random effect specifications:
# 1. Unidimensional IRT (UIRT)
# All items load on a single latent trait.
# (1 | item) = item-specific difficulty (random intercept)
# (1 | id) = person ability (random intercept)
bf_uni <- bf(
response ~ 1 + (1 | item) + (1 | id)
)
# 2. Correlated-Traits MIRT (CT-MIRT)
# Items load on one of three correlated dimensions.
# dim_m, dim_l, dim_p are binary indicators (0/1) assigning each item
# to its dimension.
# (0 + dim_m + dim_l + dim_p | id) estimates three person-level traits
# with an unrestricted covariance matrix, allowing correlations between
# dimensions.
bf_ct <- bf(
response ~ 1 + (1 | item) +
(0 + dim_m + dim_l + dim_p | id)
)
# 3. Bi-factor MIRT (BF-MIRT)
# A general factor (1 | id) plus dimension-specific factors.
# The double-bar notation (||) in brms suppresses correlations between
# random effects.
bf_bf <- bf(
response ~ 1 + (1 | item) +
(1 || id) +
(0 + dim_m || id) +
(0 + dim_l || id) +
(0 + dim_p || id)
)Prior Specification
We use weakly informative priors throughout, following recommendations in Bürkner (2021). These provide sufficient regularisation to aid convergence without imposing strong constraints on parameter estimates.
# Standard priors for person and item random effects
priors_uni <- prior("normal(0, 3)", class = "sd", group = "id") +
prior("normal(0, 3)", class = "sd", group = "item")
# For CT models: add LKJ(2) prior on the correlation matrix between
# latent dimensions. LKJ with eta > 1 weakly favours smaller
# correlations, aiding convergence.
priors_ct <- priors_uni + prior(lkj(2), class = "cor")
# For 2PL models: add prior for log-discrimination.
priors_ct_2pl <- priors_ct +
prior("normal(0, 1)", class = "sd", group = "item", dpar = "disc")
priors_bf_2pl <- priors_uni +
prior("normal(0, 1)", class = "sd", group = "item", dpar = "disc")Step 1: Preprocessing sparse response categories
Inspecting the original 5-point scale
# Load the original 5-point scale data
data_1 <- readRDS(here("data", "processed", "data_1.rds"))
# Check response distribution
print(table(data_1$response))##
## 1 2 3 4 5
## 141 913 714 1412 419# Run the plotting script to generate the descriptive plot
# This script creates the ggplot object 'fig_3'
source(here("scripts", "03_fig3.R"))
fig_3
Baseline null_rsm model
# Fit RSM on the original 5-point data
null_rsm <- brm(
formula = bf_uni,
data = data_1,
family = brmsfamily("acat", "logit"),
prior = priors_uni,
chains = 4,
cores = 4,
iter = 4000,
warmup = 2000,
seed = 1234,
file = here("models", "null_rsm")
)
null_rsm <- add_criterion(
null_rsm,
criterion = c("loo", "waic")
)
saveRDS(null_rsm, here("models", "null_rsm.rds"))Summary
summary (null_rsm)## Family: acat
## Links: mu = logit; disc = identity
## Formula: response ~ 1 + (1 | item) + (1 | id)
## Data: data_1 (Number of observations: 3599)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~id (Number of levels: 152)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.49 0.04 0.42 0.58 1.00 2577 4725
##
## ~item (Number of levels: 24)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.14 0.18 0.85 1.55 1.00 1011 2196
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -3.08 0.25 -3.57 -2.61 1.00 583 1356
## Intercept[2] -0.38 0.23 -0.86 0.06 1.00 512 1068
## Intercept[3] -0.56 0.23 -1.02 -0.12 1.00 500 950
## Intercept[4] 2.16 0.24 1.70 2.62 1.00 530 1057
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Collapsing sparse categories
Twelve items had fewer than three responses in the lowest category (“strongly disagree”), providing insufficient information for reliable threshold estimation. We collapse categories 1 and 2 into a single category, reducing the scale from 5 to 4 ordered categories.
# Load the collapsed 4-point scale data
data_2 <- readRDS(here("data", "processed", "data_2.rds"))
# Fit RSM on the collapsed 4-point data
merge_rsm <- brm(
formula = bf_uni,
data = data_2,
family = brmsfamily("acat", "logit"),
prior = priors_uni,
chains = 4,
cores = 4,
iter = 4000,
warmup = 2000,
file = here("models", "merge_rsm")
)
merge_rsm <- add_criterion(
merge_rsm,
criterion = c("loo", "waic")
)
saveRDS(merge_rsm, here("models", "merge_rsm.rds"))summary(merge_rsm)## Family: acat
## Links: mu = logit; disc = identity
## Formula: response ~ 1 + (1 | item) + (1 | id)
## Data: data_2 (Number of observations: 3599)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~id (Number of levels: 152)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.59 0.05 0.51 0.69 1.00 2591 4555
##
## ~item (Number of levels: 24)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.26 0.20 0.94 1.72 1.00 988 2083
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -0.40 0.25 -0.91 0.10 1.00 588 1140
## Intercept[2] -0.56 0.25 -1.06 -0.06 1.01 577 1091
## Intercept[3] 2.27 0.26 1.77 2.79 1.01 598 1322
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).A note on comparing models across different response scales
# You might attempt to compare the two models using LOO-CV:
loo(null_rsm, merge_rsm)## Warning: Not all models have the same y variable. ('yhash' attributes do not match)## Output of model 'null_rsm':
##
## Computed from 8000 by 3599 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -4042.8 51.1
## p_loo 143.7 4.4
## looic 8085.7 102.2
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.6]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'merge_rsm':
##
## Computed from 8000 by 3599 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -3570.8 45.7
## p_loo 147.3 4.1
## looic 7141.5 91.4
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.6]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## merge_rsm 0.0 0.0
## null_rsm -472.1 25.5# However, this comparison should be interpreted with caution.
# The two models predict different response variables (5 categories
# vs 4 categories). Reducing the number of categories mechanically
# increases per-observation predictive density — predicting among
# 4 outcomes is inherently "easier" than predicting among 5.
#
# The decision to collapse categories is instead justified by:
# (1) Descriptive evidence: 12 items had < 3 responses in category 1.
# (2) Diagnostic evidence: the initial RSM produced anomalous threshold
# ordering (see summary above), where the model could not
# meaningfully distinguish between the first two category transitions.Summary
source(here("scripts", "04_fig4.R"))
fig_4
summary (merge_rsm) # but still, the gaping issue remains.## Family: acat
## Links: mu = logit; disc = identity
## Formula: response ~ 1 + (1 | item) + (1 | id)
## Data: data_2 (Number of observations: 3599)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~id (Number of levels: 152)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.59 0.05 0.51 0.69 1.00 2591 4555
##
## ~item (Number of levels: 24)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.26 0.20 0.94 1.72 1.00 988 2083
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -0.40 0.25 -0.91 0.10 1.00 588 1140
## Intercept[2] -0.56 0.25 -1.06 -0.06 1.01 577 1091
## Intercept[3] 2.27 0.26 1.77 2.79 1.01 598 1322
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Step 2: Dimensionality — Unidimensional vs Multidimensional
We compare the three competing dimensional structures using the
collapsed dataset (data_2). All models here use the same
1PL (RSM) specification with the same 4-category response variable,
making LOO-CV comparisons valid.
Fitting the Bi-factor Model (BF-MIRT)
bf_rsm <- brm(
formula = bf_bf,
data = data_2,
family = brmsfamily("acat", "logit"),
prior = priors_uni,
chains = 4,
cores = 4,
iter = 4000,
warmup = 2000,
seed = 1234,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = here("models", "bf_rsm")
)
bf_rsm <- add_criterion(
bf_rsm,
criterion = c("loo", "waic")
)
saveRDS(bf_rsm, here("models", "bf_rsm.rds"))summary (ct_rsm)## Family: acat
## Links: mu = logit; disc = identity
## Formula: response ~ 1 + (1 | item) + (0 + dim_m + dim_l + dim_p | id)
## Data: data_2 (Number of observations: 3599)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~id (Number of levels: 152)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(dim_m) 0.77 0.08 0.61 0.93 1.00 4327 5152
## sd(dim_l) 0.81 0.07 0.69 0.96 1.00 3958 5751
## sd(dim_p) 0.72 0.09 0.55 0.90 1.00 4099 4997
## cor(dim_m,dim_l) 0.36 0.11 0.11 0.56 1.00 1475 2495
## cor(dim_m,dim_p) 0.47 0.12 0.22 0.69 1.00 2827 5135
## cor(dim_l,dim_p) 0.65 0.09 0.45 0.82 1.00 4674 5194
##
## ~item (Number of levels: 24)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.39 0.23 1.03 1.91 1.00 1535 3254
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -0.54 0.30 -1.11 0.05 1.01 1075 2321
## Intercept[2] -0.54 0.29 -1.11 0.04 1.01 1036 2266
## Intercept[3] 2.46 0.30 1.89 3.04 1.01 1088 2434
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary (bf_rsm)## Family: acat
## Links: mu = logit; disc = identity
## Formula: response ~ 1 + (1 | item) + (1 | id) + (0 + dim_m || id) + (0 + dim_l || id) + (0 + dim_p || id)
## Data: data_2 (Number of observations: 3599)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~id (Number of levels: 152)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.57 0.06 0.46 0.69 1.00 2538 3441
## sd(dim_m) 0.59 0.09 0.41 0.77 1.00 2009 3854
## sd(dim_l) 0.57 0.08 0.43 0.72 1.00 2321 4116
## sd(dim_p) 0.35 0.14 0.05 0.59 1.00 893 1250
##
## ~item (Number of levels: 24)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.39 0.23 1.03 1.92 1.00 1442 2736
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -0.55 0.28 -1.10 0.01 1.00 839 1770
## Intercept[2] -0.55 0.28 -1.10 0.01 1.00 828 1781
## Intercept[3] 2.44 0.29 1.88 3.00 1.00 857 1760
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Dimensionality comparison
loo_compare(merge_rsm, ct_rsm, bf_rsm, criterion = "loo")## elpd_diff se_diff
## bf_rsm 0.0 0.0
## ct_rsm -0.2 2.0
## merge_rsm -69.0 13.8# Both CT and BF substantially outperform the unidimensional model,
# confirming that the data have multidimensional structure.
# However, CT and BF are not clearly distinguishable at this stage.Step 3: Freeing Discrimination Parameters
The RSM constrains all item discriminations to be equal. When items differ substantially in how well they differentiate between respondents, this constraint forces the model to find compromise threshold values, which can produce anomalous threshold ordering. Freeing discrimination parameters (moving from 1PL to 2PL) serves a dual diagnostic purpose: it tests the equal-discrimination assumption and reveals items with poor measurement properties.
2PL formula specification
# CT-2pl Formula
bf_ct_2pl <- bf(
response ~ 1 + (1 | item) + (0 + dim_m + dim_l + dim_p | id),
disc ~ 1 + (1 | item)
)
priors_ct_2pl <- priors_ct +
prior("normal(0, 1)", class = "sd", group = "item", dpar = "disc")
# BF-2pl Formula
bf_bf_2pl <- bf(
response ~ 1 + (1 | item) +
(1 | id) +
(0 + dim_m || id) +
(0 + dim_l || id) +
(0 + dim_p || id),
disc ~ 1 + (1 | item)
)
priors_bf_2pl <- priors_uni +
prior("normal(0, 1)", class = "sd", group = "item", dpar = "disc")Fit 2pl (ct_grsm, bf_grsm)
ct_grsm <- brm(
formula = bf_ct_2pl,
data = data_2,
family = brmsfamily("acat", "logit"),
prior = priors_ct_2pl,
chains = 4,
cores = 4,
iter = 4000,
warmup = 2000,
seed = 1234,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = here("models", "ct_grsm")
)
ct_grsm <- add_criterion(ct_grsm, criterion = c("loo", "waic"), overwrite = TRUE)## Warning: Found 6 observations with a pareto_k > 0.7 in model 'ct_grsm'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.## Warning:
## 174 (4.8%) p_waic estimates greater than 0.4. We recommend trying loo instead.## Automatically saving the model object in '/Users/chi/Desktop/BIRT_tutorial_belief/models/ct_grsm.rds'saveRDS(ct_grsm, here("models", "ct_grsm.rds"))
bf_grsm <- brm(
formula = bf_bf_2pl,
data = data_2,
family = brmsfamily("acat", "logit"),
prior = priors_bf_2pl,
chains = 4,
cores = 4,
iter = 4000,
warmup = 2000,
seed = 1234,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = here("models", "bf_grsm")
)
bf_grsm <- add_criterion(bf_grsm, criterion = c("loo", "waic"), overwrite = TRUE)## Warning: Found 8 observations with a pareto_k > 0.7 in model 'bf_grsm'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.## Warning:
## 173 (4.8%) p_waic estimates greater than 0.4. We recommend trying loo instead.## Automatically saving the model object in '/Users/chi/Desktop/BIRT_tutorial_belief/models/bf_grsm.rds'saveRDS(bf_grsm, here("models", "bf_grsm.rds"))Comparing 1PL and 2PL specifications
loo_compare(ct_rsm, bf_rsm, ct_grsm, bf_grsm, criterion = "loo")## elpd_diff se_diff
## ct_grsm 0.0 0.0
## bf_grsm -0.7 1.4
## bf_rsm -140.8 22.6
## ct_rsm -141.0 22.4summary (bf_grsm)## Family: acat
## Links: mu = logit; disc = log
## Formula: response ~ 1 + (1 | item) + (1 | id) + (0 + dim_m || id) + (0 + dim_l || id) + (0 + dim_p || id)
## disc ~ 1 + (1 | item)
## Data: data_2 (Number of observations: 3599)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~id (Number of levels: 152)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.29 0.42 0.60 2.24 1.00 1785 2107
## sd(dim_m) 2.07 0.67 0.97 3.60 1.00 1869 2347
## sd(dim_l) 1.50 0.50 0.69 2.62 1.00 1916 2410
## sd(dim_p) 1.25 0.45 0.53 2.28 1.00 1768 2447
##
## ~item (Number of levels: 24)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 5.11 1.56 2.44 8.46 1.00 2135 2472
## sd(disc_Intercept) 0.79 0.14 0.57 1.10 1.00 2172 3296
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -0.43 0.93 -2.39 1.32 1.00 831 1739
## Intercept[2] 0.16 0.92 -1.71 1.96 1.00 800 1641
## Intercept[3] 5.95 1.99 2.73 10.35 1.00 1571 2263
## disc_Intercept -0.81 0.35 -1.45 -0.10 1.00 1688 2540
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Wright map: inspecting discrimination and difficulty
source(here("scripts", "05_fig5.R"))
fig_5
# The Wright map reveals substantial variation in discrimination
# across items, confirming that the equal-discrimination assumption
# was untenable. Threshold ordering has also improved (see paper
# Table for threshold estimates).
#
# Critically, four items (Bm_1_r to Bm_4_r) stand out with:
# - Discrimination estimates near zero
# - Extremely wide credible intervals on difficulty
# These are candidates for removal.Cross-check: Bayesian infit/outfit
The decision to flag these items came from freeing the discrimination
parameters. As a complementary check, we can also screen items with
classical infit/outfit statistics. These have a Bayesian,
posterior-predictive form
(easyRaschBayes::infit_statistic): for each posterior draw
it compares the variance-weighted standardised residuals of the observed
data against those of model-replicated data, so it is a posterior
predictive check rather than a point estimate. We run it on the RSM
stage (merge_rsm), where infit/outfit are the canonical
Rasch-family diagnostic. Misfit is read from the posterior predictive
p-value (ppp < 0.05) and from the mean-square value
(acceptable in [0.5, 1.5]).
# install once: remotes::install_github("pgmj/easyRaschBayes")
source(here("scripts", "10_itemfit.R"))## item infit_MNSQ infit_ppp outfit_MNSQ outfit_ppp
## Bm_2_r 1.651 0.026 1.484 0.120
## Bm_4_r 1.374 0.002 1.323 0.030
## Bm_1_r 1.328 0.040 1.279 0.130
## Bl_6_r 1.093 0.188 1.057 0.298
## Bm_3_r 1.092 0.308 1.174 0.232
## Bl_11_r 1.090 0.176 1.083 0.216
## Bl_10_r 1.082 0.226 1.075 0.276
## Bl_9_r 1.077 0.214 1.073 0.234
## Bl_2 1.070 0.282 1.174 0.114
## Bm_5 1.026 0.398 1.050 0.334
## Bp_2_r 1.022 0.404 0.999 0.498
## Bm_8 1.013 0.430 1.059 0.318
## Bl_7_r 1.003 0.500 1.012 0.458
## Bl_8_r 0.982 0.552 0.984 0.528
## Bp_1_r 0.970 0.598 0.984 0.530
## Bl_5 0.940 0.698 0.992 0.540
## Bl_3_r 0.939 0.718 0.929 0.736
## Bm_6 0.886 0.796 0.925 0.708
## Bp_5 0.883 0.810 0.949 0.624
## Bm_7 0.820 0.898 0.870 0.816
## Bl_4_r 0.808 0.978 0.813 0.964
## Bl_1_r 0.804 0.982 0.817 0.962
## Bp_4 0.672 0.994 0.741 0.968
## Bp_3 0.619 1.000 0.645 0.998Infit/outfit independently flag Bm_1_r, Bm_2_r and Bm_4_r. Bm_3_r, by contrast, passes item fit (infit ~1.09) and is only revealed once discrimination is freed. Item-fit and discrimination diagnostics thus have complementary sensitivities; here we proceed with the discrimination-based decision and remove all four, while a reader could reasonably treat Bm_3_r as a borderline case.
Step 4: Item removal and different link function comparisons
# This provides data3 that is without bm1r~bm4r
data_3 <- readRDS(here("data", "processed", "data_3.rds"))With dimensionality confirmed and problematic items removed, we now compare IRT models defined by different ordinal link functions.
Fitting the six models
Three link function families are compared, each crossed with two dimensional structures (CT and BF), yielding six models:
# --- GRSM: Generalised Rating Scale Model ---
# Link: adjacent-category logit
# Thresholds: shared across items (only item intercepts vary)
# In brms: family = brmsfamily("acat", "logit")
# Formula: response ~ 1 + (1 | item) + ...
#
# --- GRM: Graded Response Model ---
# Link: cumulative logit
# Thresholds: item-specific (estimated via cumulative family)
# In brms: family = brmsfamily("cumulative", "logit")
# Formula: same as GRSM; the family handles the parameterisation
#
# --- GPCM: Generalised Partial Credit Model ---
# Link: adjacent-category logit
# Thresholds: item-specific (via category-specific effects)
# In brms: family = brmsfamily("acat", "logit")
# Formula: response ~ 1 + (cs(1) | item) + ...
# cs(1) allows thresholds to vary by item
# --- 1. GRSM ---
ct_grsm_20 <- brm(
formula = bf_ct_2pl,
data = data_3,
family = brmsfamily("acat", "logit"),
prior = priors_ct_2pl,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE),
chains = 4, cores = 4, iter = 4000, warmup = 2000, seed = 1234,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = here("models", "ct_grsm_20")
)
ct_grsm_20 <- add_criterion(ct_grsm_20, "loo", overwrite = TRUE)## Warning: Found 5 observations with a pareto_k > 0.7 in model 'ct_grsm_20'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.## Automatically saving the model object in '/Users/chi/Desktop/BIRT_tutorial_belief/models/ct_grsm_20.rds'saveRDS(ct_grsm_20, here("models", "ct_grsm_20.rds"))
bf_grsm_20 <- brm(
formula = bf_bf_2pl,
data = data_3,
family = brmsfamily("acat", "logit"),
prior = priors_bf_2pl,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE),
chains = 4, cores = 4, iter = 4000, warmup = 2000, seed = 1234,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = here("models", "bf_grsm_20")
)
bf_grsm_20 <- add_criterion(bf_grsm_20, "loo", overwrite = TRUE)## Warning: Found 3 observations with a pareto_k > 0.7 in model 'bf_grsm_20'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.## Automatically saving the model object in '/Users/chi/Desktop/BIRT_tutorial_belief/models/bf_grsm_20.rds'saveRDS(bf_grsm_20, here("models", "bf_grsm_20.rds"))
# --- 2. GRM ---
ct_grm_20 <- brm(
formula = bf_ct_2pl,
data = data_3,
family = brmsfamily("cumulative", "logit"),
prior = priors_ct_2pl,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE),
chains = 4, cores = 4, iter = 4000, warmup = 2000, seed = 1234,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = here("models", "ct_grm_20")
)
ct_grm_20 <- add_criterion(ct_grm_20, "loo", overwrite = TRUE)## Warning: Found 2 observations with a pareto_k > 0.7 in model 'ct_grm_20'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.## Automatically saving the model object in '/Users/chi/Desktop/BIRT_tutorial_belief/models/ct_grm_20.rds'saveRDS(ct_grm_20, here("models", "ct_grm_20.rds"))
bf_grm_20 <- brm(
formula = bf_bf_2pl,
data = data_3,
family = brmsfamily("cumulative", "logit"),
prior = priors_bf_2pl,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE),
chains = 4, cores = 4, iter = 4000, warmup = 2000, seed = 1234,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = here("models", "bf_grm_20")
)
bf_grm_20 <- add_criterion(bf_grm_20, "loo", overwrite = TRUE)## Warning: Found 1 observations with a pareto_k > 0.7 in model 'bf_grm_20'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.## Automatically saving the model object in '/Users/chi/Desktop/BIRT_tutorial_belief/models/bf_grm_20.rds'saveRDS(bf_grm_20, here("models", "bf_grm_20.rds"))
# --- 3. GPCM ---
bf_ct_gpcm <- bf(
response ~ 1 + (cs(1) | item) + (0 + dim_m + dim_l + dim_p | id),
disc ~ 1 + (1 | item)
)
bf_bf_gpcm <- bf(
response ~ 1 + (cs(1) | item) +
(1 | id) + (0 + dim_m || id) + (0 + dim_l || id) + (0 + dim_p || id),
disc ~ 1 + (1 | item)
)
ct_gpcm_20 <- brm(
formula = bf_ct_gpcm,
data = data_3,
family = brmsfamily("acat", "logit"),
prior = priors_ct_2pl,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE),
chains = 4, cores = 4, iter = 4000, warmup = 2000, seed = 1234,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = here("models", "ct_gpcm_20")
)
ct_gpcm_20 <- add_criterion(ct_gpcm_20, "loo", overwrite = TRUE)## Warning: Found 9 observations with a pareto_k > 0.7 in model 'ct_gpcm_20'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.## Automatically saving the model object in '/Users/chi/Desktop/BIRT_tutorial_belief/models/ct_gpcm_20.rds'saveRDS(ct_gpcm_20, here("models", "ct_gpcm_20.rds"))
bf_gpcm_20 <- brm(
formula = bf_bf_gpcm,
data = data_3,
family = brmsfamily("acat", "logit"),
prior = priors_bf_2pl,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE),
chains = 4, cores = 4, iter = 4000, warmup = 2000, seed = 1234,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = here("models", "bf_gpcm_20")
)
bf_gpcm_20 <- add_criterion(bf_gpcm_20, "loo", overwrite = TRUE)## Warning: Found 7 observations with a pareto_k > 0.7 in model 'bf_gpcm_20'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.## Automatically saving the model object in '/Users/chi/Desktop/BIRT_tutorial_belief/models/bf_gpcm_20.rds'saveRDS(bf_gpcm_20, here("models", "bf_gpcm_20.rds"))Inspecting individual LOO diagnostics
Individual LOO summaries reveal Pareto-k diagnostics for each model. High Pareto-k values (> 0.7) indicate observations where the LOO approximation may be unreliable. A model with the best ELPD but many high Pareto-k values should be treated with caution.
loo(ct_grsm_20)##
## Computed from 8000 by 2993 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -2855.3 41.0
## p_loo 326.1 10.7
## looic 5710.6 81.9
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.8]).
##
## Pareto k diagnostic values:
## Count Pct. Min. ESS
## (-Inf, 0.7] (good) 2988 99.8% 312
## (0.7, 1] (bad) 5 0.2% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.loo(bf_grsm_20)##
## Computed from 8000 by 2993 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -2857.6 41.0
## p_loo 327.0 10.7
## looic 5715.3 82.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.9]).
##
## Pareto k diagnostic values:
## Count Pct. Min. ESS
## (-Inf, 0.7] (good) 2990 99.9% 265
## (0.7, 1] (bad) 3 0.1% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.loo(ct_grm_20)##
## Computed from 8000 by 2993 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -2851.2 40.7
## p_loo 342.1 10.4
## looic 5702.3 81.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
##
## Pareto k diagnostic values:
## Count Pct. Min. ESS
## (-Inf, 0.7] (good) 2991 99.9% 385
## (0.7, 1] (bad) 2 0.1% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.loo(bf_grm_20)##
## Computed from 8000 by 2993 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -2849.4 40.7
## p_loo 339.3 10.3
## looic 5698.7 81.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.9]).
##
## Pareto k diagnostic values:
## Count Pct. Min. ESS
## (-Inf, 0.7] (good) 2992 100.0% 394
## (0.7, 1] (bad) 1 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.loo(ct_gpcm_20)##
## Computed from 8000 by 2993 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -2821.3 41.8
## p_loo 364.4 12.1
## looic 5642.6 83.6
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.8]).
##
## Pareto k diagnostic values:
## Count Pct. Min. ESS
## (-Inf, 0.7] (good) 2984 99.7% 118
## (0.7, 1] (bad) 9 0.3% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.loo(bf_gpcm_20)##
## Computed from 8000 by 2993 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -2836.5 41.9
## p_loo 361.4 11.8
## looic 5673.1 83.8
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.2, 1.9]).
##
## Pareto k diagnostic values:
## Count Pct. Min. ESS
## (-Inf, 0.7] (good) 2986 99.8% 48
## (0.7, 1] (bad) 7 0.2% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.loo_compare(
ct_grsm_20, bf_grsm_20,
ct_grm_20, bf_grm_20,
ct_gpcm_20, bf_gpcm_20
)## elpd_diff se_diff
## ct_gpcm_20 0.0 0.0
## bf_gpcm_20 -15.2 5.3
## bf_grm_20 -28.1 10.1
## ct_grm_20 -29.9 9.9
## ct_grsm_20 -34.0 8.8
## bf_grsm_20 -36.3 8.9# Stacking weights across all six models
model_weights(
ct_grsm_20, bf_grsm_20,
ct_grm_20, bf_grm_20,
ct_gpcm_20, bf_gpcm_20
)## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.## ct_grsm_20 bf_grsm_20 ct_grm_20 bf_grm_20 ct_gpcm_20 bf_gpcm_20
## 3.683987e-05 2.570583e-06 6.594964e-06 2.396761e-01 7.601135e-01 1.643497e-04Clean loo
Moment matching and exact reloo refits are
computationally expensive (each refit re-runs the full sampler, ~30
min). The code that produced the calibrated LOO objects is shown below;
we then load the saved results.
# CT-GRM: moment matching + exact reloo for the few high-Pareto-k points
loo_ct_grm_20 <- loo(ct_grm_20, moment_match = TRUE, reloo = TRUE)
saveRDS(loo_ct_grm_20, here("models", "loo", "loo_ct_grm_20.rds"))
# BF-GRM: moment matching
loo_bf_grm_20 <- loo(bf_grm_20, moment_match = TRUE)
saveRDS(loo_bf_grm_20, here("models", "loo", "loo_bf_grm_20.rds"))loo_ct_grm_20 <- readRDS(here("models", "loo", "loo_ct_grm_20.rds"))
loo_bf_grm_20 <- readRDS(here("models", "loo", "loo_bf_grm_20.rds"))Clean compare
loo_compare(loo_bf_grm_20, loo_ct_grm_20)## elpd_diff se_diff
## bf_grm_20 0.0 0.0
## ct_grm_20 -1.7 1.2loo_model_weights(list(bf_grm = loo_bf_grm_20, ct_grm = loo_ct_grm_20))## Method: stacking
## ------
## weight
## bf_grm 1.000
## ct_grm 0.000Step 5: Final model Wright map
source(here("scripts", "07_fig6.R"))
fig_6
pp_check(bf_grm_20,type = "bars", ndraws = 100)
pp_check(bf_grm_20, type = "bars_grouped", group = "item", ndraws = 100)
Reliability of the final model
Finally, we check the reliability of the final model with two
complementary coefficients. Because the instrument is bi-factor, no
single latent dimension is meant to carry the measurement on its own, so
we report (i) the model-based omega coefficients
(omega_total, the reliability of the total score, with
omega_H and ECV describing the dominant
general factor), computed per posterior draw, and (ii) the per-dimension
person reliability (empirical/EAP and the Bayesian RMU of Bignardi,
Kievit & Bürkner, 2025). The two views answer different questions:
the EAP/RMU figures (a conservative lower bound under shrinkage at this
sample size) show that any single dimension is insufficient on its own,
whereas omega_total shows that the dimensions together make
the total score highly reliable.
# RMU (block C) needs: remotes::install_github("pgmj/easyRaschBayes")
source(here("scripts", "09_reliability.R"))## This is posterior version 1.6.0##
## Attaching package: 'posterior'## The following object is masked from 'package:bayesplot':
##
## rhat## The following objects are masked from 'package:stats':
##
## mad, sd, var## The following objects are masked from 'package:base':
##
## %in%, match## === (A) Model-based omega (posterior median [95% CrI]) ===
## median L95 U95
## omega_total 0.898 0.869 0.921
## omega_H 0.613 0.483 0.717
## ECV_general 0.683 0.548 0.787## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.##
## === (B) Per-dimension empirical (EAP) reliability ===
## factor reliability
## General (belief) 0.555
## Specific: mathematics 0.564
## Specific: learners 0.501
## Specific: pedagogy 0.293## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.##
## === (C) RMU reliability [95% HDCI] (easyRaschBayes) ===
## factor rmu L95 U95
## General (belief) 0.611 0.509 0.712
## Specific: mathematics 0.622 0.530 0.709
## Specific: learners 0.551 0.441 0.659
## Specific: pedagogy 0.319 0.135 0.484Step 6: Covariates and Differential Item Functioning
DIF analysis tests whether individual items function differently across demographic subgroups. In the Bayesian framework, this is implemented by adding group variables as covariates with item-level random slopes, enabling partial pooling across items.
DIF model 1: Gender only
bf_dif <- bf(
response ~ 1 + gender + (0 + gender | item) +
(1 | id) + (0 + dim_m || id) + (0 + dim_l || id) + (0 + dim_p || id),
disc ~ 1 + (1 | item)
)
bf_grm_dif <- brm(
formula = bf_dif,
data = data_3,
family = brmsfamily("cumulative", "logit"),
prior = priors_bf_2pl,
backend = "cmdstanr",
chains = 4, cores = 4, iter = 4000, warmup = 2000, seed = 1234,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = here("models", "bf_grm_dif")
)DIF model 2: Multiple covariates
bf_dif2 <- bf(
response ~ 1 + gender + District +
(1 + gender + District || item) +
(1 | id) + (0 + dim_m || id) + (0 + dim_l || id) + (0 + dim_p || id),
disc ~ 1 + (1 | item)
)
bf_grm_dif2 <- brm(
formula = bf_dif2,
data = data_3,
family = brmsfamily("cumulative", "logit"),
prior = priors_bf_2pl,
backend = "cmdstanr",
chains = 4, cores = 4, iter = 4000, warmup = 2000, seed = 1234,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = here("models", "bf_grm_dif2")
)DIF results
summary(bf_grm_dif)## Family: cumulative
## Links: mu = logit; disc = log
## Formula: response ~ 1 + gender + (0 + gender | item) + (1 | id) + (0 + dim_m || id) + (0 + dim_l || id) + (0 + dim_p || id)
## disc ~ 1 + (1 | item)
## Data: data_3 (Number of observations: 2933)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~id (Number of levels: 148)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.82 0.54 0.93 3.03 1.00 1553 2372
## sd(dim_m) 2.85 0.82 1.51 4.73 1.00 1657 2542
## sd(dim_l) 2.13 0.62 1.10 3.53 1.00 1634 2535
## sd(dim_p) 1.64 0.53 0.77 2.84 1.00 1503 2464
##
## ~item (Number of levels: 20)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(genderFemale) 3.06 0.89 1.58 5.07 1.00 1667 2542
## sd(genderMale) 3.09 0.90 1.62 5.12 1.00 1684 2566
## sd(disc_Intercept) 0.36 0.07 0.24 0.52 1.00 2317 4368
## cor(genderFemale,genderMale) 0.98 0.02 0.93 1.00 1.00 3497 4206
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -4.81 1.45 -8.07 -2.34 1.00 1555 2462
## Intercept[2] -1.33 0.80 -3.06 0.10 1.00 1397 2710
## Intercept[3] 5.85 1.80 2.96 9.93 1.00 1512 2565
## disc_Intercept -0.59 0.29 -1.12 -0.00 1.00 1631 2619
## genderMale -0.86 0.50 -1.99 -0.02 1.00 3442 4275
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(bf_grm_dif2)## Family: cumulative
## Links: mu = logit; disc = log
## Formula: response ~ 1 + gender + District + (1 + gender + District || item) + (1 | id) + (0 + dim_m || id) + (0 + dim_l || id) + (0 + dim_p || id)
## disc ~ 1 + (1 | item)
## Data: data_3 (Number of observations: 2913)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~id (Number of levels: 147)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 2.59 0.68 1.43 4.07 1.00 1346 2318
## sd(dim_m) 4.04 1.02 2.26 6.23 1.00 1404 2188
## sd(dim_l) 3.07 0.81 1.68 4.83 1.00 1414 2374
## sd(dim_p) 2.29 0.70 1.09 3.85 1.00 1161 1368
##
## ~item (Number of levels: 20)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 4.45 1.19 2.48 7.12 1.00 1515 2621
## sd(genderMale) 0.53 0.35 0.03 1.33 1.00 2443 3673
## sd(District2) 0.59 0.43 0.02 1.63 1.00 2529 3314
## sd(District3) 0.40 0.33 0.01 1.25 1.00 3083 3484
## sd(disc_Intercept) 0.35 0.07 0.24 0.52 1.00 2437 4484
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -7.21 1.93 -11.27 -3.90 1.00 1367 2465
## Intercept[2] -2.23 1.16 -4.62 -0.13 1.00 1256 2705
## Intercept[3] 8.06 2.33 4.16 13.16 1.00 1200 2173
## disc_Intercept -0.96 0.26 -1.43 -0.42 1.00 1315 2211
## genderMale -1.20 0.67 -2.62 -0.04 1.00 2985 3846
## District2 0.02 0.76 -1.51 1.55 1.00 3240 4045
## District3 -1.29 0.74 -2.89 0.03 1.00 2720 3228
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).source(here("scripts", "08_fig7.R"))
fig_7