DRModels.jl

Appearance

When variance carries signal, Part 1: location–scale ​

Status — Stable

Mirrors drmTMB's When variance carries signal, Part 1. In DRModels.jl today: the univariate Gaussian location–scale model fits by ML. Continue with Part 2: location–scale–scale when the question moves from the residual SD to the SD of a random effect.

Sometimes the science is in the spread, not the mean. A classic case: two conditions with the same average response but different variability. A mean-only model is blind to that; a location–scale model sees it.

A variance that moves while the mean stays put ​

julia
using DRModels, Random
Random.seed!(7)

n = 400
group = rand(0:1, n)
# identical mean (0), but the residual SD is larger in group 1:
y = 0.0 .+ exp.(-0.2 .+ 0.8 .* group) .* randn(n)
dat = (; y, group = Float64.(group))

fit = drm(bf(@formula(y ~ 1), @formula(sigma ~ group)), Gaussian(); data = dat)
coef(fit, :sigma)        # (Intercept), group — on log σ
2-element Vector{Float64}:
 -0.23202744783796148
  0.8594672923131855

The group coefficient on log σ is ≈ 0.8. On the response scale that is a residual-SD ratio:

julia
exp(coef(fit, :sigma)[2])     # ≈ how many times larger the SD is in group 1
exp(2 * coef(fit, :sigma)[2]) # the residual-VARIANCE ratio
5.57858179019725

Name the scale

σ coefficients are on log σ. Effects on log σ² are 2× the effects on log σ — quote the scale explicitly when comparing to DHGLM / O'Dea-style variance models.

Working in R? The drmTMB companion tutorial reports the same explicit SD and variance ratios with R syntax. The two packages have distinct APIs, fitting implementations, and route-specific evidence, so follow the diagnostics and inference guidance for the package you fit with.

Does the scale model earn its keep? ​

Compare the location–scale fit to a constant-σ fit by log-likelihood:

julia
fit0 = drm(bf(@formula(y ~ 1), @formula(sigma ~ 1)), Gaussian(); data = dat)
loglik(fit) - loglik(fit0)    # gain from letting σ depend on group
66.991378271033

The gain measures how much better the more flexible model fits these data; it does not by itself establish a biological mechanism. Compare the models with lrtest or aicc, then interpret the estimated change in spread and its uncertainty with confint (Wald or profile). See Prediction, residuals & model comparison.

Random dispersion: a scale that varies by group ​

When the spread itself varies across many groups (sites, individuals, studies), put a random intercept on σ rather than a fixed level per group. Because the random effect enters σ nonlinearly there is no closed-form marginal, so DRModels.jl integrates each group's effect out with per-group Gauss–Hermite quadrature. drmTMB uses Laplace: both target the same marginal model, but their numerical approximations need not be identical. The default (marginal = :LA) is non-adaptive 32-node quadrature: close to exact for small groups, but it loses accuracy for groups of hundreds of rows, where Laplace does well. Pass marginal = :Laplace for large groups or to match drmTMB; the capability page gives the measured differences.

julia
Random.seed!(13)
G = 40; m = 25; ng = G * m
grp = repeat(1:G, inner = m)
bg = 0.5 .* randn(G)                    # group-level log-σ deviations, SD 0.5
yr = exp.(log(0.5) .+ bg[grp]) .* randn(ng)
datre = (; y = yr, grp)

fitre = drm(bf(@formula(y ~ 1), @formula(sigma ~ 1 + (1 | grp))), Gaussian(); data = datre)
re_sd(fitre)[:grp_logsigma]  # recovered log-σ group-effect SD (≈ 0.5)
0.4556081755512807

re_sd returns the scale-RE SD; coef(fitre, :sigma) is the population (group-average) log σ, that is E[log σ], not log E[σ]. With a normal random effect of SD τ on log σ, the two differ by τ²/2.

See also ​