Coming from PCNtoolkit

PCNtoolkit users will recognise most of referent. The estimand is the same: a conditional predictive distribution of an outcome given covariates, estimated on a reference population and used to place new people. The differences are in what is exposed. referent returns the distribution itself as a first-class object and derives every score from its CDF, it separates site adaptation (shifting the distribution) from calibration (re-mapping its probabilities), and it does not define an abnormality index.

library(referent)
ref <- ref_simulate(300, seed = 1)
val <- ref_simulate(150, seed = 2)
fit <- ref_fit(ref_spec(ref_gaussian(), ~ s(age, k = 6) + sex), ref, "y")
scores <- predict(fit, newdata = val)
assessment <- ref_assess(fit, newdata = val)

Concepts

PCNtoolkit	referent	Notes
normative model (BLR, HBR)	ref_spec() + ref_fit()	one mgcv model per outcome; Gaussian or SHASH families with formulas per parameter; GPR belongs to the legacy PCNtoolkit interface, not 1.x
response variable	outcomes =	several outcomes fit as a panel, kept wide
covariates, batch effects	formula terms	s(site, bs = "re") plays the role of a batch random effect
predict() / estimate()	predict(fit, newdata)	returns scores or <distribution> columns
cross-validation (cvfolds)	ref_crossfit()	out-of-fold scores for every reference row; the deployment refit is attached
warp (SHASHb, SinArcsinh)	ref_shash()	sinh-arcsinh likelihood via mgcv::shash()
NormData batch_effects	s(site, bs = "re"); ref_calibrate(by = )	the batch variable is a random intercept in the formula; group-specific probability re-mapping is a separate step
HBR random slopes	not available; s(age, site, bs = "fs")	bs = "re" is a random intercept only; mgcv’s factor-smooth basis gives each site its own age curve shrunk toward the shared one (see below)
HBR Beta likelihood	not available	families are Gaussian and SHASH on a monotone transform of the outcome; a bounded outcome can be logit-transformed by hand before fitting
GPR	not available	smooth terms (s(age), te(age, x)) take the place of a GP mean; their uncertainty enters uncertainty = "total" through the coefficient covariance
warped BLR (WarpSinArcsinh, WarpBoxCox)	ref_spec(transform = ) and ref_shash()	a Box-Cox or log transform of the response is transform =; a sinh-arcsinh likelihood is ref_shash(); the two combine
fit_predict	ref_crossfit() + ref_fit()	cross-fit for out-of-fold reference scores, then fit on all rows for deployment
site transfer / adapt / transfer	ref_adapt()	shrunk location (and scale) offsets on the family’s own density, trajectory frozen
(no direct analogue)	ref_calibrate()	monotone re-mapping of PIT values; optional by = site
extend	refit on pooled data	ref_fit(spec, rbind(reference, new_site)) is the extend operation; the reference changes for every site
merge / federated models	none	the unit of sharing is a frozen bundle (ref_freeze()), not model parameters to be merged
harmonised outputs (Y_harmonized)	predict(type = "harmonised")	conditional quantile mapping out of the site effect (or to a chosen to level); the deviation score itself is already site-free
extreme value / abnormality	not provided	see below

Harmonised values

When a downstream analysis needs raw values on one scale, the conditional quantile map gives them: a value at the 80th centile of its own site’s predictive becomes the 80th centile of the target predictive.

sites <- ref_simulate(300, sites = 3, site_shift = c(0, 1, -1), seed = 5)
site_fit <- ref_fit(
  ref_spec(ref_gaussian(), ~ s(age, k = 5) + sex + s(site, bs = "re")),
  sites, "y"
)
harmonised <- predict(site_fit, sites, type = "harmonised")
rbind(raw = tapply(sites$y, sites$site, mean),
      harmonised = tapply(harmonised$y, harmonised$site, mean))
#>                   A        B        C
#> raw        9.570536 11.15225 9.171057
#> harmonised 9.877797 10.02819 9.988281

Random slopes by site

s(site, bs = "re") shifts the level per site. If sites also differ in the shape of the age trend, mgcv’s factor-smooth interaction, s(age, site, bs = "fs"), fits one smooth per site with a shared penalty that shrinks each curve toward the common one. It is accepted anywhere a formula is, and an unseen site still receives the population curve with support == "new_group".

multi <- ref_simulate(600, sites = 5, seed = 4)
fs_spec <- ref_spec(
  ref_gaussian(),
  location = ~ s(age, k = 6) + sex + s(age, site, bs = "fs", k = 5)
)
fs_fit <- ref_fit(fs_spec, multi, "y")
tidy(fs_fit)[, c("status", "edf_s(age)", "edf_s(age,site)")]
#> # A tibble: 1 × 3
#>   status `edf_s(age)` `edf_s(age,site)`
#>   <chr>         <dbl>             <dbl>
#> 1 ok             2.87            0.0956

