mfrmr Workflow

This vignette outlines a reproducible workflow for:

For a plot-first companion guide, see the separate mfrmr-visual-diagnostics vignette.

For a faster preliminary run without changing the final analysis target:

MML and Diagnostic Modes

mfrmr treats MML and JML differently on purpose.

For RSM and PCM, diagnostics now expose two distinct evidence paths:

Strict marginal diagnostics are screening-oriented. Use summary(diag)$diagnostic_basis to separate the legacy residual evidence from the strict marginal evidence rather than pooling them into one decision.

Load Data

library(mfrmr)

list_mfrmr_data(details = TRUE)[, c("Key", "PrimaryUse", "Design", "CountBasis")]

data("ej2021_study1", package = "mfrmr")
head(ej2021_study1)

study1_alt <- load_mfrmr_data("study1")
identical(names(ej2021_study1), names(study1_alt))

Applied Runnable Example

Start with the packaged example_operational dataset. It is intentionally compact but uses a connected two-rater assignment, unequal rater workloads, and six planned omissions represented by absent long-format rows rather than NA or sentinel scores. This makes the main tutorial closer to an applied rating design without using empirical records. The same object is also available via data("mfrmr_example_operational", package = "mfrmr"). A separate score-free assignment roster lets the pre-fit review identify those six omissions without guessing which cells should have existed. example_core remains available as an idealized complete-crossing example for fast help-page checks.

data("mfrmr_example_operational", package = "mfrmr")
data("mfrmr_example_operational_design", package = "mfrmr")
toy <- mfrmr_example_operational

data_review_toy <- describe_mfrm_data(
  data = toy,
  person = "Person",
  facets = c("Rater", "Criterion"),
  score = "Score",
  rating_min = 1,
  rating_max = 4,
  expected_design = mfrmr_example_operational_design
)
data_summary_toy <- summary(data_review_toy)
data_summary_toy$structural_missingness
data_summary_toy$design_connectivity

fit_toy <- fit_mfrm(
  data = toy,
  person = "Person",
  facets = c("Rater", "Criterion"),
  score = "Score",
  method = "MML",
  model = "RSM"
)
diag_toy <- diagnose_mfrm(
  fit_toy,
  residual_pca = "none",
  diagnostic_mode = "both",
  fit_df_method = "both"
)

# Fast fit-only summary: this does not compute diagnostics.
fit_summary_toy <- summary(fit_toy, profile = "fit", detail = "brief")

# Comprehensive review, reusing diagnostics already computed above.
facets_summary_toy <- summary(
  fit_toy,
  profile = "facets",
  detail = "brief",
  diagnostics = diag_toy
)
res_toy <- facets_summary_toy$results

fit_summary_toy$overview
fit_summary_toy$readiness
fit_summary_toy$data_review
summary(diag_toy)$overview
facets_summary_toy

# Required first fitted-scale figure.
plot(res_toy, type = "wright", preset = "publication", show_ci = TRUE, top_n = Inf)

# Optional FACETS-style ruler. Replace these examples with the study rubric.
rubric_labels <- c(
  "1" = "Level 1",
  "2" = "Level 2",
  "3" = "Level 3",
  "4" = "Level 4"
)
plot(
  res_toy,
  type = "wright",
  renderer = "facets",
  category_labels = rubric_labels,
  show_ci = FALSE,
  preset = "publication"
)
# Setting show_ci = TRUE on this FACETS-style ruler is available as a
# deliberate hybrid: FACETS ruler grammar plus mfrmr uncertainty intervals.

# Optional follow-up: Infit on x, measure on y; persons are explicit opt-in.
plot(
  res_toy,
  type = "fit_pathway",
  fit_stat = "Infit",
  include_person = TRUE,
  top_n_person = 12,
  person_labels = "none",
  facet_labels = "flagged",
  preset = "publication"
)

Optimizer code zero is only one numerical signal. fit_mfrm() also checks the terminal gradient; when the initial direct or hybrid solution stops with a larger gradient and reltol <= 1e-9, it runs a bounded sequence of warm-started polishing stages and retains the best non-worsening candidate under the recorded selection rule. Inspect fit_toy$opt$optimizer_polish$Stages when Numerical is not pass.

