A Rasch analysis workflow with easyRasch2

Descriptive plots

Before fitting any model it is worth eyeballing the response distributions:

hist(rowSums(items), col = "lightblue", main = "", xlab = "Ordinal sum scores")

Figure 1. Histogram of ordinal sum scores

RMplotBar(items, ncol = 2)

Figure 2. Faceted bar chart of response distributions

RMplotStackedbar(items, show_percent = TRUE)

Figure 3. Stacked-bar response distribution

Conditional infit MSQ

It is important to note that the RMitemInfit() function uses conditional infit, which is both robust to different sample sizes and makes ZSTD unnecessary (Müller 2020). Müller also questions the usefulness of outfit, and my simulation study (Johansson 2025) reached the same conclusion. Thus, outfit is not reported.

Conditional item infit mean-square statistics flag items whose response patterns deviate from the Rasch expectation. With RMitemInfitCutoff(), per-item highest-density intervals serve as the reference instead of rule-of-thumb cutoffs (Johansson 2025). Bootstrap p-values are also available via p_value = TRUE, with family-wise error rate (FWER; the default) or false discovery rate (FDR) correction.

NOTE: All functions that use simulation-based cutoffs (except RMitemInfitMI()) have an optional p_value = TRUE for their table outputs.

infit_cut <- RMitemInfitCutoff(items, iterations = 100, parallel = FALSE,
                           seed = 3)
RMitemInfit(items, cutoff = infit_cut)

You can also get a plot summarizing simulated and observed item infit, using RMitemInfitPlot(). Since conditional infit needs complete data, sibling functions that combine multiple imputation with conditional infit are useful when you have partial missingness: RMitemInfitMI() and RMitemInfitCutoffMI().

Based on the table, items 3, 8, and 9 underfit the Rasch model, while items 2 and 4 are overfit. Item 2 is “Feeling down, depressed, or hopeless”, which is a very general item in terms of measuring depression, so this is expected.

A low item fit value, often referred to as an item being “overfit” to the Rasch model, indicates that responses may be too predictable. This is often the case for items that are very general/broad in scope in relation to the latent variable. You will often find overfitting items to also have local dependence (residual correlations) issues with other items. Overfit may be likened to having a much stronger factor loading than other items in a confirmatory factor analysis or a higher level of discrimination in an Item Response Theory model with two or more parameters.

A high item fit value, often referred to as being “underfit” to the Rasch model, can indicate several things. Often underfit is due to multidimensionality or a question that is difficult to interpret and thus has noisy response data. The latter could for instance be caused by a question that asks about two things at the same time, or is ambiguous for other reasons.

Item-restscore

Item-restscore uses Goodman-Kruskal’s \(\gamma\) (gamma) and shows the expected and observed correlation between an item and a score based on the rest of the items (Kreiner 2011). Similarly, but inverted, to item infit, a lower observed correlation value than expected indicates underfit, that the item may not belong to the dimension. A higher than expected observed value indicates an overfitting and possibly redundant item. Overfitting items will often also show issues with local dependency.

Compared to infit, item-restscore more often flags overfit items (based on experience), and less often flags underfit items (based on a simulation study (Johansson 2025)).

RMitemRestscore(items)

Similarly to infit, item-restscore found items 2 and 4 to be overfit and 8 to be underfit. It also found item 6 to be overfit. These methods are best used together.

CFA-based cutoffs for CFI / RMSEA and item loadings

RMdimCFACutoff() simulates data from a unidimensional PCM and fits a unidimensional ordinal CFA to each simulated dataset, building a parametric-bootstrap reference distribution for the model fit indices (SRMR, CFI, RMSEA) and for the per-item standardized factor loadings. It returns this simulation object; RMdimCFA() then fits the CFA to the observed data and tabulates the observed fit indices and loadings against the simulated reference. Observed values beyond the simulated cutoffs are implausible under a unidimensional data-generating process.

cfa_cut <- RMdimCFACutoff(items, iterations = 100, parallel = FALSE,
                       seed = 2)
cfa_tbl <- RMdimCFA(items, cutoff = cfa_cut)
cfa_tbl$fit

cfa_tbl$loadings

RMdimCFAPlot(cfa_cut, data = items) returns a list of two figures: $loadings (observed loadings against their simulated expected ranges) and $fit (the simulated fit-index distributions with the observed values overlaid).

The table above agrees mostly with earlier findings. The overall model fit is clearly worse than expected under unidimensionality and item 2 is overfit (higher standardized factor loading than expected), with items 3 and 8 underfit (lower loadings). Additionally, the CFA flags item 1 as overfit. We should however take the simulation-based results with a grain of salt since we are using too few iterations to get reliable results.

Residual PCA

After fitting the Rasch model, the residuals should contain no further systematic structure. The largest eigenvalue of the residual correlation matrix can be considered the headline diagnostic; a value clearly above the simulation-based cutoff suggests a secondary dimension. However, an eigenvalue below the cutoff does not by itself support unidimensionality.

pca_cut <- RMdimResidualPCACutoff(items, iterations = 100, parallel = FALSE,
                       seed = 1)
RMdimResidualPCA(items, cutoff = pca_cut)

Also of interest is the plot of item standardized loadings on the first residual contrast and item locations. This figure can be helpful to identify clusters in data, perhaps related to local dependency and/or multidimensionality.

RMdimResidualPCA(items, output = "ggplot")

Figure 4. Standardized loadings on the first residual contrast

Martin-Löf test

This is a likelihood-ratio test of unidimensionality against an a priori specified multidimensional alternative. The p-value is obtained by parametric-bootstrap (Monte Carlo) sampling under the unidimensional null. A significant p-value constitutes strong evidence that items do not all measure one common latent variable. However, while a non-significant p-value indicates that data are consistent with one common latent variable across the partition, it is not to be considered proof of unidimensionality on its own.

