Dependent ratings within a person: testlet models

Declare what shares an effect

In this example, Rater has two distinct roles. Its fixed effect describes that rater’s severity across persons. Its testlet membership groups that rater’s criteria within each person. The local effect for Person P01 and Rater R01 is independent of the local effect for Person P02 and Rater R01. It is not one random rater effect shared across both persons. For the latter model, see vignette("mfrmr-random-raters", package = "mfrmr").

This example corresponds most directly to the distinction between fixed rater severity and Person-by-rater local effects in Wang and Wilson’s random-effects facet model (2005), equations 14–15. Their Rasch testlet model provides the broader block-dependence formulation. The implementation here uses one common local variance; it does not estimate a separate inconsistency parameter for each rater.

ratings <- load_mfrmr_data("example_core")
head(ratings[c("Person", "Rater", "Criterion", "Score")], 8)

Every row belongs to exactly one testlet. A separate Session or Station column can be used instead of Rater; it need not also be a fixed facet. Labels may be reused across persons. If two sittings within the same person have independent local effects, give them distinct testlet labels. The model does not infer membership or support overlapping groups.

Fit and check the calibration

The following call fits the model. To avoid repeating this calculation whenever the tutorial is rendered, the displayed results use the packaged synthetic example fitted with these data and settings. Run the call yourself to estimate a new calibration; reading the saved example is not a new estimation check.

testlet_fit <- fit_mfrm_testlet(
  ratings, person = "Person", score = "Score", testlet = "Rater",
  facets = c("Rater", "Criterion"), score_levels = 1:4,
  quad_points = 121
)

saved_example <- readRDS(system.file("examples", "extended-models.rds",
                                      package = "mfrmr", mustWork = TRUE))$testlet
testlet_fit <- saved_example$fit
testlet_fit$checks
summary(testlet_fit)$data_usage
testlet_fit$calibration$variance
testlet_fit$calibration$person_variance
stopifnot(testlet_fit$checks$NumericalReady,
          testlet_fit$checks$InformationPositive)

The model uses adjacent-category RSM probabilities, a mean-zero normal ability population with estimated variance, and independent zero-mean normal testlet effects with a common estimated variance. Fixed-facet effects sum to zero; free category steps set the overall location. Estimation is frequentist marginal maximum likelihood, with both ability and local effects integrated using Gaussian quadrature. There are no calibration priors and no additional fitting dependency.

person_sd = NULL, the default, estimates ability variance as well as local variance. Mean-zero ability and free step location avoid estimating two interchangeable locations. Supply person_sd = 1 to fit the earlier N(0,1) model, or another positive SD only when treating it as known is justified. With a unit Rasch slope, fixing ability variance is a population restriction, not merely choosing units. Estimating it still assumes a common normal population independent of assignment; it does not repair skewness or different ability populations across rater panels.

After integrating ability and local effects, different persons are independent under this model. Their individual ratings are not independent replications. More persons can provide information about the common variances when the design separates ability from local dependence; adding criteria in one block is not equivalent. Normal calibration intervals require an identified interior model with adequate information. Their information matrix includes both estimated variances, rather than treating fitted variances as known.

NumericalReady checks optimization, search boundaries, and agreement with a higher integration order. InformationPositive checks usable local curvature. The default order is 31; this example uses 121. Inspect the individual differences when a check fails; increasing quad_points can resolve integration errors, but cannot repair an unidentified or unsuitable model. No order is universally sufficient. Failed checks withhold intervals and block Person scoring. Passing checks does not establish model fit or frequentist coverage.

The variance can be exactly zero. Such a result does not prove that the assessment has no local dependence. Ordinary calibration intervals are withheld at an estimated zero variance; a regular variance interval is not provided. testlet_variance = 0 is a different decision: it specifies zero as known and estimates under that restricted model. A known positive variance can also be supplied; it should not be treated as known merely because another analysis estimated it.

Ability variance can also be estimated at zero. Both the population estimate and its boundary status remain visible, but predict() then returns all requested Persons as unavailable, with an explanation. It does not report zero-width ability intervals. variance_max and person_variance_max bound the searches for local and ability variance, respectively; hitting either upper bound requires numerical review.

For comparison with ordinary MFRM, preserve the observed events, category coding, fixed Rater and Criterion effects, and the ability population and location convention. With person_sd = 1, use the fixed N(0,1) ordinary model; with estimated ability variance, the ordinary model must also estimate its population and its location must be aligned. A known zero local variance gives the corresponding ordinary RSM; an estimated zero is a variance-boundary result. Standard chi-squared likelihood-ratio calibration is not automatic. compare_mfrm() checks these conditions and compares centered fixed-facet effects. Separately supplied source-roster scores add Person comparisons; no preferred model is chosen. Changes in scores or interval widths alone do not demonstrate better accuracy when the true abilities are unknown.

Review fixed-rater feedback

subset(testlet_fit$calibration_table, Facet == "Rater")
# Bounds require an explicit choice of approximation and nominal level.
confint(testlet_fit, parm = "calibration", level = .95)
plot(testlet_fit, facet = "Rater", intervals = "normal", level = .95)

Positive severity means lower expected scores at the same ability and other facet levels. Default tables retain estimates and approximate SEs with missing bounds; default plots show points only. The explicitly requested whiskers above are 95% observed-information normal approximations for the fixed effects. Their finite-sample coverage remains unresolved in the estimated-population comparison below. They are not intervals for a newly sampled rater, tests of rater quality, or simultaneous pairwise comparisons. Missing whiskers remain missing. plot(testlet_fit, facet = "Criterion") reviews criterion effects.

Score persons and interpret their intervals

For people already in the fit, use score_mfrm_persons(testlet_fit) to reuse the source ratings. This is the common entry for supported ordinary, shared-rater and testlet RSMs. The example below uses the testlet-specific predict() method to supply an explicit rating table. For a shared-rater model, predict() instead needs supplied abilities and returns response probabilities.

first <- ratings[ratings$Person %in% unique(ratings$Person)[1:4], ]
scores <- predict(testlet_fit, first)
scores$table
plot(scores)
stopifnot(all(scores$table$Status == "available_conditional"))

Estimate is the posterior mean (EAP). ConditionalSD describes the Person posterior spread. Lower and Upper delimit a continuous equal-tail posterior interval. The plot therefore labels them conditional intervals: all calibration values, including ability and local variances, are held fixed. These intervals exclude calibration-estimation uncertainty. They are neither calibration-adjusted confidence intervals nor tests of differences between persons. A 95% conditional interval does not guarantee 95% coverage at every fixed ability.

A simulation compared ordinary and testlet RSMs on the same observed ratings, with both models estimating normal ability variance. There were 120 independent datasets per condition, with 12 prespecified Persons scored from each complete calibration sample. True ability SD was 1.3; positive local variance was .8. The balanced design had two tasks with three criteria each.

The table reports marginal coverage, averaging over generated abilities and ratings, after correcting a numerical start-selection issue. Bracketed ranges are 95% Monte Carlo intervals for the simulation result, calculated with datasets as the independent units. They are not Person intervals.

All 120 panels per condition were available in the corrected replay. In the original run, 115, 113, 119 and 119 panels were available: fourteen cases selected an abnormally terminated optimization start even though a converged start had the same likelihood within floating-point precision. All fourteen passed after the selection repair, with unchanged numerical tolerances. The selection rule leaves 466 of the 480 fits unchanged; four of these successful controls were refitted and retained identical estimates and scores. This replay reuses the generated data; it is not an independent confirmation study. The original outcomes remain recorded. Fitting used 61 quadrature points, with a prespecified 121-point refit for unresolved grid agreement; the results do not establish that the default order of 31 suffices.

Accounting for local dependence improved interval coverage relative to the ordinary model in these three positive-dependence conditions. Point-score accuracy did not consistently improve. Mean squared error measures distance from the generated abilities; smaller is better. On the matched datasets, testlet minus ordinary EAP mean squared error was .00665 [95% MC interval .00190, .01140] with 120 Persons and .05174 [.02582, .07766] with 24 Persons in the balanced design. Thus point errors were larger in those two conditions even though coverage improved. Differences in the null and sparse conditions were inconclusive. Mean conditional interval widths increased by .78, .76 and .67 logits in the three positive-dependence conditions; a wider interval alone does not establish better scoring.

