ICC plot for ordinal IRT using brms

[JLC]

I found this fantastic blog post from @Solomon on plotting ICC curves for binary outcomes.

I’m working with an ordinal IRT model:

fmla_2pl <- bf(Response ~ 1 + 
                                 (1 |i| Item) + 
                                 (1 | D), 
                         disc ~ 1 + (1 |i| Item))

The response is a 4-point agreement scale and I’d like to plot the ICCs for each item.

@paul.buerkner 's IRT tutorial has been helpful for looking at “easiness” and discrimination by item:

eta <- ranef_2pl_Mod$Item[, , "Intercept"] %>%
        data.frame() %>%
        rownames_to_column("Item")

# discriminations
alpha <- ranef_2pl_Mod$Item[, , "disc_Intercept"] %>%
        exp() %>%
        data.frame() %>%
        rownames_to_column("Item")

But I’m struggling to generate a plot.

I’ve tried this:

left_join(post %>%
        select(
                starts_with("r_Item")
        ) %>%
        mutate(iter = 1:n()) %>%
        pivot_longer(starts_with("r_Item")
                     ) %>%
        mutate(Item = str_extract(name, "\\d+")) %>%
        select(-name, eta = value),
        post %>%
                select(
                        starts_with("r_Item__disc")
                ) %>%
                mutate(iter = 1:n()) %>%
                pivot_longer(starts_with("r_Item__disc")
                ) %>%
                mutate(Item = str_extract(name, "\\d+")) %>%
                select(-name, alpha = value)
        , by = c("iter", "Item")
) %>% 
        expand(nesting(iter, Item, eta, alpha),
               theta = seq(from = -10, to = 10, length.out = 10)) %>% 
        mutate(p = inv_logit_scaled(alpha * (eta + theta))) %>% 
        group_by(theta, Item) %>% 
        summarise(p = mean(p)) %>%
        ggplot(aes(x = theta, y = p, color = Item)) +
        geom_line() +
        scale_color_viridis_d(option = "H") +
        labs( 
             x = expression(theta~('Agreeability')),
             y = expression(italic(p)(y==1))) +
        theme_classic()

But the resulting plot looks odd.

Could that just be the data?

Or is it due to @Solomon 's comment:
"You should be able to generalize this workflow to IRT models for data with more than two categories. You’ll just have to be careful about juggling your thresholds. "?

I’m not sure how to include thresholds in the calculations.

[JLC]

I tried a little threshold wrangling, but am unsure I’ve calculated these probabilities correctly as a function of theta:

Fit %>%
        spread_draws(`Intercept[1]`
                     , `Intercept[2]`
                     , `Intercept[3]`
                     , r_Item__disc[Item,]
                     , r_Item[Item, ]
                     , ndraws = 250
        ) %>%
        ungroup() %>%
        select(`Intercept[1]`
               , `Intercept[2]`
               , `Intercept[3]`
               , iter = .draw
               , eta = r_Item
               , Item
               , alpha = r_Item__disc
        ) %>%
        expand(nesting(`Intercept[1]`
                       , `Intercept[2]`
                       , `Intercept[3]`
                       , iter
                       , Item
                       , eta
                       , alpha),
               theta = seq(from = -10, to = 10, length.out = 100)) %>%
        mutate(pk1 = inv_logit_scaled(`Intercept[1]` - (alpha + eta) + theta)
               , pk2 = inv_logit_scaled(`Intercept[2]` - (alpha  + eta) + theta)
               , pk3 = inv_logit_scaled(`Intercept[3]` - (alpha  + eta) + theta)
        ) %>%
        rowwise() %>%
        mutate(pk4 = 1 - (pk1 + pk2 + pk3)) %>%
        ungroup() %>%
        group_by(theta, Item) %>% 
        summarise(p_1 = mean(pk1)
                  , p_2 = mean(pk2)
                  , p_3 = mean(pk3)
                  , p_4 = mean(pk3)) %>%
        ungroup() %>%
        select(Item
               , theta
               , p_1
               , p_2
               , p_3
               , p_4) %>%
        pivot_longer(p_1:p_4
                     , names_to = "Response"
                     , values_to = "p") %>% 
        mutate(Response = factor(recode(Response
                                        , "p_1" = "Strongly Disagree"
                                        , "p_2" = "Disagree"
                                        , "p_3" = "Agree"
                                        , "p_4" = "Strongly Agree")
                                 , levels = c("Strongly Disagree"
                                              , "Disagree"
                                              , "Agree"
                                              , "Strongly Agree"))
               , Item = factor(Item)
               ) %>% 
        ggplot(aes(x = theta, y = p, colour = Item)) +
        facet_wrap(~ Response) +
        geom_line() +
        scale_color_viridis_d(option = "H") +
        labs( 
                x = expression(theta~('Agreeability')),
                y = expression(italic(p)(y==1))) +
        theme_classic()

The Item “easiness” and discrimination plot looks like this: