DRModels.jl

Appearance

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.

julia
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. By default it returns the response scale; pass type = :link for the linear predictor.

julia
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/drmTMB se.fit. You get back a NamedTuple (; prediction, se):

julia
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):

julia
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 σ:

julia
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:

julia
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:

julia
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):

julia
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:

julia
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.

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):

julia
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:

julia
(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.)

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

See also ​