---
title: "Phylogenetic mixed models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Phylogenetic mixed models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Suggested packages used only inside guarded chunks. When either is absent the
# corresponding chunk degrades to a short note instead of erroring at build.
has_ape <- requireNamespace("ape", quietly = TRUE)
has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)
```

<style>
@media (max-width: 575.98px) {
  main h1 {
    hyphens: none;
    overflow-wrap: normal;
    word-break: normal;
  }

  main table {
    display: block;
    max-width: 100%;
    overflow-x: auto;
  }
}
</style>

Closely related species are not independent observations. They share ancestry,
so a trait value carried by one species is partly inherited from the same
nodes that shaped its relatives. A phylogenetic mixed model encodes that shared
history as a structured random effect: a species-level deviation whose
covariance is read off the tree, so that sister taxa are expected to deviate
from the regression line in the same direction.

`drmTMB` exposes this through the `phylo()` term. This article explains where
the covariance comes from, how to write a `phylo()` model for a Gaussian
response, and how the same syntax carries over to non-Gaussian responses such
as counts. Everything below runs at build time on a small simulated tree, so
you can read the fitted numbers next to the values that generated them.

If you are choosing between `animal()`, `phylo()`, `spatial()`, and `relmat()`,
read the [structural-dependence overview](structural-dependence.html) first.
When observations belong to a *pair* of species drawn from two different trees,
read [two-tree phylogenetic interactions](bipartite-phylogenetic-interactions.html).

## The phylogenetic covariance

A `phylo()` term adds a vector of species-level location deviations $u$ with a
mean of zero and a covariance proportional to the phylogenetic covariance
matrix $A$,

\[
\begin{aligned}
y_{ij} &= \mathbf{x}_{ij}^\top \boldsymbol{\beta}
          + u_{s(i,j)} + \varepsilon_{ij}, \\
\mathbf{u} &\sim \mathcal{N}\!\left(\mathbf{0},\, \sigma_{\text{phylo}}^2\, A\right), \\
\varepsilon_{ij} &\sim \mathcal{N}\!\left(0,\, \sigma^2\right),
\end{aligned}
\]

where $s(i,j)$ is the species of observation $j$ in group $i$. Under a Brownian
motion model of trait evolution, $A_{kl}$ is the height of the most recent
common ancestor of species $k$ and $l$: the longer two species have shared an
evolutionary path, the more strongly their deviations covary. For an
ultrametric tree (all tips equidistant from the root), the diagonal of $A$ is
constant and the off-diagonal entries are the shared root-to-ancestor path
lengths.

Two standard deviations appear, and they answer different questions:

- $\sigma_{\text{phylo}}$ is the **phylogenetic** SD: how much species deviate
  from the fixed-effect prediction in a way that tracks the tree.
- $\sigma$ is the **residual** SD: variation among observations within a
  species, independent of ancestry.

Their ratio is the phylogenetic signal,

\[
\lambda = \frac{\sigma_{\text{phylo}}^2}{\sigma_{\text{phylo}}^2 + \sigma^2},
\]

the proportion of the random-plus-residual variance attributable to phylogeny.
`drmTMB` reports this as a derived quantity (see below). Internally `drmTMB`
does not invert the dense $A$. It uses the Hadfield and Nakagawa (2010) sparse
phylogenetic precision (an $A^{-1}$ built from the tree's branch lengths), which
keeps the augmented-state Laplace approximation fast as the number of species
grows.

## A small tree with `ape::rcoal()`

We simulate a tree small enough to fit in a fraction of a second. `ape::rcoal()`
draws a random **ultrametric** coalescent tree, which is exactly the shape
`phylo()` expects (branch lengths present, all tips contemporaneous).

```{r simulate-tree, eval = has_ape}
library(ape)
set.seed(2026)

n_species <- 16L
tree <- rcoal(n_species, tip.label = paste0("sp", seq_len(n_species)))

