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

## -----------------------------------------------------------------------------
set.seed(196)
n_onset <- 240
onset_data <- data.frame(
  canopy = factor(
    rep(c("open", "closed"), each = n_onset / 2),
    levels = c("open", "closed")
  ),
  ndvi = as.numeric(scale(runif(n_onset, 0.15, 0.85)))
)

closed <- as.numeric(onset_data$canopy == "closed")
mu_onset <- plogis(-0.85 - 0.35 * closed + 1.10 * onset_data$ndvi)
onset_data$early_onset <- rbinom(n_onset, size = 1, prob = mu_onset)

fit_onset <- drmTMB(
  bf(early_onset ~ canopy + ndvi),
  family = stats::binomial(link = "logit"),
  data = onset_data
)

coef(fit_onset, "mu")

## -----------------------------------------------------------------------------
new_onset <- data.frame(
  canopy = factor(c("open", "closed"), levels = levels(onset_data$canopy)),
  ndvi = c(0, 0)
)

data.frame(
  canopy = new_onset$canopy,
  ndvi = new_onset$ndvi,
  early_onset_probability = predict(fit_onset, newdata = new_onset, dpar = "mu")
)

## -----------------------------------------------------------------------------
set.seed(197)
n <- 360
seed_trials <- data.frame(
  treatment = factor(
    rep(c("open", "sheltered"), each = n / 2),
    levels = c("open", "sheltered")
  ),
  moisture = as.numeric(scale(runif(n, 0.1, 0.9))),
  trials = sample(18:32, n, replace = TRUE)
)

sheltered <- as.numeric(seed_trials$treatment == "sheltered")
mu_seed <- plogis(-0.55 + 0.70 * sheltered + 0.45 * seed_trials$moisture)
sigma_seed <- exp(-1.15 - 0.35 * sheltered)
phi_seed <- 1 / sigma_seed^2
tray_probability <- rbeta(
  n,
  shape1 = mu_seed * phi_seed,
  shape2 = (1 - mu_seed) * phi_seed
)

seed_trials$germinated <- rbinom(n, size = seed_trials$trials, prob = tray_probability)
seed_trials$failed <- seed_trials$trials - seed_trials$germinated

head(seed_trials)

## -----------------------------------------------------------------------------
fit_seed <- drmTMB(
  bf(cbind(germinated, failed) ~ treatment + moisture, sigma ~ treatment),
  family = beta_binomial(),
  data = seed_trials
)

## -----------------------------------------------------------------------------
check_drm(fit_seed)

## -----------------------------------------------------------------------------
coef(fit_seed, "mu")
coef(fit_seed, "sigma")

sigma_ratio <- exp(coef(fit_seed, "sigma")["treatmentsheltered"])
c(
  sigma_ratio_sheltered_vs_open = sigma_ratio,
  precision_ratio_sheltered_vs_open = sigma_ratio^(-2)
)

## -----------------------------------------------------------------------------
new_seed_trays <- data.frame(
  treatment = factor(c("open", "sheltered"), levels = levels(seed_trials$treatment)),
  moisture = c(0, 0),
  trials = c(24, 24)
)

mu_hat <- predict(fit_seed, newdata = new_seed_trays, dpar = "mu")
sigma_hat <- predict(fit_seed, newdata = new_seed_trays, dpar = "sigma")
prop_var <- mu_hat * (1 - mu_hat) *
  (1 + new_seed_trays$trials * sigma_hat^2) /
  (new_seed_trays$trials * (1 + sigma_hat^2))

data.frame(
  treatment = new_seed_trays$treatment,
  trials = new_seed_trays$trials,
  expected_probability = mu_hat,
  expected_successes = new_seed_trays$trials * mu_hat,
  sigma = sigma_hat,
  phi = 1 / sigma_hat^2,
  proportion_sd = sqrt(prop_var)
)

## ----beta-binomial-tray-figure, eval=requireNamespace("ggplot2", quietly = TRUE), fig.width=7.2, fig.height=4.4, fig.cap="Beta-binomial tray summary for the seed-germination example. Faint points are observed tray proportions; overlaid points are fitted expected germination probabilities; vertical bars show plus or minus one fitted proportion-level standard deviation, not confidence intervals.", fig.alt="Jittered point plot of observed germination proportions for open and sheltered trays. Overlaid larger points show fitted expected probabilities, and vertical bars show one fitted proportion-level standard deviation for each treatment."----
library(ggplot2)

seed_trials$observed_proportion <- seed_trials$germinated / seed_trials$trials
seed_plot_summary <- data.frame(
  treatment = new_seed_trays$treatment,
  expected_probability = mu_hat,
  proportion_sd = sqrt(prop_var)
)
seed_plot_summary$lower <- pmax(
  0,
  seed_plot_summary$expected_probability - seed_plot_summary$proportion_sd
)
seed_plot_summary$upper <- pmin(
  1,
  seed_plot_summary$expected_probability + seed_plot_summary$proportion_sd
)

ggplot(seed_trials, aes(treatment, observed_proportion, colour = treatment)) +
  geom_jitter(width = 0.12, height = 0, alpha = 0.18, size = 0.9) +
  geom_errorbar(
    data = seed_plot_summary,
    aes(
      y = expected_probability,
      ymin = lower,
      ymax = upper
    ),
    width = 0.12,
    linewidth = 0.8
  ) +
  geom_point(data = seed_plot_summary, aes(y = expected_probability), size = 3) +
  scale_colour_manual(values = c("open" = "#0072B2", "sheltered" = "#009E73")) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(
    title = "Denominator-aware proportions can still show raw trays",
    subtitle = "Bars show fitted tray-level scatter, not confidence intervals",
    x = "Microsite treatment",
    y = "Germinated proportion",
    colour = "Treatment"
  ) +
  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")
  )

## -----------------------------------------------------------------------------
set.seed(198)
n_cover <- 300
cover_data <- data.frame(
  grazing = factor(
    rep(c("ungrazed", "grazed"), each = n_cover / 2),
    levels = c("ungrazed", "grazed")
  ),
  moisture = as.numeric(scale(runif(n_cover, 0.05, 0.95)))
)

grazed <- as.numeric(cover_data$grazing == "grazed")
mu_cover <- plogis(0.35 - 0.75 * grazed + 0.35 * cover_data$moisture)
sigma_cover <- exp(-1.10 + 0.45 * grazed)
phi_cover <- 1 / sigma_cover^2
cover_data$cover <- rbeta(
  n_cover,
  shape1 = mu_cover * phi_cover,
  shape2 = (1 - mu_cover) * phi_cover
)

fit_cover <- drmTMB(
  bf(cover ~ grazing + moisture, sigma ~ grazing),
  family = beta(),
  data = cover_data
)

check_drm(fit_cover)

## -----------------------------------------------------------------------------
coef(fit_cover, "mu")
coef(fit_cover, "sigma")

new_cover <- data.frame(
  grazing = factor(c("ungrazed", "grazed"), levels = levels(cover_data$grazing)),
  moisture = c(0, 0)
)

mu_cover_hat <- predict(fit_cover, newdata = new_cover, dpar = "mu")
sigma_cover_hat <- predict(fit_cover, newdata = new_cover, dpar = "sigma")

data.frame(
  grazing = new_cover$grazing,
  expected_cover = mu_cover_hat,
  sigma = sigma_cover_hat,
  phi = 1 / sigma_cover_hat^2,
  cover_sd = sqrt(
    mu_cover_hat * (1 - mu_cover_hat) *
      sigma_cover_hat^2 / (1 + sigma_cover_hat^2)
  )
)