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

which_scale_theme <- function() {
  ggplot2::theme_minimal(base_size = 11) +
    ggplot2::theme(
      panel.grid.minor = ggplot2::element_blank(),
      panel.grid.major.x = ggplot2::element_line(
        colour = "grey90",
        linewidth = 0.3
      ),
      panel.grid.major.y = ggplot2::element_line(
        colour = "grey87",
        linewidth = 0.35
      ),
      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)
      ),
      legend.position = "bottom",
      legend.title = ggplot2::element_text(colour = "grey20")
    )
}

## ----scale-audit-data---------------------------------------------------------
set.seed(42)
n_population <- 32
n_each <- 6

population_info <- data.frame(
  population = factor(seq_len(n_population)),
  habitat = rep(c("forest", "grassland"), length.out = n_population)
)

fish <- population_info[rep(seq_len(n_population), each = n_each), ]
fish$temperature <- rnorm(nrow(fish))
fish$reliability <- ifelse(seq_len(nrow(fish)) %% 3 == 0, 2, 1)

pop_sd <- exp(-0.8 + 0.7 * (population_info$habitat == "grassland"))
b_population <- rnorm(n_population, sd = pop_sd)
fish$sigma_true <- exp(-0.7 + 0.3 * fish$temperature)
fish$growth <- 1.2 + 0.55 * fish$temperature +
  b_population[fish$population] +
  rnorm(nrow(fish), sd = fish$sigma_true)

## ----residual-scale-fit-------------------------------------------------------
fit_sigma <- drmTMB(
  bf(growth ~ temperature, sigma ~ temperature),
  family = gaussian(),
  data = fish
)

summary(fit_sigma)
round(coef(fit_sigma, "sigma"), 3)
round(range(sigma(fit_sigma)), 3)

## ----which-scale-residual-sigma-figure, fig.cap = "Fitted residual standard deviation over temperature for the `sigma ~ temperature` example; the ribbon is a 95% Wald confidence band from `predict_parameters()`.", fig.alt = "Line plot of fitted residual standard deviation over temperature with a 95 percent Wald confidence band. No raw growth points are shown on the sigma axis."----
sigma_temperature_grid <- prediction_grid(
  fit_sigma,
  focal = "temperature",
  at = list(
    temperature = seq(
      min(fish$temperature),
      max(fish$temperature),
      length.out = 80
    )
  )
)

sigma_temperature_surface <- predict_parameters(
  fit_sigma,
  newdata = sigma_temperature_grid,
  dpar = "sigma",
  conf.int = TRUE
)

unique(sigma_temperature_surface[, c(
  "dpar",
  "conf.status",
  "interval_source",
  "conf.level"
)])

if (requireNamespace("ggplot2", quietly = TRUE)) {
  plot_parameter_surface(
    sigma_temperature_surface,
    x = "temperature",
    dpar = "sigma",
    facet = NULL,
    point = FALSE
  ) +
    ggplot2::labs(
      title = "Residual scale is observation-level",
      subtitle = "Ribbon is a 95% Wald band for fitted sigma",
      x = "Temperature",
      y = "Fitted residual SD (sigma)"
    ) +
    which_scale_theme()
}

## ----likelihood-weights-fit---------------------------------------------------
fit_weighted <- drmTMB(
  bf(growth ~ temperature, sigma ~ 1),
  family = gaussian(),
  data = fish,
  weights = reliability
)

summary(fit_weighted)
head(weights(fit_weighted), 8)

## ----meta-known-v-fit---------------------------------------------------------
set.seed(101)
n_effect <- 50
meta <- data.frame(treatment = rep(c(0, 1), each = n_effect / 2))
meta$vi <- runif(n_effect, 0.02, 0.08)

mu_meta <- 0.1 + 0.25 * meta$treatment
sigma_meta <- exp(-1.1 + 0.35 * meta$treatment)
meta$yi <- rnorm(n_effect, mu_meta, sqrt(meta$vi + sigma_meta^2))

fit_meta <- drmTMB(
  bf(yi ~ treatment + meta_V(V = vi), sigma ~ treatment),
  family = gaussian(),
  data = meta
)

