########################################################################
# PLS 206 — Bonus: the GLM family
# Grey Monroe, Fall 2026
#
# Every block is tagged:
#   [CORE]   the analysis: code you would run on your own data
#   [HANDY]  useful, not essential: wrangling, plotting, small helpers
#   [DEMO]   for the lecture only: simulated data, sums by hand, and the
#            code behind slide figures
########################################################################

# [CORE]
# install.packages(c("ggplot2", "MASS"))
library(ggplot2)
library(MASS)

# ---- 1) The data: seven responses, one predictor ---------------------

# 120 wheat fields, each with one nitrogen rate (kg N/ha). Seven
# responses measured on every field (simulated):

# | column | what it is | kind of number | |---|---|---| | yield |
# grain yield, t/ha | continuous | | lodged | crop fell over: 1 = yes,
# 0 = no | 0/1 | | sown, emerged | seeds sown in a row, seedlings that
# emerged | $k$ of $n$ | | beetles | ground beetles in one pitfall
# trap, 7 days | count, same effort everywhere | | aphids | aphids on
# 10 tillers | clumped count | | leaves, lesions | leaves scored,
# septoria lesions on them | count, effort differs | | weed_biomass |
# g/m² | positive, right-skewed |

# [CORE]
survey <- read.csv("https://greymonroe.github.io/PLS_206/data/field_survey.csv")
head(survey)

# [DEMO: the seven responses against nitrogen]
long <- rbind(
  data.frame(response = "yield (t/ha)",         nitrogen = survey$nitrogen, y = survey$yield),
  data.frame(response = "lodged (0/1)",         nitrogen = survey$nitrogen, y = survey$lodged),
  data.frame(response = "emerged / sown",       nitrogen = survey$nitrogen, y = survey$emerged / survey$sown),
  data.frame(response = "beetles per trap",     nitrogen = survey$nitrogen, y = survey$beetles),
  data.frame(response = "aphids per 10 tillers", nitrogen = survey$nitrogen, y = survey$aphids),
  data.frame(response = "lesions per leaf",     nitrogen = survey$nitrogen, y = survey$lesions / survey$leaves),
  data.frame(response = "weed biomass (g/m2)",  nitrogen = survey$nitrogen, y = survey$weed_biomass))
long$response <- factor(long$response, unique(long$response))
ggplot(long, aes(nitrogen, y)) +
  geom_point(colour = "#3D7BEB", alpha = 0.6) +
  facet_wrap(~ response, scales = "free_y", nrow = 2) +
  theme_classic()

# ---- 2) The link: keep the mean where it can live --------------------

# A straight line runs from $-\infty$ to $+\infty$. The inverse link
# bends it into a range the mean can actually take. The link models the
# log (or logit) of the mean; it never transforms the data.

# [DEMO: identity, log, and logit, drawn]
eta <- seq(-4, 4, 0.05)
links <- rbind(data.frame(link = "identity: mu = eta",            eta, mu = eta),
               data.frame(link = "log: mu = exp(eta)",            eta, mu = exp(eta)),
               data.frame(link = "logit: mu = plogis(eta)",       eta, mu = plogis(eta)))
ggplot(links, aes(eta, mu)) +
  geom_line(colour = "#E0524A", linewidth = 1.2) +
  facet_wrap(~ link, scales = "free_y") +
  labs(x = "eta (the straight line)", y = "mu (the mean of y)") +
  theme_classic()

# family(model)$linkfun and family(model)$linkinv hold the two
# directions for any fitted model; plogis() and exp() are the inverse
# links you will use most.

# [HANDY]
binomial()$linkfun(0.8)   # logit: log(0.8 / 0.2)
binomial()$linkinv(1.39)  # back to a probability, same as plogis(1.39)
poisson()$linkinv(2)      # exp(2)

# ---- 3) The family: how spread grows with the mean -------------------