c(ultrametric = is.ultrametric(tree), n_tip = length(tree$tip.label))
```

The covariance $A$ implied by this tree is available from the same internal
helper the fitting code uses. We only need it here to *simulate* a trait with
known phylogenetic structure; you never compute it by hand for a real fit.

```{r tip-covariance, eval = has_ape}
A <- drmTMB:::drm_phylo_tip_covariance(tree)
dim(A)
round(A[1:4, 1:4], 3)
```

## Simulating a Gaussian trait with phylogenetic signal

We draw one species-level deviation per tip from $\mathcal{N}(0,
\sigma_{\text{phylo}}^2 A)$ using a Cholesky factor of $A$, add a fixed
covariate effect, and add independent residual noise. Several observations per
species let the model separate the phylogenetic SD from the residual SD.

```{r simulate-gaussian, eval = has_ape}
sd_phylo_true <- 0.8     # phylogenetic SD
sigma_true    <- 0.3     # residual SD
n_per_species <- 6L

# One deviation per species, correlated along the tree.
u <- as.vector(t(chol(A)) %*% rnorm(n_species, sd = sd_phylo_true))
names(u) <- tree$tip.label

species <- rep(tree$tip.label, each = n_per_species)
x <- rnorm(length(species))

# y = intercept + slope * x + phylogenetic deviation + residual noise
trait <- 0.5 - 0.4 * x + u[species] + rnorm(length(species), sd = sigma_true)

dat <- data.frame(
  trait = unname(trait),
  x = x,
  species = species
)
head(dat)
```

## Fitting with `phylo(1 | species)`

The `phylo()` term goes inside the mean (`mu`) formula. It takes a random-effect
specification, `1 | species`, and the tree as `tree = tree`. Wrap the formulas
in `bf()` (the `drmTMB` formula builder) and pass a family, exactly as for any
other `drmTMB` model.

```{r fit-gaussian, eval = has_ape}
fit <- drmTMB(
  bf(trait ~ x + phylo(1 | species, tree = tree),
     sigma ~ 1),
  family = gaussian(),
  data = dat
)
check_drm(fit)
```

Start with `check_drm()`: it reports convergence, Hessian, scale, and
phylogenetic-replication diagnostics in one public table. Then use the public
coefficient and target extractors to see the fixed slope and the two SDs:

```{r gaussian-pieces, eval = has_ape}
coef(fit, "mu")

sd_targets <- profile_targets(fit)
sd_targets[
  sd_targets$parm %in% c("sigma", "sd:mu:phylo(1 | species)"),
  c("parm", "estimate", "scale", "profile_ready", "profile_note")
]
```

The recovered slope (`x`) is close to its true value of `-0.4`, and the
phylogenetic and residual SDs are in the neighbourhood of `0.8` and `0.3`. The
fixed intercept need not match the simulated `0.5`: with only `r n_species`
species the species-level deviations carry a non-zero sample mean that the
random effect absorbs, so the intercept and the deviations trade off. This is
expected behaviour, not a fitting error -- the scientifically interpretable
quantities are the slope and the two SDs.

### Reading the phylogenetic signal

`summary()` adds a derived row for the phylogenetic signal $\lambda$, alongside
the random-effect and residual variances it is built from.

```{r gaussian-signal, eval = has_ape}
summary(fit)$derived[, c(
  "quantity", "estimate",
  "random_effect_variance", "residual_variance"
)]
```

`estimate` is $\lambda$. A value well above zero says that, after accounting for
the covariate, related species really do resemble each other more than
unrelated ones. The fitted object also contains the augmented phylogenetic
state used by the sparse-precision representation. To inspect only the
species-tip deviations, match that state by the tree's tip labels:

```{r gaussian-ranef, eval = has_ape}
phylo_dev <- ranef(fit, "phylo_mu")
tip_dev <- phylo_dev$values[tree$tip.label]
head(tip_dev)
```

### Uncertainty for the two SDs

Variance components have their own `confint()` target. The interval is on the
response (SD) scale, obtained by transforming the Wald interval for the
log-SD parameter. The table contains both the residual SD and the
phylogenetic location SD.

```{r gaussian-confint, eval = has_ape}
confint(fit, parm = "variance_components")[, c(
  "parm", "lower", "upper", "scale"
)]
```

A transformed log-SD interval is necessarily positive, so its lower bound
should **not** be read as a test of a zero variance component. Use it to
describe uncertainty in the magnitude of each SD. The earlier `check_drm(fit)`
table reports diagnostics for the recognised phylogenetic layer.

The figure uses **Confidence Eyes** rather than flat interval bars. Each pale
eye is the finite 95% Wald confidence region shaped on the log-SD scale:
compatibility is greatest near its centre and tapers towards the endpoints.
The eye is a frequentist compatibility display, not a posterior density.

```{r gaussian-figure, eval = has_ape && has_ggplot2, fig.width = 5.2, fig.height = 2.2, fig.cap = "Confidence Eyes for the two response-scale SDs. Pale shapes are the default finite 95% Wald confidence regions, constructed on the log-SD scale; hollow circles are the raw fitted SDs. The default small-sample correction shifts the phylogenetic eye slightly relative to its raw estimate. The data-generating values were 0.8 and 0.3.", fig.alt = "Two Confidence Eye rows. The phylogenetic SD has a raw fitted value of 0.78 and a broad pale confidence region from 0.50 to 1.37. The residual SD has a fitted value of 0.27 and a narrow pale confidence region from 0.23 to 0.31. Hollow circles mark the fitted values."}
vc <- confint(fit, parm = "variance_components")
sd_targets <- profile_targets(fit)
target <- c("sd:mu:phylo(1 | species)", "sigma")
interval_row <- match(target, vc$parm)
target_row <- match(target, sd_targets$parm)
stopifnot(!anyNA(interval_row), !anyNA(target_row))

