When variance carries signal, Part 2: location–scale–scale
Status — Experimental
Mirrors drmTMB's location-scale-scale vignette. In DRModels.jl today: sd(group) ~ z on the iid (1 | g) random effect (ML + REML), sd(group, phylogenetic) ~ z on the per-species phylogenetic SD (ML + REML, with dense and
A location–scale model asks whether predictors change the expected response μ and the residual SD σ. A location–scale–scale model adds a third submodel: predictors can also change the standard deviation of a latent random effect. DRModels.jl writes that third submodel as sd(group) ~ z, exactly as drmTMB.
Personality, predictability, repeatability
Suppose an exploration score is recorded repeatedly for each individual, and sex may predict three different things at once:
the mean score;
between-individual variation (how much individuals differ from each other);
within-individual variation (how repeatable each individual is).
For observation j of individual i:
The matching formula bundle — one formula per submodel:
using DRModels, Random
rng = Random.MersenneTwister(20260715)
n_id, n_each = 80, 6
sex = repeat([0.0, 1.0], inner = n_id ÷ 2) # 0 = female, 1 = male
b = randn(rng, n_id) .* [0.65, 0.40][Int.(sex) .+ 1] # between-SD differs by sex
id = repeat(1:n_id, inner = n_each)
sexl = sex[id]
y = [0.35, 0.70][Int.(sexl) .+ 1] .+ b[id] .+
randn(rng, n_id * n_each) .* [0.35, 0.60][Int.(sexl) .+ 1]
dat = (; y, sex = sexl, id)
fit = drm(bf(@formula(y ~ sex + (1 | id)),
@formula(sigma ~ sex),
@formula(sd(id) ~ sex)), Gaussian(); data = dat)
coef(fit, :sd) # log between-individual SD: (Intercept), sex2-element Vector{Float64}:
-0.3018222890807549
-0.753970804899395The same predictor appears in three formulas, and its three coefficients answer three separate questions. All SD submodels use a log link, so read them back through exp. On this simulation (truth: between-SD 0.65 vs 0.40, within-SD 0.35 vs 0.60):
(between_sd = exp.([1 0; 1 1] * coef(fit, :sd)), # ≈ 0.65 (F), 0.40 (M)
within_sd = exp.([1 0; 1 1] * coef(fit, :sigma))) # ≈ 0.35 (F), 0.60 (M)(between_sd = [0.7394694650119487, 0.34791638746540426], within_sd = [0.35818411715179976, 0.6154965068843709])sd() predictors must be constant within each group
Sex is used to model sd(id), so it must not vary within an individual. DRModels.jl checks this and errors, naming the offending predictor — exactly as drmTMB does. Do not average a genuinely within-group predictor to silence the error; that changes the scientific question.
What "repeatability" means here
The section title promises a repeatability, so the estimand has to be written down. Repeatability is the share of the variance of a single observation that is between-individual, i.e. the correlation between two observations of the same individual. Here both variance components depend on sex, so it is a function of the covariates, not a number:
where z_i is the covariate vector of individual i (the sd(id) and sigma formulas may use different predictors; the formula above uses the same z for both). No single scalar is "the" repeatability of such a model: any one number is a choice of covariate distribution, and the ratio of the average variances is not the average of the ratios. With sd(id) ~ sex and sigma ~ sex, report R by sex, never one pooled value. This is why repeatability(fit), icc(fit), heritability(fit), re_sd and vc all refuse these fits (and the refusal message states the formula above) instead of silently choosing a covariate value.
Ask for the conditional repeatability at the covariate values you care about:
r = repeatability(fit, (; sex = [0.0, 1.0])) # R for females, R for males
(sex = ["female", "male"], R = round.(r.estimate; digits = 3),
lower = round.(r.lower; digits = 3), upper = round.(r.upper; digits = 3))(sex = ["female", "male"], R = [0.81, 0.242], lower = [0.721, 0.133], upper = [0.875, 0.399])On this simulation the truth is R = 0.65² / (0.65² + 0.35²) = 0.775 for females and 0.40² / (0.40² + 0.60²) = 0.308 for males. logit R(z) = 2(α'z − γ'z) is linear in the coefficients, so the interval is an exact-linear Wald interval on the logit scale (it uses the joint vcov of the sd and sigma blocks) mapped back to (0, 1). newdata must contain every predictor of the sd(id) and sigma formulas; a categorical predictor must contain all of its training levels.
Wald and profile intervals for the coefficients themselves target the block as usual:
confint(fit; parm = :sd)2-element Vector{@NamedTuple{param::Symbol, coef::String, estimate::Float64, lower::Float64, upper::Float64}}:
(param = :sd, coef = "(Intercept)", estimate = -0.3018222890807549, lower = -0.529554047576674, upper = -0.07409053058483589)
(param = :sd, coef = "sex", estimate = -0.753970804899395, lower = -1.1609744399585336, upper = -0.34696716984025666)method = :REML is supported on this route.
The phylogenetic scale: sd(species, phylogenetic)
For comparative data the deeper question is usually about the phylogenetic SD: does among-species variation that tracks the tree change along an environmental gradient? With one trait value per species i:
where A is the Brownian-motion phylogenetic correlation, D_a = \mathrm{Diagonal}(\sigma_{a,1},\dots,\sigma_{a,n}), and each of log σ_a and log σ_e carries its own linear predictor. This is the location–scale–scale framework of the Mizuno et al. ecogeographical-rules protocol: climate in the mean, the non-phylogenetic scale, and the phylogenetic scale.
using LinearAlgebra
# a balanced 64-tip unit-height tree
function _baln(d)
node(p, k) = k == 0 ? "$(p):$(1/d)" :
(k == d ? "($(node(p*"a",k-1)),$(node(p*"b",k-1)));" :
"($(node(p*"a",k-1)),$(node(p*"b",k-1))):$(1/d)")
node("t", d)
end
phy = DRModels.augmented_phy(_baln(6))
G = phy.n_leaves
K0 = DRModels.sigma_phy_dense(phy; σ²_phy = 1.0)
dK = sqrt.(diag(K0)); K = K0 ./ (dK * dK')
rng2 = Random.MersenneTwister(11)
x = randn(rng2, G) # standardised "climate"
sda = exp.(-0.5 .+ 0.4 .* x) # phylo SD rises with climate
sde = exp.(-1.0 .- 0.3 .* x) # residual SD falls with climate
a = Diagonal(sda) * (cholesky(Symmetric(K)).L * randn(rng2, G))
yq = 1.0 .+ 0.5 .* x .+ a .+ sde .* randn(rng2, G)
datq = (y = yq, x = x, species = String.(phy.leaf_names))
fitq = drm(bf(@formula(y ~ x + phylo(1 | species)),
@formula(sigma ~ x),
@formula(sd(species, phylogenetic) ~ x)),
Gaussian(); data = datq, tree = phy)
(mu = coef(fitq, :mu), sigma = coef(fitq, :sigma), sd_phylo = coef(fitq, :sd_phylo))(mu = [0.7395755054961157, 0.47023823709398654], sigma = [-1.111655922832119, -0.4902067312323761], sd_phylo = [-1.0939684420636404, 0.6095363889181672])The estimates track the simulated truth (mean 1.0 and 0.5; the σ and σ_a slopes in the right directions). For these same data, drmTMB returns the same log-likelihood (−69.1373) and the same coefficients to seven significant figures.
When are σ_a and σ_e separately identified?
With one trait value per species there are no replicates, so the two variance surfaces are not separated by repeated measures (the personality example above has six observations per individual). They are separated only by the off-diagonal structure of A: the residual σ_e is tip-private, whereas σ_a induces covariance between related tips. For intercept-only scales V = σ_a² A + σ_e² I, and the Fisher information for (σ_a², σ_e²) is
which is singular exactly when A ∝ I (a star tree): then only σ_a² h + σ_e² is identified and any split is as good as any other. Measured on random coalescent trees (median over 20 trees per cell, true residual share 0.3, asymptotic Fisher standard errors): the correlation between the two variance estimates is −0.35 at 16 tips and −0.26 at 64 tips, and the SE of the residual share is 0.24 and 0.15. Adding a second observation per species lowers these to 0.17 and 0.11. As the tree approaches a star (mean off-diagonal tip correlation 0.1) the correlation reaches −0.96 and the SE of the share is 1.2 at 16 tips (0.65 at 64): the split is not identified. In simulation with 16 tips and a true 50 : 50 split, σ_a was estimated at zero in 51 % of data sets (64 tips: 17 %); with a true share of 0.9 residual, 73 % (64 tips: 58 %).
What the software does about it:
Both scales may depend on covariates (the three-submodel route is not removed), but separating them leans entirely on
A; expect wide intervals and useconfint(fit; method = :profile)or the bootstrap rather than Wald standard errors.When a fit with one observation per group lands on a boundary,
drmwarns and names the cause (the residual variance share ≈ 1, or the structured SD collapsed to 0, with the one-observation-per-group note), andcheck_drm(fit).variance_boundarycarries the same flags. An exactly or nearly singular information matrix (a star tree) is reported by the existing "Hessian is numerically singular" warning.The flag is intentionally silent on healthy interior fits. A warning on every one-observation-per-species fit would fire on every ordinary comparative model; the data-driven boundary check fires only when the split has actually collapsed.
Conditional on the covariates, the phylogenetic share of a tip's variance is h²(z) = σ_a(z)² / (σ_a(z)² + σ_e(z)²); as for repeatability above it is a function of z, and heritability(fit, newdata) returns it with a Wald-on-logit interval:
heritability(fitq, (; x = [-1.0, 0.0, 1.0]))(estimate = [0.10301992074861772, 0.5088428182538217, 0.9033366699465246], se_logit = [1.4973025267147162, 1.3582237853975985, 1.39820443531627], lower = [0.006067212543301157, 0.06744002755480843, 0.37623279687812017], upper = [0.6836393706200102, 0.9368756909082099, 0.9931408674972065], level = 0.95, method = :wald_logit)Species rows need not follow tree-tip order when fitting an LSS model. String labels match phy.leaf_names exactly; integer labels are positions in 1:G, not arbitrary group identifiers. Every tip must be represented in the full input, even if all responses for a tip are missing. Scale predictors must still be constant within each species.
row_order = randperm(rng2, G)
datq_shuffled = (; y = datq.y[row_order], x = datq.x[row_order],
species = datq.species[row_order])
fitq_shuffled = drm(bf(@formula(y ~ x + phylo(1 | species)),
@formula(sigma ~ x),
@formula(sd(species, phylogenetic) ~ x)),
Gaussian(); data = datq_shuffled, tree = phy)
@assert isapprox(loglik(fitq_shuffled), loglik(fitq); atol = 1e-7, rtol = 0)
@assert isapprox(coef(fitq_shuffled, :sd_phylo), coef(fitq, :sd_phylo);
atol = 4e-6, rtol = 0)
coef(fitq_shuffled, :sd_phylo)2-element Vector{Float64}:
-1.0939684420636433
0.6095363889181686Bootstrap draws follow the fitted variance model
For Gaussian models with sd() submodels, each bootstrap replicate redraws every IID and phylogenetic mean effect independently, using tree-tip names to match species. It also redraws residual error. Missing responses remain missing in each refit, and a REML seed fit is refitted with REML. These simulation checks do not establish interval coverage or large-tree performance; inspect failed-refit counts before interpreting intervals.
Grammar and solver notes
sd(species, phylogenetic)is the canonical spelling (drmTMB:sd(species, level = "phylogenetic");@formulacannot parse keyword arguments, so the level is a bare symbol). The legacysd_phylo(species)still works but is soft-deprecated on both sides.The mean formula must carry the matching
phylo(1 | species)marker, andtree = …is required.ML and REML: both
method = :ML(default) andmethod = :REMLare supported on phylogenetic and iid LSS routes.Sparse whole-tree scaling: for large trees,
sparse = true(oralgorithm = :sparse_lbfgs) invokes the exactaugmented-state GMRF engine with Takahashi selected-inverse leaf variances, avoiding dense matrix storage and dense factorization. This sparse engine is selected automatically when
Missing response handling
Like other Gaussian routes in DRModels.jl, Location–Scale–Scale models support incomplete responses (missing or NaN in y), matching response = "include" in drmTMB:
# Mask some responses
y_miss = Vector{Union{Float64, Missing}}(copy(dat.y))
y_miss[1:5] .= missing
dat_miss = (; y = y_miss, sex = dat.sex, id = dat.id)
fit_miss = drm(bf(@formula(y ~ sex + (1 | id)),
@formula(sigma ~ sex),
@formula(sd(id) ~ sex)),
Gaussian(); data = dat_miss)
nobs(fit_miss) # count of observed rows475Group-level designs (
REML on location–scale–scale models
REML accounts for estimating fixed effects when fitting covariance parameters:
fit_reml = drm(bf(@formula(y ~ sex + (1 | id)),
@formula(sigma ~ sex),
@formula(sd(id) ~ sex)),
Gaussian(); data = dat, method = :REML)
reml_loglik(fit_reml) # restricted log-likelihood-408.179000628769Both iid and phylogenetic LSS models support REML estimation. Standard errors can be unreliable when a variance approaches zero; a successful fit alone does not establish reliable uncertainty estimates.
To assess the bootstrap rather than only its interval endpoints, retain the full result:
boot = bootstrap_result(fit_reml; data = dat, B = 4,
rng = MersenneTwister(20260830), threads = true,
failures = :skip, check_converged = true)
@assert boot.attempted == boot.used + boot.failed
(attempted = boot.attempted, used = boot.used, failed = boot.failed)(attempted = 4, used = 4, failed = 0)Four replicates keep this example quick; they are too few for an interval you would report. Choose the number of replicates for your analysis and retain boot.failures, which records failed replicate seeds and messages. For a phylogenetic model, also pass the original tree. For profiles, inspect profile_result(...).stats: the generic path records nuisance termination and failed endpoints separately from limits of the searched range. These checks do not establish interval coverage or practical performance on large trees.
From R, with intervals
The same models run from R through engine = "julia", with Wald, profile, and parametric-bootstrap intervals:
fit <- drmTMB(
bf(y ~ temp + I(temp^2) + prec + I(prec^2) + phylo(1 | species, tree = tr),
sigma ~ temp + prec,
sd(species, level = "phylogenetic") ~ temp + prec), # the "M6" shape
family = gaussian(), data = d, engine = "julia"
)
confint(fit, parm = "fixef:sd_phylo:temp", method = "profile")
confint(fit, parm = "fixef:sd_phylo:temp", method = "bootstrap", R = 199,
threads = TRUE) # parallel only if JULIA_NUM_THREADS > 1threads = TRUE does not start threads by itself. Set JULIA_NUM_THREADS (for example to 4) before the first Julia call in the R session; otherwise the refits run serially and the result reports julia.threaded as FALSE. See drmTMB's Using the Julia engine.
For the ecogeographical-rules formulas assessed here, the Julia and TMB fits gave identical log likelihoods.
See also
When variance carries signal — Part 1:
μandσ.Phylogenetic structured effects — the
phylo()marker and tree handling.API stability —
sd/sd_phylosit on the Experimental tier while the surface settles.