The known-calibration reference had coverage of 94.2–95.7%. In the balanced positive-dependence conditions, fitted-calibration testlet coverage was 1.25 percentage points lower [95% MC interval -2.02, -.48] with 120 Persons and 2.50 points lower [-4.35, -.65] with 24 Persons. This compares all calibration parameters together; it does not isolate the effect of estimating ability variance. The small-sample testlet coverage MCSE was 1.01 percentage points. There was no increase in replication after seeing the outcomes.

The sparse condition retained a common three-criterion task for every Person, plus either two criteria on a second task or one on a third. It changes workload and block length together; it does not test weak/disconnected linking or informative assignment. Unassigned ratings were absent, not imputed. These results concern same-source conditional scores, not new-Person transport, calibration-adjusted intervals, or 95% coverage at every fixed ability. The small and sparse conditions did not establish nominal 95% marginal coverage. Earlier fixed-N(0,1) results remain specific to that population assumption. No sample-size threshold or automatic choice of model follows.

For each person, local-effect integration is checked at two orders, and the interval endpoints invert a continuous ability CDF. A numerical failure retains an unavailable row with a Reason; it does not silently remove the person or replace the score. Use summary(scores)$status_counts to inspect all statuses. Successful numerical integration alone does not validate the ability prior or the local-dependence assumptions.

Missing assigned scores and sparse designs

An unassigned rating is absent from the table. An assigned rating whose score was not recorded has an explicit NA. The default rejects missing assigned scores. Use missing = "omit" only after deciding that analyzing the observed ratings is appropriate; it records omissions without imputing scores or correcting informative nonresponse.

incomplete <- first
last_person <- tail(unique(incomplete$Person), 1)
incomplete$Score[incomplete$Person == last_person] <- NA_real_
incomplete$Score[1] <- NA_real_
review <- predict(testlet_fit, incomplete, missing = "omit")
review$table
review$data_usage
plot(review)
stopifnot(sum(review$table$Status == "prior_only") == 1)

The entirely missing person has status prior_only, mean zero and the SD stored in testlet_fit$calibration$person_sd. These are population values, not measured average performance. The plot uses an open circle for this case. The partially missing person is scored from the remaining ratings, with the omission count retained.

Unequal testlet sizes and absent testlets are allowed. For fitting, at least two persons must have multiple observed testlets, and at least two must have repeated ratings within a testlet. Fixed-facet levels must be observed and the fixed design must have full rank. These are initial requirements, not a proof of estimability. Sparse assignments can still provide little information about the local variance. Ability, local effects and assignment are assumed independent; choosing assignments based on unmodeled ability can violate that assumption even without missing scores. Omission and sparse assignment are therefore different issues.

Reuse and save the result

predict(testlet_fit) reuses the source assignment and its omission policy. For new data, new Person and testlet labels are allowed; fixed-facet levels must already exist in the calibration. Supply all ratings to use for each Person in one call. The function neither appends to cached responses nor conditions on a saved local-effect estimate. To score additional ratings for an existing person, supply that person’s complete set with the correct shared memberships. Do not duplicate the old rows when assembling it.

The fixed-facet restriction still applies when a column has two roles. Here Rater is both the testlet label and a fixed facet. A new Person rated by the known raters can be scored, but a new rater label has no fitted fixed severity and is refused. A separate testlet column that is not a fixed facet may use new labels. Changing labels is not a substitute for choosing which ratings actually share a local effect.

path <- tempfile(fileext = ".rds")
saveRDS(testlet_fit, path)
restored <- readRDS(path)
stopifnot(identical(restored$parameters, testlet_fit$parameters))
saved_plot <- plot(scores, draw = FALSE)
plot_data(saved_plot)$table

Connect to the common graphics and table workflow

The ordinary fixed-facet model, shared random raters and Person-specific testlets answer different questions. Use mfrmr_output_guide("models") to choose the route, then consult help("mfrmr_workflow_methods") and mfrmr_interval_guide() for its supported outputs. In particular, predict() for a shared-rater fit returns score probabilities at supplied abilities, whereas predict() for a testlet fit estimates Person abilities from ratings.

