---
title: "Animal models and additive relatedness"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Animal models and additive relatedness}
  %\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,
  dev.args = list(bg = "white")
)
library(drmTMB)

animal_guide_theme <- function() {
  ggplot2::theme_minimal(base_size = 11) +
    ggplot2::theme(
      panel.grid.minor = ggplot2::element_blank(),
      axis.title = ggplot2::element_text(colour = "grey15"),
      axis.text = ggplot2::element_text(colour = "grey25"),
      plot.title = ggplot2::element_text(
        face = "bold",
        colour = "grey10",
        margin = ggplot2::margin(b = 4)
      ),
      plot.subtitle = ggplot2::element_text(
        colour = "grey30",
        margin = ggplot2::margin(b = 8)
      ),
      plot.background = ggplot2::element_rect(fill = "white", colour = NA),
      panel.background = ggplot2::element_rect(fill = "white", colour = NA),
      legend.position = "bottom"
    )
}

animal_eye_theme <- function() {
  animal_guide_theme() +
    ggplot2::theme(
      panel.grid.major.y = ggplot2::element_blank(),
      legend.position = "none"
    )
}

relatedness_heatmap_data <- function(x) {
  out <- as.data.frame(as.table(x), stringsAsFactors = FALSE)
  names(out) <- c("row", "column", "relatedness")
  out$row <- factor(out$row, levels = rev(rownames(x)))
  out$column <- factor(out$column, levels = colnames(x))
  out
}

simulate_animal_guide_data <- function(seed = 20260523) {
  set.seed(seed)
  n_individual <- 32
  n_each <- 4
  individual_levels <- paste0("id_", seq_len(n_individual))
  individual <- factor(rep(individual_levels, each = n_each))
  age <- runif(length(individual), -1, 1)
  sex <- factor(rep(rep(c("female", "male"), each = n_each), n_individual / 2))

  A <- outer(seq_len(n_individual), seq_len(n_individual), function(i, j) {
    0.42^abs(i - j)
  })
  diag(A) <- diag(A) + 0.12
  dimnames(A) <- list(individual_levels, individual_levels)
  Ainv <- solve(A)

  animal_effect <- as.vector(t(chol(A)) %*% rnorm(n_individual)) * 0.38
  names(animal_effect) <- individual_levels
  sigma <- exp(-1.10 + 0.12 * (sex == "male"))
  body_size <- 1.6 + 0.45 * age + 0.22 * (sex == "male") +
    animal_effect[individual] + rnorm(length(individual), sd = sigma)

  list(
    data = data.frame(
      body_size = body_size,
      age = age,
      sex = sex,
      individual = individual
    ),
    A = A,
    Ainv = Ainv
  )
}