maxit is only a computational ceiling. It must not be increased until a preferred coefficient, fit statistic, or expected conclusion appears. For a final analysis, prespecify the estimator and controls and start from the package default maxit = 400 unless the analysis protocol states otherwise. If a run reaches ConvergenceStatus = "iteration_limit", keep it review-only and repeat the same data, model, method, anchors, optimizer, tolerance, and quadrature rule with the next ceiling in a prespecified sequence. Interpret only a run with Converged = TRUE, InferenceReady = TRUE, and Numerical = pass. If two separately ready runs differ materially, investigate numerical stability instead of selecting the preferred result. Report the requested ceiling, actual evaluations, convergence reason, and terminal gradient.

InferenceReady describes this numerical gate only: Data, Design, Stability, Diagnostics, and Reporting remain separate rows in fit_summary_toy$readiness. For example, a disconnected design remains a reporting hold even if its optimizer converges. Conversely, ready_for_diagnostics_and_reporting_follow_up means that fitting completed and diagnostic review is the next stage; it is neither an optimization failure nor a claim that the analysis is already manuscript-ready.

expected_design is a declared roster, not an inferred complete crossing. When no roster is available, the structural-missingness status is "not_declared"; the package does not label unassigned cells as missing.

The facets profile is FACETS-style organization, not evidence that FACETS was run and not a numerical-equivalence claim. Its brief print is selective and does not print person identifiers; full tables remain in facets_summary_toy$results$tables. For a report-oriented result set, use profile = "reporting". To inspect availability without allowing diagnostic computation, use compute = "never"; requested dependent sections are then recorded as not_computed:

reporting_summary_toy <- summary(
  fit_toy,
  profile = "reporting",
  diagnostics = diag_toy
)

availability_only <- summary(
  fit_toy,
  profile = "facets",
  compute = "never"
)
availability_only$section_status

Availability and interpretability are intentionally separate. Review res_toy$readiness and the InterpretationStatus columns in res_toy$plot_map before treating an available plot as a final result. Plots remain available during a numerical, data, design, or stability review, but they warn and carry REVIEW ONLY in the returned subtitle and drawn title until the readiness issue is resolved.

Bias/DIF, residual PCA, and anchor-drift/linking analyses are deliberately not auto-run by any summary profile. They require explicit contrasts, diagnostic settings, or multi-fit designs.

#>   Model Method MethodUsed   N Persons Facets FacetInteractions
#> 1   RSM    MML        MML 282      48      2                 0
#>   InteractionParameters InteractionCells InteractionSparseCells Categories
#> 1                     0                0                      0          4
#>              LogLik              AIC              BIC Converged InferenceReady
#> 1 -347.239803707285 712.479607414569 745.256771053012      TRUE           TRUE
#>   Iterations      IterationsBasis MMLEngineRequested MMLEngineUsed
#> 1         28 function_evaluations             direct        direct
#>                            MMLEngineDetail EMIterations EMConverged
#> 1 Direct marginal likelihood optimization.                         
#>   EMRelativeChange OptimizerMethod OptimizerInitialMethod OptimizerPolished
#> 1                         L-BFGS-B               L-BFGS-B              TRUE
#>   OptimizerPolishSucceeded OptimizerPolishStages RequestedReltol
#> 1                     TRUE                     2           1e-09
#>   EffectiveReltol  OptimizerFactr       OptimizerPgtol
#> 1           1e-13 450.35996273705 1.49011611938477e-08
#>   InitialTerminalGradientSupNorm ConvergenceCode   ConvergenceBasis
#> 1           0.000693846262430853               0 optimizer_gradient
#>   ConvergenceStatus ConvergenceReason ConvergenceSeverity
#> 1         converged     tolerance_met                pass
#>                                ConvergenceMessage
#> 1 CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH
#>                        ConvergenceDetail ReviewableWarning
#> 1 Optimizer returned convergence code 0.             FALSE
#>   GradientReviewTolerance FunctionEvaluations GradientEvaluations
#> 1                   1e-04                  28                  28
#>   TerminalGradientSupNorm  TerminalGradientRMS FacetSampleSizeFlag
#> 1    1.58249274520017e-05 8.46761078292423e-06            standard
#>   FacetMinLevelN FacetSparseCount ExtremeHighN ExtremeLowN
#> 1             38                0            0           0
#>   Observations Persons Facets Categories Subsets ResidualPCA DiagnosticMode
#> 1          282      48      2          4       1        none           both
#>   Method PrecisionTier MarginalFit
#> 1    MML   model_based   available
#>   Component
#> 1      name
#> 2      data

The same fit can then move through the recommended first reporting workflow:

report_toy <- mfrm_report(res_toy, style = "qc")

summary(res_toy)$next_actions
summary(report_toy)$overview

