--- title: "A Rasch analysis workflow" author: "Josh McGrane" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{A Rasch analysis workflow} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.2) options(digits = 4) ``` ```{r library} library(rasch) ``` This vignette follows a Rasch analysis from the overall summary to item and person fit, targeting, dependence, and differential item functioning (DIF). The order matters. A fit statistic is difficult to interpret without knowing whether the scale separates the sample or supplies information over the relevant part of the latent trait. The same analyses are available in the Shiny application. Run `rasch::run_app()` to import data, fit a model, inspect the results, and export tables or plots without writing R analysis code. ## Fit the model required by the scoring structure `rasch` fits the partial credit model by default (Rasch 1960; Andrich and Marais 2019). Dichotomous data are its one-threshold special case. Set `model = "RSM"` only when the same category threshold structure is intended to hold across items. For person $n$, item $i$, and score $x=0,\ldots,m_i$, the partial credit model is $$ P(X_{ni}=x)= \frac{\exp\left\{x\theta_n-\sum_{k=1}^{x}\delta_{ik}\right\}} {\sum_{y=0}^{m_i}\exp\left\{y\theta_n- \sum_{k=1}^{y}\delta_{ik}\right\}}. $$ The rating scale model constrains $\delta_{ik}=\beta_i+\tau_k$. Pairwise conditioning removes $\theta_n$ from the item likelihood. Person locations are subsequently estimated by Warm's weighted likelihood method (Warm 1989). The example is polytomous and contains three person groups. One item has disordered generating thresholds, one has uniform DIF, and one item pair has local response dependence. These departures make the diagnostic sequence visible without changing the commands used for observed data. Before fitting observed data, check the coding and frequency of every category. Negative scores are read as missing; valid categories begin at zero. ```{r fit} d <- simulate_rasch( n_persons = 600, n_items = 12, model = "PCM", n_categories = 4, difficulty = c(-1.5, 1.5), disordered = "I04", dependence = list(pairs = list(c("I10", "I11")), strength = 1.3), dif = list(items = "I08", uniform = 0.8), n_groups = 3, seed = 17 ) fit <- rasch(d, model = "PCM", id = "id", factors = "group") ``` ## Overall summary and reliability The summary establishes whether estimation converged and gives the overall item--trait interaction, reliability, and the package's qualitative assessment of the power of fit tests. It should be read before individual item or person results. ```{r overall-summary} fit ``` The person separation index (PSI) compares the observed variance of person locations with their mean error variance: $$ \mathrm{PSI}= \frac{\operatorname{Var}(\hat\theta)- \operatorname{mean}\{\operatorname{SE}(\hat\theta)^2\}} {\operatorname{Var}(\hat\theta)}. $$ The implementation truncates negative values at zero. `fit$psi_noext` removes persons with extreme scores and is useful when extremes inflate the observed spread. The separation ratio and number of strata express the same information on more interpretable scales. The power label is a screening judgement based on PSI, not a formal power calculation. Low reliability weakens the ordering of persons and the formation of distinct class intervals, so item--trait and DIF tests may fail to detect real departures. Power also depends on sample size, test length, targeting, category use, and trait spread. Conversely, a large sample can make a small departure statistically significant. Here the PSI is `r sprintf("%.2f", fit$psi$PSI)`, which the package classifies as `r fit$power_of_fit`; weak separation is therefore not an immediate explanation for a quiet fit test. Fit residuals, effect sizes, plots, and substantive importance remain necessary. ## Item estimates, fit, and thresholds Item locations and their standard errors should be read alongside fit residuals and the Holm-adjusted item--trait probabilities. The table shows the six items with the largest absolute residuals; the complete results are in `fit$items`. A positive residual indicates more variation than expected and a negative residual indicates responses that are more predictable than expected. The conventional $\pm2.5$ band is a screening rule rather than a separate hypothesis test. ```{r item-estimates} item_order <- order(abs(fit$items$fit_resid), decreasing = TRUE) head(fit$items[item_order, c( "item", "location", "se", "fit_resid", "p_adj" )], 6) ``` ```{r item-fit-plot, fig.alt = "Item locations plotted against item fit residuals."} plot_item_map(fit) ``` For a polytomous item, successive thresholds should normally increase on the latent scale. The threshold map shows their order directly. Item I04 has an intervening category without its own region on the scale; the category curves show the same problem in probability form. ```{r thresholds, fig.height = 5.2, fig.alt = "Estimated category thresholds for all items on the common logit scale."} plot_threshold_map(fit) ``` ```{r category-curves, fig.alt = "Category characteristic curves and observed category proportions for item I04."} plot_ccc(fit, "I04", observed = TRUE) ``` Disordering alone does not justify collapsing categories. The response-option meanings, observed use, and category curves should support any rescoring. The model must then be refitted and checked again. ## Person estimates and fit Person fit addresses the consistency of each response pattern with the fitted scale. The following table puts the largest absolute residuals first. Positive residuals indicate unexpectedly erratic patterns; negative residuals indicate patterns that are unusually predictable. ```{r person-estimates} person_order <- order(abs(fit$person$fit_resid), decreasing = TRUE, na.last = TRUE) head(fit$person[person_order, c( "id", "group", "raw", "theta", "se", "fit_resid" )], 6) ``` ```{r person-fit-plot, fig.alt = "Person locations plotted against person fit residuals."} plot_person_fit(fit) ``` An unexpected response pattern may reflect coding or data-entry errors, careless responding, a secondary trait, or a genuine but unusual person. It is not, by itself, a reason to remove the person. In this example, `r sum(abs(fit$person$fit_resid) > 2.5, na.rm = TRUE)` persons fall outside the displayed band. Fit residuals are unavailable for extreme response patterns because those patterns do not provide an interior location at which fit can be assessed. ## Targeting and information Targeting concerns the match between the person distribution and the item threshold distribution. The table reports their locations and spread, the proportions of persons beyond the threshold range, and the principal reliability indices. ```{r targeting-summary} fit$targeting ``` For a Rasch model, test information is the sum of the conditional response variances: $$ I(\theta)=\sum_i \operatorname{Var}(X_i\mid\theta), \qquad \operatorname{SE}(\hat\theta)\approx I(\theta)^{-1/2}. $$ The person--item map below places person locations and item thresholds on the same logit scale. The information curve shows where the test is most precise. ```{r targeting-map, fig.height = 5.2, fig.alt = "Person and item distributions with the test information curve on the common logit scale."} plot_pimap(fit, information = TRUE) ``` The Wright map shows the same alignment in the conventional vertical arrangement: the person distribution beside the item thresholds on one logit scale. ```{r wright, fig.height = 5.2, fig.alt = "Wright map of the person distribution and item thresholds on the common logit scale."} plot_wright(fit) ``` The optional `WrightMap` package draws the same map with greater flexibility, including several person and item panels (Torres Irribarra and Freund 2025). Polytomous thresholds are labelled `t1`, `t2`, and so on. Threshold labels are omitted for a wholly dichotomous scale; a dichotomous item in a mixed scale retains `t1`. Panels should answer a substantive question; here the person distributions are separated by the fitted group factor. ```{r wrightmap, fig.width = 8, fig.height = 5.2, fig.alt = "Wright map with one person panel per group."} if (requireNamespace("WrightMap", quietly = TRUE)) { wright_map(fit, person_panels = "group") } ``` If `WrightMap` is not installed, install it with `install.packages("WrightMap")` before running this chunk. ## Local and trait dependence Local response dependence occurs when two responses remain associated after conditioning on the latent trait. Yen's $Q3$ is the correlation between two items' standardised residuals. Because raw $Q3$ values have a negative baseline in a finite test, `q3_star` subtracts the average off-diagonal value. ```{r local-dependence} q3 <- residual_correlations(fit) q3$average head(q3$pairs[, c("item_a", "item_b", "q3", "q3_star")], 5) ``` ```{r local-dependence-plot, fig.alt = "Heatmap of adjusted residual correlations between items."} plot_resid_cor(fit) ``` There is no universal critical value for adjusted $Q3$ (Christensen, Makransky and Horton 2017). Its size, the response process, and the content of the item pair matter. When theory supports treating dependent items as one superitem, `combine_items()` refits the complete calibration and `spread_test()` compares its threshold spread with the binomial bound. The bound does not apply to a superitem containing a polytomous component. Trait dependence is examined through the residual components and by comparing person estimates from opposed item subsets (Smith 2002). For subsets $A$ and $B$, the person-level statistic is $$ t_n=\frac{\hat\theta_{nA}-\hat\theta_{nB}} {\sqrt{\operatorname{SE}(\hat\theta_{nA})^2+ \operatorname{SE}(\hat\theta_{nB})^2}}. $$ Under unidimensionality, about the nominated alpha level of these comparisons should be significant. The exact binomial interval, the residual loadings, and the score points in each subset are part of the result. ```{r trait-dependence} dimensionality <- dimensionality_test(fit) ``` ```{r trait-dependence-plot, fig.alt = "Loadings of items on the first residual component."} plot_pca(fit) ``` Here, `r sprintf("%.1f", 100 * dimensionality$prop_significant)`% of the person comparisons are significant (exact 95% interval `r sprintf("%.1f", 100 * dimensionality$ci[1])`% to `r sprintf("%.1f", 100 * dimensionality$ci[2])`%). The test does not flag trait dependence, but one opposed subset contains only `r min(dimensionality$score_points)` score points. A quiet result from a short subtest is inconclusive rather than evidence that a secondary trait is absent. ## Differential item functioning Person factors must be nominated when the model is fitted. With a person factor $G$ and trait class interval $C$, the residual model is $$ z=\mu+G+C+G\mathbin{:}C+\varepsilon. $$ The factor term tests uniform DIF; the interaction with class interval tests non-uniform DIF. `dif_anova()` fits all nominated factors together, recognises within-person factors in repeated designs, and applies Holm familywise correction over the item-by-term tests. ```{r dif} dif <- dif_anova(fit, sizes = TRUE) flagged_dif <- subset(dif$summary, uniform_DIF | nonuniform_DIF) flagged_dif[, c( "item", "term", "F_uniform", "p_uniform_adj", "eta2_uniform", "F_nonuniform", "p_nonuniform_adj", "eta2_nonuniform" )] ``` ```{r dif-plot, fig.alt = "Observed and expected item characteristic curves for item I08 by person group."} plot_icc(fit, "I08", group = "group") ``` For a significant factor with more than two levels, the follow-up should estimate the relevant differences in Rasch logits rather than rely on a generic Tukey procedure. With `sizes = TRUE`, `dif_anova()` returns the Holm-adjusted marginal pairwise comparisons for significant main effects and difference-in-differences for significant interactions. ```{r dif-follow-up} dif$posthoc[, c( "item", "contrast", "estimate", "se", "p_adj", "lower", "upper", "practical" )] ``` Statistical significance and practical magnitude answer different questions. Any split must be supported by the response process and should be applied with `resolve_dif()`, which refits the calibration and updates the item and person estimates. The revised fit then goes through the same summary, fit, targeting, dependence, and DIF sequence. ## References Andrich, D., and Marais, I. (2019). *A Course in Rasch Measurement Theory: Measuring in the Educational, Social and Health Sciences*. Springer. Christensen, K. B., Makransky, G., and Horton, M. (2017). Critical values for Yen's Q3: Identification of local dependence in the Rasch model using residual correlations. *Applied Psychological Measurement*, 41(3), 178--194. Holm, S. (1979). A simple sequentially rejective multiple test procedure. *Scandinavian Journal of Statistics*, 6(2), 65--70. Rasch, G. (1960). *Probabilistic Models for Some Intelligence and Attainment Tests*. Copenhagen: Danish Institute for Educational Research. (Expanded edition, 1980, Chicago: University of Chicago Press.) Smith, E. V. Jr. (2002). Detecting and evaluating the impact of multidimensionality using item fit statistics and principal component analysis of residuals. *Journal of Applied Measurement*, 3(2), 205--231. Torres Irribarra, D., and Freund, R. (2025). *WrightMap: IRT item-person map with ConQuest integration*. R package version 1.5. Warm, T. A. (1989). Weighted likelihood estimation of ability in item response theory. *Psychometrika*, 54(3), 427--450. Yen, W. M. (1984). Effects of local item dependence on the fit and equating performance of the three-parameter logistic model. *Applied Psychological Measurement*, 8(2), 125--145.