simulate_animal_q2_guide_data <- function(seed = 20260525) {
  set.seed(seed)
  n_individual <- 32
  n_each <- 4
  individual_levels <- paste0("id_", seq_len(n_individual))
  individual <- factor(rep(individual_levels, each = n_each))
  climate <- rnorm(length(individual))

  A <- outer(seq_len(n_individual), seq_len(n_individual), function(i, j) {
    0.40^abs(i - j)
  })
  diag(A) <- diag(A) + 0.10
  dimnames(A) <- list(individual_levels, individual_levels)
  Ainv <- solve(A)

  z_body <- rnorm(n_individual)
  z_activity <- 0.55 * z_body + sqrt(1 - 0.55^2) * rnorm(n_individual)
  u_body <- as.vector(t(chol(A)) %*% z_body) * 0.40
  u_activity <- as.vector(t(chol(A)) %*% z_activity) * 0.34
  names(u_body) <- individual_levels
  names(u_activity) <- individual_levels

  e_body <- rnorm(length(individual))
  e_activity <- 0.15 * e_body + sqrt(1 - 0.15^2) * rnorm(length(individual))

  list(
    data = data.frame(
      individual = individual,
      climate = climate,
      body_size = 1.2 + 0.30 * climate + u_body[individual] + 0.38 * e_body,
      activity = -0.3 - 0.24 * climate +
        u_activity[individual] + 0.34 * e_activity
    ),
    Ainv = Ainv
  )
}
```

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

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

Use `animal()` when individuals have known additive relatedness and the
scientific question is about among-individual location deviations after the
fixed effects have been included. The fitted first slices are Gaussian
location models. They are useful for animal-model and quantitative-genetic
workflows, but they are not yet full parity with every `phylo()` route. This
page assumes a pedigree or additive-relatedness matrix is already the right
description of your data; the [structural-dependence
overview](structural-dependence.html) is where that choice gets made.

## What is fitted today

| Question | Syntax | Status |
| --- | --- | --- |
| Does one Gaussian response have additive individual deviations? | `animal(1 | individual, pedigree = pedigree)`, `animal(1 | individual, A = A)`, or `animal(1 | individual, Ainv = Ainv)` in `mu` | Fitted first slice. Dense pedigree construction, covariance input, and precision input are supported. |
| Does one numeric fixed-effect slope also have additive individual deviations? | `animal(1 + x | individual, pedigree = pedigree)` in `mu` | Fitted first univariate one-slope slice with independent intercept and slope fields. Additional structured-slope layouts outside the exact fitted bivariate ledger cells and intercept-slope correlations remain planned. |
| Does a one-response Gaussian model have additive residual-scale deviations? | `animal(1 | individual, Ainv = Ainv)` in `sigma`, alone or matching the `mu` term; `animal(1 + x | individual, A = A)` for the exact one-slope route | The residual-scale intercept is fitted. The exact A-matrix q1 `sigma` one-slope route is fitted and inference-ready with caveats under raw uncorrected log-SD Wald-z intervals; profile is diagnostic-only at `g = 8`. Matching univariate `mu` and `sigma` animal intercepts estimate one animal mean-scale correlation. Pedigree/Ainv bridge marshalling, additional multiple or labelled sigma-slope layouts outside the exact fitted bivariate ledger cells, and structured slope correlations remain planned. |
| Do two Gaussian response means share an animal-model correlation? | matching labelled `animal(1 | p | individual, pedigree = pedigree)`, `animal(1 | p | individual, A = A)`, or `animal(1 | p | individual, Ainv = Ainv)` terms in `mu1` and `mu2` | Fitted first q=2 location-location slice. `corpairs(level = "animal")` reports the latent animal correlation separately from residual `rho12`. |
| Do animal-model deviations link means and residual scales across two responses? | the same labelled `animal(1 | p | individual, Ainv = Ainv)` term in `mu1`, `mu2`, `sigma1`, and `sigma2` | Fitted first q=4 location-scale slice. `corpairs(level = "animal")` reports six constant latent animal correlations; profile intervals for those derived correlations remain unavailable. |
| Does the animal effect change multiple slopes or a predictor-dependent correlation? | examples such as `animal(1 + x + z | individual, pedigree = pedigree)` or animal `corpair()` formulas | Planned. These routes need separate diagnostics, recovery tests, interval targets, and examples before use. |

The `pedigree` route expects a data frame whose identifiers match the grouping
variable used in the formula. The `A` route uses a precomputed additive
relationship covariance matrix, and the `Ainv` route uses the corresponding
precision matrix. Use exactly one of `pedigree`, `A`, or `Ainv` for a fitted
animal term.

## Start with the smallest useful model

For one response, start with the fitted location-intercept route:

```r
fit_animal <- drmTMB(
  y ~ treatment + animal(1 | individual, pedigree = pedigree),
  data = dat,
  family = gaussian()
)
```

For the exact Arc 1a REML route, keep `sigma ~ 1`, use an unlabelled intercept
or independent intercept-plus-one-numeric-slope shape, and set `REML = TRUE`:

```r
fit_animal_reml <- drmTMB(
  bf(
    y ~ x + animal(1 + x | individual, A = A),
    sigma ~ 1
  ),
  data = dat,
  family = gaussian(),
  REML = TRUE
)
```

The multi-seed campaign used the `A` representation shown here, with
`n_each = 20` and one fixed `M = 8` relationship matrix. Pedigree and `Ainv`
inputs have deterministic representation-parity evidence only, not multi-seed
campaign coverage. These domains do not admit slope-only, labelled,
multiple-slope, scale-side, bivariate, or non-Gaussian REML routes.

For two responses, use the same animal block label in `mu1` and `mu2`:

```r
fit_animal_q2 <- drmTMB(
  bf(
    mu1 = trait1 ~ treatment +
      animal(1 | p | individual, Ainv = Ainv),
    mu2 = trait2 ~ treatment +
      animal(1 | p | individual, Ainv = Ainv)
  ),
  data = dat,
  family = biv_gaussian()
)
corpairs(fit_animal_q2, level = "animal")
rho12(fit_animal_q2)
```

The two extractors answer different questions. `corpairs()` reports the fitted
latent animal-model correlation among individual location deviations. `rho12()`
reports the residual correlation between the two responses after the fixed
effects and random effects have been included.

For reports, keep this intercept-only q=2 animal correlation point-only. The
profile mechanism computes, but interval calibration and coverage are still
planned for this row:

```r
animal_pairs <- corpairs(fit_animal_q2, level = "animal")
plot_corpairs(animal_pairs)
```

That plot should be read as an animal-layer correlation display, not as a
residual `rho12` plot. Do not add an interval eye until this exact row has
validated interval evidence. The q=4 rows below are also point-only: current q=4
animal correlations are derived rows.

When the biological question is whether additive individual deviations in means
and residual scales covary, use the same labelled animal term in all four
bivariate endpoint formulas:

```r
fit_animal_q4 <- drmTMB(
  bf(
    mu1 = trait1 ~ treatment +
      animal(1 | p | individual, Ainv = Ainv),
    mu2 = trait2 ~ treatment +
      animal(1 | p | individual, Ainv = Ainv),
    sigma1 = ~ treatment +
      animal(1 | p | individual, Ainv = Ainv),
    sigma2 = ~ treatment +
      animal(1 | p | individual, Ainv = Ainv)
  ),
  data = dat,
  family = biv_gaussian()
)
corpairs(fit_animal_q4, level = "animal")
```

The q=4 path is constant across individuals and responses: it estimates four
latent animal endpoint scale parameters and six latent animal correlations.
Each endpoint's individual marginal SD also includes `sqrt(A[i, i])`. This is
not a direct-SD model and it does not make the residual correlation `rho12`
animal-structured.

## What to inspect

After fitting, inspect the animal layer before interpreting it:

| Output | Use |
| --- | --- |
| `check_drm(fit)` | Confirm the structured layer was recognized and review convergence diagnostics. |
| `summary(fit)$parameters` | Read the fitted latent animal scale `s`; individual `i` has marginal SD `s sqrt(A[i, i])`. |
| `ranef(fit, "animal_mu")` | Inspect conditional individual deviations on the fitted animal layer. |
| `summary(fit)$covariance` | Check how the animal SD or correlation is reported beside other covariance layers. |
| `profile_targets(fit)` | See which animal SD or q=2 correlation targets can be profiled directly; q=4 correlations are derived-only in the first slice. |
| `corpairs(fit, level = "animal")` | Read the fitted q=2 animal mean-mean correlation, or the six q=4 animal endpoint correlations when the model has the all-four block. |

If the ordinary repeatability model answers a useful sensitivity question, fit
`(1 | individual)` as a comparison. Report it as an ordinary individual
random-effect model, not as an additive animal model.

## Rendered checks

The small examples below keep three evidence grains separate. The relationship
matrix is raw structure supplied to the model. The SD plot converts the fitted
latent scale `s` to individual marginal SDs `s sqrt(A[i, i])` before comparing
them with residual `sigma`; it omits uncertainty intervals because the animal
q1 location-SD row is not interval-validated. The q=2 animal-correlation plot
is also a fitted point estimate from `corpairs()` with a dotted zero line; its
interval calibration remains planned.

```{r animal-guide-fit}
animal_example <- simulate_animal_guide_data()
animal_dat <- animal_example$data
A <- animal_example$A
Ainv <- animal_example$Ainv