# The second choice. Each family says how the variance depends on the
# mean: Normal $\sigma^2$ (constant), Poisson $\mu$, negative binomial
# $\mu + \mu^2/\theta$, Gamma $\mu^2/k$ (constant CV), binomial
# $p(1-p)/n$ for a proportion.

# [DEMO: same rising mean, five families]
set.seed(5)
x <- runif(400, 0, 10)
mv <- rbind(
  data.frame(family = "Normal",            x, y = rnorm(400, 2 + 1.6 * x, 2.5)),
  data.frame(family = "Poisson",           x, y = rpois(400, exp(0.5 + 0.28 * x))),
  data.frame(family = "negative binomial", x, y = rnbinom(400, mu = exp(0.5 + 0.28 * x), size = 2)),
  data.frame(family = "Gamma",             x, y = rgamma(400, shape = 3, rate = 3 / exp(0.5 + 0.28 * x))),
  data.frame(family = "binomial (k of 20)", x, y = rbinom(400, 20, plogis(-4 + 0.8 * x)) / 20))
mv$family <- factor(mv$family, unique(mv$family))
ggplot(mv, aes(x, y)) +
  geom_point(colour = "#3D7BEB", alpha = 0.4, size = 1) +
  facet_wrap(~ family, nrow = 1, scales = "free_y") +
  theme_classic()

########################################################################
# The seven models
########################################################################

# Every section has the same shape: fit, read the coefficients, turn
# $\beta_1$ into a statement in the response's own terms. Nitrogen runs
# 0 to 200, so the effects are reported per 100 kg N.

# ---- 4) Normal: yield ------------------------------------------------

# $\text{yield} \sim \text{Normal}(\mu, \sigma^2)$, $\mu = \beta_0 +
# \beta_1\,\text{N}$. lm() is this GLM: glm(family = gaussian) gives
# the same estimates.

# [CORE]
m1 <- glm(yield ~ nitrogen, family = gaussian, data = survey)
coef(summary(m1))
coef(lm(yield ~ nitrogen, data = survey))   # identical
100 * coef(m1)["nitrogen"]                  # t/ha per 100 kg N

# ---- 5) Binomial, 0/1: did the crop lodge? ---------------------------

# $\text{lodged} \sim \text{Bernoulli}(p)$, $\log\frac{p}{1-p} =
# \beta_0 + \beta_1\,\text{N}$. $e^{\beta_1}$ is an odds ratio.

# [CORE]
m2 <- glm(lodged ~ nitrogen, family = binomial, data = survey)
coef(summary(m2))
exp(coef(m2)["nitrogen"] * 100)                    # odds ratio, per 100 kg N
exp(confint(m2)["nitrogen", ] * 100)               # its 95% CI
predict(m2, data.frame(nitrogen = c(50, 150)), type = "response")   # probabilities

# ---- 6) Binomial, $k$ of $n$: seedlings emerged out of seeds sown ----

# $\text{emerged} \sim \text{Binomial}(\text{sown}, p)$, same logit
# link. Put cbind(successes, failures) on the left, so a row with 120
# seeds weighs more than a row with 40.

# [CORE]
m3 <- glm(cbind(emerged, sown - emerged) ~ nitrogen, family = binomial, data = survey)
coef(summary(m3))
exp(coef(m3)["nitrogen"] * 100)   # odds ratio, per 100 kg N

# The same model, written as a proportion with weights:

# [HANDY]
m3b <- glm(emerged / sown ~ nitrogen, family = binomial, weights = sown, data = survey)
coef(m3b)

# ---- 7) Poisson: beetles per trap ------------------------------------

# $\text{beetles} \sim \text{Poisson}(\mu)$, $\log\mu = \beta_0 +
# \beta_1\,\text{N}$. $e^{\beta_1}$ is a rate ratio. Every field had
# the same effort (one trap, one week).

