GPCM scope and current limitations

Does an unavailable interval mean interval estimation is impossible?

No. Distinguish the target, the method, whether this calculation returned bounds, and what supports their interpretation. A rater location, a discrimination and a new Person’s ability need different uncertainty calculations. Changing a covariance estimate does not change the target or remove estimation bias.

These distinctions also apply when an output has no CoverageQualified field: its absence does not establish coverage. A legitimately unbounded interval is different from an unresolved search. Methods with justified approximations can be useful without exact finite-sample guarantees, but the approximation and known adverse evidence must remain visible.

The interval guide locates the public routes without fitting a model:

interval_routes <- mfrmr::mfrmr_interval_guide("gpcm")
interval_routes[interval_routes$Route %in% c(
  "Fixed-facet model/sandwich intervals", "GPCM relative-slope intervals",
  "GMFRM component-slope intervals", "GPCM slope profile intervals",
  "Corrected-JML local root uncertainty"
), c("Route", "PrimaryHelper", "GPCMStatus")]

In particular, one-family MML location intervals use mfrm_facet_intervals(); they are separate from confint(fit, parm = "slopes"). The location helper supplies pointwise Wald intervals using the full joint marginal information, including estimated slopes, steps and population parameters. Their asymptotic justification requires independent Persons, a correctly specified, identified model with an interior true parameter, a fixed number of structural parameters, nonsingular limiting information and accurate integration and optimization. Adding Persons must add information in the relevant directions. The same-fit population-SD-standardized contrasts use the full delta method, including SD uncertainty and cross-covariance. These are supported approximations within the helper’s declared scope; numerical checks do not guarantee small-sample coverage or robustness to a misspecified ability distribution. GPCM sandwich and JML structural intervals remain unavailable, and the optional slope profiles and two-family location extensions retain their separate experimental status. The research existence of an alternative method does not make it an implemented public API. Conversely, an unsupported API route is not a mathematical impossibility result for every method.

Score a new cohort using the corrected calibration

Suppose next year’s learners use the same rubric and the same calibrated raters. predict_mfrm_units(fit, new_ratings) uses only those supplied ratings and the corrected facet, step and slope values. It adds a normal reference prior to obtain EAP scores. This is a different question from estimating the original learners’ profile maxima: even all-minimum or all-maximum responses can have finite EAPs. JML did not estimate this prior. Consider a substantive alternative with scoring_prior = list(mean = ..., sd = ...); sensitivity to that choice matters particularly for short or extreme response patterns.

To score in another R session without retaining the training data:

# draft <- extract_mfrm_calibration(fit, scoring_quad_points = 101)
# calibration <- freeze_mfrm_calibration(validate_mfrm_calibration(draft))
# save_mfrm_calibration(calibration, "rubric-calibration.rds")
# calibration <- load_mfrm_calibration("rubric-calibration.rds")
# scores <- score_mfrm_calibration(calibration, new_ratings)
# summary(scores)
# plot(scores, type = "interval")
# plot(scores, type = "precision")
# as_ggplot(scores)

Here 101 is an example integration order, not a universal accuracy setting. Each batch is checked against higher-order adaptive integration. Increase the extraction order and recreate the artifact if those checks fail. When a Person-by-facet cell legitimately contains independent repeated events, supply a column of distinct event IDs through event_id; duplicate cells without event IDs are refused. Missing scores cause an error by default; missing_response = "omit" records omissions and leaves missing-only Persons unscored. It does not fill unassigned design cells.

Extraction checks the corrected equation and its full Jacobian at the saved solution; it does not refit the model or demand an ordinary JML likelihood maximum. A covariance failure does not by itself prevent scoring. File format 5 preserves the correction order, source checks, scale and prior while omitting original ratings and Person estimates. The posterior SD and intervals hold calibration fixed: they do not propagate RootSE, resolve residual calibration bias or provide corrected Person ML/WLE. Calibration uncertainty and the suitability of the prior for the new cohort remain separate issues. Saved score objects support their own summaries and plots; they are not attachments to the corrected fit’s mfrm_results() report.

Describe the observed ratings at the corrected estimates

Use mfrm_response_diagnostics() to ask how the observed scores differ from the scores expected at the corrected calibration and each learner’s reprofiled ability. For example, select Judge groups to inspect rating residuals without automatically labelling anyone an inaccurate or overly predictable judge:

# response <- mfrm_response_diagnostics(fit, group_by = "Judge")
# response$measures
# response$rows
# response$probabilities
# plot(response, style = "paired")
# as_ggplot(plot(response, style = "scatter", draw = FALSE))
# res <- mfrm_results(fit, response_diagnostics = response, compute = "never")
# mfrm_report(res)

These are conditional fitted probabilities. The calculation holds the corrected facet/step/slope estimates and reprofiled abilities fixed. It does not average over a posterior distribution or account for uncertainty in those estimates. A comparison between observed and expected scores uses the same data as fitting; it is not a test of performance on unseen ratings. Sparse assignments and repeated observed events retain their original row identities. Selecting rows changes the summarized events but does not re-estimate ability.

All-minimum or all-maximum learners have infinite profile limits. Their conditional distribution concentrates on one category, so its variance is zero. The probabilities, expected score and raw residual remain available, but a residual divided by its SD is undefined at these saved values. No automatic convention replaces 0/0 with zero. Consequently:

• Infit can still be returned for a group: its sum of squared raw residuals divided by the sum of conditional variances is defined when that denominator is positive and all required values are available. Zero-variance rows remain in those sums.
• Outfit is missing if any selected standardized residual is undefined. The result records partially_available and explains why. Available counts standardizable rows; it is not an Infit exclusion count.
• A group consisting entirely of zero-variance rows has neither index. Rows are never dropped silently. An explicit rows selection changes the target group and must be reported.

The paired plot retains an available Infit when Outfit is missing; the scatter plot requires both. Both retain complete tables and reasons in plot_data(). There are no reference cutoffs, automatic overfit warnings, p-values or rater-quality decisions. Saved results and reports preserve the conditional definition and check the source slopes, Person profiles and exact observed roster before accepting an attachment. This does not enable ordinary diagnose_mfrm() tests or corrected-JML model comparisons.

The adjustment uses the MLE plug-in profile-score construction discussed by Dhaene and Weidner (2023), section 8.1. Exact owner-total expectations and a full equation Jacobian are used within the supported state limit (5,000 joint states per assignment pattern). This numerical implementation does not establish general residual-bias or sampling-coverage guarantees. See ?fit_mfrm for the full input contract.

What the paired automatic-integration study supports

The package implements the two-family response equation of Uto and Ueno (2020, equation 9) with explicit scale constraints, marginal maximum likelihood, numerical diagnostics and conditional scoring. This is a GMFRM implementation; it is not a new proposal of that response model or a reproduction of its Bayesian estimator. Numerical admission of an experimental interval and demonstrated nominal coverage are different claims.

A paired study evaluated 1,200 datasets with fixed-grid EM and adaptive direct MML: 2,400 fits of 1,200 simulated datasets. The 12 conditions crossed 60, 240 or 960 Persons, two rating assignments and generating ability SD .5 or 1, with 100 replications each. The four design/SD conditions at a given N share random draws within each replication, so the entire set is not independent. Each condition has 100 independent replications; pooled comparisons account for the shared replication blocks. This original procedure preceded the per-Person optimizer review described below; later replays are reported separately. There were three tasks, six raters and three score categories. Each source had 6N ratings. In the common-Person design, 20% of Persons were rated by all six raters and the remainder by one, on all three tasks. In the rotating-pair design, every Person was rated by two raters on all tasks. The two designs therefore had the same total rating count and different overlap. Each rater saw N/3 distinct Persons: 20, 80 or 320. The 6N rating rows are not 6N independent Persons.

Both generating populations were normal, with truth transformed to the fitted N(0,1) scale. This is not a test of robustness to a wrong population shape. The study used no anchors, weights, differential-functioning interactions, unobserved assigned scores or extra within-task response dependence. Both engines began at 31 points with a 241-point ceiling including reference rules. Direct used gpcm_mml_start = "neutral_em"; its costs include both adaptive starts and EM initialization. All attempts, including refusals, remain in the delivery denominator.

Admission means the automatic numerical checks passed, including those needed for the component-slope intervals. It does not establish a global maximum, design-wide identification or calibrated decisions. At N=60, direct succeeds where EM refuses in 68 sources, while EM succeeds where direct refuses in 44; neither engine dominates on every source. The paired delivery difference is 6 percentage points, with Monte Carlo SE about 2.41 points after averaging the paired difference across the four conditions within each of the 100 independent replication blocks.

These elapsed times describe the original full fit/check/component-interval workflow on one Mac with ten concurrent single-thread workers. They include failures and are not isolated engine benchmarks or forecasts for another computer. Later exact calculation reuse does not retroactively change this timing record. At N=960, a direct attempt took as long as about 32 minutes. max_points does not cap that duration or optimizer iterations.

Profiling those saved runs attributes about 75–77% of the fixed-EM stage cost at N=240/960 to information and precision checks. A subsequent exact probability-reuse change shortened one N=60 refusal workflow from 3.17 to 2.45 seconds; corresponding successful N=240/960 workflows changed from 10.93 to 10.59 and 65.89 to 65.19 seconds. These are means of two sequential paired repetitions per saved input, including interval computation where available, with identical output objects. They illustrate the limited scope of that optimization; small timing differences are noisy and do not establish a package-wide speedup. The numerical admission criteria were not relaxed.

After a successful automatic two-family fit, subsequent slope and location intervals reuse its complete numerical qualification in the same R session. Only the latest qualified automatic fit is retained, and the full fitted object, mfrmr.max_information_bytes and mfrmr.use_cpp11_backend settings must be unchanged. Retained source and qualification together must fit within that memory budget. Changed inputs or settings, warnings, refusals and a new R session require fresh checks. This transient reuse does not alter saved model objects or suppress the intervals’ experimental warnings. It avoids repeating completed checks; it does not make a previously refused fit eligible.

A separate paired check of this qualification reuse, using two sequential single-thread repetitions per saved input, changed the mean full workflow time from 10.75 to 7.28 seconds at N=240 and from 67.59 to 47.17 seconds at N=960. Complete fit and interval objects were identical. The N=60 refusal remained a refusal, with times of 2.50 and 2.57 seconds. These selected examples show savings from eliminating a repeated successful check, not a general runtime guarantee or improved inferential performance.

The component-slope intervals were pointwise nominal 95% log-Wald intervals. In the N=240 rotating-pair, SD=.5 condition, Task t3 contained truth in 89/100 sources for both engines (Monte Carlo 95% interval approximately 81.2–94.4%). This illustrates a coverage limitation; it is not a multiplicity-adjusted test across all correlated targets. With only 100 replications per condition, a true 95% coverage rate has Monte Carlo SE about 2.2 percentage points. The study cannot certify coverage within one percentage point of nominal. It evaluates component slopes, not all location contrasts, Person-score intervals, rater flags or decision thresholds.

Parameter bias also needs its denominator. For example, the N=60 rotating-pair, SD=.5 direct fits had mean log-slope error -0.435 for rater r3 (Monte Carlo SE .128), among only 54 admitted sources. This is conditional on admission; the other 46 do not become accurate estimates or zero errors. The result does not establish a general bias correction. A larger N improves delivery in this grid, but does not by itself qualify every slope or feedback decision.

Saved-input replay after aligning optimizer and interval checks

The revised direct procedure was applied to all 400 retained N=60 datasets, using the original starting policy, BFGS, maxit = 500, reltol = 1e-10, 31 initial points and a 241-point ceiling. These inputs had already been inspected during diagnosis: this is a matched replay, not an independent held-out confirmation study. The other 2,000 original study fits were not repeated. At these settings, the gradient threshold is unchanged at N=240/960; one retained control at each of those sample sizes reproduced the parameters, objective and complete intervals exactly.

Overall delivery rose from 277/400 (69.25%) to 337/400 (84.25%): 61 newly admitted sources and one newly refused source. The paired gain is 15 percentage points, with Monte Carlo SE 1.81 points from the 100 replication block means. All 61 gains came from the original stationarity refusals. Same-rule likelihoods did not worsen in any of 449 shared-stage comparisons. The revised procedure still refuses 63 sources: 23 at direct convergence, 33 at joint information, five at the integration ceiling and two at the stationarity check. These are stopping locations, not a proof of a unique underlying cause. The weak- information example below explains the newly refused source.

Improved delivery does not establish nominal coverage. In the rotating-pair, SD=1 condition, rater r1’s interval contained truth in 81 of 95 admitted sources: 85.3%, with an exact Monte Carlo interval of approximately 76.5–91.7%. The original result was 69/78. The 78 sources admitted by both procedures retained 69 inclusions; the 17 newly admitted sources added 12. Thus the fraction of all 100 attempts that both returned an interval and included truth increased from 69% to 81%, while coverage conditional on receiving an interval fell. This is a descriptive target-specific result, without multiplicity adjustment across the examined components. The component intervals remain experimental and coverage-unqualified. In the rotating-pair SD=.5 condition, rater r3’s conditional mean log-slope error also remained substantial: about -.422 among 69 admitted sources, compared with -.435 among the original 54.

The replay used eight concurrent single-thread workers; its elapsed times cannot be interpreted as a controlled speed comparison with the original study’s ten workers. Extra polishing can add computation. Its purpose here is a more adequately finished estimate under the existing precision criteria.

Planning a practical two-family analysis

There is no universal minimum N. The N=240 design is a useful reference when planning a comparable calibration study, not a validated minimum or a claim that N=240 is sufficient for a different rubric or rater roster. N=60 had substantial delivery losses even in correctly specified simulated models. Increasing calibration N can improve task/rater estimates; it does not remove the uncertainty in a Person scored on only a few performances.

Before collecting ratings, specify which quantity must be precise: a Person score, a task/rater slope, a severity contrast or a classification decision. Then specify an acceptable interval width or decision error and examine the proposed assignment, rather than applying a Persons-per-parameter rule. Check category use by rater, overlap across Persons and tasks, and whether the model’s conditional-independence assumption is plausible when several raters score the same performance. A connected graph is necessary for the component decomposition but does not prove marginal identification. Sparse categories or weakly linked rater/task blocks require design changes or a scientifically defensible simpler model; deleting failed fits is not a remedy.