fit_animal <- drmTMB(
  bf(
    body_size ~ age + sex + animal(1 | individual, Ainv = Ainv),
    sigma ~ 1
  ),
  family = gaussian(),
  data = animal_dat
)
```

```{r animal-relatedness-matrix, fig.width = 5.3, fig.height = 4.8, fig.cap = "Additive relatedness matrix used by the animal-model example; this heatmap is the known input structure, not an uncertainty display.", fig.alt = "Heatmap of the additive relatedness matrix for the animal-model example. Values are highest on the diagonal and fade as individuals are farther apart in the simulated pedigree-like ordering."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  animal_matrix <- relatedness_heatmap_data(A)

  ggplot2::ggplot(
    animal_matrix,
    ggplot2::aes(column, row, fill = relatedness)
  ) +
    ggplot2::geom_tile() +
    ggplot2::coord_equal() +
    ggplot2::scale_fill_gradientn(
      colours = c("#F7FBFF", "#C6DBEF", "#6BAED6", "#08519C"),
      name = "Additive\nrelatedness"
    ) +
    animal_guide_theme() +
    ggplot2::theme(
      axis.text = ggplot2::element_blank(),
      axis.ticks = ggplot2::element_blank(),
      panel.grid = ggplot2::element_blank()
    ) +
    ggplot2::labs(
      title = "Animal-model input structure",
      subtitle = "Known additive relatedness among individuals",
      x = "Individual",
      y = "Individual"
    )
}
```

```{r animal-sd-figure, fig.width = 6.6, fig.height = 3.4, fig.cap = "Fitted residual `sigma` and individual marginal animal SD point estimates `s sqrt(A[i,i])` from a univariate Gaussian model. The animal q1 location-SD interval is not validated, so this comparison is point-only.", fig.alt = "Horizontal point display comparing residual sigma with individual marginal animal standard deviations calculated as the fitted latent scale times the square root of each additive-relationship diagonal, without interval bars."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  animal_parameters <- summary(fit_animal)$parameters
  animal_scale <- animal_parameters[
    animal_parameters$parm == "sd:mu:animal(1 | individual)", "estimate"
  ]
  animal_node_sd <- animal_scale * sqrt(diag(A))
  animal_sd <- data.frame(
    estimate = c(
      unname(sigma(fit_animal)[1]),
      animal_node_sd
    ),
    label = c(
      "Residual\nsigma",
      rep("Animal node\nmarginal SD", length(animal_node_sd))
    )
  )

  ggplot2::ggplot(animal_sd, ggplot2::aes(y = label)) +
    ggplot2::geom_vline(
      xintercept = 0,
      linewidth = 0.45,
      linetype = "dashed",
      colour = "grey55"
    ) +
    ggplot2::geom_point(
      ggplot2::aes(x = estimate),
      shape = 21,
      size = 3.5,
      stroke = 1,
      fill = "white",
      colour = "#0072B2"
    ) +
    ggplot2::scale_x_continuous(
      expand = ggplot2::expansion(mult = c(0.02, 0.05))
    ) +
    animal_eye_theme() +
    ggplot2::labs(
      title = "Animal marginal SD is separate from residual sigma",
      subtitle = "Node SD = s sqrt(Aii); point estimates only",
      x = "Fitted standard deviation",
      y = NULL
    )
}
```

For two response means, the animal-model correlation row is a fitted latent
correlation among individual additive deviations. It is not the residual
`rho12` correlation and it is not a raw sample correlation among observations.

```{r animal-q2-correlation-fit}
animal_q2_example <- simulate_animal_q2_guide_data()
animal_q2_dat <- animal_q2_example$data
Ainv <- animal_q2_example$Ainv

