Package 'mfrmr'

Description

mfrmr provides estimation, diagnostics, and reporting utilities for many-facet Rasch-family models and bounded generalized partial-credit model workflows using a native R implementation.

Details

If you are new to the package, read the next four steps first and ignore the longer GPCM, simulation, and planning notes until the basic route works:

1. Fit with fit_mfrm() using method = "MML"

2. For RSM / PCM, run diagnose_mfrm() with diagnostic_mode = "both"

3. Read summary(fit) and summary(diag) before branching

4. Use plot_qc_dashboard() and reporting_checklist() as the first visual and reporting screens

Recommended workflow:

1. Fit model with fit_mfrm()

2. For RSM / PCM, compute diagnostics with diagnose_mfrm() and prefer diagnostic_mode = "both" when you want legacy residual continuity plus the newer strict marginal-fit screen

3. For RSM / PCM, run residual PCA with analyze_residual_pca() and parallel-analysis checks with check_residual_dimensionality() if needed

4. For RSM / PCM, estimate interactions with estimate_bias()

5. For RSM / PCM, choose a downstream branch: reporting_checklist() for manuscript/report preparation, or build_misfit_casebook() / build_linking_review() for operational misfit or anchor/drift review. After build_misfit_casebook(), inspect casebook$group_view_index before moving to source-specific plots.

6. For RSM / PCM, build narrative/report outputs with build_apa_outputs() and build_visual_summaries()

7. Treat GPCM, prediction, and planning helpers as advanced scope after the basic RSM / PCM route is working cleanly.

Guide pages:

• mfrmr_workflow_methods

• mfrmr_visual_diagnostics

• mfrmr_reports_and_tables

• mfrmr_reporting_and_apa

• mfrmr_linking_and_dff

• gpcm_capability_matrix

• mfrmr_compatibility_layer

Companion vignettes:

• vignette("mfrmr-workflow", package = "mfrmr")

• vignette("mfrmr-mml-and-marginal-fit", package = "mfrmr")

• vignette("mfrmr-visual-diagnostics", package = "mfrmr")

• vignette("mfrmr-reporting-and-apa", package = "mfrmr")

• vignette("mfrmr-linking-and-dff", package = "mfrmr")

A two-page landscape cheatsheet of the public API ships at system.file("cheatsheet", "mfrmr-cheatsheet.pdf", package = "mfrmr") (pre-rendered) and system.file("cheatsheet", "mfrmr-cheatsheet.Rmd", package = "mfrmr") (source). Open the PDF directly for a printable reference card, or knit the source with rmarkdown::render() when you want a customised version.

First 5-minute route

Use this order before exploring the broader feature surface:

1. fit_mfrm() with method = "MML"

2. diagnose_mfrm() with diagnostic_mode = "both" for RSM / PCM

3. summary(fit) and summary(diag)

4. plot_qc_dashboard() for first-pass triage

5. Choose the next branch: reporting_checklist() for reporting, build_weighting_audit() for Rasch-versus-GPCM weighting review, build_misfit_casebook() for operational case review, or build_linking_review() for operational linking review

Advanced scope

After the basic route above:

• the package now includes a first-version latent-regression MML branch for ordered-response RSM / PCM models with a one-dimensional conditional-normal population model and explicit one-row-per-person covariates expanded through stats::model.matrix()

• bounded GPCM support is summarized by gpcm_capability_matrix()

• bounded GPCM supports the core fit/summary/scoring/information path, direct Wright/pathway/CCC plots, residual-PCA follow-up, and the residual-based diagnostics tables/plots as exploratory tools

• posterior-predictive computation, MCMC engines, and Docker-based advanced runtimes are future extensions rather than requirements for the current bounded GPCM route

• direct GPCM data generation through build_mfrm_sim_spec(), extract_mfrm_sim_spec(), and simulate_mfrm_data() is available when the specification carries both thresholds and slopes

• fair-average, APA writer, package-native export/replay, and role-based planning/forecasting routes are available for bounded GPCM with explicit screening-tier caveats

• predict_mfrm_population() remains a scenario-level forecast helper and should not be described as the latent-regression estimator itself

• the role-based simulation/planning layer remains the PCM/GPCM route for two non-person facets, while build_mfrm_arbitrary_sim_spec(), extract_mfrm_arbitrary_sim_spec(), simulate_mfrm_arbitrary_data(), summarize_mfrm_sim_design(), plot_mfrm_sim_design(), summarize_mfrm_sim_grid(), plot_mfrm_sim_grid(), list_mfrm_sim_metrics(), plot_mfrm_sim_dashboard(), and evaluate_mfrm_bias_detection() provide a first RSM-based arbitrary-facet design, multi-metric dashboard, and bias-sensitivity branch

• latent-class mixture models and response-time / careless-rating adjustment are not estimated by mfrmr; use residual, person-fit, local-dependence, and rater-drift diagnostics as screening layers rather than as mixture-model substitutes

Equal weighting versus bounded GPCM

The package's operational reference route is still the Rasch-family RSM / PCM branch. That route enforces fixed discrimination and therefore preserves an equal-weighting scoring interpretation across observed ratings.

bounded GPCM is supported because some users want a slope-aware model- comparison or sensitivity layer inside the same many-facet workflow. However, the package does not treat bounded GPCM as a universal replacement for the Rasch-family route. A better fit under GPCM should be read as evidence about discrimination-based reweighting, not as an automatic reason to discard the equal-weighting model.

Observation weights are a different concept again. Optional Weight columns change how observed rating events enter estimation and summaries, but they do not create a free-form facet-weighting scheme and do not alter the fixed-discrimination meaning of RSM / PCM.

Function families:

• Model fitting: fit_mfrm(), summary.mfrm_fit(), plot.mfrm_fit()

• Legacy-compatible workflow wrapper: run_mfrm_facets(), mfrmRFacets()

• Diagnostics: diagnose_mfrm(), summary(diag), analyze_residual_pca(), check_residual_dimensionality(), plot_residual_pca(), plot_residual_dimensionality()

• Bias and interaction: estimate_bias(), estimate_all_bias(), summary(bias), bias_interaction_report(), plot_bias_interaction()

• Differential functioning: analyze_dff(), analyze_dif(), dif_interaction_table(), plot_dif_heatmap(), dif_report()

• Design simulation: build_mfrm_sim_spec(), extract_mfrm_sim_spec(), simulate_mfrm_data(), evaluate_mfrm_design(), evaluate_mfrm_signal_detection(), build_mfrm_arbitrary_sim_spec(), extract_mfrm_arbitrary_sim_spec(), simulate_mfrm_arbitrary_data(), summarize_mfrm_sim_design(), plot_mfrm_sim_design(), summarize_mfrm_sim_grid(), plot_mfrm_sim_grid(), list_mfrm_sim_metrics(), plot_mfrm_sim_dashboard(), evaluate_mfrm_bias_detection(), predict_mfrm_population(), predict_mfrm_units(), sample_mfrm_plausible_values() (including fit-derived empirical / resampled / skeleton-based simulation specifications; fixed-calibration unit scoring supports MML fits directly, latent-regression MML fits through the fitted population model when scored units also provide one-row-per-person background data, and JML fits through a post hoc reference-prior EAP layer; fit-derived simulation specifications, planning/forecasting helpers, curve reports, and graph-only exports are also available for bounded GPCM with the caveats documented in gpcm_capability_matrix())

• Reporting: build_apa_outputs(), build_visual_summaries(), reporting_checklist(), apa_table() for the full RSM / PCM route; bounded GPCM currently stays on the checklist / visual-summary / QC / direct-table / direct-plot side instead of the narrative layer

• Weighting review: compare_mfrm(), build_weighting_audit(), compute_information(), plot_information()

• Case review: build_misfit_casebook(), plot_unexpected(), plot_displacement(), plot_marginal_fit(), plot_marginal_pairwise()

• Linking and scale maintenance: audit_mfrm_anchors(), detect_anchor_drift(), build_equating_chain(), build_linking_review(), plot_anchor_drift()

• Dashboards and fit-direction rates: facet_quality_dashboard(), plot_facet_quality_dashboard(), fit_direction_summary(), plot_fit_direction_summary(), summarize_simulation_misfit(), plot_simulation_misfit_rates()

• Export / reproducibility: build_mfrm_manifest(), build_mfrm_replay_script(), build_conquest_overlap_bundle(), normalize_conquest_overlap_files(), normalize_conquest_overlap_tables(), audit_conquest_overlap(), export_mfrm_bundle() for the diagnostics-compatible Rasch-family route; bounded GPCM remains outside the current manifest/replay/bundle layer

• Equivalence: analyze_facet_equivalence(), plot_facet_equivalence()

• Data and anchors: describe_mfrm_data(), audit_mfrm_anchors(), make_anchor_table(), load_mfrmr_data()

Data interface:

• Input analysis data is long format (one row per observed rating).

• Required columns are one person column, one ordered score column, and one or more non-person facet columns named in facets = c(...).

• Score values should be ordered integer categories. Binary 0/1 or 1/2 input is supported as the two-category Rasch-family special case; by contrast, fractional score values should be recoded before fitting rather than relying on automatic coercion.

• If keep_original = FALSE, unused intermediate categories are collapsed to a contiguous internal scale and the mapping is stored in fit$prep$score_map.

• If the intended scale has unused boundary categories, such as a 1-5 scale with only 2-5 observed, set ⁠rating_min = 1, rating_max = 5⁠ so the zero-count boundary category remains in the fitted support. If unused intermediate categories should also remain in the original scale, set keep_original = TRUE.

• summary(describe_mfrm_data(...)) reports retained zero-count categories in Notes, printed Caveats, and ⁠$caveats⁠; summary(fit) carries full structured rows into printed Caveats and ⁠$caveats⁠, with ⁠Key warnings⁠ as a short triage subset. Summary-table exports route those rows through score_category_caveats or analysis_caveats. Treat adjacent thresholds as weakly identified when an intermediate category is unobserved.

• Optional columns such as Subset, Weight, and Group support linking, weighted analysis, and fairness-focused follow-up workflows.

• Packaged simulation data is available via load_mfrmr_data() or data().

Interpreting output

Core object classes are:

• mfrm_fit: fitted model parameters and metadata.

• mfrm_diagnostics: fit, facet-level reliability, and flag diagnostics, plus inter-rater agreement when one facet is treated as a rater facet.

• mfrm_bias: interaction bias estimates.

• mfrm_dff / mfrm_dif: differential-functioning contrasts and screening summaries.

• mfrm_population_prediction: scenario-level forecast summaries for one future design.

• mfrm_unit_prediction: posterior summaries for future or partially observed persons under the fitted scoring basis.

• mfrm_plausible_values: posterior draws for future or partially observed persons under the fitted scoring basis.

• mfrm_bundle families: summary/report bundles and plotting payloads.

Typical workflow

1. Prepare long-format data.

2. Fit with fit_mfrm().

3. For RSM / PCM, diagnose with diagnose_mfrm() and prefer diagnostic_mode = "both" for final MML runs.

4. For RSM / PCM, run analyze_dff() or estimate_bias() when fairness or interaction questions matter.

5. For RSM / PCM, report with build_apa_outputs() and build_visual_summaries().

6. For design planning, move to build_mfrm_sim_spec(), evaluate_mfrm_design(), and predict_mfrm_population(). bounded GPCM also supports direct simulation via extract_mfrm_sim_spec() / simulate_mfrm_data(), but not the broader planning helpers. Those helpers still assume two non-person facet roles even though the estimation core supports arbitrary facet counts. Treat evaluate_mfrm_design() as Monte Carlo design evaluation rather than a closed-form generalizability-theory D-study planner. predict_mfrm_population() remains the scenario-level forecast helper, not the latent-regression estimator.

7. For future-unit scoring, retain an MML calibration when you want the fitted marginal model directly, use an active latent-regression MML fit when scored units also provide one-row-per-person background data, or use a JML calibration when a post hoc fixed-calibration EAP layer is acceptable; then score with predict_mfrm_units() or sample_mfrm_plausible_values().