New two-family adaptive direct fits use the smaller of the ordinary total-gradient review tolerance and N * 1e-6. This can activate the existing bounded polishing stages when reltol <= 1e-9 (the default or stricter) and an optimizer returns code zero before the per-Person score is small enough. With a looser requested tolerance, the existing automatic-polishing restriction still applies; a failed gradient review remains a refusal. The change does not increase the requested iteration ceiling for any stage or change the starting-value policy. Each accepted polishing stage must preserve or improve the objective within roundoff. The fitted controls record this rule; older saved models retain their original convergence review. Passing the optimizer review still does not establish positive information, adequate integration, global optimality or interval coverage. A repaired score failure may reveal a later integration failure.

One retained N=60 rotating-pair example also shows why a better likelihood need not make an interval available. The improved start reduced the negative log likelihood by about 5.8e-6, while one rater slope moved from .0032 to .0027. The information matrix remained positive definite, but its smallest eigenvalue fell from about 1.26e-6 to 4.52e-7. The ill-conditioned-information safeguard then rejected a standardized Newton displacement of .000188 against its .0001 limit, at both 31 and 61 points. This is evidence of a weakly curved, insufficiently finished solution under the current numerical rule, not proof that no valid interval method can exist. Preserve the likelihood-selected solution and its refusal reason. Rater overlap, category support and model sensitivity warrant review; increasing the integration ceiling alone does not address this example.

Three retained EM refined-information failures illustrate the distinction. At the same 31-point fitted parameters, the 61-point Hessian had a negative eigenvalue in all three cases. Refitting from 61 points passed that earlier obstacle and advanced to 121 points, but the 121/241 comparisons still failed precision checks. Thus integration can contribute to an information refusal. These selected probes do not establish the cause of every refusal or justify automatically retrying every failed check.

Rater severity, rater discrimination and differential rater bias are different quantities. A low discrimination estimate is not by itself evidence of unfair scoring, nor a high severity estimate a calibrated misconduct flag. The current two-family fit excludes bias interactions and ordinary estimate_bias() tests. Use its supported descriptive response diagnostics and retain all uncertainty qualifications. A feedback process requiring validated bias screening needs an appropriate supported model and a separately justified screening rule.

The current two-family estimator also excludes fixed Person, task and rater anchors. Common Persons create overlap; they are not fixed ability anchors. Separately standardizing two cohorts to N(0,1) does not establish a common operational scale. Reusing a saved calibration to score supported new Persons is a separate conditional workflow and does not recalibrate new raters or propagate calibration uncertainty. When anchors or cohort linking are required, select a route whose anchor and uncertainty contracts support the intended operation; do not pass an anchor argument that the two-family route rejects. The accepted RSM/PCM and scoped one-family MML anchor-related outputs keep their own restrictions.

Read older performance evidence by target, design and procedure, not only by sample size. A separate development comparison used 50 datasets per cell in the historical three-task/six-assessor/three-category model. For the common-Person assignment and generating ability SD of one, the original frozen procedure gave:

The common-Person assignment includes overlap and singly assigned candidates; it is not a complete crossing. These are inspected development results, with correlated targets and substantial Monte Carlo uncertainty, not independent confirmation or a universal sample-size rule. In particular, 20/20 inclusion does not erase the 30 unavailable datasets. The separate N20 rank-policy correction returns intervals from 3/50 retained sources; the other 47 still fail numerical or category checks. Fewer Persons than free parameters no longer causes an automatic model-Hessian refusal.

Using an already computed, separately qualified source increases N60 delivery to 50/50, but Task t1 slope inclusion remains 41/50; its location inclusion is 48/50. Local profile membership at the true slope gives 42/50, which does not demonstrate a coverage repair or completion of all interval endpoints. Thus improving convergence or changing to profile intervals is not by itself an answer to the sampling question. Keep severity contrasts, component slopes and conditional Person intervals distinct when deciding what to report.

For one question chosen in advance, you can also examine a profile likelihood interval. It fixes that component slope at successive values and reestimates all remaining slopes, locations and steps. The ability population stays fixed at N(0,1). This allows the likelihood shape to depart from the quadratic shape used by Wald intervals; it does not automatically remove bias or improve coverage.

# Choose an actual level of the fitted Assessor column, for example "A1":
# profile <- confint(fit_joint, parm = "slopes", method = "profile", slope = c(Assessor = "A1"))
# plot(profile) # likelihood loss, cutoff, and same-target Wald limits
# apa_table(profile, which = "profile_endpoints")
# res <- mfrm_results(fit_joint, include = c("fit", "plots"), compute = "never",
#                     intervals = list(profile = profile))
# mfrm_report(res)

The name Assessor identifies the owner; "A1" identifies its level. Use your own column/level names. Specifying both avoids confusion when two facets share the same level labels. A first-family component retains its geometric-mean-one reference; a second-family component is free on the fixed ability scale. The interval concerns one component, not the product of the two slopes. Neither owner name nor a large slope establishes rater competence.

Profiling is an explicit, more expensive sensitivity analysis. Inspect the endpoint and numerical tables, including failed searches. The integration grid may need refinement even when it suffices for Wald intervals at the estimate: the profile also visits other coefficient values. No automatic grid change or fallback to Wald occurs. Missing endpoints mean unresolved calculations, not infinite intervals. The plot retains failed points and separates them from connected passing sections; headings and notes can be removed with title = NULL, subtitle = NULL, and caption = NULL, without deleting saved checks.

In a synthetic 240-Person example, task t3 gave profile bounds 0.899–1.227 versus Wald 0.901–1.228; assessor r3 gave 1.070–1.585 versus 1.075–1.590. These close results verify a calculation and do not demonstrate superior coverage. Profiling does not make a numerically rejected near-zero-slope fit eligible or certify that alternative boundary solutions are absent.

Sparse assignments and a narrow ability distribution can leave slopes poorly determined. A fit may converge while a constrained search for an interval endpoint does not. Obtaining endpoints by increasing numerical accuracy does not establish that the intervals have their nominal coverage. In particular, evaluating coverage only among successful searches can overstate accuracy. There is no established coverage advantage over Wald intervals for this model.

Failed trials remain saved even when both endpoints are found. An inaccurate trial beyond an endpoint need not prevent computing that endpoint: the search can try nearer values. Failures within a bracket still prevent accepting the interval. This is why a plot may contain crosses for unresolved trial points and also show two successfully computed bounds.

If a profile is unresolved, first inspect the original fit’s convergence, then the profile-check and endpoint tables. The controls address different problems: maxit in fit_mfrm() controls EM iterations, whereas profile_control = list(maxit = ...) in confint() controls each constrained BFGS stage. A finer quad_points setting requires refitting and addresses integration error; it does not repair unsuccessful optimization or agreement between starts. Keep the original failures when documenting a planned analysis; a changed numerical procedure needs its own checks.

For a two-family fit that used 121 integration points, the following illustrates how to compare with 181 points and explicitly repeat the profile calculation. These counts illustrate the workflow; they are not universal accuracy settings. Use the same rating data used for the original fit.

# review <- mml_quadrature_sensitivity(fit_joint, ratings,
#                                      quad_points = c(121, 181))
# summary(review)
# profile_finer <- confint(review$fits$q181, parm = "slopes", method = "profile",
#                          slope = c(Assessor = "A1"))
# apa_table(profile_finer, which = "profile_endpoints")

review$intervals contains Wald intervals for each grid, not recalculated profile intervals. Use the explicit confint() call above to assess the selected profile. If the numerical tables instead show that constrained optimization reached its iteration limit, increase the profile iteration budget on the same fit:

# profile_longer <- confint(fit_joint, parm = "slopes", method = "profile",
#                           slope = c(Assessor = "A1"),
#                           profile_control = list(maxit = 1600))

An increased budget allows more computation without relaxing the numerical checks. It may still leave an endpoint unresolved. Report any changed settings and the original outcome; do not select the result simply because its interval is narrower or supports a preferred conclusion.

The checks re-evaluate the saved model, score derivatives, convergence and information. They also increase the quadrature grid from q to 2q-1 at the saved parameters. Passing is a numerical result, not a coverage guarantee or proof of global identification. A converged five-point fit of a retained example failed this integration check; its 31-point fit passed. This example shows why convergence alone is insufficient, not that 31 points always suffice. See ?confint.mfrm_fit for the tolerances. Use apa_table(ci_joint, which = "checks") for pass/fail reasons and which = "numerical_checks" for measured changes. Failed checks keep estimates with missing bounds. The experimental warning, owner identities and checks follow plots, reports and saved exports; the original fit’s primary-estimate/readiness fields remain unchanged.

Empirical Person-score rank remains a diagnostic, with an additional warning when deficient. It is not a prerequisite for the model-based interval: the observed Hessian used here is different from the empirical outer product of scores. For example, fewer Persons than free parameters forces deficient empirical score rank, but does not by itself prove that the observed Hessian is singular. A nonsingular Hessian still does not establish small-sample coverage or exclude boundary alternatives. Full-parameter empirical sandwich inference is unavailable; the model-based curvature, stationarity, derivative and integration checks still apply.

A small optimization residual warning means the retained solution has a nonzero marginal score, but its local Newton displacement is at most 0.01 in joint SE units. Values above 1e-4 now retain this warning. This is an approximation at the saved estimate: confint() does not refit or move it. The other mean-score, information, derivative and quadrature checks still apply. Do not repeatedly tighten em_score_tol without inspecting the result: very strict settings can exhaust the EM iteration limit at nearly unchanged estimates. A failed quadrature sensitivity check is a separate issue; use mml_quadrature_sensitivity() to refit with more points while retaining the same model and EM stopping controls:

# fit_joint above used 31 points. Include the original point count.
integration <- mml_quadrature_sensitivity(
  fit_joint, ratings, quad_points = c(31, 61),
  adaptive_quad_points = c(31, 61)
)
summary(integration)
integration$quadrature_review # fixed-calibration integration checks per Person
integration$slopes # estimates, raw diagnostic SEs, interval availability
ci_finer <- integration$intervals$q61 # already computed; no repeated fit
apa_table(ci_finer, which = "numerical_checks")
plot(ci_finer) # retains missing intervals and cautions
# integration$fits$q61 is the corresponding fitted model.

Here the 61-point interval checks also compare 61 with 121 points at the same saved parameters; that check does not refit at 121 points. Each requested refit reconstructs the source initialization policy from the same data; this fixed-grid EM example uses a neutral start. Older adaptive fits without a recorded policy also retain their earlier neutral-only initialization. The first family’s geometric-mean reference and the second family’s fixed-population reference are retained in the slope table. adaptive_quad_points adds a separate check at each fit’s unchanged parameters. It adjusts the integration grid to each Person’s posterior while preserving the original population distribution. This helps separate inaccurate integration from differences between fitted solutions. Compare both the fixed/adaptive results and successive adaptive orders in quadrature_review. Agreement between two orders is evidence of numerical stability, not a general error bound. No single point count is guaranteed to suffice.

Two-family Person-score comparisons remain unavailable in the main comparison table; missing EAP/SD changes mean unavailable, not agreement. Posterior moments in quadrature_review are diagnostic calculations at fixed calibration; they do not supply the separately checked predict_mfrm_units() scores or Person intervals described below. Requesting the review does not change EM’s fixed-grid estimation, the saved estimates or interval availability. Neither warning handling nor additional quadrature points establishes nominal coverage in a sparse assessment design.

An initial simulation checked two assignment designs at equal cost: 240 candidates, three tasks, six assessors and 1,440 ratings. A common-candidate design had 48 candidates scored by all assessors and 192 scored by one; a rotating design gave each candidate two adjacent assessors. Every assessor scored 80 candidates. Candidate assignment was independent of ability, and there were no fixed parameter anchors or missing assigned scores. The normal ability SD was 1 or 0.5; truth was transformed to the fitted fixed-N(0,1) scale.

Using the current small-residual warning rule and a finer-grid follow-up for the common-candidate SD=1 condition gave:

These are pointwise 95% intervals for two prespecified components, not a joint coverage statement. All nine components were summarized separately. The rotating-pair SD=0.5 task result, 90/99, has a Monte Carlo 95% interval of 83.4–95.8%; it does not establish nominal coverage. In that condition its empirical log-slope SD was 0.178 versus root mean estimated variance 0.155. Some log-slope bias also remained: common-candidate SD=0.5 assessor r3 had mean error -0.111 (Monte Carlo SE 0.032), despite 99/100 coverage.

Originally, 89 common-candidate SD=1 intervals failed the 31-versus-61-point integration check. Refitting all 100 datasets, including the 11 that had passed, at 61 points resolved these failures; all also passed the 61-versus-121 check. The EM stopping tolerance, iteration limit, response data and model were unchanged. This is a numerical correction using the same sample, not an independent replication. Task t1 coverage of 93/100 has a Monte Carlo 95% interval of 86.1–97.1%; assessor r1 coverage of 98/100 has an interval of 93.0–99.8%. These results do not establish nominal coverage or prove that 61 points will suffice in other datasets. Across all nine components in this condition, coverage ranged from 90/100 to 98/100. Task t3 had 90/100 coverage (Monte Carlo 95% interval 82.4–95.1%), with empirical log-slope SD 0.130 versus root mean estimated variance 0.110. Assessor r3 retained log bias -0.058 (Monte Carlo SE 0.016). Thus, removing integration failures did not remove the observed sampling discrepancies.

The original failures must not be described as an unidentified assignment design. Equal integration-point counts did not give equivalent numerical accuracy across these designs. The remaining rotating SD=0.5 failure had a near-zero assessor slope and larger local optimization displacement. The warning change was assessed on saved fits; this is not independent confirmation of the revised rule. These results support cautious exploration and numerical review, not a general guarantee for sparse designs or a recommendation that one allocation is better.

Estimate the same response model in Stan

To inspect prior sensitivity or undertake an external Bayesian analysis, export_mfrm_model() writes a self-contained Stan bundle from fit_joint. Export does not refit the MML model or run a sampler.

stan_bundle <- export_mfrm_model(
  fit_joint, "gmfrm-stan", format = "stan", log_slope_prior_sd = 1
)
model_data <- readRDS(stan_bundle$files[["model-data.rds"]])
model_data$priors
model_data$slope_levels
# Review README.txt and run.R, then run Rscript run.R from the bundle directory.

The exported likelihood retains the product of the two slopes multiplying the entire adjacent-category predictor, the fixed N(0,1) population, source location constraints and second-owner steps. The first family’s log slopes sum to zero. Their normal-coordinate priors give each component the specified marginal log-slope SD; constrained components are dependent. Second-family log slopes have independent normal priors. Locations and step priors remain explicit as in ordinary model export. The defaults are starting choices for prior review, not a demonstrated correction for MML coverage.

