Package {easyRasch2}

Title:	Psychometric Analysis with Rasch Measurement Theory
Version:	1.0.0
Description:	Streamlines reproducible Rasch measurement theory analyses for ordinal item-response data, combining estimation routines from 'eRm', 'psychotools', 'mirt', 'iarm', and 'lavaan' with consistent diagnostic, plotting, and reporting layers. Covers the four basic psychometric criteria summarised by Christensen et al. (2021) <doi:10.1111/sms.13908> – unidimensionality, local independence, ordered response category thresholds, and invariance across subgroups – together with item fit, targeting, reliability, category functioning, and descriptive item-response plots. A distinguishing feature is the use of simulation-based critical values to replace rule-of-thumb cutoffs for conditional infit mean-square, Yen's Q3 local-dependence statistic, the largest residual-PCA eigenvalue, ordinal CFA fit indices, and partial-gamma DIF and local-dependence coefficients, optionally augmented with multiplicity-corrected bootstrap p-values. Outputs are knitr::kable() tables and 'ggplot2' figures suitable for direct inclusion in 'Quarto' and 'R Markdown' reports.
License:	GPL (≥ 3)
Encoding:	UTF-8
LazyData:	true
RoxygenNote:	7.3.3
URL:	https://github.com/pgmj/easyRasch2, https://pgmj.github.io/easyRasch2/
BugReports:	https://github.com/pgmj/easyRasch2/issues
Depends:	R (≥ 4.1.0)
Imports:	graphics, knitr, mirt, psychotools (≥ 0.7-3), stats, utils, rlang
Suggests:	difR, eRm, geomtextpath, ggdist, ggtext, iarm, mirai, ggplot2 (≥ 3.4.0), partykit, psychotree, stablelearner, testthat (≥ 3.0.0), rmarkdown, patchwork, scales, mice, ggrepel, lavaan, withr
Config/testthat/edition:	3
VignetteBuilder:	knitr
NeedsCompilation:	no
Packaged:	2026-07-05 12:52:56 UTC; magnus.johansson.3
Author:	Magnus Johansson [aut, cre], Nicklas Korsell [ctb] (PCM simulation code), Mirka Henninger [ctb] (MH / partial-gamma effect-size and ETS-classification algorithms in dif_tree.R, adapted under MIT licence from the raschtreeMH and effecttree packages), Jan Radek [ctb] (partial-gamma effect-size and ETS-classification algorithms in dif_tree.R, adapted under MIT licence from the effecttree package)
Maintainer:	Magnus Johansson <pgmj@pm.me>
Repository:	CRAN
Date/Publication:	2026-07-05 15:50:08 UTC

easyRasch2: Psychometric Analysis with Rasch Measurement Theory

Description

Streamlines reproducible Rasch measurement theory analyses for ordinal item-response data, combining estimation routines from 'eRm', 'psychotools', 'mirt', 'iarm', and 'lavaan' with consistent diagnostic, plotting, and reporting layers. Covers the four basic psychometric criteria summarised by Christensen et al. (2021) doi:10.1111/sms.13908 – unidimensionality, local independence, ordered response category thresholds, and invariance across subgroups – together with item fit, targeting, reliability, category functioning, and descriptive item-response plots. A distinguishing feature is the use of simulation-based critical values to replace rule-of-thumb cutoffs for conditional infit mean-square, Yen's Q3 local-dependence statistic, the largest residual-PCA eigenvalue, ordinal CFA fit indices, and partial-gamma DIF and local-dependence coefficients, optionally augmented with multiplicity-corrected bootstrap p-values. Outputs are knitr::kable() tables and 'ggplot2' figures suitable for direct inclusion in 'Quarto' and 'R Markdown' reports.

Author(s)

Maintainer: Magnus Johansson pgmj@pm.me (ORCID)

Other contributors:

• Nicklas Korsell (PCM simulation code) [contributor]

• Mirka Henninger (ORCID) (MH / partial-gamma effect-size and ETS-classification algorithms in dif_tree.R, adapted under MIT licence from the raschtreeMH and effecttree packages) [contributor]

• Jan Radek (ORCID) (partial-gamma effect-size and ETS-classification algorithms in dif_tree.R, adapted under MIT licence from the effecttree package) [contributor]

See Also

Useful links:

• https://github.com/pgmj/easyRasch2

• https://pgmj.github.io/easyRasch2/

• Report bugs at https://github.com/pgmj/easyRasch2/issues

Format a human-readable label for the cutoff method

Description

Format a human-readable label for the cutoff method

Usage

.format_cutoff_method_label(cutoff_method, hdci_width)

Arguments

Value

A character label, or NULL if the method is unknown/unset.

Format a human-readable label for the gamma cutoff method

Description

Format a human-readable label for the gamma cutoff method

Usage

.format_gamma_cutoff_method_label(cutoff_method, hdci_width)

Arguments

Value

A character label, or NULL if the method is unknown/unset.

Relative Measurement Uncertainty (RMU)

Description

Bayesian-style reliability estimate (Bignardi, Kievit & Bürkner, 2025) computed from a matrix of posterior or plausible-value draws. The columns of input_draws are split at random into two halves; reliability is the Pearson correlation across persons of paired columns from the two halves, summarised across pairs as a posterior mean with HDCI.

Usage

RMUreliability(input_draws, level = 0.95, verbose = FALSE)

Arguments

Details

Adapted (with permission, GPL-2/3) from https://github.com/giac01/gbtoolbox/blob/main/R/reliability.R.

The function silently returns 0 for any column pair where either side has zero variance (the correlation is undefined there).

Requires the ggdist package (Suggests).

Value

A 1-row data.frame with columns rmu_estimate, hdci_lowerbound, hdci_upperbound, plus the .width/.point/.interval metadata columns added by ggdist::mean_hdci().

References

Bignardi, G., Kievit, R., & Bürkner, P. C. (2025). A general method for estimating reliability using Bayesian Measurement Uncertainty. PsyArXiv. doi:10.31234/osf.io/h54k8_v1

See Also

RMreliability()

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 = FALSE,
  correction = c("fwer", "fdr_bh", "fdr_by", "none"),
  alpha = 0.05,
  output = "kable"
)

Arguments

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).

Bootstrap p-values. When p_value = TRUE, each item's observed partial gamma is compared against its simulated null distribution (from cutoff$results, where the DIF variable is random 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). correction = "fwer" uses the Westfall-Young studentised-max step-down, which exploits the bootstrap dependence among items (Ferreira, 2024); it is liberal when the simulation is small, so at least 1000 iterations in RMdifGammaCutoff() are recommended (a warning 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.

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." (a star-string indicator from iarm::partgam_DIF()). 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.

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)

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

Simulation-Based Partial Gamma DIF Cutoff Determination

Description

Uses parametric bootstrap simulation to determine appropriate cutoff values for partial gamma DIF analysis via partgam_DIF. Under a correctly fitting Rasch model where the DIF variable is unrelated to item responses (i.e., no true DIF), this function generates the expected distribution of absolute partial gamma values per item, providing empirical critical values.

Usage

RMdifGammaCutoff(
  data,
  dif_var,
  iterations = 250,
  parallel = TRUE,
  n_cores = NULL,
  verbose = FALSE,
  seed = NULL,
  cutoff_method = "hdci",
  hdci_width = 0.99
)

Arguments

Details

For each simulation iteration the function:

1. Resamples person parameters (thetas) with replacement from the WLE person locations.

2. Simulates item response data under a Rasch model (dichotomous via psychotools::rrm() or polytomous via an internal partial credit simulator).