# This is a controlled analysis archive, not a deidentified shareable export.
export_dir <- file.path(tempdir(), "mfrmr-workflow-export")
export_toy <- export_mfrm_results(
  res_toy,
  output_dir = export_dir,
  include = c("default", "report"),
  overwrite = TRUE,
  acknowledge_sensitive = TRUE
)
head(export_toy$written_files)

The acknowledgement suppresses the warning only; it does not redact person identifiers, person-level results, local paths, or the complete RDS object. Review every exported file under the study’s data-handling policy before sharing.

#>   Priority               Area
#> 1        1           Overview
#> 2        2             Triage
#> 3        2         Wright map
#> 4        3        Diagnostics
#> 5        4 Visual diagnostics
#> 6        5        Fit pathway
#> 7        5          Precision
#> 8        6          Reporting
#> 9       11             Tables
#>                                                                      Action
#> 1                                         Read the compact results summary.
#> 2                            Read the first-screen triage before branching.
#> 3                   Create and inspect the required shared-logit scale map.
#> 4                    Review diagnostic key warnings before report drafting.
#> 5                     Open the QC dashboard after reviewing the Wright map.
#> 6 Review Infit against measure, including selected person rows when useful.
#> 7        Inspect fit, separation, reliability, and ZSTD wording boundaries.
#> 8        Use the reporting checklist as a guide for preparing a manuscript.
#> 9                            Create an appendix-ready summary-table bundle.
#>                                                                                                                                                                     Route
#> 1                                                                                                                                                            summary(res)
#> 2                                                                                                                                                     summary(res)$triage
#> 3                                                                                         plot(res, type = "wright", preset = "publication", show_ci = TRUE, top_n = Inf)
#> 4                                                                                                                                   summary(res$diagnostics)$key_warnings
#> 5                                                                                                                          plot(res, type = "qc", preset = "publication")
#> 6 plot(res, type = "fit_pathway", fit_stat = "Infit", include_person = TRUE, top_n_person = 12, person_labels = "none", facet_labels = "flagged", preset = "publication")
#> 7                                                                                                                                summary(res$components$precision_review)
#> 8                                                                                                                             summary(res$components$reporting_checklist)
#> 9                                                                                                                                         build_summary_table_bundle(res)
#>                                                                                                                                                 Reason
#> 1                                                           Confirms input mode, model, method, section status, table coverage, and available figures.
#> 2                             Triage orders unavailable, review, information, and OK signals across diagnostics, tables, plots, and reporting outputs.
#> 3 The Wright map is the primary fitted-scale figure: compare person targeting with facet levels and step thresholds before branching into diagnostics.
#> 4                                            Diagnostic warnings identify the highest-priority fit, precision, residual, or category follow-up checks.
#> 5                                                            The QC dashboard gives a focused follow-up view of fit, residual, and category summaries.
#> 6                                          This follow-up separates measure uncertainty from fit displacement while keeping person inclusion explicit.
#> 7                                         Precision review keeps fit-size, standardized fit, and separation evidence in separate reporting categories.
#> 8                                                                                Checklist rows identify report-ready, missing, and caveated sections.
#> 9                                                                   The bundle exposes table roles, plot readiness, and conservative appendix presets.
#>   Style OverallStatus     FirstAction ReviewAreas NotComputedAreas CaveatAreas
#> 1    qc        review Start with Fit.           1                0           0
#>   OptionalAreas UnavailableAreas OkAreas
#> 1             3                0       1
#>                                                       SourceInclude
#> 1 fit, diagnostics, tables, precision, reporting, categories, plots
#>                 Component Format                                      Path Note
#> 1        summary_overview    csv        mfrmr_results_summary_overview.csv     
#> 2          summary_status    csv          mfrmr_results_summary_status.csv     
#> 3 summary_component_index    csv mfrmr_results_summary_component_index.csv     
#> 4     summary_table_index    csv     mfrmr_results_summary_table_index.csv     
#> 5        summary_plot_map    csv        mfrmr_results_summary_plot_map.csv     
#> 6          summary_triage    csv          mfrmr_results_summary_triage.csv     
#>            DataHandling
#> 1 review_before_sharing
#> 2 review_before_sharing
#> 3 review_before_sharing
#> 4 review_before_sharing
#> 5 review_before_sharing
#> 6 review_before_sharing

Diagnostics and Reporting

