Prediction, residuals & model comparison

Status — Stable

The post-fit surface a fitted drm model exposes: prediction (point + standard error, per parameter, on a grid), quantile residuals, and model-comparison helpers. Mirrors the predict / residuals / anova verbs you reach for after a drmTMB fit.

A fit is the beginning, not the end. Once drm returns, you usually want to predict at new covariate values, check the model with residuals, and compare it against a simpler one. This page walks the whole surface on one small location–scale model.

using DRM, Random, Statistics
Random.seed!(11)

n = 500
x = randn(n)
# mean μ = 1 + 0.5x ; log-scale logσ = -0.2 + 0.6x  (the spread grows with x)
y = 1.0 .+ 0.5 .* x .+ exp.(-0.2 .+ 0.6 .* x) .* randn(n)
dat = (; y, x)

fit = drm(bf(@formula(y ~ 1 + x), @formula(sigma ~ 1 + x)), Gaussian(); data = dat)
coef(fit, :mu), coef(fit, :sigma)

([0.9978887393874984, 0.5008810755477906], [-0.17403196560290976, 0.5525691658061693])

Predicting the mean at new data

predict evaluates the fitted mean at a newdata column table. By default it returns the response scale; pass type = :link for the linear predictor.

nd = (; x = [-1.0, 0.0, 1.0])
predict(fit, nd)                 # μ̂ at x = -1, 0, 1  (response scale)

3-element Vector{Float64}:
 0.49700766383970785
 0.9978887393874984
 1.498769814935289

Standard errors (the delta method)

Add se = true for delta-method standard errors — the glmmTMB/drmTMB se.fit. You get back a NamedTuple (; prediction, se):

ps = predict(fit, nd; type = :link, se = true)
ps.se                            # SE of the linear predictor at each new row

3-element Vector{Float64}:
 0.02823637531426556
 0.042090491806048215
 0.06614441643335006

On the response scale the SE is chained through the inverse-link derivative automatically (for the Gaussian identity link the two coincide):

predict(fit, nd; type = :response, se = true).se

3-element Vector{Float64}:
 0.02823637531426556
 0.042090491806048215
 0.06614441643335006

Predicting every distributional parameter

A distributional model has more than a mean. predict_parameters returns each parameter at the new data — here μ and σ:

predict_parameters(fit, nd)      # Dict(:mu => …, :sigma => …), response scale

Dict{Symbol, Vector{Float64}} with 2 entries:
  :mu    => [0.497008, 0.997889, 1.49877]
  :sigma => [0.48355, 0.84027, 1.46015]

With se = true each entry becomes a (; value, se) pair:

pp = predict_parameters(fit, nd; se = true)
pp[:sigma].value, pp[:sigma].se  # fitted σ and its SE at x = -1, 0, 1

([0.48354972215712677, 0.8402700374907776, 1.4601471235574912], [0.022384472882962025, 0.026625026130919105, 0.06344325301657909])

Sweeping a predictor: prediction_grid

To trace a parameter across a covariate, build a grid with prediction_grid — it sweeps the named predictor(s) over a range and holds everything else at a reference value, then composes straight into predict_parameters:

grid = prediction_grid((; x = dat.x), x = range(-2, 2; length = 5))
predict_parameters(fit, grid)[:sigma]   # σ̂ rising across the x sweep

5-element Vector{Float64}:
 0.27826808450349017
 0.48354972215712677
 0.8402700374907776
 1.4601471235574912
 2.5373148241723613

In-sample fitted parameters

For the per-observation fitted parameters at the training data there is a cheap accessor that reads straight from the fit (no recomputation) — marginal_parameters. In-sample it equals predict_parameters(fit, data):

mp = marginal_parameters(fit)
first(mp[:mu], 3), first(mp[:sigma], 3)

([0.8137954867791451, 0.6540498515418572, 0.8175584898843813], [0.6858319862414286, 0.5750184916308269, 0.6886850146546453])

Residuals: response and quantile

residuals defaults to response residuals (y − μ̂). For a distribution-aware diagnostic — the DHARMa/glmmTMB randomized quantile residual — pass type = :quantile. Under a correctly specified model these are standard normal, which makes a QQ-plot or a simple moment check interpretable even for non-Gaussian families:

rq = residuals(fit; type = :quantile)
(mean = mean(rq), sd = std(rq))         # ≈ (0, 1) when the model fits

(mean = -0.004829463709049791, sd = 1.0009898288987558)

Why quantile residuals

Raw y − μ̂ residuals are misleading when the variance changes with the mean (exactly the location–scale case here) or for discrete responses. Quantile residuals fold each observation through its own fitted CDF, so a well-specified model always yields ≈ N(0, 1) — the scale is the same across families. (Implemented for Gaussian and Poisson today; other families are tracked in issue #183.)

Does the extra structure earn its keep?

Fit a simpler model with a constant scale and compare. lrtest (alias anova) runs the nested likelihood-ratio test; it returns a NamedTuple (; statistic, dof, pvalue):

reduced = drm(bf(@formula(y ~ 1 + x), @formula(sigma ~ 1)), Gaussian(); data = dat)
lrtest(reduced, fit)             # is σ ~ x worth the one extra parameter?

(statistic = 301.74324091067547, dof = 1, pvalue = 1.373983432463021e-67)

A small p-value says the moving scale is real signal. The information criteria agree — lower is better, and aicc is the small-sample-corrected AIC:

(aic = aic(fit),   bic = bic(fit),   aicc = aicc(fit)),
(aic = aic(reduced), bic = bic(reduced), aicc = aicc(reduced))

((aic = 1287.7430095377558, bic = 1304.6014419314447, aicc = 1287.823817618564), (aic = 1587.4862504484313, bic = 1600.1300747436978, aicc = 1587.5346375452054))

Refitting: update

update re-fits with a new formula, reusing the original family — handy for building a comparison ladder. (A DrmFit does not retain its data, so pass data again.)

mean_only = update(fit, bf(@formula(y ~ 1 + x), @formula(sigma ~ 1)); data = dat)
length(coef(mean_only)) == length(coef(reduced))

true

See also

• When variance carries signal — the model fit here, in depth.

• Rosetta (R ↔ Julia) — the verb-by-verb drmTMB mapping.

• Model fitting & post-fit reference — every accessor's docstring.