3. Creates a random DIF variable by sampling group labels with the same proportions as the observed dif_var, so there is no true DIF by construction.

4. Computes partial gamma DIF statistics via iarm::partgam_DIF().

The distribution of partial gamma values across iterations provides empirical critical values per item. Values from real data that fall outside these bounds suggest DIF that exceeds what would be expected by chance under a correctly fitting Rasch model. Failed iterations (e.g., due to convergence issues or degenerate data) are silently discarded.

The generating model uses CML item thresholds via psychotools::pcmodel() (a dichotomous item is a 2-category PCM) and WLE person locations, consistent with the rest of the package; responses are simulated with psychotools::rrm() (dichotomous) or an internal partial credit score simulator (polytomous).

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, 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.

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 estimated person parameters.

item_names

Character vector of item names from data.

dif_group_sizes

Named integer vector of group sizes used in the simulation (matches proportions in the observed dif_var).

cutoff_method

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

hdci_width

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

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

See Also

partgam_DIF

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)
  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
}

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 via partgam_DIF.

Usage

RMdifGammaPlot(simfit, data, dif_var)

Arguments

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 via 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, ggdist, and optionally iarm packages must be installed (they are in Suggests, not Imports).

Value

A ggplot object.

See Also

RMdifGammaCutoff, RMdifGamma

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)
  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

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"),
  ...
)

Arguments

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.

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

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

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

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

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. For dichotomous data the fast algorithm of Christensen & Kreiner (2007, p. 23) is used (sample without replacement weighted by item easinesses). For polytomous data the recursive \gamma-function approach is used, with each item's response sampled conditional on the remaining items' joint score distribution (computed via psychotools::elementary_symmetric_functions()).

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.

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.

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.

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

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

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

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

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.

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., 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

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

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

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

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("cut", "score"),
  class_intervals = 4,
  ci = TRUE,
  error_band = FALSE,
  conf_level = 0.95,
  min_n = 8L,
  items = NULL,
  output = c("patchwork", "list")
)

Arguments

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. With method = "cut", total scores are split into class_intervals equal-count bins using common boundaries (so groups share the x-axis); each bin contributes one observed point at its mean total score. With method = "score", every observed total score is a point. If the score distribution is too sparse to form class_intervals distinct bins, the function falls back to score-level points.

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.

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

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 = FALSE,
  correction = c("fwer", "fdr_bh", "fdr_by", "none"),
  alpha = 0.05,
  output = "kable",
  sort
)

Arguments

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 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), and the studentised-max (FWER) correction is liberal when the simulation is small (the bootstrap mean/SD used for studentisation are then too noisy). At least 1000 iterations in RMitemInfitCutoff() are recommended – in simulations the family-wise error rate is then controlled at the nominal level – and a warning is issued when the simulation is smaller. 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).

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

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.

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")
  }
}

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 = 250,
  parallel = TRUE,
  n_cores = NULL,
  verbose = FALSE,
  seed = NULL,
  cutoff_method = "hdci",
  hdci_width = 0.999,
  dgp = c("resample", "conditional")
)

Arguments

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.

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").

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

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.

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

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.

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

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.999).

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

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, output = "kable", sort, p_adj = "BH")

Arguments

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).

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).

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).

References

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

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 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

Details

Useful with large samples, where the asymptotic test underlying RMitemRestscore can flag items that are not practically misfitting; bootstrapping gives a more nuanced view of the probability of an item actually being misfit.

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))
}
}

Partial Gamma Local Dependence Analysis

Description

Computes partial gamma coefficients for Local Dependence (LD) assessment using iarm::partgam_LD(). Each pair of items is tested for residual association, controlling for the rest score (total score minus one of the items in the pair).

Usage

RMlocdepGamma(
  data,
  cutoff = NULL,
  p_value = FALSE,
  correction = c("fwer", "fdr_bh", "fdr_by", "none"),
  alpha = 0.05,
  output = "kable",
  n_pairs = NULL
)

Arguments

Details

Partial gamma (Christensen, Kreiner & Mesbah, 2013) measures the residual association between pairs of items after controlling for the rest score (total score minus one item). Because it matters which item is subtracted, calculations are done for each pair in both directions, yielding two data.frames.

Values near 0 indicate no local dependence. Large positive values suggest positive LD (items share variance beyond the latent trait), while large negative values suggest negative LD.

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

Bootstrap p-values. When p_value = TRUE, each pair's observed partial gamma (canonical direction) is compared against its simulated null distribution (from cutoff$results, simulated under local independence). The per-pair statistic is the residual studentised by the bootstrap mean and SD; the marginal p-value is the one-sided Monte-Carlo p-value ⁠(1 + #\{t* >= t\}) / (B + 1)⁠ for excess positive LD (redundancy, the diagnostic target — matching RMlocdepQ3), so it can be no smaller than 1 / (B + 1). The band still shows both bounds for reference. correction = "fwer" uses the Westfall-Young studentised-max step-down over the family of all pairs, which exploits the bootstrap dependence among them (Ferreira, 2024); it is liberal when the simulation is small, so at least 1000 iterations in RMlocdepGammaCutoff() are recommended (a warning is issued below that). Unlike the asymptotic p-values from iarm::partgam_LD(), 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.

Value

• If output = "kable": an object of class "RMlocdepGamma". Internally a list with two knitr_kable elements, ⁠$direction1⁠ and ⁠$direction2⁠. In both, the rest score is the total score minus Item 2 (the second column); the two elements list each item pair in the two possible orders, so together they cover both rest-score directions for every pair. Each has columns "Item 1", "Item 2", "Partial gamma", "Adj. p-value (BH)", and "p-value sign." (a star-string indicator from iarm::partgam_LD()). When cutoff is provided, additional columns "Gamma low", "Gamma high", and "Flagged" are included.

  The object has custom print() and knitr::knit_print() methods: in the R console it prints the two tables stacked vertically; in a Quarto / R Markdown chunk it renders as two distinct pipe tables. Access the individual tables explicitly as result$direction1 and result$direction2 if needed.

• If output = "dataframe": a named list of two data.frames (⁠$direction1⁠, ⁠$direction2⁠) with columns Item1, Item2, gamma, 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 (identical for a pair in both directions).

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

Christensen, K. B., Kreiner, S. & Mesbah, M. (Eds.) (2013). Rasch Models in Health, pp. 133–135. ISTE & Wiley. doi:10.1002/9781118574454

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

