I’m trying to extract item-specific thresholds (i.e., difficulty/easiness parameters) from an item response model in brms. In the following example, the simulated data have 3 response categories (1, 2, and 3), so the model is similar to a graded response model, and I would like to extract/calculate the threshold parameter for the boundary from 1-2 (i.e., b1) and the threshold parameter for the boundary from 2-3 (i.e., b2).
First, I simulate the data for the IRT model:
library("mirt")
library("brms")
library("tidyverse")
# Simulate data for graded response model with 20 items and 3 response categories
set.seed(52242)
a <- matrix(rlnorm(20,.2,.3))
diffs <- t(apply(matrix(runif(20*2, .3, 1), 20), 1, cumsum))
diffs <- -(diffs - rowMeans(diffs))
d <- diffs + rnorm(20)
dat_wide <- simdata(a, d, 500, itemtype = "graded")
dat_wide <- dat_wide + 1
dat_wide <- data.frame(dat_wide)
dat_wide$ID <- 1:nrow(dat_wide)
dat_wide <- dat_wide %>%
select(ID, everything())
dat_long <- dat_wide %>%
as.data.frame() %>%
pivot_longer(
cols = Item_1:Item_20,
names_to = "item",
values_to = "response")
Here is the graded response model fitted in mirt:
# IRT Model in mirt
mirtModel <- mirt(
data = dat_wide %>% select(starts_with("Item_")),
model = 1,
itemtype = "graded",
SE = TRUE)
Here are the item parameters (a = discrimination, b1 = difficulty for boundary from 1-2; b2 = difficulty for boundary from 2-3):
> coef(mirtModel, simplify = TRUE, IRTpars = TRUE)
$items
a b1 b2
Item_1 1.549 -0.844 -0.340
Item_2 0.637 -0.513 -0.001
Item_3 1.071 -0.958 -0.582
Item_4 1.067 -0.215 0.324
Item_5 2.219 0.555 0.763
Item_6 1.331 0.079 0.357
Item_7 1.881 0.163 0.499
Item_8 0.785 -1.001 -0.086
Item_9 0.996 0.343 1.050
Item_10 0.966 0.180 0.756
Item_11 0.978 0.645 0.985
Item_12 1.405 0.113 0.393
Item_13 1.762 1.508 1.717
Item_14 1.089 -2.364 -1.688
Item_15 0.774 -0.771 -0.301
Item_16 0.687 0.872 1.517
Item_17 1.057 -0.159 0.131
Item_18 0.744 0.046 0.748
Item_19 1.739 -0.470 0.037
Item_20 0.839 -0.502 0.208
Here are the equivalent item parameters in intercept-slope form:
> coef(mirtModel, simplify = TRUE, IRTpars = FALSE)
$items
a1 d1 d2
Item_1 1.549 1.307 0.527
Item_2 0.637 0.327 0.001
Item_3 1.071 1.026 0.623
Item_4 1.067 0.229 -0.346
Item_5 2.219 -1.232 -1.693
Item_6 1.331 -0.105 -0.475
Item_7 1.881 -0.307 -0.938
Item_8 0.785 0.786 0.068
Item_9 0.996 -0.342 -1.046
Item_10 0.966 -0.174 -0.730
Item_11 0.978 -0.631 -0.963
Item_12 1.405 -0.159 -0.552
Item_13 1.762 -2.657 -3.024
Item_14 1.089 2.575 1.839
Item_15 0.774 0.597 0.233
Item_16 0.687 -0.599 -1.042
Item_17 1.057 0.168 -0.138
Item_18 0.744 -0.034 -0.557
Item_19 1.739 0.817 -0.064
Item_20 0.839 0.421 -0.175
Here is a brms model (using intercept-slope parameterization, where beta is item easiness, alpha is item discrimination, and theta is the person parameter):
# IRT Model in brms
formula_grm2PL <- bf(
response ~ beta + alpha * theta,
theta ~ 0 + (1 | ID),
beta ~ 1 + (1 | item),
alpha ~ 1 + (1 | item),
nl = TRUE
)
# Model priors
get_prior(
formula_grm2PL,
data = dat_long,
family = cumulative(link = "logit", threshold = "flexible")
)
model_prior <-
prior("", class = "b", nlpar = "alpha", lb = 0) +
prior("constant(1)", class = "sd", group = "ID", nlpar = "theta")
fit_grm2PL <- brm(
formula = formula_grm2PL,
data = dat_long,
family = cumulative(link = "logit", threshold = "flexible"),
prior = model_prior,
seed = 52242,
cores = 4)
Here is the model output:
> summary(fit_grm2PL)
Family: cumulative
Links: mu = logit
Formula: response ~ beta + alpha * theta
theta ~ 0 + (1 | ID)
beta ~ 1 + (1 | item)
alpha ~ 1 + (1 | item)
Data: dat_long (Number of observations: 10000)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~ID (Number of levels: 500)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(theta_Intercept) 1.00 0.00 1.00 1.00 NA NA NA
~item (Number of levels: 20)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(beta_Intercept) 1.03 0.19 0.73 1.46 1.01 828 1402
sd(alpha_Intercept) 0.66 0.45 0.28 1.72 1.53 7 35
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept[1] -0.34 2.00 -4.45 3.57 1.00 2801 2291
Intercept[2] 0.18 2.00 -3.92 4.11 1.00 2800 2292
beta_Intercept -0.22 2.01 -4.29 3.68 1.00 2683 2181
alpha_Intercept 0.90 0.49 0.01 1.38 1.54 7 30
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).
Putting aside that the model hasn’t fully converged, how can I extract or calculate b1 and b2 from this model? As far as I can tell, the model doesn’t estimate b1/b2 directly, so I’m wondering whether/how I can calculate b1 and b2 from the available information. I can extract a given parameter (e.g., beta) for each item like this:
newData <- data.frame(
item = unique(dat_long$item)
)
itemEasiness <- posterior_epred(
fit_grm2PL,
nlpar = "beta",
newdata = newData,
re_formula = ~ (1 | item))
Here is the beta (item easiness) for each item:
> colMeans(itemEasiness, na.rm = TRUE)
[1] 0.77273459 0.09772734 0.77438866 -0.12948904 -1.44226197 -0.36620740 -0.64832907 0.32993252 -0.73160887
[10] -0.51793267 -0.90402998 -0.42919358 -2.78840530 2.01113821 0.35868236 -0.91325270 -0.05345497 -0.37049963
[19] 0.26500245 0.04817425
In case it’s helpful, this is what ChatGPT suggests for calculating b1 and b2:
# Get posterior draws
post <- as_draws_df(fit_grm2PL)
# Example: compute thresholds for one item
item <- "Item_1"
beta_i <- post[[paste0("r_item__beta[", item, ",Intercept]")]] + post[["b_beta_Intercept"]]
alpha_i <- post[[paste0("r_item__alpha[", item, ",Intercept]")]] + post[["b_alpha_Intercept"]]
# Inspect threshold variable names
names(post)[grep(item, names(post))]
# For example, if thresholds are stored as "Intercept[1]", "Intercept[2]", etc.
int1_i <- post[["Intercept[1]"]]
int2_i <- post[["Intercept[2]"]]
# Compute item thresholds (b1, b2)
b1 <- (int1_i - beta_i) / alpha_i
b2 <- (int2_i - beta_i) / alpha_i
# Summarize
> median(b1)
[1] -0.7493667
> median(b2)
[1] -0.3870398
Is that approach correct? In any case, I prefer to model beta, alpha, and theta explicitly so I can include predictors/covariates of those parameters, so it would be helpful for any solution to generalize to cases where I have predictors of those parameters (and are thus not relying merely on the intercepts of those parameters). The above approach relies on the intercepts of those parameters, so it would not meet my goals.
Thanks in advance!