Partial Credit Model Analysis with easyRaschBayes

Overview

This vignette demonstrates a Partial Credit Model (PCM) Rasch analysis workflow using easyRaschBayes. The analysis uses the eRm::pcmdat2 dataset, a small polytomous item response data set included in the eRm package.

All functions in this package work with a Bayesian brms model object fitted with the acat (adjacent categories) family, which parameterises the PCM. Dichotomous Rasch models can also be fit using brms and analyzed with the functions in this package. A code example is available here.

There is brief text explaining the output and interpretation in this vignette. For a more extensive treatment of the Rasch analysis steps, please see the easyRasch vignette.

Data Preparation

library(easyRaschBayes)
library(brms)
library(dplyr)
library(tidyr)
library(tibble)
library(ggplot2)

eRm::pcmdat2 is in wide format with item responses coded 0, 1, 2, …. The brms acat family expects response categories starting at 1, so we add 1 to every response before reshaping to long format.

df_pcm <- eRm::pcmdat2 %>%
  mutate(across(everything(), ~ .x + 1)) %>%
  rownames_to_column("id") %>%
  pivot_longer(!id, names_to = "item", values_to = "response")

head(df_pcm)
#> # A tibble: 6 × 3
#>   id    item  response
#>   <chr> <chr>    <dbl>
#> 1 1     I1           2
#> 2 1     I2           2
#> 3 1     I3           2
#> 4 1     I4           2
#> 5 2     I1           1
#> 6 2     I2           1

Fitting the PCM

The model is fitted once and saved to disk. The code chunk below shows the fitting call (not evaluated during R CMD check). A pre-fitted model stored at fits/fit_pcm.rds is loaded instead.

prior_pcm <- prior("normal(0, 5)", class = "Intercept") +
  prior("normal(0, 3)", class = "sd", group = "id")
fit_pcm <- brm(
  response | thres(gr = item) ~ 1 + (1 | id),
  data    = df_pcm,
  family  = acat,
  prior   = prior_pcm,
  chains  = 4,
  cores   = 4,
  iter    = 2000
)
saveRDS(fit_pcm, "fits/fit_pcm.rds")
fit_pcm <- readRDS("fits/fit_pcm.rds")

Item Fit: Infit and Outfit Statistics

infit_statistic() computes posterior predictive infit and outfit statistics for each item. Values near 1.0 indicate good fit; values substantially above 1 suggest underfit (unexpected responses), values below 1 suggest overfit (too predictable).

The ndraws_use argument limits the number of posterior draws used, which speeds up computation during exploration. For final reporting, use all draws (set ndraws_use = NULL or omit it).

fit_stats <- infit_statistic(fit_pcm, ndraws_use = 500)

# Posterior predictive p-values summarised per item
fit_stats %>%
  group_by(item) %>%
  summarise(
    infit_obs = round(mean(infit),3),
    infit_rep = round(mean(infit_rep),3),
    infit_ppp = round(mean(infit_rep > infit),3)
  )
#> # A tibble: 4 × 4
#>   item  infit_obs infit_rep infit_ppp
#>   <chr>     <dbl>     <dbl>     <dbl>
#> 1 I1        1.04      0.998     0.262
#> 2 I2        1.08      0.998     0.11 
#> 3 I3        0.922     0.996     0.856
#> 4 I4        1.03      0.995     0.312

infit_obs indicates the observed conditional infit, which can be compared to infit_rep, which is akin to the model expected value. Posterior predictive p-values (*_ppp) close to 0.5 indicate that the observed statistic falls near the centre of the posterior predictive distribution, implying good fit. Values near 0 or 1 warrant further investigation.

Item–Rest Score Association

item_restscore_statistic() computes Goodman-Kruskal’s gamma between each item’s observed responses and the rest score (total score minus the focal item). In a well-fitting Rasch model, gamma should be positive and moderate; very high values may indicate redundancy, very low or negative values suggest the item does not relate well to the latent trait.