RMlocdepGammaCutoff, RMlocdepGammaPlot

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)

  # Default kable output
  RMlocdepGamma(sim_data)

  # Return as data.frame list
  RMlocdepGamma(sim_data, output = "dataframe")

  # Simulation-based cutoffs (slow): 100+ Monte-Carlo iterations
  if (requireNamespace("ggdist", quietly = TRUE)) {
    cutoff_res <- RMlocdepGammaCutoff(sim_data, iterations = 100,
                                      parallel = FALSE, seed = 42)
    RMlocdepGamma(sim_data, cutoff = cutoff_res)

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

Simulation-Based Partial Gamma LD Cutoff Determination

Description

Uses parametric bootstrap simulation to determine appropriate cutoff values for partial gamma Local Dependence analysis via partgam_LD. Under a correctly fitting Rasch model where items are locally independent, this function generates the expected distribution of partial gamma values per item pair, providing empirical critical values.

Usage

RMlocdepGammaCutoff(
  data,
  iterations = 250,
  parallel = TRUE,
  n_cores = NULL,
  verbose = FALSE,
  seed = NULL,
  cutoff_method = "hdci",
  hdci_width = 0.99
)

Arguments

Details

For each simulation iteration the function:

1. Resamples person parameters (thetas) with replacement from the WLE person locations.

2. Simulates item response data under a Rasch model (dichotomous via psychotools::rrm() or polytomous via an internal partial credit simulator).

3. Computes partial gamma LD statistics via iarm::partgam_LD().

Because the data are simulated under the Rasch model, items are locally independent by construction. The distribution of partial gamma values across iterations provides empirical critical values per item pair. Values from real data that fall outside these bounds suggest local dependence that exceeds what would be expected by chance. Failed iterations (e.g., due to convergence issues or degenerate data) are silently discarded.

The generating model uses CML item thresholds via psychotools::pcmodel() (a dichotomous item is a 2-category PCM) and WLE person locations, consistent with the rest of the package; responses are simulated with psychotools::rrm() (dichotomous) or an internal partial credit score simulator (polytomous).

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, Item1, Item2, and gamma (one row per item pair per successful iteration). Contains results from direction 1 only (rest score = total - Item2), which is the conventional direction.

pair_cutoffs

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

actual_iterations

Number of successful 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.

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").

References

Christensen, K. B., Kreiner, S. & Mesbah, M. (Eds.) (2013). Rasch Models in Health, pp. 133–135. ISTE & Wiley. doi:10.1002/9781118574454

See Also

partgam_LD, RMlocdepGamma, RMlocdepGammaPlot

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 <- RMlocdepGammaCutoff(sim_data, iterations = 100,
                                    parallel = FALSE, seed = 42)
  cutoff_res$pair_cutoffs
}

Plot Distribution of Simulated Partial Gamma LD Values

Description

Visualises the distribution of simulation-based partial gamma LD values from RMlocdepGammaCutoff, optionally overlaying observed partial gamma values computed from real data via partgam_LD.

Usage

RMlocdepGammaPlot(simfit, data, items = NULL, n_pairs = NULL)

Arguments

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).

The plot shows one row per item pair (labelled as "Item1 - Item2"). Only direction 1 (rest score = total - Item2) is plotted, matching the convention used in the simulation.

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, the function:

1. Computes observed partial gamma values via iarm::partgam_LD().

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

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

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

Value

A ggplot object.

See Also

RMlocdepGammaCutoff, RMlocdepGamma

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 <- RMlocdepGammaCutoff(sim_data, iterations = 100,
                                    parallel = FALSE, seed = 42)

  # Simulated distribution only
  RMlocdepGammaPlot(cutoff_res)

  # With observed partial gamma overlaid
  RMlocdepGammaPlot(cutoff_res, data = sim_data)

  # Plot only a subset of items
  RMlocdepGammaPlot(cutoff_res, data = sim_data,
                    items = c("Item1", "Item2", "Item3"))
}

Q_3 Residual Correlations for Local Dependence Assessment

Description

Computes Yen's Q_3 residual correlations between item pairs. By default the Rasch model is fitted by conditional maximum likelihood with WLE person locations (estimator = "CML"); marginal ML via mirt is available with estimator = "MML". High correlations (above the dynamic cut-off) indicate potential local dependence between items. See RMlocdepQ3Cutoff for how to determine the appropriate dynamic cut-off for your data.

Usage

RMlocdepQ3(
  data,
  cutoff = NULL,
  output = "kable",
  n_pairs = NULL,
  p_value = FALSE,
  correction = c("fwer", "fdr_bh", "fdr_by", "none"),
  alpha = 0.05,
  estimator = c("CML", "MML")
)

Arguments

Details

The Q_3 statistic (Yen, 1984) is the correlation between residuals of pairs of items after accounting for the latent trait. Under local independence, Q_3 values are expected to be around -1/(k-1) where k is the number of items. When cutoff is supplied, the dynamic cut-off is the mean of all off-diagonal Q_3 values plus cutoff, following the approach of Christensen et al. (2017). Use RMlocdepQ3Cutoff to obtain a simulation-based cutoff recommendation.

Q_3 is the column-wise correlation matrix of the model standardized residuals (x - E)/\sqrt{Var}. By default (estimator = "CML") item parameters are estimated by conditional maximum likelihood (psychotools) and person locations by Warm's weighted likelihood; estimator = "MML" instead uses mirt's marginal-ML model and its built-in Q_3 residuals. The two estimators give very similar Q_3 values (off-diagonal correlation typically > 0.95); what matters for inference is that the observed Q_3 and the simulated cut-off in RMlocdepQ3Cutoff use the same estimator, which the functions enforce.

Two views of local dependence. Given the full RMlocdepQ3Cutoff object, two complementary tables are returned. The ⁠$matrix⁠ applies a single global cut-off (the Christensen et al. approach: the 99th percentile of the simulated max-minus-mean Q_3) – a family-wise "is there any local dependence" overview. The ⁠$pairs⁠ table is the per-comparison view: each observed Q_3 against its own simulated expected range (the Low/High bounds), so individual dependent pairs can be read off and ranked.

Bootstrap p-values. When p_value = TRUE, the ⁠$pairs⁠ table also tests each observed Q_3 against its simulated null (from cutoff$pair_results) with a one-sided (upper-tail) test for excess local dependence. The pair statistic is studentised by the bootstrap mean and SD; the marginal p-value is ⁠(1 + #{Q3* >= Q3}) / (B + 1)⁠, and correction applies the family-wise (Westfall-Young step-down) or FDR adjustment across the k(k-1)/2 pairs. As for item fit, the family-wise correction is liberal when the simulation is small, so >= 1000 iterations in RMlocdepQ3Cutoff() are recommended (a warning is issued otherwise).

Value

With cutoff = NULL or a bare numeric cut-off, a single object (kable or data.frame) holding the lower triangle of the Q_3 matrix; a numeric cut-off adds a Flagged row flag and a caption describing the dynamic cut-off.

With the full RMlocdepQ3Cutoff() object, a named list of two:

⁠$matrix⁠

the Q_3 lower-triangle matrix with the global dynamic cut-off (mean off-diagonal Q_3 + suggested cut-off), as above.

⁠$pairs⁠

one row per item pair: Item1, Item2, Observed (Q_3), Low/High (the per-pair expected range, i.e. the simulated bounds), and Flagged – "above" (Q_3 above the upper bound, indicating local dependence), "below" (below the lower bound), or "". Sorted by absolute departure from the per-pair simulated median and truncated to n_pairs. With p_value = TRUE, columns p_q3 and padj_q3 are added and Flagged reflects padj_q3 < alpha and only flags "above".

The Q_3 tile heatmap that earlier versions returned as ⁠$plot⁠ is now produced by RMlocdepQ3Plot (as its ⁠$matrix⁠ element), so the table and plot outputs share the same ⁠$matrix⁠/⁠$pairs⁠ structure.

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

Yen, W. M. (1984). Effects of local item dependence on the fit and equating performance of the three-parameter logistic model. Applied Psychological Measurement, 8(2), 125–145. doi:10.1177/014662168400800201

Christensen, K. B., Makransky, G., & Horton, M. (2017). Critical values for Yen's Q_3: Identification of local dependence in the Rasch model. Applied Psychological Measurement, 41(3), 178–194. doi:10.1177/0146621616677520

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

Examples

# Simulate binary item response data (10 items, 200 persons)
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)

# Raw Q3 matrix (no cutoff)
RMlocdepQ3(sim_data)

# Get the underlying data.frame
q3_df <- RMlocdepQ3(sim_data, output = "dataframe")