This is not HBR: the site curves are penalised deviations from a shared smooth, not draws from a hierarchical prior with its own estimated hyperparameters, and ref_adapt() transports only the level and scale, not the slope. Use it when the reference contains enough rows per site to support a curve each, and cross-fit with strata = site to check that it earns its effective degrees of freedom.

Metrics

ref_assess() reports related evaluation quantities alongside proper scoring rules, but shared names do not guarantee the same estimand. The table labels each relationship rather than silently treating near-matches as identical. The PCNtoolkit side follows its 1.x evaluation-metrics definitions, not the legacy 0.35 interface.

PCNtoolkit	referent column	Equivalence	Relationship
Z	z	Conceptual analogue	Both place an observation in a predictive distribution. referent always computes qnorm(centile) from the distribution CDF, including non-Gaussian families.
yhat	median	Exact when yhat is the predictive median	For a Gaussian the median also equals the mean. Check the model family’s definition before comparing other families.
legacy BLR s2	variance(dist)	Model-specific analogue	Both can include observation and coefficient uncertainty, but PCNtoolkit 1.x does not expose s2 as a general output and HBR posterior averaging need not match referent’s mixture distribution.
MSLL	no exact column	Not numerically equivalent	PCNtoolkit subtracts a Gaussian baseline fitted to the sample being scored and uses lower-is-better loss. standardized_log_score uses the frozen reference-sample baseline and higher-is-better log score. Rescore both against one declared baseline before comparing.
SMSE	smse	Near-equivalent	Both divide MSE by outcome variance, but PCNtoolkit uses a population variance denominator while R’s var() uses a sample denominator. On the same complete rows, referent’s value is PCNtoolkit’s value times (n - 1) / n.
EV (EXPV)	ev	Algebraically equivalent on nonconstant complete data	Both compute 1 - var(residual) / var(y). Constant outcomes and filtering of non-finite rows can produce different edge behaviour.
Rho	no exact column	Different statistic	PCNtoolkit 1.x reports Spearman rank correlation. cor is Pearson correlation of the observed value and predictive median.
RMSE	rmse	Exact on the same included rows	Missing and non-finite row handling must also match.
MACE	mace	Exact only without batch effects	Both use the mean absolute calibration gap on a centile grid. With batch combinations, PCNtoolkit 1.3.0 averages the batch-level values equally; referent currently pools rows.
ShapiroW	shapiro_w	Same statistic	Shapiro-Wilk W of Z; referent subsamples 5000 scores with a fixed seed when n is larger, because shapiro.test() accepts no more.
Skew, Kurtosis	skew_z, excess_kurtosis_z	Related, not exact	PCNtoolkit 1.3.0 uses SciPy’s bias-corrected skewness and excess kurtosis. referent reports unadjusted standardised moments (and also uses excess kurtosis). The difference is material for small samples.
(none)	mean_log_score, crps	referent additions	Proper scores of the full predictive distribution; larger log score and smaller CRPS are better.

One implementation detail matters for external rescoring: with outcome standardisation enabled, PCNtoolkit 1.3.0 maps predictions and centiles back to the outcome scale but returns logp on the standardised scale. Subtract log(sd(y_train)) (using the population standard deviation used by the scaler), plus any response-transform Jacobian, before comparing response-scale log scores.

assessment$overall[, c("standardized_log_score", "smse", "ev", "cor", "rmse", "crps")]
#> # A tibble: 1 × 6
#>   standardized_log_score  smse    ev   cor  rmse  crps
#>                    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1                  0.455 0.416 0.609 0.794  1.32 0.743

If you need PCNtoolkit’s Rho, compute it explicitly from the held-out score rows. Pearson and Spearman answer different questions and can disagree even on the same predictions.

metric_rows <- assessment$scores
c(
  pearson = cor(metric_rows$observed, metric_rows$median, method = "pearson"),
  spearman_rho = cor(metric_rows$observed, metric_rows$median, method = "spearman")
)
#>      pearson spearman_rho
#>    0.7927178    0.7831264

The normality metrics are in the marginal table. The chunk after it recomputes them from the same held-out scores (assessment$scores; the scores object above was computed under uncertainty = "total" and differs slightly), which is also how to get them on a different centile grid or without the kurtosis offset.

assessment$marginal[, c("mean_z", "var_z", "skew_z", "excess_kurtosis_z",
                        "shapiro_w", "mace")]