8. For bounded GPCM, use summary.mfrm_fit(), diagnose_mfrm(), analyze_residual_pca(), check_residual_dimensionality(), predict_mfrm_units(), sample_mfrm_plausible_values(), compute_information(), plot_qc_dashboard(), plot.mfrm_fit(), category_structure_report(), category_curves_report(), fair_average_table(), estimate_bias(), build_visual_summaries(), run_qc_pipeline(), graph-only facets_output_file_bundle(), direct simulation-spec generation/data generation, caveated APA/export/replay bundles, caveated role-based planning/forecasting, and the residual-based table helpers while FACETS compatibility score exports remain blocked. Use gpcm_capability_matrix() as the formal boundary statement.

Model formulation

The many-facet Rasch model (MFRM; Linacre, 1989) extends the basic Rasch model by incorporating multiple measurement facets into a single linear model on the log-odds scale.

General MFRM equation

For an observation where person $n$ with ability $\theta_n$ is rated by rater $j$ with severity $\delta_j$ on criterion $i$ with difficulty $\beta_i$, the probability of observing category $k$ (out of $K$ ordered categories) is:

$P(X_{nij} = k \mid \theta_n, \delta_j, \beta_i, \tau) = \frac{\exp\bigl[\sum_{s=1}^{k}(\theta_n - \delta_j - \beta_i - \tau_s)\bigr]} {\sum_{c=0}^{K}\exp\bigl[\sum_{s=1}^{c}(\theta_n - \delta_j - \beta_i - \tau_s)\bigr]}$

where $\tau_s$ are the Rasch-Andrich threshold (step) parameters and $\sum_{s=1}^{0}(\cdot) \equiv 0$ by convention. Additional facets enter as additive terms in the linear predictor $\eta = \theta_n - \delta_j - \beta_i - \ldots$.

This formulation generalises to any number of facets; the facets argument to fit_mfrm() accepts an arbitrary-length character vector.

Rating Scale Model (RSM)

Under the RSM (Andrich, 1978), all levels of the step facet share a single set of threshold parameters $\tau_1, \ldots, \tau_K$.

Partial Credit Model (PCM)

Under the PCM (Masters, 1982), each level of the designated step_facet has its own threshold vector on the package's common observed score scale. In the current implementation, threshold locations may vary by step-facet level, but the fitted score range is still defined by one global category set taken from the observed data.

Ordered-response scope

The current public response-model scope is ordered categorical only. Binary responses are the $K = 1$ special case of the same formulation, so they are handled through the ordinary ordered-score interface. This means mfrmr supports ordered binary and ordered polytomous data under RSM and PCM, plus a narrow bounded GPCM branch with one designated slope_facet that currently must equal step_facet. Unordered nominal/multinomial response models are not yet implemented.

Estimation methods

Marginal Maximum Likelihood (MML)

MML integrates over the person ability distribution using Gauss-Hermite quadrature, in the broader marginal-likelihood framework introduced by Bock & Aitkin (1981) for IRT:

$L = \prod_{n} \int P(\mathbf{X}_n \mid \theta, \boldsymbol{\delta}) \, \phi(\theta) \, d\theta \approx \prod_{n} \sum_{q=1}^{Q} w_q \, P(\mathbf{X}_n \mid \theta_q, \boldsymbol{\delta})$

where $\phi(\theta)$ is the assumed normal prior and $(\theta_q, w_q)$ are quadrature nodes and weights. Person estimates are obtained post-hoc via Expected A Posteriori (EAP):

$\hat{\theta}_n^{\mathrm{EAP}} = \frac{\sum_q \theta_q \, w_q \, L(\mathbf{X}_n \mid \theta_q)} {\sum_q w_q \, L(\mathbf{X}_n \mid \theta_q)}$

MML avoids the incidental-parameter problem and is generally preferred for smaller samples.

Note: Bock & Aitkin (1981) is the canonical citation for the Gauss-Hermite-quadrature MML framework. The default mfrmr engine (mml_engine = "direct") optimises this marginal log-likelihood by direct gradient methods (BFGS / L-BFGS-B), not by Bock & Aitkin's signature EM algorithm. The "em" and "hybrid" engines do follow the EM template but use a BFGS M-step rather than B&A's probit IRLS, because the target is the polytomous Rasch family rather than B&A's 2PL probit model.

Joint Maximum Likelihood (JML)

JML estimates all person and facet parameters simultaneously as fixed effects by maximising the joint log-likelihood $\ell(\boldsymbol{\theta}, \boldsymbol{\delta} \mid \mathbf{X})$ directly. It does not assume a parametric person distribution, which can be advantageous when the population shape is strongly non-normal, but parameter estimates are known to be biased when the number of persons is small relative to the number of items (Neyman & Scott, 1948). The package still accepts "JMLE" as a backward-compatible alias, but user-facing summaries and documentation use "JML" as the public label.

See fit_mfrm() for practical guidance on choosing between the two.

Strict marginal diagnostics

For RSM / PCM, diagnose_mfrm(..., diagnostic_mode = "both") returns two complementary targets: the legacy residual / EAP diagnostics and a marginal_fit layer whose expected counts and pairwise summaries are integrated over the posterior quadrature bundle rather than plugged in at the EAP point. The screen is structured as limited-information evidence (Orlando & Thissen, 2000; Haberman & Sinharay, 2013; Sinharay & Monroe, 2025), not as an omnibus accept / reject test, and it complements rather than replaces separation / reliability and inter-rater agreement summaries. The full derivation, with notation and pairwise local-dependence events, lives in vignette("mfrmr-mml-and-marginal-fit", package = "mfrmr").

Statistical background

Key statistics reported throughout the package:

Infit (Information-Weighted Mean Square)

Weighted average of squared standardized residuals, where weights are the model-based variance of each observation:

$\mathrm{Infit}_j = \frac{\sum_i Z_{ij}^2 \, \mathrm{Var}_i \, w_i} {\sum_i \mathrm{Var}_i \, w_i}$

Expected value is 1.0 under model fit. Values below 0.5 suggest overfit (Mead-style responses); values above 1.5 suggest underfit (noise or misfit). Infit is most sensitive to unexpected patterns among on-target observations (Wright & Masters, 1982).

Note: The 0.5–1.5 range is the general "productive for measurement" band given by Linacre (2002, RMT 16(2), 878). Context-specific bands come from Wright & Linacre (1994, RMT 8(3), 370): 0.8–1.2 for high-stakes MCQ, 0.7–1.3 for run-of-the-mill MCQ, 0.6–1.4 for rating-scale surveys, 0.5–1.7 for clinical observation, and 0.4–1.2 for judged performance. See also Bond & Fox (2015) for textbook summaries of these conventions.

Outfit (Unweighted Mean Square)

Simple average of squared standardized residuals:

$\mathrm{Outfit}_j = \frac{\sum_i Z_{ij}^2 \, w_i}{\sum_i w_i}$

Same expected value and flagging thresholds as Infit, but more sensitive to extreme off-target outliers (e.g., a high-ability person scoring the lowest category).

ZSTD (Standardized Fit Statistic)

Wilson-Hilferty (1931) cube-root transformation that converts the mean-square chi-square ratio to an approximate standard normal deviate:

$\mathrm{ZSTD} = \frac{\mathrm{MnSq}^{1/3} - (1 - 2/(9\,\mathit{df}))} {\sqrt{2/(9\,\mathit{df})}}$

Values near 0 indicate expected fit; $|\mathrm{ZSTD}| > 2$ flags potential misfit at the 5\ 1\ ZSTD is reported alongside every Infit and Outfit value.

PTMEA (Point-Measure Correlation)

Pearson correlation between observed scores and estimated person measures within each facet level. Positive values indicate that scoring aligns with the latent trait dimension; negative values suggest reversed orientation or scoring errors.

Separation

Package-reported separation is the ratio of adjusted true standard deviation to root-mean-square measurement error:

$G = \frac{\mathrm{SD}_{\mathrm{adj}}}{\mathrm{RMSE}}$

where $\mathrm{SD}_{\mathrm{adj}} = \sqrt{\mathrm{ObservedVariance} - \mathrm{ErrorVariance}}$. Higher values indicate the facet discriminates more statistically distinct levels along the measured variable. In mfrmr, Separation is the model-based value and RealSeparation provides a more conservative companion based on RealSE.

Reliability

$R = \frac{G^2}{1 + G^2}$

Analogous to Cronbach's alpha or KR-20 for the reproducibility of element ordering. In mfrmr, Reliability is the model-based value and RealReliability gives the conservative companion based on RealSE. For MML, these are anchored to observed-information ModelSE estimates for non-person facets; JML keeps them as exploratory summaries. This is a Rasch/FACETS-style separation reliability on the fitted logit scale, not an intra-class correlation. Use compute_facet_icc() only when you want the complementary random-effects variance-share view on the observed-score scale; for non-person facets, large ICC values indicate systematic facet variance rather than desirable measurement reliability.

Strata

Number of statistically distinguishable groups of elements:

$H = \frac{4G + 1}{3}$

Three or more strata are commonly used as a practical target (Wright & Masters, 1982), but in this package the estimate inherits the same approximation limits as the separation index.

Key references

• Andrich, D. (1978). A rating formulation for ordered response categories. Psychometrika, 43, 561–573.

• Bond, T. G., & Fox, C. M. (2015). Applying the Rasch model (3rd ed.). Routledge.

• Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation of item parameters: Application of an EM algorithm. Psychometrika, 46, 443–459.

• Burnham, K. P., & Anderson, D. R. (2002). Model selection and multimodel inference: A practical information-theoretic approach (2nd ed.). Springer. (AIC / BIC weights and Delta-IC bands used by compare_mfrm().)

• 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. (Source for the lz person-fit statistic implemented in compute_person_fit_indices().)

• Haberman, S. J., & Sinharay, S. (2013). Generalized residuals for general models for contingency tables with application to item response theory. Journal of the American Statistical Association, 108, 1435–1444.

• Eckes, T. (2005). Examining rater effects in TestDaF writing and speaking performance assessments: A many-facet Rasch analysis. Language Assessment Quarterly, 2, 197–221.

• Muraki, E. (1992). A generalized partial credit model: Application of an EM algorithm. Applied Psychological Measurement, 16(2), 159–176. (Source for the bounded GPCM extension used in fit_mfrm(model = "GPCM"), fair_average_table(), and estimate_bias().)