# Simulation-based cutoff (use 500+ iterations in real analyses)
if (requireNamespace("ggdist", quietly = TRUE)) {
  cutoff_res <- RMlocdepQ3Cutoff(sim_data, iterations = 50, parallel = FALSE)

  # Bare numeric cutoff -> just the matrix
  RMlocdepQ3(sim_data, cutoff = cutoff_res$suggested_cutoff)

  # Full object -> list of two tables: $matrix and $pairs
  res <- RMlocdepQ3(sim_data, cutoff = cutoff_res, output = "dataframe")
  res$pairs

  # Top 5 pairs, with bootstrap p-values (use iterations >= 1000 in practice)
  RMlocdepQ3(sim_data, cutoff = cutoff_res, n_pairs = 5, p_value = TRUE,
             output = "dataframe")$pairs
}

Simulation-Based Q_3 Cutoff Determination

Description

Uses parametric bootstrap simulation to determine an appropriate cutoff value for RMlocdepQ3. Under a correctly fitting Rasch model, Q_3 residuals have an unknown distribution; this function simulates that distribution and returns empirical percentiles.

Usage

RMlocdepQ3Cutoff(
  data,
  iterations = 500,
  parallel = TRUE,
  n_cores = NULL,
  verbose = FALSE,
  seed = NULL,
  cutoff_method = "hdci",
  hdci_width = 0.99,
  estimator = c("CML", "MML"),
  dgp = c("resample", "conditional")
)

Arguments

Details

The generating model is fitted once: CML item parameters (via psychotools) and WLE person locations. For each simulation iteration, those WLE thetas are resampled with replacement, response data are simulated under the Rasch / Partial Credit model, the model is refitted, and Q_3 residuals are computed under estimator. The distribution of max(Q3) - mean(Q3) across iterations provides empirical critical values. Failed iterations (e.g., due to convergence issues) are silently discarded.

Supports both dichotomous data (simulated via psychotools::rrm()) and polytomous data (via an internal partial credit score simulator).

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 mean, max, diff (one row per successful iteration).

pair_results

Long data.frame with columns Item1, Item2, Q3, iteration — one row per item pair per successful iteration. Used by RMlocdepQ3Plot.

pair_cutoffs

data.frame with per-pair cutoff summaries: Item1, Item2, Q3_low, Q3_high. Boundaries are computed via the method specified by cutoff_method.

actual_iterations

Number of successful iterations.

sample_n

Number of persons in the original data.

sample_n_total

Equal to sample_n: no respondents are dropped (incomplete responses are retained). Stored for consistency with the other ⁠*Cutoff()⁠ objects.

sample_has_na

Logical. Whether the data contained any missing values.

sample_summary

Summary statistics of estimated person parameters.

item_names

Character vector of item names from data.

max_diff, sd_diff

Max and SD of the diff distribution.

p95, p99, p995, p999

Empirical percentiles of diff.

suggested_cutoff

The 99th percentile (p99) — recommended scalar cutoff for RMlocdepQ3.

cutoff_method

The method used for pair_cutoffs ("hdci" or "quantile").

hdci_width

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

estimator

The estimator used for the simulated Q_3 ("CML" or "MML"); reused by RMlocdepQ3 and RMlocdepQ3Plot.

dgp

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

See Also

RMlocdepQ3

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)

  # Few iterations for a fast example; use 500+ in real analyses
  cutoff_res <- RMlocdepQ3Cutoff(sim_data, iterations = 50, parallel = FALSE,
                                 seed = 42)
  cutoff_res$suggested_cutoff  # 99th percentile

  # Use the cutoff in RMlocdepQ3()
  RMlocdepQ3(sim_data, cutoff = cutoff_res$suggested_cutoff)
}

Plot Distribution of Simulated Q_3 Residual Correlations

Description

Visualises the distribution of simulation-based Yen's Q_3 residual correlations per item pair from RMlocdepQ3Cutoff, optionally overlaying observed Q_3 values computed from real data via mirt::residuals(..., type = "Q3").

Usage

RMlocdepQ3Plot(simfit, data, items = NULL, n_pairs = NULL)

Arguments

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).

The ⁠$pairs⁠ plot shows one row per item pair (labelled as "Item1 - Item2"). Only the upper triangle of the Q_3 matrix is plotted (pairs are unordered under symmetric Q_3, unlike partial gamma which is direction-dependent).

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

When data is supplied, the function:

1. Computes observed Q_3 residual correlations under the same estimator used to build simfit (its ⁠$estimator⁠: CML/WLE by default, or MML via mirt).

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

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

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

Value

