Distributional outputs and adequacy

A fitted mu formula answers “what is the expected response, given the covariates?” It does not answer two other questions an applied paper usually needs: “did I model the whole conditional distribution, not just its average?” and “how extreme, or how likely to exceed a threshold, is this observation for its covariates?” A model can have a perfectly reasonable mean structure and still misrepresent the spread, tails, or boundary behaviour of the response – and coefficients alone never reveal that. This article follows the check_drm() step in the post-fit workflow described in Checking and using fitted models.

This article covers drmTMB’s distributional-output and adequacy layer, built on one shared foundation: every fitted family exposes its own density, CDF, and quantile function at the fitted parameters (fitted_distribution()). Two families of outputs are built on it:

• Adequacy diagnostics – residuals(fit, type = "quantile"), worm_plot(), qq_plot() – randomized quantile residuals (Dunn & Smyth
  1996. that should look like a standard normal sample under a correctly specified fixed-effect distributional form.
• Distributional outputs – predict(fit, type = "quantile"), exceedance(), centile_chart() – conditional quantiles, exceedance probabilities, and centile curves at the fitted model.

The GAMLSS Primer (Merder et al., Nature Reviews Methods Primers, 2026, 6:49) provides background on distributional regression and model checking. The contrast below is a package-generated Gaussian example with heteroscedastic data: a location-only fit has constant sigma, while the matching location-scale fit lets sigma vary with x.

Setup

set.seed(20260712)
n <- 300
x <- stats::rnorm(n)
sigma_true <- exp(0.1 + 0.9 * x)
y <- 0.3 + 0.6 * x + stats::rnorm(n, sd = sigma_true)
dat <- data.frame(y = y, x = x)

y is generated with a residual SD that grows with x: a Gaussian location-scale problem, not a location-only one.

Quantile residuals and worm/QQ plots on a correctly specified fit

Fit the model that matches the data-generating process – mu and sigma both depend on x – and compute randomized quantile residuals.

fit_true <- drmTMB(bf(y ~ x, sigma ~ x), family = gaussian(), data = dat)
r_true <- residuals(fit_true, type = "quantile")
stats::ks.test(r_true, "pnorm")
#> 
#>  Asymptotic one-sample Kolmogorov-Smirnov test
#> 
#> data:  r_true
#> D = 0.031089, p-value = 0.9339
#> alternative hypothesis: two-sided

The residuals are indistinguishable from a standard normal sample: a Kolmogorov-Smirnov test does not reject. worm_plot() (a detrended QQ plot) and qq_plot() show the same thing visually – a flat scatter around the reference line is “no detectable departure from N(0,1),” never proof the model is correct.

worm_plot(fit_true) +
  ggplot2::labs(
    subtitle = "Fixed-effect adequacy: no detectable departure from N(0,1)"
  )

Worm plot of quantile residuals for the correctly specified location-scale fit; the smooth stays near the zero reference across the bulk, with only sparse-tail drift (the KS test does not reject).

qq_plot(fit_true)

Normal QQ plot of quantile residuals for the correctly specified fit; points on the diagonal are no detectable departure from N(0,1).

A location-only fit to heteroscedastic data

Now fit a location-only model to the same data – sigma held constant, so the growing spread in y has nowhere to go but into the residuals.

fit_mis <- drmTMB(bf(y ~ x, sigma ~ 1), family = gaussian(), data = dat)
r_mis <- residuals(fit_mis, type = "quantile")
stats::ks.test(r_mis, "pnorm")
#> 
#>  Asymptotic one-sample Kolmogorov-Smirnov test
#> 
#> data:  r_mis
#> D = 0.19845, p-value = 1.093e-10
#> alternative hypothesis: two-sided

The KS test now rejects sharply. The worm plot shows a systematic bend: at extreme covariate values, observations are consistently too close to or too far from mu because the constant sigma cannot track the growing variance.

worm_plot(fit_mis) +
  ggplot2::labs(
    subtitle = "Fixed-effect adequacy: bend flags a detectable departure from N(0,1)"
  )

Worm plot of quantile residuals for the location-only fit to heteroscedastic data; the systematic bend flags a mis-specified distributional shape.

Both worm_plot() and qq_plot() apply this package’s Confidence-Eye figure contract: hollow points, no filled dots, no interval bars around the reference line – the reference line is a null-hypothesis benchmark, not an uncertainty interval.

The two fits have the same mu formula and produce similar coefficient tables; only the adequacy check separates them. That is the point of a distributional diagnostic: it catches a mis-specification a mean-only residual plot would miss.

Distributional outputs: quantiles, exceedance, and centiles

Beyond checking the fit, the same fitted_distribution() foundation answers conditional-quantile and threshold questions directly, using the correctly specified fit_true.

Conditional quantiles

newdata <- data.frame(x = c(-1, 0, 1))
predict(fit_true, newdata = newdata, type = "quantile", prob = c(0.025, 0.5, 0.975))
#>           2.5%        50%    97.5%
#> [1,] -1.192822 -0.2513602 0.690102
#> [2,] -1.802150  0.3907378 2.583625
#> [3,] -4.074916  1.0328358 6.140588
#> attr(,"calibrated")
#> [1] FALSE
#> attr(,"prob")
#> [1] 0.025 0.500 0.975
#> attr(,"label")
#> [1] "distributional (plug-in) interval"

The output widens at extreme x because sigma grows with x in this model – exactly the heteroscedasticity the location-only fit above could not represent.

Exceedance probability

exceedance() returns Pr(Y > threshold | x) at the fitted model, a thin wrapper over the shared CDF.

exceedance(fit_true, threshold = 2, newdata = newdata)
#> [1] 1.386542e-06 7.517022e-02 3.552728e-01
#> attr(,"calibrated")
#> [1] FALSE
#> attr(,"lower.tail")
#> [1] FALSE
#> attr(,"threshold")
#> [1] 2 2 2

Centile chart

centile_chart() draws fitted response centiles against one covariate, holding every other predictor at a reference value.

centile_chart(fit_true, covariate = "x") +
  ggplot2::scale_colour_viridis_d()

Model-conditional centile chart for the correctly specified location-scale fit; centiles fan out with x because sigma grows with x.

The centiles fan out as x increases – a direct visual reading of the fitted sigma ~ x relationship that a mean-only prediction plot would not show.

What this layer detects – and what it does not

Detects. The adequacy diagnostic reliably catches distributional shape or atom mis-specification that a family cannot reabsorb through its own free parameters: heavy tails fit as Gaussian, ignored overdispersion or zero-inflation in a family with no free dispersion parameter, ignored truncation, a missing zero/one atom, and the heteroscedasticity example reproduced above. In these situations, a clear departure in the residual plots is a reason to reconsider the family or its distributional formulas; it is not, by itself, proof of the replacement model.

Does not detect – a structural blind spot, not evidence of adequacy. When a fitted family’s own free nuisance, dispersion, or inflation parameter can absorb a mis-specification, the marginal quantile residual stays close to N(0,1) even though the mean or variance structure is wrong. Simulation examples show heteroscedasticity absorbed by Student-t nu (power 0.035, versus 1.0 for the same heteroscedasticity under Gaussian, which has no absorbing parameter); missing zero-inflation absorbed by nbinom2’s sigma (power 0.035, versus 0.96 under Poisson); and a constant-vs-covariate zero-inflation/hurdle/zero-one-inflation mechanism, which splits in two: for hurdle_nbinom2/zero_one_beta power stays flat near the type-I rate regardless of sample size (a genuine structural blind spot); for zi_poisson/zi_nbinom2 power rises with sample size but stays far below 0.8 even at n = 3000, so it is impractical to detect at realistic sample sizes rather than permanently undetectable. A mean-structure diagnostic, not this one, is what catches an absorbed mis-specification.

Wording and scope, throughout this package:

• Passing this diagnostic means “no detectable departure from N(0,1),” and is worded that way everywhere – never “the model is adequate,” “valid,” or “correct.”
• residuals(type = "quantile"), worm_plot(), and qq_plot() are fixed-effect adequacy: for a random-effect or structured fit, they are conditional on the fixed-effect prediction, not a marginal check of the full model.
• predict(type = "quantile") and exceedance() results carry attr(., "calibrated") <- FALSE: they are distributional (plug-in) outputs at the fitted parameter point estimate, and do not propagate parameter uncertainty. centile_chart() is a model-conditional curve at that same point estimate, not a WHO-style or other population reference standard.
• Quantile predictions and residual diagnostics do not establish that a confidence interval is appropriate. For family-specific limits and inference guidance, see Capability and limits.

See also Checking and using fitted models for the post-fit workflow that starts with these diagnostics.