The shared plot accessors work directly on the testlet results. They return the stored results without drawing or refitting:

plot_data_components(scores)
score_table <- plot_data(scores, component = "table")
plot_data(testlet_fit, component = "table", facet = "Rater")
stopifnot(identical(score_table, scores$table))

Optional ggplot2 conversion retains all labeled rows, prior-only symbols, missing whiskers and the interpretation note. It does not change the model or recalculate an interval. The incomplete-data result makes those distinctions visible:

if (requireNamespace("ggplot2", quietly = TRUE)) {
  figure <- as_ggplot(review)
  print(figure + ggplot2::theme_minimal(base_size = 12))
  # For a file: ggplot2::ggsave("person-scores.png", figure,
  #                            width = 7, height = 5, dpi = 300)
}

as_ggplot(testlet_fit, facet = "Rater") uses the same route for fixed-rater effects. For another question, plot(scores, style = "precision") shows estimates against conditional interval width; style = "distribution" shows the empirical cumulative distribution of response-based point estimates, excluding prior-only and unavailable scores. This is not the latent population distribution. Shrinkage and Person selection affect it.

The extended-model plots share sort = "estimate" or "uncertainty", palette = "mono", show_labels, reference, text_scale and point_size. Replace title or caption, use "" to remove either, or set show_title = FALSE, show_notes = FALSE when preparing a figure with a separate caption. These controls also pass through plot(res, type = "scores", ...). Explain the conditional interval target in that caption or the surrounding text. plot_data(scores)$alt_text and $display_data provide a text alternative and row inclusion/reasons; supply them alongside images for accessible output. The shared-rater tutorial shows these views and presentation controls with saved estimates. plot_data(scores)$settings and plot_data(scores)$notes retain the basis for custom displays. Plot data saved by earlier versions can be converted without refitting; recreating a plot from the saved fit or scores also refreshes its display metadata.

For table export, keep Status and Reason rather than selecting only finite scores. An RDS object preserves settings that a CSV table alone cannot:

csv_path <- tempfile(fileext = ".csv")
utils::write.csv(review$table, csv_path, row.names = FALSE, na = "NA")
score_path <- tempfile(fileext = ".rds")
saveRDS(review, score_path)
stopifnot(identical(readRDS(score_path), review))

Compare fixed effects with an ordinary RSM

For the same rating events, compare the usual rater and criterion effects with an ordinary model that estimates a common normal ability SD. This asks how modeling local dependence changes the fitted fixed effects; it does not test whether the local variance is zero or select a preferred model.

ordinary_fit <- fit_mfrm(
  ratings, person = "Person", facets = c("Rater", "Criterion"), score = "Score",
  model = "RSM", method = "MML", rating_min = 1, rating_max = 4,
  keep_original = TRUE, population_formula = ~1,
  person_data = unique(ratings["Person"]), quad_points = 121
)

model_comparison <- compare_mfrm(
  ordinary_fit, testlet_fit, labels = c("Ordinary RSM", "Testlet RSM")
)
model_comparison$checks
model_comparison$effects
plot(model_comparison, facet = "Rater", style = "difference")

The comparison retains both models’ population and numerical-check metadata. Within each facet it subtracts the unweighted mean over the same complete set of levels; the raw values, centers and differences remain available. Missing level estimates or failed numerical checks remain unavailable. No difference SEs, limits of agreement, automatic AIC/BIC preference or chi-squared LRT are computed. The equality or zero line is an orientation aid. Person scores, step locations and predictive Infit/Outfit are not compared by this facet output. Use style = "paired" for an equality-line display, or as_ggplot() for further styling; title/caption and accessibility controls are shared.

Observed events, category coding and omitted identities must agree. New ordinary fits retain prep$omitted_data; old fits that omitted scores without this provenance need refitting from the complete assigned roster. Do not match datasets by counts alone. Zero known testlet variance retains the fixed rater effects, unlike zero shared-rater variance, which removes rater differences.

Inspect posterior predictive residuals