t4_toy <- unexpected_response_table(
  fit_toy,
  diagnostics = diag_toy,
  abs_z_min = 1.5,
  prob_max = 0.4,
  top_n = 10
)
t12_toy <- fair_average_table(fit_toy, diagnostics = diag_toy)
t13_toy <- bias_interaction_report(
  estimate_bias(fit_toy, diag_toy,
                facet_a = "Rater", facet_b = "Criterion",
                max_iter = 2),
  top_n = 10
)

class(summary(t4_toy))
class(summary(t12_toy))
class(summary(t13_toy))

names(plot(t4_toy, draw = FALSE))
names(plot(t12_toy, draw = FALSE))
names(plot(t13_toy, draw = FALSE))

chk_toy <- reporting_checklist(fit_toy, diagnostics = diag_toy)
subset(
  chk_toy$checklist,
  Section == "Visual Displays",
  c("Item", "DraftReady", "NextAction")
)
#>                      Object        SummaryClass
#> 1 unexpected_response_table summary.mfrm_bundle
#> 2        fair_average_table summary.mfrm_bundle
#> 3   bias_interaction_report summary.mfrm_bundle
#>                      Object Components
#> 1 unexpected_response_table name, data
#> 2        fair_average_table name, data
#> 3   bias_interaction_report name, data
#>                                  Item DraftReady
#> 1                          Wright map       TRUE
#> 2                QC / facet dashboard       TRUE
#> 3                Residual PCA visuals      FALSE
#> 4 Connectivity / design-matrix visual       TRUE
#> 5  Inter-rater / displacement visuals       TRUE
#> 6             Strict marginal visuals      FALSE
#> 7                  Bias / DIF visuals      FALSE
#> 8      Precision / information curves       TRUE
#> 9                Fit/category visuals       TRUE
#>                                                                                                                      NextAction
#> 1                                      Include a Wright map when the manuscript benefits from a shared-scale targeting display.
#> 2                     Use the dashboard as a first-pass triage view, then move to the specific follow-up plot behind each flag.
#> 3                                         Run residual PCA if you want scree/loadings visuals for residual-structure follow-up.
#> 4                                                       Use the design-matrix view to support linkage and comparability claims.
#> 5                                       Use displacement and inter-rater views to localize QC issues after dashboard screening.
#> 6 Treat strict marginal plots as exploratory corroboration screens, then corroborate with design review and legacy diagnostics.
#> 7                                                        Run bias or DIF screening before discussing interaction-level visuals.
#> 8                                Use information curves to describe precision across theta when that is the reporting question.
#> 9                                        Use category curves and fit visuals as local descriptive follow-up after QC screening.

Fit and Diagnose with Full Data

For a larger sparse synthetic illustration, use the packaged Study 1 dataset:

fit <- fit_mfrm(
  data = ej2021_study1,
  person = "Person",
  facets = c("Rater", "Criterion"),
  score = "Score",
  method = "MML",
  model = "RSM",
  quad_points = 7
)

diag <- diagnose_mfrm(
  fit,
  residual_pca = "none",
  diagnostic_mode = "both",
  fit_df_method = "both"
)

summary(fit, profile = "fit", detail = "brief")
summary(diag)

# Keep the final figure flow explicit: fit -> Wright map -> follow-up plots.
s <- summary(fit, profile = "facets", diagnostics = diag)
res <- s$results
s
plot(res, type = "wright", preset = "publication", show_ci = TRUE, top_n = Inf)

# Optional closest FACETS-style asterisk ruler (without mfrmr CI overlays).
plot(res, type = "wright", renderer = "facets",
     category_labels = rubric_labels, show_ci = FALSE,
     preset = "publication")

plot(
  res,
  type = "fit_pathway",
  fit_stat = "Infit",
  include_person = TRUE,
  top_n_person = 12,
  person_labels = "none",
  facet_labels = "flagged",
  preset = "publication"
)

This full-data figure caps the displayed person layer at 12 and suppresses routine point labels to keep the first screen legible. The selected person IDs and every retained facet row remain in plot(..., draw = FALSE)$data$table; use person_labels = "all" or facet_labels = "all" for a point-identification figure.

If you need residual-structure evidence for a final report, you can add residual PCA after the initial diagnostic pass. Treat this as an exploratory screen, not as a standalone unidimensionality test or as a DIMTEST/UNIDIM substitute. In MFRM reporting, a cautious claim should combine global residual fit, element-level fit, residual PCA, and local-dependence screens, for example: “evidence consistent with essential unidimensionality under the specified facet structure.”

diag_pca <- diagnose_mfrm(
  fit,
  residual_pca = "both",
  pca_max_factors = 6
)