fit_animal_q2_example <- drmTMB(
  bf(
    mu1 = body_size ~ climate +
      animal(1 | p | individual, Ainv = Ainv),
    mu2 = activity ~ climate +
      animal(1 | p | individual, Ainv = Ainv),
    sigma1 = ~ 1,
    sigma2 = ~ 1,
    rho12 = ~ 1
  ),
  family = c(gaussian(), gaussian()),
  data = animal_q2_dat
)

animal_q2_pairs <- corpairs(
  fit_animal_q2_example,
  level = "animal"
)
animal_q2_pairs
```

```{r animal-q2-confidence-eye, fig.width = 6.4, fig.height = 2.7, fig.cap = "Animal-model intercept-only q=2 location-location point estimate from `corpairs()`; the dotted vertical line marks zero correlation and no interval is shown because calibration remains planned.", fig.alt = "Single-row point plot for the animal-model mean-mean correlation, with a hollow point estimate to the right of the dotted zero reference line."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  animal_q2_display <- animal_q2_pairs
  animal_q2_display$display_label <- "Animal\nmu1-mu2"

  plot_corpairs(
    animal_q2_display,
    colour = "level",
    label = "display_label",
    facet = NULL
  ) +
    ggplot2::scale_colour_manual(values = c("animal" = "#0072B2")) +
    ggplot2::scale_fill_manual(values = c("animal" = "#0072B2")) +
    animal_eye_theme() +
    ggplot2::labs(
      title = "Animal-model latent mean correlation",
      subtitle = "Point estimate only; dotted line marks zero",
      x = "Correlation estimate"
    )
}
```

## Boundaries

The following animal-model routes remain planned:

- sparse large-pedigree precision construction;
- multiple animal-model slopes such as
  `animal(1 + x + z | individual, pedigree = pedigree)`;
- animal-model intercept-slope correlations;
- multiple or labelled residual-scale structured slopes beyond the exact
  A-matrix q1 `animal(1 + x | individual, A = A)` sigma route;
- predictor-dependent animal `corpair()` regressions;
- generic direct-SD grammar for animal-model standard deviations;
- richer non-Gaussian animal structured effects: labelled covariance blocks,
  multiple structured slopes, structured scale routes beyond the NB2 one-slope
  gate, and q=2 or q=4 non-Gaussian animal covariance blocks. Exact q=1
  animal gates already fit today for `poisson()`, `nbinom2()`, and `beta()`:
  Poisson, NB2, and beta `mu` accept an unlabelled intercept plus one independent
  slope; NB2 `sigma` accepts the same q1 intercept-plus-one-slope form; beta
  `sigma` accepts an intercept only. These rows retain their recorded
  point-fit/recovery tiers -- trust the point estimate, not the interval.

Use the [structural-dependence overview](structural-dependence.html) when you
are choosing among `animal()`, `phylo()`, `spatial()`, and `relmat()`. Use the
detailed [structural-dependence tutorial](https://itchyshin.github.io/drmTMB/articles/phylogenetic-spatial.html) when you
need the current fitted examples, equations, and broader parity ladder.