When variance carries signal, Part 2: location-scale-scale models

A location-scale model asks whether predictors change the expected response mu and the residual standard deviation sigma. A location-scale-scale model adds a third submodel: predictors can also change the standard deviation of a latent random effect. drmTMB writes this third submodel as sd(group) ~ predictors.

This is Part 2 of the variance sequence. Return to Part 1: location-scale models if mu and sigma are new to you or if the scientific question concerns only mean and residual variation. Read Which scale are you modelling? when you need to distinguish residual SD, group-level SD, known sampling variance, and likelihood weights.

Personality, predictability, and repeatability

Suppose an exploration score is recorded repeatedly for each individual. We want to ask whether sex predicts three different features of behaviour:

the mean exploration score;
between-individual variation in expected scores; and
within-individual variation around each expected score.

For observation \(j\) from individual \(i\), write

\[ \begin{aligned} y_{ij} &\sim \operatorname{Normal}(\mu_{ij}, \sigma_{e,i}^2),\\ \mu_{ij} &= \beta_0 + \beta_1\operatorname{sex}_i + b_i,\\ b_i &\sim \operatorname{Normal}(0, \sigma_{b,i}^2),\\ \log(\sigma_{e,i}) &= \gamma_0 + \gamma_1\operatorname{sex}_i,\\ \log(\sigma_{b,i}) &= \alpha_0 + \alpha_1\operatorname{sex}_i. \end{aligned} \]

Here, \(\operatorname{sex}_i=0\) denotes female and \(\operatorname{sex}_i=1\) denotes male. \(\sigma_{b,i}\) is the between-individual SD and \(\sigma_{e,i}\) is the within-individual residual SD. The matching drmTMB formula is

bf(
  exploration_score ~ sex + (1 | individual),
  sigma ~ sex,
  sd(individual) ~ sex
)

The same predictor appears in three formulas, but its coefficients answer three separate questions. Sex must be constant within each individual because it is used to model sd(individual).

If drmTMB() reports that an sd(individual) predictor varies within an individual, check the grouping variable and source data. Correct miscoded values or remove that predictor from the sd() formula. Do not average a genuinely within-individual predictor to silence the error, because doing so changes the scientific question. See Errors, warnings, and convergence for the next diagnostic steps.

Simulate and fit the three submodels

The example below gives females and males different means, between-individual SDs, and within-individual SDs. It is intentionally simple enough that the three sources of signal remain visible.

set.seed(20260715)
n_individual <- 80L
n_each <- 6L

individual_info <- data.frame(
  individual = factor(seq_len(n_individual)),
  sex = factor(
    rep(c("female", "male"), each = n_individual / 2),
    levels = c("female", "male")
  )
)

mean_by_sex <- c(female = 0.35, male = 0.70)
between_sd_by_sex <- c(female = 0.65, male = 0.40)
within_sd_by_sex <- c(female = 0.35, male = 0.60)

individual_effect <- stats::rnorm(
  n_individual,
  sd = between_sd_by_sex[individual_info$sex]
)

personality <- individual_info[rep(seq_len(n_individual), each = n_each), ]
personality$exploration_score <-
  mean_by_sex[personality$sex] +
  individual_effect[as.integer(personality$individual)] +
  stats::rnorm(
    nrow(personality),
    sd = within_sd_by_sex[personality$sex]
  )

fit_personality <- drmTMB(
  bf(
    exploration_score ~ sex + (1 | individual),
    sigma ~ sex,
    sd(individual) ~ sex
  ),
  family = gaussian(),
  data = personality
)

