mfrmr Linking and DFF

Limits that matter before using real data

• Whole-facet centering includes anchors. Unless a facet is noncenter_facet, centering_scope = "facet" makes all its levels sum to zero. With R01 = 0.70, the other ungrouped raters sum to -0.70; a single remaining rater is fixed at -0.70. Group totals also contribute. Incompatible totals stop before fitting, even when anchor_policy = "silent". centering_scope = "remainder" reproduces the earlier restriction that the ungrouped, unanchored levels sum to zero. Saved results are not converted by reloading them; fitting under a different convention requires re-estimation. For ordinary additive RSM/PCM JML, centering_scope = "anchor_origin" instead lets direct or group anchors replace their facet’s centering, as FACETS does. Choose an unanchored, active noncenter facet and repeat the same arithmetic mean target for each group member. This does not change the default or account for estimated-anchor error.
• The helper does not copy every setting. The model, estimation method, centering convention and MML integration mode are inherited automatically. Older baselines without convention metadata retain "remainder". Repeat the intended category policy, rubric, orientation, PCM/GPCM roles, uncentered facet and population model. RSM/PCM MML without a population model fixes ability to a standard-normal distribution; it does not estimate a changed cohort mean/variance. An explicitly requested population_formula = ~ 1 estimates normal-population moments and needs person_data; see help("fit_mfrm") for its assumptions.
• Person anchors require JML. Current MML integrates ability and rejects fixed Person values or Person group means. Common Persons can be present in the design without being anchored. Choosing JML changes the estimator and should follow the analysis question, not merely avoid an error.
• Anchor uncertainty is conditional. The new fit treats supplied values as exact. It does not propagate their source estimation error; an unavailable SE for a fixed value does not mean its source was perfectly known. Compare defensible anchor sets when conclusions depend on that choice.
• Locations are not the entire calibration. To hold thresholds and other supported calibration components fixed while scoring, use extract_mfrm_calibration(), validate and freeze the saved calibration, then call score_mfrm_calibration(). See help("mfrm_calibration_workflow") for its supported models and required checks.

The FACETS centering rules help distinguish an origin constraint from a fixed external value, but are not mfrmr’s partial-anchor convention. Uto (2021, pp. 1442–1444) explains common-rater/task linking and its design dependence. That simulation fixes generated common parameters and removes a redundant rater constraint; its Bayesian estimation does not validate mfrmr’s different constraints or supply uncertainty propagation for estimated anchors.

Identify the performance, judge and rating occasion

Keep one row per individual rating on one criterion. The following IDs illustrate identities to preserve; they are not supplied example data.

Also retain criterion, score, session, training phase and video/live mode. The same speaker giving a new speech needs a new performance ID when the measurement target is that speech. Preserve individual ratings before discussion or disclosure of reference scores. Consensus scores alone cannot recover differences among judges. Include performances spanning relevant levels and criterion profiles; there is no universal adequate recording count. The five-common-element helper warning is not such a rule.

Choose what is shared and what is allowed to change

In joint calibration, analyze years together, sharing an unknown measure for each common recording while allowing judges’ severities to vary by year. A judge-year effect represents a particular judge in a particular year; retain the original judge ID so these effects can later be compared. The recording measures can be estimated rather than supplied as numerical anchors.

In fixed-reference calibration, use previously estimated recording measures as fixed values for a later year and estimate that year’s judge effects. This preserves a chosen reference scale but conditions on its estimated values. When the recordings take the Person role, current mfrmr requires JML for direct Person anchors; MML rejects them. Sharing unknown performance measures is a different design and does not require these anchors. Formal JML structural intervals remain unqualified.

Returning judges should not automatically be fixed to their previous severities when their changes are the question. A panel mean fixed separately in every year also prevents estimating a common shift. For example, centering away a +0.40-logit increase for all judges would leave no panel-wide increase to report. FACETS describes group anchoring as an assumption about average severity, not proof that the panel stayed stable (Linacre, 2012, p. 9).

Specify scale constraints, criterion/category stability and ability-population assumptions together. Unrestricted year effects and judge-year effects would overlap; freely changing criteria can also absorb a common scoring shift. Separate fixed-standard-normal MML fits do not estimate changing cohort means. In a joint model, do not assign conflicting Person-level year covariates to a common recording rated in several years. Its rating occasion is not a different underlying performance.

Read the existing helper results for their stated purpose

When judges are the linking set, detect_anchor_drift() subtracts their common shift. If all selected differences are +0.40 logits, LinkOffset is +0.40 and every residual Drift is zero. This describes relative alignment, not evidence that the panel was stable. Interpreting LinkOffset itself as panel change is also unjustified without a separate scale reference.

Agreement after training, operational stability and correctness against an expert benchmark are different outcomes. Check familiar-score memory, for example with additional recordings not used to teach the reference scores. If the contest is rated live but the calibration uses video, paired live/video ratings of some performances can help assess whether the estimated severity transfers. Repeated ratings can be dependent even when their IDs are correct.

Annual-difference uncertainty needs the covariance between the relevant estimates, or propagation of reference-calibration uncertainty where applicable. Adding two squared SEs assumes away covariance; marginal interval overlap is not a test of change. Report the reference, assumptions, number of common performances and unavailable results alongside any point comparisons.

Reusing an analysis saved before the centering change

Start with what you want to preserve or change. Saving a new copy does not convert the statistical model.

For a runnable illustration, create an analysis with the legacy convention using the packaged ratings already loaded above. If you already have an old fit file, start by reading that file instead. The example fixes a Rater location, uses RSM with a fixed standard-normal population and keeps Person as the uncentered facet. These choices illustrate migration; choose the population assumption for your assessment question.