Do some raters have more residual variability under the testlet model? Compute the predicted category distribution by integrating both ability and the Person-specific local effect. Its variance includes uncertainty in the conditional mean, not just the conditional rating variance. The resulting Infit and Outfit are descriptive summaries of the same data used for fitting. An expected value of one and ordinary-model cutoffs are not established; a lower value is not evidence that a rater should be removed.

response_review <- mfrm_response_diagnostics(testlet_fit, group_by = "Rater")
response_review$measures
plot(response_review, style = "paired")

Use style = "scatter" to compare the two indices on separate axes. palette = "mono", show_title = FALSE, show_notes = FALSE, show_labels = FALSE and custom title/caption control the presentation. as_ggplot(response_review, style = "scatter") preserves the definition and text alternative. plot_data() retains unavailable groups and their reasons.

response_review$rows contains original row numbers, predicted means, predictive variances and residuals. response_review$probabilities contains category probabilities on the declared score scale. Assigned missing scores have no residual and are not imputed; absent assignments are never added. Groups with any unresolved observed row have unavailable summaries. rows = ... selects original rows for output, while the full observed roster still informs the posterior. This is not held-out validation or a formal test; ordinary plug-in Infit/Outfit cannot be compared directly with these values.

Compare predictions on the same definition

Does accounting for Person-local dependence change expected scores and residual summaries? Calculate the ordinary RSM’s posterior predictive quantities, then compare them with the saved extension results. Both models use the same observed ratings and selected events. The ordinary plug-in indices returned by diagnose_mfrm() have a different definition and cannot replace this step.

ordinary_response_review <- mfrm_response_diagnostics(ordinary_fit,
  group_by = "Rater")
model_comparison <- compare_mfrm(ordinary_fit, testlet_fit,
  labels = c("Ordinary RSM", "Testlet RSM"),
  response_diagnostics = list(ordinary_response_review, response_review))
model_comparison$responses$measures
plot(model_comparison, metric = "infit", style = "difference")

These differences are descriptive and hold each model’s calibration fixed. The zero line means the summaries agree. A smaller Infit does not show that the model is better or the rater is more accurate: these predictions reuse the observations being checked, and their mean squares have no established reference value of one.

as_ggplot(model_comparison, metric = "probability", show_labels = FALSE,
  palette = "mono", show_title = FALSE, show_notes = FALSE)

Use metric = "expected_score", "variance" or "outfit" to change the quantity, and category = 4 with metric = "probability" to focus on one score. Expected scores remain on the original scale; variances are in squared score units. responses$rows retains both original row numbers and responses$events identifies each matched rating. Repeated events retain their multiplicity. If one prediction is unavailable, its difference remains missing with a reason. Group summaries never silently discard an unresolved observed row. Row numbers can differ between fits; matching uses event contents.

The report below saves these comparisons and supports plot(res, type = "response_comparison", metric = "infit"). To collect the ordinary model’s summaries separately without running its plug-in diagnostics, use mfrm_results(ordinary_fit, response_diagnostics = ordinary_response_review, compute = "never"). Plotting, collecting and exporting saved comparisons performs no new integration.

Compare Person scores and inspect model-aware maps

How do individual scores change when the model accounts for local dependence? Use conditional scores computed from the same complete source roster. Selecting returned Persons does not remove the other observed responses. The ordinary model’s scoring retains its fitted normal mean and variance.

source_scores <- score_mfrm_persons(testlet_fit, persons = scores$table$Person)
ordinary_person_scores <- score_mfrm_persons(ordinary_fit,
  persons = source_scores$table$Person)
model_comparison <- compare_mfrm(ordinary_fit, testlet_fit,
  response_diagnostics = list(ordinary_response_review, response_review),
  person_scores = list(ordinary_person_scores, source_scores))
model_comparison$persons$table
plot(model_comparison, metric = "person", style = "difference", show_labels = TRUE)

The comparison subtracts each model’s fitted population mean to remove its arbitrary origin; it preserves the Rasch logit unit. It retains raw scores, origins, posterior SDs and conditional endpoints. Shrinkage and population variance can still differ. Separate score intervals do not constitute an interval for their difference, a ranking of models or a test between Persons. Prior-only and unavailable differences remain withheld.

map_results <- mfrm_results(testlet_fit, scores = source_scores,
  response_diagnostics = response_review, comparison = model_comparison)