sd_tab <- data.frame(
  label = factor(
    c("Phylogenetic SD", "Residual SD"),
    levels = c("Residual SD", "Phylogenetic SD")
  ),
  estimate = c(
    sd_targets$estimate[target_row]
  ),
  lower = vc$lower[interval_row],
  upper = vc$upper[interval_row]
)
stopifnot(
  all(is.finite(unlist(sd_tab[c("estimate", "lower", "upper")]))),
  all(sd_tab$lower > 0),
  all(sd_tab$lower <= sd_tab$estimate),
  all(sd_tab$estimate <= sd_tab$upper)
)

sd_eye <- do.call(rbind, lapply(seq_len(nrow(sd_tab)), function(i) {
  log_lower <- log(sd_tab$lower[i])
  log_upper <- log(sd_tab$upper[i])
  log_centre <- 0.5 * (log_lower + log_upper)
  log_value <- seq(log_lower, log_upper, length.out = 401L)
  half_width <- 0.5 * (log_upper - log_lower)
  height <- pmax(1 - ((log_value - log_centre) / half_width)^2, 0)
  data.frame(
    label = as.character(sd_tab$label[i]),
    value = exp(log_value),
    height = height
  )
}))
sd_eye$label <- factor(sd_eye$label, levels = levels(sd_tab$label))
sd_eye$y <- as.numeric(sd_eye$label)
sd_tab$y <- as.numeric(sd_tab$label)

ggplot2::ggplot() +
  ggplot2::geom_ribbon(
    data = sd_eye,
    ggplot2::aes(
      x = value,
      ymin = y - 0.20 * height,
      ymax = y + 0.20 * height,
      group = label
    ),
    fill = "#0072B2",
    alpha = 0.24,
    colour = NA
  ) +
  ggplot2::geom_point(
    data = sd_tab,
    ggplot2::aes(x = estimate, y = y),
    shape = 21,
    fill = "white",
    colour = "#0072B2",
    size = 3.0,
    stroke = 1.0
  ) +
  ggplot2::scale_y_continuous(
    breaks = seq_along(levels(sd_tab$label)),
    labels = levels(sd_tab$label),
    expand = ggplot2::expansion(add = 0.38)
  ) +
  ggplot2::scale_x_continuous(
    limits = c(0, NA),
    expand = ggplot2::expansion(mult = c(0, 0.04))
  ) +
  ggplot2::labs(
    x = "Standard deviation (response scale)",
    y = NULL
  ) +
  ggplot2::theme_minimal(base_size = 12.5) +
  ggplot2::theme(
    axis.line.x = ggplot2::element_line(colour = "grey40", linewidth = 0.35),
    axis.ticks.x = ggplot2::element_line(colour = "grey40", linewidth = 0.35),
    panel.grid.major.y = ggplot2::element_blank(),
    panel.grid.minor = ggplot2::element_blank(),
    axis.text.y = ggplot2::element_text(colour = "grey15")
  )
