Package 'easyRasch2'

PHQ-9 Depression Screener (NHANES Subsample)

Description

A processed subsample of the Patient Health Questionnaire 9-item (PHQ-9) depression screener from the U.S. National Health and Nutrition Examination Survey (NHANES), September 2024 release. Six hundred respondents were drawn at random from the cycle's PHQ-9 module subject to having complete responses on all nine items, while retaining a realistic share of respondents with a sum-score of zero (n = 8) so that floor behaviour can be illustrated in a Rasch analysis.

Usage

phq9

Format

A data frame with 600 rows and 12 variables:

q1

Little interest or pleasure in doing things. Integer 0–3.

q2

Feeling down, depressed, or hopeless. Integer 0–3.

q3

Trouble falling/staying asleep, or sleeping too much. Integer 0–3.

q4

Feeling tired or having little energy. Integer 0–3.

q5

Poor appetite or overeating. Integer 0–3.

q6

Feeling bad about yourself — or that you are a failure or have let yourself or your family down. Integer 0–3.

q7

Trouble concentrating on things, such as reading the newspaper or watching television. Integer 0–3.

q8

Moving or speaking so slowly that other people could have noticed — or the opposite, being so fidgety or restless that you have been moving around a lot more than usual. Integer 0–3.

q9

Thoughts that you would be better off dead, or of hurting yourself in some way. Integer 0–3.

gender

Self-reported gender, factor with levels "Female" and "Male" (31 respondents with missing values).

age

Age in years (integer, range 15–85).

edu

Highest educational attainment, factor with levels "Elementary School", "High school", "University".

Each PHQ-9 item uses a four-point ordinal response scale, scored 0 ("Not at all"), 1 ("Several days"), 2 ("More than half the days") and 3 ("Nearly every day").

Details

The dataset is a processed subsample intended for teaching and for the package's worked example; it should not be treated as a canonical NHANES microdata file. Users wishing to validate against NCHS-published figures should download the original public-use microdata directly from the NHANES website (see Source).

Source

U.S. Centers for Disease Control and Prevention, National Center for Health Statistics. National Health and Nutrition Examination Survey, September 2024 release. https://wwwn.cdc.gov/nchs/nhanes/search/datapage.aspx?Component=Questionnaire&CycleBeginYear=2024. NHANES data are released to the public domain by the U.S. federal government (https://www.cdc.gov/nchs/policy/data-release-policy.html).

References

Kroenke, K., Spitzer, R. L., & Williams, J. B. W. (2001). The PHQ-9: Validity of a brief depression severity measure. Journal of General Internal Medicine, 16(9), 606–613. doi:10.1046/j.1525-1497.2001.016009606.x

Examples

data(phq9)
str(phq9)
summary(rowSums(phq9[, 1:9]))

---

Partial Gamma DIF Analysis

Description

Computes partial gamma coefficients for Differential Item Functioning (DIF) using iarm::partgam_DIF(). Each item is tested for association with a single categorical DIF variable, controlling for the total score.

Usage

RMdifGamma(
  data,
  dif_var,
  cutoff = NULL,
  p_value = NULL,
  correction = c("fwer", "fdr_bh", "fdr_by", "none"),
  alpha = 0.05,
  output = "kable"
)

Arguments

data	A data.frame or matrix of item responses. Items must be scored starting at 0 (non-negative integers). Missing values (NA) are allowed, but at least one complete case must exist after combining data and dif_var.
dif_var	A vector (factor or character) of the same length as nrow(data), representing the grouping variable for DIF analysis.
cutoff	Optional. Default NULL (no cutoff applied). Can be: The return value of RMdifGammaCutoff (a list with ⁠$item_cutoffs⁠): the data.frame is extracted automatically and simulation metadata is included in the kable caption. The ⁠$item_cutoffs⁠ data.frame from RMdifGammaCutoff directly: must have columns Item, gamma_low, gamma_high. When provided, adds columns Gamma_low, Gamma_high, and Flagged (logical; TRUE when the observed partial gamma falls outside the credible range) to the result.
p_value	Logical or NULL. When TRUE, adds two-sided bootstrap p-values (p_gamma, padj_gamma) comparing each item's observed partial gamma against its simulated null distribution, and flagged reflects padj_gamma < alpha instead of the credible range. The asymptotic BH-adjusted p-value and star columns from iarm::partgam_DIF() are dropped in this mode (two p-value families in one table would invite double-reading); the simulated gamma_low / gamma_high band is kept as the effect-size reference. TRUE requires the full RMdifGammaCutoff object as cutoff (it carries the simulated distributions in ⁠$results⁠). NULL, the default, means TRUE when cutoff is that full object and FALSE otherwise, so calling with no cutoff keeps the asymptotic table unchanged. Pass FALSE to flag against the interval instead. The interval makes a poor decision rule, since its width sets a family-wise error rate of 1 - width^k over all $k$ items at once (Johansson, 2026).
correction	Character. Multiplicity correction for the bootstrap p-values: "fwer" (default; Westfall-Young studentised-max step-down), "fdr_bh", "fdr_by", or "none". Ignored when p_value = FALSE.
alpha	Numeric in (0, 1). Significance level used to flag items on the corrected p-value. Default 0.05. Ignored when p_value = FALSE.
output	Character string controlling the return value. Either "kable" (default) for a formatted knitr::kable() table, or "dataframe" for the underlying data.frame.

Details

Partial gamma (Bjorner et al., 1998) measures the association between item response and an exogenous grouping variable, controlling for the total score. Values near 0 indicate no DIF. Recommended interpretive thresholds (Bjorner et al., 1998):

• No or negligible DIF: gamma within $[-0.21, 0.21]$, or gamma not significantly different from 0.

• Slight to moderate DIF: gamma within $[-0.31, 0.31]$ (and outside $[-0.21, 0.21]$), or not significantly outside $[-0.21, 0.21]$.

• Moderate to large DIF: gamma outside $[-0.31, 0.31]$, and significantly outside $[-0.21, 0.21]$.

The iarm package must be installed (it is in Suggests, not Imports). The asymptotic adjusted p-value is a Benjamini-Hochberg adjustment over the items of the p-values from iarm::partgam_DIF(), applied by easyRasch2. iarm's own adjusted column is a Bonferroni correction whatever method is named, because iarm adjusts one p-value at a time, so it is not used.

Bootstrap p-values. When p_value = TRUE, each item's observed partial gamma is compared against its simulated null distribution (from cutoff$results, simulated with no DIF by construction). The per-item statistic is the residual studentised by the bootstrap mean and SD; the marginal p-value is the two-sided Monte-Carlo p-value ⁠(1 + #\{|t*| >= |t|\}) / (B + 1)⁠, so it can be no smaller than 1 / (B + 1). The test is two-sided because DIF has no privileged direction: an item can favour either group. This differs from RMlocdepGamma(), which tests one-sided for excess positive local dependence. correction = "fwer" uses the Westfall-Young studentised-max step-down, which exploits the bootstrap dependence among items (Ferreira, 2024); it is mildly liberal below 400 iterations in RMdifGammaCutoff() (a message is issued below that). Unlike the asymptotic p-values from iarm::partgam_DIF(), these are calibrated against the simulated Rasch null rather than the asymptotic SE; they are model-conditional and sample-size-sensitive, and are reported alongside the simulated effect-size band, not in place of it. The observed partial gamma they test is computed by the same internal routine as the simulated values, which reproduces iarm::partgam_DIF() exactly.

Value

• If output = "kable": a knitr_kable object with columns "Item", "Partial gamma", "SE", "Lower CI", "Upper CI", "Adj. p-value (BH)", and "p-value sign." (stars for the BH-adjusted p-value, as iarm shows them). When cutoff is provided, additional columns "Gamma low", "Gamma high", and "Flagged" are included. With p_value = TRUE, the asymptotic p-value columns are replaced by bootstrap "p" and "p (adj)".

• If output = "dataframe": a data.frame with columns Item, gamma, se, lower, upper, padj_bh, Significance. When cutoff is provided, columns gamma_low, gamma_high, and flagged are also included. With p_value = TRUE, padj_bh and Significance are replaced by p_gamma and padj_gamma.

Multiple comparisons

The marginal p-value controls the error rate of a single comparison: for one item (or item pair) decided on in advance it is the relevant value. But scanning all k comparisons and flagging whichever fall below alpha tests k hypotheses at once, so the chance of at least one false flag inflates to roughly $1 - (1 - \alpha)^k$ (e.g. about 34% for k = 8 at alpha = 0.05) – even when every marginal p-value is correctly calibrated. The corrected (adjusted) p-value controls this: correction = "fwer" bounds the probability of any false flag (strict, lower power), while "fdr_bh" / "fdr_by" bound the expected proportion of false flags among those raised (a more lenient middle ground). Rule of thumb: use the marginal p-value for a single pre-specified comparison, and a corrected p-value when screening the whole table – the usual workflow.

References

Bjorner, J. B., Kreiner, S., Ware, J. E., Damsgaard, M. T., & Bech, P. (1998). Differential item functioning in the Danish translation of the SF-36. Journal of Clinical Epidemiology, 51(11), 1189–1202. doi:10.1016/S0895-4356(98)00111-5

Ferreira, J. A. (2024). Methods of testing a 'small' or 'moderate' number of hypotheses simultaneously. Journal of Statistical Theory and Practice, 19(6). doi:10.1007/s42519-024-00412-4

Westfall, P. H., & Young, S. S. (1993). Resampling-Based Multiple Testing. Wiley.

Johansson, M. (2026). Simulation-based cutoffs for conditional item fit in Rasch models: Iterations, multiplicity correction, and decision stability. PsyArXiv. doi:10.31234/osf.io/7pqz4_v2

See Also

RMdifGammaCutoff

Examples

if (requireNamespace("iarm", quietly = TRUE)) {
  set.seed(42)
  sim_data <- as.data.frame(
    matrix(sample(0:1, 200 * 10, replace = TRUE), nrow = 200, ncol = 10)
  )
  colnames(sim_data) <- paste0("Item", 1:10)
  dif_group <- factor(sample(c("A", "B"), 200, replace = TRUE))

  # Default kable output
  RMdifGamma(sim_data, dif_group)

  # Return as data.frame
  RMdifGamma(sim_data, dif_group, output = "dataframe")

  # Simulation-based cutoffs (100 Monte-Carlo iterations)
  if (requireNamespace("ggdist", quietly = TRUE)) {
    cutoff_res <- RMdifGammaCutoff(sim_data, dif_var = dif_group,
                                   iterations = 100, parallel = FALSE,
                                   seed = 42)
    RMdifGamma(sim_data, dif_group, cutoff = cutoff_res)

    # The full cutoff object switches on bootstrap p-values with the
    # family-wise (Westfall-Young) correction; 100 iterations keeps the
    # example quick, use the default 400 or more in real analyses
    RMdifGamma(sim_data, dif_group, cutoff = cutoff_res,
               output = "dataframe")
  }
}

---

Simulation-Based Partial Gamma DIF Cutoff Determination

Description

Simulates the distribution of the partial gamma DIF coefficient for each item when there is no DIF, to give empirical critical values and the simulated null behind the bootstrap p-values in RMdifGamma. Every simulated dataset keeps each respondent's observed group membership and total score, so the null reflects the observed group sizes and any difference between the groups' latent distributions.

Usage

RMdifGammaCutoff(
  data,
  dif_var,
  iterations = 400,
  parallel = TRUE,
  n_cores = NULL,
  verbose = FALSE,
  seed = NULL,
  cutoff_method = "hdci",
  hdci_width = 0.95,
  dgp = c("conditional", "permutation")
)

Arguments