A named list of two ggplot objects (mirroring the ⁠$matrix⁠ / ⁠$pairs⁠ structure of RMlocdepQ3's table output):

⁠$pairs⁠

the per-pair plot described below (always returned).

⁠$matrix⁠

a lower-triangle tile heatmap of the observed Q_3 matrix, with pairs above the global dynamic cut-off outlined. This needs the observed data, so it is NULL (with a message) when data is not supplied. When items is given, the heatmap is subset to those items; n_pairs does not apply to it.

See Also

RMlocdepQ3, RMlocdepQ3Cutoff, RMlocdepGammaPlot

Examples

if (requireNamespace("ggplot2", 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 simulation (use more iterations, e.g. 500+, in real analyses)
  cutoff_res <- RMlocdepQ3Cutoff(sim_data, iterations = 50,
                                 parallel = FALSE, seed = 42)

  # Simulated distribution only
  RMlocdepQ3Plot(cutoff_res)

  # With observed Q3 overlaid
  RMlocdepQ3Plot(cutoff_res, data = sim_data)

  # Top 10 pairs by departure from null
  RMlocdepQ3Plot(cutoff_res, data = sim_data, n_pairs = 10)
}

Person-Fit Statistics for a Rasch / Partial Credit Model

Description

Computes per-respondent person-fit statistics and resampling-based p-values. Three statistics are reported: conditional infit and outfit mean-squares (MSQ) – the magnitude (effect-size) measures familiar from Rasch analysis – and the standardized log-likelihood statistic lz. The conditional MSQ statistics use response probabilities conditional on the total score, which require no person estimate and are therefore unbiased (Kreiner & Christensen, 2011); they are computed on each person's observed response pattern, so partial missingness is handled directly.

Usage

RMpersonFit(
  data,
  statistics = c("infit", "outfit", "lz"),
  estimator = c("CML", "MML"),
  theta_method = c("WLE", "EAP"),
  iterations = 500L,
  flag_alpha = 0.05,
  flag = c("both", "underfit"),
  zstd = FALSE,
  parallel = FALSE,
  n_cores = NULL,
  seed = NULL,
  output = c("kable", "dataframe", "ggplot")
)

Arguments

Details

Statistical significance is assessed by Monte-Carlo resampling under the fitted model rather than by an assumed asymptotic distribution. The asymptotic null distributions of MSQ and lz are known to be unreliable – the Wilson-Hilferty (ZSTD) transformation of MSQ in particular adds nothing once conditional estimation is used (Müller, 2020) – so resampling provides the valid reference (Sinharay, 2016).

Conditional MSQ. For person v with answered items A_v and total score r_v, the conditional residual is z_{vi} = (x_{vi} - E(X_{vi} \mid R_v = r_v)) / \sqrt{\mathrm{Var}(X_{vi} \mid R_v = r_v)}, with conditional moments obtained from elementary symmetric functions of the answered items' parameters. Outfit is the unweighted mean of z_{vi}^2; infit is the information-weighted mean. Because the conditional moments use only item parameters, no (biased) person estimate enters the residual.

Interpreting MSQ direction. Values near 1 indicate fit. MSQ > 1 (underfit) means the responses are noisier / more erratic than the model expects – the validity-relevant direction, indicating careless or aberrant responding. MSQ < 1 (overfit) means responses are overly deterministic (Guttman-like); usually benign for score validity, though a suspiciously perfect pattern can occasionally signal copied or fabricated data. Use flag = "underfit" to flag only the former. lz is one-sided: low (negative) values flag the aberrant (underfit-like) direction.

lz. The standardized log-likelihood of the response pattern evaluated at the estimated person location (Drasgow, Levine & Williams, 1985); small (negative) values indicate misfit. Because the location is estimated rather than known, the variance of lz is below 1 and its asymptotic standard-normal null is invalid, which makes a naive test conservative (Snijders, 2001; Sinharay, 2016). Instead of applying the analytic lz\* standardization – derived by Snijders (2001) for dichotomous items and extended to polytomous / mixed-format items by Sinharay (2016, 2026) – easyRasch2 obtains the reference distribution by resampling with the person location re-estimated for every simulated pattern (see Resampling). This reproduces the ability-estimation effect and is the resampling analogue of lz\*, so the resampled p-value is well calibrated.

Resampling. Each statistic is referenced against the null that matches it, removing the need to choose a scheme. The conditional MSQ statistics use patterns sampled conditional on the person's total score – the exact Rasch-native null, requiring no person estimate and fully consistent with the (conditional) statistic. lz, which is defined at the estimated location, uses patterns simulated at that location and then scored with the location re-estimated from each simulated pattern, so the null carries the same ability-estimation effect as the observed lz (for a Rasch / PCM model the location is a function of the total score, so this re-estimation is a cheap score-based lookup). Both schemes are computed per person over the items actually answered, so partial missingness is handled by either. The p-value is the proportion of the iterations replicates at least as extreme as the observed value. This follows the resampling-based person-fit approach (Sinharay, 2016) and the bootstrap recommendation of Müller (2020).

Value

For output = "dataframe", a data.frame with one row per respondent (input order): id, n_answered, sum_score, the requested statistics (infit_msq, outfit_msq, lz), their resampled p-values (p_infit, p_outfit, p_lz) when iterations > 0, flagged, and – if zstd = TRUE – infit_zstd, outfit_zstd. The flagged column is TRUE when the marginal (uncorrected) p-value of any requested statistic is below flag_alpha for that respondent; it is therefore a per-person screening flag, not corrected for the number of respondents tested (so under fit expect about flag_alpha of respondents to be flagged by chance). Each output = "ggplot" map instead colours and counts respondents flagged by its own statistic alone. Extreme scorers (minimum or maximum possible given the items they answered) receive NA statistics and are not assessed. For output = "kable" the same content as a knitr_kable; for output = "ggplot" the named list of maps described under that argument.

References

de la Torre, J., & Deng, W. (2008). Improving person-fit assessment by correcting the ability estimate and its reference distribution. Journal of Educational Measurement, 45(2), 159-177. doi:10.1111/j.1745-3984.2008.00058.x

Drasgow, F., Levine, M. V., & Williams, E. A. (1985). Appropriateness measurement with polychotomous item response models and standardized indices. British Journal of Mathematical and Statistical Psychology, 38(1), 67-86. doi:10.1111/j.2044-8317.1985.tb00817.x

Kreiner, S., & Christensen, K. B. (2011). Exact evaluation of bias in Rasch model residuals. Advances in Mathematics Research, 12, 19-40.

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

Sinharay, S. (2016). Assessment of person fit using resampling-based approaches. Journal of Educational Measurement, 53(1), 63-85. doi:10.1111/jedm.12101

Sinharay, S. (2016). Asymptotically correct standardization of person-fit statistics beyond dichotomous items. Psychometrika, 81(4), 992-1013. doi:10.1007/s11336-015-9465-x

Sinharay, S. (2026). Refining the asymptotically correct standardization of person-fit statistics for mixed-format tests. British Journal of Mathematical and Statistical Psychology. doi:10.1111/bmsp.70049

Snijders, T. A. B. (2001). Asymptotic null distribution of person fit statistics with estimated person parameter. Psychometrika, 66(3), 331-342. doi:10.1007/BF02294440

See Also

RMpersonParameters(), RMitemInfit(), RMitemParameters()

Examples

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

# Conditional infit/outfit MSQ + lz with resampled p-values
RMpersonFit(dat, iterations = 200, output = "dataframe") |> head()

# Person-fit maps: a named list with one plot per statistic
if (requireNamespace("ggplot2", quietly = TRUE)) {
  plots <- RMpersonFit(dat, iterations = 200, output = "ggplot")
  plots$infit    # infit map (each plot's caption reports the % flagged)
}

# Flag only underfit (noisy responding), the validity-relevant direction
RMpersonFit(dat, iterations = 200, flag = "underfit",
            output = "dataframe") |> head()

Person Locations for a Rasch / Partial Credit Model

Description

Estimates a person location (theta) and its standard error of measurement (SEM) for every respondent. Estimation is performed directly on each person's observed response pattern, so partial missingness is handled correctly: two respondents with the same sum score on different subsets of items receive different estimates. This avoids the sum-score lookup used by RMscoreSE(), which assumes complete data.

Usage

RMpersonParameters(
  data,
  method = c("WLE", "EAP"),
  item_params = NULL,
  estimator = c("CML", "MML"),
  theta_range = c(-10, 10),
  prior_mean = 0,
  prior_sd = NULL,
  output = c("kable", "dataframe", "ggplot", "file"),
  filename = NULL
)

Arguments

Details

Item parameters are obtained once (by CML or MML, or taken from item_params) and treated as fixed. Each person's theta is then estimated from the items they answered.

WLE solves Warm's weighted-likelihood score equation ⁠l'(theta) + J(theta) / (2 I(theta)) = 0⁠ by bracketed root finding, where the bias-correction term ⁠J / (2 I)⁠ keeps the equation solvable at the boundaries. The SEM is 1 / sqrt(I(theta)). Unlike the MLE, Warm's estimator therefore yields finite locations for the minimum and maximum possible scores (a sensible step beyond the next-most extreme score rather than +/-Inf; cf. Warm, 1989), matching established implementations such as catR and TAM. Such scores are still marked in the extreme column because they are extrapolated and carry large standard errors. Only if the root lies outside theta_range is the boundary returned with NA SEM.

EAP integrates the pattern likelihood against a normal prior over a quadrature grid; the estimate is the posterior mean and the SEM is the posterior standard deviation. Unlike WLE, EAP is finite at extreme scores because the prior shrinks them inward, at the cost of depending on the assumed prior. The prior actually used – including the marginal-ML estimate of prior_sd when it is left NULL – is reported in the kable caption and attached to the result as attr(result, "prior").

Value

For output = "dataframe", a data.frame with one row per respondent (in input order; respondents with no responses are dropped with a message) and columns theta (the person location), sem (standard error of measurement), sum_score, n_answered (number of non-missing responses), and extreme (logical: a minimum or maximum possible score given the items answered). For output = "kable", the same content as a knitr_kable object. For output = "ggplot", a histogram of theta.

References

Warm, T. A. (1989). Weighted likelihood estimation of ability in item response theory. Psychometrika, 54(3), 427-450. doi:10.1007/BF02294627

Bock, R. D., & Mislevy, R. J. (1982). Adaptive EAP estimation of ability in a microcomputer environment. Applied Psychological Measurement, 6(4), 431-444. doi:10.1177/014662168200600405

Magis, D. (2015). A note on weighted likelihood and Jeffreys modal estimation of proficiency levels in polytomous item response models. Psychometrika, 80(1), 200-204. doi:10.1007/s11336-013-9378-5

See Also

RMitemParameters(), RMscoreSE(), RMreliability()

Examples

set.seed(1)
dat <- as.data.frame(
  matrix(sample(0:2, 200 * 6, replace = TRUE), nrow = 200, ncol = 6)
)
colnames(dat) <- paste0("Item", 1:6)
# Introduce some missingness
dat[cbind(sample(200, 30), sample(6, 30, replace = TRUE))] <- NA

# Default: WLE person locations
RMpersonParameters(dat, output = "dataframe") |> head()

# EAP with a data-estimated prior SD
eap <- RMpersonParameters(dat, method = "EAP", output = "dataframe")
attr(eap, "prior")

# EAP with a fixed N(0, 1) prior
RMpersonParameters(dat, method = "EAP", prior_sd = 1, output = "dataframe") |>
  head()

# Write the person-location table to a CSV (also returned invisibly)
RMpersonParameters(dat, output = "file",
                   filename = tempfile(fileext = ".csv"))

Item Response Distribution Bar Chart

Description

Creates a faceted bar chart showing the response distribution for each item, with counts and percentages displayed on each bar. Each item gets its own panel, with response categories on the x-axis and percentage of responses on the y-axis. This is a descriptive data visualization tool intended for use before model fitting.

Usage

RMplotBar(
  data,
  item_labels = NULL,
  category_labels = NULL,
  ncol = 1L,
  label_wrap = 25L,
  text_y = 6,
  viridis_option = "A",
  viridis_end = 0.9,
  font = "sans"
)

Arguments

Details

Each item is displayed as a separate facet panel with the item label in the strip on the left side. Bars are coloured by response category using the viridis palette. Each bar shows the count (n = X) as text.

The plot caption reports the sample in the standard ⁠n = X respondents (policy)⁠ form; item-level NAs are retained – each bar counts the non-missing responses for that item.

Input requirements:

• All columns must be numeric (integer-valued).

• The data frame must contain at least 2 columns (items) and at least 1 row (person).

Value

A ggplot2::ggplot object.

See Also

RMplotStackedbar(), RMplotTile()

Examples

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

  # Basic response distribution plot
  RMplotBar(pcmdat2)

  # With custom item labels
  RMplotBar(
    pcmdat2,
    item_labels = c("Mood", "Sleep", "Appetite", "Energy")
  )

  # Two-column layout with wrapped labels
  RMplotBar(
    pcmdat2,
    item_labels = c(
      "General mood and emotional wellbeing",
      "Quality of sleep at night",
      "Appetite and eating habits",
      "Overall energy level during the day"
    ),
    ncol = 2, label_wrap = 20
  )

  # With custom category labels
  RMplotBar(
    pcmdat2,
    category_labels = c("Never", "Sometimes", "Often")
  )
}

Stacked Bar Chart of Item Response Distributions

Description

Creates a horizontal stacked bar chart showing the response distribution for all items. Each bar represents one item, with segments coloured by response category. Counts are displayed as text labels within each segment. This is a descriptive data visualization tool intended for use before model fitting.

Usage

RMplotStackedbar(
  data,
  item_labels = NULL,
  category_labels = NULL,
  show_n = TRUE,
  show_percent = FALSE,
  text_color = "sienna1",
  text_size = 3,
  min_label_n = 0L,
  viridis_option = "D",
  viridis_end = 0.99,
  title = "Item responses"
)

Arguments

Details

Items are displayed on the y-axis in the same order as the columns in data (first column at the top). Each bar is divided into segments representing response categories, with the lowest category on the left and the highest on the right. The total bar length equals the number of non-missing responses for that item.

Categories with zero responses still appear in the legend but produce no visible bar segment, which helps identify gaps in the response distribution.

The plot caption reports the sample in the standard ⁠n = X respondents (policy)⁠ form; item-level NAs are retained – each bar counts the non-missing responses for that item.

Input requirements:

• All columns must be numeric (integer-valued).

• The data frame must contain at least 2 columns (items) and at least 1 row (person).

Value

A ggplot2::ggplot object.

See Also

RMplotBar(), RMplotTile()

Examples

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

  # Basic stacked bar chart
  RMplotStackedbar(pcmdat2)

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

  # Show percentages, suppress small segments
  RMplotStackedbar(
    pcmdat2,
    show_percent = TRUE,
    show_n       = FALSE,
    min_label_n  = 5
  )
}