migration_anchor <- data.frame(Facet = "Rater", Level = "R01", Anchor = 0.70)
legacy_fit <- fit_mfrm(
  ratings, "Person", c("Rater", "Criterion"), "Score",
  model = "RSM", method = "MML", anchors = migration_anchor,
  noncenter_facet = "Person", centering_scope = "remainder",
  rating_min = 1, rating_max = 4, category_policy = "preserve",
  quad_points = 31, maxit = 400, reltol = 1e-9
)
fit_file <- tempfile(fileext = ".rds")
saveRDS(legacy_fit, fit_file)
restored_fit <- readRDS(fit_file)
unlink(fit_file)
spec <- specifications_report(restored_fit)$data_spec
spec[spec$Setting %in% c("NonCenteredFacet", "CenteringScope"), ]

CenteringScope = "remainder" describes the fitted constraint, including for native files saved before this field existed. It is not changed to the new default merely because a newer package reads the file. Do not assign a new value to fit$config$centering_scope or recenter its coefficient table: the likelihood, parameter mapping and uncertainty still belong to the original fit.

To change the convention, fit a separate model from the original responses. Here only the declared centering changes. Inspect its numerical and design status; the number of quadrature points is not an accuracy guarantee.

updated_fit <- fit_mfrm(
  ratings, "Person", c("Rater", "Criterion"), "Score",
  model = "RSM", method = "MML", anchors = migration_anchor,
  noncenter_facet = "Person", centering_scope = "facet",
  rating_min = 1, rating_max = 4, category_policy = "preserve",
  quad_points = 31, maxit = 400, reltol = 1e-9
)
summary(updated_fit, compute = "never", detail = "brief")
old_raters <- subset(restored_fit$facets$others, Facet == "Rater")
new_raters <- subset(updated_fit$facets$others, Facet == "Rater")
data.frame(
  Convention = c("remainder", "facet"),
  FixedR01 = c(old_raters$Estimate[old_raters$Level == "R01"],
               new_raters$Estimate[new_raters$Level == "R01"]),
  OtherRaterSum = c(sum(old_raters$Estimate[old_raters$Level != "R01"]),
                    sum(new_raters$Estimate[new_raters$Level != "R01"]))
)

Both fits keep R01 at 0.70. The remaining raters sum to zero in the legacy fit and -0.70 in the whole-facet fit. This changes R01’s comparison with the other raters; it is not just a renamed column. Retain the two analyses under different filenames, and report the convention with the other model settings. Complete facet anchors can also become incompatible with a whole-facet zero sum. Resolve that model specification explicitly rather than rounding anchors or silently switching the population model.

For more complex analyses, build_mfrm_replay_script(restored_fit, data_file = "ratings.csv") provides a recipe with stored anchors, group means and the original convention. It requires the original analysis data at the specified path. Review its generated fit_mfrm() call before deliberately changing the convention; preserve population covariates, category handling, orientation, uncentered facets and PCM/GPCM owners. This recipe refits. By contrast, the replay written by export_mfrm_results() reopens its saved results bundle without reconstructing the model.

If only a frozen calibration survives, retain its original scoring meaning. Generic old identification labels are insufficient to infer its centering convention. See vignette("mfrmr-portable-calibration", package = "mfrmr") for re-extraction and validation; changing a label is not recalibration. Earlier MML fits containing Person anchors require special review because those values were ignored. Current fitting rejects them; choosing centering_scope = "remainder" does not repair that earlier analysis.

Comparing an anchor with TAM or ConQuest

An anchor is a fixed value in all three programs, but the surrounding constraints determine what that value means. Copying a number from another program is insufficient when its parameterization or population model differs.

The TAM MML documentation, TAM JML documentation and ConQuest command reference describe these controls. For ConQuest, read import together with set; importing anchors alone does not specify all location and scale constraints. Parameter indices must be checked against the actual model design and retained levels. A changed index can fix the wrong coefficient without a syntax error.

For example, let four rater severities be \(\gamma_1,\ldots,\gamma_4\), with R01 fixed at \(\gamma_1=0.70\), and no group constraints:

The second row applies, for example, to TAM’s tam.mml.mfr(formulaA = ~ item + rater + step, constraint = "cases") and ConQuest’s model task+rater+step with lconstraints=cases. Those settings also leave the first task term uncentered, unlike mfrmr’s default. These are specific model configurations, not universal software conventions. The different rater contrasts show why applying the same anchor and then shifting the final locations is not generally an equivalent analysis. Group-mean anchors need corresponding linear constraints; fixing every member to the group mean would impose a different model.

For unit-slope RSM/PCM with all facet locations centered and an estimated population mean \(\mu\), another equivalent origin uses a zero population mean and uncentered task locations: set \(\theta'=\theta-\mu\) and \(b'_{task}=b_{task}-\mu\). Rater severities, centered step ladders and population variance stay the same, so every response probability is unchanged. This mapping applies to the configuration with Rater anchors above. If a task value is also anchored, its constraint must be transformed with it. In particular, estimating population variance and fixing it to one are different models; moving the origin does not reconcile that difference. Checking this mapping is distinct from running and comparing both programs’ optimizers.

To compare programs, first match the ratings, category coding, response probabilities, free and fixed parameters, and population assumptions. A coordinate transformation must carry the anchors and population parameters with it. Then check numerical integration and optimization. Finally compare SEs for the same parameter or contrast, transforming the covariance matrix when coordinates differ. A point-estimate match alone does not validate an uncertainty method or include source-anchor uncertainty.

A bounded check using all 48 Persons, four raters and four criteria in example_core illustrates both outcomes. Under the whole-rater constraint above, TAM 4.3.25 and ConQuest 5.47.5 agreed to within \(3.1\times10^{-6}\) in their free response-design coefficients. A separate explicit TAM design representing mfrmr’s legacy remainder constraints agreed with that mode to within \(1.8\times10^{-7}\), when mfrmr used 81-point Gauss–Hermite integration. Both comparisons fixed ability to \(N(0,1)\). For the latter comparison, 31-point integration in that legacy mode left a maximum coefficient difference of about 0.00067; increasing integration precision resolved it. Equal node counts across programs do not guarantee equal quadrature rules or accuracy.