data	A data.frame or matrix of item responses. Items must be scored starting at 0 (non-negative integers). Only complete cases (rows without any NA) are used.
dif_var	A vector (factor, character, or integer) defining group membership for DIF analysis. Must have the same length as nrow(data). Group membership is held fixed in every simulated dataset, so the group sizes and each group's distribution of total scores are those observed.
iterations	Integer. Number of simulation iterations (default 400).
parallel	Logical. Use parallel processing via mirai if available (default TRUE).
n_cores	Integer or NULL. Number of parallel workers. When NULL, getOption("mc.cores") is checked first. If neither is set and parallel = TRUE, a warning is issued and execution falls back to sequential (single core) processing.
verbose	Logical. Show a progress bar (default FALSE).
seed	Integer or NULL. Random seed for reproducibility. See easyRasch2-reproducibility for what this guarantees and how it interacts with parallel.
cutoff_method	Character string specifying how cutoff intervals are computed. Either "hdci" (default) for the Highest Density Interval via ggdist::hdci(), or "quantile" for the 2.5th/97.5th percentiles via stats::quantile().
hdci_width	Numeric. Width of the HDCI when cutoff_method = "hdci". Default is 0.95 (95\ Ignored when cutoff_method = "quantile".
dgp	Character. How the null datasets are generated. "conditional" (default) draws each respondent's response pattern from the Rasch conditional distribution given their observed total score, with the CML item thresholds fixed, as the conditional option of RMitemRestscoreCutoff does. "permutation" keeps the observed responses and permutes group membership among respondents with the same total score, the Monte Carlo form of an exact conditional test of item-group independence given the score (Kreiner, 1987). It needs no item parameters, keeps any misfit elsewhere in the observed responses, and is 25 to 60 times faster per iteration. Both held the nominal error rate in simulation, including when the groups differed by 1 SD.

Details

Partial gamma conditions on the total score, and under the Rasch model the response to an item is independent of group membership given the total score. Both data-generating processes preserve each respondent's group and total score, so the simulated null has the observed distribution of groups across score strata. An earlier version assigned simulated respondents to groups at random, which gave both groups the same latent distribution. When a minority group differed from the majority by 1 SD, that null was too narrow and flagged at least one item in about 14 percent of datasets with no DIF, against a nominal 5 percent.

For each iteration the function draws one null dataset (see dgp) and computes partial gamma for every item. The computation reproduces iarm::partgam_DIF() exactly but is vectorised, so iarm is not needed here. Iterations that fail are discarded and reported through actual_iterations.

The conditional data-generating process uses CML item thresholds from psychotools::pcmodel() (a dichotomous item is a 2-category PCM). The group order follows the levels of dif_var (alphabetical for character vectors), which sets the sign of gamma as in iarm::partgam_DIF().

Parallel processing is provided by the mirai package (optional). Install it with install.packages("mirai") to enable parallelisation.

Value

A list with components:

results

data.frame with columns iteration, Item, and gamma (one row per item per successful iteration).

item_cutoffs

data.frame with per-item cutoff summaries: Item, gamma_low, gamma_high. Bounds are computed using the method specified by cutoff_method.

actual_iterations

Number of successful iterations.

requested_iterations

Number of iterations asked for.

sample_n

Number of complete cases used.

sample_n_total

Number of respondents in the raw input data, before removing rows with NA in data or dif_var.

sample_has_na

Logical. Whether data or dif_var contained any missing values.

sample_summary

Summary statistics of the WLE person locations, or NULL when the model could not be fitted with dgp = "permutation".

item_names

Character vector of item names from data.

dif_group_sizes

Integer vector of group sizes, held fixed in every simulated dataset.

cutoff_method

The method used to compute cutoffs ("hdci" or "quantile").

hdci_width

The HDCI width used (only meaningful when cutoff_method = "hdci").

dgp

The data-generating process used.

References

Bjorner, J. B., Kreiner, S., Ware, J. E., Damsgaard, M. T., & Bech, P. (1998). Differential item functioning in the Danish translation of the SF-36. Journal of Clinical Epidemiology, 51(11), 1189–1202. doi:10.1016/S0895-4356(98)00111-5

Henninger, M., Radek, J., Debelak, R., & Strobl, C. (2025). Partial credit trees meet the partial gamma coefficient for quantifying DIF and DSF in polytomous items. Behaviormetrika, 52, 221–257. doi:10.1007/s41237-024-00252-3

Kreiner, S. (1987). Analysis of multidimensional contingency tables by exact conditional tests: Techniques and strategies. Scandinavian Journal of Statistics, 14(2), 97–112.

See Also

partgam_DIF, RMdifGamma, RMdifGammaPlot

Examples

if (requireNamespace("ggdist", quietly = TRUE)) {
  set.seed(42)
  sim_data <- as.data.frame(
    matrix(sample(0:1, 200 * 10, replace = TRUE), nrow = 200, ncol = 10)
  )
  colnames(sim_data) <- paste0("Item", 1:10)
  dif_sex <- sample(c("male", "female"), 200, replace = TRUE)

  # Run 100 iterations sequentially for a quick demo
  cutoff_res <- RMdifGammaCutoff(sim_data, dif_var = dif_sex,
                                 iterations = 100, parallel = FALSE,
                                 seed = 42)
  cutoff_res$item_cutoffs

  # The permutation null needs no item parameters and is much faster
  perm_res <- RMdifGammaCutoff(sim_data, dif_var = dif_sex,
                               iterations = 100, parallel = FALSE,
                               seed = 42, dgp = "permutation")
  perm_res$item_cutoffs
}

---

Plot Distribution of Simulated Partial Gamma DIF Values

Description

Visualises the distribution of simulation-based partial gamma DIF values from RMdifGammaCutoff, optionally overlaying observed partial gamma values computed from real data with the same routine the simulation uses, which reproduces partgam_DIF exactly.

Usage

RMdifGammaPlot(simfit, data, dif_var)

Arguments

simfit	The return value of RMdifGammaCutoff (a list with components results, item_cutoffs, actual_iterations, sample_n, and item_names).
data	Optional. A data.frame or matrix of item responses for computing and overlaying observed partial gamma values. Items must be scored starting at 0 (non-negative integers). When provided, the plot includes orange diamond markers for the observed partial gamma alongside the simulated distribution, plus segment summaries from the cutoff intervals.
dif_var	Required when data is supplied. A vector (factor, character, or integer) defining group membership for the DIF analysis. Must have the same length as nrow(data).

Details

Uses ggdist::stat_dotsinterval() (when data is not supplied) or ggdist::stat_dots() (when data is supplied) with point_interval = "median_hdci" and .width = c(0.66, 0.95, 0.99).

When data is not supplied, the function plots the simulated partial gamma distributions as dot-interval plots using ggdist::stat_dotsinterval() with median and Highest Density Continuous Interval (HDCI) summaries.

When data is supplied (along with dif_var), the function:

1. Computes observed partial gamma values (identical to iarm::partgam_DIF()).

2. Overlays observed gamma values as orange diamond markers on the simulated distributions.

3. Shows per-item cutoff intervals (from simfit$item_cutoffs) as black line segments, with thicker segments for the 66\ black dots for the median.

The ggplot2 and ggdist packages must be installed (they are in Suggests, not Imports).

Value

A ggplot object.

See Also

RMdifGammaCutoff, RMdifGamma

Examples

if (requireNamespace("ggdist", quietly = TRUE) &&
    requireNamespace("ggplot2", quietly = TRUE)) {
  set.seed(42)
  sim_data <- as.data.frame(
    matrix(sample(0:1, 200 * 10, replace = TRUE), nrow = 200, ncol = 10)
  )
  colnames(sim_data) <- paste0("Item", 1:10)
  dif_group <- factor(sample(c("A", "B"), 200, replace = TRUE))

  # Run simulation
  cutoff_res <- RMdifGammaCutoff(sim_data, dif_var = dif_group,
                                 iterations = 100, parallel = FALSE,
                                 seed = 42)

  # Simulated distribution only
  RMdifGammaPlot(cutoff_res)

  # With observed partial gamma overlaid
  RMdifGammaPlot(cutoff_res, data = sim_data, dif_var = dif_group)
}

---

DIF analysis via Andersen's likelihood-ratio test

Description

Splits a Rasch model by an external grouping variable using eRm::LRtest() and reports per-group item locations (or per-group threshold locations) together with their standard errors. A single function replaces the four legacy helpers (RIdifTableLR, RIdifThreshTblLR, RIdifFigureLR, RIdifThreshFigLR) by exposing the two underlying axes – level (item or threshold) and output (data.frame, kable, or ggplot) – as arguments. The same data preparation pipeline feeds all six combinations.

Usage

RMdifLR(
  data,
  dif_var,
  model = c("auto", "PCM", "RM"),
  level = c("item", "threshold"),
  output = c("kable", "dataframe", "ggplot"),
  cutoff = 0.5,
  conf = 0.95,
  sort = FALSE
)

Arguments

data	A data.frame or matrix of item responses (non-negative integers, 0-based). One column per item, one row per person. Person IDs and grouping variables must not be included – pass the grouping variable separately via dif_var.
dif_var	Vector of length nrow(data) (factor, character, or numeric) defining the DIF grouping variable. Coerced to factor; unused levels are dropped. Rows where dif_var is NA are dropped with a message. Must result in at least 2 groups after cleaning.
model	One of "auto" (default), "PCM", or "RM". "auto" fits eRm::RM() when the data are dichotomous (max response = 1) and eRm::PCM() otherwise. "RM" errors on polytomous data.
level	One of "item" (default) or "threshold". "item" reports each item's mean threshold location per group; "threshold" reports each individual threshold per group. For dichotomous (RM) data the two views are equivalent (one threshold per item).
output	One of "kable" (default), "dataframe", or "ggplot". The data.frame view always carries a Flagged logical column; the kable view bolds flagged rows; the ggplot view shows confidence intervals.
cutoff	Numeric or NULL. Threshold (in logits) for the Flagged column: a row is flagged when MaxDiff > cutoff, where MaxDiff is the difference between the largest and smallest per-group location for that item (or threshold). Set to NULL to suppress flagging. Default 0.5.
conf	Numeric in (0, 1). Confidence level used for the ggplot error bars. Default 0.95.
sort	Logical. kable output only: sort rows by MaxDiff (descending). Default FALSE.

Details

The Partial Credit Model (PCM) is fitted by default for polytomous data and the dichotomous Rasch Model (RM) is fitted when all responses are 0/1; this can be overridden via model.

For the data.frame and kable outputs, locations are reported on the centred eRm parameterisation returned by eRm::thresholds(). Per-group fits come from eRm::LRtest(..., splitcr = dif_var); the unsplit fit (All column) is the model fitted to the full dataset. The Andersen LR statistic, df, and p-value reported as the lr_test attribute / caption come directly from LRtest()'s return value.

cell_spec()-style HTML cell colouring used in the legacy easyRasch package has been dropped in favour of a logical Flagged column (and bold rendering in the kable output), so the kable renders correctly in HTML, LaTeX, and pipe/markdown.

Value

A data.frame, a knitr_kable object, or a ggplot object, depending on output.

The data.frame has one row per item (level = "item") or per item x threshold (level = "threshold"), with columns Item (and Threshold at threshold level), one numeric column per group level, an All column for the unsplit fit, MaxDiff, Flagged (when cutoff is non-NULL), and matching SE_* columns.

The Andersen LR test result is attached as attr(result, "lr_test") on the data.frame, in the kable footnote, and in the ggplot caption (LR $\chi^2$, df, p-value).

Examples

if (requireNamespace("eRm", quietly = TRUE)) {
  set.seed(1)
  data("pcmdat2", package = "eRm")
  grp <- factor(sample(c("A", "B"), nrow(pcmdat2), replace = TRUE))

  # Default: kable of per-group item locations
  RMdifLR(pcmdat2, dif_var = grp)

  # ggplot panel of item locations with 95% CIs
  RMdifLR(pcmdat2, dif_var = grp, output = "ggplot")

  # Threshold-level kable, sorted by MaxDiff
  RMdifLR(pcmdat2, dif_var = grp, level = "threshold", sort = TRUE)

  # Tidy data.frame for downstream use
  df <- RMdifLR(pcmdat2, dif_var = grp, output = "dataframe")
  attr(df, "lr_test")
  df[df$Flagged, ]
}

---

Tree-based DIF analysis with effect-size classification

Description

Detects Differential Item Functioning (DIF) by recursively splitting the sample on one or more covariates using psychotree::raschtree() (dichotomous data) or psychotree::pctree() (polytomous data), then computes per-split effect-size measures for every item – the Mantel-Haenszel odds-ratio (Holland & Thayer, 1986), on the Delta scale developed at the Educational Testing Service (ETS; Zwick, 2012), for dichotomous data, or the partial gamma coefficient for polytomous data – and classifies them into ETS A/B/C categories (Bjorner et al., 1998). Optionally, an iterative purification step (Henninger et al., 2025) is applied, and tree stability across resamples can be assessed via stablelearner::stabletree() (Philipp et al., 2018).

Usage

RMdifTree(
  data,
  covariates,
  model = c("auto", "PCM", "RM"),
  effect_size = c("auto", "MH", "pgamma"),
  purification = c("none", "iterative"),
  p_adj = c("none", "fdr", "bonferroni"),
  thresholds = c(0.21, 0.31),
  alpha = 0.05,
  prune_negligible = FALSE,
  stability = FALSE,
  stability_B = 100L,
  stability_sampler = c("subsampling", "bootstrap"),
  min_n_per_level = 20L,
  on_rescale = c("message", "warning", "stop"),
  output = c("kable", "dataframe", "tree", "plot", "list"),
  ...
)

Arguments

data	A data.frame or matrix of item responses (non-negative integers, 0-based). One column per item. Items only – pass covariates separately via covariates.
covariates	A data.frame with nrow(data) rows giving the DIF covariates. Columns may be numeric (continuous), factor, ordered factor, or logical. At least one column is required.
model	One of "auto" (default), "PCM", or "RM". "auto" picks RM when the data are dichotomous (max response = 1) and PCM otherwise.
effect_size	One of "auto" (default), "MH" (Mantel-Haenszel), or "pgamma" (partial gamma). "auto" uses MH for RM and partial gamma for PCM.
purification	"none" (default) or "iterative". Iterative purification recomputes the effect size while excluding items already classified as DIF from the matching score (Bjorner et al., 1998).
p_adj	Multiple-testing adjustment for the partial-gamma classification: "none" (default), "fdr" (Benjamini & Hochberg), or "bonferroni". Ignored for MH (which uses the ETS test directly).
thresholds	Numeric length-2 vector with the partial-gamma B/C boundaries (default c(0.21, 0.31)). Ignored for MH (the ETS Delta-scale boundaries 1.0 / 1.5 are used).
alpha	Significance level for the A/C tests. Default 0.05.
prune_negligible	Logical. If TRUE, splits whose every item is classified as A are pruned from the tree. Default FALSE.
stability	Logical. If TRUE, runs stablelearner::stabletree() on the fitted tree to assess variable-selection and cutpoint stability across resamples. Default FALSE.
stability_B	Integer. Number of resamples for stability assessment. Default 100.
stability_sampler	One of "subsampling" (default) or "bootstrap". Subsampling is recommended for tree stability because bootstrap resamples can drop levels of categorical covariates, breaking the model fit on small subsamples.
min_n_per_level	Integer. Minimum count required for any factor level in covariates. Default 20. Set to 0 to disable.
on_rescale	One of "message" (default), "warning", or "stop". Controls how the function reports the situation in which one or more items have no responses in their lowest category within a terminal node of the tree (psychotree silently rescales those items per node, which leaves the per-node item parameters on a different metric). The MH and partial-gamma effect sizes reported here are computed on the raw responses and are not affected. The diagnostic concerns only the per-node item parameters that are visible via plot(tree).
output	One of "kable" (default), "dataframe", "tree", "plot", or "list". "kable" renders the per-split per-item effect-size table as a knitr::kable(); "dataframe" returns the underlying tidy data.frame; "tree" returns the augmented partykit tree object; "plot" returns the partykit tree plot with item names on the terminal-node x-axis (tp_args = list(names = TRUE)); "list" returns both views from a single fit, as list(table = <data.frame>, tree = <tree object>). For full control over the plot, use output = "tree" and call plot() on the result with your own tp_args. When stability = TRUE, the stability summary is attached as attr(result, "stability") (a small data.frame) and attr(result, "stability_kable") (pre-rendered kable).
...	Additional arguments forwarded to psychotree::raschtree() or psychotree::pctree(), e.g. minsize, alpha for the parameter-instability test. Note: this alpha (the tree-fitting argument) is independent of the alpha for ETS classification above; pass it via ....

Details

Continuous covariates (e.g., age in years) and interactions among multiple covariates are handled natively by the model-based recursive-partitioning machinery of partykit / psychotree.

ETS Delta-scale Mantel-Haenszel. The MH common odds ratio $\hat{\alpha}_{MH}$ is mapped to the ETS Delta scale via $\Delta_{MH} = -2.35 \log \hat{\alpha}_{MH}$. The classification rules (Holland & Thayer, 1988; Zwick, 2012) are:

• Class A: $|\Delta_{MH}| < 1$ or test of $\Delta_{MH} = 0$ not rejected at alpha.

• Class C: $|\Delta_{MH}| \ge 1.5$ and test of $|\Delta_{MH}| \le 1$ rejected at alpha.

• Class B: otherwise.

easyRasch2 uses the sign convention of effecttree: positive $\Delta$ indicates that the item is more difficult for the second (reference) group.

Partial gamma. iarm::partgam_DIF() provides the gamma estimate and SE per item. ETS-style classification follows Bjorner et al. (1998): A if $|\gamma| < 0.21$ or the test of $\gamma = 0$ is not rejected at alpha; C if $|\gamma| > 0.31$ and the test of $|\gamma| \le 0.21$ is rejected at alpha; B otherwise.

Rule-of-thumb caveat. The A/B/C boundaries (Delta 1.0 / 1.5; gamma 0.21 / 0.31) are conventions carried over from large-scale educational testing, not values calibrated to the sample and items at hand – in contrast to the simulation-based cutoffs used elsewhere in this package – so the classification is best read as a rough magnitude guide rather than a calibrated test. For a sample-calibrated partial-gamma DIF test, see RMdifGamma() with an RMdifGammaCutoff() object (optionally with p_value = TRUE).

Continuous and interaction effects. The recursive partitioning step automatically handles continuous covariates (the parameter-instability test searches for the optimal cutpoint) and multivariate interactions (later splits are conditional on earlier ones). To assess sensitivity of variable selection and cutpoints to the particular sample, set stability = TRUE.

Stability assessment. When stability = TRUE, the same tree-fitting call is replayed on stability_B resamples of the data. The returned stabletree object reports, per covariate, the proportion of resamples in which it was selected at any split, and (for continuous covariates) the empirical distribution of cutpoints. Stability is independent of effect-size classification – a stable split with a negligible effect is still negligible, and a large effect on an unstable split should be interpreted with care. Cost is roughly stability_B times the original fitting time.

Acknowledgement. The effect-size machinery – MH and partial gamma per node, iterative purification, ETS A/B/C classification – follows the implementations by Henninger & Radek in the GitHub packages raschtreeMH and effecttree. The relevant code has been adapted here under the MIT licence (see file header).

Value

Depending on output:

"kable"

A knitr_kable object with one row per (split node x item), grouped by node via pack_rows-style section headers. The caption summarises model, effect-size measure, purification, and (if requested) stability.

"dataframe"

A data.frame with one row per (split node x item): columns NodeID, Split (human-readable description of the split), Variable, Direction, Item, EffectSize, SE, Class (A/B/C), Flagged (TRUE for B or C), n_left, n_right.

"tree"

The fitted partykit tree object with class c("RMdifTree", ...), with effect-size results stored at tree$info$effectsize.

"plot"

A plotted partykit tree.

"list"

A list with both views from a single fit: table (the "dataframe" result, with its attributes) and tree (the "tree" result).

Stability results (when stability = TRUE) are attached to the return value as attr(result, "stability") (data.frame) and attr(result, "stability_kable") (pre-rendered kable).

References

Bjorner, J. B., Kreiner, S., Ware, J. E., Damsgaard, M. T., & Bech, P. (1998). Differential item functioning in the Danish translation of the SF-36. Journal of Clinical Epidemiology, 51(11), 1189-1202. doi:10.1016/S0895-4356(98)00111-5

Henninger, M., Debelak, R., & Strobl, C. (2023). A new stopping criterion for Rasch trees based on the Mantel-Haenszel effect size measure for DIF. Educational and Psychological Measurement, 83, 181-212. doi:10.1177/00131644221077135

Henninger, M., Radek, J., Debelak, R., & Strobl, C. (2025). Partial credit trees meet the partial gamma coefficient for quantifying DIF and DSF in polytomous items. Behaviormetrika, 52, 221-257. doi:10.1007/s41237-024-00252-3

Holland, P. W., & Thayer, D. T. (1986). Differential item performance and the Mantel-Haenszel procedure. ETS Research Report Series(2). doi:10.1002/j.2330-8516.1986.tb00186.x

Philipp, M., Rusch, T., Hornik, K., & Strobl, C. (2018). Measuring the stability of results from supervised statistical learning. Journal of Computational and Graphical Statistics, 27, 685-700. doi:10.1080/10618600.2018.1473779

Asamoah, N. A. B., Turner, R. C., Lo, W.-J., Crawford, B. L., & Jozkowski, K. N. (2025). Impacts of DIF Item Balance and Effect Size Incorporation With the Rasch Tree. Educational and Psychological Measurement. doi:10.1177/00131644251370605

Zwick, R. (2012). A review of ETS differential item functioning assessment procedures: Flagging rules, minimum sample size requirements, and criterion refinement. ETS Research Report Series(1). doi:10.1002/j.2333-8504.2012.tb02290.x

See Also

RMdifLR, RMdifGamma

Examples

if (requireNamespace("psychotree", quietly = TRUE) &&
    requireNamespace("partykit", quietly = TRUE) &&
    requireNamespace("difR", quietly = TRUE) &&
    requireNamespace("iarm", quietly = TRUE)) {
  data("DIFSimPC", package = "psychotree")

  items <- as.data.frame(as.matrix(DIFSimPC$resp))
  covs  <- DIFSimPC[, c("age", "gender", "motivation")]

  # Default: kable of per-split effect sizes
  RMdifTree(items, covariates = covs)

  # Tidy data.frame -- one row per (split node x item)
  df <- RMdifTree(items, covariates = covs, output = "dataframe")
  df[df$Flagged, ]

  # Tree object (for plotting via partykit::plot)
  tree <- RMdifTree(items, covariates = covs, output = "tree")
  plot(tree)

  # Stability assessment refits the tree on B resamples.
  # (use more resamples, e.g. 100+, in real analyses)
  if (requireNamespace("stablelearner", quietly = TRUE)) {
    kbl <- RMdifTree(items, covariates = covs,
                     purification = "iterative", p_adj = "fdr",
                     stability = TRUE, stability_B = 25)
    kbl                          # main effect-size kable
    attr(kbl, "stability_kable") # pre-rendered stability kable
    attr(kbl, "stability")       # raw stability data.frame
  }
}

---

Observed one-factor CFA fit and loadings vs a simulated reference

Description

Fits the observed one-factor categorical CFA to data and compares its fit indices and per-item standardized loadings against the simulated null distribution produced by RMdimCFACutoff. Returns a list of two tables: model-fit indices and per-item loadings, each with the observed value, the expected reference from the simulation, and a flag.

Usage

RMdimCFA(
  data,
  cutoff,
  p_value = FALSE,
  correction = c("fwer", "fdr_bh", "fdr_by", "none"),
  alpha = 0.05,
  output = c("kable", "dataframe")
)

Arguments

data	A data.frame or matrix of item responses (non-negative integers, 0-based), the same items used for the cutoff simulation.
cutoff	The list returned by RMdimCFACutoff. Required: observed CFA fit indices are not interpretable without the simulated reference, so the function errors if it is missing.
p_value	Logical. When TRUE, adds bootstrap p-values from the simulated null distributions: one-sided in the unfavourable direction for the fit indices (CFI low; RMSEA / SRMR high), two-sided for the per-item loadings. The Flagged columns then reflect padj < alpha instead of the percentile cutoffs. The fit indices and the loadings are corrected as two separate families. Default FALSE.
correction	Character. Multiplicity correction for the p-values: "fwer" (default; Westfall-Young studentised-max step-down), "fdr_bh" (Benjamini-Hochberg), "fdr_by" (Benjamini-Yekutieli), or "none". Ignored when p_value = FALSE.
alpha	Numeric in (0, 1). Significance level used to flag comparisons on the corrected p-value. Default 0.05. Ignored when p_value = FALSE.
output	Character. "kable" (default) returns each table as a knitr::kable(); "dataframe" returns plain data.frames.

Details

Bootstrap p-values. The per-comparison statistic is the residual studentised by the bootstrap mean and SD. Marginal p-values are Monte-Carlo, (1 + count) / (B + 1), so they can be no smaller than 1 / (B + 1). correction = "fwer" uses the Westfall-Young studentised-max step-down, which exploits the bootstrap dependence among the statistics (Ferreira, 2024); it is liberal when the simulation is small, so at least 1000 iterations in RMdimCFACutoff() are recommended (a warning is issued below that). These p-values are model-conditional and sample-size-sensitive and are reported alongside the simulated expected ranges, not in place of them.

Value

A named list with two elements, fit and loadings:

fit

CFI / RMSEA / SRMR with columns Index, Observed, Cutoff, Direction, Flagged (one-sided, in the unfavourable direction). With p_value = TRUE, columns p and padj are added and Flagged reflects padj < alpha.

loadings

One row per item with columns Item, Observed, Expected_low, Expected_high, Flagged ("below" / "above" / ""). With p_value = TRUE, columns p_loading and padj_loading are added and Flagged reflects padj_loading < alpha (direction from the sign of the deviation from the simulated mean).

Each element is a knitr_kable (when output = "kable") or a data.frame (when output = "dataframe").

Multiple comparisons

The marginal p-value controls the error rate of a single comparison: for one item (or item pair) decided on in advance it is the relevant value. But scanning all k comparisons and flagging whichever fall below alpha tests k hypotheses at once, so the chance of at least one false flag inflates to roughly $1 - (1 - \alpha)^k$ (e.g. about 34% for k = 8 at alpha = 0.05) – even when every marginal p-value is correctly calibrated. The corrected (adjusted) p-value controls this: correction = "fwer" bounds the probability of any false flag (strict, lower power), while "fdr_bh" / "fdr_by" bound the expected proportion of false flags among those raised (a more lenient middle ground). Rule of thumb: use the marginal p-value for a single pre-specified comparison, and a corrected p-value when screening the whole table – the usual workflow.

References

Ferreira, J. A. (2024). Methods of testing a 'small' or 'moderate' number of hypotheses simultaneously. Journal of Statistical Theory and Practice, 19(6). doi:10.1007/s42519-024-00412-4

Westfall, P. H., & Young, S. S. (1993). Resampling-Based Multiple Testing. Wiley.

See Also

RMdimCFACutoff, RMdimCFAPlot

Examples

if (requireNamespace("lavaan", quietly = TRUE) &&
    requireNamespace("eRm", quietly = TRUE)) {
  data("raschdat1", package = "eRm")
  sim <- RMdimCFACutoff(raschdat1[, 1:8], iterations = 50,
                        parallel = FALSE, seed = 1)
  tabs <- RMdimCFA(raschdat1[, 1:8], cutoff = sim)
  tabs$fit
  tabs$loadings

  # Bootstrap p-values with family-wise (Westfall-Young) correction
  # (use iterations >= 1000 in real analyses for stable p-values)
  RMdimCFA(raschdat1[, 1:8], cutoff = sim, p_value = TRUE)
}

---

Simulated null distribution for one-factor CFA fit and loadings under PCM unidimensionality

Description

Generates a parametric-bootstrap null distribution against which observed one-factor categorical-CFA results can be compared. The simulation draws iterations datasets from the fitted PCM (or RM, for dichotomous data) using the observed item parameters and a resampled person distribution; each simulated dataset is fitted with lavaan::cfa(..., ordered = TRUE, estimator = "WLSMV") and both the three fit indices (CFI, RMSEA, SRMR) and the per-item standardized factor loadings are recorded. Because the simulated data satisfy the PCM unidimensional assumption by construction, the resulting distributions are the "expected" reference for what a correctly fitting unidimensional model produces at this sample size and item structure.

Usage

RMdimCFACutoff(
  data,
  iterations = 250L,
  percentile = 99,
  output = c("list", "kable"),
  parallel = TRUE,
  n_cores = NULL,
  verbose = FALSE,
  seed = NULL,
  estimator = "WLSMV"
)

Arguments

data	A data.frame or matrix of item responses (non-negative integers, 0-based). One column per item, one row per person.
iterations	Integer. Number of parametric-bootstrap iterations. Default 250.
percentile	Numeric in (50, 100). The strictness of the cutoffs. Default 99. Fit indices use a one-sided cutoff in the unfavourable direction (CFI from below; RMSEA and SRMR from above); standardized loadings use a two-sided central interval covering percentile% of the simulated distribution (an item is flagged when its observed loading falls in the outer 100 - percentile%, split across the two tails).
output	Character. Only "list" (the default) is supported and the function always returns the simulation object. "kable" is retained only to raise an informative error: tables now come from RMdimCFA.
parallel	Logical. If TRUE (default), uses parallel processing via mirai. Falls back to sequential if mirai is not installed or n_cores cannot be resolved.
n_cores	Integer or NULL. Number of parallel workers. When NULL, getOption("mc.cores") is consulted; if neither is set, sequential is used.
verbose	Logical. Show a progress bar (default FALSE).
seed	Integer or NULL. Master seed for reproducibility. See easyRasch2-reproducibility for what this guarantees and how it interacts with parallel.
estimator	Character. The lavaan estimator passed to lavaan::cfa(). Default "WLSMV". Other limited-information estimators that produce robust/scaled fit indices (e.g., "DWLS", "ULSMV") are also accepted; full-information ML is rejected (incompatible with ordered = TRUE).

Details

This function only generates the simulated reference. To obtain the observed-vs-expected tables, pass its result to RMdimCFA; for the figures, pass it to RMdimCFAPlot.

Generative model. The data-generating process for each simulated dataset is the PCM (or RM) fitted to the observed data, with persons drawn from the empirical theta distribution (resampled with replacement). This means the simulated data perfectly satisfy the PCM unidimensional assumption.

Estimation model. The CFA on each simulated dataset uses a single-factor model with all items as ordinal indicators (F1 =~ I1 + I2 + ...), fitted with WLSMV by default. Reported CFI / RMSEA are the Satorra-Bentler-scaled variants (cfi.scaled, rmsea.scaled) for consistency across iterations; SRMR is reported unchanged. Standardized loadings are the est.std of the ⁠=~⁠ paths from lavaan::standardizedSolution().

Why a null distribution. A perfectly PCM-unidimensional dataset will typically not yield CFA fit indices at their ideal values (CFI = 1, RMSEA = 0), nor identical loadings across items: PCM uses a logistic threshold structure while WLSMV uses a probit link via the polychoric correlation matrix, and finite samples add sampling variability. The simulated distributions capture both, giving a more honest reference than rule-of-thumb cutoffs derived under continuous-data ML.

Iteration failures. Some simulated datasets cause WLSMV to fail (non-positive-definite polychoric matrix, boundary thresholds, empty categories). Failed iterations are dropped; actual_iterations reflects the number that succeeded.

Value

A list (the simulation object), with components:

simulated

data.frame with one row per successful iteration and columns iteration, cfi, rmsea, srmr.

simulated_loadings

data.frame with one row per successful iteration: an iteration column followed by one column per item holding the simulated standardized loading.

percentile

Numeric: the strictness setting used.

cutoffs

Named numeric vector (cfi, rmsea, srmr) of one-sided fit-index cutoffs at the chosen percentile.

loading_cutoffs

data.frame Item, low, high — the two-sided expected loading interval per item.

actual_iterations

Number of successful MC iterations.

sample_n

Number of complete cases used.

sample_n_total

Number of respondents in the raw input data, before the complete-case filter.

sample_has_na

Logical. Whether the raw input data contained any missing values.

n_items

Number of items.

item_names

Character vector of item names.

is_polytomous

Logical: was a PCM (vs RM) fitted?

estimator

The lavaan estimator used.

References

Yuan, K.-H., & Bentler, P. M. (2000). Three likelihood-based methods for mean and covariance structure analysis with nonnormal missing data. Sociological Methodology, 30(1), 165-200. doi:10.1111/0081-1750.00078

Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling. Journal of Statistical Software, 48(2), 1-36. doi:10.18637/jss.v048.i02

See Also

RMdimCFA, RMdimCFAPlot, RMdimResidualPCA, RMdimMartinLof

Examples

if (requireNamespace("lavaan", quietly = TRUE) &&
    requireNamespace("eRm", quietly = TRUE)) {
  data("raschdat1", package = "eRm")

  # Few iterations for a fast example; use 250+ in real analyses
  sim <- RMdimCFACutoff(raschdat1[, 1:8], iterations = 50,
                        parallel = FALSE, seed = 1)

  # Observed-vs-expected tables
  RMdimCFA(raschdat1[, 1:8], cutoff = sim)

  if (requireNamespace("ggplot2", quietly = TRUE)) {
    plots <- RMdimCFAPlot(sim, data = raschdat1[, 1:8])
    plots$loadings
    plots$fit
  }
}

---

Plot observed CFA fit and loadings against the simulated null

Description

Returns two figures (in a list) comparing the observed one-factor CFA to the simulated null distribution from RMdimCFACutoff: a per-item standardized-loadings plot (observed marker against each item's simulated distribution and expected range, in the style of RMitemInfitPlot), and a faceted plot of the CFI / RMSEA / SRMR distributions with the observed value overlaid.

Usage

RMdimCFAPlot(simfit, data, percentile = NULL)

Arguments

simfit	The list returned by RMdimCFACutoff.
data	The item-response data the CFA was run on (the same items used for the cutoff). Required: the observed values are computed from it.
percentile	Numeric in (50, 100) or NULL. When supplied, the cutoffs and flags are recomputed at this percentile from the stored simulated distributions (no re-simulation). When NULL (default), the percentile from the original RMdimCFACutoff() call is reused.

Value

A named list of two ggplot objects:

loadings

Per-item standardized loadings: simulated distribution (dots), expected interval, and the observed loading as a diamond (red when flagged).

fit

Faceted CFI / RMSEA / SRMR simulated distributions with the observed value and cutoff overlaid.

See Also

RMdimCFACutoff, RMdimCFA

Examples

if (requireNamespace("lavaan", quietly = TRUE) &&
    requireNamespace("ggplot2", quietly = TRUE) &&
    requireNamespace("eRm", quietly = TRUE)) {
  data("raschdat1", package = "eRm")
  sim <- RMdimCFACutoff(raschdat1[, 1:8], iterations = 50,
                        parallel = FALSE, seed = 1)
  plots <- RMdimCFAPlot(sim, data = raschdat1[, 1:8])
  plots$loadings
  plots$fit
}

---

Martin-Löf Test of Unidimensionality

Description

Likelihood-ratio test of unidimensionality against an a priori specified multidimensional alternative, generalised to polytomous Rasch / partial credit models (Christensen, Bjorner, Kreiner, & Petersen, 2002). The p-value is obtained by parametric-bootstrap (Monte Carlo) sampling under the unidimensional null, following Christensen & Kreiner (2007), because the asymptotic chi-square approximation is biased toward conservatism for realistic sample sizes – especially with polytomous items, where the degrees of freedom can be very large.

Usage

RMdimMartinLof(
  data,
  partition,
  iterations = 1000L,
  stopping = c("none", "sequential"),
  h = 50L,
  alpha = 0.05,
  parallel = TRUE,
  n_cores = NULL,
  verbose = FALSE,
  seed = NULL
)

Arguments

data	A data.frame or matrix of item responses (0-based, non-negative integers). Complete cases only: rows with any NA on the items named in partition are dropped. Missingness on items outside partition is ignored.
partition	The hypothesised partition of items into subscales. One of: a list of column-name or column-index vectors, e.g. list(c("I1","I2","I3"), c("I4","I5","I6")); a vector of length ncol(data) indicating each item's subscale (factor, character, or integer), e.g. c(1,1,1,2,2,2). Each subscale must contain at least two items. Subscales must not overlap; items not assigned to any subscale are dropped with a warning.
iterations	Integer. Maximum number of Monte Carlo iterations (default 1000).
stopping	Character. "none" (default) runs all iterations. "sequential" uses Besag & Clifford's (1991) sequential rule: stop as soon as h simulated statistics have exceeded the observed value. The sequential strategy substantially reduces compute time when H0 holds but cannot be parallelised.
h	Integer. Sequential-stopping count threshold (default 50). Ignored when stopping = "none".
alpha	Numeric in (0, 1). Nominal significance level used only for the rejected flag in the result; default 0.05.
parallel	Logical. Use parallel processing via mirai (default TRUE). Ignored when stopping = "sequential".
n_cores	Integer or NULL. Number of parallel workers. When NULL, getOption("mc.cores") is checked first; if neither is set, falls back to sequential with a warning.
verbose	Logical. Show a progress bar (default FALSE).
seed	Integer or NULL. Random seed for reproducibility. Items are processed internally in a fixed (alphabetical) order, so the same seed reproduces the same p-value regardless of how the data's columns are arranged and regardless of the parallel setting. See easyRasch2-reproducibility for what this guarantees and how it interacts with parallel.

Details

This is not a routine screening tool. The test requires an a priori partition of items into subscales; using it post-hoc on, e.g., the partition suggested by RMdimResidualPCA()'s PC1 sign would inflate the Type-I error rate. Both source papers state this explicitly.

Test statistic. With items partitioned into D subscales, total score $t$ and subscores $(t_1, \ldots, t_D)$ (Christensen et al. 2002, eq. 22):

$T = 2\Bigl[\sum_{t_1, \ldots, t_D} n_{t_1, \ldots, t_D}\log(n_{t_1, \ldots, t_D}/N) - \sum_t n_t\log(n_t/N) - \ell_C(\hat{\epsilon}) + \sum_d \ell_C(\hat{\epsilon}^{(d)})\Bigr]$

where $\ell_C$ is the conditional log-likelihood and the $\hat{\epsilon}^{(d)}$ are CML estimates on the d-th subscale alone. CML fits use psychotools::raschmodel() (RM) or psychotools::pcmodel() (PCM) for speed.

Monte Carlo sampling under H0. Following Christensen & Kreiner (2007): (a) sample N total scores from the empirical score distribution $n_t/N$; (b) for each sampled score, sample an item-response vector from the conditional distribution $p(x \mid t, \hat{\epsilon})$ given by eq. 4 of the paper. Step (b) uses the recursive $\gamma$-function algorithm (Andersen, 1995, eq. 15.22 and 15.27) for dichotomous and polytomous items alike: items are filled in from the last backwards, each conditional on the score still to be distributed over the items before it. The nested $\gamma^{(1)}, \ldots, \gamma^{(n)}$ depend only on the item parameters, which are fixed across the run, so the recursion is evaluated once rather than per person.

The faster dichotomous shortcut described by Christensen & Kreiner (2007, p. 23), choosing $s$ items one at a time with probabilities $\phi_i / \sum_{\mathrm{remaining}} \phi$, is deliberately not used: successive sampling of that kind follows a Wallenius-type distribution, whereas the Rasch conditional distribution weights a set of items by the product of their easinesses. See the note on validation below.

Iterations that fail (e.g., simulated dataset has an empty category for an item) are silently dropped.

Item parameters are estimated once on the observed data and held fixed across MC iterations. Christensen & Kreiner (2007) use the extended likelihood function (Tjur, 1982) with the empirical score distribution as a non-parametric estimate of the latent distribution, so no distributional assumption about $\theta$ is needed.

Missing data. Complete cases only. The statistic is built from the joint table of subscores, so a respondent needs a defined subscore on every subscale; there is no partial-missingness path, and the SAS macro this implementation was validated against has the same requirement. Rows with NA are dropped after items outside partition have been removed, so missingness on items the test does not use costs no respondents. sample_n reports the respondents analysed and sample_n_total the raw input rows.

Validation. The statistic, the conditional log-likelihoods, the conditional sampler and the expected counts and residuals of RMdimMartinLofResiduals have been checked against the pml SAS macro (Christensen, 2004), kindly shared by Karl Bang Christensen. The statistic and the residual table agree with the macro to numerical precision. The sampler follows the macro's recursion rather than the dichotomous shortcut published in Christensen & Kreiner (2007).

Value

A list with components:

T_obs

Observed Martin-Löf likelihood-ratio statistic.

p_value

Monte Carlo p-value with (n_exceed + 1) / (n + 1) correction. The attainable p-values are k / (n + 1) for ⁠k = 1, ..., n + 1⁠, so the p-value's resolution is limited by the number of iterations: with 100 iterations the smallest attainable value is 1/101 = 0.0099. A reported p-value equal to the floor (see p_value_floor) means no simulated statistic reached the observed one and should be read as "p < floor" – the true p-value may be much smaller; increase iterations for finer resolution.

p_value_floor

The smallest attainable p-value, 1 / (actual_iterations + 1).

actual_iterations

Number of successful MC iterations completed.

rejected

Logical: is p_value < alpha?

partition

Normalised partition (list of integer indices).

n_subscales

Number of subscales.

is_polytomous

Whether a PCM was fitted.

sample_n

Number of complete cases analysed, counted over the items in partition only.

sample_n_total

Number of respondents in the raw input data, before the complete-case filter.

sample_has_na

Logical. Whether the partitioned items contained any missing values.

n_items

Number of items.

stopping

The stopping strategy used.

h

The sequential-stopping count, or NA for stopping = "none".

T_rep

Numeric vector of successful MC test statistics.

wle_scores

data.frame with one row per person and one column per subscale (subscale_1_wle, ..., subscale_D_wle), giving Warm's Weighted Likelihood Estimate of theta from a CML fit on each subscale alone. Persons whose subscore equals the minimum or maximum on a subscale produce non-finite WLEs (Inf / -Inf) and are excluded from wle_correlation pairwise.

wle_correlation

data.frame of pairwise Pearson correlations between subscale WLEs, with columns subscale_a, subscale_b, r, ci_lower, ci_upper (95% CI from stats::cor.test), p_value, and n (number of persons with finite WLEs on both subscales). One row per pair; for D = 2, a single row. Useful as an effect-size companion to p_value – a rejected test with r near 1 indicates a small effect; r clearly below 1 indicates substantive multidimensionality.

References

Andersen, E. B. (1995). Polytomous Rasch models and their estimation. In G. H. Fischer & I. W. Molenaar (Eds.), Rasch models: Foundations, recent developments, and applications (pp. 271-291). Springer-Verlag.

Besag, J., & Clifford, P. (1991). Sequential Monte Carlo p-values. Biometrika, 78(2), 301-304. doi:10.1093/biomet/78.2.301

Christensen, K. B. (2004). pml: A SAS macro for testing unidimensionality in polytomous Rasch models (Technical note). National Institute of Occupational Health, Denmark, and Department of Biostatistics, University of Copenhagen.

Christensen, K. B., Bjorner, J. B., Kreiner, S., & Petersen, J. H. (2002). Testing unidimensionality in polytomous Rasch models. Psychometrika, 67(4), 563-574. doi:10.1007/BF02295132

Christensen, K. B., & Kreiner, S. (2007). A Monte Carlo approach to unidimensionality testing in polytomous Rasch models. Applied Psychological Measurement, 31(1), 20-30. doi:10.1177/0146621605286204

See Also

RMdimResidualPCA, RMdimResidualPCACutoff

Examples

set.seed(1)
# Build 2-dimensional polytomous data: 4 items per subscale, 5 categories
n     <- 400
theta1 <- rnorm(n)
theta2 <- 0.6 * theta1 + sqrt(1 - 0.6^2) * rnorm(n)
make_pcm <- function(theta, n_items, taus) {
  sapply(seq_len(n_items), function(j) {
    # ... toy simulation here
    sample(0:4, n, replace = TRUE)
  })
}
dat <- cbind(make_pcm(theta1, 4, NULL), make_pcm(theta2, 4, NULL))
colnames(dat) <- paste0("I", 1:8)

# Few iterations for a fast example; use 1000+ in real analyses
RMdimMartinLof(dat,
            partition = list(c("I1","I2","I3","I4"),
                             c("I5","I6","I7","I8")),
            iterations = 100, parallel = FALSE, seed = 1)

# Sequential stopping: stop as soon as h = 25 simulated statistics exceed
# the observed one (cuts compute time under H0).
RMdimMartinLof(dat,
            partition = c(1,1,1,1,2,2,2,2),
            iterations = 200, stopping = "sequential", h = 25,
            seed = 1)

---

Standardised Residuals from the Joint Subscore Distribution

Description

Diagnostic accompanying RMdimMartinLof: per-cell standardised residuals from the joint distribution of subscores under unidimensionality (Christensen, Bjorner, Kreiner, & Petersen, 2002, eq. 13). Useful for identifying where a partition deviates from the unidimensional null rather than just whether it does (which RMdimMartinLof() answers).

Usage

RMdimMartinLofResiduals(
  data,
  partition,
  output = c("kable", "dataframe", "ggplot"),
  flag_threshold = 2,
  color_by = c("residual", "n"),
  color_limits = NULL,
  min_expected = NULL
)

Arguments

data	A data.frame or matrix of item responses (0-based, non-negative integers). Complete cases only: rows with any NA on the items named in partition are dropped. Missingness on items outside partition is ignored.
partition	Same format as in RMdimMartinLof: a list of item-name/index vectors, or a length-ncol(data) vector of group labels. Each subscale must contain at least 2 items.
output	Character. "kable" (default) for a 2-D pipe-format residual table when D = 2 (long-format kable when D > 2), "dataframe" for the underlying long-format data.frame, or "ggplot" for a diverging-fill heatmap (with facet_wrap over t3 for D = 3, error for D > 3).
flag_threshold	Numeric. Cells with ⁠|residual| > flag_threshold⁠ are flagged: marked as ⁠**bold**⁠ in the kable, shown in the flagged column of the dataframe. Default 2.
color_by	Character. For output = "ggplot", what the tile fill colour encodes. "residual" (default) – diverging red-white-blue scale centred at 0, the most directly diagnostic. "n" – sequential blue scale on the observed cell count, useful for spotting whether large-magnitude residuals are driven by sparse cells. Either way, the numeric residual is printed inside each cell.
color_limits	Numeric length-2 vector or NULL. Caps the colour scale for output = "ggplot". Default c(-5, 5) when color_by = "residual" (sparse-cell residuals from (o-e)/sqrt(n*p*(1-p)) can be enormous when expected counts are tiny; capping the scale stops outliers from compressing the rest of the plot). The unclipped residual values still appear as cell labels. NULL for color_by = "n", which uses the natural data range.
min_expected	Numeric or NULL. If set, cells with expected < min_expected have their residual set to NA (they appear grey in the heatmap and are not flagged). Default NULL (no filtering). Setting min_expected = 1 (or 5) is the analogue of the Cochran rule for sparse-cell chi-square contributions and removes residuals whose asymptotic standard normal approximation is unreliable.

Details

For each cell of the joint subscore table (indexed by $(t_1, \ldots, t_D)$), the conditional probability under H0 given the total score $t = \sum_d t_d$ is

$p(t_1, \ldots, t_D \mid t) = \prod_d \gamma^{(d)}_{t_d} / \gamma_t,$

the expected count is $e = n_t \cdot p$, and the residual is $(o - e) / \sqrt{n_t \cdot p \cdot (1 - p)}$. CML estimates from the unidimensional model are used for the $\gamma$-functions.

Reading the table (D = 2): under the unidimensional null, residuals should be patternless and roughly N(0, 1). Multidimensionality with positively correlated dimensions typically shows up as positive residuals at the corners of each antidiagonal (high t1 + low t2, low t1 + high t2) and negative residuals near the antidiagonal centre (matched subscores). Negatively correlated dimensions show positive residuals at the table corners (high/low and low/high) and negative residuals at high/high and low/low. See Christensen et al. (2002, section7) for a worked example.

Cells where the total score has no observed cases (n_t = 0) are uninformative and are dropped from the output.

Complete cases only, as in RMdimMartinLof, and rows with NA are dropped only after items outside partition have been removed. The reported ⁠n = X of Y respondents⁠ counts respondents complete on the partitioned items against the raw input rows.

Expected counts and residuals agree to numerical precision with the pml SAS macro (Christensen, 2004), kindly shared by Karl Bang Christensen. The one deliberate difference is at the extreme totals t = 0 and t = max, where the conditional probability is 1 and the residual is undefined: the macro prints 0, this function returns NA.

Value

• output = "kable": a knitr_kable object. For D = 2, a wide table with rows = t1, columns = t2, cells = standardised residual (⁠**bold**⁠ if flagged, em-dash if NA). For D > 2, long-format with one row per cell.

• output = "dataframe": a long-format data.frame with columns t1, ..., tD, total, observed, expected, residual, flagged.

• output = "ggplot": a geom_tile() heatmap (D = 2 or 3 only).

References

Christensen, K. B. (2004). pml: A SAS macro for testing unidimensionality in polytomous Rasch models (Technical note). National Institute of Occupational Health, Denmark, and Department of Biostatistics, University of Copenhagen.

Christensen, K. B., Bjorner, J. B., Kreiner, S., & Petersen, J. H. (2002). Testing unidimensionality in polytomous Rasch models. Psychometrika, 67(4), 563-574. doi:10.1007/BF02295132

See Also

RMdimMartinLof

Examples

set.seed(1)
dat <- as.data.frame(matrix(sample(0:1, 400 * 8, replace = TRUE),
                            nrow = 400, ncol = 8))
colnames(dat) <- paste0("I", 1:8)

# Wide kable table for D = 2
RMdimMartinLofResiduals(dat,
                     partition = list(c("I1","I2","I3","I4"),
                                      c("I5","I6","I7","I8")))

# Heatmap
if (requireNamespace("ggplot2", quietly = TRUE)) {
  RMdimMartinLofResiduals(dat,
                       partition = c(1,1,1,1,2,2,2,2),
                       output = "ggplot")
}

# Underlying data.frame for custom analysis
df <- RMdimMartinLofResiduals(dat,
                           partition = c(1,1,1,1,2,2,2,2),
                           output = "dataframe")
df[df$flagged, ]

---

PCA of Standardized Rasch Residuals

Description

Fits a Rasch model by CML via psychotools (a dichotomous item is a 2-category PCM), extracts the standardized residuals $(x - E)/\sqrt{Var}$ at WLE person locations, and runs an unrotated principal-component analysis on those residuals via stats::prcomp(). The function reports the top n_components eigenvalues and their proportions of unexplained variance, and optionally compares the first-contrast eigenvalue against a simulation-based bound from RMdimResidualPCACutoff.

Usage

RMdimResidualPCA(
  data,
  cutoff = NULL,
  p_value = FALSE,
  n_components = 5L,
  output = "kable"
)

Arguments

data	A data.frame or matrix of item responses. Items must be scored starting at 0 (non-negative integers). Rows with any NA are dropped before PCA, since prcomp() does not accept missing values.
cutoff	Optional. The list returned by RMdimResidualPCACutoff (its suggested_cutoff is used), or a single numeric value to use as the cutoff directly. When provided, the result includes a Flagged column (logical: is the eigenvalue above the simulated bound?) and the kable caption notes the cutoff.
p_value	Logical. When TRUE, adds a one-sided bootstrap p-value for the first-contrast eigenvalue: the proportion of simulated first-contrast eigenvalues at least as large as the observed one, ⁠(1 + #\{lambda* >= lambda\}) / (B + 1)⁠. Requires the full RMdimResidualPCACutoff object as cutoff (it carries the simulated eigenvalues in ⁠$results⁠); a bare numeric cutoff is not sufficient. The simulated null is the distribution of the largest eigenvalue, so the p-value applies to PC1 only (NA for the other components). This is a single test — no multiplicity correction is involved. Default FALSE.
n_components	Integer. Number of eigenvalues to report. Capped at the number of items. Default 5.
output	Character. "kable" (default) for a formatted knitr::kable() table, "dataframe" for the underlying data.frame, or "ggplot" for a ggplot of PC1 loadings against item locations. "loadings" is accepted as a backward-compatible alias for "ggplot".

Details

Rule-of-thumb thresholds for the first-contrast eigenvalue (e.g., the "> 2" heuristic occasionally cited from Winsteps documentation) are not reliable indicators of multidimensionality; the first-contrast eigenvalue under a correctly fitting unidimensional model varies systematically with sample size, test length, and item-parameter spread. Empirical (simulated) bounds tailored to the data structure should be used instead — see RMdimResidualPCACutoff, and Chou & Wang (2010) for the underlying simulation argument.

The PCA is performed on the standardized residuals $(x - E)/\sqrt{Var}$ from the shared CML/WLE engine (CML item parameters via psychotools, WLE person locations). The reported eigenvalues are unrotated; rotation is appropriate for interpreting a multidimensional solution but obscures the dominant first contrast that dimensionality assessment is concerned with.

Item locations on the loadings plot are the per-item mean of the CML Andrich thresholds.

The variance partition follows Linacre's convention: per-item observed variance is compared to per-item expected variance under the fitted model, summed across items. Expected scores are computed from the CML item parameters and WLE person locations. WLE is finite at extreme scores, so all persons are retained (the previous MLE partition dropped extreme-score cases).

Bootstrap p-value. When p_value = TRUE, the observed first-contrast eigenvalue is compared against the simulated null distribution of largest eigenvalues (from cutoff$results), giving the one-sided Monte-Carlo p-value ⁠(1 + #\{lambda* >= lambda\}) / (B + 1)⁠. Because the maximum eigenvalue is a single family-wise statistic, no multiplicity correction applies. The p-value is model-conditional and sample-size-sensitive; it is reported alongside the simulated cutoff, not in place of it, and can be no smaller than 1 / (B + 1).

Value

• If output = "kable": a knitr_kable object with columns Component, Eigenvalue, Proportion of variance (and Flagged when cutoff is provided; p when p_value = TRUE). The caption gives the variance partition (% of total observed variance explained by measures vs. unexplained), the model fitted, sample size, and cutoff/p-value metadata if applicable.

• If output = "dataframe": a data.frame with columns Component, Eigenvalue, Proportion_of_variance (and Flagged when cutoff is provided; p when p_value = TRUE, non-NA for PC1 only). The variance partition is attached as the "variance_partition" attribute — a list with elements total, explained, unexplained, pct_explained, pct_unexplained, n_persons. Access via attr(result, "variance_partition").

• If output = "ggplot": a ggplot showing each item's PC1 loading on the x-axis and Rasch item location on the y-axis, with dashed reference lines at zero, and the variance partition in the figure caption. Item names are labelled via ggrepel::geom_text_repel() when ggrepel is installed; otherwise plain geom_text().

References

Chou, Y.-T., & Wang, W.-C. (2010). Checking dimensionality in item response models with principal component analysis on standardized residuals. Educational and Psychological Measurement, 70(5), 717-731. doi:10.1177/0013164410379322

See Also

RMdimResidualPCACutoff

Examples

set.seed(1)
dat <- as.data.frame(
  matrix(sample(0:1, 200 * 12, replace = TRUE), nrow = 200, ncol = 12)
)
colnames(dat) <- paste0("I", 1:12)

# Default kable output
RMdimResidualPCA(dat)

# PC1 loadings vs item location plot
if (requireNamespace("ggplot2", quietly = TRUE) &&
    requireNamespace("ggrepel", quietly = TRUE)) {
  RMdimResidualPCA(dat, output = "ggplot")
}

# Simulation-based cutoff (use 250+ iterations in real analyses)
bound <- RMdimResidualPCACutoff(dat, iterations = 50, parallel = FALSE, seed = 1)
RMdimResidualPCA(dat, cutoff = bound)

# With the one-sided bootstrap p-value for the first contrast
RMdimResidualPCA(dat, cutoff = bound, p_value = TRUE)

---

Simulation-based Cutoff for First-Contrast Eigenvalue

Description

Parametric bootstrap producing an empirical upper percentile for the largest eigenvalue from a PCA of standardized residuals under a correctly fitting Rasch model. For each iteration the function resamples theta values from the data-based estimates, simulates response data, refits the appropriate Rasch model, and extracts the first-contrast eigenvalue. Several upper-tail percentiles of the resulting distribution are returned; the 99th percentile is reported as the suggested cutoff (matching the convention used by RMlocdepQ3Cutoff).

Usage

RMdimResidualPCACutoff(
  data,
  iterations = 250,
  parallel = TRUE,
  n_cores = NULL,
  verbose = FALSE,
  seed = NULL
)

Arguments

data	A data.frame or matrix of item responses (0-based, non-negative integers).
iterations	Integer. Number of simulation iterations. Default 250.
parallel	Logical. Use parallel processing via mirai. Default TRUE.
n_cores	Integer or NULL. Number of parallel workers. When NULL, getOption("mc.cores") is checked first; if neither is set, parallel = TRUE falls back to sequential with a warning.
verbose	Logical. Show a progress bar. Default FALSE.
seed	Integer or NULL. Random seed for reproducibility. See easyRasch2-reproducibility for what this guarantees and how it interacts with parallel.

Details

Rule-of-thumb cutoffs for the first-contrast eigenvalue depend strongly on sample size, test length, and item-parameter spread; simulation-based cutoffs tailored to the data are more defensible (Chou & Wang, 2010).

The generating model uses CML item parameters (psychotools) and WLE person locations. Per iteration: theta values are sampled with replacement from the WLE estimates; response data are simulated under the model (psychotools::rrm for dichotomous, partial-credit simulator for polytomous); CML/WLE standardized residuals are computed (the same engine as the observed analysis); prcomp() is run; the first-contrast eigenvalue is recorded.

Iterations that fail (e.g., due to a degenerate simulated dataset where some category isn't represented) are silently dropped. The iarm package is not required.

Value

A list with components:

results

data.frame: iteration, eigenvalue.

p95, p99, p995, p999

Empirical percentiles of eigenvalue.

max

The largest simulated eigenvalue.

suggested_cutoff

The 99th percentile (p99) — pass this list back into RMdimResidualPCA via ⁠cutoff = ⁠, or use suggested_cutoff directly.

suggested_cutoff_percentile

The percentile used for suggested_cutoff, currently always 99.

actual_iterations

Number of successful iterations.

sample_n

Number of complete cases used.

item_names

Character vector of item names.

References

Chou, Y.-T., & Wang, W.-C. (2010). Checking dimensionality in item response models with principal component analysis on standardized residuals. Educational and Psychological Measurement, 70(5), 717-731. doi:10.1177/0013164410379322

See Also

RMdimResidualPCA

Examples

set.seed(1)
dat <- as.data.frame(
  matrix(sample(0:1, 200 * 12, replace = TRUE), nrow = 200, ncol = 12)
)
colnames(dat) <- paste0("I", 1:12)

# Few iterations for a fast example; use 250+ in real analyses
bound <- RMdimResidualPCACutoff(dat, iterations = 50, parallel = FALSE, seed = 1)
bound$suggested_cutoff

RMdimResidualPCA(dat, cutoff = bound)

---

Item Category Probability Curves

Description

Plots model-implied response-category probability curves for each item as a function of the latent trait $\theta$. Item parameters are estimated by conditional maximum likelihood via psychotools::pcmodel() (a dichotomous item is a 2-category PCM). Each item gets its own facet panel, with one curve per response category coloured from low to high using the viridis palette. Comparable in scope to eRm::plotICC() and mirt's trace plots, with a ggplot2 / viridis output and optional descriptive labels for items and categories.

Usage

RMitemCatProb(
  data,
  item_labels = NULL,
  category_labels = NULL,
  theta_range = NULL,
  n_points = 200L,
  viridis_option = "D",
  viridis_end = 0.95,
  facet_ncol = NULL,
  label_wrap = 25L,
  line_width = 0.9,
  font = "sans",
  output = c("ggplot", "dataframe"),
  label_curves = c("legend", "path"),
  item = NULL,
  text_size = 4
)

Arguments

data	A data.frame or matrix of item responses. Items must be scored starting at 0 (non-negative integers). Missing values (NA) are allowed — the model fit handles them.
item_labels	Optional character vector of descriptive item labels (facet strip titles). Must be the same length as ncol(data). If NULL (the default), column names are used.
category_labels	Optional character vector of labels for the response categories (legend). Must be the same length as the number of categories spanning from 0 to the maximum observed value. If NULL (the default), numeric category values are used.
theta_range	Numeric length 2. Range of the latent trait $\theta$ (logits) plotted along the x-axis. When NULL (default) the range is c(-4, 4) for label_curves = "legend" and a wider c(-5, 5) for label_curves = "path" — the wider range gives geomtextpath more horizontal room to place per-curve labels at their peaks. Pass an explicit c(lo, hi) to override.
n_points	Integer. Number of evenly-spaced $\theta$ values at which the curves are evaluated. Default 200L.
viridis_option	Character. Viridis palette identifier. One of "A" through "H". Default "D" (viridis green).
viridis_end	Numeric in (0, 1]. Upper end of the viridis palette range; lower values keep the palette inside its mid-tones (avoids the very bright yellow at 1.0). Default 0.95.
facet_ncol	Optional integer. Number of columns in the facet layout. Default NULL (ggplot2::facet_wrap() auto-layout).
label_wrap	Integer. Characters per line for facet-strip label wrapping. Default 25L.
line_width	Numeric. Line width for the probability curves. Default 0.9.
font	Character. Font family for all text. Default "sans".
output	Character. Either "ggplot" (default) for a ggplot object, or "dataframe" for the underlying long-format probability table.
label_curves	Character. How response categories are identified. "legend" (default) draws lines and a separate colour legend with one swatch per category — the standard multi-item presentation. "path" uses geomtextpath::geom_textpath() to write each category's label along its own curve (a la classic IRT trace plots) and suppresses the legend; valid for a single item at a time because narrow facets do not give path labels enough horizontal room to render legibly. The model is still fit on the full multi-item data (PCM thresholds require multi-item input for CML estimation); the item argument then selects which item's curves to plot. Requires the optional geomtextpath package.
item	Character or integer. Used only when label_curves = "path". Either an item name (matching a column in data) or a column index, identifying which item's category curves to plot. Required for path mode.
text_size	Numeric. Used only when label_curves = "path". Size of the path labels in mm. Default 4.

Details

For each polytomous item i with response categories $0, 1, \ldots, K_i$ and threshold parameters $\delta_{i,1}, \ldots, \delta_{i,K_i}$ (CML, psychotools), the PCM category probability is

$P(X_i = k \mid \theta) = \frac{\exp(\sum_{j=1}^{k} (\theta - \delta_{i,j}))}{\sum_{k'=0}^{K_i} \exp(\sum_{j=1}^{k'} (\theta - \delta_{i,j}))},$

with the empty sum (when $k = 0$) taken as zero. For dichotomous items the function fits a Rasch model and treats the item difficulty $\delta_i = -\beta_i$ as the single threshold, recovering the standard two-category logistic ICC.

The colour mapping uses scale_color_viridis_c() against the integer category value, so the natural ordering of response categories is preserved visually — low categories at one end of the palette, high categories at the other. When category_labels is provided, the legend uses those labels (e.g., "Never" / "Sometimes" / "Often") while the colour mapping stays on the integer category value.

Items with fewer response categories than the maximum (e.g., an otherwise four-category scale with one three-category item) contribute only the categories they actually have to their own facet — the y-axis still spans $[0, 1]$.

Value

• If output = "ggplot": a ggplot2::ggplot object, one facet per item.

• If output = "dataframe": a long-format data.frame with columns Item (factor in column order), Category (integer), and Theta, Probability (numeric), one row per item × category × theta gridpoint.

See Also

RMitemICCPlot() for conditional ICCs binned by total score, RMitemHierarchy() for item threshold locations on the logit scale.

Examples

if (requireNamespace("eRm", quietly = TRUE) &&
    requireNamespace("ggplot2", quietly = TRUE)) {
  data(pcmdat2, package = "eRm")

  # Default plot
  RMitemCatProb(pcmdat2)

  # Custom item and category labels
  RMitemCatProb(
    pcmdat2,
    item_labels     = c("Mood", "Sleep", "Appetite", "Energy"),
    category_labels = c("Never", "Sometimes", "Often")
  )

  # Underlying probability data
  df <- RMitemCatProb(pcmdat2, output = "dataframe")
  head(df)

  # Single-item plot with labels written along each curve
  # (classic IRT trace-plot style). Model is still fit on all four
  # items; `item` picks which one to plot.
  if (requireNamespace("geomtextpath", quietly = TRUE)) {
    RMitemCatProb(
      pcmdat2,
      category_labels = c("Never", "Sometimes", "Often"),
      label_curves    = "path",
      item            = "I1"
    )
  }
}

---

Item-Threshold Hierarchy Plot

Description

Visualises item and threshold locations on the logit scale for a Partial Credit Model. Items are sorted by their (mean-threshold) location; each item shows its location as a black diamond and its individual thresholds as coloured dots with confidence-interval error bars.

Usage

RMitemHierarchy(
  data,
  show_numbers = TRUE,
  sem_multiplier = 1.405,
  item_labels = NULL,
  output = c("ggplot", "dataframe")
)

Arguments

data	A data.frame or matrix of polytomous item responses (non-negative integers, 0-based, max value > 1). One column per item, one row per person.
show_numbers	Logical. When TRUE (default), prints the numerical item location and threshold values next to the points on the plot. When FALSE, only the threshold labels (T1, T2, ...) are shown.
sem_multiplier	Numeric multiplier for the threshold SE used to draw the error bars. Default 1.405 (84% CI). Common alternatives: 1.96 (95% CI), 2.576 (99% CI).
item_labels	Optional character vector of length ncol(data) providing descriptive labels for items. Default NULL uses the column names of data. Labels are appended to the column names on the y-axis (name - label) and wrapped at 36 characters via base-R strwrap().
output	One of "ggplot" (default) – a faceted hierarchy figure – or "dataframe" – the long-format underlying data with one row per (item × threshold).

Details

Threshold locations are centred at the grand mean of all thresholds across items, so the dashed reference line at 0 represents the mean threshold location across the scale.

Confidence intervals around thresholds are 84% by default (sem_multiplier = 1.405), following Payton, Greenstone, & Schenker (2003) – non-overlap of 84% intervals approximately corresponds to a two-sample significance test at $\alpha = 0.05$. Use sem_multiplier = 1.96 for 95% intervals.

Polytomous only. Dichotomous items have a single threshold that coincides with the item location; the hierarchy plot is visually degenerate in that case. For dichotomous data use RMtargeting() or RMscoreSE() instead.

Centring convention. The CML PCM thresholds (estimated via psychotools::pcmodel()) are shifted so that their grand mean is zero; each item's location is then the mean of its centred thresholds. The dashed horizontal reference line on the plot marks this zero – i.e., the average threshold across all items.

Value

Either a ggplot (default) or a data.frame with columns Item, ItemLabel, Threshold, ThresholdLocation, ThresholdSE, and ItemLocation (the per-item mean of the centred thresholds).

References

Payton, M. E., Greenstone, M. H., & Schenker, N. (2003). Overlapping confidence intervals or standard error intervals: What do they mean in terms of statistical significance? Journal of Insect Science, 3(34), 1-6. doi:10.1093/jis/3.1.34

See Also

RMtargeting(), RMscoreSE(), RMitemICCPlot()

Examples

if (requireNamespace("eRm", quietly = TRUE)) {
  data("pcmdat2", package = "eRm")
  RMitemHierarchy(pcmdat2)

  # 95% CI instead of 84%
  RMitemHierarchy(pcmdat2, sem_multiplier = 1.96)

  # Underlying data.frame
  RMitemHierarchy(pcmdat2, output = "dataframe")
}

---

Conditional Item Characteristic Curves

Description

Faceted panel of Conditional Item Characteristic Curves (cICCs), one per item. Each panel shows the model-expected conditional item score $E[X_i \mid R = r]$ against the total score $r$, together with the average observed item score within class intervals of the total score and their confidence intervals. When a dif_var is supplied, observed averages are computed separately per group, turning each panel into a visual DIF check, and the partial-gamma DIF magnitude is reported per item.

Usage

RMitemICCPlot(
  data,
  dif_var = NULL,
  method = c("quantile", "width", "score", "manual", "cut"),
  class_intervals = 4,
  score_breaks = NULL,
  ci = TRUE,
  error_band = FALSE,
  conf_level = 0.95,
  min_n = 8L,
  items = NULL,
  output = c("patchwork", "list")
)

Arguments

data	A data.frame or matrix of item responses (non-negative integers, 0-based). One column per item, one row per person.
dif_var	Optional vector of length nrow(data) (factor or coercible to factor) defining a DIF grouping variable. When supplied, observed averages are shown separately per group and the partial-gamma DIF coefficient (iarm::partgam_DIF(), as in RMdifGamma()) is annotated per panel. Default NULL (no DIF).
method	How the observed means are grouped along the total score (see Class intervals in Details). One of "quantile" (default): class_intervals groups of approximately equal numbers of respondents; "width": class_intervals equal-width intervals over the observed total-score range; "score": every observed total score is its own group; or "manual": groups defined by score_breaks. "cut" is accepted as a legacy alias for "quantile".
class_intervals	Integer >= 2. Number of class intervals for method = "quantile" and method = "width". Default 4. Ignored (without an error) by method = "score" and method = "manual".
score_breaks	Integer vector for method = "manual": the total scores at which a new group starts, in increasing order. For example, score_breaks = c(3, 6, 8) on a 0-9 scale defines the groups 0-2, 3-5, 6-7, and 8-9 (the endpoints are added automatically). Serves the same purpose as lower.groups in RASCHplot's CICCplot(), without the leading zero. Default NULL.
ci	Logical. Draw confidence intervals (error bars) on the observed class-interval means. Default TRUE.
error_band	Logical. Add a shaded band around the model-expected curve: the model's interval for the observed mean at each total score, $E \pm z \sqrt{\mathrm{Var}/n_r}$ (n_r = number of respondents at that total score), as in RASCHplot's CICCplot. Complementary to ci, not a replacement – see Details. Default FALSE.
conf_level	Numeric in (0, 1). Confidence level for the observed error bars and the model band. Default 0.95.
min_n	Integer. A group-by-interval cell needs at least this many respondents to contribute an observed point + CI; sparser cells are skipped. Default 8 (the package's per-cell stability floor).
items	Optional character or integer vector selecting which items to plot. The model is always fitted on all items; only the rendering is filtered. Default NULL (all items).
output	One of "patchwork" (default) – composite patchwork figure – or "list" – a named list of per-item ggplot objects.

Details

The model curve and its conditional variance come from the shared CML engine (CML item thresholds via psychotools and the exact conditional distribution of each item score given the total score, via elementary symmetric functions). The conditional-ICC approach follows Buchardt, Christensen & Jensen (2023) and their RASCHplot package; the implementation here is native to easyRasch2's CML/WLE engine.

Conditioning. The expected curve is the exact conditional expectation of the item score given the total score (it accounts for the item being part of the total), not a marginal ICC.

Class intervals. Following Buchardt et al. (2023), the empirical item means can be shown for each value of the total score or for each value of the grouped total score; all grouping happens on the total-score scale (never on an estimated latent score). Four grouping rules are available. method = "quantile" (the default) splits the observed total scores into class_intervals bins of approximately equal numbers of respondents, using common boundaries so DIF groups share the x-axis – the same approach as iarm::ICCplot()'s class intervals (Hmisc::cut2(g = ...)). method = "width" instead splits the observed score range into class_intervals intervals of equal width, which reads naturally on the score axis but can leave sparse intervals (pruned by min_n). method = "score" uses every observed total score as its own group – the maximal-resolution display, sensible for short scales, with min_n pruning thinly-populated scores. method = "manual" places the group boundaries exactly where you say via score_breaks (RASCHplot's lower.groups concept). Each bin contributes one observed point at its mean total score. If the score distribution is too sparse to form class_intervals distinct bins, the quantile and width methods fall back to score-level points.

The figure caption names the grouping that was actually used, which is not always the one requested. Quantile grouping forms fewer bins than asked for when total scores tie at the boundaries, equal-width and manual grouping can define intervals no respondent falls into, and both binned methods fall back to score level on a degenerate score distribution. The caption reports the realised number of groups, how many of them contain respondents when that is fewer, and any fall back.

Confidence intervals. Observed error bars use the normal approximation $\bar{x}_l \pm z \sqrt{\mathrm{var}(x_l) / n_l}$ within each (group, interval) cell, clamped to the item's score range; cells with fewer than min_n respondents are dropped. In DIF mode this makes sparse group differences visibly uncertain rather than over-interpreted.

The model band (error_band). The shaded band is drawn around the model-expected curve at $E_{ri} \pm z \sqrt{V_{ri}/n_r}$, where $E_{ri}$ and $V_{ri}$ are the exact conditional mean and variance of the item score given total score $r$ (model quantities), and $n_r$ is the observed number of respondents at that total score. If the Rasch model is true, this is where the empirical mean item score at total score $r$ should fall, given how many people actually sit at that score; an observed diamond outside the band is localized graphical misfit. The band is deliberately the standard error of the mean ($\sqrt{V/n_r}$), not the far wider individual-response SD band, which would flag nothing. Its width therefore tracks the data: narrow where many respondents sit (a strict test), wide in sparse score regions (a lenient one), with gaps where $n_r = 0$ and collapsing to the curve at the deterministic extremes (total score 0 and the maximum). It is complementary to, not a replacement for, the ci error bars: the band is the model-implied uncertainty of the observed mean per raw score; the error bars are the sample-based uncertainty of the observed group means.

DIF magnitude. The annotated partial gamma (Bjorner et al., 1998) is the association between the item score and the group conditional on the rest score – a non-parametric effect size, complementary to the total-score-conditional visual. It is the same statistic as RMdifGamma().

Value

Either a patchwork/ggplot composite (default) or a named list of per-item ggplot objects (output = "list"). In DIF mode the per-item partial-gamma table is attached as attr(., "dif_gamma").

References

Buchardt, A.-S., Christensen, K. B., & Jensen, S. N. (2023). Visualizing Rasch item fit using conditional item characteristic curves in R. Psychological Test and Assessment Modeling, 65(2), 206-219.

Andersen, E. B. (1995). Polytomous Rasch models and their estimation. In G. H. Fischer & I. W. Molenaar (Eds.), Rasch Models: Foundations, Recent Developments, and Applications (pp. 271-291). Springer. (Conditional expectation of item scores given the total score; formula 15.22.)

Bjorner, J. B., Kreiner, S., Ware, J. E., Damsgaard, M. T., & Bech, P. (1998). Differential item functioning in the Danish translation of the SF-36. Journal of Clinical Epidemiology, 51(11), 1189-1202.

See Also

RMdifGamma, RMitemRestscore

Examples

if (requireNamespace("ggplot2", quietly = TRUE) &&
    requireNamespace("patchwork", quietly = TRUE) &&
    requireNamespace("iarm", quietly = TRUE)) {
  set.seed(42)
  sim_data <- as.data.frame(
    matrix(sample(0:2, 200 * 4, replace = TRUE), nrow = 200, ncol = 4)
  )
  colnames(sim_data) <- paste0("Item", 1:4)

  # Default: quantile groups (~equal numbers of respondents per group)
  RMitemICCPlot(sim_data)

  # Equal-width intervals on the total-score scale, with the model band
  RMitemICCPlot(sim_data, method = "width", error_band = TRUE)

  # Every total score as its own group (short scales)
  RMitemICCPlot(sim_data, method = "score")

  # Manual grouping: new groups start at total scores 3, 5, and 7
  RMitemICCPlot(sim_data, method = "manual", score_breaks = c(3, 5, 7))
}

---

Conditional Item Infit MSQ

Description

Computes conditional infit mean-square (MSQ) statistics for each item using iarm::out_infit(), enriched with item locations relative to the sample mean person location.

Usage

RMitemInfit(
  data,
  cutoff = NULL,
  p_value = NULL,
  correction = c("fwer", "fdr_bh", "fdr_by", "none"),
  alpha = 0.05,
  output = "kable",
  sort,
  statistic = "infit"
)

Arguments

data	A data.frame or matrix of item responses. Items must be scored starting at 0 (non-negative integers). Missing values (NA) are allowed, but at least one complete case (row with no NA) must be present.
cutoff	Optional. Default NULL (no cutoff applied, behaviour is identical to the current version). Can be: The return value of RMitemInfitCutoff (a list with ⁠$item_cutoffs⁠): the data.frame is extracted automatically and the number of simulation iterations, cutoff_method, and hdci_width are included in the kable caption. The ⁠$item_cutoffs⁠ data.frame from RMitemInfitCutoff directly: must have columns Item, infit_low, and infit_high (or outfit_low and outfit_high when statistic = "outfit"). When provided, adds columns Infit_low, Infit_high, and Flagged to the result. Flagged labels the misfit direction: "overfit" (infit below the range), "underfit" (above the range), or "" (within range, no misfit).
p_value	Logical or NULL. Whether to compute bootstrap p-values from the simulated null distribution and flag on them. NULL (the default) means use them when they are available: TRUE when cutoff is the full RMitemInfitCutoff object, which carries the per-item simulated values in its ⁠$results⁠ element, and FALSE when cutoff is NULL or the summarised ⁠$item_cutoffs⁠ data.frame, neither of which can support a p-value. Explicit TRUE with an insufficient cutoff is an error rather than a silent downgrade. With FALSE, items are flagged against the interval instead, and a one-time message reports the family-wise error rate that implies. The interval describes where a fitting item's statistic is expected to fall. Using it as a decision rule tests all k items at once, so its width sets a family-wise error rate of 1 - width^k, which is 37% for the default 95% interval over nine items. The corrected p-value targets alpha directly and needs far fewer iterations to do it (Johansson, 2026). Up to and including version 1.1.1 the default was FALSE, so Flagged came from the interval. Existing scripts that pass the full cutoff object will now flag on the corrected p-value instead and may report different items.
correction	Character. Multiple-comparison correction applied across items when p_value = TRUE: "fwer" (default) for the Westfall-Young studentised-max step-down (family-wise error rate), "fdr_bh" / "fdr_by" for Benjamini-Hochberg / Benjamini-Yekutieli false discovery rate control, or "none" for uncorrected per-item p-values.
alpha	Numeric in (0, 1). Significance level for the Flagged column when p_value = TRUE (an item is flagged when its corrected p-value is below alpha). Default 0.05.
output	Character string controlling the return value. Either "kable" (default) for a formatted knitr::kable() table, or "dataframe" for the underlying data.frame.
sort	Optional character string. When sort matches statistic ("infit" by default), rows are sorted by the reported MSQ column in descending order before output.
statistic	Character string. Which conditional fit statistic the table reports: "infit" (default) or "outfit". Both come from the same iarm::out_infit() call and the same parametric bootstrap, so the choice only selects which one is displayed and tested. Column names, cutoff bounds, p-values and Flagged all follow it, and when p_value = TRUE the multiplicity correction is applied across items within the selected statistic (the family is the k items, not 2k item-by-statistic combinations). Testing both statistics on the same items and reporting whichever flags is a larger family than either call corrects for.

Details

Infit MSQ is a weighted fit statistic that emphasises deviations near the item location. Values close to 1.0 indicate good fit. Values substantially above 1.0 suggest underfit (unexpected responses), while values substantially below 1.0 suggest overfit (overly predictable responses). The definition of "substantially" depends on several factors such as sample size, and needs to be determined by simulation using RMitemInfitCutoff. There is no general rule-of-thumb value that is correct.

Conditional infit MSQ statistics are computed via iarm::out_infit(), which uses the conditional distribution of the sufficient statistics (Müller, 2020). Only complete cases (rows without any NA) are used in the conditional fit calculation.

Item parameters are estimated by conditional maximum likelihood via psychotools::pcmodel() (a dichotomous item is a 2-category PCM); the conditional infit/outfit MSQ comes from iarm::out_infit() and is invariant to the estimation engine. Per-item average locations are the means of the CML thresholds, and the person-location reference is the mean of the Warm WLE estimates.

Relative item location is defined as the item's average location minus the sample mean person location, providing a measure of item targeting.

The iarm package must be installed (it is in Suggests, not Imports).

Bootstrap p-values. When p_value = TRUE, each item's observed infit (or outfit, per statistic) is compared against its simulated null distribution (from cutoff$results). The per-item statistic is the residual studentised by the bootstrap mean and SD – deliberately the empirical SD rather than the Wilson-Hilferty / ZSTD transform, which is uninformative for conditional MSQ (Müller, 2020). The marginal p-value is the two-sided Monte-Carlo p-value ⁠(1 + #{|t*| >= |t|}) / (B + 1)⁠. For correction = "fwer" the family-wise adjustment uses the Westfall-Young studentised-max step-down, which exploits the bootstrap dependence among items and is more powerful than Bonferroni/Holm (Ferreira, 2024); its validity rests on subset pivotality. "fdr_bh"/"fdr_by" apply Benjamini-Hochberg / Benjamini- Yekutieli instead. These are model-conditional, sample-size-sensitive p-values and are reported alongside the simulated effect-size band, not in place of it. p-values can be no smaller than 1 / (B + 1). Below 400 iterations in RMitemInfitCutoff() the studentised-max (FWER) correction is mildly liberal, because the bootstrap mean and SD used for studentisation are too noisy, and a message says so. From 400 iterations the family-wise error rate is at the nominal level, and more iterations improve reproducibility instead: two analysts using different seeds disagree about at least one item roughly 10 percent of the time at 400 iterations against 4 percent at 2000 (Johansson, 2026). The table caption notes this below 1000 iterations, so use 1000 to 2000 for a final analysis. The marginal p-values are well calibrated even at a few hundred iterations.

Value

• If output = "kable": a knitr_kable object (plain text table via format = "pipe") with columns "Item", "Infit MSQ", and "Relative location", and a caption showing the number of complete cases. When cutoff is provided, columns "Infit low", "Infit high", and "Flagged" are also included, and the caption notes the simulation-based cutoffs.

• If output = "dataframe": a data.frame with columns Item, Infit_MSQ, and Relative_location. When cutoff is provided, columns Infit_low, Infit_high, and Flagged are also included (inserted after Infit_MSQ, before Relative_location). Flagged is a character column with values "overfit", "underfit", or "" (not the previous logical). When p_value = TRUE, columns p_infit (marginal two-sided p-value) and padj_infit (corrected p-value) are added and Flagged reflects items with padj_infit < alpha (direction from Infit_MSQ relative to 1).

With statistic = "outfit" the layout is identical and the column names follow the statistic: Outfit_MSQ, Outfit_low, Outfit_high, p_outfit, padj_outfit.

Multiple comparisons

The marginal p-value controls the error rate of a single comparison: for one item (or item pair) decided on in advance it is the relevant value. But scanning all k comparisons and flagging whichever fall below alpha tests k hypotheses at once, so the chance of at least one false flag inflates to roughly $1 - (1 - \alpha)^k$ (e.g. about 34% for k = 8 at alpha = 0.05) – even when every marginal p-value is correctly calibrated. The corrected (adjusted) p-value controls this: correction = "fwer" bounds the probability of any false flag (strict, lower power), while "fdr_bh" / "fdr_by" bound the expected proportion of false flags among those raised (a more lenient middle ground). Rule of thumb: use the marginal p-value for a single pre-specified comparison, and a corrected p-value when screening the whole table – the usual workflow.

Flags depend on the other items

Item infit and item-restscore both judge each item against expectations from a Rasch model fitted to all items, so misfitting items change what is expected of the others. A flag in one direction can therefore produce flags in the other. Overfitting items can push fitting items into the underfit range, and underfitting items, for example from a second dimension, can make fitting items look overfit. How large the effect is depends on how strong the misfit is and how many items share it. In simulation with seven polytomous items, two overfitting items led to about 6 percent of the fitting items being flagged as underfit, and one or two mildly multidimensional items to 1 to 2 percent being flagged as overfit. When items are flagged in both directions, consider removing the clearest misfit and testing again before interpreting the rest.

References

Müller, M. (2020). Item fit statistics for Rasch analysis: Can we trust them? Journal of Statistical Distributions and Applications, 7(5). doi:10.1186/s40488-020-00108-7

Ferreira, J. A. (2024). Methods of testing a 'small' or 'moderate' number of hypotheses simultaneously. Journal of Statistical Theory and Practice, 19(6). doi:10.1007/s42519-024-00412-4

Westfall, P. H., & Young, S. S. (1993). Resampling-Based Multiple Testing. Wiley.

Johansson, M. (2026). Simulation-based cutoffs for conditional item fit in Rasch models: Iterations, multiplicity correction, and decision stability. PsyArXiv. doi:10.31234/osf.io/7pqz4_v2

See Also

RMitemInfitCutoff

Examples

if (requireNamespace("iarm", quietly = TRUE)) {
  # Simulate binary item response data (5 items, 40 persons)
  set.seed(42)
  sim_data <- as.data.frame(
    matrix(sample(0:1, 40 * 5, replace = TRUE), nrow = 40, ncol = 5)
  )
  colnames(sim_data) <- paste0("Item", 1:5)

  # Default kable output
  RMitemInfit(sim_data)

  # Sorted by infit MSQ descending
  RMitemInfit(sim_data, sort = "infit")

  # Return as data.frame for further processing
  df <- RMitemInfit(sim_data, output = "dataframe")

  # Simulation-based cutoffs (100 Monte-Carlo iterations)
  if (requireNamespace("ggdist", quietly = TRUE)) {
    cutoff_res <- RMitemInfitCutoff(sim_data, iterations = 100,
                                    parallel = FALSE, seed = 42)
    RMitemInfit(sim_data, cutoff = cutoff_res)
    RMitemInfit(sim_data, cutoff = cutoff_res, output = "dataframe")

    # Bootstrap p-values with family-wise (Westfall-Young) correction
    # (use iterations >= 1000 in real analyses for stable p-values)
    RMitemInfit(sim_data, cutoff = cutoff_res, p_value = TRUE,
                output = "dataframe")

    # The same for conditional outfit MSQ
    RMitemInfit(sim_data, cutoff = cutoff_res, p_value = TRUE,
                statistic = "outfit", output = "dataframe")
  }
}

---

Simulation-Based Infit MSQ Cutoff Determination

Description

Uses parametric bootstrap simulation to determine appropriate cutoff values for RMitemInfit. This function simulates data from a correctly fitting Rasch model that mimics your data and returns per-item empirical cutoffs.

Usage

RMitemInfitCutoff(
  data,
  iterations = 400,
  parallel = TRUE,
  n_cores = NULL,
  verbose = FALSE,
  seed = NULL,
  cutoff_method = "hdci",
  hdci_width = 0.95,
  dgp = c("resample", "conditional")
)

Arguments

data	A data.frame or matrix of item responses. Items must be scored starting at 0 (non-negative integers). Only complete cases (rows without any NA) are used.
iterations	Integer. Number of simulation iterations (default 400).
parallel	Logical. Use parallel processing via mirai if available (default TRUE).
n_cores	Integer or NULL. Number of parallel workers. When NULL, getOption("mc.cores") is checked first. If neither is set and parallel = TRUE, a warning is issued and execution falls back to sequential (single core) processing.
verbose	Logical. Show a progress bar (default FALSE).
seed	Integer or NULL. Random seed for reproducibility. See easyRasch2-reproducibility for what this guarantees and how it interacts with parallel.
cutoff_method	Character string specifying how cutoff intervals are computed. Either "hdci" (default) for the Highest Density Interval via ggdist::hdci(), or "quantile" for the 2.5th/97.5th percentiles via stats::quantile().
hdci_width	Numeric. Width of the HDCI when cutoff_method = "hdci". Default is 0.95 (95% HDCI). Ignored when cutoff_method = "quantile". The interval is a description of where a fitting item's statistic is expected to fall, not a decision rule. Flagging every item outside a width-w interval tests all k items at once, so the family-wise error rate is 1 - w^k, which is 37% for a 95% interval over nine items. Decisions should come from the corrected p-value instead (RMitemInfit with p_value = NULL and the full object returned here). The default was 0.999 up to and including version 1.1.1, which needs roughly 5000 iterations before the interval reaches its stated width; 0.95 reaches it by about 1000, so the band shown to readers means close to what it says (Johansson, 2026).
dgp	Character. Data-generating process for the parametric bootstrap. "resample" (default) resamples WLE person locations with replacement and simulates responses under the model (a marginal null). "conditional" simulates each respondent's pattern from the exact Rasch conditional distribution given their observed total score, item parameters fixed (a conditional null). Because the conditional infit/outfit statistic is itself conditional on the total score, "conditional" is its naturally matched null. Experimental.

Details

The generating model is CML item parameters (via psychotools) with WLE person locations. For each iteration a dataset is simulated under the chosen dgp, the model is refitted by CML (psychotools::pcmodel(), which handles dichotomous and polytomous data and is accepted by iarm), and conditional infit and outfit MSQ are computed via iarm::out_infit(). The distribution of these statistics across iterations provides empirical critical values per item. Failed iterations (e.g., degenerate simulated data) are silently discarded.

Parallel processing is provided by the mirai package (optional). Install it with install.packages("mirai") to enable parallelisation.

The iarm package must be installed (it is in Suggests, not Imports).

Value

A list with components:

results

data.frame with columns iteration, Item, InfitMSQ, OutfitMSQ (one row per item per successful iteration).

item_cutoffs

data.frame with per-item cutoff summaries: Item, infit_low, infit_high, outfit_low, outfit_high. Bounds are computed using the method specified by cutoff_method.

actual_iterations

Number of successful iterations. Everything downstream rests on this rather than on iterations, so it is the number to report.

requested_iterations

The iterations argument, kept so callers can tell how many simulated datasets were discarded.

sample_n

Number of complete cases used.

sample_n_total

Number of respondents in the raw input data, before the complete-case filter.

sample_has_na

Logical. Whether the raw input data contained any missing values.

sample_summary

Summary statistics of estimated person parameters.

item_names

Character vector of item names from data.

cutoff_method

The method used to compute cutoffs ("hdci" or "quantile").

hdci_width

The HDCI width used (only meaningful when cutoff_method = "hdci").

dgp

The data-generating process used ("resample" or "conditional").

References

Johansson, M. (2025). Detecting item misfit in Rasch models. Educational Methods & Psychometrics, 3(18). doi:10.61186/emp.2025.5

Johansson, M. (2026). Simulation-based cutoffs for conditional item fit in Rasch models: Iterations, multiplicity correction, and decision stability. PsyArXiv. doi:10.31234/osf.io/7pqz4_v2

See Also

RMitemInfit

Examples

if (requireNamespace("iarm", quietly = TRUE) &&
    requireNamespace("ggdist", quietly = TRUE)) {
  set.seed(42)
  sim_data <- as.data.frame(
    matrix(sample(0:1, 200 * 10, replace = TRUE), nrow = 200, ncol = 10)
  )
  colnames(sim_data) <- paste0("Item", 1:10)

  # Run 100 iterations sequentially for a quick demo
  cutoff_res <- RMitemInfitCutoff(sim_data, iterations = 100,
                                  parallel = FALSE, seed = 42)
  cutoff_res$item_cutoffs

  # Use the cutoffs in RMitemInfit()
  RMitemInfit(sim_data)
}

---

Simulation-Based Infit MSQ Cutoff Determination for Multiply Imputed Data

Description

Extends RMitemInfitCutoff to work with multiply imputed datasets produced by the mice package. Runs the parametric bootstrap simulation on each imputed dataset and stacks the resulting distributions, so that the final cutoff intervals reflect both sampling variability and imputation uncertainty.

Usage

RMitemInfitCutoffMI(
  mids_object,
  iterations = 500,
  parallel = TRUE,
  n_cores = NULL,
  verbose = FALSE,
  seed = NULL,
  cutoff_method = "hdci",
  hdci_width = 0.999
)

Arguments

mids_object	A mids object (multiply imputed dataset) as returned by mice::mice(). Each completed dataset must contain only the item response columns to be analysed (i.e., no ID or grouping variables). Items must be scored starting at 0 (non-negative integers).
iterations	Integer. Total number of simulation iterations to run across all imputations. These are distributed approximately evenly across the m imputed datasets (default 500).
parallel	Logical. Use parallel processing via mirai within each imputed dataset (default TRUE). Passed to RMitemInfitCutoff.
n_cores	Integer or NULL. Number of parallel workers. Passed to RMitemInfitCutoff.
verbose	Logical. Show progress messages (default FALSE).
seed	Integer or NULL. Master random seed for reproducibility. A unique per-imputation seed is derived from this value. See easyRasch2-reproducibility for what this guarantees and how it interacts with parallel.
cutoff_method	Character string specifying how cutoff intervals are computed from the stacked distribution. Either "hdci" (default) for the Highest Density Interval via ggdist::hdci(), or "quantile" for the 2.5th/97.5th percentiles via stats::quantile().
hdci_width	Numeric. Width of the HDCI when cutoff_method = "hdci". Default is 0.999 (99.9% HDCI). Ignored when cutoff_method = "quantile".

Details

The function completes each of the m imputed datasets via mice::complete(), then calls RMitemInfitCutoff on each one. The total number of iterations is split approximately evenly across imputations (i.e., each imputed dataset receives ceiling(iterations / m) or floor(iterations / m) iterations). The per-imputation simulation results are stacked into a single distribution from which cutoff intervals are computed, naturally incorporating imputation uncertainty.

Imputed datasets that cause model convergence failures are dropped with a warning. If all imputations fail, the function stops with an error.

The mice package must be installed (it is in Suggests, not Imports).

Value

A list with the same structure as RMitemInfitCutoff, so that the result can be passed directly to RMitemInfit, RMitemInfitMI, and RMitemInfitPlot:

results

data.frame with columns iteration, imputation, Item, InfitMSQ, OutfitMSQ — the stacked simulation results from all imputed datasets.

item_cutoffs

data.frame with per-item cutoff summaries: Item, infit_low, infit_high, outfit_low, outfit_high. Computed from the stacked distribution.

actual_iterations

Total number of successful iterations across all imputations.

sample_n

Number of rows (respondents) per imputed dataset.

sample_summary

Summary statistics of estimated person parameters from the first imputed dataset.

item_names

Character vector of item names.

cutoff_method

The method used to compute cutoffs.

hdci_width

The HDCI width used.

n_imputations

Number of imputed datasets used.

iterations_per_imputation

Integer vector of requested iterations per imputed dataset.

actual_iterations_per_imputation

Integer vector of successful iterations per imputed dataset.

Interval width under multiple imputation

hdci_width defaults to 0.999 here, not to the 0.95 used by RMitemInfitCutoff. RMitemInfitMI has no corrected-p-value path yet, so under imputation the interval is still the decision rule and keeps a width chosen for that job. Combining bootstrap p-values across imputations is planned but needs its own calibration, since Johansson (2026) covers complete data only.

See Also

RMitemInfitCutoff, RMitemInfitMI, RMitemInfitPlot

Examples

if (requireNamespace("mice", quietly = TRUE) &&
    requireNamespace("iarm", quietly = TRUE) &&
    requireNamespace("ggdist", quietly = TRUE)) {
  # Create example data with ~10% MCAR missingness
  set.seed(42)
  mat <- matrix(sample(0:1, 200 * 8, replace = TRUE), nrow = 200, ncol = 8)
  mat[sample(length(mat), round(0.10 * length(mat)))] <- NA
  sim_data <- as.data.frame(mat)
  colnames(sim_data) <- paste0("Item", 1:8)

  # mice's ordinal method (`polr`) requires the items to be ordered
  # factors, so code them as such before imputing. RMitemInfitCutoffMI()
  # converts the completed factors back to numeric internally.
  sim_data[] <- lapply(sim_data, function(x) factor(x, ordered = TRUE))

  # Impute (use more imputations, e.g. m = 5+, in real analyses)
  imp <- mice::mice(sim_data, m = 2, method = "polr", seed = 123,
                    printFlag = FALSE)

  # Compute simulation-based cutoffs across imputations
  # (use more iterations, e.g. 250+, in real analyses)
  cutoff_mi <- RMitemInfitCutoffMI(imp, iterations = 50, parallel = FALSE,
                                seed = 42)
  cutoff_mi$item_cutoffs

  # Use with RMitemInfitMI()
  RMitemInfitMI(imp, cutoff = cutoff_mi)
}

---

Conditional Item Infit MSQ for Multiply Imputed Data

Description

Extends RMitemInfit to work with multiply imputed datasets produced by the mice package. Computes conditional infit MSQ on each imputed dataset and pools the results using Rubin's rules.

Usage

RMitemInfitMI(mids_object, cutoff = NULL, output = "kable", sort)

Arguments

mids_object	A mids object (multiply imputed dataset) as returned by mice::mice(). Each completed dataset must contain only the item response columns to be analysed (i.e., no ID or grouping variables). Items must be scored starting at 0 (non-negative integers).
cutoff	Optional. Default NULL (no cutoff applied). Can be: The return value of RMitemInfitCutoff or RMitemInfitCutoffMI (a list with ⁠$item_cutoffs⁠): the data.frame is extracted automatically and metadata is included in the kable caption. The ⁠$item_cutoffs⁠ data.frame directly: must have columns Item, infit_low, and infit_high. When provided, adds columns Infit_low, Infit_high, and Flagged to the result. Flagged is a character column labelling the misfit direction: "overfit" (pooled infit below the range), "underfit" (above), or "" (within range).
output	Character string controlling the return value. Either "kable" (default) for a formatted knitr::kable() table, or "dataframe" for the underlying data.frame.
sort	Optional character string. When sort = "infit", rows are sorted by Infit_MSQ in descending order before output.

Details

For each of the m imputed datasets, the function:

1. Fits a Rasch model by CML via psychotools::pcmodel() (a dichotomous item is a 2-category partial credit model), consistent with RMitemInfit and the rest of the package.

2. Computes conditional infit MSQ and its standard error via iarm::out_infit().

3. Computes item locations (mean of the grand-mean-centred CML Andrich thresholds) and the mean WLE person location.

The per-imputation estimates are then pooled using Rubin's rules:

Pooled MSQ

The mean of the m infit MSQ point estimates.

Within-imputation variance

The mean of the m squared standard errors.

Between-imputation variance

The sample variance of the m point estimates.

Total variance

Within + (1 + 1/m) * Between.

Pooled SE

The square root of the total variance.

Relative item location is the mean of per-imputation relative locations (item location minus sample mean person location).

Caveat on the pooled SE. The within-imputation variance is the squared conditional infit SE from iarm::out_infit(). Müller (2020) showed that this asymptotic SE is an unreliable measure of uncertainty for the conditional infit statistic; Rubin's pooled SE inherits that limitation, so the Infit_SE/⁠Infit SE⁠ column should be read as an approximate indication of imputation-related variability rather than a trustworthy inferential standard error. For item misfit decisions, prefer the simulation-based cutoffs from RMitemInfitCutoffMI.

Imputed datasets that cause model convergence failures are dropped with a warning. If all imputations fail, the function stops with an error. At least two successful imputations are required to estimate between-imputation variance.

The mice and iarm packages must be installed (they are in Suggests, not Imports).

Value

• If output = "kable": a knitr_kable object (plain text table via format = "pipe") with columns "Item", "Infit MSQ", "Infit SE", "Relative location", and a caption noting the number of imputations and complete cases. When cutoff is provided, columns "Infit low", "Infit high", and "Flagged" are also included.

• If output = "dataframe": a data.frame with columns Item, Infit_MSQ, Infit_SE, and Relative_location. When cutoff is provided, columns Infit_low, Infit_high, and Flagged are also included (inserted after Infit_SE, before Relative_location). Flagged is a character column ("overfit" / "underfit" / ""), not the previous logical.

Flagging differs from the complete-data function

RMitemInfit flags on the Westfall-Young corrected p-value by default. This function has no p-value path and flags against the interval, so the two are not directly comparable. Combining bootstrap p-values across imputations is planned but needs calibration of its own, since Johansson (2026) studied complete data. Until then RMitemInfitCutoffMI keeps hdci_width = 0.999, a width suited to a decision rule, and the family-wise error rate of that rule is 1 - 0.999^k over k items.

References

Müller, M. (2020). Item fit statistics for Rasch analysis: Can we trust them? Journal of Statistical Distributions and Applications, 7(5). doi:10.1186/s40488-020-00108-7

See Also

RMitemInfit, RMitemInfitCutoffMI

Examples

if (requireNamespace("mice", quietly = TRUE) &&
    requireNamespace("iarm", quietly = TRUE) &&
    requireNamespace("ggdist", quietly = TRUE)) {
  # Create example data with ~10% MCAR missingness
  set.seed(42)
  mat <- matrix(sample(0:1, 200 * 8, replace = TRUE), nrow = 200, ncol = 8)
  mat[sample(length(mat), round(0.10 * length(mat)))] <- NA
  sim_data <- as.data.frame(mat)
  colnames(sim_data) <- paste0("Item", 1:8)

  # mice's ordinal method (`polr`) requires the items to be ordered
  # factors, so code them as such before imputing. RMitemInfitMI()
  # converts the completed factors back to numeric internally.
  sim_data[] <- lapply(sim_data, function(x) factor(x, ordered = TRUE))

  # Impute (use more imputations, e.g. m = 5+, in real analyses)
  imp <- mice::mice(sim_data, m = 2, method = "polr", seed = 123,
                    printFlag = FALSE)

  # Pooled infit table (no cutoffs)
  RMitemInfitMI(imp)

  # With simulation-based cutoffs
  # (use more iterations, e.g. 250+, in real analyses)
  cutoff_mi <- RMitemInfitCutoffMI(imp, iterations = 50, parallel = FALSE,
                                seed = 42)
  RMitemInfitMI(imp, cutoff = cutoff_mi)

  # As data.frame
  df <- RMitemInfitMI(imp, cutoff = cutoff_mi, output = "dataframe")
}

---

Plot Distribution of Simulated Infit and Outfit MSQ Values

Description

Visualises the distribution of simulation-based conditional item fit values from RMitemInfitCutoff, optionally overlaying observed item fit from the original data.

Usage

RMitemInfitPlot(simfit, data, statistic = "infit")

RMitemInfitCutoffPlot(...)

Arguments

simfit	The return value of RMitemInfitCutoff (a list with components results, item_cutoffs, actual_iterations, sample_n, and item_names).
data	Optional. A data.frame or matrix of item responses for computing and overlaying observed conditional item fit values. Items must be scored starting at 0 (non-negative integers). When provided, the plot includes orange diamond markers for the observed infit/outfit MSQ alongside the simulated distribution, plus segment summaries from the cutoff intervals.
statistic	Character string. Either "infit" (default) to show only the infit panel, "outfit" to show only the outfit panel, or "both" to show infit and outfit side by side (requires the patchwork package when data is supplied).
...	Arguments passed on to RMitemInfitPlot().

Details

Uses ggdist::stat_dotsinterval() (when data is not supplied) or ggdist::stat_dots() (when data is supplied) with point_interval = "median_hdci". The outer .width follows simfit$hdci_width, so the shaded interval matches the one RMitemInfit() tabulates.

When data is not supplied, the function plots the simulated MSQ distributions as dot-interval plots using ggdist::stat_dotsinterval() with median and Highest Density Continuous Interval (HDCI) summaries, faceted by statistic (InfitMSQ / OutfitMSQ).

When data is supplied, the function:

1. Fits a Rasch / Partial Credit model by CML via psychotools::pcmodel() (a dichotomous item is a 2-category PCM) and computes observed conditional infit and outfit MSQ via iarm::out_infit().

2. Overlays observed fit values as orange diamond markers on the simulated distributions.

3. Shows per-item cutoff intervals (from simfit$item_cutoffs) as black line segments, with thicker segments for the 66% HDCI range and black dots for the median.

The ggplot2, ggdist, and optionally patchwork packages must be installed (they are in Suggests, not Imports).

RMitemInfitCutoffPlot() is a deprecated alias for RMitemInfitPlot(), retained for backward compatibility with code written against easyRasch2 0.8.0. It warns and forwards to RMitemInfitPlot().

Value

A ggplot object (or a patchwork object when statistic = "both" and data is supplied).

See Also

RMitemInfitCutoff, RMitemInfit

Examples

if (requireNamespace("iarm", quietly = TRUE) &&
    requireNamespace("ggdist", quietly = TRUE) &&
    requireNamespace("ggplot2", quietly = TRUE)) {
  set.seed(42)
  sim_data <- as.data.frame(
    matrix(sample(0:1, 200 * 10, replace = TRUE), nrow = 200, ncol = 10)
  )
  colnames(sim_data) <- paste0("Item", 1:10)

  # Run simulation
  cutoff_res <- RMitemInfitCutoff(sim_data, iterations = 100,
                                  parallel = FALSE, seed = 42)

  # Simulated distribution only (infit + outfit faceted)
  RMitemInfitPlot(cutoff_res)

  # With observed fit overlaid (infit only, the default)
  RMitemInfitPlot(cutoff_res, data = sim_data)

  # Both infit and outfit panels side by side
  if (requireNamespace("patchwork", quietly = TRUE)) {
    RMitemInfitPlot(cutoff_res, data = sim_data, statistic = "both")
  }
}

---

Item Parameters for a Rasch / Partial Credit Model

Description

Estimates item difficulty (dichotomous) or item-category threshold (polytomous) parameters and returns them in long or wide format, with optional standard errors and Wald confidence intervals. Item parameters are estimated by conditional maximum likelihood (CML, via psychotools) by default, with marginal maximum likelihood (MML, via mirt) available for sparse data where CML can be unstable.

Usage

RMitemParameters(
  data,
  estimator = c("CML", "MML"),
  format = c("long", "wide"),
  se = TRUE,
  ci_level = 0.95,
  center = TRUE,
  output = c("kable", "dataframe", "file"),
  filename = NULL
)

Arguments

data	A data.frame or matrix of item responses. Items must be scored starting at 0 (non-negative integers). Missing values (NA) are allowed; both estimators handle them.
estimator	Character. "CML" (default) estimates item parameters with psychotools::pcmodel() (a dichotomous item is a 2-category PCM); "MML" uses mirt::mirt() with itemtype = "Rasch". CML is preferred for Rasch measurement; MML can be more robust when data are sparse. When estimator = "CML" and sparse response categories are detected, a warning suggests switching to MML.
format	Character. "long" (default) returns one row per item (dichotomous) or per item-threshold (polytomous); "wide" returns one row per item with threshold columns t1, t2, ... plus a mean location column.
se	Logical. If TRUE (default), standard-error and confidence-interval columns are added.
ci_level	Numeric in (0, 1). Confidence level for the Wald interval (⁠estimate +/- z * SE⁠). Default 0.95.
center	Logical. If TRUE (default), all thresholds are shifted so their grand mean is zero, the usual Rasch identification. CML estimates are already centred; the shift mainly affects the MML path, keeping the two estimators on a common scale.
output	Character. "kable" (default) for a formatted knitr::kable() table, "dataframe" for the underlying data.frame, or "file" to write that data.frame to a CSV at filename (the data.frame is also returned invisibly).
filename	Character. Path to the CSV file to write when output = "file". Required in that case; ignored otherwise.

Details

Items are detected as dichotomous (maximum score 1) or polytomous (maximum score > 1), and the Rasch or Partial Credit model is chosen accordingly. Thresholds are reported as Andrich thresholds (the person locations at which adjacent response categories are equally probable) on the logit difficulty scale, matching RMtargeting().

Standard errors. For the CML path, threshold SEs are the square roots of the diagonal of the threshold-parameter covariance from psychotools::threshpar(vcov = TRUE). For the MML path, SEs come from the mirt parameter covariance, propagated by the delta method through the linear threshold map. Confidence intervals are Wald intervals and are symmetric on the logit scale.

Value

For output = "dataframe", a data.frame. In long format the columns are item, threshold (integer; 1 for dichotomous items), location, and – when se = TRUE – se, ci_lower, ci_upper. In wide format the columns are item, the threshold locations (t1, t2, ... or location for dichotomous items), and a mean location; when se = TRUE, matching se_t1, se_t2, ... columns are appended. For output = "kable", the same content as a knitr_kable object.

References

Andrich, D. (1978). A rating formulation for ordered response categories. Psychometrika, 43(4), 561-573. doi:10.1007/BF02293814

Mair, P., & Hatzinger, R. (2007). Extended Rasch modeling: The eRm package for the application of IRT models in R. Journal of Statistical Software, 20(9), 1-20. doi:10.18637/jss.v020.i09

See Also

RMpersonParameters(), RMscoreSE(), RMtargeting()

Examples

set.seed(1)
poly <- as.data.frame(
  matrix(sample(0:2, 250 * 5, replace = TRUE), nrow = 250, ncol = 5)
)
colnames(poly) <- paste0("Item", 1:5)

# Default: long-format kable with SE and 95% CI
RMitemParameters(poly)

# Wide format, point estimates only
RMitemParameters(poly, format = "wide", se = FALSE, output = "dataframe")

# Dichotomous data
dich <- as.data.frame(
  matrix(sample(0:1, 250 * 6, replace = TRUE), nrow = 250, ncol = 6)
)
colnames(dich) <- paste0("Item", 1:6)
RMitemParameters(dich, output = "dataframe")

# Write the parameter table to a CSV (also returned invisibly)
RMitemParameters(poly, output = "file",
                 filename = tempfile(fileext = ".csv"))

---

Item Restscore Analysis

Description

Computes observed and model-expected item-restscore correlations using iarm::item_restscore(), and enriches the output with the absolute difference between observed and expected values, item average locations, and item locations relative to the sample mean person location.

Usage

RMitemRestscore(
  data,
  cutoff = NULL,
  p_value = NULL,
  correction = c("fwer", "fdr_bh", "fdr_by", "none"),
  alpha = 0.05,
  output = "kable",
  sort,
  p_adj = "BH"
)

Arguments

data	A data.frame or matrix of item responses. Items must be scored starting at 0 (non-negative integers). Missing values (NA) are allowed, but at least one complete case (row with no NA) must be present.
cutoff	Optional. Default NULL, which tests each item against the asymptotic normal reference distribution from iarm::item_restscore(), adjusted by p_adj. Can be: The return value of RMitemRestscoreCutoff: the parametric bootstrap null distribution of the observed minus expected gamma. Adds the columns Diff_low and Diff_high (the descriptive interval), and, with p_value resolving to TRUE, p_restscore and padj_restscore. The ⁠$item_cutoffs⁠ data.frame from RMitemRestscoreCutoff directly (columns Item, diff_low, diff_high). Items are then flagged against the interval. With a cutoff, the asymptotic p_adjusted column is dropped and p_adj is ignored.
p_value	Logical or NULL. Whether to compute bootstrap p-values from the simulated null distribution and flag on them. NULL (the default) means use them when they are available: TRUE when cutoff is the full RMitemRestscoreCutoff object, which carries the per-item simulated values in its ⁠$results⁠ element, and FALSE otherwise. Explicit TRUE without the full object is an error rather than a silent downgrade. With FALSE and a cutoff, items are flagged against the interval instead, and a one-time message reports the family-wise error rate that implies.
correction	Character. Multiple-comparison correction applied across items when p_value = TRUE: "fwer" (default) for the Westfall-Young studentised-max step-down (family-wise error rate), "fdr_bh" / "fdr_by" for Benjamini-Hochberg / Benjamini-Yekutieli false discovery rate control, or "none" for uncorrected per-item p-values.
alpha	Numeric in (0, 1). Significance level for the Flagged column when p_value = TRUE. Default 0.05.
output	Character string controlling the return value. Either "kable" (default) for a formatted knitr::kable() table, or "dataframe" for the underlying data.frame.
sort	Optional character string. When sort = "diff", rows are sorted by the absolute magnitude of Difference in descending order, so that both over- and underfitting items appear near the top.
p_adj	Character string specifying the p-value adjustment method passed to iarm::item_restscore(). Default "BH" (Benjamini-Hochberg); use "none" for unadjusted p-values. Run ?stats::p.adjust for the list of available methods. Only used when cutoff = NULL.

Details

Item-restscore correlations using Goodman-Kruskal's gamma (Kreiner, 2011) measure the association between a person's score on a single item and their total score on the remaining items (the "restscore"). Under a correctly fitting Rasch model, observed and model-expected correlations should agree closely.

Item parameters are estimated by conditional maximum likelihood via psychotools::pcmodel() (a dichotomous item is a 2-category PCM); the item-restscore statistic itself comes from iarm::item_restscore() and is conditional on the total score, so it is invariant to the estimation engine. Per-item average locations are the means of the CML thresholds, and the person-location reference is the mean of the Warm WLE estimates.

Relative item location is defined as the item's average location minus the sample mean person location, providing a measure of item targeting.

The iarm package must be installed (it is in Suggests, not Imports).

The asymptotic p-value is miscalibrated. Without a cutoff, each item is tested with (observed - expected) / SE against the standard normal, where the SE is that of the observed gamma alone. The expected gamma is estimated from the same data and correlates with the observed one, so the SE is too large for the difference. The difference is also biased upwards in small samples. Under a true Rasch model the test therefore flags too many items as overfit in small samples and too few as underfit at every sample size, and the adjustment chosen in p_adj does not correct either. In simulation, 20 dichotomous items at n = 150 gave at least one BH-flagged item in 13 percent of datasets (30 percent when 1.5 logits off target), almost all of them overfit, while 9 polytomous items at n = 1000 gave about 1 percent. Pass the object from RMitemRestscoreCutoff() as cutoff for p-values from a parametric bootstrap null instead.

Bootstrap p-values. When p_value = TRUE, each item's observed Difference is compared against its simulated null distribution (from cutoff$results), studentised by the bootstrap mean and SD. Because the bootstrap mean is subtracted, the small-sample upward bias of the difference is removed, and because the bootstrap SD is used, the correlation between observed and expected gamma is accounted for. The marginal p-value is the two-sided Monte-Carlo p-value ⁠(1 + #{|t*| >= |t|}) / (B + 1)⁠, and correction = "fwer" applies the Westfall-Young studentised-max step-down across items (Ferreira, 2024).

The two directions do not get equal shares of alpha. Gamma is bounded at 1, so the null of Difference is left-skewed and a two-sided test on ⁠|t|⁠ rejects more often in the long (underfit) tail. In simulation under a true Rasch model the total rate was nominal, but the per-item marginal rate was about 3 percent for underfit and 2 percent for overfit with 20 dichotomous items at n = 150 and 1.5 logits off target, narrowing toward equal shares with polytomous items and larger samples. With correction = "fwer", over four such conditions, at least one item was flagged as underfit in 3.7 percent of datasets and as overfit in 1.6 percent, 5.0 percent in total. An equal-tailed test, with alpha/2 for each direction, was evaluated and not adopted: it detected 2 to 3 percentage points fewer underfitting items and no more overfitting ones.

Reproducibility. The bootstrap p-values depend on the simulated null, so two analyses of the same data with different seeds can disagree about items near the decision boundary. In simulation, in conditions chosen to include such items, two seeds disagreed about the flag of at least one item in about 10 percent of analyses at 400 iterations and about 6 percent at 1000 with correction = "fwer", and in 13 and 11 percent with correction = "fdr_bh". Set seed in RMitemRestscoreCutoff() for a reproducible analysis, and use 1000 or more iterations for a final one.

The direction in Flagged follows the studentised value, not the sign of Difference. At small sample sizes the null mean of Difference is slightly positive, so an item can be flagged "underfit" while its Difference is still just above zero. This is rare.

Value

• If output = "kable": a knitr_kable object (plain text table via format = "pipe") with columns for item name, observed and expected restscore correlations, the signed difference (observed minus expected), adjusted p-value, the Flagged misfit label, and item location relative to the sample mean person location.

• If output = "dataframe": a data.frame with columns Item, Observed, Expected, Difference, p_adjusted, Flagged, and Relative_location. Flagged is "overfit" (observed above expected, adj. p < .05), "underfit" (below, adj. p < .05), or "" (not flagged).

With a cutoff, p_adjusted is replaced by Diff_low and Diff_high, followed by p_restscore (marginal two-sided bootstrap p-value) and padj_restscore (corrected p-value) when p_value resolves to TRUE. Flagged then reflects padj_restscore < alpha, or the interval when p_value = FALSE.

The Difference column is signed (observed minus expected): positive values indicate that the item correlates more strongly with the rest-score than the Rasch model predicts (over-discrimination / overfit, often associated with local dependence), and negative values indicate weaker-than-expected association (under-discrimination / underfit, often associated with multidimensionality or noise).

Multiple comparisons

The marginal p-value controls the error rate of a single comparison: for one item (or item pair) decided on in advance it is the relevant value. But scanning all k comparisons and flagging whichever fall below alpha tests k hypotheses at once, so the chance of at least one false flag inflates to roughly $1 - (1 - \alpha)^k$ (e.g. about 34% for k = 8 at alpha = 0.05) – even when every marginal p-value is correctly calibrated. The corrected (adjusted) p-value controls this: correction = "fwer" bounds the probability of any false flag (strict, lower power), while "fdr_bh" / "fdr_by" bound the expected proportion of false flags among those raised (a more lenient middle ground). Rule of thumb: use the marginal p-value for a single pre-specified comparison, and a corrected p-value when screening the whole table – the usual workflow.

Flags depend on the other items

Item infit and item-restscore both judge each item against expectations from a Rasch model fitted to all items, so misfitting items change what is expected of the others. A flag in one direction can therefore produce flags in the other. Overfitting items can push fitting items into the underfit range, and underfitting items, for example from a second dimension, can make fitting items look overfit. How large the effect is depends on how strong the misfit is and how many items share it. In simulation with seven polytomous items, two overfitting items led to about 6 percent of the fitting items being flagged as underfit, and one or two mildly multidimensional items to 1 to 2 percent being flagged as overfit. When items are flagged in both directions, consider removing the clearest misfit and testing again before interpreting the rest.

References

Kreiner, S. (2011). A Note on Item–Restscore Association in Rasch Models. Applied Psychological Measurement, 35(7), 557–561. doi:10.1177/0146621611410227

Ferreira, J. A. (2024). Methods of testing a 'small' or 'moderate' number of hypotheses simultaneously. Journal of Statistical Theory and Practice, 19(6). doi:10.1007/s42519-024-00412-4

See Also

RMitemRestscoreCutoff, RMitemRestscorePlot

Examples

if (requireNamespace("iarm", quietly = TRUE)) {
  # Simulate binary item response data (8 items, 200 persons)
  set.seed(42)
  sim_data <- as.data.frame(
    matrix(sample(0:1, 200 * 8, replace = TRUE), nrow = 200, ncol = 8)
  )
  colnames(sim_data) <- paste0("Item", 1:8)

  # Default kable output
  RMitemRestscore(sim_data)

  # Sorted by absolute difference
  RMitemRestscore(sim_data, sort = "diff")

  # Return as data.frame for further processing
  df <- RMitemRestscore(sim_data, output = "dataframe")

  # Bootstrap null distribution, flagging on Westfall-Young corrected
  # p-values (use 1000 or more iterations in a final analysis)
  if (requireNamespace("ggdist", quietly = TRUE)) {
    cutoff_res <- RMitemRestscoreCutoff(sim_data, iterations = 100,
                                        parallel = FALSE, seed = 42)
    RMitemRestscore(sim_data, cutoff = cutoff_res)
  }
}

---

Bootstrap Item-Restscore Misfit Detection

Description

Non-parametric bootstrap of item-restscore fit using iarm::item_restscore(). For each iteration, a sample of size samplesize is drawn from data with replacement, the appropriate Rasch model is refitted, and item-restscore results are classified as "overfit", "underfit", or "no misfit" based on the BH-adjusted p-value (< .05) and the sign of expected - observed. The function returns the percentage of iterations in which each item is flagged.

Usage

RMitemRestscoreBoot(
  data,
  iterations = 200,
  samplesize = 600,
  parallel = TRUE,
  n_cores = NULL,
  cutoff = 5,
  verbose = FALSE,
  seed = NULL,
  output = "kable"
)

Arguments

data	A data.frame or matrix of item responses. Items must be scored starting at 0 (non-negative integers).
iterations	Integer. Number of bootstrap samples (default 200).
samplesize	Integer. Size of each bootstrap sample (default 600). Must not exceed nrow(data).
parallel	Logical. Use parallel processing via mirai if available (default TRUE).
n_cores	Integer or NULL. Number of parallel workers. When NULL, getOption("mc.cores") is checked first. If neither is set and parallel = TRUE, a warning is issued and execution falls back to sequential processing.
cutoff	Numeric. Items flagged in fewer than this percentage of iterations are excluded from the kable output (default 5). Has no effect when output = "dataframe" or output = "raw".
verbose	Logical. Show a progress bar (default FALSE).
seed	Integer or NULL. Random seed for reproducibility. See easyRasch2-reproducibility for what this guarantees and how it interacts with parallel.
output	Either "kable" (default) for a formatted knitr::kable() table, "dataframe" for the per-item summary data.frame, or "raw" for the per-iteration long data.frame (useful for custom plotting).

Details

Each iteration uses the asymptotic test from iarm::item_restscore() with Benjamini-Hochberg adjustment, so the result inherits that test's calibration at the subsample size. Under a true Rasch model it flags too many items as overfit in small or mistargeted samples of dichotomous items, too few items of any kind in large polytomous samples, and too few underfitting items at every sample size (see RMitemRestscore()). The percentage is therefore how often that test flags an item in samples of size samplesize, not the probability that the item misfits, and it changes with samplesize. A low underfit percentage is weak evidence of fit. For a test of item fit, use RMitemRestscoreCutoff() with RMitemRestscore(), whose bootstrap p-values are calibrated at the observed sample size. This function is best read as a description of how stable the asymptotic flags are across resamples.

The full-sample model is fitted by CML via psychotools::pcmodel() (a dichotomous item is a 2-category partial credit model), and item locations (mean of the grand-mean-centred CML Andrich thresholds) and the mean WLE person location are computed from it — consistent with RMitemRestscore and the rest of the package. Each bootstrap iteration draws a sample of size samplesize with replacement and refits via psychotools::pcmodel(..., hessian = FALSE) for speed; iarm::item_restscore() accepts the fitted model.

Conditional infit MSQ (computed once on the full sample via iarm::out_infit()) and relative item locations are reported alongside the bootstrap percentages for context. The item-restscore classification and the infit statistic are conditional and engine-invariant; only the relative-item location shifts slightly relative to the previous eRm implementation, because it now uses the WLE person mean (finite at extreme scores) rather than the eRm MLE mean.

Iterations that fail (e.g., due to convergence issues on a degenerate bootstrap sample) are silently discarded; the caption / actual_iterations reflects only successful runs.

Parallel processing is provided by the mirai package (optional). Install it with install.packages("mirai") to enable parallelisation.

The iarm package must be installed (it is in Suggests, not Imports).

Value

• If output = "kable": a knitr_kable object listing items flagged in more than cutoff% of iterations, with columns Item, Item-restscore result, % of iterations, Conditional MSQ infit, and Relative item location, and a caption noting iteration count, bootstrap size, and number of complete cases.

• If output = "dataframe": a data.frame with one row per item × classification combination (Item, item_restscore, n, percent, Infit_MSQ, Relative_location), including "no misfit" rows.

• If output = "raw": a long data.frame with one row per item × successful iteration (iteration, Item, item_restscore, diff, diff_abs), where diff = expected - observed.

References

Kreiner, S. (2011). A Note on Item-Restscore Association in Rasch Models. Applied Psychological Measurement, 35(7), 557-561. doi:10.1177/0146621611410227

See Also

RMitemRestscore, RMitemInfit

Examples

if (requireNamespace("iarm", quietly = TRUE)) {
set.seed(42)
sim_data <- as.data.frame(
  matrix(sample(0:1, 400 * 8, replace = TRUE), nrow = 400, ncol = 8)
)
colnames(sim_data) <- paste0("Item", 1:8)

# Few iterations for a fast example; use 100+ in real analyses
# Default kable output (only items flagged > cutoff%)
RMitemRestscoreBoot(sim_data, iterations = 50, samplesize = 300,
                parallel = FALSE, seed = 1)

# Per-item summary data.frame (all classifications, including "no misfit")
summary_df <- RMitemRestscoreBoot(sim_data, iterations = 50, samplesize = 300,
                              parallel = FALSE, seed = 1,
                              output = "dataframe")

# Per-iteration long data for custom plotting
raw_df <- RMitemRestscoreBoot(sim_data, iterations = 50, samplesize = 300,
                          parallel = FALSE, seed = 1, output = "raw")

# Distribution of (expected - observed) across iterations, per item
if (requireNamespace("ggplot2", quietly = TRUE)) {
  library(ggplot2)
  ggplot(raw_df, aes(x = Item, y = diff)) +
    geom_hline(yintercept = 0, linetype = "dashed", colour = "grey50") +
    geom_violin(fill = "grey90", colour = NA) +
    geom_jitter(aes(colour = item_restscore),
                width = 0.15, alpha = 0.5, size = 0.8) +
    scale_colour_manual(values = c("overfit"   = "#377eb8",
                                   "underfit"  = "#e41a1c",
                                   "no misfit" = "grey60")) +
    labs(y = "Expected - observed restscore correlation",
         colour = NULL) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
}
}

---

Simulation-Based Item-Restscore Null Distribution

Description

Uses parametric bootstrap simulation to build the null distribution of the item-restscore statistic for RMitemRestscore. This function simulates data from a correctly fitting Rasch model that mimics your data and returns, per item, the simulated difference between observed and expected item-restscore gamma.

Usage

RMitemRestscoreCutoff(
  data,
  iterations = 400,
  parallel = TRUE,
  n_cores = NULL,
  verbose = FALSE,
  seed = NULL,
  cutoff_method = "hdci",
  hdci_width = 0.95,
  dgp = c("conditional", "resample")
)

Arguments

data	A data.frame or matrix of item responses. Items must be scored starting at 0 (non-negative integers). Only complete cases (rows without any NA) are used.
iterations	Integer. Number of simulation iterations (default 400).
parallel	Logical. Use parallel processing via mirai if available (default TRUE).
n_cores	Integer or NULL. Number of parallel workers. When NULL, getOption("mc.cores") is checked first. If neither is set and parallel = TRUE, a warning is issued and execution falls back to sequential (single core) processing.
verbose	Logical. Show a progress bar (default FALSE).
seed	Integer or NULL. Random seed for reproducibility. See easyRasch2-reproducibility for what this guarantees and how it interacts with parallel.
cutoff_method	Character string specifying how the intervals are computed. Either "hdci" (default) for the Highest Density Interval via ggdist::hdci(), or "quantile" for the 2.5th/97.5th percentiles via stats::quantile().
hdci_width	Numeric. Width of the HDCI when cutoff_method = "hdci". Default is 0.95 (95% HDCI). Ignored when cutoff_method = "quantile". The interval is a description of where a fitting item's difference is expected to fall, not a decision rule. Flagging every item outside a width-w interval tests all k items at once, so the family-wise error rate is 1 - w^k. Decisions should come from the corrected p-value instead (RMitemRestscore with p_value = NULL and the full object returned here).
dgp	Character. Data-generating process for the parametric bootstrap. "conditional" (default) simulates each respondent's pattern from the exact Rasch conditional distribution given their observed total score, item parameters fixed (a conditional null). The expected gamma is computed from the observed score distribution, which the conditional null holds fixed. "resample" resamples WLE person locations with replacement and simulates responses under the model (a marginal null). In simulation under a true Rasch model, the conditional null gave a family-wise rate of 5.0 percent with correction = "fwer" and the resample null 6.4 percent, pooled over four designs. The conditional null takes about 1.2 to 1.6 times as long. Experimental.

Details

The asymptotic test in iarm::item_restscore() divides the observed minus expected gamma by the standard error of the observed gamma alone. The expected gamma is estimated from the same data and correlates with the observed one, so that standard error is too large for the difference, and the difference is also biased upwards in small samples. Under a true Rasch model the resulting test is liberal for dichotomous items in small or mistargeted samples and conservative for polytomous items, in every case flagging too few underfitting items. The bootstrap null replaces the asymptotic reference distribution and absorbs both problems.

The generating model is CML item parameters (via psychotools) with WLE person locations. For each iteration a dataset is simulated under the chosen dgp, the model is refitted by CML (psychotools::pcmodel()), and the observed and expected item-restscore gamma are computed as in iarm::item_restscore(), by a faster internal routine that skips the standard errors and the rounding of the printed values. The refit matters: the expected gamma varies from sample to sample because the thresholds do, and holding them fixed would reproduce the problem the bootstrap exists to solve. Failed iterations (e.g., degenerate simulated data) are silently discarded.

Parallel processing is provided by the mirai package (optional). Install it with install.packages("mirai") to enable parallelisation.

The iarm package must be installed (it is in Suggests, not Imports).

Value

A list with components:

results

data.frame with columns iteration, Item, Observed, Expected, Difference (one row per item per successful iteration). Difference is Observed - Expected, the statistic RMitemRestscore tests.

item_cutoffs

data.frame with per-item interval bounds for the difference: Item, diff_low, diff_high. Bounds are computed using the method specified by cutoff_method.

actual_iterations

Number of successful iterations. Everything downstream rests on this rather than on iterations, so it is the number to report.

requested_iterations

The iterations argument, kept so callers can tell how many simulated datasets were discarded.

sample_n

Number of complete cases used.

sample_n_total

Number of respondents in the raw input data, before the complete-case filter.

sample_has_na

Logical. Whether the raw input data contained any missing values.

sample_summary

Summary statistics of estimated person parameters.

item_names

Character vector of item names from data.

cutoff_method

The method used to compute the intervals ("hdci" or "quantile").

hdci_width

The HDCI width used (only meaningful when cutoff_method = "hdci").

dgp

The data-generating process used ("resample" or "conditional").

References

Kreiner, S. (2011). A Note on Item-Restscore Association in Rasch Models. Applied Psychological Measurement, 35(7), 557-561. doi:10.1177/0146621611410227

Johansson, M. (2026). Simulation-based cutoffs for conditional item fit in Rasch models: Iterations, multiplicity correction, and decision stability. PsyArXiv. doi:10.31234/osf.io/7pqz4_v2

See Also

RMitemRestscore, RMitemRestscorePlot

Examples

if (requireNamespace("iarm", quietly = TRUE) &&
    requireNamespace("ggdist", quietly = TRUE)) {
  set.seed(42)
  sim_data <- as.data.frame(
    matrix(sample(0:1, 200 * 10, replace = TRUE), nrow = 200, ncol = 10)
  )
  colnames(sim_data) <- paste0("Item", 1:10)

  # Run 100 iterations sequentially for a quick demo
  cutoff_res <- RMitemRestscoreCutoff(sim_data, iterations = 100,
                                      parallel = FALSE, seed = 42)
  cutoff_res$item_cutoffs

  # Flag on bootstrap p-values in RMitemRestscore()
  RMitemRestscore(sim_data, cutoff = cutoff_res)
}

---