Tile Plot of Item Response Distributions

Description

Creates a tile (heat map) plot showing the distribution of responses across all items and response categories. Each cell displays the count (or percentage) of responses, with optional conditional highlighting for cells with low counts. Optional faceting by a grouping variable is provided for inspecting subgroup response distributions before DIF analyses – particularly useful for spotting empty categories or under-represented subgroups before fitting Rasch models per group.

Usage

RMplotTile(
  data,
  group = NULL,
  cutoff = 10,
  highlight = TRUE,
  percent = FALSE,
  text_color = "orange",
  item_labels = NULL,
  category_labels = NULL,
  group_labels = NULL,
  facet_ncol = NULL,
  output = c("ggplot", "dataframe")
)

Arguments

Details

Adapted from easyRaschBayes::plot_tile() and extended with the group parameter for faceted display.

Items are placed on the y-axis (in the same order as the columns of data, top to bottom) and response categories on the x-axis. Cell shading represents the count of responses (darker = more responses). Categories with zero responses are explicitly shown (n = 0), which helps identify gaps in the response distribution – one of the primary purposes of the plot, especially before DIF analyses where under-represented categories within a subgroup can break model fitting on that subgroup.

When group is supplied, percentages and the highlight cutoff are applied within each group, so a cell labelled "5" in the group-A facet contains the count for group A only.

The plot caption reports the sample in the standard ⁠n = X of Y respondents (policy)⁠ form: rows with NA group are dropped (and counted in Y only), while item-level NAs are retained – each cell simply counts the non-missing responses for that item.

Value

Either a ggplot object or a data.frame, depending on output.

Examples

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

  # Basic tile plot
  RMplotTile(pcmdat2)

  # With percentages
  RMplotTile(pcmdat2, percent = TRUE)

  # Faceted by an external grouping variable
  set.seed(1)
  grp <- sample(c("A", "B"), nrow(pcmdat2), replace = TRUE)
  RMplotTile(pcmdat2, group = grp)

  # With custom labels and tighter cutoff
  RMplotTile(pcmdat2,
             group = grp,
             group_labels = c("Female", "Male"),
             cutoff = 5,
             facet_ncol = 2)

  # Underlying counts as a data.frame
  RMplotTile(pcmdat2, group = grp, output = "dataframe")
}

Reliability metrics for a Rasch model

Description

Computes three reliability indices for a Rasch / partial credit model: the Person Separation Index (PSI) – WLE-based separation reliability – the marginal reliability (native, CML test information integrated over the estimated normal latent density), and Relative Measurement Uncertainty (RMU) via RMUreliability() applied to plausible values from mirt::fscores().

Usage

RMreliability(
  data,
  conf_int = 0.95,
  draws = 1000,
  rmu_iter = 50,
  estim = "WLE",
  boot = FALSE,
  boot_iter = 200,
  parallel = TRUE,
  n_cores = NULL,
  seed = NULL,
  verbose = FALSE,
  theta_range = c(-10, 10),
  output = "kable"
)