After sampling, component_slopes and log_component_slopes map to the owners and levels in slope-levels.csv. run.R saves draws in fit.rds before reporting. Its parameter-summary.csv contains explicit equal-tailed 95% credible bounds, endpoint Monte Carlo SEs, R-hat and bulk/tail ESS for both reported parameters and free sampling coordinates. sampler-diagnostics.csv records each chain’s divergences, tree-depth hits and E-BFMI; posterior-summary.rds retains these tables and interpretation notes. Read summary-decision.csv first for the sampling review and next action; sampling-review.csv records the checks and criteria. The parameter table joins original labels and reference/direction notes, with value_role and interval_status beside the numbers. Fixed values have no estimated SD or interval; their uncertainty fields are NA, not evidence of perfect precision.

From the exported directory, run Rscript run.R --summarize to rebuild these reports without compiling or sampling. It replaces the report files and preserves the saved draws. Bounds remain visible when diagnostics warn: interval_computed means finite endpoints for a varying parameter, not qualified inference. Endpoint MCSE measures simulation error in a reported bound, whereas posterior SD describes model-based uncertainty. endpoint_mcse_fraction is the larger endpoint MCSE divided by interval width; its acceptable size depends on the required precision. Good precision for one parameter cannot override divergent transitions or poor exploration elsewhere. Constant draws yield NA diagnostics; check whether a model constraint fixes the component or sampling is stuck.

Inspect all coordinates and chains, prior sensitivity and model fit before interpreting posterior intervals. These are Bayesian credible intervals, not MML confidence intervals or proof of frequentist coverage. Predictions condition on the sampled abilities of the existing Persons.

This uses the response equation of Uto and Ueno (2020), equation 9, while preserving mfrmr’s logistic metric (L-metric), with D = 1. It does not copy their Appendix’s extra 1.7 multiplier or its prior construction. Fixed N(0,1) abilities do not imply a normal-ogive link. When ability, locations, steps and the first slope family are held fixed, an otherwise matching D = 1.7 parameterization needs second-family slopes equal to the mfrmr slopes divided by 1.7 to preserve response probabilities. Dividing both families by 1.7 would not preserve them. Bayesian equivalence additionally requires transforming the priors; retaining the same normal log-slope priors changes the model. One-family GPCM exports to Stan or brms with its own geometric-mean-one slope constraint; two-family brms remains unsupported. Supported fixed facet anchors and group means are retained when exporting a one-family fit. These anchors carry no external calibration uncertainty. Native two-family fitting still excludes anchors. The files contain ratings and original labels and need review before sharing.

Inspect observed ratings with posterior residuals

Suppose you want to describe where a judge’s observed scores differ from the model’s predictions. mfrm_response_diagnostics() integrates each candidate’s ability posterior, using all of that candidate’s observed ratings and holding the fitted calibration fixed. Both slope families enter the response equation. It does not simply plug in the candidate’s EAP score. Its predictive variance includes both rating variation at a given ability and variation in the expected rating across the ability posterior.

# After the two-family fit has passed its engine's convergence checks:
# response <- mfrm_response_diagnostics(fit_joint, group_by = "Assessor")
# response$rows     # review unavailable rows and integration differences first
# response$measures
# plot(response, style = "paired")
# as_ggplot(plot(response, style = "scatter", draw = FALSE))
# res <- mfrm_results(fit_joint, include = c("fit", "plots"), compute = "never",
#                     response_diagnostics = response)
# mfrm_report(res)

Use your fitted identifier name for group_by. For adaptive direct MML, this uses a grid adapted to each candidate’s posterior; fixed-grid EM keeps its fixed grid. In both cases, all of the candidate’s observed ratings condition the prediction, even when rows returns only a subset. If response integration is unresolved, explicitly increase quad_points and inspect the new check; the fitting grid need not suffice for every response prediction. This recomputes probabilities at the saved calibration. Returned probabilities use the higher order (2 * quad_points + 1), and unresolved rows and groups remain visible. In contrast, mml_quadrature_sensitivity() refits the calibration across specified grids. Neither numerical check establishes model adequacy.

With rows, select original input row numbers for display and grouping. All observed ratings still inform each candidate’s posterior. Incomplete assignments stay incomplete; the function adds no unassigned ratings. A group with any unresolved observed row has unavailable summaries rather than silently dropping that row. Available output connects to paired/scatter plots, ggplot conversion, saved reports and exports without refitting or reintegration.

These are same-data descriptive summaries: calibration and checking use the same observations. Infit and Outfit have no established expectation-one reference here. A value above one is not an automatic misfit flag, and a value below one is not evidence that a judge is too predictable. No formal fit test, confidence interval or rater-quality classification is supplied. These residuals do not test unidimensionality or local independence and exclude calibration uncertainty. No slope confidence interval is needed to compute the fixed-point probabilities; its availability answers a different question.

Scoring new Persons with the fitted two-family calibration

predict_mfrm_units(fit_joint, new_ratings) supplies experimental EAP scores, posterior SDs and continuous equal-tail posterior intervals. It uses only the new response rows and holds both slope components, locations and category steps fixed. The scoring prior remains N(0,1); a prior override or population covariates are not supported. Non-Person levels must be known, scores must use the original zero-based coding, and observation weights must be one. Known levels may form new combinations under the product model, but these combinations have not thereby been checked empirically. Repeated Person-facet combinations are rejected; missing responses are omitted with a warning and row counts, and a Person without valid rows receives no score.

# new_ratings has the original Person, facet and score columns.
scores <- predict_mfrm_units(fit_joint, new_ratings, scoring_quad_points = 121)
summary(scores)
saveRDS(scores, "two-family-new-person-scores.rds")

The saved calibration must reproduce its specification, point estimates, likelihood and numerical convergence. A finer integration rule must change the NLL per calibration Person by at most 1e-6, with maximum absolute gradient per Person at most 1e-6. A separate check compares each reported EAP and SD against two adaptive reference orders (tolerance 1e-5). Source checks and scoring checks remain separately recorded. A numerical refinement failure requires explicit readiness_policy = "review" and keeps the affected result flagged; invalid saved calibration or failed convergence is rejected even in review mode. These checks do not require slope confidence intervals.

The posterior intervals exclude calibration uncertainty. They do not establish frequentist coverage, population transport, global identification, boundary absence or decision accuracy. Saved scores retain the calibration roles, category coding, numerical checks and limitations without recalculation. Their 95% posterior probability also differs from containing the fixed true ability in 95% of repeated samples at every ability value. With a known correct calibration and population law, exact posterior probability averages to 95% over that population; substituting an estimated calibration need not preserve that property. Tail-specific shrinkage and population-average performance are therefore separate questions.

For example, the retained common-Person SD1 study at calibration N=480 returned all 3,000 held-out population intervals and included truth in 94.2% (Monte Carlo SE 0.50 percentage points, based on 50 calibration replicates). Separate fixed-ability panels gave 84.1%, 98.0% and 86.1% inclusion at z=-2, 0 and 2. At calibration N=60 the population result was 91.2% (MCSE 0.64 percentage points), again with all 3,000 intervals returned under the prespecified integration-refinement rule. These development results do not identify calibration uncertainty as the sole cause, establish population transport, or justify treating Persons sharing a calibration as independent replicates. sample_mfrm_plausible_values() uses the same checked route for discrete posterior draws; these do not establish a downstream population-analysis method. Attach saved scores for results tables, reports and exports:

res <- mfrm_results(fit_joint, scores = scores)
res$tables$person_scores
res$tables$scoring_integration
report <- mfrm_report(res)

The same attachment accepts the format-6 scores below when the source fit is available. It checks saved source identity against the fit’s calibration, data, category/owner coding, integration and recorded status, without repeating numerical checks or scoring. Reports and CSV retain unrounded estimates and their conditional interpretation, separate source/batch checks and review-only labels. Portable results also retain the not-scored roster and row dispositions; native results record omitted row counts, without enumerating not-scored Persons. The attachment provides tables and reports, not a new score plot route. Older scores without source identity remain readable on their standalone route but must be regenerated before attachment; re-extract older portable calibrations from the saved fit. No refitting is required.

The same conditional target is available from a portable format 6 calibration:

draft <- extract_mfrm_calibration(fit_joint, scoring_quad_points = 121)
calibration <- freeze_mfrm_calibration(validate_mfrm_calibration(draft))
save_mfrm_calibration(calibration, "two-family-calibration.rds")

# In a later R session: the original fit and training data are not needed.
calibration <- load_mfrm_calibration("two-family-calibration.rds")
scores <- score_mfrm_calibration(calibration, new_ratings)
summary(scores)

Extraction requires passing source checks; review-only native scores cannot be turned into a frozen calibration. The artifact retains both slope components, location/step constraints, category codes, fixed prior and observed facet combinations, without training responses or Person estimates. Loading checks the stored identity and evidence; it cannot re-evaluate the absent source data. Every scored batch must pass integration checks, with no review override. row_dispositions$ObservedContext marks combinations absent during calibration. Prior overrides and repeated-event extensions remain unsupported. Missing responses require explicit missing_response = "omit"; Persons without valid responses are reported as not scored. File validation and freezing do not add statistical qualification or calibration uncertainty.

plot(fit_joint), ordinary fit diagnostics, Wright/Pathway maps, location/curve intervals and model ranking are not yet available. Unsupported fitting arguments produce an error rather than being ignored. Saving and reopening res preserves this analysis; it does not create a portable calibration artifact.

Which uncertainty claims are supported?

All the inferential routes below are approximations. A numerical eligibility check establishes whether a calculation can be returned under its stated rules; it does not certify the nominal confidence level for your assessment. The available studies support different claims for different targets:

For all routes, keep unavailable or unbounded outcomes and saved cautions in your report. Bonferroni concerns the finite family requested in one call; sandwich covariance requires the declared independent clusters. Neither option removes model bias or guarantees small-sample calibration. The APIs remain available with these limits; no data-dependent widening or arbitrary sample-size cutoff is used to conceal the adverse findings.

How uncertain is a relative slope?

Suppose three rubric criteria have relative slopes 0.8, 1.0 and 1.25. Their geometric mean is one. The last criterion has the steepest response curve on this model’s common ability scale. An interval describes sampling uncertainty in that relative discrimination; it does not say that the criterion should receive more scoring weight or that its raters agree more closely.

# After fitting a native GPCM with method = "MML":
# slope_ci <- confint(fit, parm = "slopes", level = 0.95)
# slope_ci
# attr(slope_ci, "diagnostics")[, c("SlopeFacet", "Estimate",
#   "CIEligible", "InferenceReview")]

The calculation inverts the full joint observed information, including estimated population parameters. It then maps its free log-slope covariance to all slope levels using the sum-zero constraint. The interval for slope \(a_g\) is \(\exp[\log(a_g) \pm z SE(\log(a_g))]\). Working on the log scale keeps bounds positive and accounts for the dependent slope coordinates. Inverting only the slope block of the information matrix would ignore estimation of the other parameters and generally understate uncertainty.

Scale constraints alone do not guarantee that the data identify the model. For example, suppose every Person answers the same two binary items once. The four response patterns have only three independent probabilities. The default GPCM still estimates four quantities: one item-location contrast, one relative log-slope contrast, and the ability population mean and variance. Increasing the number of Persons cannot resolve this deficit. The package retains descriptive optimizer points but withholds joint covariance, regular intervals and inferential model comparison with an explanation. A simpler model or additional informative response designs may resolve the problem; merely increasing quadrature or optimization effort cannot.

This dimension check uses each distinct rating roster and population-covariate row, including the number of independent repeated responses. It does not require more Persons than parameters and does not enumerate every response pattern. Passing this necessary check is not a global identification proof. A recorded complete-pattern first-order rank failure also prevents regular information-based inference; a tiny positive numerical Hessian eigenvalue cannot override it. An observed score subset with deficient rank is a different, inconclusive diagnostic and is not used as an automatic model-information veto.

The current likelihood and gradient must match a numerically adequate fit, and the information must be positive without regularization. Unit weights, supported score categories and at least 31 quadrature points are required. Missing bounds retain a reason; increasing iterations cannot repair every failure. Examine different starts and a denser grid for consequential results. diagnose_mfrm(fit)$parameter_uncertainty$slopes supplies the same 95% intervals; attach_diagnostics = TRUE also makes them available in the fitted slope table and subsequent weighting reviews. Ordinary Wright/Pathway maps display ability, facet location or fit, not these discriminations. Existing information plots remain plug-in curves; mfrm_curve_intervals() adds calibration intervals.

These defaults are pointwise model-based approximations. The optional methods below explicitly change the target, covariance or multiplicity rule. An interval for PopulationSD * Estimate needs the variance and covariance of the estimated population SD. A slope interval containing one is not a test that all slopes are equal; use the matched PCM/GPCM test for that question.

The following targeted check used saved fits from an earlier implementation. Later numerical repairs can change fitted solutions or interval availability; these percentages are not an evaluation of refitting every dataset with the current estimator. The check used 100 independent datasets per condition, with three criteria, three raters and three score categories. It crossed the slope owner (Criterion or Rater), equal versus unequal slopes (approximately 0.67, 1 and 1.49), and two designs. With 100 Persons rated by all three raters, intervals were available in every dataset; pointwise coverage across levels was 92–98%. With 40 Persons each rated by two of the three raters, equal-slope availability was 100% and coverage was 95–97%. Unequal-slope availability fell to 96–97%, with coverage among available intervals from 90.6% to 97.9%; the lowest covered-and-returned fraction over all planned datasets was 87%. Some intervals were extremely wide. Four datasets were withheld for regularized information and three for category-support review. These results support an approximate procedure with explicit limits, not a general 95% coverage claim. Each coverage estimate has only 96–100 independent repetitions; the three intervals within a dataset are dependent. Sample size and incomplete assignment changed together, so their effects cannot be isolated from this check.

Choosing the interval for the question

For rubric criteria, a comparison of two slopes asks whether their response curves differ in steepness, conditional on this model. It is not a comparison of rater agreement or a reason to change scoring weights. Prespecify the comparison before looking at which estimates differ most.

The following fictional rubric data give a reproducible example. The simulated model is correctly specified; this illustrates the API, not field performance.

library(mfrmr)
ratings <- simulate_mfrm_data(
  n_person = 80, n_rater = 3, n_criterion = 3, score_levels = 3,
  model = "GPCM", step_facet = "Criterion", slope_facet = "Criterion",
  slopes = c(C01 = 0.8, C02 = 1, C03 = 1.25), seed = 924090
)
fit <- fit_mfrm(
  ratings, "Person", c("Rater", "Criterion"), "Score",
  model = "GPCM", method = "MML", step_facet = "Criterion",
  slope_facet = "Criterion", quad_points = 31, maxit = 400, reltol = 1e-10
)