check_drm(fit_personality)
#> <drm_check: 15 checks>
#> ok: 14; notes: 1; warnings: 0; errors: 0
#>                         check status
#>         optimizer_convergence     ok
#>              optimizer_budget     ok
#>              finite_objective     ok
#>         logsigma_clamp_active     ok
#>                fixed_gradient     ok
#>               sdreport_status     ok
#>     hessian_positive_definite     ok
#>        standard_errors_finite     ok
#>      standard_errors_inflated     ok
#>                  dropped_rows     ok
#>                positive_scale     ok
#>     random_effect_sd_boundary     ok
#>    interval_reliability_scope   note
#>      fixed_effect_design_size     ok
#>  mu_random_effect_replication     ok
#>                                                                                   value
#>                                                                                       0
#>                                                 iterations=27; function=39; gradient=28
#>                                                                                   424.2
#>                                                                                    <NA>
#>                                                  max=0.0005111; component=beta_sigma[2]
#>                                                                                      ok
#>                                                                                    TRUE
#>                                                                  range=[0.05000,0.1961]
#>                                           n_inflated=0; max_se=0.1961; median_se=0.1069
#>                                                                     nobs=480; dropped=0
#>                                                                              min=0.3689
#>                   min=0.4393; boundary=0.0001000; term=sd(individual).sd(individual):41
#>                                                          sd_targets=80; assessed_here=0
#>  total_mb=0.08308; max_cols=2; largest=mu; largest_class=matrix; largest_density=0.7500
#>                                                                      (1 | individual)=6
#>                                                                                                                                                                                                                                                                                                                                    message
#>                                                                                                                                                                                                                                                                                                              nlminb convergence code is 0.
#>                                                                                                                                                                                                                                                        Optimizer evaluation counts recorded; no eval.max or iter.max control was supplied.
#>                                                                                                                                                                                                                                                                                                   Objective and log-likelihood are finite.
#>                                                                                                                                                                                                                                                                                         The log(sigma) clamp is not active at the optimum.
#>                                                                                                                                                                                                                                                           Maximum absolute fixed gradient is <= 0.001; largest component is beta_sigma[2].
#>                                                                                                                                                                                                                                                                                                    TMB::sdreport() completed successfully.
#>                                                                                                                                                                                                                                                                                              sdreport reports a positive-definite Hessian.
#>                                                                                                                                                                                                                                                                                               All fixed-effect standard errors are finite.
#>                                                                                                                                                                                                                                                                         No fixed-effect standard error is inflated relative to the others.
#>                                                                                                                                                                                                                                                                         No rows were dropped by model-frame or known-covariance filtering.
#>                                                                                                                                                                                                                                                                                           All fitted scale values are finite and positive.
#>                                                                                                                                                                                                               All fitted random-effect standard deviations are finite, positive, and above the requested lower-boundary warning threshold.
#>  This fit has 80 random-effect standard-deviation targets. `check_drm()` assesses the fit, not interval reliability: a target can pass every check above and still return an interval that `confint()` warns about at a variance boundary. Before reporting an interval, call `confint()` and read `conf.status` and any boundary warning.
#>                                                                                                                                                                                                                                                                                Dense fixed-effect design matrices are modest for this fit.
#>                                                                                                                                                                                                                                                                            Every random-effect level has at least two fitted observations.
round(coef(fit_personality, "mu"), 3)
#> (Intercept)     sexmale 
#>       0.353       0.355
round(coef(fit_personality, "sigma"), 3)
#> (Intercept)     sexmale 
#>      -0.997       0.587
round(coef(fit_personality, "sd(individual)"), 3)
#> (Intercept)     sexmale 
#>      -0.549      -0.273

Inspect interval availability before requesting an interval. The two slope rows below are coefficients on log-SD scales; profile_ready says whether the fitted object retained what a direct profile interval needs.

personality_targets <- profile_targets(fit_personality)
personality_targets[
  personality_targets$parm %in% c(
    "fixef:sigma:sexmale",
    "fixef:sd(individual):sexmale"
  ),
  c("parm", "estimate", "profile_ready", "profile_note")
]
#>                           parm   estimate profile_ready profile_note
#> 4          fixef:sigma:sexmale  0.5872034          TRUE        ready
#> 6 fixef:sd(individual):sexmale -0.2730528          TRUE        ready