Arguments

Details

Confidence intervals for PSI and Marginal reliability are obtained by non-parametric bootstrap (resampling respondents; all three indices are recomputed natively per resample, no model is refitted by mirt). The RMU interval is the HDCI of correlations across plausible-value draws, averaged over rmu_iter random splits of the draws.

Marginal reliability is the native Green (1984) coefficient, 1 - \overline{1/I(\theta)}/\sigma^2, where the test information I(\theta) is summed from the CML item parameters and the average error variance is taken over the estimated normal latent density N(0, \sigma^2) (\sigma from marginal ML). Unlike mirt::marginal_rxx(), which assumes N(0,1), this integrates over the estimated latent variance, so it is correct on the Rasch logit scale (where \sigma is typically well above 1, and the N(0,1) assumption underestimates reliability). It is the model-based complement to the sample-based PSI; a large gap between the two flags an off-target or non-normal sample.

PSI is the WLE-based separation reliability, 1 - \overline{SEM^2} / \mathrm{Var}(\hat\theta), computed from CML item thresholds (psychotools) and Warm's WLE person locations / analytic SEMs. Respondents with extreme (min/max) raw scores are excluded – their boundary estimates would inflate the person variance and overstate reliability. (Earlier versions used eRm::SepRel() with MLE; the values can differ, most noticeably for scales with many extreme scorers, e.g. dichotomous items.)

RMU is from Bignardi, Kievit, & Bürkner (2025), modified here to use mirt plausible values rather than fully Bayesian posterior draws (see Mislevy, 1991, for the plausible-values framework).

Bootstrap iterations that fail to converge are silently dropped.

Value

• If output = "kable": a knitr_kable object with one row per metric.

• If output = "dataframe": a data.frame with columns metric, estimate, lower, upper, notes.

References

Bignardi, G., Kievit, R., & Bürkner, P. C. (2025). A general method for estimating reliability using Bayesian Measurement Uncertainty. PsyArXiv. doi:10.31234/osf.io/h54k8_v1

Green, B. F., Bock, R. D., Humphreys, L. G., Linn, R. L., & Reckase, M. D. (1984). Technical Guidelines for Assessing Computerized Adaptive Tests. Journal of Educational Measurement, 21(4), 347–360. doi:10.1111/j.1745-3984.1984.tb01039.x

Mislevy, R. J. (1991). Randomization-Based Inference about Latent Variables from Complex Samples. Psychometrika, 56(2), 177-196. doi:10.1007/BF02294457

Adams, R. J. (2005). Reliability as a measurement design effect. Studies in Educational Evaluation, 31(2), 162-172. doi:10.1016/j.stueduc.2005.05.008

See Also

RMUreliability()

Examples

if (requireNamespace("ggdist", quietly = TRUE) &&
    requireNamespace("eRm", quietly = TRUE)) {
  set.seed(1)
  RMreliability(eRm::raschdat1[, 1:20], draws = 1000)

  # Bootstrap CI for PSI and Marginal
  # (use more bootstrap iterations, e.g. 200+, in real analyses)
  RMreliability(eRm::raschdat1[, 1:20], draws = 1000,
                boot = TRUE, boot_iter = 25, parallel = FALSE, seed = 42)
}

Raw-Score to Logit Score Transformation Table

Description

For a given set of items, returns the score-to-theta lookup that maps each possible raw sum score to a person-location estimate (in logits) and its standard error. Useful when reporting a scale's measurement properties or converting raw totals to interval-scaled scores for downstream analysis.

Usage

RMscoreSE(
  data,
  method = "WLE",
  output = "kable",
  ci_multiplier = 1.96,
  point_size = 3,
  error_width = 0.5,
  theta_range = c(-10, 10)
)

Arguments

Details

The function automatically detects whether the data is dichotomous (max score 1) or polytomous (max score > 1) and selects the appropriate Rasch / Partial Credit model.

method = "WLE" fits the model by CML with psychotools::pcmodel(), centres the item thresholds to grand-mean-zero, and solves Warm's weighted-likelihood equation for each raw score with the same engine used by RMpersonParameters(); the two functions therefore report identical locations and standard errors. Warm's bias correction yields finite locations even at the minimum and maximum scores (only a root outside theta_range is clamped to the boundary with NA SE). The reported logit_se is the information-based standard error 1 / sqrt(I(theta)), matching catR, TAM and most Rasch software.

method = "EAP" fits the model with mirt::mirt(..., itemtype = "Rasch") (MML) and obtains sum-score-based EAP estimates and posterior SDs via mirt::fscores(method = "EAPsum", full.scores = FALSE, full.scores.SE = TRUE). EAP estimates are finite at all score boundaries (the prior shrinks them inward), but they depend on the assumed normal prior on theta. Item parameters from MML differ slightly from the CML values used by the WLE path; for well-behaved data the difference is small.

Value

• If output = "kable": a knitr_kable object with columns "Ordinal sum score", "Logit score", and "Logit std.error", and a caption noting the estimation method.

• If output = "dataframe": a data.frame with columns raw_score, logit_score, and logit_se (one row per possible raw sum score from 0 to the theoretical maximum).

• If output = "ggplot": a ggplot object — points at each (logit_score, raw_score) with horizontal error bars at ± ⁠ci_multiplier × logit_se⁠.

References

Warm, T. A. (1989). Weighted likelihood estimation of ability in item response theory. Psychometrika, 54(3), 427-450. doi:10.1007/BF02294627

Bock, R. D., & Mislevy, R. J. (1982). Adaptive EAP estimation of ability in a microcomputer environment. Applied Psychological Measurement, 6(4), 431-444. doi:10.1177/014662168200600405

Examples

set.seed(42)
sim_data <- as.data.frame(
  matrix(sample(0:3, 200 * 6, replace = TRUE), nrow = 200, ncol = 6)
)
colnames(sim_data) <- paste0("Item", 1:6)

# Default kable output, WLE
RMscoreSE(sim_data)

# Underlying data.frame
RMscoreSE(sim_data, output = "dataframe")

# ggplot figure
if (requireNamespace("ggplot2", quietly = TRUE)) {
  RMscoreSE(sim_data, output = "ggplot")
}

# EAP via mirt
RMscoreSE(sim_data, method = "EAP")

Person-Item Targeting Plot (Wright Map)

Description

Produces a three-panel targeting plot with a shared logit scale x-axis:

1. Top: Histogram of person location estimates, with a reference line for the mean (or median) and shading for ±1 SD (or ±1 MAD).

2. Middle: Inverted histogram of item threshold locations, with the same summary annotations.

3. Bottom: Dot-and-whisker plot of individual item thresholds with confidence intervals based on threshold standard errors.

Usage

RMtargeting(
  data,
  robust = FALSE,
  sort_items = c("data", "location"),
  bins,
  xlim = c(-4, 4),
  ci_level = 0.95,
  person_fill = "#0072B2",
  threshold_fill = "#D55E00",
  height_ratios = c(3, 2, 5),
  output = "patchwork"
)

Arguments

Details

Together, the top and middle panels form a back-to-back histogram that makes it easy to assess whether the test is well-targeted to the sample.

Estimation method selection. The function checks whether any item response category has fewer than 3 observations. If all categories have at least 3 responses, item threshold locations and their standard errors are estimated via Conditional Maximum Likelihood (CML) using psychotools::pcmodel() (a dichotomous item is a 2-category PCM). If any category has fewer than 3 responses, the function falls back to Marginal Maximum Likelihood (MML) estimation via mirt::mirt() with itemtype = "Rasch" and SE = TRUE, which is more numerically stable under sparse-category conditions. A message is emitted when the MML fallback is used.