confint(fit, parm = "slopes", scale = "standardized")
slope_levels <- fit$slopes$SlopeFacet
comparison <- matrix(0, 1, length(slope_levels),
                     dimnames = list("first versus second", slope_levels))
comparison[1, 1:2] <- c(1, -1)
ratio <- confint(fit, parm = "slopes", contrasts = comparison)
difference <- confint(fit, parm = "slopes", contrasts = comparison,
                      contrast_scale = "difference")
attr(ratio, "diagnostics")[, c("Estimate", "CI_Lower", "CI_Upper", "PValue")]

The ratio null is one; the difference null is zero. More rows specify a finite family. simultaneous = "bonferroni" adjusts both intervals and the corresponding Wald p-values for all requested rows, not for other analyses or data-selected comparisons. Standardized slopes multiply the relative slope by estimated population SD, propagating its uncertainty and cross-covariances. With latent regression this is the residual population SD, not the marginal spread across different covariate groups. Ratios cancel a common scale; differences do not. JML inference is outside these MML methods.

confint(fit, parm = "slopes", method = "sandwich", simultaneous = "bonferroni")
# If Persons belong to independent schools, supply a complete Person/Cluster
# mapping and use clusters = schools. Do not cluster separate rows of one Person.

Sandwich covariance uses the marginal-likelihood score vector for each Person, or sums those vectors within a supplied cluster. It needs many independent clusters with enough score-vector rank. Its target is the fitted working-model parameter. It does not repair bias from a wrong population model, informative assignment or dependence between clusters. The optional adjust = TRUE only multiplies the covariance by G/(G-1); it is not a small-sample correction with guaranteed coverage. The earlier 800-dataset relative-interval study does not qualify these additional procedures by itself. An earlier reanalysis recalculated intervals on all 800 saved fits, without new simulations or optimizer runs, and found 90.7–99% pointwise coverage for standardized model intervals and 90.7–100% for standardized sandwich intervals. The minimum was 88/97 available intervals in a small incomplete, unequal-rater-slope condition (Monte Carlo 95% interval 83.1–95.7%; 88/100 covered-and-returned). For three standardized slopes, Bonferroni family coverage ranged from 91.7–98% with model covariance and 90–97% with sandwich covariance. Each cell has only 100 independent datasets. These correctly specified normal-population results show that neither sandwich covariance nor multiplicity adjustment guarantees nominal coverage in small samples; they do not test robustness under misspecification.

Examine one relative discrimination with a profile likelihood

The default slope interval uses the local curvature near the fitted solution. If the likelihood is asymmetric, that approximation can differ substantially from an interval obtained by fixing a candidate discrimination and refitting all remaining parameters. Development mfrmr offers the latter as an explicit, experimental option for one prespecified relative slope. It does not establish better small-sample coverage.

# Choose a level from fit$slopes$SlopeFacet before looking for a favorable result.
ci <- confint(fit, parm = "slopes", method = "profile", slope = "R01")
print(ci)
plot(ci, type = "profile")
# Hide text when composing a larger figure:
as_ggplot(ci, type = "profile", title = NULL, subtitle = NULL, caption = NULL)
apa_table(ci, which = "profile_endpoints")
apa_table(ci, which = "profile_checks", digits = 8)

# Keep the calculated profile; reporting does not fit the model again.
res <- mfrm_results(fit, intervals = list(profile = ci), compute = "never")
report <- mfrm_report(res)
saveRDS(ci, "one-slope-profile.rds")

The curve shows twice the log-likelihood loss. The horizontal dashed line is the requested chi-square(1) cutoff; its crossings define the connected local interval. Dotted vertical lines show the same-target Wald limits. The slope axis is logarithmic. Crosses identify unresolved calculations; lines stop at failed points, and evaluations without finite coordinates remain in the tables. A steeper discrimination is not a measure of rater agreement or a scoring weight recommendation.

The relative slopes still have geometric mean one. Fixing even the dependent last slope preserves that constraint, while the population mean and variance and other nuisance parameters remain free. Two starting values and a finer quadrature are checked at each candidate. A better source likelihood, failed numerical checks, observed nonmonotonicity or an exhausted search produce explicitly unresolved bounds. Increase profile_control$max_steps only to extend a stated search; failing to find a crossing does not prove an infinite interval. Successful finite searches cannot exclude other, disconnected regions.

The initial scope is an eligible GPCM MML fit with an estimated intercept-only normal population, unit weights, and no anchors or interactions, with shared or separate slope/step owners. Standardized targets, ratios/differences, JML, latent regression and simultaneous profile intervals are not included. The source fit and its global audit states remain unchanged. This method requires additional nuisance fits and is not assumed cheaper than bootstrap. Regular likelihood-ratio asymptotics still matter near boundaries; a profile is not a correction for a wrong model or a guarantee of interval coverage.

A comparison using all 100 saved calibrations from one small-sample condition illustrates why successful profiles alone are insufficient. Each dataset had 40 people, three criteria and a linked rotating assignment of two of three raters. Rater slopes and steps shared their owner. The target was the first rater’s relative slope, with true value exp(-0.4); both methods used nominal 95% intervals on the same stored estimate.

The exact 95% Monte Carlo interval for profile availability is 42.8%–63.1%; for coverage among available profiles it is 93.3%–100%. The latter describes the selected 53 datasets, not a reliable procedure for all datasets. Among the same 53 datasets, Wald covered 49; their median interval widths were 1.010 for Wald and 1.035 for profile. This comparison does not establish general profile superiority. Wald also remains approximate: its available-interval coverage has a Monte Carlo interval of 83.3%–95.7%, and some intervals were extremely wide.

Of the 47 unresolved profiles, 11 failed before profiling (source integration, information or category checks) and 36 failed during constrained searches. Those 36 calls did return result objects: 32 retain one computed endpoint (12 lower-only, 20 upper-only), and four retain neither endpoint. Thus 53 is the count of complete finite intervals, not the count of returned objects. Computed endpoints and warnings remain available in the result; apa_table(ci, which = "profile_endpoints") shows each side’s status. A missing opposite endpoint does not turn the computed side into a separately calibrated one-sided confidence interval. Of the 47 incomplete profiles, 45 have a separately available log-Wald interval; retain its method and cautions when displaying it, rather than labelling it as the missing profile interval. These are limitations of this numerical procedure on the stored calibrations; they do not show that a statistical interval cannot exist. No failed case was replaced, and no tolerance was relaxed. Times include the interval method’s checks, exclude original fitting, and were measured with two local workers.

These previously analysed data provide a diagnostic comparison, not independent confirmation of the current estimator. They do not assess separate owners, probability intervals or all sparse designs. Retain profile intervals as an explicit experimental sensitivity analysis; do not use them as an automatic Wald replacement or report their successful subset as a coverage guarantee.

A fitted-model bootstrap alternative

The parametric bootstrap asks what variation occurs if independent Persons and scores follow the fitted model on the same observed rating assignment. It preserves covariates and facet labels, and reestimates the free population parameters. Unassigned ratings stay unassigned. It does not test whether the assignment mechanism or the normal population assumption is correct.

boot <- bootstrap_mfrm_gpcm(fit, nsim = 499, seed = 92401)
confint(boot, scale = "standardized")
confint(boot, contrasts = comparison, simultaneous = "bonferroni")
boot$trials  # Includes every planned trial, its last stage and any reason/warning.
apa_table(boot, which = "checks")  # Returned fits: categories, convergence and information.
# boot$refit_draws retains returned parameters for diagnosis, not interval calculation.
# For an eligible matched PCM fit with the same population model:
lr_boot <- bootstrap_mfrm_gpcm(fit, null_fit = pcm_fit,
                              nsim = 499, seed = 92402)
lr_boot$test

The interval uses basic bootstrap error quantiles, on a log scale for positive slopes/ratios and on the original scale for differences. Failed refits and returned estimates rejected by eligibility checks are retained as unknown outcomes. The reported limits enclose all completions, possibly reaching zero or infinity; deleting failures could instead conceal the least stable cases. The LRT reports Monte Carlo precision and p-value bounds when any refit is unresolved. Its null-model draws cannot supply slope intervals around the alternative fit. Saved results support new confidence levels without refitting. More repetitions improve tail resolution, not model validity; inspect attr(confint(boot), "expected_tail_draws"), especially with many simultaneous comparisons. This is an approximate bootstrap method, not a profile-likelihood implementation or an exact small-sample test.

For basic slope bootstrapping, a category observed once is now a warning rather than an automatic exclusion. The source fit and each refit must still pass model-identity and numerical checks, including unregularized joint information. The exception requires a fresh audit confirming that every category is observed in each step scope and no step coordinate or category contrast is unsupported. It does not change the fit’s saved readiness or permission for Wald intervals, IC comparisons or LRTs. checks$BootstrapEligible and checks$WaldEligible distinguish these decisions. Source decisions are in source_checks. Accepted singleton warnings remain in the interval diagnostics, printed output, default plot subtitle and reports, even when limits are finite. A supplied subtitle, including NULL, overrides the default plot text.

The basic interval formula itself does not require a variance estimate in every replicate; variance estimates are required for studentized intervals, as explained in the R boot documentation. Separating those formulas does not establish that an unstable estimate is suitable for basic-bootstrap inference. Use trials$Stage, trials$AlternativeReturned and checks to distinguish the reasons. Missing check fields mean that the quantity was not recorded. Older saved results may have only a combined reason and no rejected parameter vectors.

A small-sample example illustrates why unresolved replicates matter. Its original 499-replicate run accepted 408 estimates. Rechecking 62 category-support cases under a revised acceptance rule produced an assembled result with 470 accepted and 29 unresolved estimates. That assembled result is a comparison of procedures, not a complete bootstrap run under the current estimator. All three standardized 95% intervals still extended from zero to infinity. This is one dataset, not evidence of repeated-dataset coverage or a failure rate for other designs.

Preserve an original bootstrap when changing estimation or acceptance rules. A complete new run is needed to evaluate the changed procedure; replacing selected failed draws cannot establish its performance. Saved reanalysis history, when present, follows printing, slope intervals and reporting tables. Missing history does not establish that an older result used the current procedure. Neither more replicates nor selective replacement automatically resolves poor identification or validates interval coverage.

Why an information warning needs a parameter-specific review

Inspect boot$checks (or apa_table(boot, which = "checks")) for PopulationSD, MinimumStandardizedSlope and MaximumStandardizedSlope. These describe the returned optimizer estimates, including rejected refits. For an estimated population, standardized slopes multiply relative slopes by the residual population SD; for a fixed standard-normal population that SD is one. Missing entries indicate missing retained estimates. These fields are diagnostics, not thresholds for judging raters or permission to use a CI. APA check tables use significant digits for these values and the gradient/ information diagnostics, preserving small positive values in table conversion. The original boot$checks columns remain unrounded numeric values.

The 29 unresolved trials illustrate different reasons for withholding inference. Independent likelihood and derivative calculations helped separate finite numerical stopping points from evidence for boundary approaches. These checks diagnose the selected examples; they do not certify global maxima or repeated-sample interval accuracy.

These are primary findings; numerical and boundary problems can overlap. An unused category in a replicate is a sampling outcome, not evidence that the scoring category should be removed. A near-zero standardized slope can coexist with other finite standardized slopes. Keeping the geometric mean of relative slopes at one can then make relative slopes and location coordinates extreme. Conversely, a stable discrimination estimate does not establish a finite threshold. The assembled example retains all 29 unresolved trials. The two improved finite solutions are separate numerical evidence; they have not been inserted into its draws or counted as a current-procedure bootstrap result.

For small fixed-grid GPCM MML fits (up to 64 free parameters, with reltol <= 1e-9), the direct optimizer now reviews numerical curvature even after a small raw gradient. Negative curvature triggers up to three BFGS restarts in rescaled search coordinates. Each uses the requested iteration ceiling. Only a non-worsening, numerically converged candidate without detected negative curvature can replace the selected point. Otherwise the estimate is retained with a numerical warning. Inspect fit$opt$optimizer_polish$Stages: StageLabel == "curvature_rescale" identifies these attempts, and SmallestCurvature, CurvatureScale and CurvatureReviewError explain the curvature review. Missing curvature entries mean it was not evaluated. The search rescaling does not provide a regularized covariance estimate. An improved optimizer result still needs separate information and integration checks; positive numerical curvature alone is insufficient.

Poor conditioning does not by itself mean that an information matrix must be regularized or discarded. It can reflect parameter units, weak information, or numerical differentiation error. The shared MML calculation now compares two refined Hessians, checks their unregularized inverse and reviews the gradient in the curvature metric. If these checks pass, intervals can be returned with a saved caution. Otherwise they remain unavailable. The information workspace budget covers the additional calculation too.

In the remaining finite-solution example, the standardized log-slope SE changed from about 30 with the default numerical differentiation to about 27 after refinement, agreeing with the reference parameterization. Nevertheless, the weakest discrimination’s 95% log-Wald interval spans roughly \(10^{-25}\) to \(10^{21}\). Returning that interval makes the lack of precision visible; it does not make the approximation trustworthy near zero. Examine attr(intervals, "cautions") and InferenceReview, which also appear in printing, default plot subtitles and reports. This numerical review does not establish sampling coverage or turn a boundary limit into a finite estimate.

Further requirements include integration and parameter-specific boundary estimates and inference. The examples support different treatment of different targets; they do not validate ordinary Wald intervals or demonstrate repeated-dataset coverage.

Uncertainty in probability and information curves

For example, examine a single rater/criterion context across known ability values. The existing probability and information plots still describe fitted points. This API propagates the full calibration covariance into each curve point, using logit intervals for probabilities and log intervals for information.

grid <- expand.grid(Theta = seq(-3, 3, length.out = 41),
                    Rater = "R01", Criterion = "C01")
curves <- mfrm_curve_intervals(fit, grid)
plot(curves, title = NULL, subtitle = NULL)

information <- mfrm_curve_intervals(fit, grid, type = "information")
plot(information)

Crosses identify fitted points whose intervals are unavailable. The caption counts these points; a missing ribbon does not mean negligible uncertainty. Ribbons stop at unavailable grid points. Inspect plot_data(curves, "table") for the saved InferenceReview reasons. Use caption = NULL to hide the caption; the crosses and saved reasons remain. Long default warnings wrap within the figure.

For slope or ratio intervals spanning many orders of magnitude, a logarithmic axis can make the estimates easier to compare:

# When estimates and finite interval bounds are strictly positive:
# as_ggplot(ci) + ggplot2::scale_x_log10()

Keep a linear axis for signed differences. A logarithmic display does not reduce uncertainty or establish interval accuracy.

