Troubleshooting reference-model workflows

This article starts where Getting started with referent (vignette("getting-started", package = "referent")) ends. You have a fitted reference model, but a score is missing, a transported site behaves unexpectedly, or a longitudinal extension refuses to make a claim. Each section follows the same path: recognize the symptom, identify the cause, and make the smallest defensible fix.

library(referent)

reference <- ref_simulate(240, seed = 10)
held_out <- ref_simulate(100, seed = 11)
spec <- ref_spec(ref_gaussian(), location = ~ s(age, k = 5) + sex)
fit <- ref_fit(spec, data = reference, outcomes = "y")

Why is a score missing for an ordinary-looking row?

Symptom

A used covariate is missing, or the row lies outside the range represented by the reference data. referent keeps the fitted median where it can, but masks the probability-scale quantities that would otherwise look more trustworthy than they are.

in_rows <- which(ref_support(fit, held_out)$support == "in")[1:4]
incomplete <- held_out[in_rows, ]
incomplete$age[2] <- NA_real_

missing_scores <- predict(fit, incomplete, uncertainty = "conditional")
missing_scores[, c(".row", "z", "support", "status")]
#> # A tibble: 4 × 4
#>    .row      z support status
#>   <int>  <dbl> <chr>   <chr>
#> 1     1 -1.29  in      ok
#> 2     2 NA     unknown missing_predictor
#> 3     3 -0.969 in      ok
#> 4     4  0.509 in      ok

An out-of-range covariate is different from a missing one, but it receives the same conservative treatment on the score scale.

outside <- held_out[in_rows[1:3], ]
outside$age[1] <- max(reference$age) + 20

outside_scores <- predict(fit, outside, uncertainty = "conditional")
outside_scores[, c(".row", "median", "z", "support")]
#> # A tibble: 3 × 4
#>    .row median      z support
#>   <int>  <dbl>  <dbl> <chr>
#> 1     1  13.3  NA     out
#> 2     2  10.4   0.416 in
#> 3     3   7.73 -0.969 in

Cause

The formulas in ref_spec() define the covariates the model actually needs. An NA in one of those columns gives status == "missing_predictor" and support == "unknown". A finite value beyond the observed reference range gives support == "out". Neither condition is repaired by imputing or silently extending the fitted curve.

Fix

Supply the missing predictor from a defensible data source, or restrict the analysis to the declared reference domain. Use ref_support(fit, newdata) before a large scoring run to audit support without fitting or scoring again. Do not use allow_extrapolation = TRUE merely to turn an unsupported age into a number.

Why is a transported site’s score still marked new_group?

Symptom

The shared model was fitted on sites A-D. Site E has local reference controls, and ref_adapt() estimates a site-E offset, but the original support record still says that E was absent from the training reference.

site_spec <- ref_spec(
  ref_gaussian(),
  location = ~ s(age, k = 5) + sex + s(site, bs = "re")
)
site_fit <- ref_fit(site_spec, reference, outcomes = "y")

site_e <- ref_simulate(
  600, sites = 5,
  site_shift = c(0, 0, 0, 0, 1),
  seed = 12
)
site_e <- site_e[site_e$site == "E", ]
local_controls <- site_e[1:20, ]
site_e_target <- site_e[21:50, ]
adapted <- ref_adapt(site_fit, local_controls, by = site)

By default, unseen-group probability columns remain masked. After fitting an explicit transport model from local controls, you can deliberately request the transported scores while retaining support == "new_group" as provenance.

masked <- predict(adapted, site_e_target, uncertainty = "conditional")
transported <- predict(
  adapted,
  site_e_target,
  uncertainty = "conditional",
  allow_extrapolation = TRUE
)

data.frame(
  default_z = head(masked$z, 3),
  transported_z = head(transported$z, 3),
  support = head(transported$support, 3)
)
#>   default_z transported_z   support
#> 1        NA     0.6037310 new_group
#> 2        NA     0.4638874 new_group
#> 3        NA     0.4062116 new_group

Cause

Adaptation changes the predictive location or scale; it does not rewrite the historical fact that site E was outside the original site’s factor levels. allow_extrapolation = TRUE exposes the adapted distribution, while the support column preserves that boundary.

Fix

Use this override only after fitting a deliberate transport step from an appropriate local reference sample. Keep support in every exported score table, report how many local controls estimated the transport, and validate on different local observations. See Brain charts across sites for the full comparison among population scoring, adaptation, calibration, and refitting.

Why does a small local correction look weak or fail to calibrate?

Symptom

An adaptation offset is smaller than the raw local mean difference, or a calibration object built from a tiny sample leaves calibrated == FALSE.

site_e_offset <- adapted$adaptation$offsets$y$E
data.frame(
  n = site_e_offset$n,
  location = site_e_offset$location,
  location_se = site_e_offset$location_se
)
#>    n  location location_se
#> 1 20 0.9672202   0.2613119

A single probability value cannot define the package’s monotone calibration map, so no row is reported as calibrated.

tiny_calibration <- ref_calibrate(
  site_fit,
  data = local_controls[1, , drop = FALSE],
  by = site,
  uncertainty = "conditional"
)
tiny_scores <- predict(
  tiny_calibration,
  site_e_target,
  uncertainty = "conditional",
  allow_extrapolation = TRUE
)

data.frame(
  calibration_n = tiny_calibration$calibration$n,
  any_mapped = any(tiny_scores$calibrated)
)
#>   calibration_n any_mapped
#> 1             1      FALSE

Cause