summary(diag_pca)

Strict Diagnostics for RSM and PCM

For RSM and PCM, the package can now keep the legacy residual path and the strict marginal path side by side:

fit_rsm_strict <- fit_mfrm(
  data = toy,
  person = "Person",
  facets = c("Rater", "Criterion"),
  score = "Score",
  method = "MML",
  model = "RSM",
  quad_points = 7,
  maxit = 30
)
diag_rsm_strict <- diagnose_mfrm(
  fit_rsm_strict,
  diagnostic_mode = "both",
  residual_pca = "none"
)

fit_pcm_strict <- fit_mfrm(
  data = toy,
  person = "Person",
  facets = c("Rater", "Criterion"),
  score = "Score",
  method = "MML",
  model = "PCM",
  step_facet = "Criterion",
  quad_points = 7,
  maxit = 30
)
diag_pcm_strict <- diagnose_mfrm(
  fit_pcm_strict,
  diagnostic_mode = "both",
  residual_pca = "none"
)

summary(diag_rsm_strict)$diagnostic_basis[, c("DiagnosticPath", "Status", "Basis")]
summary(diag_pcm_strict)$diagnostic_basis[, c("DiagnosticPath", "Status", "Basis")]

When you want a compact simulation-based screening check for the strict branch, use evaluate_mfrm_diagnostic_screening() on a small design:

screen_rsm <- evaluate_mfrm_diagnostic_screening(
  design = list(person = 18, rater = 3, criterion = 3, assignment = 3),
  reps = 1,
  scenarios = c("well_specified", "local_dependence"),
  model = "RSM",
  maxit = 30,
  quad_points = 7,
  seed = 123
)
screen_pcm <- evaluate_mfrm_diagnostic_screening(
  design = list(person = 18, rater = 3, criterion = 3, assignment = 3),
  reps = 1,
  scenarios = c("well_specified", "step_structure_misspecification"),
  model = "PCM",
  maxit = 30,
  quad_points = 7,
  seed = 123
)

screen_rsm$performance_summary[, c("Scenario", "EvaluationUse", "LegacyAnyFlagRate", "StrictAnyFlagRate")]
screen_pcm$performance_summary[, c("Scenario", "EvaluationUse", "LegacySensitivityProxy", "StrictSensitivityProxy", "DeltaStrictMinusLegacyFlagRate")]

The same strict branch is now reflected in the reporting router:

chk_rsm_strict <- reporting_checklist(fit_rsm_strict, diagnostics = diag_rsm_strict)
subset(
  chk_rsm_strict$checklist,
  Section == "Visual Displays" &
    Item %in% c("QC / facet dashboard", "Strict marginal visuals", "Precision / information curves"),
  c("Item", "Available", "DraftReady", "NextAction")
)

Residual PCA and Reporting

pca <- analyze_residual_pca(diag_pca, mode = "both")
plot_residual_pca(pca, mode = "overall", plot_type = "scree")
data("mfrmr_example_bias", package = "mfrmr")
bias_df <- mfrmr_example_bias
fit_bias <- fit_mfrm(
  bias_df,
  person = "Person",
  facets = c("Rater", "Criterion"),
  score = "Score",
  method = "MML",
  model = "RSM",
  quad_points = 7
)
diag_bias <- diagnose_mfrm(fit_bias, residual_pca = "none")
bias <- estimate_bias(fit_bias, diag_bias, facet_a = "Rater", facet_b = "Criterion")
fixed <- build_fixed_reports(bias)
apa <- build_apa_outputs(fit_bias, diag_bias, bias_results = bias)

mfrm_threshold_profiles()
vis <- build_visual_summaries(fit_bias, diag_bias, threshold_profile = "standard")
vis$warning_map$residual_pca_overall

The same example_bias dataset also carries a Group variable so DIF-oriented examples can show a non-null pattern instead of a fully clean result. It can be loaded either with load_mfrmr_data("example_bias") or data("mfrmr_example_bias", package = "mfrmr").

Human-Readable Reporting API

spec <- specifications_report(fit, title = "Study run")
data_qc <- data_quality_report(
  fit,
  data = ej2021_study1,
  person = "Person",
  facets = c("Rater", "Criterion"),
  score = "Score"
)
iter <- estimation_iteration_report(fit, max_iter = 8)
subset_rep <- subset_connectivity_report(fit, diagnostics = diag)
facet_stats <- facet_statistics_report(fit, diagnostics = diag)
cat_structure <- category_structure_report(fit, diagnostics = diag)
cat_curves <- category_curves_report(fit, theta_points = 101)
bias_rep <- bias_interaction_report(bias, top_n = 20)
plot_bias_interaction(bias_rep, plot = "scatter")