This checks one RSM configuration and its point estimates, not general compatibility for PCM, GPCM, group constraints, JML or SEs. It does not make 81 points a universal requirement. More importantly, being able to reproduce mfrmr’s restricted model in another program does not establish that its partial-anchor constraints answer a particular linking question. Decide which locations and population moments should remain free before choosing the model.

Design the link, not only the percentage

No anchor percentage or count is universally adequate. In sparse designs, link quality also depends on coverage, the distribution and redundancy of the links, whether one articulation point or critical edge holds the network together, the location and model fit of the linking set, and the intended use of ranks or classifications. Inspect both Rater-centered and Person/task- centered assignment graphs rather than treating a minimum count as proof of a common scale. This distinction is consistent with work on sparse linking networks by Myford and Wolfe (2000), Wind and Jones (2018), and Uto (2021).

Fixed anchor values also carry uncertainty from their source calibration. Current mfrmr drift and chain summaries do not propagate every source-fit, offset, and cross-fit covariance component. When results depend on the anchor set, compare defensible alternatives and report that sensitivity. Robitzsch (2024, pp. 1181–1182, 1187–1188) examines a different source of linking error: random uniform DIF in fixed-item 2PL calibration. Its derivation neglects Person sampling error and does not provide a variance estimator for fitted population parameters. It supports checking invariance, not claiming that mfrmr propagates all linking uncertainty.

Review a particular school comparison before fitting

The question is whether the rating assignments identify the requested interaction contrast, rather than merely connect the panels. This small constructed assignment isolates that question: R4 rates School A only, while the other raters rate both origins. The scores only supply the declared three-category format; they are not evidence for a school effect, and eight performances are not a recommended sample size.

school_ratings <- expand.grid(
  Person = paste0("P", 1:8), Rater = paste0("R", 1:4),
  Criterion = c("C1", "C2"), stringsAsFactors = FALSE
)
school_ratings$School <- ifelse(
  school_ratings$Person %in% paste0("P", 1:4), "A", "B"
)
school_ratings$Score <- (as.integer(factor(school_ratings$Person)) +
  as.integer(factor(school_ratings$Rater)) +
  as.integer(factor(school_ratings$Criterion))) %% 3
school_ratings <- subset(school_ratings, !(Rater == "R4" & School == "B"))

school_questions <- list(
  facet = "Rater", group = "School", groups = c("A", "B"),
  model = "RSM", method = "JML",
  contrasts = rbind(
    R1_vs_panel = c(R1 = .75, R2 = -.25, R3 = -.25, R4 = -.25),
    R1_vs_R2 = c(R1 = 1, R2 = -1, R3 = 0, R4 = 0)
  )
)
school_review <- describe_mfrm_data(
  school_ratings, "Person", c("Rater", "Criterion", "School"), "Score",
  rating_min = 0, rating_max = 2, category_policy = "preserve",
  include_agreement = FALSE, comparisons = school_questions
)
school_review$comparison_review$table[, c(
  "Comparison", "Group1", "Group2", "State", "N_Group1", "N_Group2"
)]
school_review$comparison_review$cells

Each row applies its weights to the raters’ A-minus-B interaction effects. R1_vs_panel subtracts the equal-weight average across all four raters; R1_vs_R2 subtracts R2’s school difference from R1’s. The first is target_unidentified, while the second is target_identified_only. R4’s missing B ratings matter to the whole-panel reference. Dropping R4 and renormalizing weights would answer a different question.

The IdentificationBasis column makes the role of assumptions explicit. For the same ratings and target weights, compare a common normal ability population with a normal regression allowing school-specific means. These reviews do not estimate either model:

population_reviews <- lapply(
  list(common_normal = ~1, school_means = ~School),
  function(population) {
    question <- school_questions
    question$method <- "MML"
    question$population_formula <- population
    describe_mfrm_data(
      school_ratings, "Person", c("Rater", "Criterion", "School"), "Score",
      rating_min = 0, rating_max = 2, category_policy = "preserve",
      include_agreement = FALSE, comparisons = question
    )$comparison_review
  }
)
lapply(population_reviews, function(review) review$table[, c(
  "Comparison", "ConditionalState", "FreePersonState", "IdentificationBasis", "State"
)])

For R1_vs_panel, the common normal model reports population_restriction: the unrestricted-Person ambiguity disappears only after imposing that model’s population restrictions. Allowing school-specific means restores the ambiguity, so the basis is unidentified. No additional rating was observed in either review. The common population assumption is substantive; it is not a harmless change of coordinates or evidence that the missing R4/B ratings exist.

For R1_vs_R2, both models report free_person_response_model: the declared response map identifies this target even with unrestricted Person locations. That statement still assumes the response model, including its slope and step restrictions. It does not make an MML estimate distribution-free. Both positive MML conditional results keep State = "not_evaluated", because the review does not prove global marginal identification. Read IdentificationDetail and the review’s notes before selecting a model; a more convenient result is not a reason to assume equal ability populations.

The declared response model includes a zero-valued School main facet and a zero-margin Rater-by-School interaction, alongside free JML Person locations. Its settings are saved in school_review$comparison_review$fit_arguments. The review does not fit this rank-deficient model, estimate the identified contrast, or supply an interval. The ordinary full-fit guard still refuses the model. Identification of one target is not permission to set other effects to arbitrary values and report the resulting estimates.

Adding genuinely shared recordings from both origins can supply missing comparisons. Keep each recording’s original Person and School labels across raters. Neutral recordings may connect the measurement scale without filling the missing A/B interaction information. PCM reviews also require an explicit step owner. GPCM does not currently support this interaction model and is reported as unsupported; its separate slope/product qualifications cannot be substituted for this school review.

Preserve that question alongside a residual screen

For the bundled data used above, explicitly request an MML prospective Rater-by-Group review, then attach it to a residual screen of the same ratings. The model in this review is a prospective interaction model; the fitted screening model remains the earlier additive RSM. No interaction fit is silently added.