Replace these labels with levels in the fitted model and include every non-Person facet. Information is per rating: slope squared times the conditional category variance. It is not total information summed over the observed exposures. Theta is fixed on the fitted native ability scale; uncertainty in a Person score is not included. Lines and ribbons retain their numerical table, use color plus line type, and accept custom palettes and omitted titles. Sandwich covariance is optional. Bonferroni adjustment is applied over the finite output grid, including all categories; connecting points does not produce a simultaneous confidence band for every value between them.

What the additional interval checks found

Do the slope-interval results also describe uncertainty in differences and curves? To examine that question, a further reanalysis used all 800 saved fits, with the current interval calculations and their availability checks. It did not generate new datasets or refit the historical estimates. Numerical repairs to estimation can change those estimates, so this is evidence about the saved-fit analysis, not an independent evaluation of the current estimator.

The generating facet effects were saved for every dataset. They were centered to the fitted reference before comparing curves at native ability values Theta = -2, 0, 2. Each slope-owner level was evaluated with the other facet at its second level. This gives nine rating contexts, 27 category probabilities and nine information values. Standardized differences compare the first with the second slope and the second with the third; they include population-SD uncertainty. Model and independent-Person sandwich covariance used the same fits. The confidence level was 95%, without a small-cluster adjustment.

The table gives ranges across the eight original conditions and the specified targets. Pointwise coverage concerns one target at a time. Family coverage requires every interval in the specified call to cover its truth; a family contains two differences, 27 probabilities or nine information values.

These percentages condition on available intervals or complete available families. Each condition has 100 independent datasets, not one independent replication per curve point. Depending on the target and condition, 96–100 pointwise intervals were available. In the small incomplete unequal-rater-slope condition, the sandwich probability family covered in 82 of 98 available datasets: 83.7%, with a Monte Carlo 95% interval of 74.8–90.4%. It was covered and returned in 82 of all 100 planned datasets. That Monte Carlo interval describes this cell, not simultaneous uncertainty for the minimum across cells.

Bonferroni does not fix an inaccurate marginal approximation. The curve results show a limitation even under the generating model’s assumptions; sandwich covariance is not an automatic solution. Read a returned interval as an asymptotic approximation, especially in small incomplete designs. Do not interpret numerical eligibility as demonstrated 95% coverage, turn curve differences into rater-quality verdicts, or assume that the bootstrap has resolved this limitation without its own repeated-dataset evidence. The two source designs also change sample size and assignment together, so the results cannot isolate an effect of sparsity. They do not establish performance for other grids, misspecified models or the latest refitted solutions.

Checking the updated estimator on the affected design

The entire small incomplete, unequal-rater-slope condition was then refitted: all 100 datasets, including cases whose intervals missed the truth or were unavailable. The updated numerical implementation used ordinary initialization and the same model, 61-point quadrature, iteration ceiling and convergence tolerance. All fitted parameters, objective values, probability estimates and interval endpoints were identical to the earlier results. Model and sandwich Bonferroni families still covered 83/98 and 82/98 available datasets, respectively. Updating the optimizer therefore did not resolve this shortfall in the selected condition.

These are the same 100 datasets already included in the 800-fit study, not another 100 independent replications or an independent confirmation after selecting a condition. Analytic probability derivatives and full-covariance delta intervals were also checked on two prespecified refits. The largest endpoint discrepancy from the API was below 2e-9. This checks the probability transformation; it does not establish global optimization, quadrature accuracy or finite-sample coverage, nor identify approximation error as the only cause.

For a small incomplete assessment, use these bands as approximate descriptions of calibration uncertainty. Do not rely on their nominal 95% label alone to claim that an entire curve is covered or that two curves clearly differ. Such a conclusion needs an uncertainty method justified for the study’s design and model. A larger sample, a simpler model or bootstrap may be worth assessing, but none is an automatic correction established by this check. Printing and default figure subtitles now identify the intervals as approximate; CIEligible continues to describe numerical availability.

Carry the selected inference into reports and figures

A rater’s severity, discrimination and residual fit answer different questions. Wright maps locate fitted abilities and facet effects; Pathway maps relate location to fit statistics. A slope interval instead describes discrimination under its stated scale and covariance method. Adding a slope interval must not replace severity intervals, Infit/Outfit, or a feedback decision with a different target. In particular, a larger slope does not automatically mean a better rater or a larger recommended scoring weight.

ci <- confint(fit, parm = "slopes", scale = "standardized")
results <- mfrm_results(fit, intervals = list(slopes = ci, curves = curves),
                        include = c("fit", "plots"), compute = "never")
plot(results, type = "gpcm_slopes", title = NULL)
apa_table(ci)
head(plot_data(results, type = "gpcm_curves", component = "table"))

The reproduction code in summary(results)$reproducible_code saves and reloads the complete result, preserving attached inference. In that code, res denotes your result object (res <- results here); choose the file path before saving. The starter export index also links the saved GPCM inference figures alongside the ordinary maps, with separate interpretation guidance.

For separately computed one-family Person scores, use the same collector:

# new_ratings contains response rows for the Persons being scored, using known
# facet levels and the fitted category coding.
# scores <- predict_mfrm_units(fit, new_ratings)
# scored_results <- mfrm_results(fit, scores = scores,
#   include = "fit", compute = "never")
# scored_results$tables$person_scores
# scored_results$tables$scoring_priors
# scored_report <- mfrm_report(scored_results)
# export_mfrm_results(scored_results, output_dir = "person_results",
#   include = c("tables", "report", "replay"), acknowledge_sensitive = TRUE)

Here Estimate is EAP and SD is the Person’s posterior SD. It is different from population SD, the SD of EAPs across Persons, and a calibration-parameter standard error. Lower/Upper are continuous equal-tail posterior intervals at the saved level; the ordinary diagnostic normal bands remain separate. Both condition on the point calibration and scoring prior and exclude their estimation uncertainty. Posterior probability is not a guarantee of repeated- sampling coverage for each fixed ability.

The attachment checks source calibration, population coding, prior and numerical records without rescoring. scoring_priors preserves retained and selected priors by Person, including an explicit scoring_prior override or differing covariate-based means. Review-only scores and omission records remain visible. Native scoring does not enumerate an entire not-scored roster, so NotScored is unknown rather than zero. Reports and saved replay retain the complete scoring object. This does not add a score plot route or change diagnostic maps. One-family portable calibration scores still use their standalone methods.

The same stored values are available through mfrm_report(results) and export_mfrm_results(results, ...). Named results define lowercase plot routes. Use as_ggplot(results, type = "gpcm_curves") for further ggplot customization; title = NULL or subtitle = NULL removes the corresponding annotation. Plots distinguish unavailable intervals by crosses and unbounded intervals by open circles, retaining their full limits and reasons in the tables. Color and line type distinguish probability curves. No reference cutoff is chosen for a slope plot; reference is an explicit optional plotting choice.

Source matching checks fitted data, parameters, population and integration settings before attachment. Reusing a saved bundle does not change the selected method, recompute intervals or rerun a bootstrap. For bootstrap slopes, pass confint(boot, scale = "standardized", level = .90) to preserve that choice; passing the raw boot object selects its default relative 95% pointwise intervals. The export replay reloads the saved RDS. The full archive contains the fitted data and identifiers, so review it before sharing.

Check that the integration grid does not change the decision

mml_quadrature_sensitivity(fit, data, quad_points = c(41, 61)) preserves explicit population formulas, covariates and factor coding. Its slope table contains both interval endpoints and eligibility; the summary reports their changes. For a PCM/GPCM comparison, apply the helper to each model and rerun compare_mfrm() on their fits at each common grid. Compare the IC/LRT decisions as well as estimates; increasing the GPCM grid alone changes the comparison basis. A small observed change is evidence for that fit and grid, not a universal integration certificate.

Information criteria and parameter intervals need different evidence. The R documentation for AIC requires maximized likelihoods on comparable data. Computing an AIC number does not require a slope confidence interval. In mfrmr, the comparison uses the declared parameter count, numerical solution and integration settings. Users should assess competing starts and integration sensitivity, particularly when a model-choice decision is close. An interval procedure also needs a justified covariance or other uncertainty method, the correct constrained-parameter transformation, and evidence about its repeated-sampling behavior. Neither decision requires a theorem covering every conceivable GPCM likelihood path.

An unavailable slope interval does not automatically invalidate a separately eligible IC comparison or matched PCM/GPCM test. Older saved fits must not gain intervals from stored flags alone; rerun confint() or diagnose_mfrm() to check the retained solution and refresh its uncertainty tables.

Why intervals, information criteria and tests need separate checks

Slope confidence intervals. The package uses the joint covariance in free coordinates, maps it to sum-zero log slopes and exponentiates the log-Wald limits. Its separate qualification checks preserve missing results for singular, regularized or unstable information. The result is an asymptotic approximation, not exact coverage for every possible dataset. The mirt documentation also provides GPCM fitting and information-matrix SE methods; that establishes a methodological precedent, not numerical validation of mfrmr.

Information-criterion comparison. The package already calculates \(\mathrm{AIC}=-2\widehat\ell+2p\) and \(\mathrm{BIC}=-2\widehat\ell+p\log N\), with \(N\) equal to the number of Persons for the ordinary MML route and \(p\) the free-parameter count. For GPCM, the comparison reevaluates the retained likelihood, gradient and observed information without rerunning the optimizer. It checks agreement with the stored objective, the existing gradient tolerance (at most \(10^{-4}\)), and positive unregularized information using the same relative eigenvalue tolerance as covariance inversion. IC/LRT and slope intervals share the full joint information calculation. The old IC/LRT-specific 80-coordinate cutoff is replaced by a dense-matrix workspace budget, configurable with options(mfrmr.max_information_bytes = 256 * 1024^2). This estimates eight double-precision square matrices; it excludes data, integration arrays and other process memory. An unavailable check returns its reason. This local check is separate from InferenceReady and does not qualify a slope interval. Same-data/likelihood, free-parameter-count, unit-weight and integration checks still apply. There is no requirement that the candidates be nested. Local curvature does not prove global optimality, and sorting candidate IC values is not a test that the winning model is true.

PCM/GPCM likelihood-ratio test. For the current relative-slope parameterization, let \(\lambda_g=\log\alpha_g\), with \(\sum_g\lambda_g=0\). The PCM null is \(\lambda_1=\cdots=\lambda_G=0\). It imposes \(G-1\) independent restrictions when the population model, steps, other facets and constraints are identical. The null value \(\alpha_g=1\) is inside the positive slope space; it is not a zero-variance boundary. Under the usual identified, interior and regular MML conditions, the candidate reference is \(2(\widehat\ell_{GPCM}-\widehat\ell_{PCM})\sim\chi^2_{G-1}\) asymptotically. The implemented test verifies the matched null/alternative population model, step and facet constraints, fixed interactions and free dimensions. It reuses the likelihood, gradient and unregularized-information checks described above. A negative likelihood gain is a numerical failure and does not become p = 1. The reference assumes independent Persons and correctly specified response and population models, with regular interior solutions. It is an asymptotic approximation, especially uncertain in small or weakly linked designs. Default calls can differ in population assumptions, so changing only model = "PCM" to model = "GPCM" does not necessarily test slopes alone. The mirt model-comparison documentation distinguishes ordinary likelihood-ratio comparisons from boundary-null tests. That distinction does not require a mixture chi-square test merely because GPCM slopes are positive.

A simulation illustration examined the approximation under a true PCM, using three raters, three criteria, three score categories and normally distributed Person abilities. In the crossed design, each of 100 Persons was rated by all three raters on all three criteria. In the rotating design, each of 40 Persons was rated by two of the three raters on all three criteria. The latter changes both sample size and assignment density, so it does not isolate the effect of sparsity. Rater and criterion locations varied across replications; the slope null was always true. Both fitted models estimated the population mean and variance on the same 61-point integration grid.

All 400 pairs supplied a test in these conditions. The intervals above are exact binomial intervals describing uncertainty in the simulated rejection rates, not parameter confidence intervals. The point estimates do not show a large departure from 5%, but 100 replications per condition give limited precision. They do not guarantee calibration under other sample sizes, category support, dependence, assignment or population distributions. In the first replication of each condition, increasing both fits to 121 integration points changed the likelihood-ratio statistic by less than \(10^{-6}\); this illustrates a sensitivity check, not a universal integration-error bound.

These explanations separate missing mfrmr functionality from failures of a particular fit. They do not imply that every model needs every downstream workflow, or that JML inherits the MML reference distribution.

When would two slope families help in practice?

Consider a speaking assessment: each candidate completes several tasks, and several assessors score overlapping sets of performances. Difficulty concerns which tasks tend to receive lower scores. Severity concerns which assessors tend to give lower scores. Discrimination asks a different question: how strongly do the probabilities of higher categories change with ability? Two slope families let that change depend on both the task and the assessor. The provisional fit_mfrm() route supports exploring this formulation. Experimental component-slope intervals are available when their numerical checks pass. The separate mfrm_facet_intervals(fit, "Rater") route supplies experimental normal location intervals or named within-facet contrasts from the inverse full marginal information. It retains the fitted location constraints and numerical refusal reasons; its sampling coverage is not qualified by the component-slope evidence. Attach it through mfrm_results(fit, intervals = list(locations = ci)) and select plot(res, type = "gpcm_locations"). Location differences do not imply uniform expected-rating differences when slopes or category steps also vary.

For a comparison in population SD units, supply an explicit zero-sum contrast:

levels <- fit$config$facet_levels$Rater
C <- matrix(c(1, -1, rep(0, length(levels) - 2)), nrow = 1,
            dimnames = list("First minus second rater", levels))
ci <- mfrm_facet_intervals(fit, "Rater", contrasts = C,
                          scale = "standardized")
ci$table
plot_data(plot(ci, draw = FALSE), "notes")

The standardized target is \(d/\sigma\), where \(d\) is the native severity contrast and \(\sigma\) is the population SD. When SD is estimated, the full joint delta method includes both its uncertainty and its covariance with the contrast. Dividing native limits by the fitted SD would omit these terms. This option supports the existing model-based RSM/PCM and GPCM routes with fixed N(0,1) or an estimated intercept-only normal population. A regression’s residual SD is not a marginal population SD, so latent regression is excluded. In the fixed-N(0,1) two-family model the known SD is one and the contrast is numerically unchanged. The saved result records the scale, native estimate, Jacobian and SD. Plot axes identify population SD units; explanations remain returned notes. Numerical eligibility does not establish finite-sample coverage, and a near-zero SD can make the ratio unstable. These contrasts align ability units; they do not make different models’ expected-rating differences equal.