Design Simulation and Prediction

The package also supports a separate simulation/prediction layer. The key distinction is:

if (requireNamespace("lme4", quietly = TRUE)) {
  gt <- mfrm_generalizability(fit)
  gt$coefficients[, c("G", "Phi", "GStatus", "PhiStatus",
                      "IdentificationStatus")]

  ds <- mfrm_d_study(
    gt,
    data.frame(Rater = c(2, 3, 4), Criterion = 4),
    residual_scaling = "sensitivity"
  )
  ds[, c("n_Rater", "n_Criterion", "ResidualScaling",
         "G", "Phi", "GStatus", "PhiStatus", "IdentificationStatus")]
}
sim_spec <- build_mfrm_sim_spec(
  n_person = 30,
  n_rater = 4,
  n_criterion = 4,
  raters_per_person = 2,
  assignment = "rotating"
)

recovery <- suppressWarnings(
  evaluate_mfrm_recovery(
    sim_spec = sim_spec,
    reps = 2,
    maxit = 30,
    include_diagnostics = TRUE,
    diagnostic_fit_df_method = "both",
    seed = 2
  )
)

summary(recovery)$recovery_summary[, c("ParameterType", "Facet", "RMSE", "Bias")]
plot(recovery, type = "summary", metric = "rmse", draw = FALSE)$data$plot_table

recovery_review <- assess_mfrm_recovery(
  recovery,
  min_reps = 2,
  min_se_available = NULL,
  max_mcse_rmse_ratio = NULL,
  max_rmse = c(facet = 1, step = 1, default = 1.5),
  max_abs_bias = c(default = 0.75)
)

summary(recovery_review)$checklist[, c("Section", "Item", "Status")]
summary(recovery_review)$reading_order
summary(recovery_review)$condition_reporting_notes
summary(recovery_review)$condition_review
summary(recovery_review)$diagnostic_reporting_notes
summary(recovery_review)$diagnostic_review

status_plot <- plot(recovery_review, type = "status", draw = FALSE)
status_plot$data$section_status
status_plot$data$reading_order

metric_plot <- plot(recovery_review, type = "metrics", metric = "rmse", draw = FALSE)
metric_plot$data$plot_table
metric_plot$data$guidance

recovery_bundle <- build_summary_table_bundle(
  recovery_review,
  appendix_preset = "recommended"
)
recovery_bundle$table_index[, c("Table", "Rows", "Role")]

pred_pop <- predict_mfrm_population(
  sim_spec = sim_spec,
  reps = 2,
  maxit = 30,
  seed = 1
)

summary(pred_pop)$forecast[, c("Facet", "MeanSeparation", "McseSeparation")]

keep_people <- unique(toy$Person)[1:18]
toy_mml <- suppressWarnings(
  fit_mfrm(
    toy[toy$Person %in% keep_people, , drop = FALSE],
    person = "Person",
    facets = c("Rater", "Criterion"),
    score = "Score",
    method = "MML",
    quad_points = 5,
    maxit = 30
  )
)

new_units <- data.frame(
  Person = c("NEW01", "NEW01"),
  Rater = unique(toy$Rater)[1],
  Criterion = unique(toy$Criterion)[1:2],
  Score = c(2, 3)
)

pred_units <- predict_mfrm_units(toy_mml, new_units, n_draws = 0)
pv_units <- sample_mfrm_plausible_values(toy_mml, new_units, n_draws = 2, seed = 1)

summary(pred_units)$estimates[, c("Person", "Estimate", "Lower", "Upper")]
summary(pv_units)$draw_summary[, c("Person", "Draws", "MeanValue")]

For a report or appendix handoff, pass the recovery objects through the same summary-table export route used by the rest of the package:

export_summary_appendix(
  list(recovery = recovery, recovery_review = recovery_review),
  output_dir = tempdir(),
  prefix = "mfrmr_recovery_appendix",
  preset = "recommended",
  include_html = FALSE,
  overwrite = TRUE
)

For an initial exploratory run, reps = 2 or another very small value is useful only to check the data-generating setup and refit path. For a study report, increase reps, keep the ADEMP-style metadata in the exported tables, and set substantive RMSE/Bias thresholds so that assess_mfrm_recovery() can mark those rows as ok, review, or concern rather than not_assessed.