## ----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
  )
}

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

## ----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"
    )
}

## ----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
    )
}

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

## ----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"
    )
}