```

## Non-Gaussian responses

The same `phylo(1 | species, tree = tree)` term works for ordinary Poisson and
negative-binomial (NB2) location models. The phylogenetic deviation now acts on
the linear predictor of the count mean through the log link, so related species
share a baseline abundance. Two differences from the Gaussian case:

- Poisson has no residual `sigma` formula. NB2 instead has a modelled
  overdispersion `sigma`; its exact q1 phylogenetic `sigma` gate accepts an
  intercept plus one independent slope at recovery grade, separately from the
  `mu` field;
- the simulated deviation is built from the **correlation** form of $A$ (unit
  diagonal) scaled by the phylogenetic SD, so the SD is interpretable on the
  log-mean scale.

```{r simulate-count, eval = has_ape}
set.seed(11)

# Standardise A to a correlation matrix, then scale by the phylogenetic SD.
A_cor <- A / outer(sqrt(diag(A)), sqrt(diag(A)))
sd_phylo_count <- 0.5
u_count <- as.vector(t(chol(A_cor)) %*% rnorm(n_species)) * sd_phylo_count
names(u_count) <- tree$tip.label

species_c <- rep(tree$tip.label, each = n_per_species)
x_c <- rep(seq(-1, 1, length.out = n_per_species), times = n_species)

eta <- log(3) - 0.3 * x_c + u_count[species_c]   # log mean
count <- rpois(length(eta), lambda = exp(eta))

dat_count <- data.frame(count = count, x = x_c, species = species_c)
range(dat_count$count)
```

A Poisson fit uses the same call with `family = poisson()`:

```{r fit-poisson, eval = has_ape}
fit_pois <- drmTMB(
  bf(count ~ x + phylo(1 | species, tree = tree)),
  family = poisson(link = "log"),
  data = dat_count
)

check_drm(fit_pois)
coef(fit_pois, "mu")     # log-mean intercept near log(3) ~ 1.10, slope near -0.3
summary(fit_pois)$parameters # phylogenetic SD on the log-mean scale
```

If the counts are overdispersed relative to a Poisson, swap in `nbinom2()`. The
NB2 family adds an overdispersion (scale) parameter while keeping the identical
phylogenetic location term.

```{r fit-nbinom2, eval = has_ape}
fit_nb <- drmTMB(
  bf(count ~ x + phylo(1 | species, tree = tree)),
  family = nbinom2(),
  data = dat_count
)

check_drm(fit_nb)
summary(fit_nb)$parameters # phylogenetic SD, NB2 mean model
```

The phylogenetic SD is reported on the log-mean (link) scale for count
families, so it is not directly comparable to the Gaussian response-scale SD;
compare it instead to other log-scale effects in the same model.

## A larger fit, for reference

The fits above are tiny by design. For intuition about runtime on a more
realistic tree, the chunk below sketches a 200-species fit. It is marked
`eval = FALSE` so the vignette never blocks on it; the numbers in the comments
are illustrative of the shape of the output, not a benchmarked claim.

```{r large-sketch, eval = FALSE}
set.seed(99)
big_tree <- ape::rcoal(200, tip.label = paste0("t", 1:200))
A_big <- drmTMB:::drm_phylo_tip_covariance(big_tree)
u_big <- as.vector(t(chol(A_big)) %*% rnorm(200, sd = 0.7))
names(u_big) <- big_tree$tip.label

sp <- rep(big_tree$tip.label, each = 4L)
xb <- rnorm(length(sp))
yb <- 0.2 + 0.5 * xb + u_big[sp] + rnorm(length(sp), sd = 0.3)
big <- data.frame(y = yb, x = xb, species = sp)