For a quick first interval, request the named coefficients and retain the status and method columns returned with the endpoints:

personality_ci <- confint(
  fit_personality,
  parm = c(
    "fixef:sigma:sexmale",
    "fixef:sd(individual):sexmale"
  ),
  method = "wald"
)
personality_ci
#>                           parm level      lower     upper scale
#> 1          fixef:sigma:sexmale  0.95  0.4486130 0.7257939  link
#> 2 fixef:sd(individual):sexmale  0.95 -0.6573551 0.1112496  link
#>     transformation tmb_parameter index method profile.engine conf.status
#> 1 linear_predictor    beta_sigma     2   wald           <NA>        wald
#> 2 linear_predictor    beta_sd_mu     2   wald           <NA>        wald
#>   profile.boundary profile.message
#> 1               NA            <NA>
#> 2               NA            <NA>

The sigma and sd(individual) coefficients use log-SD links. For example, exponentiating the male coefficient from sigma gives the male-to-female ratio of within-individual SDs. Exponentiating the corresponding sd(individual) coefficient gives the male-to-female ratio of between-individual SDs.

When an interval row has a successful conf.status, exponentiating both log-SD endpoints gives an interval for the corresponding male-to-female SD ratio. Keep the returned interval method and status with the reported result. If the status marks a boundary or failure, do not transform missing or unreliable endpoints; inspect the profile target and follow the suggested profile or diagnostic route instead.

interval_ok <- personality_ci$conf.status == "wald" &
  stats::complete.cases(personality_ci[c("lower", "upper")])
data.frame(
  parm = personality_ci$parm[interval_ok],
  sd_ratio_lower = exp(personality_ci$lower[interval_ok]),
  sd_ratio_upper = exp(personality_ci$upper[interval_ok]),
  method = personality_ci$method[interval_ok],
  conf.status = personality_ci$conf.status[interval_ok]
)
#>                           parm sd_ratio_lower sd_ratio_upper method conf.status
#> 1          fixef:sigma:sexmale      1.5661385       2.066371   wald        wald
#> 2 fixef:sd(individual):sexmale      0.5182201       1.117674   wald        wald

Before reporting either contrast, inspect the diagnostics printed by check_drm(). A warning about weak replication, a non-positive-definite Hessian, or a large terminal gradient means the fitted SD surface may be unstable. Simplify the SD formula or increase group-level information, refit, and use the errors, warnings, and convergence guide rather than interpreting the coefficient anyway.

Two panels show repeated exploration scores for female and male individuals. Grey observations cluster around blue individual-mean ticks, while a vermillion horizontal line marks the fitted sex-specific mean. The female panel has wider spread among individual means, while the male panel has more scatter within individuals.

Repeated exploration scores for females and males. Grey points are observations, blue ticks are individual means, and the vermillion line is the fitted sex-specific mean. The spread of individual means represents between-individual variation; scatter around each individual mean represents within-individual variation. Both panels use the same vertical scale.

Read the fitted scales and calculate repeatability

Predictions with newdata are population-level values. They give the three fitted quantities for each sex directly on their natural scales.

sex_grid$mean_score <- predict(
  fit_personality, newdata = sex_grid, dpar = "mu"
)
sex_grid$between_individual_sd <- predict(
  fit_personality, newdata = sex_grid, dpar = "sd(individual)"
)
sex_grid$within_individual_sd <- predict(
  fit_personality, newdata = sex_grid, dpar = "sigma"
)

sex_grid$repeatability <- with(
  sex_grid,
  between_individual_sd^2 /
    (between_individual_sd^2 + within_individual_sd^2)
)

repeatability_table <- sex_grid[c(
  "sex",
  "mean_score",
  "between_individual_sd",
  "within_individual_sd",
  "repeatability"
)]
repeatability_table[-1] <- lapply(
  repeatability_table[-1],
  round,
  digits = 3
)
repeatability_table
#>      sex mean_score between_individual_sd within_individual_sd repeatability
#> 1 female      0.353                 0.577                0.369         0.710
#> 2   male      0.708                 0.439                0.664         0.305