• Muraki, E. (1993). Information functions of the generalized partial credit model. Applied Psychological Measurement, 17(4), 351–363. (Companion paper to Muraki 1992 that derives the GPCM item information identity $I_j(\theta) = D^2 a_j^2 \mathrm{Var}(T \mid \theta)$ via Samejima's (1974) polytomous information formula. This is the canonical reference for compute_information() under bounded GPCM.)

• Samejima, F. (1974). Normal ogive model on the continuous response level in the multidimensional latent space. Psychometrika, 39, 111–121. (General polytomous information formula that Muraki 1993 specializes to the GPCM.)

• Snijders, T. A. B. (2001). Asymptotic null distribution of person fit statistics with estimated person parameter. Psychometrika, 66(3), 331–342. (compute_person_fit_indices() reports a Snijders-style score-projection lz_star for JML fits, conditional on the fitted non-person parameters. For MML/EAP fits lz_star is unavailable because EAP posterior means do not satisfy the ML person-score estimating equation; the old finite-N screen is reported separately as lz_finite_n.)

• Linacre, J. M. (1989). Many-facet Rasch measurement. MESA Press.

• Linacre, J. M. (2002). What do Infit and Outfit, mean-square and standardized mean? Rasch Measurement Transactions, 16(2), 878.

• Masters, G. N. (1982). A Rasch model for partial credit scoring. Psychometrika, 47, 149–174.

• Orlando, M., & Thissen, D. (2000). Likelihood-based item-fit indices for dichotomous item response theory models. Applied Psychological Measurement, 24, 50–64.

• Orlando, M., & Thissen, D. (2003). Further investigation of the performance of S-X2: An item fit index for use with dichotomous item response theory models. Applied Psychological Measurement, 27, 289–298.

• Sinharay, S., Johnson, M. S., & Stern, H. S. (2006). Posterior predictive assessment of item response theory models. Applied Psychological Measurement, 30, 298–321.

• Sinharay, S., & Monroe, S. (2025). Assessment of fit of item response theory models: A critical review of the status quo and some future directions. British Journal of Mathematical and Statistical Psychology, 78, 711–733.

• Wright, B. D., & Masters, G. N. (1982). Rating scale analysis. MESA Press.

• Wright, B. D., & Linacre, J. M. (1994). Reasonable mean-square fit values. Rasch Measurement Transactions, 8(3), 370.

• Wilson, E. B., & Hilferty, M. M. (1931). The distribution of chi-square. Proceedings of the National Academy of Sciences of the United States of America, 17(12), 684-688.

Model selection

RSM vs PCM

The Rating Scale Model (RSM; Andrich, 1978) assumes all levels of the step facet share identical threshold parameters. The Partial Credit Model (PCM; Masters, 1982) allows each level of the step_facet to have its own set of thresholds on the package's shared observed score scale. Use RSM when the rating rubric is identical across all items/criteria; use PCM when category boundaries are expected to vary by item or criterion. In the current implementation, PCM still assumes one common observed score support across the fitted data, so it should not be described as a fully mixed-category model with arbitrary item-specific category counts.

MML vs JML

Marginal Maximum Likelihood (MML) integrates over the person ability distribution using Gauss-Hermite quadrature and does not directly estimate person parameters; person estimates are computed post-hoc via Expected A Posteriori (EAP). Joint Maximum Likelihood (JML) estimates all person and facet parameters simultaneously as fixed effects; "JMLE" remains a backward-compatible alias.

MML is generally preferred for smaller samples because it avoids the incidental-parameter problem of JML. JML does not assume a normal person distribution and can be lighter computationally in some settings, which may be an advantage when the population shape is strongly non-normal.

See fit_mfrm() for usage.

Fixed-calibration scoring after fitting

predict_mfrm_units() and sample_mfrm_plausible_values() score future or partially observed persons on a quadrature grid under the fitted scoring basis. For ordinary MML fits, these summaries inherit the fitted marginal calibration directly. For latent-regression MML fits, they use the fitted one-dimensional conditional normal population model and therefore require one-row-per-person background data for the scored units when the fitted population model includes covariates. Intercept-only latent-regression fits (population_formula = ~ 1) can reconstruct that minimal person table from the scored person IDs. For JML fits, mfrmr uses the fitted facet and step parameters together with a standard normal reference prior introduced only for the post hoc scoring layer. This is useful for practical fixed-scale scoring, but it should still be described as a limited approximation rather than as full ConQuest-style population modeling.

Current ConQuest overlap

The package now includes a first-version latent-regression MML branch, but the overlap with ConQuest should still be described conservatively. The documented overlap is: ordered-response RSM / PCM, one latent dimension, a conditional-normal person population model, and person covariates supplied through an explicit one-row-per-person table and expanded through the package-built model matrix. Categorical person covariates carry fitted levels and contrasts into scoring. This is a scoped overlap, not a claim of broad ConQuest numerical equivalence for arbitrary imported design matrices, multidimensional models, imported design specifications, or the full plausible-values workflow.

Author(s)

Maintainer: Ryuya Komuro [email protected] (ORCID) [copyright holder]

See Also

Useful links:

• https://github.com/Ryuya-dot-com/mfrmr

• Report bugs at https://github.com/Ryuya-dot-com/mfrmr/issues

Examples

mfrm_threshold_profiles()
list_mfrmr_data()


toy <- load_mfrmr_data("example_core")
fit <- fit_mfrm(
  toy,
  person = "Person",
  facets = c("Rater", "Criterion"),
  score = "Score",
  method = "MML",
  model = "RSM",
  quad_points = 7
)
diag <- diagnose_mfrm(fit, diagnostic_mode = "both", residual_pca = "none")
summary(diag)

Description

Tests whether the difficulty of facet levels differs across a grouping variable (e.g., whether rater severity differs for male vs. female examinees, or whether item difficulty differs across rater subgroups).

analyze_dif() is retained for compatibility with earlier package versions. In many-facet workflows, prefer analyze_dff() as the primary entry point.

Usage

analyze_dff(
  fit,
  diagnostics,
  facet,
  group,
  data = NULL,
  focal = NULL,
  method = c("residual", "refit"),
  min_obs = 10,
  p_adjust = "holm"
)

analyze_dif(...)

Arguments

Details

Differential facet functioning (DFF) occurs when the difficulty or severity of a facet element differs across subgroups of the population, after controlling for overall ability. In an MFRM context this generalises classical DIF (which applies to items) to any facet: raters, criteria, tasks, etc.

Differential functioning is a threat to measurement fairness: if Criterion 1 is harder for Group A than Group B at the same ability level, the measurement scale is no longer group-invariant.

Two methods are available:

Residual method (method = "residual"): Uses the existing fitted model's observation-level residuals. For each facet-level $\times$ group cell, the observed and expected score sums are aggregated and a standardized residual is computed as:

$z = \frac{\sum (X_{obs} - E_{exp})}{\sqrt{\sum \mathrm{Var}}}$

Pairwise contrasts between groups compare the mean observed-minus-expected difference for each facet level, with uncertainty summarized by a Welch/Satterthwaite approximation. This method is fast, stable with small subsets, and does not require re-estimation. Because the resulting contrast is not a logit-scale parameter difference, the residual method is treated as a screening procedure rather than an ETS-style classifier.

Refit method (method = "refit"): Subsets the data by group, refits the MFRM model within each subset, anchors all non-target facets back to the baseline calibration when possible, and compares the resulting facet-level estimates using a Welch t-statistic:

$t = \frac{\hat{\delta}_1 - \hat{\delta}_2} {\sqrt{SE_1^2 + SE_2^2}}$

This provides group-specific parameter estimates on a common scale when linking anchors are available, but is slower and may encounter convergence issues with small subsets. ETS categories are reported only for contrasts whose subgroup calibrations retained enough linking anchors to support a common-scale interpretation and whose subgroup precision remained on the package's model-based MML path.

When facet refers to an item-like facet (for example Criterion), this recovers the familiar DIF case. When facet refers to raters or prompts/tasks, the same machinery supports DRF/DPF-style analyses.

For the refit method only, effect size is classified following the ETS (Educational Testing Service) DIF guidelines when subgroup calibrations are both linked and eligible for model-based inference:

• A (Negligible): $|\Delta| <$ 0.43 logits

• B (Moderate): 0.43 $\le |\Delta| <$ 0.64 logits

• C (Large): $|\Delta| \ge$ 0.64 logits

Multiple comparisons are adjusted using Holm's step-down procedure by default, which controls the family-wise error rate without assuming independence. Alternative methods (e.g., "BH" for false discovery rate) can be specified via p_adjust.

Value

An object of class mfrm_dff (with compatibility class mfrm_dif) with:

• dif_table: data.frame of differential-functioning contrasts.

• cell_table: (residual method) per-cell detail table.

• summary: counts by screening or ETS classification.

• group_fits: (refit method) per-group facet estimates.

• config: list with facet, group, method, min_obs, p_adjust settings.

Choosing a method

In most first-pass DFF screening, start with method = "residual". It is faster, reuses the fitted model, and is less fragile in smaller subsets. Use method = "refit" when you specifically want group-specific parameter estimates and can tolerate extra computation. Both methods should yield similar conclusions when sample sizes are adequate ($N \ge 100$ per group is a useful guideline for stable differential-functioning detection).

Interpreting output

• ⁠$dif_table⁠: one row per facet-level x group-pair with contrast, SE, t-statistic, p-value, adjusted p-value, effect metric, and method-appropriate classification. Includes Method, N_Group1, N_Group2, EffectMetric, ClassificationSystem, ContrastBasis, SEBasis, StatisticLabel, ProbabilityMetric, DFBasis, ReportingUse, PrimaryReportingEligible, and sparse columns.

• ⁠$cell_table⁠: (residual method only) per-cell detail with N, ObsScore, ExpScore, ObsExpAvg, StdResidual.

• ⁠$summary⁠: counts by screening result (method = "residual") or ETS category plus linked-screening and insufficient-linking rows (method = "refit").

• ⁠$group_fits⁠: (refit method only) list of per-group facet estimates and subgroup linking diagnostics.

Typical workflow

1. Fit a model with fit_mfrm(). For RSM / PCM fairness review, prefer method = "MML".

2. Run diagnose_mfrm() and, for RSM / PCM, prefer diagnostic_mode = "both" so legacy and strict marginal screens remain visible together.

3. Run analyze_dff(fit, diagnostics, facet = "Criterion", group = "Gender", data = my_data).

4. Inspect ⁠$dif_table⁠ for flagged levels and ⁠$summary⁠ for counts.

5. Use dif_interaction_table() when you need cell-level diagnostics.

6. Use plot_dif_heatmap() or dif_report() for communication.

See Also

fit_mfrm(), estimate_bias(), compare_mfrm(), dif_interaction_table(), plot_dif_heatmap(), dif_report(), subset_connectivity_report(), mfrmr_linking_and_dff

Examples

toy <- load_mfrmr_data("example_bias")

fit <- fit_mfrm(toy, "Person", c("Rater", "Criterion"), "Score",
                 method = "MML", model = "RSM", maxit = 200)
diag <- diagnose_mfrm(fit, residual_pca = "none", diagnostic_mode = "both")
dff <- analyze_dff(fit, diag, facet = "Rater", group = "Group", data = toy)
dff$summary
# Look for: a small `FlaggedPairs` count relative to `Pairs`. Under
#   method = "residual", `ClassificationSystem` is "screening", not
#   ETS. "Screen positive" rows are prompts for substantive review.
head(dff$dif_table[, c("Level", "Group1", "Group2", "Contrast",
                       "Classification", "ClassificationSystem")])
# The residual contrast is an observed-minus-expected average contrast
# between groups. It is useful for screening, but it is not an ETS
# A/B/C logit-delta classification.
dff_refit <- analyze_dff(fit, diag, facet = "Rater", group = "Group",
                         data = toy, method = "refit")
unique(dff_refit$dif_table$ClassificationSystem)
# Look for: "ETS" only when subgroup calibration, linking, and precision
#   checks all support a common-scale model-based contrast.
sc <- subset_connectivity_report(fit, diagnostics = diag)
plot(sc, type = "design_matrix", draw = FALSE)
if ("ScaleLinkStatus" %in% names(dff_refit$dif_table)) {
  unique(dff_refit$dif_table$ScaleLinkStatus)
}
# Look for: "linked" in `ScaleLinkStatus` confirms the focal and
#   reference groups share enough common elements for a comparable
#   contrast; "demoted_*" rows lose linking under the refit branch
#   and should be read as exploratory.

Description

Runs limited classical DIF screens on a long-format score table. The Mantel-Haenszel route uses a generalized Cochran-Mantel-Haenszel test over ordered score categories and total-score strata. The logistic route is a binary logistic-regression screen; for polytomous data it runs only when logistic_threshold is supplied, making the dichotomization explicit.

Usage

analyze_dif_classical(
  x,
  facet,
  group,
  data = NULL,
  score = NULL,
  person = NULL,
  methods = c("mantel_haenszel", "logistic"),
  focal = NULL,
  min_obs = 10L,
  match_bins = 5L,
  p_adjust = "holm",
  logistic_threshold = NULL
)

Arguments

Details

Rows are first collapsed to one mean score per person x group x facet level. The matching variable for a target level is the person's total observed score across the screened facet minus the target-level score, so the target response is not used to condition on itself. When the matching score has more distinct values than match_bins, quantile bins are used as score strata.

The Mantel-Haenszel option forms a $Group \times Score \times MatchStratum$ table and calls stats::mantelhaen.test() without continuity correction. This is a generalized Cochran-Mantel-Haenszel screening p value for ordered score categories. The reported Contrast remains the simple Group2-minus-Group1 mean score difference for direction; it is not a Mantel-Haenszel common odds ratio.

The logistic option fits binary models

$Y \sim MatchScore,$

$Y \sim MatchScore + Group,$

$Y \sim MatchScore * Group.$

The logistic_uniform row is the likelihood-ratio comparison of the first two models. The logistic_nonuniform row is the likelihood-ratio comparison of the second and third models. For polytomous scores, Y is defined as $1(Score \ge logistic_threshold)$; no implicit dichotomization is used.

Value

An object of class mfrm_dff / mfrm_dif with dif_table, cell_table, summary, and config fields.

Scope

This is a classical screening helper, not a replacement for analyze_dff() or SIBTEST. It does not estimate MFRM subgroup parameters, does not use anchors, and does not claim ETS A/B/C classifications.

References

• Mantel, N., & Haenszel, W. (1959). Statistical aspects of the analysis of data from retrospective studies of disease. Journal of the National Cancer Institute, 22, 719-748.

• Holland, P. W., & Thayer, D. T. (1988). Differential item performance and the Mantel-Haenszel procedure. In H. Wainer & H. I. Braun (Eds.), Test validity.

• Swaminathan, H., & Rogers, H. J. (1990). Detecting differential item functioning using logistic regression procedures. Journal of Educational Measurement, 27(4), 361-370.

Examples

toy <- load_mfrmr_data("example_bias")
cls <- analyze_dif_classical(
  toy, facet = "Criterion", group = "Group",
  person = "Person", score = "Score",
  methods = "mantel_haenszel"
)
cls$summary

Description

Analyze practical equivalence within a facet

Usage

analyze_facet_equivalence(
  fit,
  diagnostics = NULL,
  facet = NULL,
  equivalence_bound = 0.5,
  ci_level = 0.95,
  conf_level = NULL
)

Arguments

Details

This function tests whether facet elements (e.g., raters) are similar enough to be treated as practically interchangeable, rather than merely testing whether they differ significantly. This is the key distinction from a standard chi-square heterogeneity test: absence of evidence for difference is not evidence of equivalence.

The function uses existing facet estimates and their standard errors from diagnostics$measures; no re-estimation is performed.

The bundle combines four complementary views:

1. Fixed chi-square test: tests $H_0$: all element measures are equal. A non-significant result is necessary but not sufficient for interchangeability. It is reported as context, not as direct evidence of equivalence.

2. Pairwise TOST (Two One-Sided Tests): for each pair of elements, tests whether the difference falls within $\pm$equivalence_bound. The TOST procedure (Schuirmann, 1987) rejects the null hypothesis of non-equivalence when both one-sided tests are significant at level $\alpha$. A pair is declared "Equivalent" when the TOST p-value < 0.05.

3. BIC-based Bayes-factor heuristic: an approximate screening tool (not full Bayesian inference) that compares the evidence for a common-facet model (all elements equal) against a heterogeneity model (elements differ) via $\mathrm{BF}_{01} \approx \exp((\mathrm{BIC}_{H_1} - \mathrm{BIC}_{H_0}) / 2)$ (Kass & Raftery, 1995). Values > 3 favour the common-facet model; < 1/3 favour heterogeneity.

4. ROPE-style grand-mean proximity: the proportion of each element's normal-approximation confidence distribution that falls within $\pm$equivalence_bound of the weighted grand mean. This is a descriptive proximity summary, not a Bayesian ROPE decision rule around a prespecified null value.

Choosing equivalence_bound: the default of 0.5 logits is a moderate criterion. For high-stakes certification, 0.3 logits may be appropriate; for exploratory or low-stakes contexts, 1.0 logits may suffice. The bound should reflect the smallest difference that would be practically meaningful in your application.

Value

A named list with class mfrm_facet_equivalence.

What this analysis means

analyze_facet_equivalence() is a practical-interchangeability screen. It asks whether facet levels are close enough, under a user-defined logit bound, to be treated as practically similar for the current use case.

What this analysis does not justify

• A non-significant chi-square result is not evidence of equivalence.

• Forest/ROPE displays are descriptive and do not replace the pairwise TOST decision rule.

• The BIC-based Bayes-factor summary is a heuristic screen, not a full Bayesian equivalence analysis.

Interpreting output

Start with summary$Decision, which is a conservative summary of the pairwise TOST results. Then use the remaining tables as context:

• chi_square: is there broad heterogeneity in the facet?

• pairwise: which specific pairs meet the practical-equivalence bound?

• rope / forest: how close is each level to the facet grand mean?

Smaller equivalence_bound values make the criterion stricter. If the decision is "partial_pairwise_equivalence", that means some pairwise contrasts satisfy the practical-equivalence bound but not all of them do.

Decision rule

The final Decision is a pairwise TOST summary rather than a global equivalence proof. If all pairwise contrasts satisfy the practical- equivalence bound, the facet is labeled "all_pairs_equivalent". If at least one, but not all, pairwise contrasts are equivalent, the facet is labeled "partial_pairwise_equivalence". If no pairwise contrasts meet the practical-equivalence bound, the facet is labeled "no_pairwise_equivalence_established". The chi-square, Bayes-factor, and grand-mean proximity summaries are reported as descriptive context.

How to read the main outputs

• summary: one-row pairwise-TOST decision summary and aggregate context.

• pairwise: pair-level TOST detail; use this for the primary inferential read.

• chi_square: broad heterogeneity screen.

• rope / forest: level-wise proximity to the weighted grand mean.

Recommended next step

If the result is borderline or high-stakes, re-run the analysis with a tighter or looser equivalence_bound, then inspect pairwise and plot_facet_equivalence() before deciding how strongly to claim interchangeability.

Typical workflow

1. Fit a model with fit_mfrm().

2. Run analyze_facet_equivalence() for the facet you want to screen.

3. Read summary and chi_square first.

4. Use plot_facet_equivalence() to inspect which levels drive the result.

Output

The returned bundle has class mfrm_facet_equivalence and includes:

• summary: one-row overview with convergent decision

• chi_square: fixed chi-square / separation summary

• pairwise: pairwise TOST detail table

• rope: element-wise ROPE probabilities around the weighted grand mean

• forest: element-wise estimate, confidence interval, and ROPE status

• settings: applied facet and threshold settings

References

Kass, R. E., & Raftery, A. E. (1995). Bayes factors. Journal of the American Statistical Association, 90(430), 773-795.

Schuirmann, D. J. (1987). A comparison of the two one-sided tests procedure and the power approach for assessing the equivalence of average bioavailability. Journal of Pharmacokinetics and Biopharmaceutics, 15(6), 657-680.

See Also

facets_chisq_table(), fair_average_table(), plot_facet_equivalence()

Examples

toy <- load_mfrmr_data("example_core")
fit <- fit_mfrm(toy, "Person", c("Rater", "Criterion"), "Score",
                method = "JML", maxit = 25)
eq <- analyze_facet_equivalence(fit, facet = "Rater")
eq$summary[, c("Facet", "Elements", "Decision", "MeanROPE")]
head(eq$pairwise[, c("ElementA", "ElementB", "Equivalent")])

Description

One-stop audit that combines the nesting, cross-tabulation, ICC, and design-effect reports into a single object. Designed to be reused by the publication-workflow surface: its summary feeds into reporting_checklist(), and its tables are picked up by build_mfrm_manifest() for reproducibility bundles.

Usage

analyze_hierarchical_structure(
  data,
  facets = NULL,
  person = "Person",
  score = "Score",
  compute_icc = TRUE,
  ci_method = c("none", "profile", "boot"),
  ci_level = 0.95,
  ci_boot_reps = 1000L,
  ci_boot_seed = NULL,
  igraph_layout = TRUE,
  icc_ci_method = NULL,
  icc_ci_level = NULL,
  icc_ci_boot_reps = NULL,
  icc_ci_boot_seed = NULL
)

Arguments

Value

A list of class mfrm_hierarchical_structure with:

• nesting: output of detect_facet_nesting().

• crosstabs: list of pairwise observation-count data.frames (long format, suitable for heatmap plotting).

• icc: output of compute_facet_icc() when requested.

• design_effect: output of compute_facet_design_effect() when requested.

• connectivity: named list with bipartite-graph component summary when igraph is available.

• summary: one-row summary used by downstream reporting helpers.

• facets: character vector of facet names that were audited (echoed for downstream reporting helpers that need to label rows by audit scope).

Interpreting output

• nesting: a detect_facet_nesting() object with every facet pair classified as Crossed / Partially / Near-perfectly / Fully nested.

• crosstabs: list of ⁠(LevelA, LevelB, N)⁠ long-format tables, one per facet pair. Plot via plot(x, type = "crosstab", pair = "FacetA__FacetB").

• icc: per-facet variance shares. See compute_facet_icc() for the two-scale interpretation.

• design_effect: Kish (1965) Deff and EffectiveN.

• connectivity: number of bipartite components linking Person x facet levels. A single component is required for a common measurement scale; multiple components indicate a disconnected design.

Typical workflow

1. Optional: fit the MFRM with fit_mfrm().

2. Call analyze_hierarchical_structure(fit) (or on the raw data).

3. Read summary(x) for the condensed view.

4. Feed the object to reporting_checklist() and build_mfrm_manifest() to record the audit in publication bundles. build_apa_outputs() uses the fit-level FacetSampleSizeFlag to add a Methods sentence automatically.

References

McEwen, M. R. (2018). The effects of incomplete rating designs on results from many-facets-Rasch model analyses (Doctoral thesis, Brigham Young University). https://scholarsarchive.byu.edu/etd/6689/

Linacre, J. M. (2026). A User's Guide to FACETS, Version 4.5.0. Winsteps.com. https://www.winsteps.com/facets.htm

Kish, L. (1965). Survey Sampling. New York: Wiley.

Koo, T. K., & Li, M. Y. (2016). A guideline of selecting and reporting intraclass correlation coefficients for reliability research. Journal of Chiropractic Medicine, 15(2), 155-163.

See Also

detect_facet_nesting(), facet_small_sample_audit(), compute_facet_icc(), compute_facet_design_effect(), reporting_checklist(), build_mfrm_manifest(), fit_mfrm().

Examples

toy <- load_mfrmr_data("example_core")
hs <- analyze_hierarchical_structure(toy,
                                     facets = c("Rater", "Criterion"),
                                     compute_icc = FALSE,
                                     igraph_layout = FALSE)
summary(hs)


# Full audit when lme4 and igraph are available.
if (requireNamespace("lme4", quietly = TRUE) &&
    requireNamespace("igraph", quietly = TRUE)) {
  hs_full <- analyze_hierarchical_structure(toy,
                                            facets = c("Rater", "Criterion"))
  summary(hs_full)
  plot(hs_full, type = "icc")
}

Description

Legacy-compatible residual diagnostics can be inspected in two ways:

1. overall residual PCA on the person x combined-facet matrix

2. facet-specific residual PCA on person x facet-level matrices

Usage

analyze_residual_pca(
  diagnostics,
  mode = c("overall", "facet", "both"),
  facets = NULL,
  pca_max_factors = 10L
)

Arguments

Details

The function works on standardized residual structures derived from diagnose_mfrm(). When a fitted object from fit_mfrm() is supplied, diagnostics are computed internally.

Conceptually, this follows the Rasch residual-PCA tradition of examining structure in model residuals after the primary Rasch dimension has been extracted. In mfrmr, however, the implementation is an exploratory many-facet adaptation: it works on standardized residual matrices built as person x combined-facet or person x facet-level layouts, rather than reproducing FACETS/Winsteps residual-contrast tables one-to-one.

Output tables use:

• Component: principal-component index (1, 2, ...)

• Eigenvalue: eigenvalue for each component

• Proportion: component variance proportion

• Cumulative: cumulative variance proportion

For mode = "facet" or "both", by_facet_table additionally includes a Facet column.

summary(pca) is supported through summary(). plot(pca) is dispatched through plot() for class mfrm_residual_pca. Available types include "overall_scree", "facet_scree", "overall_loadings", and "facet_loadings".

Value

A named list with:

• mode: resolved mode used for computation

• facet_names: facets analyzed

• overall: overall PCA bundle (or NULL)

• by_facet: named list of facet PCA bundles

• overall_table: variance table for overall PCA

• by_facet_table: stacked variance table across facets

• errors: named list of any per-facet PCA errors that were caught and turned into NA_real_ rows in the variance tables (e.g., psych::principal() failure on a near-singular residual matrix). The list is empty when every facet PCA succeeded.

Interpreting output

Use overall_table first:

• early components with noticeably larger eigenvalues or proportions suggest stronger residual structure that may deserve follow-up.

Then inspect by_facet_table:

• helps localize which facet contributes most to residual structure.

Finally, inspect loadings via plot_residual_pca() to identify which variables/elements drive each component. For a simulation-calibrated null threshold, use check_residual_dimensionality() and plot_residual_dimensionality().

References

The residual-PCA idea follows the Rasch residual-structure literature, especially Linacre's discussions of principal components of Rasch residuals. The current mfrmr implementation should be interpreted as an exploratory extension for many-facet workflows rather than as a direct reproduction of a single FACETS/Winsteps output table.

• Horn, J. L. (1965). A rationale and test for the number of factors in factor analysis. Psychometrika, 30, 179-185.

• Glorfeld, L. W. (1995). An improvement on Horn's parallel analysis methodology. Educational and Psychological Measurement, 55, 377-393.

• Linacre, J. M. (1998). Structure in Rasch residuals: Why principal components analysis (PCA)? Rasch Measurement Transactions, 12(2), 636.

• Linacre, J. M. (1998). Detecting multidimensionality: Which residual data-type works best? Journal of Outcome Measurement, 2(3), 266-283.

• 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, 717-731.

Typical workflow

1. Fit model and run diagnose_mfrm() with residual_pca = "none" or "both".

2. Call analyze_residual_pca(..., mode = "both").

3. Review summary(pca), then plot scree/loadings.

4. Optionally call check_residual_dimensionality() for a parallel-analysis null threshold.

5. Cross-check with fit/misfit diagnostics before conclusions.

See Also

diagnose_mfrm(), plot_residual_pca(), check_residual_dimensionality(), mfrmr_visual_diagnostics

Examples

toy <- load_mfrmr_data("example_core")
fit <- fit_mfrm(toy, "Person", c("Rater", "Criterion"), "Score", method = "JML", maxit = 25)
diag <- diagnose_mfrm(fit, residual_pca = "both")
pca <- analyze_residual_pca(diag, mode = "both")
pca2 <- analyze_residual_pca(fit, mode = "both")
summary(pca)
p <- plot_residual_pca(pca, mode = "overall", plot_type = "scree", draw = FALSE)
p$data$plot
head(p$data)
head(pca$overall_table)

Description

Re-estimates a many-facet Rasch model on new data while holding selected facet parameters fixed at the values from a previous (baseline) calibration. This is the standard workflow for placing new data onto an existing scale, linking test forms, or carrying a baseline calibration across administration windows.

Usage

anchor_to_baseline(
  new_data,
  baseline_fit,
  person,
  facets,
  score,
  anchor_facets = NULL,
  include_person = FALSE,
  weight = NULL,
  model = NULL,
  method = NULL,
  anchor_policy = "warn",
  ...
)

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

## S3 method for class 'mfrm_anchored_fit'
summary(object, ...)

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

Arguments

Details

This function automates the baseline-anchored calibration workflow:

1. Extracts anchor values from the baseline fit using make_anchor_table().

2. Re-estimates the model on new_data with those anchors fixed via fit_mfrm(..., anchors = anchor_table).

3. Runs diagnose_mfrm() on the anchored fit.

4. Computes element-level differences (new estimate minus baseline estimate) for every common element.

The model and method arguments default to the baseline fit's settings so the calibration framework remains consistent. Elements present in the anchor table but absent from the new data are handled according to anchor_policy: "warn" (default) emits a message, "error" stops execution, and "silent" ignores silently.

The returned drift table is best interpreted as an anchored consistency check. When a facet is fixed through anchor_facets, those anchored levels are constrained in the new run, so their reported differences are not an independent drift analysis. For genuine cross-wave drift monitoring, fit the waves separately and use detect_anchor_drift() on the resulting fits.

Element-level differences are calculated for every element that appears in both the baseline and the new calibration:

$\Delta_e = \hat{\delta}_{e,\text{new}} - \hat{\delta}_{e,\text{base}}$

An element is flagged when $|\Delta_e| > 0.5$ logits or $|\Delta_e / SE_{\Delta_e}| > 2.0$, where $SE_{\Delta_e} = \sqrt{SE_{\mathrm{base}}^2 + SE_{\mathrm{new}}^2}$.

Value

Object of class mfrm_anchored_fit with components:

fit

The anchored mfrm_fit object.

diagnostics

Output of diagnose_mfrm() on the anchored fit.

baseline_anchors

Anchor table extracted from the baseline.

drift

Tibble of element-level drift statistics.

Which function should I use?

• Use anchor_to_baseline() when you have one new dataset and want to place it directly on a baseline scale.

• Use detect_anchor_drift() when you already have multiple fitted waves and want to compare their stability.

• Use build_equating_chain() when you need cumulative offsets across an ordered series of waves.

Interpreting output

• ⁠$drift⁠: one row per common element with columns Facet, Level, Baseline, New, Drift, SE_Baseline, SE_New, SE_Diff, Drift_SE_Ratio, and Flag. Read this as an anchored consistency table. Small absolute differences indicate that the anchored re-fit stayed close to the baseline scale. Flagged rows warrant review, but they are not a substitute for a separate drift study on unanchored common elements.

• ⁠$fit⁠: the full anchored mfrm_fit object, usable with diagnose_mfrm(), measurable_summary_table(), etc.

• ⁠$diagnostics⁠: pre-computed diagnostics for the anchored calibration.

• ⁠$baseline_anchors⁠: the anchor table fed to fit_mfrm(), useful for auditing which elements were constrained.

Typical workflow

1. Fit the baseline model: fit1 <- fit_mfrm(...).

2. Collect new data (e.g., a later administration).

3. Call res <- anchor_to_baseline(new_data, fit1, ...).

4. Inspect summary(res) to confirm the anchored run remains close to the baseline scale.

5. For multi-wave drift monitoring, fit waves separately and pass the fits to detect_anchor_drift() or build_equating_chain().

See Also

fit_mfrm(), make_anchor_table(), detect_anchor_drift(), diagnose_mfrm(), build_equating_chain(), mfrmr_linking_and_dff

Examples

d1 <- load_mfrmr_data("study1")
keep1 <- unique(d1$Person)[1:15]
d1 <- d1[d1$Person %in% keep1, , drop = FALSE]
fit1 <- fit_mfrm(d1, "Person", c("Rater", "Criterion"), "Score",
                 method = "JML", maxit = 15)
d2 <- load_mfrmr_data("study2")
keep2 <- unique(d2$Person)[1:15]
d2 <- d2[d2$Person %in% keep2, , drop = FALSE]
res <- anchor_to_baseline(d2, fit1, "Person",
                          c("Rater", "Criterion"), "Score",
                          anchor_facets = "Criterion")
summary(res)
head(res$drift[, c("Facet", "Level", "Drift", "Flag")])
res$baseline_anchors[1:3, ]

Description

Build APA-style table output using base R structures

Usage

apa_table(
  x,
  which = NULL,
  diagnostics = NULL,
  digits = 2,
  caption = NULL,
  note = NULL,
  bias_results = NULL,
  context = list(),
  whexact = FALSE,
  branch = c("apa", "facets")
)

Arguments

Details

This helper avoids styling dependencies and returns a reproducible base data.frame plus metadata.

Supported which values:

• For mfrm_fit: "summary", "person", "facets", "steps"

• For summary() outputs or mfrm_summary_table_bundle: names listed in build_summary_table_bundle(x)$table_index

• For diagnostics list: "overall_fit", "measures", "fit", "reliability", "facets_chisq", "bias", "interactions", "interrater_summary", "interrater_pairs", "obs"

• For bias-result list: "table", "summary", "chi_sq"

Value

A list of class apa_table with fields:

• table (data.frame)

• which

• caption

• note

• digits

• branch, style

Interpreting output

• table: plain data.frame ready for export or further formatting.

• which: source component that produced the table.

• caption/note: manuscript-oriented metadata stored with the table.

• For fit-based which = "summary", the automatic caption describes the model overview table; use which = "facets" or which = "person" for Table 1-style measurement summaries.

Typical workflow

1. Build table object with apa_table(...).

2. Inspect quickly with summary(tbl).

3. Render with as_kable(tbl) for R Markdown / Quarto or as_flextable(tbl) for Word when those packages are installed.

4. Render base preview via plot(tbl, ...) or export tbl$table.

See Also

fit_mfrm(), diagnose_mfrm(), build_apa_outputs(), reporting_checklist(), mfrmr_reporting_and_apa

Examples

toy <- load_mfrmr_data("example_core")
fit <- fit_mfrm(toy, "Person", c("Rater", "Criterion"), "Score", method = "JML", maxit = 25)
tbl <- apa_table(fit, which = "summary", caption = "Model summary", note = "Toy example")
tbl_facets <- apa_table(fit, which = "summary", branch = "facets")
fit_bundle <- build_summary_table_bundle(summary(fit))
tbl_from_summary <- apa_table(fit_bundle, which = "facet_overview")
summary(tbl)
p <- plot(tbl, draw = FALSE)
p_facets <- plot(tbl_facets, type = "numeric_profile", draw = FALSE)
p$data$plot
p_facets$data$plot
if (interactive()) {
  plot(
    tbl,
    type = "numeric_profile",
    main = "APA Table Numeric Profile (Customized)",
    palette = c(numeric_profile = "#2b8cbe", grid = "#d9d9d9"),
    label_angle = 45
  )
}
tbl$note

Description

Post-hoc shrinkage helper that augments an mfrm_fit with James-Stein / empirical-Bayes shrunk estimates for each non-person facet. The shrinkage variance $\hat{\tau}^2$ is estimated by method of moments from the facet-level point estimates and their standard errors:

$\hat{\tau}^2 = \max\!\left(0, \frac{1}{K}\sum_{j=1}^{K}\hat{\delta}_j^{2} - \overline{\mathrm{SE}^2}\right),$

where the first term is the population variance of the facet point estimates around their known mean of zero (the mfrmr sum-to-zero identification pins the facet mean exactly at 0, so no degree of freedom is consumed by mean estimation). The shrinkage factor is $B_j = \mathrm{SE}_j^2 / (\hat{\tau}^2 + \mathrm{SE}_j^2)$, and the shrunk point / standard error are $\hat{\delta}_j^{EB} = (1 - B_j)\hat{\delta}_j$ and $\mathrm{SE}_j^{EB} = \sqrt{(1 - B_j)\mathrm{SE}_j^2}$. The posterior SE form treats $\hat{\tau}^2$ as known; it omits the Morris (1983, eqs. 4.1-4.2, p. 51) confidence-interval correction $v \cdot \hat{\delta}_j^{2}$ with $v = 2 B_j^2 / (K - r - 2)$, where $r$ is the number of regression coefficients used to model the prior mean (under mfrmr's sum-to-zero pinning, $r = 0$, so the divisor is $K - 2$). This correction adds variance proportional to the squared deviation $\hat{\delta}_j^{2}$, accounting for uncertainty in $\hat{\tau}^2$. Under the equal-variance assumption $\hat{\delta}_j^{2} \approx \hat{\tau}^2$, the omitted variance is on the order of $2 / (K - 2)$ times the reported posterior variance $V(1 - B_j)$, so the true SE is approximately $\sqrt{1 + 2/(K - 2)}$ times the reported ShrunkSE. Magnitudes: SE understated by ~73\ at $K = 8$, ~7\ ShrunkSE as a lower bound rather than a calibrated posterior SE.

Usage

apply_empirical_bayes_shrinkage(
  fit,
  facet_prior_sd = NULL,
  shrink_person = FALSE
)

Arguments

Details

fit$facets$others gains ShrunkEstimate, ShrunkSE, and ShrinkageFactor columns, and fit$shrinkage_report records the per-facet $\hat{\tau}^2$, mean shrinkage, and effective degrees of freedom ($\mathrm{EffectiveDF}_f = \sum_j (1 - B_j)$, which matches the "effective number of parameters" defined by Efron & Morris, 1973). The original Estimate / SE columns are preserved.

Value

The same mfrm_fit, with augmented columns and a new shrinkage_report list entry, and with fit$config$facet_shrinkage set to "empirical_bayes".

Typical workflow

1. Fit the model as usual with fit_mfrm().

2. Call apply_empirical_bayes_shrinkage(fit) when small-N facets are present (see facet_small_sample_audit()).

3. Report both the original and shrunk estimates in the manuscript, citing Efron & Morris (1973). build_apa_outputs() will add the sentence automatically when fit$config$facet_shrinkage is set.

References

Efron, B., & Morris, C. (1973). Combining possibly related estimation problems. Journal of the Royal Statistical Society: Series B, 35(3), 379-402.

Efron, B. (2021). Empirical Bayes: Concepts and methods (Technical report). Department of Statistics, Stanford University. https://efron.ckirby.su.domains/papers/2021EB-concepts-methods.pdf

Morris, C. N. (1983). Parametric empirical Bayes inference: Theory and applications. Journal of the American Statistical Association, 78(381), 47-55.

See Also

fit_mfrm() (which accepts facet_shrinkage directly), facet_small_sample_audit(), compute_facet_icc().

Examples

toy <- load_mfrmr_data("example_core")
fit <- fit_mfrm(toy, "Person", c("Rater", "Criterion"), "Score",
                method = "JML", maxit = 25)
fit_eb <- apply_empirical_bayes_shrinkage(fit)
fit_eb$shrinkage_report
# Look for:
# - `Tau2` is the estimated between-level prior variance per facet.
#   `Tau2 = 0` means the data did not justify any pooling and the
#   shrunken estimates equal the raw estimates (`MeanShrinkage = 0`).
# - `MeanShrinkage` near 0 = little movement, near 1 = heavy pooling
#   toward 0. Small-N facets typically pull values further than
#   well-identified ones.
# - `EffectiveDF` is the implied "effective number of parameters"
#   (Efron & Morris 1973); EffectiveDF much smaller than the row
#   count of the facet means most levels were pooled together.
head(fit_eb$facets$others[, c("Facet", "Level", "Estimate",
                               "ShrunkEstimate", "ShrinkageFactor")])
# Look for: rows where `ShrinkageFactor` is large (close to 1) had
#   their estimates pulled most strongly toward the facet mean (0).

Description

Generic for converting objects to a flextable

Usage

as_flextable(x, ...)

Arguments

Value

A flextable object (concrete return type from the underlying method, e.g. ⁠[as_flextable.apa_table()]⁠ returns a flextable ready for flextable::save_as_docx()).

See Also

as_flextable.apa_table() for the apa_table method; as_kable() for a knitr::kable-targeted alternative; apa_table() for constructing an apa_table in the first place.

Description

Produces a Word / PowerPoint-friendly flextable with the caption and note wired in. Requires flextable (in Suggests).

Usage

## S3 method for class 'apa_table'
as_flextable(x, ...)

Arguments

Value

A flextable object, or a message when flextable is unavailable.

See Also

as_kable.apa_table(), apa_table().

Description

Generic for converting objects to a knitr::kable

Usage

as_kable(x, ...)

Arguments

Value

A knitr::kable object (concrete return type from the underlying method, e.g. ⁠[as_kable.apa_table()]⁠ returns a kableExtra object when the package is installed).

See Also

as_kable.apa_table() for the apa_table method; as_flextable() for a flextable-targeted alternative; apa_table() for constructing an apa_table in the first place.

Description

Renders the table payload for direct inclusion in RMarkdown, Quarto, or HTML reports, wiring the caption and note slots into the standard APA placement (caption above, note below). When kableExtra is installed the note is attached as a footer; otherwise the note is appended as a knitr::asis_output() block.

Usage

## S3 method for class 'apa_table'
as_kable(x, format = c("pipe", "html", "latex"), digits = 3L, ...)

Arguments

Value

A knitr_kable object ready to be printed inline in a report, or a message when knitr is unavailable.

See Also

as_flextable.apa_table(), apa_table().

Description

Convert simulation evaluation objects to data frames.

Usage

## S3 method for class 'mfrm_design_evaluation'
as.data.frame(
  x,
  row.names = NULL,
  optional = FALSE,
  component = c("results", "rep_overview", "design_summary", "overview"),
  ...
)

## S3 method for class 'summary.mfrm_design_evaluation'
as.data.frame(
  x,
  row.names = NULL,
  optional = FALSE,
  component = c("design_summary", "overview"),
  ...
)

## S3 method for class 'mfrm_signal_detection'
as.data.frame(
  x,
  row.names = NULL,
  optional = FALSE,
  component = c("results", "rep_overview", "detection_summary", "overview"),
  ...
)

## S3 method for class 'summary.mfrm_signal_detection'
as.data.frame(
  x,
  row.names = NULL,
  optional = FALSE,
  component = c("detection_summary", "overview"),
  ...
)

## S3 method for class 'mfrm_bias_detection'
as.data.frame(
  x,
  row.names = NULL,
  optional = FALSE,
  component = c("results", "estimates", "reliability", "fit_statistics",
    "fit_summary", "pair_results", "rep_overview", "target_summary",
    "pair_summary", "design_grid"),
  ...
)

## S3 method for class 'summary.mfrm_bias_detection'
as.data.frame(
  x,
  row.names = NULL,
  optional = FALSE,
  component = c("target_summary", "pair_summary", "fit_summary"),
  ...
)

Arguments

Details

The simulation evaluators already store their core outputs as ordinary data frames. These methods make that contract explicit and provide a stable route for write.csv(as.data.frame(x, component = "results"), ...) and custom graphics workflows.

For evaluate_mfrm_design(), the "results" component includes facet-level separation, strata, reliability, fit summaries, and recovery metrics for each design and replication. For evaluate_mfrm_bias_detection(), use "estimates" for fitted measure estimates and "reliability" or "fit_summary" for simulation-derived reliability coefficients.

Value

A base data.frame.

See Also

evaluate_mfrm_design(), evaluate_mfrm_signal_detection(), evaluate_mfrm_bias_detection()

Examples

spec <- build_mfrm_arbitrary_sim_spec(
  n_person = 10,
  facets = c(Rater = 3, Criteria = 2, Task = 2),
  facets_per_person = c(Rater = 2),
  score_levels = 4
)
targets <- data.frame(Rater = "Rater03", Task = "Task02", Effect = -0.5)
eval <- suppressWarnings(evaluate_mfrm_bias_detection(
  spec,
  bias_targets = targets,
  reps = 1,
  fit_method = "JML",
  maxit = 20,
  bias_max_iter = 1,
  seed = 1
))
head(as.data.frame(eval, component = "estimates"))
head(as.data.frame(eval, component = "reliability"))

Description

Returns all facet-level estimates (person and others) in a single tidy data.frame. Useful for quick interactive export: write.csv(as.data.frame(fit), "results.csv").

Usage

## S3 method for class 'mfrm_fit'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

Arguments

Details

This method returns four columns (Facet, Level, Estimate, Extreme) so that the result is easy to inspect, join, or write to disk.

Value

A data.frame with columns Facet, Level, Estimate, and Extreme. The Extreme column is populated for person rows from the extreme-score flag added in 0.1.6 ("Min" / "Max" / NA); non-person facet rows carry NA in that column by design.

Interpreting output

Person estimates are returned with Facet = "Person". All non-person facets are stacked underneath in the same schema.

Typical workflow

1. Fit a model with fit_mfrm().

2. Convert with as.data.frame(fit) for a compact long-format export.

3. Join additional diagnostics later if you need SE or fit statistics.

See Also

fit_mfrm, export_mfrm

Examples

toy <- load_mfrmr_data("example_core")
fit <- fit_mfrm(toy, "Person", c("Rater", "Criterion"), "Score",
                method = "JML", model = "RSM", maxit = 25)
head(as.data.frame(fit))

Description

Coerce residual dimensionality output to a data frame

Usage

## S3 method for class 'mfrm_residual_dimensionality'
as.data.frame(
  x,
  row.names = NULL,
  optional = FALSE,
  component = c("comparison", "observed", "null_distribution"),
  ...
)

Arguments

Value

A data frame.

Description

Audit an exact-overlap ConQuest comparison against an mfrmr overlap bundle

Usage

audit_conquest_overlap(
  bundle,
  conquest_population = NULL,
  conquest_item_estimates = NULL,
  conquest_case_eap = NULL,
  conquest_population_term = "auto",
  conquest_population_estimate = "auto",
  conquest_item_id = "auto",
  conquest_item_estimate = "auto",
  item_id_source = c("auto", "response_var", "level"),
  conquest_case_person = "auto",
  conquest_case_estimate = "auto"
)

Arguments

Details

This helper compares normalized ConQuest output tables against the exact- overlap bundle produced by build_conquest_overlap_bundle(). It is intentionally conservative:

• it does not parse raw ConQuest text output automatically;

• it expects already normalized data frames or output from normalize_conquest_overlap_tables();

• and it reports numerical differences and missing elements without claiming that any fixed tolerance implies software equivalence.

This is the package's external-table audit path. It is distinct from reference_case_benchmark(cases = "synthetic_conquest_overlap_dry_run"), which only round-trips package-native tables through the same normalization and audit contract without executing ConQuest.

The intended workflow is:

1. export an exact-overlap bundle with build_conquest_overlap_bundle();

2. run the narrow matching case in ConQuest;

3. normalize the resulting ConQuest outputs into data frames;

4. pass those tables here to inspect direct differences, centered item agreement, and case-level EAP agreement.

Value

A named list with class mfrm_conquest_overlap_audit.

Output

The returned object has class mfrm_conquest_overlap_audit and includes:

• overall: one-row comparison summary with missing/duplicate/non-numeric attention-item counts and worst-row labels

• population_comparison: parameter-by-parameter comparison table

• item_comparison: centered item-estimate comparison table

• case_comparison: case-level EAP comparison table

• attention_items: missing, malformed, or unmatched elements

• settings: audit settings

• notes: interpretation notes

Interpretation

• Read summary(audit)$audit_scope first to confirm that the result is a supplied-table audit, not raw ConQuest text parsing or a software- equivalence claim.

• Population slopes and sigma2 are intended for direct comparison.

• Item estimates should be interpreted after centering.

• Case estimates should be interpreted as posterior EAP summaries under the fitted population model.

• The overall table reports both mean and maximum absolute differences for compared population, centered item, and case rows. The PopulationMaxAbsParameter, ItemCenteredMaxAbsItem, and CaseMaxAbsPerson columns identify the row where each maximum absolute difference occurs.

• Missing or non-numeric rows in attention_items indicate that the external tables do not yet align cleanly with the exported overlap bundle.

See Also

build_conquest_overlap_bundle(), normalize_conquest_overlap_files(), normalize_conquest_overlap_tables(), reference_case_benchmark()

Examples

bundle <- build_conquest_overlap_bundle()
raw_pop <- data.frame(
  Term = bundle$mfrmr_population$Parameter,
  Est = bundle$mfrmr_population$Estimate
)
raw_item <- data.frame(
  Item = bundle$mfrmr_item_estimates$ResponseVar,
  Est = bundle$mfrmr_item_estimates$Estimate
)
raw_case <- data.frame(
  PID = bundle$mfrmr_case_eap$Person,
  EAP = bundle$mfrmr_case_eap$Estimate
)
normalized <- normalize_conquest_overlap_tables(
  conquest_population = raw_pop,
  conquest_item_estimates = raw_item,
  conquest_case_eap = raw_case,
  conquest_population_term = "Term",
  conquest_population_estimate = "Est",
  conquest_item_id = "Item",
  conquest_item_estimate = "Est",
  conquest_case_person = "PID",
  conquest_case_estimate = "EAP"
)
audit <- audit_conquest_overlap(bundle, normalized)
summary(audit)$summary

Description

Audit and normalize anchor/group-anchor tables

Usage

audit_mfrm_anchors(
  data,
  person,
  facets,
  score,
  anchors = NULL,
  group_anchors = NULL,
  weight = NULL,
  rating_min = NULL,
  rating_max = NULL,
  keep_original = FALSE,
  missing_codes = NULL,
  min_common_anchors = 5L,
  min_obs_per_element = 30,
  min_obs_per_category = 10,
  noncenter_facet = "Person",
  dummy_facets = NULL
)

Arguments

Details

Anchoring (also called "fixing" or scale linking) constrains selected parameter estimates to pre-specified values, placing the current analysis on a previously established scale. This is essential when comparing results across administrations, linking test forms, or monitoring rater drift over time.

This function applies the same preprocessing and key-resolution rules as fit_mfrm(), but returns an audit object so constraints can be checked before estimation. Running the audit first helps avoid estimation failures caused by misspecified or data-incompatible anchors.

Anchor types:

• Direct anchors fix individual element measures to specific logit values (e.g., Rater R1 anchored at 0.35 logits).

• Group anchors constrain the mean of a set of elements to a target value, allowing individual elements to vary freely around that mean.

• When both types overlap for the same element, the direct anchor takes precedence.

Design checks verify that each anchored element has at least min_obs_per_element weighted observations (default 30) and each score category has at least min_obs_per_category (default 10). These thresholds follow standard Rasch sample-size recommendations (Linacre, 1994).

Value

A list of class mfrm_anchor_audit with:

• anchors: cleaned anchor table used by estimation

• group_anchors: cleaned group-anchor table used by estimation

• facet_summary: counts of levels, constrained levels, and free levels

• design_checks: observation-count checks by level/category

• thresholds: active threshold settings used for recommendations

• issue_counts: issue-type counts

• issues: list of issue tables

• recommendations: package-native anchor guidance strings

Interpreting output

• issue_counts/issues: concrete data or specification problems.

• facet_summary: constraint coverage by facet.

• design_checks: whether anchor targets have enough observations.

• recommendations: action items before estimation.

Typical workflow

1. Build candidate anchors (e.g., with make_anchor_table()).

2. Run audit_mfrm_anchors(...).

3. Resolve issues, then fit with fit_mfrm().

See Also

fit_mfrm(), describe_mfrm_data(), make_anchor_table()

Examples

toy <- load_mfrmr_data("example_core")

anchors <- data.frame(
  Facet = c("Rater", "Rater"),
  Level = c("R1", "R1"),
  Anchor = c(0, 0.1),
  stringsAsFactors = FALSE
)
aud <- audit_mfrm_anchors(
  data = toy,
  person = "Person",
  facets = c("Rater", "Criterion"),
  score = "Score",
  anchors = anchors
)
aud$issue_counts
summary(aud)
p_aud <- plot(aud, draw = FALSE)
p_aud$data$plot

Description

Build a bias-cell count report

Usage

bias_count_table(
  bias_results,
  min_count_warn = 10,
  branch = c("original", "facets"),
  fit = NULL
)

Arguments

Details

This helper summarizes how many observations contribute to each bias-cell estimate and flags sparse cells.

Branch behavior:

• "facets": keeps legacy manual-aligned column labels (Sq, ⁠Observd Count⁠, ⁠Obs-Exp Average⁠, ⁠Model S.E.⁠) for side-by-side comparison with external workflows.

• "original": keeps compact field names (Count, BiasSize, SE) for custom QC workflows and scripting.

Value

A named list with:

• table: cell-level counts with low-count flags

• by_facet: named list of counts aggregated by each interaction facet

• by_facet_a, by_facet_b: first two facet summaries (legacy compatibility)

• summary: one-row summary

• thresholds: applied thresholds

• branch, style: output branch metadata

• fit_overview: optional one-row fit metadata when fit is supplied

Interpreting output

• table: cell-level contribution counts and low-count flags.

• by_facet: sparse-cell structure by each interaction facet.

• summary: overall low-count prevalence.

• fit_overview: optional run context (when fit is supplied).

Low-count cells should be interpreted cautiously because bias-size estimates can become unstable with sparse support.

Typical workflow

1. Estimate bias with estimate_bias().

2. Build bias_count_table(...) in desired branch.

3. Review low-count flags before interpreting bias magnitudes.

Further guidance

For a plot-selection guide and a longer walkthrough, see mfrmr_visual_diagnostics and vignette("mfrmr-visual-diagnostics", package = "mfrmr").

Output columns

The table data.frame contains, in the legacy-compatible branch:

FacetA, FacetB

Interaction facet level identifiers; placeholder names for the two interaction facets.

Sq

Sequential row number.

Observd Count

Number of observations for this cell.

Obs-Exp Average

Observed minus expected average for this cell.

Model S.E.

Standard error of the bias estimate.

Infit, Outfit

Fit statistics for this cell.

LowCountFlag

Logical; TRUE when count < min_count_warn.

The summary data.frame contains:

InteractionFacets

Names of the interaction facets.

Cells, TotalCount

Number of cells and total observations.

LowCountCells, LowCountPercent

Number and share of low-count cells.

See Also

estimate_bias(), unexpected_after_bias_table(), build_fixed_reports(), mfrmr_visual_diagnostics

Examples

toy <- load_mfrmr_data("example_bias")
fit <- fit_mfrm(toy, "Person", c("Rater", "Criterion"), "Score", method = "JML", maxit = 25)
diag <- diagnose_mfrm(fit, residual_pca = "none")
bias <- estimate_bias(fit, diag, facet_a = "Rater", facet_b = "Criterion", max_iter = 2)
t11 <- bias_count_table(bias)
t11_facets <- bias_count_table(bias, branch = "facets", fit = fit)
summary(t11)
p <- plot(t11, draw = FALSE)
p2 <- plot(t11, type = "lowcount_by_facet", draw = FALSE)
if (interactive()) {
  plot(
    t11,
    type = "cell_counts",
    draw = TRUE,
    main = "Bias Cell Counts (Customized)",
    palette = c(count = "#2b8cbe", low = "#cb181d"),
    label_angle = 45
  )
}

Description

Bundles the ranked flagged-cells view of a bias-interaction run for downstream printing and plotting. The three sibling reports in this family are intentionally distinct:

• bias_interaction_report() (this one) = FACETS Table 13: a ranked list of interaction cells with t, ⁠bias size⁠, and screening tail area – use when reviewing which ⁠(facet_a, facet_b)⁠ cells deserve follow-up.

• bias_iteration_report() = iteration history / convergence trace for the bias recalibration (FACETS Table 9 territory) – use when diagnosing whether the bias run itself stabilised.

• bias_pairwise_report() = pairwise contrast table for a target facet (FACETS Table 14 territory) – use when comparing levels within a facet while controlling for the other.

Usage

bias_interaction_report(
  x,
  diagnostics = NULL,
  facet_a = NULL,
  facet_b = NULL,
  interaction_facets = NULL,
  max_abs = 10,
  omit_extreme = TRUE,
  max_iter = 4,
  tol = 0.001,
  top_n = 50,
  abs_t_warn = 2,
  abs_bias_warn = 0.5,
  p_max = 0.05,
  sort_by = c("abs_t", "abs_bias", "prob")
)

Arguments

Details

Preferred bundle API for interaction-bias diagnostics. The function can:

• use a precomputed bias object from estimate_bias(), or

• estimate internally from mfrm_fit + facet specification.

Value

A named list with bias-interaction plotting/report components. Class: mfrm_bias_interaction.

Interpreting output

Focus on ranked rows where multiple screening criteria converge:

• large absolute t statistic

• large absolute bias size

• small screening tail area

The bundle is optimized for downstream summary() and plot_bias_interaction() views.

Typical workflow

1. Run estimate_bias() (or provide mfrm_fit here).

2. Build bias_interaction_report(...).

3. Review summary(out) and visualize with plot_bias_interaction().

See Also

estimate_bias(), build_fixed_reports(), plot_bias_interaction()

Examples

toy <- load_mfrmr_data("example_bias")
fit <- fit_mfrm(toy, "Person", c("Rater", "Criterion"), "Score", method = "JML", maxit = 25)
diag <- diagnose_mfrm(fit, residual_pca = "none")
bias <- estimate_bias(fit, diag, facet_a = "Rater", facet_b = "Criterion", max_iter = 2)
out <- bias_interaction_report(bias, top_n = 10)
summary(out)
p_bi <- plot(out, draw = FALSE)
p_bi$data$plot

Description

This report is NOT an alias of bias_interaction_report() despite the similar name. It focuses on the recalibration path of a bias run: iteration table, convergence summary, and orientation audit. Use this to confirm that the bias recalibration itself converged; use bias_interaction_report() to review the ranked flagged cells from the converged run.

Usage

bias_iteration_report(
  x,
  diagnostics = NULL,
  facet_a = NULL,
  facet_b = NULL,
  interaction_facets = NULL,
  max_abs = 10,
  omit_extreme = TRUE,
  max_iter = 4,
  tol = 0.001,
  top_n = 10
)

Arguments

Details

This report focuses on the recalibration path used by estimate_bias(). It provides a package-native counterpart to legacy iteration printouts by exposing the iteration table, convergence summary, and orientation audit in one bundle.

Value

A named list with:

• table: iteration history

• summary: one-row convergence summary

• orientation_audit: interaction-facet sign audit

• settings: resolved reporting options

• direction_note: one-line interpretive note describing which direction the iteration moved (carried from the bias estimator; empty string when the underlying estimator does not emit one)

• recommended_action: one-line recommended action label (e.g. "converged", "increase max_iter"); empty string when the underlying estimator does not emit one

See Also

estimate_bias(), bias_interaction_report(), build_fixed_reports()

Examples

toy <- load_mfrmr_data("example_bias")
fit <- fit_mfrm(toy, "Person", c("Rater", "Criterion"), "Score", method = "JML", maxit = 25)
diag <- diagnose_mfrm(fit, residual_pca = "none")
out <- bias_iteration_report(fit, diagnostics = diag, facet_a = "Rater", facet_b = "Criterion")
summary(out)

Description

Build a pairwise contrast table that, for a chosen target facet (e.g. raters), compares each pair of target-facet levels while holding a context facet (e.g. items / criteria) constant. This is the FACETS Table 14 view: it answers "is rater A consistently more severe than rater B on the same items?" rather than "which (rater, item) cell has the largest local bias?" – the latter is covered by bias_interaction_report().

Usage

bias_pairwise_report(
  x,
  diagnostics = NULL,
  facet_a = NULL,
  facet_b = NULL,
  interaction_facets = NULL,
  max_abs = 10,
  omit_extreme = TRUE,
  max_iter = 4,
  tol = 0.001,
  target_facet = NULL,
  context_facet = NULL,
  top_n = 50,
  p_max = 0.05,
  sort_by = c("abs_t", "abs_contrast", "prob")
)

Arguments

Details

This helper exposes the pairwise contrast table that was previously only reachable through fixed-width output generation. It is available only for 2-way interactions. The pairwise contrast statistic uses a Welch/Satterthwaite approximation and is labeled as a Rasch-Welch comparison in the output metadata.

Value

A named list with:

• table: pairwise contrast rows

• summary: one-row contrast summary

• orientation_audit: interaction-facet sign audit

• settings: resolved reporting options

• direction_note: one-line interpretive note describing the dominant pairwise-contrast direction (carried from the underlying bias estimator; empty string when not applicable)

• recommended_action: one-line recommended-action label (e.g. routing the user to follow-up review of the largest flagged pairs); empty string when the underlying estimator does not emit one

Interpreting output

• table: one row per ordered (target_level_1, target_level_2) pair, with Bias_diff, SE_diff, t_diff, df_diff, p_diff, and the underlying per-level bias rows. Rows are sorted so that the largest-magnitude ⁠|t_diff|⁠ rises to the top.

• summary: one-row screening summary with MaxAbsBiasDiff, MaxAbsT, Significant (count of flagged pairs at p_max), BonferroniSignificant, and HolmSignificant.

• orientation_audit carries the same facet-orientation sign audit as the parent estimate_bias() run.

• The SE caveat below applies: read Significant / BonferroniSignificant as a screening triage, not as formal inferential tests.

Typical workflow

1. Fit and diagnose the model.

2. Run estimate_bias() to get the underlying interaction effects.

3. Pass that result to bias_pairwise_report() for the rater-pair contrast table.

4. Use summary(out)$MaxAbsT and the top rows of out$table to flag rater-pair systematic differences for follow-up review.

5. For the ranked flagged-cells view (which (rater, item) pairs have the largest local bias), use bias_interaction_report() on the same estimate_bias() output.

Standard-error caveat

The contrast standard error is computed as SE(b_i - b_j) = sqrt(SE_i^2 + SE_j^2) – the independence approximation. For same-facet bias values that share a sum-to-zero identification, Cov(b_i, b_j) < 0, so the true contrast variance is SE_i^2 + SE_j^2 - 2 * Cov(b_i, b_j), which is smaller than the reported value. The reported t-statistics and p-values are therefore conservative for same-facet contrasts (the true significance is higher than reported). For across-facet contrasts the covariance term is approximately zero and the approximation is appropriate. Use the report as a screening / triage table; for inferential claims that hinge on a marginally-significant same-facet contrast, follow up with a contrast that uses the full parameter covariance.

References

• Linacre, J. M. (1989). Many-Facet Rasch Measurement. MESA Press.

• Eckes, T. (2005). Examining rater effects in TestDaF writing and speaking performance assessments: A many-facet Rasch analysis. Language Assessment Quarterly, 2(3), 197-221.

• Myford, C. M., & Wolfe, E. W. (2003). Detecting and measuring rater effects using many-facet Rasch measurement: Part I. Journal of Applied Measurement, 4(4), 386-422.

• Myford, C. M., & Wolfe, E. W. (2004). Detecting and measuring rater effects using many-facet Rasch measurement: Part II. Journal of Applied Measurement, 5(2), 189-227.

See Also

estimate_bias(), bias_interaction_report(), build_fixed_reports()

Examples

toy <- load_mfrmr_data("example_bias")
fit <- fit_mfrm(toy, "Person", c("Rater", "Criterion"), "Score", method = "JML", maxit = 25)
diag <- diagnose_mfrm(fit, residual_pca = "none")
out <- bias_pairwise_report(fit, diagnostics = diag, facet_a = "Rater", facet_b = "Criterion")
s <- summary(out)
s$summary
# Look for: `MaxAbsBiasDiff` < ~0.5 logits and `Significant = 0` mean
#   no rater pair contrasts above the screen. The `BonferroniSignificant`
#   / `HolmSignificant` columns count pairs that survive multiple-
#   testing correction; both being 0 is a stronger "no rater-pair
#   inconsistency" signal than the raw screen-positive count alone.
head(out$table)
# Look for: top rows with `|t_diff|` > 2 and |Bias_diff| > 0.5 logits
#   warrant content-review of the two raters' scoring conventions on
#   the conditioning context facet (e.g. compare their item-level
#   marks for systematic strictness/leniency patterns).

Description

Build an APA reporting bundle from model results

Usage

build_apa_outputs(
  fit,
  diagnostics,
  bias_results = NULL,
  context = list(),
  whexact = FALSE
)

Arguments

Details

context is an optional named list for narrative customization. Frequently used fields include:

• assessment, setting, scale_desc

• rater_training, raters_per_response

• rater_facet (used for targeted reliability note text)

• line_width (optional text wrapping width for report_text; default = 92)

Output text includes residual-PCA screening commentary if PCA diagnostics are available in diagnostics.

build_apa_outputs() is the single front-door helper for concise manuscript reporting output. It intentionally reuses the same facts exposed by summary(fit), summary(diagnostics), and companion reporting helpers, but the object returned here is the paper-facing bundle: printing it shows the compact Method / Results narrative, whereas summary(apa) is a QA checklist for completeness, convergence/precision readiness, and wording alignment.

For bounded GPCM, this helper returns a caveated APA scaffold. It uses the GPCM-specific model wording and carries a support_status table plus a caveat field. Keep fair-average and bias language at the screening tier, and do not describe conditional fair-average SEs as full joint-uncertainty intervals. Use gpcm_capability_matrix() as the formal boundary statement for that branch.

By default, report_text includes:

• model/data design summary (N, facet counts, scale range)

• optimization/convergence metrics (Converged, Iterations, LogLik, AIC, BIC)

• anchor/constraint summary (noncenter_facet, anchored levels, group anchors, dummy facets)

• latent-regression population-model wording when fit has an active population_formula

• category/threshold diagnostics (including disordered-step details when present)

• overall fit, misfit count, and top misfit levels

• facet reliability/separation, residual PCA summary, and bias-screen counts

Value

An object of class mfrm_apa_outputs with:

• report_text: APA-style Method/Results draft prose

• table_figure_notes: consolidated draft notes for tables/visuals

• table_figure_captions: draft caption candidates without figure numbering

• section_map: package-native section table for manuscript assembly

• contract: structured APA reporting contract used for downstream checks

Interpreting output

• report_text: manuscript-draft narrative covering Method (model specification, estimation, convergence) and Results (global fit, facet separation/reliability, misfit triage, category diagnostics, residual-PCA screening, bias screening). Written in third-person past tense following APA 7th edition conventions, but still intended for human review.

• table_figure_notes: reusable draft note blocks for table/figure appendices.

• table_figure_captions: draft caption candidates aligned to generated outputs.

• active latent-regression fits add a population-model section and Table 5 notes/captions that distinguish conditional-normal coefficient reporting from post hoc regression on EAP/MLE scores.

When bias results or PCA diagnostics are not supplied, those sections are omitted from the narrative rather than producing placeholder text.

Typical workflow

1. Build diagnostics (and optional bias results). For RSM / PCM reporting runs, prefer an MML fit and diagnose_mfrm(..., diagnostic_mode = "both").

2. Run build_apa_outputs(...).

3. Check summary(apa) for completeness and analysis-readiness flags.

4. Print the returned object to view the concise manuscript narrative (cat(apa$report_text) is equivalent for scripted output).

5. Insert apa$report_text, apa$section_map, and note/caption fields into manuscript drafts after checking the listed cautions.

Context template

A minimal context list can include fields such as:

• assessment: name of the assessment task

• setting: administration context

• scale_desc: short description of the score scale

• rater_facet: rater facet label used in narrative reliability text

Input validation

fit must be an mfrm_fit object from fit_mfrm(). diagnostics must be an mfrm_diagnostics object from diagnose_mfrm(). context must be a list (use NULL or list() for no extra context). If supplied, bias_results must come from estimate_bias() or another package-native bias helper that provides a table component.

See Also

build_visual_summaries(), estimate_bias(), reporting_checklist(), mfrmr_reporting_and_apa

Examples

# Fast smoke run: a JML fit and a legacy diagnostic let us build the
# APA bundle and confirm `report_text` is non-empty in well under
# a second.
toy <- load_mfrmr_data("example_core")
fit_quick <- fit_mfrm(toy, "Person", c("Rater", "Criterion"), "Score",
                      method = "JML", maxit = 15)
diag_quick <- diagnose_mfrm(fit_quick, residual_pca = "none",
                             diagnostic_mode = "legacy")
apa_quick <- build_apa_outputs(fit_quick, diag_quick)
nchar(apa_quick$report_text) > 0


fit <- fit_mfrm(toy, "Person", c("Rater", "Criterion"), "Score",
                method = "MML", maxit = 200)
diag <- diagnose_mfrm(fit, residual_pca = "both", diagnostic_mode = "both")
apa <- build_apa_outputs(
  fit,
  diag,
  context = list(
    assessment = "Toy writing task",
    setting = "Demonstration dataset",
    scale_desc = "0-2 rating scale",
    rater_facet = "Rater"
  )
)
s_apa <- summary(apa)
s_apa$overview[, c("DraftContractPass", "AnalysisReady")]
apa$section_map[, c("SectionId", "Available", "SentenceCount")]
# Look for: `DraftContractPass = TRUE` before using the generated prose.
#   `AnalysisReady = FALSE` means convergence or formal precision still
#   needs review before submission.
chk <- reporting_checklist(fit, diagnostics = diag)
head(chk$checklist[, c("Section", "Item", "DraftReady", "NextAction")])
# Look for: rows with `DraftReady = "yes"` are ready to draft from
#   under the documented caveats. `"no"` rows tell you which helper / setting
#   needs to run before that paragraph can be drafted, via
#   `NextAction`. Aim for every Visual Displays / Reliability /
#   Diagnostics row to be `"yes"` before submitting.
apa
apa$section_map[, c("SectionId", "Available")]

Description

Build a scoped ConQuest-overlap bundle

Usage

build_conquest_overlap_bundle(
  fit = NULL,
  case = c("synthetic_latent_regression"),
  output_dir = NULL,
  prefix = "conquest_overlap",
  overwrite = FALSE,
  quad_points = 7L,
  maxit = 40L,
  reltol = 1e-06
)

Arguments

Details

This helper prepares a narrow ConQuest comparison bundle for an RSM / PCM latent-regression MML fit and records the mfrmr-side tables to compare after an external ConQuest run. The supported overlap is intentionally narrow:

• ordered-response RSM / PCM only;

• binary responses only;

• exactly one non-person facet, treated as the item facet;

• active latent-regression MML;

• exactly one numeric person covariate beyond the intercept;

• complete person-by-item rectangular data.

The returned bundle standardizes the responses to ⁠{0, 1}⁠, pivots them to a one-row-per-person wide CSV, stores the corresponding person covariates, and records the mfrmr estimates that should be compared externally.

The conquest_command component is a conservative starting template, not a guaranteed version-invariant automation. The conquest_output_contract component records which requested external output should feed each normalized audit table. Use normalize_conquest_overlap_files() or normalize_conquest_overlap_tables() and then audit_conquest_overlap() only after the matching ConQuest run has been executed externally and the relevant output tables have been extracted. The bundle and command template alone are not external validation evidence.

Value

A named list with class mfrm_conquest_overlap_bundle.

Comparison targets

• regression slope: compare directly;

• residual variance sigma2: compare directly;

• item estimates: compare after centering because the Rasch location origin remains constraint-dependent;

• case EAP estimates: compare as posterior summaries under the fitted population model.

Output

The returned object has class mfrm_conquest_overlap_bundle and includes:

• summary: one-row scope summary with posterior-basis and population-model audit fields

• comparison_targets: comparison rules for the exported tables

• conquest_output_contract: requested ConQuest outputs and audit handoff

• response_long: long-format binary response data used by the bundle

• response_wide: wide CSV-ready response matrix for the ConQuest template

• person_data: one-row-per-person covariate table

• item_map: mapping from exported response columns to original item levels

• mfrmr_population: fitted population-model coefficients plus sigma2

• mfrmr_item_estimates: fitted item estimates with centered values

• mfrmr_case_eap: posterior EAP summaries for the fitted persons

• conquest_command: conservative ConQuest command template

• written_files: file inventory when output_dir is supplied

• settings: bundle settings

• notes: interpretation notes

See Also

normalize_conquest_overlap_files(), normalize_conquest_overlap_tables(), audit_conquest_overlap(), reference_case_benchmark(), build_mfrm_replay_script(), export_mfrm_bundle()

Examples

bundle <- build_conquest_overlap_bundle()
bundle$summary[, c("Case", "Facet", "Covariate", "Persons", "Items")]
summary(bundle)$conquest_command_scope
summary(bundle)$conquest_output_contract
cat(substr(bundle$conquest_command, 1, 120))