A retained comparison illustrates why model complexity and interval usefulness need separate checks. It used eight PCM-truth conditions, fifty datasets each, three raters, three criteria and two raters per Person, comparing PCM with Criterion-slope GPCM while both estimated the normal population. The equal-slope truth belongs to both models. At N20 and generating SD1, Criterion-difference RMSE on the same 47 interval-returning datasets was .3113 for PCM and .4006 for GPCM; the paired MSE increase was .06357 (Monte Carlo SE .02815). Rater-difference RMSE instead decreased from .4142 to .3846. Extra slopes did not consistently improve these targets.

At N20 and generating SD.5, Criterion intervals contained truth in 43/44 PCM and 39/41 GPCM returns, but mean widths were 6.98 and 8.24 population SD units. The corresponding returned-and-containing counts were 43/50 and 39/50. Already refused near-boundary fits also produced extreme finite ratio points; these were retained in the all-attempt results. Thus high coverage among returned intervals does not by itself establish useful precision. These are descriptive saved-source results, not independent coverage confirmation, current-source readmission, or a rule for selecting a model. They do not show whether simplifying a model with genuinely unequal slopes is beneficial.

Neither interval route establishes a rule for feedback decisions; curve intervals remain unavailable. Individual two-family sheets can now collect the saved descriptive evidence through mfrm_report(..., style = "rater", facet = "Rater", rater = "r1"). Select the actual fitted facet and level; the sheet does not classify competence or demonstrate a training effect.

For a simple numerical illustration, let a task slope be 0.8 and an assessor slope be 1.2. Their effective slope is 0.96. Combining the same task with an assessor slope of 0.7 gives 0.56. These are multipliers in the adjacent-category log-odds equation, not score weights, accuracy percentages or observed score increases. The actual expected-score change also depends on ability, locations and category steps. In particular, a large slope can coexist with little information where almost everyone receives an extreme category. Read the probability and information curves over the ability range relevant to the assessment, not just a ranked slope table.

The multiplication is also a substantive restriction. The task/assessor sensitivity factors combine as a product; they are not a freely estimated slope for every task-by-assessor pair. A particular assessor’s unusual reaction to only one task may require a different model. Nor do two slopes create two latent abilities: the formulation discussed here has one ability dimension.

A useful analysis starts with the decision and rating design:

1. Define the ability being assessed, the meaning of each facet and the proposed feedback or task-review decision. Keep the column labels distinct from the chosen slope reference and step owner.
2. Inspect who rated which performances. Overlapping performances and multiple ratings per person help separate effects. A connected assignment graph alone does not establish that both slope families are identified; inspect weak links and systematic differences in the assigned candidate groups too.
3. Compare equal-slope and simpler one-family explanations on matched data, response equations and scales before attributing a pattern to two families. Supported one-family PCM/GPCM comparisons are available now; the current comparison API does not qualify a two-family test or ranking.
4. For a qualified two-family analysis, inspect uncertainty, starting-value and integration sensitivity, curves and substantive evidence together. Mark predictions for combinations not observed together. Use the findings to investigate tasks or feedback needs, not as an automatic removal or competence rule.

Start from an actual rating table

Suppose an assessment coordinator wants to review how raters use a writing rubric. The question is whether differences concern consistently higher/lower scores, category use, or different sensitivity to writing ability. A two-slope model is a candidate explanation when both criteria and raters may differ in sensitivity. It is not automatically preferable to an equal-slope or one-slope-family model, and a high rater slope is not evidence of accurate scoring. Reference performances and a substantive review of the rubric are needed before deciding what feedback to give.

The separately installed sirt package supplies the empirical data.ratings1 table from a 2009 Austrian grade-8 German writing survey (sirt data documentation). In sirt 4.2.133 it contains 135 students, seven observed raters, five criteria (k1–k5) and scores 0–3. The stored rater factor has additional unused levels; these are not additional observed raters. The documentation does not provide descriptive names for the five criteria, so retain their labels instead of inventing rubric meanings. Using this dataset does not make mfrmr’s response equation equivalent to sirt::rm.facets().

The following conversion preserves every supplied score and identifier. It requires sirt to be installed; it does not download data during analysis.

library(mfrmr)
data("data.ratings1", package = "sirt")
writing_ratings <- do.call(rbind, lapply(paste0("k", 1:5), function(k) {
  data.frame(
    Person = as.character(data.ratings1$idstud),
    Rater = as.character(data.ratings1$rater),
    Criterion = k, Score = data.ratings1[[k]]
  )
}))
writing_review <- describe_mfrm_data(
  writing_ratings, person = "Person", facets = c("Criterion", "Rater"),
  score = "Score", rating_min = 0, rating_max = 3,
  category_policy = "preserve", rater_facet = "Rater"
)
writing_review$design_connectivity
writing_review$structural_missingness$summary
# Count distinct raters for each student, independently of criterion rows.
writing_pairs <- unique(writing_ratings[c("Person", "Rater")])
table(table(writing_pairs$Person))
xtabs(~ Rater + Score, writing_ratings)

There are 1,370 observed ratings: 89 students have one rater, 27 have two, two have six, and 17 have all seven. The rater network is connected, but its information is concentrated unevenly across students. No score is missing within the supplied rows. That does not establish that every intended rating was collected: a separate planned assignment roster is unavailable. The unobserved combinations must not become zeros or automatically be labelled missing at random. Each criterion appears with each rater, but five rater-by-criterion-by-category counts are zero. These observations matter more than a single overall sample-size recommendation.

To examine criterion and rater slopes, use the preceding two-family call with data = writing_ratings, person = "Person", facets = c("Criterion", "Rater"), score = "Score", slope_facet = c("Criterion", "Rater"), step_facet = "Rater" and noncenter_facet = "Rater". Keep the declared 0–3 categories. This assumes one writing ability, a common standard-normal population, conditional independence of ratings, and rater-owned steps shared across criteria. Five analytic scores on the same performance may challenge independence; the fitted slope table does not test that assumption.

If your records distinguish Task, Criterion and Rater, preserve all three roles during design review. The current two-family route admits exactly two non-Person facets. Pooling task and criterion labels or dropping a facet changes the model and can conceal dependence; such a dataset needs an explicit model-scope decision before fitting.

An analysis of the complete table at 61 and 121 quadrature points converged under the same EM stopping rule, yet the maximum category-probability difference on a common ability grid from -4 to 4 was about 0.17. This compares fitted curves, not empirical predictive accuracy; the discrepancy also occurs within -2 to 2. Component-slope intervals were unavailable for both fits. At the 61-point calibration, the response-diagnostic check at 121 versus 243 points retained only 318 of 1,370 rows under its probability tolerance; all grouped Infit/Outfit values remained unavailable. Thus this example does not support issuing rater feedback from the returned estimates. Integration and solution sensitivity need resolution first. More convergence messages, a low in-sample residual, or a ranked slope plot cannot replace that review.

Holding these two saved calibrations fixed, an adaptive 61-point review was checked against separate continuous integration for all 135 Persons. The maximum differences were below 2e-8 for log marginal likelihood and 8e-8 for posterior mean/SD. The original fixed-grid likelihoods contained enough integration error to reverse which saved calibration had the higher likelihood. This identifies a numerical approximation problem; it does not establish that either calibration is the maximum of the accurately integrated likelihood, nor does it repair unavailable intervals. The same generic adaptive_quad_points option above supports this distinction for other data.

A subsequent direct adaptive fit on the same data was checked against continuous integration and local observed information. Its posterior residuals now use the same adaptive integration method. At this saved calibration, response checks with quad_points = 61 return 1,369 of 1,370 rows; quad_points = 121 returns all rows and all 12 criterion/rater summaries. The returned probabilities at the latter setting agree with separate continuous integrals within 1e-12. This is numerical agreement at fixed calibration, not evidence of predictive accuracy, calibrated fit cutoffs or the correctness of a one-ability model. Keep the unresolved row from the lower-order run visible; these results do not make 121 a universally sufficient order or establish that feedback would improve rating practice.

When these saved diagnostics are attached to mfrm_results(), the response_overview table retains the selected and available row counts, unresolved calculations, missing scores, zero-variance rows and unselected source rows. The report marks the 61-point result as requiring review and the complete 121-point result as descriptive with caveats. These labels do not change the saved calculations or establish model adequacy.

For a different operational setting, Uto et al. (2024) publish medical interview OSCE scores and their data description on Dryad: two raters cover every examinee and three additional raters cover different subsets. This asks whether rubric use and rater effects differ in a small assessment cohort. Their published model includes rater-by-item location interactions and item-specific steps, which the current two-family route does not implement. A native fit would be a restricted analysis, not a reproduction of that study. The writing example does not qualify OSCE, music or sport applications by renaming facets.

What if the analysis could affect a championship?

First identify the decision. The winner under the competition’s scoring rules, the athlete with the highest modelled ability, and the winner of a future performance are different targets. A model can help review judging or plan future assessments without supplying a replacement competition result.

The official 2026 Olympic women’s figure-skating results provide a concrete example. Adding the published short-program and free-skating scores gives the following decomposition for the first two finishers:

Both have zero deductions. The 1.89-point total difference is an official score difference, not a latent-ability estimate or its uncertainty. The component-score ordering differs from the total-score ordering because the competition also rewards technical elements. Fitting only the component marks would answer a narrower question; it would not determine who should have won.

The official short-program and free-skating protocols contain 29 and 24 performances, respectively, with nine judges and three program components per performance. Reconstructing the trimmed means, component factors and rounding reproduces all 53 program-component totals; combining them with the published technical scores and signed deductions reproduces all 24 finalists’ totals. This verifies data interpretation, not the accuracy of a latent model. Technical GOE conversions were not rederived.

Preserve the following information before fitting:

• Judge identity and performance occasion. The two published panels have 13 distinct judges, five shared. Join actual names/IDs from the short-program and free-skating rosters; a slot such as J1 is not the same person across segments.
• Advancement and assignment. The five non-finalists have no free-skating performance in this event. Those absent records are not missing assigned ratings. Advancement also selects the population observed in the final.
• Score meaning. Component marks use quarter-point increments from .25 to 10; integer recoding preserves a 40-category ladder. Sparse category use still needs review. GOE grades, difficulty/base values and deductions are different quantities and cannot be pooled as one ordinal response.
• Model scope. Athlete, segment, component and judge are distinct roles. The current two-family route admits only two non-Person facets. Discarding the segment or treating each performance as an independent athlete changes the model. Shared-performance dependence, one-dimensionality and the selected athlete population also require justification.

For winner selection, average parameter recovery and overall rank correlation are insufficient. Validation must ask how often the procedure selects the wrong winner, whether it wrongly presents that choice as certain, and how often it cannot distinguish the leaders. Include near ties, exact ties, a clearly separated leader, weak panel links and performance dependence. Keep unavailable results in the denominator. Changing a judge or a model in a sensitivity analysis is not a random draw from future competitions.

An individual 95% interval does not give a 95% confidence statement about being first among many athletes. Even a difference between two estimates needs their covariance: its variance is Var(A) + Var(B) - 2*Cov(A,B). Selecting the apparent leaders after seeing the data also requires inference that accounts for that selection. Shared calibration uncertainty cannot be recovered by independently sampling normal values from printed SEs.

mfrmr currently provides no winner-probability or simultaneous athlete-rank confidence-set API. One-family and experimental two-family EAP intervals condition on point calibration and a specified prior; they do not establish championship-decision accuracy. Multivariate G/D studies can examine dependability of supported score composites and designs, but a high G coefficient is not the probability of selecting the correct champion. Their linear composite weights do not implement the score-dependent trimming of competition panels.

Sport applications have precedents. Looney (1997)¹ applied many-facet Rasch analysis to 1994 Olympic figure skating under the then-used rank aggregation system; that study does not validate current scoring rules or this GMFRM implementation. In gymnastics, Mercier and Heiniger (2019) examined performance-dependent judging error and limitations of ranking-based judge evaluation. Their comparison scores approximate performance quality; agreement with a panel is not independent proof of truth. These studies motivate checks of score meaning, judging variability and decision consequences, rather than automatic exclusion of a judge whose marks change the podium.

Can EM estimate two slope families?

Yes. With fixed task and rater effects and one latent ability, adding a second slope family does not add another latent integration dimension. MML–EM can combine all ratings of each person in the E-step and jointly update task slopes, rater slopes, locations and steps in a constrained numerical M-step. This extends the estimation approach used for GPCM by Muraki (1992); it does not turn the model into JML. Numerical EM can converge slowly or to a local solution, and small likelihood changes alone do not establish convergence. Uncertainty must account for unobserved ability: the curvature of a single M-step is not the observed marginal information matrix.

In the two-family formulation discussed here, fixing the ability distribution to N(0,1) and the geometric mean of task slopes to one leaves the rater slopes free. A task slope of one is the task-family geometric reference. A rater slope of one is a unit multiplier on this model’s identified scale, not the average rater’s slope. Neither is a cutoff for rater competence. Tables and plots must identify both the owning facet and its level, even when task and rater levels use the same labels. The effective slope for a rating is their product, not either component alone.

“Task” and “Rater” describe the example’s roles, not required column names. The same roles could be labelled “Criterion” and “Assessor”. Keep the column names separate from the choice of slope reference and step owner: changing a label should not change the model, whereas assigning steps to a different facet generally does. Neither the order of columns in a data frame nor words such as “rater” should silently select the model’s constraints. The public fit_mfrm() accepts chosen person, facet and score column names. For two families, the order of slope_facet = c(first_owner, second_owner) declares their roles explicitly; it is independent of data-column or facets order.

sirt::rm.facets() provides an existing MML–EM route with item and rater slopes. Its documented slopes multiply ability, with item category intercepts and unscaled rater severity. The Uto–Ueno equation instead multiplies the whole adjacent-category predictor and uses rater steps. They are generally different models, so their estimates should not be compared as alternative optimizers without first matching the response equation and identification. Availability of both fitting routes does not establish software equivalence or qualify uncertainty for mfrmr’s provisional two-family route.

Assignment also matters. Having every task rated by every rater is not enough to distinguish two slope families. For example, with two tasks, two raters, binary scores and a fixed standard-normal ability distribution, the two-family model above has six free location/slope parameters. If each person supplies only one rating, the data provide just four marginal success probabilities, one per task-rater pair. More people with that same assignment cannot resolve all six parameters. Several ratings of the same person provide additional information through their joint response patterns. This does not guarantee precise estimates: category use, sample size, the population assumption and numerical integration still matter. Optimizer convergence alone cannot resolve an insufficient assignment design.

What does changing the ability SD test?