fit_big <- drmTMB(
  bf(y ~ x + phylo(1 | species, tree = big_tree), sigma ~ 1),
  family = gaussian(),
  data = big
)
summary(fit_big)$parameters
# The sparse-precision path keeps this on the order of a second on a laptop;
# cost grows roughly linearly in the number of species rather than cubically.
```

## What `phylo()` fits today

This is a reader-level summary of the current capability ledger. The tier words
matter: **diagnostic-only** means a fit/extractor smoke, **point-fit recovery**
means simulation or oracle evidence for fitted values, and **inference-ready
with caveats** means interval evidence exists only inside the stated design.
The [model map](model-map.html) gives the wider package-level boundary.

<style>
.phylo-capability-table {
  overflow-x: auto;
  -webkit-overflow-scrolling: touch;
}
.phylo-capability-table table {
  min-width: 760px;
}
</style>

<div class="phylo-capability-table" markdown="1">

| Question | Syntax | Status |
| --- | --- | --- |
| Gaussian location SD | `phylo(1 | species, tree = tree)` or one independent slope in `mu` | ML intercept: inference-ready with caveats. ML intercept-plus-slope and REML intercept: point-fit recovery. |
| Gaussian residual-scale phylogeny | `sigma ~ phylo(1 | species, tree = tree)` or one independent slope | ML intercept: point-fit recovery. ML intercept-plus-slope: inference-ready with caveats. REML pure scale route: point-fit recovery. |
| Gaussian location and residual scale together | matching labelled `phylo()` terms in `mu` and `sigma` | Constant and first one-slope q=2 blocks have ML point-fit recovery; the constant REML q=2 block also has point-fit recovery. |
| Gaussian predictor-dependent phylogenetic SD | `phylo(1 | species, tree = tree)` in `mu` plus `sd(species, level = "phylogenetic") ~ x` | Fitted direct-SD route; REML admission is included in the q=1 point-fit-recovery cell. This is not a blanket interval/coverage claim. |
| Poisson/NB2 q1 phylogenetic `mu` intercept-plus-one-slope | `phylo(1 + x | species, tree = tree)` in `mu` | The exact unlabelled intercept plus one independent slope at recovery grade has point-fit recovery. |
| NB2 overdispersion deviations | NB2 q1 phylogenetic `sigma` via `phylo(1 + x | species, tree = tree)` | The exact unlabelled intercept-plus-one-slope route has point-fit recovery; no interval or coverage promotion. |
| Gamma or lognormal location SD | `phylo(1 | species, tree = tree)` in `mu` | Exact ML intercept routes have point-fit recovery; slopes, structured `sigma`, REML, and intervals are outside those cells. |
| Beta phylogenetic location-SD regression | `phylo(1 | spp_id, tree = tree)` in `mu` plus `sd(spp_id, level = "phylogenetic") ~ x` | Exact unlabelled q=1 ML route is inference-ready with caveats only for the validated interior-response, 1,024-species, four-replicate designs. Family `sigma` remains fixed-effect and $\phi = \sigma^{-2}$. |
| Student-t tail-weight deviations | Student-t q1 phylogenetic `nu` via `phylo(1 | species, tree = tree)` | Diagnostic-only fit/extractor evidence; not point-fit recovery. |
| cumulative-logit location deviations | cumulative-logit q1 phylogenetic `mu` via `phylo(1 | species, tree = tree)` | Diagnostic-only fit/extractor evidence; not point-fit recovery. |
| Two Gaussian means sharing phylogeny | matching `phylo(1 | p | species, tree = tree)` in `mu1` and `mu2` | Intercept q=2 has ML and REML point-fit recovery; the exact ML slope-only q=2 cell is inference-ready with caveats. Read correlations with `corpairs(level = "phylogenetic")`. |
| Bivariate phylogenetic location-scale blocks | matching labelled `phylo()` terms across `mu1`, `mu2`, `sigma1`, and `sigma2` | Exact q=2+q=2 and broader all-four cells are fitted, but their evidence ranges from diagnostic-only to point-fit recovery; they are not blanket interval/coverage claims. |

</div>

When a predictor changes the phylogenetic SD itself, or when the model contains
location plus two scale submodels, continue with [When variance carries signal,
Part II: location-scale-scale models](location-scale-scale.html). For the
bivariate covariance layouts and the broader parity ladder, use the
[structural-dependence tutorial](https://itchyshin.github.io/drmTMB/articles/phylogenetic-spatial.html).

Anything outside these exact rows remains unpromoted, including other
non-Gaussian families, additional slope/label combinations, public matrix input
such as `phylo(1 | species, A = A)`, simultaneous `phylo()` and `spatial()`
layers, and structured random effects inside residual `rho12`.

## References

Hadfield, J. D. and Nakagawa, S. (2010). General quantitative genetic methods
for comparative biology: phylogenies, taxonomies and multi-trait models for
continuous and categorical characters. *Journal of Evolutionary Biology*,
23(3), 494-508.