rater_levels <- sort(unique(as.character(bias_df$Rater)))
panel_weights <- matrix(-1 / length(rater_levels), 1, length(rater_levels),
  dimnames = list("First rater versus panel", rater_levels))
panel_weights[1, 1] <- panel_weights[1, 1] + 1
prospective <- describe_mfrm_data(
  bias_df, "Person", c("Rater", "Criterion", "Group"), "Score",
  rating_min = 1, rating_max = 4, category_policy = "preserve",
  include_agreement = FALSE,
  comparisons = list(facet = "Rater", group = "Group", groups = c("A", "B"),
    contrasts = panel_weights, model = "RSM", method = "MML")
)
rater_screen <- analyze_dff(
  fit, facet = "Rater", group = "Group", data = bias_df,
  comparison_review = prospective
)
summary(rater_screen)$comparison_review$table
support <- plot(rater_screen, draw = FALSE)
plot_data(support, "comparison_review")$weights
report <- dif_report(rater_screen)
report$comparison_review$table
saveRDS(rater_screen, file.path(tempdir(), "rater-screen-with-review.rds"))

Positive ConditionalState values concern the conditional linear response map. Under MML, State = "not_evaluated" deliberately withholds a positive global marginal-identification claim; it does not mean the residual calculation failed. A proven same-conditional-law ambiguity can still produce target_unidentified after mixing. FreePersonState and Basis expose the population assumption separately. Neither rank nor a finite residual difference qualifies a DFF test. The attached review keeps all requested rows, weights, counts and limitations; it does not replace the residual table with structural parameter estimates. Changed rating subsets or origin assignments require a new review before attachment.

Estimate the declared contrast in one joint calibration

When the assignment supports the whole model, the next question is the size of the specified interaction contrast in a joint fit. Use the bundled example_bias ratings from the screen above. Its Group labels illustrate the origin role; this synthetic example is not evidence about actual schools. Allow Group-specific population means through ~Group, so an equal-population restriction is not silently used to explain the group difference. School/Group remains a zero main facet; the population mean and rater-by-group interaction have separate roles.

The first comparison is the first rater’s A-minus-B interaction relative to the equal-weight whole panel. The second subtracts the second rater’s A-minus-B interaction. This compares school differences, not the raters’ overall severity.

joint_questions <- prospective$comparison_review$specification
joint_questions$population_formula <- ~Group
pair_weights <- setNames(rep(0, length(rater_levels)), rater_levels)
pair_weights[1:2] <- c(1, -1)
joint_questions$contrasts <- rbind(
  panel_weights, "First versus second rater" = pair_weights
)
joint_review <- describe_mfrm_data(
  bias_df, "Person", c("Rater", "Criterion", "Group"), "Score",
  rating_min = 1, rating_max = 4, category_policy = "preserve",
  include_agreement = FALSE, comparisons = joint_questions
)
joint_fit <- do.call(fit_mfrm, c(
  list(data = bias_df, quad_points = 63L),
  joint_review$comparison_review$fit_arguments
))
joint_results <- summary(joint_fit, comparison_review = joint_review)
joint_comparisons <- joint_results$interaction_comparisons
joint_comparisons[, c(
  "Comparison", "Group1", "Group2", "Estimate", "EstimateStatus",
  "State", "ConditionalState", "QuadratureState", "FormalInferenceEligible"
)]
attr(joint_comparisons, "fit_checks")$numerical
attr(joint_comparisons, "comparison_review")$weights
saveRDS(list(fit = joint_fit, review = joint_review, results = joint_results),
  file.path(tempdir(), "joint-school-comparisons.rds"))

The two fitted contrasts are approximately 0.58 and 0.92 logits. A positive value indicates the weighted relative shift toward higher coded scores for A versus B under this model. It is not an absolute advantage for School A, a comparison of population means, or proof of unfairness. The values differ from the earlier residual-screen statistic because that statistic asks a different question.

interaction_effect_table(joint_fit, comparison_review = joint_review) returns the same comparison table directly. Without comparison_review, it continues to return individual interaction cells. The numerical checks reevaluate the saved joint likelihood and full gradient without optimization. For this example, 31-node fixed quadrature did not meet the finer-grid check; 63 nodes were selected after checking agreement with 127. This order is not a universal accuracy recommendation. QuadratureState = "stable_at_retained_point" means the retained point passed the specified integration check, not that the model or its uncertainty has been validated.

EstimateStatus = "descriptive_only" keeps the fitted value separate from identification and inference. Positive conditional MML ranks still retain State = "not_evaluated"; this example does not prove global marginal identification. Inspect FitReadiness, BoundaryState, NumericalState, EstimateReason and the saved checks. A failed check leaves Estimate missing; OptimizerEstimate is then only a numerical trace. This extraction table supplies no SEs, confidence limits or formal tests. The separate interval example below uses mfrm_facet_intervals(comparison_review = ...) with additional joint-information and finer-grid checks. Keep the full saved object to preserve weights, empty cells, population assumptions and checks; a CSV does not retain these attributes.

Estimate an identified comparison when the full model is deficient

Suppose the last rater has no ratings from Group B. The first-versus-second rater comparison can still be identified, but the comparison with the whole panel includes an unobserved rater/group cell. An ordinary full fit correctly refuses this deficient design. Explicit comparison-only estimation can retain the observed likelihood and estimate its invariant targets without reporting arbitrary values for the unidentifiable interaction cells.

This example uses JML with free Person locations. It does not assume equal normal populations across groups. It uses the same declared interaction target as before, on the reduced assignment, but changes the estimator and the observed ratings; its numerical value is therefore not a replication of the preceding full-design MML value.

incomplete_ratings <- bias_df[!(bias_df$Rater == tail(rater_levels, 1) &
                                bias_df$Group == "B"), ]
