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

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

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

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

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

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

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

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

## ----gaussian-confint, eval = has_ape-----------------------------------------
confint(fit, parm = "variance_components")[, c(
  "parm", "lower", "upper", "scale"
)]

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

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

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

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

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