In both cases, item threshold locations are centered (shifted so the grand mean of all thresholds equals zero).

Person estimates are obtained by Warm's weighted likelihood (WLE) from the fitted item thresholds, consistent with the rest of the package. WLE is finite at extreme scores, so all-zero and perfect responders are located rather than dropped.

Confidence intervals for item thresholds are based on Wald-type intervals: threshold estimate ± z × SE, where z is the standard normal quantile corresponding to ci_level.

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

Value

• If output = "patchwork": a patchwork object (combined ggplot).

• If output = "list": a named list with elements p1 (person histogram), p2 (threshold histogram), and p3 (item threshold dot-whisker plot).

References

Wright, B. D. & Stone, M. H. (1979). Best Test Design. MESA Press.

See Also

psychotools::pcmodel(), mirt::mirt()

Examples

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

  # Default: mean/SD, data order, 95% CI
  RMtargeting(sim_data)

  # Robust (median/MAD), sorted by location, 84% CI
  RMtargeting(sim_data, robust = TRUE, sort_items = "location",
              ci_level = 0.84)

  # Get list of sub-plots for customisation
  plots <- RMtargeting(sim_data, output = "list")
  plots$p1 + ggplot2::ggtitle("My custom title")

  # Dichotomous example
  sim_bin <- as.data.frame(
    matrix(sample(0:1, 200 * 10, replace = TRUE), nrow = 200, ncol = 10)
  )
  colnames(sim_bin) <- paste0("Item", 1:10)
  RMtargeting(sim_bin)
}

Function renaming in easyRasch2 0.8.0

Description

In version 0.8.0, 22 exported functions were renamed under a consistent domain-prefix \to method \to variant-suffix scheme. The renames have no semantic effect — only names changed. Use the table below to migrate existing scripts.

Domain prefixes

• RMdif — differential item functioning

• RMlocdep — local dependence

• RMitem — per-item statistics (fit, restscore, ICC, hierarchy)

• RMdim — dimensionality assessment

• RMplot — descriptive response-distribution plots

Variant suffixes attached at the end:

• Cutoff — simulation-based critical values

• Plot — visualisation

• MI — multiple-imputation variant (replaces former ⁠_mi⁠)

Old to new name map

Unchanged: RMreliability(), RMUreliability(), RMitemHierarchy(), RMdifLR(), RMdifTree(), RMlocdepQ3(), RMtargeting(), RMscoreSE().

No deprecation aliases

No deprecation shims are shipped: each old name is simply gone. This keeps the export list clean ahead of the first CRAN submission. Run a search-and-replace against the table above to migrate.

knitr knit_print method for RMlocdepGamma kable output

Description

Inside a knitr / Quarto / R Markdown chunk, returns the pre-combined two-table asis string so pandoc renders them as two distinct pipe tables. Outside knitr, R's normal dispatch falls back to print.RMlocdepGamma().

Usage

## S3 method for class 'RMlocdepGamma'
knit_print(x, ...)

Arguments

Value

A knit_asis character object.

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]))

Print method for RMlocdepGamma kable output

Description

Prints the two rest-score direction tables stacked vertically with a blank line between them. Each table renders via knitr_kable's own print method as a clean pipe-markdown table.

Usage

## S3 method for class 'RMlocdepGamma'
print(x, ...)

Arguments

Value

Invisibly returns x.

Run item-restscore bootstrap iterations in parallel using mirai

Description

Run item-restscore bootstrap iterations in parallel using mirai

Usage

run_boot_restscore_parallel(
  iterations,
  boot_seeds,
  boot_data_list,
  n_cores,
  verbose = FALSE
)

Arguments

Value

List of raw results (one element per iteration).

Run item-restscore bootstrap iterations sequentially

Description

Run item-restscore bootstrap iterations sequentially

Usage

run_boot_restscore_sequential(
  iterations,
  boot_seeds,
  boot_data_list,
  verbose = FALSE
)

Arguments

Value

List of raw results (one element per iteration).

Run infit simulations in parallel using mirai

Description

Run infit simulations in parallel using mirai

Usage

run_infit_sim_parallel(
  iterations,
  sim_seeds,
  sim_data_list,
  n_cores,
  verbose = FALSE
)

Arguments

Value

List of raw results (one element per iteration).

Run infit simulations sequentially

Description

Run infit simulations sequentially

Usage

run_infit_sim_sequential(iterations, sim_seeds, sim_data_list, verbose = FALSE)

Arguments

Value

List of raw results (one element per iteration).

Run partial gamma LD simulations in parallel using mirai

Description

Run partial gamma LD simulations in parallel using mirai

Usage

run_partgam_LD_sim_parallel(
  iterations,
  sim_seeds,
  sim_data_list,
  n_cores,
  verbose = FALSE
)

Arguments

Value

List of raw results (one element per iteration).

Run partial gamma LD simulations sequentially

Description

Run partial gamma LD simulations sequentially

Usage

run_partgam_LD_sim_sequential(
  iterations,
  sim_seeds,
  sim_data_list,
  verbose = FALSE
)

Arguments

Value

List of raw results (one element per iteration).

Run partial gamma DIF simulations in parallel using mirai

Description

Run partial gamma DIF simulations in parallel using mirai

Usage

run_partgam_sim_parallel(
  iterations,
  sim_seeds,
  sim_data_list,
  n_cores,
  verbose = FALSE
)

Arguments

Value

List of raw results (one element per iteration).

Run partial gamma DIF simulations sequentially

Description

Run partial gamma DIF simulations sequentially

Usage

run_partgam_sim_sequential(
  iterations,
  sim_seeds,
  sim_data_list,
  verbose = FALSE
)

Arguments

Value

List of raw results (one element per iteration).

Run Q3 simulations in parallel using mirai

Description

Run Q3 simulations in parallel using mirai

Usage

run_q3_sim_parallel(
  iterations,
  sim_seeds,
  sim_data_list,
  n_cores,
  verbose = FALSE
)

Arguments

Value

List of raw results (one element per iteration).

Run Q3 simulations sequentially

Description

Run Q3 simulations sequentially

Usage

run_q3_sim_sequential(iterations, sim_seeds, sim_data_list, verbose = FALSE)

Arguments

Value

List of raw results (one element per iteration).

Run a single item-restscore bootstrap iteration

Description

Run a single item-restscore bootstrap iteration

Usage

run_single_boot_restscore(seed, data_list)

Arguments

Value

A data.frame with columns Item, item_restscore, diff, diff_abs, or a character string on failure.

Run a single infit simulation iteration

Description

Run a single infit simulation iteration

Usage

run_single_infit_sim(seed, data_list)

Arguments

Value

A data.frame with columns Item, InfitMSQ, OutfitMSQ, or a character string on failure.

Run a single partial gamma LD simulation iteration

Description

Run a single partial gamma LD simulation iteration

Usage

run_single_partgam_LD_sim(seed, data_list)

Arguments

Value

A data.frame with columns Item1, Item2, and gamma, or a character string on failure.

Run a single partial gamma DIF simulation iteration

Description

Run a single partial gamma DIF simulation iteration

Usage

run_single_partgam_sim(seed, data_list)

Arguments

Value

A data.frame with columns Item and gamma, or a character string on failure.

Run a single Q3 simulation iteration

Description

Run a single Q3 simulation iteration

Usage

run_single_q3_sim(seed, data_list)

Arguments

Value

A list with mean and max Q3, or a character string on failure.