partial_questions <- joint_questions
partial_questions$method <- "JML"
partial_questions$population_formula <- NULL
partial_review <- describe_mfrm_data(
  incomplete_ratings, "Person", c("Rater", "Criterion", "Group"), "Score",
  rating_min = 1, rating_max = 4, category_policy = "preserve",
  include_agreement = FALSE, comparisons = partial_questions
)
partial_fit <- do.call(fit_mfrm, c(
  list(data = incomplete_ratings, fit_scope = "comparisons",
       comparison_review = partial_review, interaction_policy = "silent"),
  partial_review$comparison_review$fit_arguments
))
partial_comparisons <- interaction_effect_table(partial_fit)
partial_comparisons[, c(
  "Comparison", "ConditionalState", "Estimate", "EstimateStatus",
  "BoundaryState", "NumericalState", "EstimateReason"
)]
partial_fit$checks$numerical
saveRDS(partial_fit, file.path(tempdir(), "partial-school-comparisons.rds"))

The whole-panel comparison remains missing; the first-versus-second rater contrast is approximately 1.08 logits, with target_identified_only and descriptive_only. This is an A-minus-B difference between those two raters’ interaction effects. It does not recover the missing last-rater/B effect. summary(partial_fit) returns the same table. This dedicated result cannot be used as an ordinary full model for Person scoring or individual-cell inference.

Identifiability alone is insufficient. For JML, the comparison route uses lpSolve to check whether any parameter direction can increase the likelihood without reaching a finite maximum. A detected or unchecked recession withholds all estimates, even if some contrasts might have limits in an extended model; that extension is outside this route. Failed convergence, insufficient reduced curvature or unsupported categories also retain missing values and reasons.

For MML, declare its population model in the review before fitting. The same coordinate reduction preserves every conditional probability and the marginal likelihood. Multiple starting variances and finer quadrature assess numerical adequacy; they do not prove global marginal identification or a unique global maximum. Positive conditional results still keep State = "not_evaluated". Neither route supplies SEs, confidence intervals or formal DFF tests. The saved representation and optimization are audit coordinates, not reportable cell parameters. See ?fit_mfrm for the supported scope and acceptance criteria.

Obtain CML uncertainty for the identified comparison

For the same incomplete assignment, CML conditions on each Person’s total score. It needs no ability distribution and avoids estimating one incidental ability parameter per Person. The inferential target remains the declared rater-by-origin interaction contrast. Start with the bundled synthetic ratings and specify the two comparisons explicitly:

cml_data <- load_mfrmr_data("example_bias")
cml_raters <- sort(unique(as.character(cml_data$Rater)))
cml_ratings <- cml_data[!(cml_data$Rater == tail(cml_raters, 1) &
                          cml_data$Group == "B"), ]
cml_weights <- rbind(
  panel = c(1, rep(0, length(cml_raters) - 1)) - 1 / length(cml_raters),
  pair = c(1, -1, rep(0, length(cml_raters) - 2))
)
colnames(cml_weights) <- cml_raters
cml_review <- describe_mfrm_data(
  cml_ratings, "Person", c("Rater", "Criterion", "Group"), "Score",
  rating_min = 1, rating_max = 4, category_policy = "preserve",
  include_agreement = FALSE,
  comparisons = list(facet = "Rater", group = "Group", groups = c("A", "B"),
                     contrasts = cml_weights, model = "RSM", method = "CML")
)
cml_fit <- do.call(fit_mfrm, c(
  list(data = cml_ratings, fit_scope = "comparisons",
       comparison_review = cml_review, interaction_policy = "silent"),
  cml_review$comparison_review$fit_arguments
))
cml_table <- interaction_effect_table(cml_fit)
cml_table[, c("Comparison", "ConditionalState", "Estimate", "SE", "Lower",
              "Upper", "FormalInferenceEligible", "EstimateReason")]
cml_fit$covariance
saveRDS(cml_fit, file.path(tempdir(), "cml-school-comparisons.rds"))

Read the pair row: the estimated contrast is approximately 0.925 logits, SE 0.518, with a pointwise 95% Wald interval of −0.090 to 1.939. It is the first rater’s A-minus-B interaction minus the second rater’s A-minus-B interaction. The interval includes zero. The panel estimate, SE and limits remain NA, because its weights involve the missing last-rater/B cell. Zero-filling that row would change the inferential claim. The full contrast covariance is retained in cml_fit$covariance, including missing entries for unavailable targets; summary(cml_fit) returns the same table.

Reuse that fit through the same interval interface used for full models. No estimation or information calculation is repeated. The missing panel comparison stays in the display as “not estimated”, with no point or interval:

cml_intervals <- mfrm_facet_intervals(cml_fit)
plot(cml_intervals)
apa_table(cml_intervals, which = "intervals")$table
saveRDS(cml_intervals, file.path(tempdir(), "cml-school-intervals.rds"))

Use mfrm_facet_intervals(cml_fit, level = .90) for pointwise 90% limits; only the normal critical value changes. The saved comparison_review retains the group order, weights and empty cells, and fit_checks retains the numerical qualification. apa_table(cml_intervals, which = "comparison_cells") shows which crossings were observed. These dedicated CML objects are not ordinary mfrm_results() fits and do not enable Person scoring or full-model diagnostics.

CML identification uses all possible patterns at the observed totals, not just the unconditional response design. Extreme-total Persons contribute no information. The fit removes unidentified directions only from its internal coordinates, then checks finite-MLE existence, convergence and the full conditional information in the identified subspace. A recession or failed check withholds inference with a reason. For PCM, declare model = "PCM" and its step_facet in the review before fitting.

The intervals are asymptotic under a correct response model, independent Persons, fixed structural dimension and increasing information. Numerical agreement with independent calculations is not a finite-sample coverage study for every school design. These estimates do not establish a calibrated DFF classification or recover Person abilities or between-group ability means. CML estimates no population SD, so it cannot supply a contrast divided by that fit’s estimated ability SD. The separate MML standardization uses its joint delta method, including uncertainty and cross-covariance for the fitted SD. Save the whole CML object as RDS to retain weights, checks and missing-cell information; CSV exports only the displayed columns.

