## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.1,
  dpi = 144
)
if (!"package:drmTMB" %in% search()) {
  library(drmTMB)
}

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

location_scale_habitat_scales <- function() {
  list(
    ggplot2::scale_colour_manual(
      values = c("forest" = "#0072B2", "grassland" = "#009E73")
    ),
    ggplot2::scale_fill_manual(
      values = c("forest" = "#0072B2", "grassland" = "#009E73")
    )
  )
}

## ----location-scale-simulate--------------------------------------------------
set.seed(42)
n <- 240
dat <- data.frame(
  habitat = factor(rep(c("forest", "grassland"), each = n / 2)),
  temperature = rnorm(n)
)

mu <- 8 + 1.2 * (dat$habitat == "grassland") + 0.6 * dat$temperature
sigma_true <- exp(log(0.6) + 0.55 * (dat$habitat == "grassland"))
dat$growth <- rnorm(n, mean = mu, sd = sigma_true)

## ----location-scale-fit-------------------------------------------------------
fit_growth <- drmTMB(
  drm_formula(growth ~ habitat + temperature, sigma ~ habitat),
  family = gaussian(),
  data = dat
)

## ----location-scale-check-----------------------------------------------------
check_drm(fit_growth)

## ----location-scale-profile-targets-------------------------------------------
profile_targets(fit_growth)[
  ,
  c("parm", "estimate", "profile_ready", "profile_note")
]

## ----location-scale-summary---------------------------------------------------
summary(fit_growth)

## ----location-scale-interpret-sigma-------------------------------------------
sigma_habitat <- coef(fit_growth, "sigma")["habitatgrassland"]
data.frame(
  coefficient = sigma_habitat,
  residual_sd_ratio = exp(sigma_habitat),
  residual_variance_ratio = exp(2 * sigma_habitat)
)

## ----location-scale-response-table--------------------------------------------
newdat <- data.frame(
  habitat = factor(c("forest", "grassland"), levels = levels(dat$habitat)),
  temperature = 0
)

growth_report <- data.frame(
  habitat = newdat$habitat,
  fitted_mean_growth = predict(fit_growth, newdata = newdat, dpar = "mu"),
  fitted_residual_sd = predict(fit_growth, newdata = newdat, dpar = "sigma")
)
growth_report$fitted_residual_variance <- growth_report$fitted_residual_sd^2
growth_report

## ----location-scale-growth-figure, fig.cap = "Raw growth observations and fitted response-scale `mu` surfaces for the Gaussian location-scale example; ribbons are 95% Wald confidence bands from `predict_parameters()`.", fig.alt = "Scatterplot of simulated growth over temperature for forest and grassland habitats, with fitted mean growth lines and 95 percent Wald confidence bands for each habitat."----
growth_mu_grid <- prediction_grid(
  fit_growth,
  focal = c("temperature", "habitat"),
  at = list(
    temperature = seq(
      min(dat$temperature),
      max(dat$temperature),
      length.out = 80
    )
  )
)

growth_mu_surface <- predict_parameters(
  fit_growth,
  newdata = growth_mu_grid,
  dpar = "mu",
  conf.int = TRUE
)

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

if (requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::ggplot(
    dat,
    ggplot2::aes(x = temperature, y = growth, colour = habitat)
  ) +
    ggplot2::geom_point(alpha = 0.38, size = 1.25) +
    ggplot2::geom_ribbon(
      data = growth_mu_surface,
      ggplot2::aes(
        x = temperature,
        ymin = conf.low,
        ymax = conf.high,
        fill = habitat
      ),
      inherit.aes = FALSE,
      alpha = 0.18,
      colour = NA
    ) +
    ggplot2::geom_line(
      data = growth_mu_surface,
      ggplot2::aes(y = estimate),
      linewidth = 0.85
    ) +
    location_scale_habitat_scales() +
    ggplot2::labs(
      title = "Mean growth and observed scatter",
      subtitle = "Points are raw growth; ribbons are 95% Wald bands for mu",
      x = "Temperature",
      y = "Growth",
      colour = "Habitat",
      fill = "Habitat"
    ) +
    location_scale_theme() +
    ggplot2::guides(fill = "none")
}

## ----location-scale-sigma-contrast-figure, fig.width = 5.8, fig.height = 2.6, fig.cap = "Fitted residual standard deviations by habitat at average temperature; horizontal intervals are 95% Wald confidence intervals requested from `predict_parameters()`.", fig.alt = "Horizontal point interval plot comparing fitted residual standard deviation for forest and grassland habitats at temperature zero, with 95 percent Wald confidence intervals."----
growth_sigma_grid <- prediction_grid(
  fit_growth,
  focal = "habitat",
  condition = list(temperature = 0)
)

growth_sigma_surface <- predict_parameters(
  fit_growth,
  newdata = growth_sigma_grid,
  dpar = "sigma",
  conf.int = TRUE
)

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

if (requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::ggplot(
    growth_sigma_surface,
    ggplot2::aes(x = estimate, y = habitat, colour = habitat)
  ) +
    ggplot2::geom_segment(
      ggplot2::aes(x = conf.low, xend = conf.high, yend = habitat),
      linewidth = 0.8
    ) +
    ggplot2::geom_point(size = 2.8) +
    location_scale_habitat_scales() +
    ggplot2::labs(
      title = "Residual SD differs by habitat",
      subtitle = "Points are fitted sigma values; bars are 95% Wald intervals",
      x = "Fitted residual SD (sigma)",
      y = NULL,
      colour = "Habitat"
    ) +
    location_scale_theme() +
    ggplot2::guides(colour = "none")
}

## ----location-scale-translation-table-----------------------------------------
growth_translation <- data.frame(
  model_piece = c(
    "fixed mean slope",
    "fixed residual-SD contrast",
    "fixed residual-variance contrast"
  ),
  fitted_term = c(
    "mu:temperature",
    "exp(sigma:habitatgrassland)",
    "exp(2 * sigma:habitatgrassland)"
  ),
  response_scale_value = c(
    unname(coef(fit_growth, "mu")["temperature"]),
    unname(exp(sigma_habitat)),
    unname(exp(2 * sigma_habitat))
  ),
  interpretation = c(
    "mean growth change per one-unit temperature increase",
    "grassland residual SD divided by forest residual SD",
    "grassland residual variance divided by forest residual variance"
  )
)
growth_translation$response_scale_value <-
  round(growth_translation$response_scale_value, 3)
growth_translation