A narrow ability distribution may give different category exposure and less information about a response curve than a broad distribution. But parameter recovery must first be compared on the same identified scale. For the two-owner equation above, write ability as theta = mu + sigma*z, with standard-normal z. The same response probabilities result when task slopes are unchanged, rater slopes are multiplied by sigma, task locations and rater steps are divided by sigma, and rater locations become (location - mu)/sigma. This preserves the task-slope geometric mean of one. A change in raw slope values across these representations is not by itself estimator bias.

A single correctly transformed normal population is different from rater panels assigned different ability means or SDs, or a skewed/mixture population fitted with a normal model. Common rated performances help connect panels, but are not fixed ability anchors unless their values are explicitly supplied. Neither overlap counts nor fixed severity anchors alone establish the ability scale and both slope families. The one-family GPCM and provisional two-family formulation have distinct identification constraints; do not transfer parameter values or anchor settings without matching them.

Do two slopes handle testlets or local dependence?

No. Two fixed slope families allow the response curve to differ across two facets. They still assume independent ratings after conditioning on ability and the specified fixed parameters. Two scores from the same performance may remain related beyond ability: for example, fluency and accuracy ratings may both reflect an unusually successful response to one speaking task. Merely adding task and assessor slopes does not model this shared fluctuation.

Local independence is a conditional assumption. A testlet model adds a latent effect for a block of related ratings and assumes independence after conditioning on ability and that effect. Integrating out the effect leaves dependence among the block’s ratings. This is compatible with one substantive ability, but it introduces additional latent variables for dependence. It should not be described as having only one latent variable or as measuring several substantive abilities automatically.

A testlet effect is itself a random effect; the distinction from random rater severity is which ratings share it, not whether the model uses random effects. Block identity must follow the assessment process. A Person-by-Task effect shared across assessors represents something different from a Person-by-Task-by-Assessor impression shared only across that assessor’s rubric scores. With testlet = "Rater", the current testlet API assigns a separate effect to each Person/rater pair; it does not create one rater effect shared across persons. Repeated ratings and assignment overlap are needed to separate such effects. Nonzero testlet variance is not proof of halo, and zero estimated variance in a weak design is not proof of independence.

The current testlet and random-rater routes are separate RSM models. Neither adds dependence effects to a GPCM fit, and they cannot yet be combined. A future combined model must specify how effects are shared, whether a local effect is multiplied by discrimination, and how all shared effects are integrated. Independently redrawing a common rater for each person would change the model. EM cannot resolve that problem by itself.

For feedback and planning, specify whether a prediction concerns an observed performance/rater or a new one. An observed block’s fitted effect is not the predictive distribution of a new block. Read model-specific residuals and uncertainty alongside the response curves; ordinary infit thresholds do not automatically become calibrated tests in the extended models. Do not sum per-rating marginal information as if dependent ratings were independent, or interpret testlet variance as a G-theory coefficient. A testlet model also does not automatically make differently sized tasks equally weighted.

The Rasch testlet model of Wang and Wilson (2005) and their random-effects facet model provide the local-effect background. Rijmen (2009) explains how testlet and bifactor formulations can be related under particular loading constraints. These relationships do not mean every latent effect or multilevel model is interchangeable. For available operations, consult ?fit_mfrm_testlet, ?fit_mfrm_random_rater and their prediction methods.

The one-ability assumption can be examined before a multidimensional GMFRM is available. The visual diagnostics guide compares residual PCA, Q3-style and marginal screens, model-generated reference checks and explicit alternative-model comparisons. It also distinguishes residual networks/EGA from assignment networks. Current residual screens are exploratory; neither a residual component count nor a network community count is a certified number of abilities.

How does this relate to multivariate G theory?

The two approaches answer different questions about the same assessment. A many-facet GPCM describes category responses in terms of ability, task location, rater severity, steps and discrimination. Multivariate G theory instead asks how observed scores on several criteria, such as fluency and accuracy, vary together across persons, raters and tasks. It can then compare assessment plans, for example adding raters or tasks for a prespecified weighted score. Two slope families still describe one latent ability; they do not create two score dimensions or estimate their G-study covariances.

Start the G-study from numeric score columns and the declared design, not a matrix of GPCM parameter estimates. The current G-study supports one or two common random measurement facets, including its documented nested structure. Balanced ANOVA and identifiable incomplete-design MINQUE(0) analyses are available. The latter does not correct informative missingness. D-study plans use the documented balanced future designs; their help describes the distinct restrictions on point estimates and plan-difference intervals. G/D-study objects have their own summary(), plot() and plot_data() routes and are not mfrm_results() inputs.

For design decisions, begin with the existing two-dimensional D-study plots: condition count on the horizontal axis, G/Phi or SEM on the vertical axis, with separate lines for the other count. Inspect one criterion or declared composite at a time. Multiple score criteria do not require a 3D surface. An increase in G/Phi across populations need not mean smaller measurement error: coefficients depend on person variation as well as error variance. Read SEM in score units alongside them. The plot help explains the distinction between point projections and intervals for prespecified plan differences.

Individual sheets support native additive RSM/PCM and experimental two-family GPCM MML through mfrm_report(res, style = "rater", facet = ..., rater = ...). The two-family sheet separates location, component slope, category use and assignment exposure. Its saved component-slope intervals and same-data posterior residuals retain their own limitations. Missing intervals and unresolved response rows remain visible; location and step intervals are not displayed in this sheet. Separately requested location/contrast intervals belong in the analyst report. The standalone sheet excludes source identifiers and the fitted object. Review it before distribution; the full results RDS remains an analyst artifact. One-family GPCM sheets remain unsupported.

The separate mfrm_generalizability(fit) convenience route does accept an MFRM fit, but uses its stored ratings (or supplied data) for an observed-score main-effects analysis. It does not convert latent model parameters to G-theory components or fit the multivariate model.

There is direct precedent for complementary use. Hirai and Koizumi (2013) used multivariate G theory with mGENOVA to examine criterion-specific dependability across stories and MFRM with Facets to study ratings and scales. Their G-study excluded raters; that design does not by itself establish generalization to new raters. The paper also set negative variance estimates to zero, whereas mfrmr retains raw component estimates and reports when a requested coefficient is unavailable. Using the same broad methods is not a numerical reproduction of that analysis. Jiang et al. (2020) explain estimation of multivariate G-theory parameters through a linear mixed model framework; this is not an automatic conversion of GPCM slopes.

There are also joint IRT/G-theory approaches, including Briggs and Wilson (2007). Such integration requires specifying the population of tasks/raters, the score or latent scale and the quantity to generalize. It is not implemented by passing the sampling covariance of estimated GPCM slopes to a G-study: that covariance describes estimation error, not sources of variation in observed scores. G-study score weights express the intended assessment composite, not estimated GPCM identification or discrimination weights.

For complementary analyses, keep the person, rater and task identifiers and state which ratings and score criteria each analysis uses, including any aggregation or exclusions. Do not silently replace unassigned ratings with predictions or analyze EAP scores as error-free observations. Those choices change the target and its uncertainty. A joint GMFRM/G-theory model with propagated calibration uncertainty and a declared generalization population remains a separate extension; the current APIs do not claim that integration.

What can the related facets literature add?

Wang and Liu (2007) combine item discrimination and item steps with latent regression: person background variables explain the latent ability distribution while response uncertainty remains in the model. The existing MML arguments population_formula and person_data provide a conditional-normal latent regression route in mfrmr. Regressing fitted EAP scores as if they were observed without error is a different procedure. Their model does not itself add a second, multiplicative rater-slope family, and its item steps differ from the rater steps in the Uto–Ueno example.

Wang, Wu and Qiu (2025, preprint) propose a derivative-based rater capability index for binary ratings. The useful practical question is how a rater’s response changes across the ability range. In a GPCM context with effective slope A, expected-score sensitivity is A * Var(Score | Theta, context), whereas information is A^2 * Var(Score | Theta, context). A steep curve and a large amount of information are related but not identical. Existing probability and information curves can help examine this distinction; mfrmr does not currently provide the paper’s normalized capability index.

Population-averaged feedback also depends on the ability distribution and the tasks used for averaging. Comparing raters on different assigned populations can change the comparison; a common reference population answers a different question from each rater’s actual assignment. Such summaries describe a fitted response model, not agreement with an independently established correct score. The preprint’s binary normalization and Laplace approximation should not be treated as validated defaults for polytomous GPCM or sparse rating panels. EM and Laplace address different numerical choices: EM is an optimization strategy, whereas Laplace approximates integration over the latent ability. A Laplace-based marginal likelihood is not ordinary JML, and its accuracy needs checking separately from optimizer convergence.

Which facet receives different slopes?

The slope is not one common number for the entire fitted GPCM. The selected slope_facet contributes one positive slope for each of its levels. If it has \(G\) levels, the returned table contains \(G\) slopes and the model estimates \(G-1\) free relative log-slope contrasts because their geometric mean is fixed to one.

Every facet that is not selected remains an additive location term inside \(\eta\). For example, a criterion-owned GPCM can still estimate rater severity, but it does not give each rater a separate discrimination parameter. The criterion slope nevertheless scales that rater-severity contribution because the contribution is inside \(\eta\). “No separate rater slope” and “rater severity is unscaled” are therefore not the same statement.

Criterion discrimination with rater-specific category use

Suppose several raters score each performance on several rubric criteria. Use step_facet for whose category transitions may differ and slope_facet for whose responses may distinguish ability more sharply. These are distinct questions from overall rater severity, which remains a location effect.

fit <- fit_mfrm(ratings, "Person", c("Rater", "Criterion"), "Score",
  model = "GPCM", method = "MML",
  step_facet = "Rater", slope_facet = "Criterion")
slopes <- confint(fit, parm = "slopes")
grid <- expand.grid(Theta = seq(-2, 2, length.out = 21),
                    Rater = "R01", Criterion = "C01") # use your fitted labels
curves <- mfrm_curve_intervals(fit, grid)
plot(curves, title = NULL, subtitle = NULL)
results <- mfrm_results(fit, intervals = list(slopes = slopes, curves = curves),
                        include = c("fit", "plots"), compute = "never")
saveRDS(results, "criterion-slopes-rater-steps.rds")
mfrm_report(readRDS("criterion-slopes-rater-steps.rds"))

Relative slopes have geometric mean one. Step contrasts are centered within each rater, separately from severity. Adequate crossing and category use are needed to distinguish these effects; a sparse or confounded design can leave intervals unavailable. A large slope or unusual step is not, by itself, a rater-quality or training diagnosis. The numerical workflow checks do not establish finite-sample coverage for this new structure.

For a PCM comparison, keep the same data, rater steps, constraints and population model; the GPCM adds one fewer slope contrast than there are criteria. The default GPCM estimates a population mean and variance, whereas default PCM fixes its population, so defaults alone do not form that matched comparison. Follow the matched population example earlier in this vignette. Existing likelihood and local-information checks still apply. bootstrap_mfrm_gpcm() preserves both owners and the observed assignment for parametric resampling; this is not a new-design simulation. build_weighting_review(), extract_mfrm_sim_spec() and the simulation/design workflows still require a shared owner. Portable GPCM calibration and simultaneous rater/criterion slope families remain unavailable.

An explanation for a high-school reader

Imagine that several teachers grade presentations on several criteria.

• PCM gives everyone rulers with the same magnification. One criterion may be harder and one teacher may be stricter, but a one-unit ability difference changes the rating probabilities at the same basic sharpness everywhere.
• The current mfrmr GPCM lets you choose one box of rulers. If the box is Criterion, each criterion can have a different magnification. If the box is Rater, each rater can have a different magnification. You may choose only one box for magnification. With MML, another box may supply the category steps: for example, criteria can differ in magnification while teachers differ in how they use the score categories. With JML, the same box must supply both magnification and category steps.
• A broader generalized MFRM can give both the task/criterion and the rater their own magnification knobs. In the Uto–Ueno example, the effective sharpness is the product of the task and rater settings. It can also model rater severity and rater-specific use of rating categories.

Thus the current GPCM is a smaller branch inside the broader family: it adds a carefully controlled slope vector to a many-facet PCM, but it is not the full two-knob generalized MFRM. Also, once discrimination is allowed to vary, the model is generalized from the Rasch baseline; it no longer retains the strict equal-discrimination Rasch restriction.

Comparing PCM and GPCM

The package provides three complementary comparison layers:

1. compare_mfrm() returns AIC, Person-based BIC and SABIC comparisons when comparison$table$ICComparable is true. Its local-solution check is independent of slope-interval availability. Add nested = TRUE to test equal relative slopes under the matched PCM/GPCM specifications below; inspect $lrt and $comparison_basis$lrt_reason.
2. build_weighting_review() shows what the GPCM slopes changed: the slope profile, facet-measure shifts, and redistribution of information across levels.
3. build_model_choice_review(..., run_weighting_review = TRUE) bundles the statistical comparison, operational model roles, support boundaries, and the weighting review.

For example, dat is a long-format rating table with Person, Rater, Criterion and Score columns. Fit the same estimated-normal population in both models; ~1 estimates its mean and variance without covariates.

persons <- data.frame(Person = unique(dat$Person))
fit_pcm <- fit_mfrm(
  dat, "Person", c("Rater", "Criterion"), "Score",
  method = "MML", model = "PCM",
  step_facet = "Criterion", quad_points = 61,
  population_formula = ~1, person_data = persons, person_id = "Person"
)

fit_gpcm <- fit_mfrm(
  dat, "Person", c("Rater", "Criterion"), "Score",
  method = "MML", model = "GPCM",
  step_facet = "Criterion", slope_facet = "Criterion",
  quad_points = 61,
  population_formula = ~1, person_data = persons, person_id = "Person"
)

comparison <- compare_mfrm(fit_pcm, fit_gpcm,
                           labels = c("PCM", "GPCM"), nested = TRUE)
comparison$lrt
comparison$comparison_basis$lrt_reason
comparison$table[, c(
  "Label", "LogLik", "AIC", "BIC", "SABIC", "ICComparable"
)]

weighting <- build_weighting_review(fit_pcm, fit_gpcm, nested = TRUE)
summary(weighting)
weighting$comparison_contract

review <- build_model_choice_review(
  PCM = fit_pcm,
  GPCM = fit_gpcm,
  run_weighting_review = TRUE
)
summary(review)
review$model_roles[, c(
  "Model", "StepStructure", "StepCoordinates", "FreeStepParameters",
  "SlopeStructure", "SlopeCoordinates", "FreeSlopeParameters",
  "FitReadiness", "FormalInference", "Interpretation"
)]

StepCoordinates and SlopeCoordinates count the values shown in fitted tables. FreeStepParameters and FreeSlopeParameters count the independent coordinates optimized after identification constraints. For example, with four score categories and four levels of the step facet, PCM and GPCM display 12 step coordinates but optimize 8 free step parameters because each three-step ladder is sum-zero. GPCM additionally displays four relative slopes but optimizes three free log slopes because their geometric mean is one.