ref_adapt() deliberately shrinks local offsets toward the shared reference; the amount depends on the local sample size and the prior strength. ref_calibrate() instead estimates an empirical probability map. With fewer than two finite calibration values there is no map, and with only a modest number its tail resolution is necessarily coarse.

Fix

Choose the correction that matches the failure. Adapt when a scanner or domain plausibly shifts location or scale while preserving the shared trajectory. Calibrate when held-out diagnostics show probability miscalibration that a location/scale shift cannot repair. Reserve separate local rows for evaluation, report the local sample size, and inspect centre, spread, coverage, and tails; do not judge a correction on the rows that estimated it.

Why did one outcome fail while the rest of the panel fitted?

Symptom

ref_fit() returns a panel object, but tidy() reports unsupported_type for one outcome.

panel <- reference
panel$group_label <- factor(ifelse(panel$y > median(panel$y), "high", "low"))

panel_fit <- ref_fit(
  spec,
  data = panel,
  outcomes = c("y", "group_label")
)
tidy(panel_fit)[, c("outcome", "status", "message")]
#> # A tibble: 2 × 3
#>   outcome     status           message
#>   <chr>       <chr>            <chr>
#> 1 y           ok               <NA>
#> 2 group_label unsupported_type Outcome is factor; only numeric outcomes are sup…

Cause

The current model families define continuous predictive distributions for numeric outcomes. A factor, character, or logical outcome is recorded as a failed panel member instead of aborting every other outcome.

Fix

Select the intended numeric measurements explicitly. Convert an encoded numeric column only when its scale has a defensible quantitative meaning; do not turn categories into arbitrary integer codes to make the status disappear. Inspect tidy(fit) before deployment and decide whether any partial panel is acceptable for the study.

Why is an assessment too optimistic?

Symptom

The model looks well calibrated when assessed on the same rows used to fit it. ref_assess() now rejects this reuse by default. An explicitly requested training diagnostic remains labelled in_sample; it is not validation.

training_assessment <- ref_assess(fit, newdata = reference, allow_in_sample = TRUE)
held_out_assessment <- ref_assess(fit, newdata = held_out)

data.frame(
  sample = c("training", "held out"),
  in_sample = c(training_assessment$in_sample, held_out_assessment$in_sample),
  var_z = c(training_assessment$marginal$var_z, held_out_assessment$marginal$var_z),
  cover_95 = c(training_assessment$marginal$cover_95,
               held_out_assessment$marginal$cover_95)
)
#>     sample in_sample     var_z  cover_95
#> 1 training      TRUE 0.9706020 0.9625000
#> 2 held out     FALSE 0.7216211 0.9895833

Cause

The fitted rows influenced the curve against which they are being scored. Smoothness selection and parameter estimation can make their deviations look smaller than deviations for genuinely new observations.

Fix

Pass a held-out reference sample to ref_assess(). When a separate sample is not available, use out-of-fold scores from ref_crossfit() and keep whole subjects together with cluster = when rows repeat within person.

oof <- ref_crossfit(spec, reference, outcomes = "y", folds = 3)
table(oof$.in_sample)
#>
#> FALSE
#>   240

Why are longitudinal change scores all NA?

Symptom

The marginal model fits, but ref_dynamics() reports an unidentified process and ref_transition() returns no innovation or change Z-scores.

sparse <- ref_simulate(80, seed = 20)
sparse$participant_id <- seq_len(nrow(sparse))
sparse$participant_id[1:3] <- sparse$participant_id[4:6]

sparse_fit <- ref_fit(
  spec,
  data = sparse,
  outcomes = "y",
  id = participant_id
)
sparse_dynamics <- ref_dynamics(
  sparse_fit,
  data = sparse,
  id = participant_id,
  time = age,
  crossfit = 0
)
sparse_dynamics$components[, c("identified", "stable", "dynamic", "measurement")]
#> # A tibble: 1 × 4
#>   identified stable dynamic measurement
#>   <lgl>       <dbl>   <dbl>       <dbl>
#> 1 FALSE          NA      NA          NA

sparse_transitions <- ref_transition(
  sparse_dynamics,
  data = sparse,
  id = participant_id,
  time = age
)
sparse_transitions[, c(".id", "innovation_z", "change_z", "support")]
#> # A tibble: 3 × 4
#>     .id innovation_z change_z support
#>   <int>        <dbl>    <dbl> <chr>
#> 1     4           NA       NA unidentified
#> 2     5           NA       NA unidentified
#> 3     6           NA       NA unidentified

Cause

Only three subjects have repeat observations here. That is not enough to identify within-person dependence, so the package fails closed instead of turning starting values into longitudinal evidence. Even when change is identifiable, separating stable rank, dynamic variation, and measurement noise requires informative visit timing, such as three or more visits for some subjects or short-interval repeats.

Fix

Treat visit design as part of the estimand. Collect repeated observations from enough subjects, include varied lags when estimating a time-decaying process, and inspect components, identified, ell_identified, and support before interpreting change. Longitudinal change explains the stable fixed-lag fallback, innovation Z, change Z, velocity, and forecasting.

Where should you go next?

Return to Getting started with referent for the ordinary fit-score-assess path. Use Brain charts across sites for site transport, joint scores, and reporting; use Longitudinal change for repeated visits; and use Coming from PCNtoolkit when translating an existing PCNtoolkit workflow. The function reference for predict.ref_fit(), ref_support(), ref_assess(), and ref_dynamics() gives the exact status and return-value contracts used here.