Standardize an MML interaction contrast by that fit’s population SD

A contrast in logits and a contrast in ability-population SD units answer different questions. For the latter, define the denominator using the same fit that estimates the interaction. The following example uses the complete example_bias assignment and explicitly assumes one common normal ability population through population_formula = ~1. This differs from the preceding MML example with ~Group, which allowed group-specific means, and from the incomplete CML example. The change is a modeling assumption, not just a display option; do not replace ~Group with ~1 solely to obtain a standardized result.

ratings <- load_mfrmr_data("example_bias")
raters <- sort(unique(as.character(ratings$Rater)))
weights <- rbind(
  panel = c(1, rep(0, length(raters) - 1)) - 1 / length(raters),
  pair = c(1, -1, rep(0, length(raters) - 2))
)
colnames(weights) <- raters
sd_review <- describe_mfrm_data(
  ratings, "Person", c("Rater", "Criterion", "Group"), "Score",
  rating_min = 1, rating_max = 4, category_policy = "preserve",
  include_agreement = FALSE,
  comparisons = list(facet = "Rater", group = "Group", groups = c("A", "B"),
    contrasts = weights, model = "RSM", method = "MML", population_formula = ~1)
)
sd_fit <- do.call(fit_mfrm, c(
  list(data = ratings, quad_points = 63, interaction_policy = "silent"),
  sd_review$comparison_review$fit_arguments
))
sd_intervals <- mfrm_facet_intervals(
  sd_fit, comparison_review = sd_review, scale = "standardized"
)
sd_intervals$settings$population_sd
summary(sd_intervals)[, c("Target", "NativeEstimate", "Estimate", "SE", "Lower", "Upper")]
saveRDS(sd_intervals, file.path(tempdir(), "standardized-school-intervals.rds"))

Here the fitted SD is approximately 0.976. The first-minus-second rater’s A-minus-B interaction is 0.927 logits, or 0.950 population SD units, with SE 0.494 and a pointwise 95% interval of −0.018 to 1.918. The interval includes zero. NativeEstimate keeps the original scale alongside the standardized result, and covariance retains covariance between the two requested contrasts. The pointwise intervals do not adjust for selecting or examining multiple contrasts.

For contrast \(d\) and log population variance \(\lambda=\log(\sigma^2)\), \(t=d/\sigma\) has derivatives \(J_d/\sigma\) and \(-t/2\). The implemented variance is

\[ \operatorname{Var}(\hat t)\approx \frac{\operatorname{Var}(\hat d) +\hat d^2\operatorname{Var}(\hat\lambda)/4 -\hat d\operatorname{Cov}(\hat d,\hat\lambda)}{\hat\sigma^2}. \]

It uses the full joint marginal-information inverse. Dividing the native SE or interval endpoints by \(\hat\sigma\) leaves out denominator uncertainty and its covariance with the contrast. The ratio can be unstable when the fitted SD is near zero; standardization does not cure weak information.

plot(sd_intervals)

The usual saved-results path carries the same targets, scale, checks and review into plots, reports and table exports, without fitting again:

sd_results <- mfrm_results(
  sd_fit, intervals = list(school = sd_intervals),
  include = c("fit", "plots"), compute = "never"
)
sd_results$tables$facet_school_intervals[, c("Target", "Estimate", "ScaleReference")]
saveRDS(sd_results, file.path(tempdir(), "standardized-school-results.rds"))

For a supported full MML fit with ~Group, request scale = "native" instead; its fitted residual SD is not substituted for a marginal population SD. The fixed N(0,1) case uses known SD one. CML estimates no such SD, and this MML consumer does not extend inference to partial MML fits or GPCM interactions. Numerical eligibility means that source, information and finer-grid checks passed. CoverageQualified = FALSE and the retained marginal-identification state keep finite-sample performance and global identification separate from those checks. These intervals require the declared response/population model and independent Persons; they are not a calibrated DFF classification.

Comparing with MH or ordinal logistic DIF

mantel_haenszel_table() now supplies a separate binary MH comparison with explicit matching strata and one response per Person and target. Its default conditional exact inference uses stats::mantelhaen.test(); an asymptotic option is also available. analyze_dff() retains its residual and refit routes. ordinal_logistic_table() supplies the independent-Person cumulative-logit comparison through the optional ordinal package. Pooled repeated-rating comparisons remain planned work. First specify what each result describes:

First declare the intended response roster and rating plan before running the methods. Passing that roster through residual analyze_dff(data = ...) now selects the fitted observations by Person, all facets, score, weight and repeated-response multiplicity. It uses the original fitted expectations, without recalibrating on the selected subset. The same contract applies to dif_interaction_table().

Inspect $cases in each result: residual output maps InputRow in the supplied data to FitRow in fit$prep$data; MH/ordinal output identifies input rows with Row. Missing groups, invalid responses and responses absent from the fit cannot silently acquire another observation’s residual. Requesting more copies than the fit contains is rejected. Identical fitted copies can be exchanged for this descriptive calculation only when their model moments agree; their FitRow assignment does not recover an occasion identity.

Residual N counts response rows, while NPersons counts distinct Persons. For example, 12 Persons with three criteria give 36 rows for a rater, not 36 independent Persons. Descriptive residual summaries can retain those rows; ordinary MH/ordinal inference cannot treat them as independent target responses. Earlier saved residual summaries must be recomputed from the saved fit and source data to obtain this matching audit; a model refit is unnecessary.

Method-specific exclusions can still differ: MH removes uninformative strata, whereas ordinal models may retain their observations. Report these differences and the matching rule; do not automatically intersect outcome-dependent inclusion masks and treat that as a prespecified common cohort. Matching rows is necessary for understanding a comparison, but does not make the methods’ quantities or null hypotheses identical.