plot(map_results, type = "wright", show_labels = TRUE)

The Person column shows the selected source EAPs. Facet positions equal the severity coefficient plus the mean step location: the average adjacent category boundary when other effects are zero. The category column shows each step at that zero reference. The Person-local testlet effect is set to zero for these reference locations; it has no global facet column. These positions are conditional references, not posterior-averaged category crossings. A higher reference location means more ability is needed to cross the corresponding boundary. Adding several plotted facet locations together would count the mean step more than once; use the model equation for a specific rating profile.

Whiskers show only Person intervals conditional on calibration. Translating an ordinary facet interval by the step mean would omit the covariance and uncertainty of that mean, so no such composite interval is drawn. The separate calibration/rater plots retain their original uncertainty targets.

as_ggplot(plot(map_results, type = "fit_pathway", facet = "Rater",
  palette = "mono", show_title = FALSE, draw = FALSE))

This pathway connects position with the previously defined descriptive posterior residual index. Here the summaries include every observed source rating. Use fit_stat = "Outfit" to change the horizontal index. A large or small value is a prompt to inspect the selected data, not a rater-quality classification. No ordinary fit cutoff, ZSTD or automatic exclusion is supplied.

Use facet, persons, show_steps, show_intervals and show_labels to select the display. show_title = FALSE and show_notes = FALSE hide annotations while retaining their meaning in plot_data(). Shapes accompany colour; palette = "mono" supports grayscale output. Tables retain unavailable rows and original/transformed coordinates. Older saved testlet scores lacking a complete scoring_data roster must be regenerated from the saved fit before these maps or score comparisons; no refit is needed. See ?mfrmr_model_maps.

Build a report from saved results

The common reporting route accepts the fit and separately computed scores. It never runs predict(), refits the model or adds ordinary-model diagnostics. The scoring table retains every requested Person, including prior-only and unavailable rows with their reasons.

res <- mfrm_results(testlet_fit, scores = source_scores, comparison = model_comparison,
                    diagnostics = response_review)
res$tables$interval_basis
res$tables$scoring_status
report <- mfrm_report(res)
summary(report, view = "reader")
archive_dir <- tempfile("testlet-report-")
archive <- export_mfrm_results(res, archive_dir, preset = "starter",
  acknowledge_sensitive = TRUE) # Synthetic example; real exports retain IDs
stopifnot(nrow(archive$plot_errors) == 0)
restored_results <- readRDS(file.path(archive_dir, "mfrmr_results_results.rds"))
stopifnot(identical(restored_results$tables$person_scores, source_scores$table))

Open index.html in archive_dir for the model-specific figures and report. CSV files keep all rows; RDS preserves the fit, scores and their settings. Run the replay script from that folder to reload the report without fitting or scoring. Scores made before prediction-source metadata was introduced must be regenerated with predict() from their saved fit before attachment; refitting is unnecessary. Matching checks prevent mixing different column roles, facet levels, calibrations or model settings.

Saved fits can score in a new R session without a live optimizer. Scores and plot data can also be saved separately. The result is a distinct mfrm_testlet object: ordinary mfrm_fit diagnostics, response-MI pooling, portable calibration objects and random-rater methods do not accept it. Numerical checks are not model-fit diagnostics. The Shiny viewer, ordinary Wright map and ordinary plug-in Infit/Outfit diagnostics are unavailable for this class. The separate posterior predictive summaries above do not supply formal fit tests or ordinary-model thresholds. An earlier ordinary model needs a new fit with explicit testlet membership; changing its saved class does not convert the model.

Earlier saved testlet fits retain their fixed N(0,1) population when rescored or reported. A new fit estimates ability variance by default; specify person_sd = 1 to reproduce the earlier model. Keep a new fit and its newly computed scores together when changing this assumption; do not attach scores from a different population fit to a report.

The supported scope is one ability, RSM, unit weights, additive fixed facets and one common variance for independent Person-specific testlets. PCM, heterogeneous or correlated local effects, shared random raters and joint multidimensional abilities need other models. The testlet modeling idea is described by Wang and Wilson (2005), The Rasch Testlet Model.