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 DRModels, 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.9978887393881511, 0.5008810755480976], [-0.1740319656027385, 0.5525691658062271])

Predicting the mean at new data

predict evaluates the fitted mean at a newdata column table. type = :response (the default) returns invlink(Xβ̂) with every random effect set to 0; pass type = :link for the linear predictor Xβ̂. Those random or structured effects are not integrated out, so on a random-effect fit this is the prediction for an average group, not the mean averaged over groups.

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

3-element Vector{Float64}:
 0.49700766384005357
 0.9978887393881511
 1.4987698149362487

Standard errors (the delta method)

Add se = true for delta-method standard errors — the glmmTMB se.fit convention (drmTMB's predict() has no standard-error option). They carry the uncertainty in β̂ only. 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.02823637531426892
 0.04209049180605554
 0.06614441643336108

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.02823637531426892
 0.04209049180605554
 0.06614441643336108

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.48354972215718167, 0.8402700374909216, 1.4601471235578256], [0.022384472882966716, 0.02662502613092813, 0.06344325301660712])

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.27826808450350565
 0.48354972215718167
 0.8402700374909216
 1.4601471235578256
 2.5373148241730887

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) for the mean and sigma of the GLM-type families (Gaussian, Poisson, NegBinomial2, Gamma, Beta, Binomial), as here (the predict_parameters docstring lists the exceptions):

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

([0.813795486779685, 0.6540498515422992, 0.8175584898849235], [0.6858319862415315, 0.5750184916309026, 0.6886850146547489])

Residuals: response and quantile

residuals defaults to response residuals (y − μ̂). For a distribution-aware diagnostic — the DHARMa/glmmTMB randomized quantile residual — pass type = :quantile. For a correctly specified model without random effects, like this one, they are approximately 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.004829463709711091, sd = 1.0009898288985777)

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 without random effects yields ≈ N(0, 1), on the same scale for every family. They are implemented for every family except Tweedie and SkewNormal; for zero-inflated and hurdle count fits the CDF is that of the whole mixture, zero part included. Random-effect and temporal() fits judge each observation against a different reference (a random intercept integrated out, or the fitted temporal path), so what "≈ N(0, 1)" means changes; residuals lists what each kind of fit returns.

Accuracy limit when the random-effect variance is large (#1103)

Where a random intercept is integrated out by the 32-node Gauss–Hermite mixture (non-Gaussian (1 | g) on the mean, a Gaussian sigma ~ 1 + (1 | g) intercept, and structured effects that keep the mixture), the conditional CDF becomes a near-step in the node variable once σ̂_b ≫ σ̂, and 32 nodes stop resolving it. Measured for a Gaussian mean intercept (max |Φ⁻¹(mixture) − exact z| over |z| ≤ 4): about 1e-6 at σ_b/σ = 2, 6e-3 at 4, 0.02 at 5, 0.12 at 10 and 0.21 at 20. Read the residuals from those routes as approximate when σ̂_b/σ̂ is above about 3. Gaussian (1 | g) on the mean and Gaussian / LogNormal phylo(1 | g) already use the closed-form marginal and are not affected. Separately, every continuous family computes Φ⁻¹(clamp(Φ(z), eps, 1 − eps)), which loses upper-tail precision (error about 6e-6 at z = 7 and 8e-3 at z = 8) and caps |r| at about 8.13.

Quantile residuals with a structured mean: DRModels.jl vs drmTMB

For a fit whose mean has a phylo(), relmat() or animal() random intercept, the two packages do not return the same numbers, although the fits agree (same log-likelihood to about 6e-8 on the checked Gaussian phylo(1 | species) model).

Convention	Residual	Who returns it
marginal	(y − Xβ̂) / √(σ̂² + σ̂_b² C_ii)	DRModels.jl (default Gaussian phylo(1 | g) route)
conditional	(y − Xβ̂ − û) / σ̂	drmTMB's residuals(type = "quantile") code
fixed-effect only	(y − Xβ̂) / σ̂	what drmTMB's residual help text describes

On that model the two engines agree to within 3.5e-9 under either convention, so the gap is the convention, not the fit. Marginal residuals stay near N(0, 1) for each row, but rows that share ancestry are correlated, so read the whole set (a QQ plot, the SD) and do not hold each row to N(0, 1). Do not compare the two packages' quantile residuals number by number for these fits. residuals lists the convention for every other kind of fit. Which convention drmTMB should use is an open drmTMB-side decision (#1101).

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 supports the moving-scale term within these nested, checked candidate models. Also inspect the estimated scale change, residual diagnostics, and whether that effect answers the biological question. 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.