The MH result keeps the MH point estimate separately from the conditional maximum-likelihood estimate used by exact inference. Its strata and cases tables retain group/category counts, opposing effects, missing responses and strata without overlap. Supply multiple matching columns, such as total and rating plan, through strata; continuous values are not silently binned. The function rejects duplicate Person-target cells. Pointwise intervals are not simultaneous across targets; Holm adjustment applies to p-values only. Exact probability-ordered two-sided p-values need not invert equal-tail confidence intervals. See ?mantel_haenszel_table for a runnable example.

The ordinal result fits matching-only (M1), matching plus group (M2), and matching plus group and their interaction (M3) to identical complete rows. Its table separates uniform (M2/M1), nonuniform (M3/M2) and joint (M3/M1) likelihood-ratio tests, with Holm adjustment across all three tests and targets. DeltaNagelkerkeR2 names the pseudo-R-squared definition explicitly. The contrasts table reports M3’s reference/focal upper-tail odds ratios at the endpoints and midpoint of the overlapping matching ranges; values above one favor higher reference scores. Fits use all complete responses. An interaction means this odds ratio changes with matching, so one constant group coefficient does not summarize the full difference.

The separate proportional_odds table tests whether slopes can differ across category cuts, using a global nonproportional alternative. This diagnostic has its own multiplicity family and does not select which group tests are run. For binary responses it is not applicable. models, matching and cases preserve convergence information, group ranges and row exclusions. Empty declared categories, absent range overlap and numerical failures remain unavailable; repeated rows are rejected. See ?ordinal_logistic_table for a runnable example and the exact numerical checks. All reported p-values are asymptotic: convergence or a nonsignificant diagnostic does not establish adequate small-sample inference or proportional odds.

A bounded known-model check illustrates this limit. With three categories, balanced groups, a fixed matching grid and total N=40, all 1,000 simulated data sets had computable group tests. Under no group difference and true proportional odds, the unadjusted joint test rejected in 65/1,000 cases (6.5%; 95% Monte Carlo Wilson interval 5.1%–8.2%). The proportional-odds diagnostic was unavailable in 42 cases and rejected in 114/958 available cases (11.9%; interval 10.0%–14.1%), at nominal 5%. Treating its numerical availability as assurance of a valid 5% test would therefore be misleading. For the default Holm-adjusted family of three group tests, at least one rejection occurred in 52/1,000 cases (5.2%; interval 4.0%–6.8%). This family is separate from the proportional-odds diagnostic and does not repair it. At N=200 and N=1,000 under the same mechanism, joint rejection was 4.7% and 4.6%, and diagnostic rejection was 5.8% and 4.5%, with all results available. These particular conditions do not establish a universal sample-size rule. The function offers asymptotic comparisons, not qualified small-sample tests. No significance cutoff or numerical tolerance was changed in response to these results; the diagnostic does not determine whether group comparisons are run.

Matching is a substantive choice. In the binary Rasch case with a common item set and invariant other items, conditioning on the total including the studied item removes the group ability distribution from that item’s conditional response. Removing the studied item can leave group differences even without IRT DIF (Holland and Thayer, 1986, pp. 14–18). This is not a general endorsement of any total score. Contamination by other DIF items, missing responses and different rater/task plans require separate attention. In a many-facet assignment, two persons with the same total may have faced different raters or tasks. Specify common exposure or an explicit plan adjustment with overlap before treating their totals as comparable.

For each target, retain the persons, group coding, category order, matching definition, exclusions and rating plan. Repeated ratings from a Person are not independent new persons in ordinary MH tables or logistic likelihoods. A comparison must account for that dependence; silently averaging ratings changes the response being analyzed. Connectedness alone does not resolve either problem.

MH can miss effects that reverse direction across matching strata. Rogers and Swaminathan’s study also used different purification procedures for MH and logistic regression, so its comparison is not solely a change of test statistic (1993, pp. 111–115). Report uniform and nonuniform terms separately when fitting a logistic comparison. In Zumbo’s three-model sequence, add the group term to the matching-only model, then add the matching-by-group interaction; the joint test compares the last model with the first (1999, pp. 23–27).

The cumulative proportional-odds assumption in that ordinal model is different from the adjacent-category structure of RSM/PCM and is not generally satisfied by it. Check the matching-response relationship and model adequacy before interpreting a difference between procedures as DFF. Pseudo-R-squared definitions and their changes must be named: the historical .130 screen in Zumbo’s application is not a universal cutoff or a default for mfrmr. A nonsignificant result is not evidence that differential functioning is absent, and statistical DIF alone does not establish bias.

Zumbo, B. D. (1999). A handbook on the theory and methods of differential item functioning (DIF): Logistic regression modeling as a unitary framework for binary and Likert-type (ordinal) item scores. Ottawa, ON: Directorate of Human Resources Research and Evaluation, Department of National Defense.

Compare an explicitly shared response roster

The following self-contained example puts the three procedures beside each other. It requires the optional ordinal package. The binary Rasch design has 240 independent Persons answering the same eight items once each. Both groups have the same item difficulties, but different ability distributions. These synthetic data illustrate how to assemble and audit a comparison; this one sample is not a validation study or a recommended sample size.

Here Total includes the studied item and is computed only after checking the complete common item roster. MH conditions on each total separately; the logistic models instead assume a linear relationship with total on the logit scale. Even in this Rasch example, those are different matching models. The residual calculation uses the fitted Person and item parameters.

library(mfrmr)
set.seed(241010)
group_order <- c("Reference", "Focal")
person_group <- rep(group_order, each = 120)
ability <- rnorm(240, mean = ifelse(person_group == "Focal", 0.5, -0.5))
shared <- expand.grid(Person = seq_len(240), Item = paste0("I", 1:8),
                      KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE)