In behavioural ecology, this intraclass correlation is usually called repeatability. For sex \(s\),

\[ R_s = \frac{\sigma_{b,s}^2} {\sigma_{b,s}^2 + \sigma_{e,s}^2}. \]

Repeatability is derived from the two fitted scale submodels; it is not a fourth formula. The calculation above is a point estimate. It does not supply an uncertainty interval for this nonlinear ratio.

Do not create a repeatability interval by combining the separate endpoints of the two SD intervals. That ignores their covariance and does not produce an interval for the ratio itself. Report the repeatability calculation as a point estimate unless the analysis includes a validated joint uncertainty method.

Three small panels compare fitted female and male values. The expected score and within-individual standard deviation are higher for males, while the between-individual standard deviation is higher for females. No uncertainty intervals are shown.

Model-implied values for females and males. Panels show the expected exploration score, between-individual SD from sd(individual), and within-individual residual SD from sigma. Lines aid comparison and are not uncertainty intervals.

A short phylogenetic extension

The same three-part logic applies when the latent location effect follows a phylogeny. Start from the familiar constant-scale model

\[ \begin{aligned} y_i &= \mu_i + a_i + e_i,\\ \mu_i &= \beta_0 + \beta_1 T_i,\\ \mathbf a &\sim \operatorname{MVN}(\mathbf 0, \sigma_a^2 A),\\ e_i &\sim \operatorname{Normal}(0, \sigma_e^2), \end{aligned} \]

where \(A\) is the phylogenetic correlation matrix and \(T_i\) is temperature. A location-scale-scale extension lets both SDs change linearly with temperature:

\[ \log(\sigma_{e,i}) = \gamma_0 + \gamma_1 T_i, \qquad \log(\sigma_{a,i}) = \alpha_0 + \alpha_1 T_i. \]

The matching syntax is deliberately simple:

fit_phylo_lss <- drmTMB(
  bf(
    trait ~ temperature + phylo(1 | species, tree = tree),
    sigma ~ temperature,
    sd(species, level = "phylogenetic") ~ temperature
  ),
  family = gaussian(),
  data = dat
)

Formula	Model quantity	Interpretation
`trait ~ temperature + phylo(...)`	\(\mu_i\) and \(a_i\)	expected trait and phylogenetically correlated location deviation
`sigma ~ temperature`	\(\sigma_{e,i}\)	independent SD
`sd(species, level = "phylogenetic") ~ temperature`	\(\sigma_{a,i}\)	SD of the phylogenetic location deviation

The scalar covariance \(\sigma_a^2 A\) is the constant-SD starting point. The third formula generalizes it by allowing the phylogenetic-effect SD to vary among species with temperature. The older formula spelling sd_phylo(species) ~ temperature is soft-deprecated; fitted objects retain the output label sd_phylo(species) for extractor compatibility.

With repeated observations within species, sigma is the within-species residual SD. With one response row per species, it is the independent non-phylogenetic species-level deviation. It should not automatically be described as measurement error.

What this example supports

Ordinary sd(group) ~ predictors is implemented for distinct unlabelled Gaussian mu random intercepts. The matching (1 | group) term must appear in the location formula, and scale predictors must be constant within group.
sd(group, level = "phylogenetic") ~ predictors targets the location phylogenetic effect introduced by phylo() in mu; predictors must be constant within species.
This article teaches the Gaussian route. Check narrow non-Gaussian cells in What can I fit today? rather than generalizing from it.
Random effects on the right-hand side of an sd() formula and generic direct-SD levels for spatial(), animal(), and relmat() remain separate implementation and validation questions.
Use check_drm(), inspect the terminal gradient and Hessian diagnostics, and retain failed fits. A clean optimizer code alone does not establish that every SD surface is well identified.