In the use case demonstrated below, we have hypothesized that the five psychosomatic symptoms are a separate subscale from the other items. Looking at the PCA plot above, we can see tendencies of this clustering based on the loadings on the first residual contrast factor, with items 3-5 most clearly deviating from 1, 2, 6, and 9. Using the PCA plot to determine item partitioning for the M-L test is not recommended; see ?RMdimMartinLof for details and references.

mlof <- RMdimMartinLof(items, iterations = 100,
                       partition = list(
                         c("q1","q2","q6","q9"),
                         c("q3","q4","q5","q7","q8")
                       )
)

mlof$p_value
#> [1] 0.00990099
mlof$wle_correlation
#>   subscale_a subscale_b     r ci_lower ci_upper   p_value   n
#> 1          1          2 0.717    0.676    0.754 5.994e-96 600

The p-value rejects unidimensionality across the partition, while the subscale WLE correlation indicates the dimensions are strongly related. The M-L test results can be further investigated using the diagnostic function RMdimMartinLofResiduals().

Partial gamma LD

A second perspective on LD is the partial gamma coefficient (Kreiner and Christensen 2004; Kreiner 2007) between observed item pairs, conditional on the rest-score. Note that this function evaluates both directions of LD, thus the output is two tables. We’ll restrict the output to the 6 item-pairs with largest LD deviations.

RMlocdepGamma(items, n_pairs = 6)

You can also get simulation-based thresholds for partial gamma LD, using RMlocdepGammaCutoff(), which can be used with RMlocdepGamma() and also to plot the results with RMlocdepGammaPlot()

Item pairs flagged by both \(Q_3\) and partial gamma are the strongest candidates for further inspection or possible item revision. Some argue for creating testlets by combining a locally dependent pair into a single polytomous super-item rather than eliminating one in an LD pair. I suggest looking closely at the item content before taking any action. Often, one will see very similarly worded items, where one item in an LD pair is clearly redundant. Another solution to LD can be the Graphical Loglinear Rasch Model by Kreiner and Christensen (Kreiner 2007), but that is beyond the scope of this vignette.

Andersen LR-test

RMdifLR(items, dif_var = phq9$gender, level = "threshold", output = "ggplot")

Figure 9. Andersen LR-test DIF locations by gender

The plot shows the item threshold locations estimated in each gender group with the corresponding confidence band. The global LR test is also reported in the plot caption, with a statistically significant test indicating problems with DIF between groups. However, this test is sensitive to sample size and number of items and should, like all global tests, be interpreted with care.

Partial-gamma DIF

This plot annotates each item with its partial-\(\gamma\) DIF coefficient and shows group differences within class intervals (respondents grouped by their total/latent score).

RMitemICCPlot(items, dif_var = phq9$gender, error_band = TRUE)

Figure 10. Partial Gamma DIF Conditional Item Characteristic Curves

To demonstrate the FWER-adjusted p-values, we’ll run the partial-\(\gamma\) DIF analysis with parametric bootstrap (only 100 here for rendering speed, which limits the reliability of the p-values).

difpg <- RMdifGammaCutoff(items, dif_var = phq9$gender, iterations = 100,
                          parallel = FALSE, seed = 5)
RMdifGamma(items, cutoff = difpg, dif_var = phq9$gender, p_value = TRUE)

Rasch tree (model-based recursive partitioning)

The DIF methods above compare groups that we specify in advance. The Rasch tree approach (Strobl et al. 2015) instead searches for DIF: the model is fitted to the full sample, covariates are tested for parameter instability, and whenever instability is detected the sample is split at the covariate value that maximizes it. The procedure then repeats within each subgroup, so several useful things come for free. Continuous covariates such as age need no arbitrary pre-categorization — the optimal cutpoint is estimated from the data. Interactions are handled naturally, since later splits are conditional on earlier ones (e.g. a gender split appearing only among older respondents). And if no instability is found, the tree simply does not split.

Statistically significant splits are not necessarily important splits, particularly in large samples. RMdifTree() therefore pairs the tree with an effect-size classification for every item at every split (Henninger et al. 2025): the Mantel-Haenszel odds ratio on the Delta scale developed at the Educational Testing Service (ETS) for dichotomous data, or the partial gamma coefficient for polytomous data, classified into the ETS categories A (negligible), B (slight to moderate), and C (moderate to large) (Holland and Thayer 1986; Zwick 2012). For partial gamma the B/C boundaries are the familiar 0.21 / 0.31 used by RMdifGamma(), so the two functions read on the same scale. Note that these boundaries are conventions carried over from large-scale educational testing, not values calibrated to your sample and items — in contrast to the simulation-based cutoffs used elsewhere in this package — so the A/B/C labels are best read as a rough magnitude guide rather than a calibrated test.

RMdifTree(items, covariates = phq9$gender)

Node 1 – gender: Female vs Male (n left = 426, n right = 143)

The tree splits on gender, and the effect sizes match what the partial gamma analysis showed above: q2 is classified as C and q1 and q5 as B, while the remaining items are negligible (A). With a single binary covariate the tree reduces to the familiar two-group comparison — its advantages appear when you supply several covariates at once (covariates = phq9[, c("gender", "age")]), where it tests all of them, picks split points for continuous ones, and uncovers interactions.

Three options are worth knowing about: purification = “iterative” recomputes effect sizes while excluding already-flagged items from the matching score; stability = TRUE refits the tree on resamples to report how often each covariate is selected and where the cutpoints land — a guard against overinterpreting a single sample’s tree structure; and output = “plot” draws the tree with per-node item parameters.