summary(fit_meta)

meta_report <- data.frame(
  treatment = meta$treatment,
  known_sampling_variance = meta$vi,
  extra_heterogeneity_sd = sigma(fit_meta)
)
meta_report$extra_heterogeneity_variance <-
  meta_report$extra_heterogeneity_sd^2
meta_report$total_observation_variance <-
  meta_report$known_sampling_variance +
  meta_report$extra_heterogeneity_variance

meta_summary <- aggregate(
  meta_report[c(
    "extra_heterogeneity_sd",
    "extra_heterogeneity_variance",
    "total_observation_variance"
  )],
  by = list(treatment = meta_report$treatment),
  FUN = mean
)
round(meta_summary, 3)

## ----random-effect-scale-fit--------------------------------------------------
fit_sd <- drmTMB(
  bf(
    growth ~ temperature + (1 | population),
    sigma ~ temperature,
    sd(population) ~ habitat
  ),
  family = gaussian(),
  data = fish
)

round(coef(fit_sd, "sd(population)"), 3)
round(tapply(
  predict(fit_sd, dpar = "sd(population)"),
  population_info$habitat,
  mean
), 3)

## ----which-scale-population-sd-figure, fig.width = 5.8, fig.height = 2.6, fig.cap = "Fitted among-population standard deviations by habitat for the `sd(population) ~ habitat` example; points are fitted random-effect SDs, with no interval drawn because the prediction table marks this random-effect-SD surface as interval-unavailable.", fig.alt = "Point plot comparing fitted among-population standard deviation for forest and grassland habitats, with no interval bars because this random-effect standard-deviation surface has no supported interval in the current prediction table."----
sd_population_rows <- predict_parameters(
  fit_sd,
  dpar = "sd(population)",
  conf.int = TRUE
)
sd_population_rows$habitat <- population_info$habitat[sd_population_rows$row]

unique(sd_population_rows[, c(
  "dpar",
  "component",
  "conf.status",
  "interval_source"
)])

sd_population_display <- aggregate(
  estimate ~ habitat,
  data = sd_population_rows,
  FUN = mean
)

if (requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::ggplot(
    sd_population_display,
    ggplot2::aes(x = estimate, y = habitat)
  ) +
    ggplot2::geom_point(size = 3, colour = "#D55E00") +
    ggplot2::labs(
      title = "Group-level scale is not residual sigma",
      subtitle = "Points are fitted sd(population); no supported interval is drawn",
      x = "Fitted among-population SD",
      y = NULL
    ) +
    which_scale_theme()
}

## ----residual-coscale-fit-----------------------------------------------------
set.seed(12)
behaviour <- data.frame(treatment = rep(c(0, 1), each = 50))
Sigma0 <- matrix(c(0.6^2, 0.2 * 0.6 * 0.5,
                   0.2 * 0.6 * 0.5, 0.5^2), 2, 2)
Sigma1 <- matrix(c(0.6^2, 0.65 * 0.6 * 0.5,
                   0.65 * 0.6 * 0.5, 0.5^2), 2, 2)
Y <- matrix(NA_real_, nrow(behaviour), 2)

for (i in seq_len(nrow(behaviour))) {
  Sigma_i <- if (behaviour$treatment[i] == 0) Sigma0 else Sigma1
  mu_i <- c(1 + 0.2 * behaviour$treatment[i],
            0.5 + 0.1 * behaviour$treatment[i])
  Y[i, ] <- as.numeric(mu_i + t(chol(Sigma_i)) %*% rnorm(2))
}

behaviour$activity <- Y[, 1]
behaviour$boldness <- Y[, 2]

fit_rho12 <- drmTMB(
  bf(
    mu1 = activity ~ treatment,
    mu2 = boldness ~ treatment,
    sigma1 = ~ treatment,
    sigma2 = ~ treatment,
    rho12 = ~ treatment
  ),
  family = c(gaussian(), gaussian()),
  data = behaviour
)

summary(fit_rho12)
round(coef(fit_rho12, "rho12"), 3)
round(tapply(rho12(fit_rho12), behaviour$treatment, mean), 3)