#> # A tibble: 1 × 6
#>   mean_z var_z skew_z excess_kurtosis_z shapiro_w   mace
#>    <dbl> <dbl>  <dbl>             <dbl>     <dbl>  <dbl>
#> 1  0.240 0.811 -0.166             0.251     0.989 0.0605
held <- assessment$scores
z <- held$z[is.finite(held$z)]
u <- held$centile[is.finite(held$centile)]
grid <- c(0.05, 0.25, 0.5, 0.75, 0.95)
c(
  MACE = mean(vapply(grid, function(q) abs(q - mean(u <= q)), numeric(1))),
  ShapiroW = unname(shapiro.test(z)$statistic),
  Skew = mean((z - mean(z))^3) / sd(z)^3,
  Kurtosis = mean((z - mean(z))^4) / sd(z)^4
)
#>        MACE    ShapiroW        Skew    Kurtosis
#>  0.06054054  0.98930171 -0.16584463  3.25101784

d <- predict(fit, newdata = val[1:3, ], type = "distribution")$y
data.frame(
  yhat = median(d),
  s2 = variance(d),
  Z = scores$z[1:3]
)
#>        yhat       s2          Z
#> 1  8.017298 2.014648 0.50918165
#> 2 10.645450 2.015517 0.04345949
#> 3 10.667744 2.016018 0.86435217

Z here includes uncertainty in the fitted curves (uncertainty = "total", the default). That is closest to legacy BLR’s coefficient-aware predictive variance, but it is not a general claim of numerical equivalence across PCNtoolkit model families. In particular, PCNtoolkit HBR averages draw-specific Z scores and quantiles, while referent derives them from a predictive mixture. Those are different estimands even with identical draws. uncertainty = "conditional" gives the plug-in score.

What numerical agreement means

The executable comparison target is PCNtoolkit 1.3.0. Fixed Gaussian distributions and transformed-distribution algebra are checked as exact; SHASH distributions are checked after an explicit parameter conversion; deliberately matched fitted models use prediction-level tolerances. When the estimators differ, agreement is assessed by held-out proper scores and calibration rather than by forcing their parameters or predictions to match.

These are independent claims: better held-out performance cannot excuse a failure of shared distribution algebra, and numerical agreement does not by itself establish superiority. The complete classifications, tolerances, and CI tiers are recorded in the repository’s docs/design/pcntoolkit-validation-contract.md, with current machine-readable results under docs/evidence/pcntoolkit/v1.3.0/. Both paths are linked from the repository README.

Site adaptation

PCNtoolkit’s adaptation re-estimates site-specific parameters of a pre-trained model from a small local sample. ref_adapt() does the same with the shared trajectory frozen, estimating a location offset (and optionally a log-scale offset) per site by penalised maximum likelihood, with shrinkage set in pseudo-observations.

local <- ref_simulate(40, site_shift = 1, seed = 3)
adapted <- ref_adapt(fit, data = local, parameters = c("location", "scale"))
adapted$adaptation
#> <ref_adaptation> parameters: location, scale
#> local n = 40

PCNtoolkit has no separate calibration step. referent keeps the two ideas apart because they fix different problems: adaptation moves the distribution when a new site’s level or spread differs; calibration re-maps probabilities when the distributional family is wrong in a way that is common to everyone. vignette("brain-charts-across-sites") compares both against a refit on held-out data.

What is deliberately missing

PCNtoolkit offers extreme-value and abnormality statistics that summarise how far a profile lies in the tails. referent does not provide an abnormality index. A large |z| is a threshold exceedance whose expected count under calibration is known; ref_flag() reports that count and an FDR-adjusted tail probability, and ref_joint() gives a calibrated joint centile across outcomes. Whether an exceedance means anything is a question for the study design, not the scoring function.

To state the rest plainly, referent does not do the following.

• Hierarchical Bayesian regression. There is no MCMC, no posterior over site slopes, and no Beta or other bounded likelihood. Site effects are penalised random intercepts (or factor smooths); uncertainty in the curves is Gaussian in the coefficients.
• Gaussian-process regression. Smooth terms play that role.
• Harmonisation as a separate model. predict(type = "harmonised") maps each value through its own site’s predictive CDF and the target site’s quantile function, so it removes exactly the site effect the reference model estimated and nothing else. There is no ComBat-style empirical-Bayes step, and the deviation score, which is site-free by construction, is the quantity meant for downstream use.
• Model merging or federated fitting. A reference is fitted on data in one place and shipped as a bundle. Extending a reference means refitting on the pooled data.
• Classification or clinical thresholds. No function maps a score to a label.

Continue with Troubleshooting reference-model workflows (vignette("troubleshooting", package = "referent")), or return to Getting started with referent for the default end-to-end workflow.