---
title: "When variance carries signal, Part 2: location-scale-scale models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{When variance carries signal, Part 2: location-scale-scale models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.1,
  dpi = 144
)
if (!"package:drmTMB" %in% search()) {
  library(drmTMB)
}
```

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](location-scale.html) 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?](which-scale.html) 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:

1. the mean exploration score;
2. between-individual variation in expected scores; and
3. 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

```r
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](convergence.html) 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.

```{r personality-simulate}
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]
  )
```

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

check_drm(fit_personality)
round(coef(fit_personality, "mu"), 3)
round(coef(fit_personality, "sigma"), 3)
round(coef(fit_personality, "sd(individual)"), 3)
```

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.

```{r personality-interval-targets}
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")
]
```

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

```{r personality-component-intervals}
personality_ci <- confint(
  fit_personality,
  parm = c(
    "fixef:sigma:sexmale",
    "fixef:sd(individual):sexmale"
  ),
  method = "wald"
)
personality_ci
```

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.

```{r personality-ratio-intervals}
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]
)
```

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](convergence.html) rather
than interpreting the coefficient anyway.

```{r personality-data-figure, echo = FALSE, fig.width = 7, fig.height = 4.1, fig.cap = "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.", fig.alt = "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."}
individual_summary <- stats::aggregate(
  exploration_score ~ individual + sex,
  data = personality,
  FUN = mean
)
individual_summary <- individual_summary[
  order(individual_summary$sex, individual_summary$exploration_score),
]
individual_summary$display_id <- ave(
  individual_summary$exploration_score,
  individual_summary$sex,
  FUN = seq_along
)

plot_data <- merge(
  personality,
  individual_summary[c("individual", "sex", "display_id")],
  by = c("individual", "sex"),
  sort = FALSE
)

sex_grid <- data.frame(
  sex = factor(c("female", "male"), levels = levels(personality$sex))
)
sex_grid$fitted_mean <- predict(
  fit_personality,
  newdata = sex_grid,
  dpar = "mu"
)

ggplot2::ggplot(
  plot_data,
  ggplot2::aes(x = display_id, y = exploration_score)
) +
  ggplot2::geom_point(
    position = ggplot2::position_jitter(width = 0.08, height = 0, seed = 1),
    colour = "grey55",
    alpha = 0.42,
    size = 0.8
  ) +
  ggplot2::geom_segment(
    data = individual_summary,
    ggplot2::aes(
      x = display_id - 0.26,
      xend = display_id + 0.26,
      y = exploration_score,
      yend = exploration_score
    ),
    inherit.aes = FALSE,
    colour = "#0072B2",
    linewidth = 0.65
  ) +
  ggplot2::geom_hline(
    data = sex_grid,
    ggplot2::aes(yintercept = fitted_mean),
    inherit.aes = FALSE,
    colour = "#D55E00",
    linewidth = 0.8
  ) +
  ggplot2::geom_text(
    data = sex_grid,
    ggplot2::aes(
      x = Inf,
      y = fitted_mean,
      label = "fitted sex mean"
    ),
    inherit.aes = FALSE,
    colour = "#D55E00",
    hjust = 1.05,
    vjust = -0.6,
    size = 3
  ) +
  ggplot2::facet_wrap(
    ~sex,
    nrow = 1,
    labeller = ggplot2::as_labeller(c(female = "Female", male = "Male"))
  ) +
  ggplot2::labs(
    x = "Individuals, ordered within sex by their observed mean",
    y = "Exploration score"
  ) +
  ggplot2::theme_minimal(base_size = 11) +
  ggplot2::theme(
    panel.grid.minor = ggplot2::element_blank(),
    panel.grid.major.x = ggplot2::element_blank(),
    axis.text.x = ggplot2::element_blank(),
    axis.ticks.x = ggplot2::element_blank()
  )
```

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

```{r personality-components}
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
```

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.

```{r personality-component-figure, echo = FALSE, fig.width = 7, fig.height = 2.9, fig.cap = "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.", fig.alt = "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."}
component_surface <- rbind(
  data.frame(
    sex = sex_grid$sex,
    component = "Expected score",
    estimate = sex_grid$mean_score
  ),
  data.frame(
    sex = sex_grid$sex,
    component = "Between-individual SD",
    estimate = sex_grid$between_individual_sd
  ),
  data.frame(
    sex = sex_grid$sex,
    component = "Within-individual SD",
    estimate = sex_grid$within_individual_sd
  )
)
component_surface$component <- factor(
  component_surface$component,
  levels = c(
    "Expected score",
    "Between-individual SD",
    "Within-individual SD"
  )
)

ggplot2::ggplot(
  component_surface,
  ggplot2::aes(x = sex, y = estimate, group = 1)
) +
  ggplot2::geom_line(colour = "grey65", linewidth = 0.55) +
  ggplot2::geom_point(colour = "#0072B2", size = 2.4) +
  ggplot2::geom_text(
    ggplot2::aes(label = sprintf("%.2f", estimate)),
    vjust = -0.8,
    colour = "grey20",
    size = 3.1
  ) +
  ggplot2::facet_wrap(~component, scales = "free_y", nrow = 1) +
  ggplot2::expand_limits(y = 0) +
  ggplot2::scale_y_continuous(
    expand = ggplot2::expansion(mult = c(0.02, 0.18))
  ) +
  ggplot2::scale_x_discrete(labels = c(female = "Female", male = "Male")) +
  ggplot2::labs(x = NULL, y = "Fitted value") +
  ggplot2::theme_minimal(base_size = 11) +
  ggplot2::theme(
    panel.grid.minor = ggplot2::element_blank(),
    panel.grid.major.x = ggplot2::element_blank(),
    strip.text = ggplot2::element_text(face = "bold")
  )
```

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

```r
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?](model-map.html) 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.

Return to [Part 1: location-scale models](location-scale.html) for residual
variability without a model for the group-level SD. For broader phylogenetic
syntax, continue to [Phylogenetic mixed models](phylogenetic-models.html). For
prediction tables and the distinction between `sigma` and `sd(group)`, use
[Which scale are you modelling?](which-scale.html).

## Reference

Nakagawa, S., Mizuno, A., Williams, C., Lagisz, M., Yang, Y., and Drobniak,
S. M. (2025). Quantifying macro-evolutionary patterns of trait mean and
variance with phylogenetic location-scale models. *Methods in Ecology and
Evolution*. [doi:10.1111/2041-210X.70160](https://doi.org/10.1111/2041-210X.70160).