## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")

## -----------------------------------------------------------------------------
library(drmTMB)

set.seed(101)
n <- 180
seedlings <- data.frame(
  drought = factor(rep(c("ambient", "dry"), each = n / 2))
)
dry <- as.numeric(seedlings$drought == "dry")
mu <- 1.2 - 0.45 * dry
sigma <- exp(-1 + 0.35 * dry)
seedlings$growth <- mu + sigma * rt(n, df = 5)

## -----------------------------------------------------------------------------
fit_gaussian <- drmTMB(
  bf(growth ~ drought, sigma ~ drought),
  family = gaussian(),
  data = seedlings
)

## -----------------------------------------------------------------------------
fit_student <- drmTMB(
  bf(growth ~ drought, sigma ~ drought, nu ~ 1),
  family = student(),
  data = seedlings
)

## -----------------------------------------------------------------------------
student_checks <- check_drm(fit_student)
student_checks
student_checks[
  student_checks$check == "student_nu",
  c("status", "value", "message")
]

## -----------------------------------------------------------------------------
coef(fit_student, "mu")
coef(fit_student, "sigma")
coef(fit_student, "nu")

## -----------------------------------------------------------------------------
head(predict(fit_student, dpar = "nu"))

## -----------------------------------------------------------------------------
AIC(fit_gaussian, fit_student)

## -----------------------------------------------------------------------------
coef(fit_gaussian, "mu")
coef(fit_student, "mu")

coef(fit_gaussian, "sigma")
coef(fit_student, "sigma")

## ----robust-student-tail-figure, eval=requireNamespace("ggplot2", quietly = TRUE), fig.width=7.2, fig.height=4.4, fig.cap="Robust-model check for the seedling example. Faint points show observed growth values; overlaid points compare fitted Gaussian and Student-t expected growth by drought treatment. No interval bars are drawn because this fixture is a raw-data and fitted-point comparison, not an interval summary.", fig.alt="Jittered point plot of seedling growth by ambient and dry drought treatments. Faint raw observations show heavy-tailed residual variation, and overlaid points compare Gaussian and Student-t fitted expected growth for each treatment."----
library(ggplot2)

student_plot_grid <- data.frame(
  drought = factor(c("ambient", "dry"), levels = levels(seedlings$drought))
)
student_plot_means <- rbind(
  data.frame(
    student_plot_grid,
    model = "Gaussian",
    fitted_mu = predict(fit_gaussian, newdata = student_plot_grid, dpar = "mu")
  ),
  data.frame(
    student_plot_grid,
    model = "Student-t",
    fitted_mu = predict(fit_student, newdata = student_plot_grid, dpar = "mu")
  )
)

ggplot(seedlings, aes(drought, growth, colour = drought)) +
  geom_jitter(width = 0.12, height = 0, alpha = 0.22, size = 1.1) +
  geom_point(
    data = student_plot_means,
    aes(y = fitted_mu, shape = model),
    position = position_dodge(width = 0.35),
    size = 3.2,
    stroke = 1.1
  ) +
  scale_colour_manual(values = c("ambient" = "#0072B2", "dry" = "#D55E00")) +
  scale_shape_manual(values = c("Gaussian" = 16, "Student-t" = 1)) +
  labs(
    title = "Robust fits should be checked against the raw tails",
    subtitle = "Points are observed growth values; symbols are fitted expected growth",
    x = "Drought treatment",
    y = "Growth",
    colour = "Treatment",
    shape = "Model"
  ) +
  guides(colour = "none") +
  theme_minimal(base_size = 11) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(colour = "grey30")
  )