# [CORE]
m4 <- glm(beetles ~ nitrogen, family = poisson, data = survey)
coef(summary(m4))
exp(coef(m4)["nitrogen"] * 100)        # rate ratio, per 100 kg N
deviance(m4) / df.residual(m4)         # about 1: Poisson spread fits

# ---- 8) Overdispersion: is Poisson enough for the aphids? ------------

# Poisson forces variance = mean. Residual deviance / df well above 1
# means the counts are more spread out than that, and the Poisson
# standard errors are too small.

# [CORE]
m5p <- glm(aphids ~ nitrogen, family = poisson, data = survey)
deviance(m5p) / df.residual(m5p)

# [DEMO: the Poisson 95% band is far too narrow for the aphids]
g <- data.frame(nitrogen = 0:200)
g$mu <- predict(m5p, g, type = "response")
g$lo <- qpois(0.025, g$mu); g$hi <- qpois(0.975, g$mu)
ggplot(survey, aes(nitrogen, aphids)) +
  geom_ribbon(data = g, aes(nitrogen, ymin = lo, ymax = hi), inherit.aes = FALSE, fill = "red", alpha = 0.15) +
  geom_point(colour = "#3D7BEB") +
  geom_line(data = g, aes(nitrogen, mu), colour = "red") +
  theme_classic()
outside <- with(survey, aphids < qpois(0.025, fitted(m5p)) | aphids > qpois(0.975, fitted(m5p)))
mean(outside)   # share of fields outside the band (should be about 0.05)

# ---- 9) Negative binomial: aphids ------------------------------------

# $\text{aphids} \sim \text{NegBin}(\mu, \theta)$, $\log\mu = \beta_0 +
# \beta_1\,\text{N}$, $\text{Var} = \mu + \mu^2/\theta$. Same link and
# same reading of $\beta_1$ as Poisson; one more parameter for the
# spread.

# [CORE]
m5 <- glm.nb(aphids ~ nitrogen, data = survey)
coef(summary(m5))
m5$theta                               # smaller = more clumped
exp(coef(m5)["nitrogen"] * 100)        # rate ratio, per 100 kg N

# [HANDY]
# the standard error of the nitrogen slope, three ways
m5q <- glm(aphids ~ nitrogen, family = quasipoisson, data = survey)
c(poisson = coef(summary(m5p))[2, 2],
  quasipoisson = coef(summary(m5q))[2, 2],
  negbin = coef(summary(m5))[2, 2])

# ---- 10) Counts per effort: the offset -------------------------------

# Scouts scored different numbers of leaves per field, and more on
# high-nitrogen fields. The question is lesions per leaf. An offset
# puts $\log(\text{leaves})$ in the linear predictor with its
# coefficient fixed at 1:

# $$\log\mu = \log(\text{leaves}) + \beta_0 + \beta_1\,\text{N}
# \quad\Longleftrightarrow\quad \log\frac{\mu}{\text{leaves}} = \beta_0
# + \beta_1\,\text{N}$$

# [CORE]
m6 <- glm(lesions ~ nitrogen + offset(log(leaves)), family = poisson, data = survey)
coef(summary(m6))
exp(coef(m6)["nitrogen"] * 100)        # lesions per leaf, ratio per 100 kg N
deviance(m6) / df.residual(m6)

# Leave the offset out and nitrogen also gets the credit for the extra
# leaves scored:

# [CORE]
m6x <- glm(lesions ~ nitrogen, family = poisson, data = survey)
exp(coef(m6x)["nitrogen"] * 100)       # no offset: too steep
cor(survey$nitrogen, survey$leaves)    # because effort rises with nitrogen

# Predictions from an offset model need the effort column. leaves = 1
# gives the rate per leaf:

# [HANDY]
predict(m6, data.frame(nitrogen = c(50, 150), leaves = 1), type = "response")

# The negative binomial takes an offset the same way: glm.nb(y ~ x +
# offset(log(effort))).

# ---- 11) Gamma: weed biomass -----------------------------------------