The comparison_contract describes the evidence available from the two fits. Its FormalModelSelectionAvailable field follows the verified IC comparison, not numerical convergence alone:

For a JML question, follow that table with recovery conditions rather than an IC or LRT claim. simulate_mfrm_recovery() already distinguishes unit_slopes, near_flat, moderate, and high_dispersion generator conditions. The unit-slope condition asks how often harmless sampling and optimizer variation can create an apparent slope spread; the non-unit conditions ask whether the declared heterogeneity can be recovered. FACETS can support the aligned PCM/JML side of this review, but its Table 7 discrimination remains a post-fit diagnostic and not a free-slope GPCM fit.

PCM is obtained from the aligned GPCM response kernel by setting all slopes to one. Calling compare_mfrm(..., nested = TRUE) records PCM_in_GPCM when the structural restrictions are verified. The result in $lrt is returned only when the numerical and likelihood checks also pass; otherwise read $comparison_basis$lrt_reason. Older saved comparisons must be rebuilt to obtain the new decision. Weighting consequences still need substantive review.

Save the comparison and carry its limitations into reporting

Start with comparison$comparison_basis, then read the criteria and test. For example, two fits can have the same number of Persons and ratings but different response data. In that case same_data and ICComparable are false: raw criteria may remain visible, but criterion differences, weights, preferences and the LRT are withheld. nested = TRUE does not override this decision. A structurally nested response kernel alone is insufficient.

The ordinary PCM/GPCM reporting route uses build_model_choice_review() and export_summary_appendix(). The review retains comparison warnings alongside its tables. It does not request the LRT; retain the separate comparison object above to preserve that explicitly requested test and its reason. The comparison argument to mfrm_results() is for the separate descriptive testlet/shared-rater workflow, not this ordinary IC/LRT result.

tables <- build_summary_table_bundle(review)
tables$tables$comparison_table
tables$tables$comparison_warnings
export_summary_appendix(tables, output_dir = "model-comparison",
                        prefix = "pcm-gpcm", include_html = TRUE)

curves <- plot(fit_pcm, comparison = fit_gpcm, type = "ccc", view = "difference", labels = c("PCM", "GPCM"), draw = FALSE)
curves$data$notes
as_ggplot(curves)

saveRDS(list(comparison = comparison, review = review, curves = curves),
        "pcm-gpcm-comparison.rds")
saved <- readRDS("pcm-gpcm-comparison.rds")
saved$comparison$comparison_basis
as_ggplot(saved$curves)

The saved comparison preserves unrounded values; appendix summaries may round them for presentation. Reloading these objects does not refit either model. Rebuilding compare_mfrm() checks the retained solution again, including its likelihood, gradient and information; it does not rerun the optimizer.

Paired plots show descriptive differences on the fitted coordinates. Their checks concern score coding, constraints, estimator and coordinate basis; they do not require identical response data or establish IC/LRT eligibility. A drawable comparison therefore does not authorize model ranking. By default, titles and notes are omitted from the figure and retained in its returned data. Read data$notes and data$fit_readiness before interpreting the plot; visual differences have no difference intervals or equivalence test.

Why total-score conditioning does not remove ability

For two binary ratings with adjacent log odds \(a_j(\theta-d_j)\), the two patterns \((1,0)\) and \((0,1)\) have the same total, but their log probability ratio is

\[ \log\frac{P(1,0\mid\theta)}{P(0,1\mid\theta)} = (a_1-a_2)\theta-a_1d_1+a_2d_2. \]

The ability term disappears when the slopes are equal. With unknown unequal slopes it remains, so the ordinary total is not sufficient for ability. The analogous polytomous joint likelihood contains \(\theta\sum_j a_j y_j\). If the slopes were known, that weighted total would be sufficient for ability; substituting estimated slopes does not give the RSM/PCM conditional-likelihood argument.

This does not mean that every possible conditional method is impossible. For example, when ratings share slope owners, conditioning on the vector of totals within each owner makes \(\sum_j a_j y_j\) constant and removes ability. It also conditions away more information: if each owner supplies one rating, that vector fixes the complete response and leaves no conditional variation. Identification and usefulness of any remaining comparisons require their own analysis. mfrmr implements total-score CML for additive RSM/PCM, not this different GPCM conditional procedure. Use the supported MML route for GPCM structural inference, retaining its population-model assumptions.

Model structure: which effect does discrimination change?

TAM and ConQuest support flexible slope designs. That does not make every many-facet GPCM in these programs the same model.

The TAM construction is documented in Example 14, Model 14c. ConQuest’s Note 8, pages 1–5 describes its scoring design, population model and identification. It updates the earlier ConQuest 3 note; limitations described only in that older note should not be assigned to current ConQuest.

For a rubric criterion \(c\), rater \(r\) and category transition \(k\), compare these adjacent log odds, using minus signs for difficulty and severity:

\[ \begin{aligned} \text{Current mfrmr:}\quad &\alpha_c(\theta_p-d_c-\rho_r-\tau_{ck}),\\ \text{Additive-intercept design:}\quad &a_c z_p-b_c-u_r-v_{ck}. \end{aligned} \]

The second expression represents the TAM example and a ConQuest design with criterion-shared scores and additive rater intercepts. It is not every model those programs can fit. The comparison is an algebraic consequence of their documented designs: the first expression contains the product \(\alpha_c\rho_r\); independent linear intercept and slope designs do not impose that product of two estimated parameters.

For example, on the mfrmr scale, two raters differing by 0.5 severity units produce an absolute log-odds difference of 0.4 for a criterion with slope 0.8 and 0.625 for one with slope 1.25. In the additive-intercept model, the shared rater contrast has the same log-odds effect at every criterion. These are different assumptions about how rater severity acts. Changing which facet supplies slopes and steps does not by itself change this assumption.

With two slope families and Rater-owned steps, the experimental native predictor is \(a_i c_r(\theta_p-d_i-\rho_r-\tau_{rk})\). A generalized item for each Task–Rater context can represent its probabilities at fixed parameter values. Free estimation must additionally preserve the product structure of the slopes and the shared, slope-dependent intercepts. Giving every context its own free slope and steps removes these restrictions and fits a different model. Consequently, a fixed-calibration probability or scoring check is useful but does not establish free-estimator equivalence.

sirt::rm.facets() provides a useful second MML reference: it estimates item and/or rater slopes under product-one centering and an estimated normal trait spread. Its item-only GPCM and equal-discrimination many-facet reductions are candidate matched comparisons after category, threshold, quadrature, and centering conventions are aligned. Its general free-slope rater kernel is not the current mfrmr kernel: sirt uses the product of item and rater slopes on the trait term while rater severity remains a separate location term, whereas mfrmr multiplies the complete adjacent-category predictor in its one-family route; the experimental two-family route uses the same full action with a product slope. Such results are different-model sensitivity evidence, not automatic numerical validation. immer provides PCM-design JML/CML/CCML routes and a hierarchical rater model, not a direct free-GPCM many-facet match. Its role is consequently a unit-slope reduction or alternative-model sensitivity analysis rather than free-slope numerical validation.

Exact overlap and what has been checked

With items only, free item-specific steps can absorb the slope multiplication. For an estimated normal ability population, write \(\theta=\mu+X\beta+\sigma\epsilon\), with standard-normal \(\epsilon\), and \(\prod_i\alpha_i=1\). Set \(z=(\theta-\mu)/\sigma\), \(a_i=\sigma\alpha_i\), \(b_i=\alpha_i(d_i-\mu)\), \(v_{ik}=\alpha_i\tau_{ik}\) and \(\beta_z=\beta/\sigma\). Then \(\alpha_i(\theta-d_i-\tau_{ik})=a_i z-b_i-v_{ik}\). Here \(\sigma\) is the residual population SD given any covariates \(X\); without covariates it is the population SD. This derives an identical response model on a unit residual-variance scale; the intercept/step signs and centering must still match each program. Fixing both population variance and the geometric mean of slopes to one instead imposes an additional restriction, so it is not the same conversion. The overlap is for positive slopes; ConQuest also exposes a positivescores setting, so the permitted parameter space must be checked.

A retained validation comparison exercises the exact item-only TAM overlap on TAM’s data.gpcm fixture at q=31 and q=41. TAM fixes the latent variance and estimates absolute slopes; mfrmr fixes the geometric mean of its relative slopes and estimates the latent variance. After applying the analytic scale and location transformation, maximum q-specific differences were about 1.8e-5 for relative slopes, 1.3e-5 for transition thresholds, 6e-6 for fitted category probabilities, and 1.4e-6 for deviance. The transformed TAM probability formula itself agreed to binary64 precision. This is external evidence for the item-only GPCM kernel, coordinate map, and continuous MML target; it is not a comparison of the rater-plus-slope model, a cross-engine SE validation, or permission to override mfrmr readiness. These results use TAM 4.3-25. A separate q=31 many-facet re-expression was retained as a different-model comparison after the algebraic distinction above was found.

An earlier local ConQuest 5.47.5 MML run checked the item-only overlap with one population covariate, 120 Persons, five items and four categories. The deviances differed by about 1.3e-5 after matched identification. Its exact coordinate transformation reproduces probabilities within 1e-14. This is historical, version-specific evidence for that model and mapping, not a fresh check of the current source, interval coverage or a multifacet slope model. The public ConQuest overlap bundle does not automate this GPCM comparison.

Inference: same estimates do not imply the same intervals

The saved TAM 4.3-25 fits retain marginal slope and intercept SEs, but not the joint slope covariance or slope–intercept cross-covariance required by the nonlinear transformation to mfrmr’s relative slopes, thresholds, and free population scale. Two positive- definite covariance matrices can preserve exactly those reported marginal SEs while producing different transformed SEs. This non-identification means the comparison must be withheld rather than completed by assuming zero covariance. TAM’s documented tam.se() also omits parameter covariances and labels loading SEs for GPCM.design experimental. This describes that route and the saved objects; it does not rule out obtaining a joint matrix through additional methods.

ConQuest’s command reference provides estimatecovariances with stderr=empirical and matrixout, and show table 10 reports parameter-estimate error covariances. In contrast, export covariance exports the latent-population covariance. It is not the sampling covariance needed for transformed calibration SEs. Parameter ordering, constraints and the inclusion of free score parameters must be checked in the actual exported matrix before using it. No such current-source joint-covariance comparison is claimed here. ConQuest’s estimated-score route is an MML comparator; its JML route does not estimate free scores.

For example, let \(s=\log\sigma\) and \(\ell_i=\log\alpha_i\). A population-SD-standardized slope has \(\log a_i=s+\ell_i\), so its delta-method variance requires

\[ \operatorname{Var}(\log a_i)=\operatorname{Var}(s)+ \operatorname{Var}(\ell_i)+2\operatorname{Cov}(s,\ell_i). \]

Multiplying the current relative-slope interval endpoints by the estimated \(\sigma\) omits two of these terms. Use confint(fit, parm = "slopes", scale = "standardized") to include all three; the default still reports relative slopes. Conversely, converting absolute slopes to geometric-mean-one slopes requires their joint covariance.

TAM’s anova() method reports information criteria and likelihood-ratio comparisons. These methods and mfrmr’s current IC/LRT operations still require comparable likelihoods. Replicating another program’s fit requires the same model. IC selection can instead compare different, nonnested models when their likelihoods refer to the same observed data with compatible likelihood definitions and constants. Verify category coding, population specification, each model’s free-parameter count and integration accuracy. An LRT additionally requires a nested null; different many-facet slope actions are not automatically nested. The simple PCM/GPCM null with \(G-1\) relative-slope restrictions does not establish a reference distribution for arbitrary cross-program design changes.

These distinctions also govern diagnostics. FACETS’ discrimination, mfrmr residual PCA or bias screens, and hierarchical-rater variance components answer different questions. Agreement among them can be supporting evidence; a difference is not a software error until the fitted likelihood, estimand, identification, data rows, and output transformation have first been shown to match.

A population-level validation also asked whether the current complete-predictor slope and a loading-only slope are merely two coordinate systems. They reduce exactly when slopes are all one or when the non-owner Rater contrast is zero. In a fully crossed four-Criterion/four-Rater design, however, their identifying adjacent-logit difference in differences is zero for loading-only and \(-\{\alpha_c-\alpha_d\}\{\rho_r-\rho_s\}\) for the current kernel. After all available slopes, severities, and transition boundaries were re-estimated, moderate and strong fixed conditions retained population KL losses of about 0.0016 and 0.0085–0.0088 per response, with maximum common-grid probability differences of about 0.047–0.184. q=31 and q=41 conclusions agreed closely. These are fixed model-form demonstrations, not practical cutoffs or finite-sample operating characteristics. They support treating a loading-only kernel as a separate possible future family rather than as an interpretation switch for an existing GPCM fit.

Sparsity changes this distinction structurally, not merely by reducing sample size. On an observed Rater-by-Criterion graph with no cycles, the products \(\alpha_c\rho_r\) can be represented exactly by new additive Rater and Criterion coordinates. Thus even a connected tree cannot distinguish the two slope actions. A known-ability validation compared four connected graphs. At 250 Persons, truth-selection rates were about 0.96 under complete crossing, 0.76–0.78 for a balanced eight-edge cycle, and only 0.58 for an eight-edge cycle localized to two Raters and two Criteria; the connected tree remained observationally equivalent. Each Person was observed on every edge, and abilities and both parameter vectors were supplied to the comparison, so these are optimistic fixed-parameter results rather than JML or MML model- selection performance. The design lesson is to inspect which slope and severity contrasts are traversed by cycles, not to infer identifiability from connectedness, density, or an anchor percentage alone.

A subsequent finite-sample validation removed the known-ability advantage and refitted both slope actions by direct MML. The common identification estimated geometric-mean-one relative slopes, sum-zero Rater severities, mean-zero transition boundaries, and a normal population mean and standard deviation. Across 12 replications per truth and design at 250 training Persons, all 96 candidate fits converged and all full 19-coordinate Hessians were positive definite. Two starts differed by at most about 7.9e-6 log-likelihood units, and q=31 versus q=41 changed none of the 48 paired family decisions. Complete crossing selected the true family in every training and independent holdout dataset. Under the balanced cycle, training selected the truth in 7–9 of 12 datasets and median relative-slope RMSE increased from about 0.06 to as much as 0.12. A stricter retained-solution polish reduced the largest gradient by roughly an order of magnitude without materially changing likelihoods. These small fixed counts clarify optimizer, curvature, and design contributions; they do not calibrate a universal selection, standard-error, or readiness rule.