rest_stats <- item_restscore_statistic(fit_pcm, ndraws_use = 500)

rest_stats[,1:5] %>% 
  mutate(across(where(is.numeric), ~ round(.x, 3)))
#> # A tibble: 4 × 5
#>   item  gamma_obs gamma_rep gamma_diff   ppp
#>   <chr>     <dbl>     <dbl>      <dbl> <dbl>
#> 1 I3        0.643     0.532      0.11  0.968
#> 2 I4        0.535     0.542     -0.007 0.466
#> 3 I1        0.473     0.543     -0.07  0.1  
#> 4 I2        0.441     0.548     -0.107 0.032

The output is limited to columns 1 through 5 in the output above.

Dimensionality: Residual PCA

plot_residual_pca() performs a principal components analysis on the person-item residuals and plots the loadings on the first contrast together with the item locations and the uncertainty of both. Substantial loadings on the first contrast suggest multidimensionality.

pca <- plot_residual_pca(fit_pcm, ndraws_use = 500)
pca$plot
plot of chunk pca-plot

plot of chunk pca-plot

Items with positive loadings cluster on one end, negative loadings on the other. If the observed largest eigenvalue is smaller than the replicated, unidimensionality is supported. The ppp should not be close to 1.

Local Dependence: Q3 Residual Correlations

q3_statistic() computes Yen’s Q3 statistic — the correlation between person-item residuals for every item pair. After conditioning on the latent trait, residuals should be uncorrelated; elevated Q3 values indicate local dependence (LD). Our primary metric here is the ppp, that should not be close to 1. The output is filtered on ppp values above 0.95.

q3_stats <- q3_statistic(fit_pcm, ndraws_use = 500)

q3_stats %>% 
  filter(ppp > 0.95)
#> # A tibble: 1 × 14
#>   item_1 item_2 q3_obs  q3_rep q3_diff   ppp q3_obs_q025 q3_obs_q975
#>   <chr>  <chr>   <dbl>   <dbl>   <dbl> <dbl>       <dbl>       <dbl>
#> 1 I3     I4      0.342 0.00208   0.340     1       0.285       0.393
#> # ℹ 6 more variables: q3_diff_q025 <dbl>, q3_diff_q975 <dbl>,
#> #   q3_obs_q005 <dbl>, q3_obs_q995 <dbl>, q3_diff_q005 <dbl>,
#> #   q3_diff_q995 <dbl>

Pairs flagged as locally dependent should be examined for substantive overlap in item content.

Item Category Probabilities

This plot shows the probability of using a response category on the y axis and the latent score on the x axis. The crossover points, where lines meet, are the item category threshold locations. Uncertainty is shown with the shaded area around each line.

plot_ipf(fit_pcm, theta_range = c(-6,5))
plot of chunk ipf-plot

plot of chunk ipf-plot

Person–Item Targeting

plot_targeting() produces a Wright map (person–item targeting plot) showing the distribution of person locations alongside the item threshold locations on the same logit scale. Good targeting occurs when person and item distributions overlap substantially.

plot_targeting(fit_pcm)
plot of chunk targeting

plot of chunk targeting

If the person distribution is systematically above or below the item thresholds, the test may be too easy or too hard for the sample.

Reliability: Relative Measurement Uncertainty

RMUreliability() provides a Bayesian reliability estimate via Relative Measurement Uncertainty (RMU, see Bignardi et al., 2025). It requires a matrix of person location draws with dimensions \[persons × draws\]. The output is a point estimate and lower/upper 95% highest density continuous intervals (HDCI).

person_draws <- fit_pcm %>%
  as_draws_df() %>%
  as_tibble() %>% 
  select(starts_with("r_id")) %>%
  t()
rmu <- RMUreliability(person_draws)
rmu
#>   rmu_estimate hdci_lowerbound hdci_upperbound .width .point .interval
#> 1    0.6728529       0.6132611       0.7265942   0.95   mean      hdci

RMU values range from 0 to 1, with higher values indicating higher reliability, similarly to traditional reliability metrics.