# Positive, continuous, right-skewed, spread growing with the mean.
# $\text{weeds} \sim \text{Gamma}$, $\log\mu = \beta_0 +
# \beta_1\,\text{N}$, $\text{Var} = \mu^2/k$.

# [CORE]
m7 <- glm(weed_biomass ~ nitrogen, family = Gamma(link = "log"), data = survey)
coef(summary(m7))
exp(coef(m7)["nitrogen"] * 100)        # ratio of mean biomass, per 100 kg N
1 / summary(m7)$dispersion             # the shape k (about 2 in the simulation)

# The usual alternative, lm(log(y) ~ x), models the mean of $\log y$ (a
# geometric mean), not the mean itself, and fails on zeros:

# [HANDY]
exp(coef(lm(log(weed_biomass) ~ nitrogen, data = survey))["nitrogen"] * 100)

########################################################################
# Reading and plotting any GLM
########################################################################

# ---- 12) One rule per link -------------------------------------------

# | link | model | $x + 1$ does | |---|---|---| | identity | $\mu =
# \beta_0 + \beta_1 x$ | adds $\beta_1$ to the mean | | log | $\mu =
# e^{\beta_0} e^{\beta_1 x}$ | multiplies the mean by $e^{\beta_1}$ | |
# logit | $\text{odds} = e^{\beta_0} e^{\beta_1 x}$ | multiplies the
# odds by $e^{\beta_1}$ |

# [CORE]
new <- data.frame(nitrogen = c(50, 150), leaves = 1)
predict(m2, new, type = "link")        # log-odds
predict(m2, new, type = "response")    # probability

# ---- 13) The fitted curve with a confidence band, for any family -----

# Build the band on the link scale (where it is symmetric), then
# back-transform.

# [HANDY]
glm_curve <- function(model, data, x = "nitrogen", ...) {
  g <- data.frame(seq(min(data[[x]]), max(data[[x]]), length.out = 200), ...)
  names(g)[1] <- x
  p <- predict(model, g, type = "link", se.fit = TRUE)
  inv <- family(model)$linkinv
  g$fit <- inv(p$fit)
  g$lo <- inv(p$fit - 1.96 * p$se.fit)
  g$hi <- inv(p$fit + 1.96 * p$se.fit)
  g
}
band <- glm_curve(m4, survey)
ggplot(survey, aes(nitrogen, beetles)) +
  geom_ribbon(data = band, aes(nitrogen, ymin = lo, ymax = hi), inherit.aes = FALSE, alpha = 0.2) +
  geom_point(colour = "#3D7BEB") +
  geom_line(data = band, aes(nitrogen, fit), colour = "red", linewidth = 1) +
  theme_classic()

# ---- 14) All seven at once -------------------------------------------

# [HANDY]
models <- list(yield = m1, lodged = m2, emerged = m3, beetles = m4,
               aphids = m5, lesions = m6, weeds = m7)
data.frame(
  family = sapply(models, function(m) family(m)$family),
  link   = sapply(models, function(m) family(m)$link),
  beta1  = sapply(models, function(m) coef(m)["nitrogen"]),
  per100 = sapply(models, function(m) {
    b <- 100 * coef(m)["nitrogen"]
    if (family(m)$link == "identity") b else exp(b)
  }),
  p = sapply(models, function(m) coef(summary(m))["nitrogen", 4]))

# per100 is a difference for yield (t/ha), an odds ratio for lodged and
# emerged, and a ratio of means for the rest.

# ---- 15) Beyond these seven ------------------------------------------

# | response | function | |---|---| | ordered categories (none < mild <
# severe) | MASS::polr() | | three or more unordered categories |
# nnet::multinom() | | too many zeros | glmmTMB::glmmTMB(..., ziformula
# = ~ 1) | | continuous proportions (% cover) | glmmTMB::glmmTMB(...,
# family = beta_family()) | | any of these with plots nested in farms |
# lme4::glmer(..., family = ...) |