shared$Group <- person_group[shared$Person]
difficulty <- setNames(seq(-1, 1, length.out = 8), paste0("I", 1:8))
shared$Score <- rbinom(nrow(shared), 1,
  plogis(ability[shared$Person] - difficulty[shared$Item]))
stopifnot(!anyNA(shared), all(table(shared$Person, shared$Item) == 1L))
shared$Total <- ave(shared$Score, shared$Person, FUN = sum)

shared_fit <- fit_mfrm(shared, person = "Person", facets = "Item",
                       score = "Score", method = "JML", model = "RSM")
residual <- analyze_dff(shared_fit, facet = "Item", group = "Group",
  data = shared, focal = group_order[2], method = "residual")
mh <- mantel_haenszel_table(shared, person = "Person", target = "Item",
  score = "Score", group = "Group", strata = "Total", groups = group_order)
ologit <- ordinal_logistic_table(shared, person = "Person", target = "Item",
  score = "Score", group = "Group", matching = "Total",
  groups = group_order, categories = 0:1)

First retain each method’s response audit. The join below uses original row numbers from the same shared data, not the position of a row in a result table. It retains exclusions rather than choosing the intersection of the methods’ retained rows. With one response per Person-item, the counts below are also distinct-Person counts within each item and group.

case_audit <- rbind(
  data.frame(Method = "Residual", Row = residual$cases$InputRow,
    Included = residual$cases$Included, Reason = residual$cases$Reason),
  data.frame(Method = "MH", mh$cases[c("Row", "Included", "Reason")]),
  data.frame(Method = "Logistic", ologit$cases[c("Row", "Included", "Reason")])
)
case_audit$Target <- shared$Item[case_audit$Row]
case_audit$Person <- shared$Person[case_audit$Row]
case_audit$Group <- shared$Group[case_audit$Row]
case_audit$Count <- 1L
case_counts <- aggregate(Count ~ Target + Method + Group + Included + Reason,
                         data = case_audit, FUN = sum)
case_counts[case_counts$Target == "I2", ]

Next align effect direction explicitly. focal above makes the residual contrast Reference minus Focal. All odds ratios below are Reference/Focal odds of a higher score. Thus a positive residual difference and an odds ratio above one point in the same score direction, but their magnitudes cannot be compared. The MH column uses the conditional maximum-likelihood estimate from the default exact analysis; the distinct MH estimator remains in mh$table$MHOddsRatio. Logistic values come from M3 and can change direction over total scores. Keep their matching values and numerical status visible.

r <- residual$dif_table
stopifnot(all(r$Group1 == group_order[1]), all(r$Group2 == group_order[2]))
effects <- data.frame(Target = r$Level, ResidualDifference = r$Contrast,
  ResidualStatus = ifelse(is.finite(r$Contrast), "descriptive", "unavailable"))
m <- mh$table[c("Target", "OddsRatio", "Status")]
names(m)[2:3] <- c("MH_ExactOR", "MH_Status")
effects <- merge(effects, m, by = "Target", all = TRUE)
mid <- subset(ologit$contrasts, Position == "OverlapMid",
              select = c(Target, Matching, OddsRatio, Status))
names(mid)[2:4] <- c("TotalAtMidpoint", "Logistic_MidOR", "Logistic_Status")
effects <- merge(effects, mid, by = "Target", all = TRUE)
effects[effects$Target == "I2", ]

# Preserve endpoint effects; a midpoint alone can conceal a crossing effect.
subset(ologit$contrasts, Target == "I2")
subset(mh$table, Target == "I2", select = c(Target, NReference, NFocal,
  Lower, Upper, PValue, PAdjusted, Status, IntervalStatus))
subset(ologit$table, Target == "I2", select = c(Target, Comparison,
  DeltaNagelkerkeR2, PValue, PAdjusted, Status))

In the displayed sample, residual and logistic calculations use 120 Persons per group for I2. MH uses 96 Reference and 97 Focal Persons; its audit accounts for the other 47 responses in constant-response strata. This is a difference in usable conditional information, not an unnoticed change of the input roster. At total 0 or 8, the item response is determined by the total itself. Logistic endpoint odds ratios are predictions from its smooth response model; they are not empirical odds ratios in those degenerate strata. Overlapping matching ranges alone do not establish response support or a correct linear logit relationship. A binary proportional-odds diagnostic cannot check that relationship.

All eight items were analyzed before displaying I2, chosen in advance as the example target. MH adjusts across eight item tests; logistic adjusts across 24 item-by-comparison tests. These are separate families, not an adjustment for selecting among methods. Residual differences have no p-values. MH’s intervals are pointwise; logistic contrasts here have no intervals. For these binary data the proportional-odds diagnostic is not_applicable_binary, not evidence that proportional odds holds for a future ordinal analysis.

Read the case counts before interpreting different results: MH can omit constant-response or single-group strata that logistic retains. Logistic Included = TRUE means a row was supplied to a model, not that estimation succeeded; inspect the effect/test statuses and ologit$models. A missing effect stays NA, with its reason, rather than being replaced by zero or removed from the comparison. The descriptive JML residuals do not acquire inferential validity from the other methods’ p-values.

To retain this analysis, save the input, matching definition, calibration and complete results together. Later displays can be rebuilt from this object without rerunning the models.

comparison <- list(data = shared, groups = group_order, categories = 0:1,
  matching_rule = "Complete common eight-item total, including the target",
  fit = shared_fit, residual = residual, mh = mh, ordinal = ologit,
  effects = effects, case_audit = case_audit, case_counts = case_counts,
  session = sessionInfo())
saveRDS(comparison, "dff-comparison.rds")

For actual ordinal scores, retain the declared categories and omit ordinary binary MH; do not dichotomize merely to fill this table. For a Rater target with several Criteria per Person, the repeated rows do not satisfy this example’s independence contract. Differing rater/task plans also require justified adjustment. Those extensions remain open. Agreement of signs or nonsignificant tests in this display is not a robustness, invariance or fairness classification.