--- title: "Planting misfit and watching the diagnostics fire" author: "rasch" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Planting misfit and watching the diagnostics fire} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.2, dpi = 96 ) options(digits = 4) ``` ```{r library} library(rasch) ``` ## 1. Why simulate A diagnostic is only trustworthy if it fires when it should and stays quiet when it should not. The surest way to check that is to *manufacture* the fault it is meant to detect and confirm that the statistic responds. The package ships a family of simulators (`simulate_rasch`, `simulate_btl`, `simulate_mfrm`, `simulate_efrm`) built around exactly this plant-and-detect loop: each generates data from the appropriate member of the Rasch family and lets a *known* departure from the model be dialled in, so the matching diagnostic can be watched as it responds. Every simulator returns data ready for its fit function, with the generating parameters and the planted departures carried on the object as `attr(x, "truth")`. The print method reports what was planted, so a dataset never becomes anonymous once it leaves the call that made it. ```{r whatplanted} d <- simulate_rasch(n_persons = 400, n_items = 10, seed = 101) d # the print method reports the plant names(attr(d, "truth")) # the truth travels with the data ``` With no departures requested the data conform to the model, and the print method says so. The rest of this vignette plants one fault at a time and runs the diagnostic built to find it. The diagnostics themselves follow the conventions of Andrich and Marais (2019); the simulators exist to test them. ## 2. A clean baseline Before planting anything, it is worth seeing what recovery looks like when the model holds exactly, so that later departures are read against an honest null. `sim_recovery` compares the parameters a fit recovers against the ones the simulator planted (mean-centring the locations, which the model identifies only up to an origin), and `plot_recovery` draws each true-versus-estimated panel with its correlation and root-mean-square error. ```{r baseline} fit <- rasch(d, id = "id") rec <- sim_recovery(fit, d) rec ``` ```{r baseline-plot, fig.alt = "Recovery scatter plots of planted against recovered item difficulty and person ability."} plot_recovery(rec) ``` The two panels tell different stories, and the difference is not a defect of the estimator but the arithmetic of information. Each item difficulty pools evidence from all 400 persons, so it is recovered tightly: the correlation is near unity and the error is around a tenth of a logit. Each person ability, by contrast, rests on only the ten items that person answered, so its weighted likelihood estimate (Warm 1989) carries a standard error near seven tenths of a logit; the recovery correlation is correspondingly lower and the root-mean-square error much larger. The bias is negligible in both. A short test measures persons imprecisely *even when the model is exactly true* --- a point worth holding onto before reading any person-level result from a ten-item scale. ## 3. Planting item misfit The `discrimination` argument scales an item's slope. A value above one makes an item over-discriminate --- its responses become more deterministic than the model expects, a Guttman-like pattern --- while a value below one makes it under-discriminate, adding noise. The fit residuals should separate both from the well-behaved items. Here one central item is given a slope of 2.5 and another a slope of 0.4. ```{r discrim} disc <- rep(1, 10); disc[5] <- 2.5; disc[6] <- 0.4 d2 <- simulate_rasch(400, 10, discrimination = disc, seed = 21) fit2 <- rasch(d2, id = "id") fit2$items[, c("item", "location", "infit_ms", "outfit_ms")] ``` The two planted items depart from the pack in opposite directions, as the log-of-mean-square fit statistics (Andrich and Marais 2019, ch. 23) require. The over-discriminating item (I05) produces standardised residuals that are *too small* --- responses more predictable than the model allows --- so its infit and outfit mean squares fall below one. The under-discriminating item (I06) produces residuals that are too large, so its mean squares rise above one. The remaining items sit near the expected value of one. With a 400-by-10 draw these mean squares are themselves estimates; what matters is the direction of departure and that the planted items separate from the rest. ## 4. Planting DIF Differential item functioning is planted with `dif`, naming the affected items and the size of the group shift, with `n_groups = 2` so there is a group to differ. Here item I06 is given a uniform shift of one logit for the second group; the other nine items are invariant. `dif_anova` analyses the residuals by group over automatically sized trait class intervals (Hagquist and Andrich 2017), and should flag I06 alone. ```{r dif} d3 <- simulate_rasch(500, 10, dif = list(items = "I06", uniform = 1), n_groups = 2, seed = 303) fit3 <- rasch(d3, id = "id", factors = "group") da <- dif_anova(fit3) da$summary[, c("item", "term", "F_uniform", "p_uniform_adj", "uniform_DIF")] ``` Only I06 carries a significant group term after the Benjamini-Hochberg adjustment across items; the other nine are quiet, as invariant items should be. The analysis of variance answers *whether* there is DIF. To ask *how much*, and whether it matters on the measurement scale, `dif_size` resolves the item into one copy per group, refits, and reports the distance between the resolved locations in logits (Andrich and Marais 2019, ch. 16). ```{r difsize} dif_size(fit3, "I06", by = "group") ``` The resolved gap recovers the planted shift of one logit and clears the half-logit practical criterion. Flagging and magnitude are complementary: a significant `dif_anova` term says the item is not invariant; the `dif_size` logit difference says whether the non-invariance is large enough to act on. ## 5. Planting local dependence Response dependence --- one item's answer partly following another's --- is planted with `dependence`, naming the item pairs and a strength. It should raise the residual correlation of the planted pair above the small negative value expected under local independence. `residual_correlations` returns Yen's (1984) Q3 and, following Christensen, Makransky and Horton (2017), the adjusted Q3\* (each Q3 less the average off-diagonal value), flagging pairs whose excess passes a threshold. ```{r dependence} d4 <- simulate_rasch(500, 10, dependence = list(pairs = list(c("I04", "I05")), strength = 1.8), seed = 41) fit4 <- rasch(d4, id = "id") rc <- residual_correlations(fit4) rc$average # near -1/(L-1) under independence head(rc$pairs, 3) # the planted pair leads the table rc$flagged ``` The planted I04--I05 pair sits at the top of the sorted table and is the only one flagged; the average off-diagonal value is close to the $-1/(L-1)$ expected when the items are locally independent. As with DIF, flagging is one thing and magnitude another. `dependence_magnitude` puts the effect on the logit scale by the resolution method of Andrich and Kreiner (2010): the dependent item is split by the category of the item it follows, refitted, and the threshold displacement read off. ```{r depmag} dependence_magnitude(fit4, dependent = "I05", independent = "I04") ``` The estimated dependence is around three quarters of a logit and highly significant --- a substantial displacement of I05's threshold by the response to I04, not a marginal correlation. (The residual principal components of the same fit, via `plot_scree`, are the complementary tool for a *second dimension* rather than a dependent pair.) ## 6. Paired comparisons The paired-comparison simulator plants faults of its own. `erratic_judges` sets a proportion of judges who choose at random; their disorder should surface in the judge fit residuals and in the transitivity of their choices. Here two of eight judges are erratic. ```{r btl} b <- simulate_btl(n_objects = 7, n_judges = 8, reps_per_pair = 30, erratic_judges = 0.25, seed = 61) bt <- btl(b, "object_a", "object_b", winner = "winner", judge = "judge") bt$judges[order(-bt$judges$fit_resid), ] ``` The two erratic judges (J1 and J2) carry the largest infit and outfit mean squares and the largest positive fit residuals; the well-behaved judges sit below one. The judge fit residual is the paired-comparison counterpart of person fit, and here it is the sharpest instrument. `btl_transitivity` asks the single-dimension question in its native paired-comparison form. A Bradley-Terry-Luce scale implies that preferences stack into one consistent order: if A beats B and B beats C then A should beat C. A *circular triad* (A beats B, B beats C, C beats A) is a local contradiction, and their rate is compared with the one quarter expected from pure guessing (Kendall and Babington Smith 1940). ```{r transitivity} tr <- btl_transitivity(bt) tr$summary[, c("n_objects", "n_triples", "n_circular", "circular_rate", "consistency")] tr$judges[, c("judge", "n_triples", "circular_rate", "consistency")] ``` Pooled across all judges the object scale is nearly transitive --- majority verdicts stack into essentially one order --- so the overall consistency is high. The per-judge index separates the individuals, with the erratic judges falling toward the bottom of the table. With only a handful of complete triples per judge this index is noisier than the fit residual, but the two agree in direction. Where the concern is not judge noise but a second attribute steering some contests, `btl_dimensionality` provides the residual-"swirl" analogue of residual principal components. ## 7. A small Monte Carlo power estimate The plant-and-detect loop scales up. `sim_replicate` calls a simulator many times with successive seeds, returning the datasets as a list, so a diagnostic can be run across all of them and its behaviour summarised. Asking how often a planted departure is flagged turns that into a detection-power estimate. Here a uniform DIF of 0.8 logits is planted on one item across ten datasets, and each is fitted and passed to `dif_anova`. ```{r power} batch <- sim_replicate(simulate_rasch, 10, n_persons = 400, n_items = 8, dif = list(items = "I04", uniform = 0.8), n_groups = 2, seed = 700) flagged <- vapply(batch, function(dd) { s <- dif_anova(rasch(dd, id = "id", factors = "group"))$summary isTRUE(s$uniform_DIF[s$item == "I04"]) }, logical(1)) mean(flagged) # proportion of runs that flagged I04 ``` The estimate is the proportion of the ten runs in which the planted item was flagged after adjustment. A shift of 0.8 logits split over two groups of 200 is detected in a majority of runs but by no means all: a real departure of moderate size against a modest sample is caught only some of the time. Ten replicates suffice to demonstrate the loop, not to pin the number down --- the Monte Carlo error on ten draws is large --- but the same six lines, run with a few hundred replicates and swept over sample size and effect, are a complete power study. ## References Andrich, D., and Kreiner, S. (2010). Quantifying response dependence between two dichotomous items using the Rasch model. *Applied Psychological Measurement*, 34, 181--192. Andrich, D., and Marais, I. (2019). *A Course in Rasch Measurement Theory*. Springer. Christensen, K. B., Makransky, G., and Horton, M. (2017). Critical values for Yen's Q3. *Applied Psychological Measurement*, 41, 178--194. Hagquist, C., and Andrich, D. (2017). Recent advances in analysis of differential item functioning in health research using the Rasch model. *Health and Quality of Life Outcomes*, 15, 181. Kendall, M. G., and Babington Smith, B. (1940). On the method of paired comparisons. *Biometrika*, 31, 324--345. Warm, T. A. (1989). Weighted likelihood estimation of ability in item response theory. *Psychometrika*, 54, 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, 125--145.