Why Generalized Linear Models?

Non-normal outcomes: binary, counts, proportions, rates.

GLM Framework: the idea first

Two plain sentences before any density:

• \(\theta_i\) is a re-parameterization of the mean, called the natural parameter.
• \(\phi\) is a scale / dispersion constant that sets the variance level.

Everything you need follows from one cumulant function:

GLM Framework: the exponential family

Each response \(Y_i\) has a density in the exponential family. Using FLW’s notation (eq. 11.3, p. 328):

\[ f(y_i \mid \theta_i,\phi) = \exp\!\left[ \frac{y_i\theta_i - a(\theta_i)}{\phi} + b(y_i,\phi) \right] \]

with \(\mathbb{E}(Y_i)=\mu_i=a'(\theta_i)\) and \(\mathrm{Var}(Y_i)=\phi\,a''(\theta_i)\).

Exponential Family: a concrete example

For Poisson with mean \(\mu\) (here \(\phi=1\)):

• Natural parameter: \(\theta = \log(\mu)\)
• Cumulant function: \(a(\theta) = e^\theta = \mu\)
• \(a'(\theta) = e^\theta = \mu\) (the mean)
• \(a''(\theta) = e^\theta = \mu\) (the variance)

Insight: for Poisson, mean \(=\) variance. This falls out of the family automatically.

GLM Components

Here \(X_{ij}\) is the row of covariates (\(1\times p\)), so the mean model is \(X_{ij}\boldsymbol\beta\).

Canonical Links

Estimation: maximum likelihood

GLM coefficients \(\boldsymbol\beta\) are estimated by maximum likelihood (FLW Section 11.7). FLW states only this; the algorithm below is standard GLM theory.

IRLS: The Intuition

Think of IRLS as “smart” weighted regression that learns as it goes:

Key insight: higher-variance observations get lower weights, so uncertain points have less influence.

IRLS: Why Iteration?

Unlike OLS, GLMs have a nonlinear link between parameters and predictions.

OLS: one-shot solution \((\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{Y}\).

GLM: weights depend on \(\mu_i\), which depends on \(\boldsymbol\beta\), which depends on the weights…

Solution: iterate to convergence. Start somewhere reasonable, improve step by step. In R, glm() handles this automatically (it warns if there are convergence problems).

Interpreting Logistic Coefficients

Model on the log-odds scale:

\[ \log\!\left(\frac{\pi}{1-\pi}\right) = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p \]

Case Study: TLC Lead-Exposure Trial (real data)

tlc <- tlc |>
  dplyr::mutate(elevated_w6 = as.integer(w6 >= 20))   # 1 = elevated lead at wk 6
table(Treatment = tlc$trt, Elevated_wk6 = tlc$elevated_w6)

Logistic GLM: Fit

fit_logit <- glm(elevated_w6 ~ trt + w0, data = tlc, family = binomial)
summary(fit_logit)

Logistic: Odds Ratios and CIs

or_tab <- broom::tidy(fit_logit, conf.int = TRUE) |>
  dplyr::mutate(OR = exp(estimate),
                OR_low = exp(conf.low),
                OR_high = exp(conf.high)) |>
  dplyr::select(term, estimate, OR, OR_low, OR_high, p.value)
or_tab

Plain-Language Interpretation: Logistic Output

In a sentence: adjusting for baseline lead, the succimer odds ratio for remaining elevated at week 6 is 0.18 (below 1, i.e. lower odds than placebo), and the baseline-lead odds ratio is 1.33 (above 1, i.e. higher baseline raises the odds).

From Coefficients to Predicted Probabilities

Logistic inverse link (row form):

\[ \hat\pi_{ij} = \frac{e^{X_{ij}\hat{\boldsymbol\beta}}}{1 + e^{X_{ij}\hat{\boldsymbol\beta}}} \]

In R, predict(fit, type = "response") returns \(\hat\pi_{ij}\).

tlc$pred_prob <- predict(fit_logit, type = "response")
head(tlc[, c("trt", "w0", "elevated_w6", "pred_prob")])

Visualize Fitted Probabilities

newdat <- expand.grid(
  trt = levels(tlc$trt),
  w0  = seq(min(tlc$w0), max(tlc$w0), length = 100)
)
newdat$pred_prob <- predict(fit_logit, newdata = newdat, type = "response")
ggplot(newdat, aes(w0, pred_prob, color = trt)) +
  geom_line(linewidth = 1.2) +
  labs(y = "Predicted P(elevated at wk 6)", x = "Baseline lead (wk 0)",
       title = "Logistic Fitted Probabilities: TLC Trial",
       color = "Treatment") +
  theme_minimal()

Static SAS: PROC LOGISTIC (reference only)

The same TLC logistic model in SAS (FLW Section 11.6 uses PROC GENMOD/LOGISTIC):

proc logistic data=tlc descending;
   class trt (ref="Placebo") / param=ref;
   model elevated_w6 = trt w0;
   /* odds ratios printed automatically; oddsratio for custom contrasts */
run;

Using Probabilities for Classification

Choose a threshold \(t\) and predict:

\[ \widehat Y = \begin{cases} 1, & \hat\pi \ge t \\ 0, & \text{otherwise} \end{cases} \]

Default \(t=0.5\), but it can be tuned to the application’s costs.

tlc$pred_class_05 <- ifelse(tlc$pred_prob >= 0.5, 1, 0)
table(Predicted = tlc$pred_class_05, Observed = tlc$elevated_w6)

Sensitivity, Specificity, and FPR

Varying the threshold traces an ROC curve (TPR on the \(y\)-axis vs FPR on the \(x\)-axis).

ROC Curve and AUC

suppressPackageStartupMessages(library(pROC))
roc_obj <- roc(tlc$elevated_w6, tlc$pred_prob, quiet = TRUE)
plot(roc_obj, main = paste0("TLC Logistic ROC (AUC = ", round(auc(roc_obj), 3), ")"))

Interpreting AUC

Calibration: Predicted vs Observed

dd <- tlc |>
  dplyr::mutate(dec = cut(pred_prob,
                          breaks = quantile(pred_prob, probs = seq(0, 1, 0.1), na.rm = TRUE),
                          include.lowest = TRUE)) |>
  dplyr::group_by(dec) |>
  dplyr::summarise(mean_pred = mean(pred_prob),
                   obs_rate  = mean(elevated_w6),
                   n = dplyr::n(), .groups = "drop")
ggplot(dd, aes(mean_pred, obs_rate, size = n)) +
  geom_point(alpha = 0.8) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  labs(title = "TLC Logistic Calibration (deciles of predicted risk)",
       x = "Mean predicted probability", y = "Observed event rate") +
  theme_minimal()

Logistic Diagnostics: Residuals

res_df <- data.frame(
  fitted = fitted(fit_logit),
  dev    = residuals(fit_logit, type = "deviance"),
  pear   = residuals(fit_logit, type = "pearson")
)
p1 <- ggplot(res_df, aes(fitted, dev)) + geom_point(alpha = .6) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Deviance Residuals vs Fitted") + theme_bw()
p2 <- ggplot(res_df, aes(fitted, pear)) + geom_point(alpha = .6) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Pearson Residuals vs Fitted") + theme_bw()
print(p1); print(p2)

Influence and Leverage

lev <- hatvalues(fit_logit)
cd  <- cooks.distance(fit_logit)
par(mfrow = c(1, 2))
plot(lev, type = "h", main = "Leverage (hat values)", ylab = "hat")
abline(h = 2 * mean(lev), col = "red", lty = 2)
plot(cd, type = "h", main = "Cook's Distance", ylab = "Cook's D")
abline(h = 4 / length(cd), col = "red", lty = 2)

Goodness-of-Fit: Hosmer-Lemeshow

if (requireNamespace("ResourceSelection", quietly = TRUE)) {
  hl <- ResourceSelection::hoslem.test(tlc$elevated_w6, fitted(fit_logit), g = 10)
  hl
} else {
  cat("Install 'ResourceSelection' for the Hosmer-Lemeshow test.\n")
}

Logistic Models: Quick Checklist

* enrichment beyond FLW Ch. 11.

Transition to Count Outcomes

Counts (number of events) call for the Poisson family.

Case Study: Epilepsy Seizure-Count Trial (real data)

str(ep[, c("id", "trtf", "age", "base", "total")])

Fit Poisson GLM (with offset)

ep$logexp <- log(8)   # 8 post-randomization weeks of observation
fit_pois <- glm(total ~ trtf + base + age + offset(logexp),
                data = ep, family = poisson)
summary(fit_pois)

Poisson GLM: Interpret Rate Ratios

• \(e^{\hat\beta_{\text{Progabide}}} =\) 0.82
• \(e^{\hat\beta_{\text{base}}} =\) 1.023 per additional baseline seizure

Interpretation: adjusting for baseline count and age, progabide is associated with a seizure rate 0.82 times that of placebo (a 18% lower rate point estimate). Each extra baseline seizure multiplies the rate by 1.023.

Plain-Language Interpretation: Poisson Output

Note on the offset: the offset makes us model rates (events per unit time), not raw counts, so subjects with different follow-up are compared fairly.

Visualize Fitted Rates

newd <- expand.grid(
  age   = seq(min(ep$age), max(ep$age), length.out = 100),
  trtf  = levels(ep$trtf),
  base  = median(ep$base),
  logexp = log(8)
)
newd$pred <- predict(fit_pois, newdata = newd, type = "response")
ggplot(newd, aes(age, pred, color = trtf)) +
  geom_line(linewidth = 1.1) +
  labs(y = "Predicted seizures / 8 weeks", x = "Age",
       title = "Poisson GLM: Predicted Seizure Count (at median baseline)",
       color = "Treatment") +
  theme_minimal()

Poisson Diagnostics: Pearson Residuals

pp <- tibble(fitted = fitted(fit_pois),
             rpear  = residuals(fit_pois, type = "pearson"))
ggplot(pp, aes(fitted, rpear)) +
  geom_point(alpha = .6) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Epilepsy Poisson: Pearson Residuals vs Fitted") +
  theme_bw()

Checking for Overdispersion

phi_hat <- sum(residuals(fit_pois, type = "pearson")^2) / fit_pois$df.residual
phi_hat

If \(\hat\phi > 1\), the variance exceeds the mean (overdispersion). Here \(\hat\phi\) is far above 1, so the Poisson SEs are too small (invalid inference).

Remedy: quasi-Poisson or negative binomial.

Quasi-Poisson Model

fit_qp <- glm(total ~ trtf + base + age + offset(logexp),
              data = ep, family = quasipoisson)
summary(fit_qp)

Compare Poisson vs Quasi-Poisson

Negative Binomial

fit_nb <- glm.nb(total ~ trtf + base + age + offset(logexp), data = ep)
summary(fit_nb)

Static SAS: PROC GENMOD with offset (reference only)

The epilepsy Poisson rate model in SAS (FLW Section 11.6, Table 11.9 style):

data ep; set ep; logexp = log(8); run;

proc genmod data=ep;
   class trtf (ref="Placebo") / param=ref;
   model total = trtf base age / dist=poisson link=log offset=logexp;
   /* for overdispersion, add  dscale  to scale the SEs by sqrt(Pearson/df) */
run;

/* negative binomial alternative */
proc genmod data=ep;
   class trtf (ref="Placebo") / param=ref;
   model total = trtf base age / dist=negbin link=log offset=logexp;
run;

Ordinal Response Setup

Ordered categories (e.g. symptom severity: Low \(<\) Med \(<\) High).

Proportional-odds model (POM), R/polr convention:

\[ \log\!\left(\frac{P(Y\le k)}{P(Y>k)}\right)=\kappa_k - X_{ij}\boldsymbol{\beta}, \qquad k=1,\dots,K-1 \]

Simulate Ordinal Data

set.seed(667)
N <- 250
dat_ord <- tibble(id = 1:N, trt = rbinom(N, 1, 0.5), time = rnorm(N)) %>%
  mutate(eta = -0.5 + 1.1 * trt + 0.6 * time,
         p_low = plogis(1 - eta),
         p_med = plogis(3 - eta) - plogis(1 - eta),
         u = runif(N),
         ord = case_when(u < p_low ~ "L",
                         u < p_low + p_med ~ "M",
                         TRUE ~ "H"),
         trt = factor(trt, labels = c("Ctl", "Tx")),
         ord = factor(ord, ordered = TRUE, levels = c("L", "M", "H")))
head(dat_ord)

Fit Ordinal Logistic Model

fit_polr <- MASS::polr(ord ~ trt * time, data = dat_ord, Hess = TRUE)
summary(fit_polr)

What the Coefficients Mean

Models the log-odds of lower categories (\(\le k\)) vs higher (\(>k\)). Under the minus-sign (polr) convention:

Example: if \(\hat\beta_{\text{Tx}}=0.8\), then OR \(=e^{0.8}=2.23\): treatment has 2.23x the odds of being in a higher category, and one \(\beta\) applies across all cutpoints (the proportional-odds assumption).

Comparing Ordinal vs Multinomial

suppressPackageStartupMessages(library(nnet))
fit_multi <- multinom(ord ~ trt * time, data = dat_ord, trace = FALSE)
# AIC comparison: POM (shared beta) vs multinomial (separate beta per category)
data.frame(model = c("Ordinal POM", "Multinomial"),
           params = c(length(coef(fit_polr)) + length(fit_polr$zeta),
                      length(as.vector(coef(fit_multi)))),
           AIC = c(AIC(fit_polr), AIC(fit_multi)))

Is Proportional Odds Reasonable?

A quick check: fit separate binary logits at each cutpoint and compare slopes (if PO holds, the slopes should be similar). A formal Brant test is also available.

# Cutpoint 1: P(Y > L);  Cutpoint 2: P(Y > M)
ge_M <- as.integer(dat_ord$ord >= "M")   # not Low
ge_H <- as.integer(dat_ord$ord >= "H")   # High
s1 <- coef(glm(ge_M ~ trt + time, data = dat_ord, family = binomial))
s2 <- coef(glm(ge_H ~ trt + time, data = dat_ord, family = binomial))
round(rbind(`cut: >L` = s1, `cut: >M` = s2), 3)

Predicted Category Probabilities

newo <- expand.grid(trt = c("Ctl", "Tx"), time = seq(-2, 2, length = 100))
preds <- predict(fit_polr, newdata = newo, type = "probs")
pred_df <- cbind(newo, preds) %>%
  pivot_longer(cols = c(L, M, H), names_to = "level", values_to = "prob")
ggplot(pred_df, aes(time, prob, color = level)) +
  geom_line(linewidth = 1) + facet_wrap(~trt) +
  labs(title = "POM Predicted Category Probabilities", y = "Probability") +
  theme_minimal()

Diagnostic Goals for GLMs

Single reference table (used by both the logistic and Poisson diagnostics above):

Deviance and Pearson Residuals

\[ r_i^{(P)} = \frac{y_i - \hat\mu_i}{\sqrt{\hat V(\hat\mu_i)}}, \qquad r_i^{(D)} = \operatorname{sign}(y_i-\hat\mu_i)\sqrt{2\left[\ell(y_i)-\ell(\hat\mu_i)\right]} \]

Assessing Model Fit with Deviance

Residual deviance measures fit relative to the saturated model (the model with one parameter per observation, i.e. a perfect fit – the deviance is twice the log-likelihood gap between your model and that perfect fit).

Compare nested models via the likelihood-ratio test:

\[ D = 2(\ell_{\text{full}} - \ell_{\text{reduced}}) \sim \chi^2_{df} \]

For nested models, smaller deviance indicates better fit. (This is FLW-native, Section 11.7.)

Choosing the Link Function

Zero-Inflated and Hurdle Models

When counts have far more zeros than Poisson allows, use a two-part model:

# Illustrative ZIP on the epilepsy counts (note: these are not zero-inflated,
# so this is a demonstration of syntax, not a recommended model here).
if (requireNamespace("pscl", quietly = TRUE)) {
  zi <- pscl::zeroinfl(total ~ trtf + base | 1, data = ep, dist = "poisson")
  round(coef(zi), 3)
} else {
  cat("Install 'pscl' for zero-inflated models.\n")
}

Quasi-Likelihood and Robust (Sandwich) SEs

For misspecified-variance or correlated data, robust (sandwich) SEs give valid inference even if the variance model is wrong.

if (requireNamespace("sandwich", quietly = TRUE) &&
    requireNamespace("lmtest", quietly = TRUE)) {
  lmtest::coeftest(fit_pois, vcov = sandwich::sandwich)
} else {
  cat("Install 'sandwich' and 'lmtest' for robust SEs.\n")
}

GEE Connection (Preview)

Single-level GLMs assume independent observations. Longitudinal data have within-subject correlation; Ch. 12-13 relax independence using GEE.

The GEE working covariance (NOT the true \(\mathrm{Var}(\mathbf{Y}_i)\)) is

\[ V_i = \mathbf{A}_i^{1/2}\,\operatorname{Corr}(\mathbf{Y}_i)\,\mathbf{A}_i^{1/2}, \]

where \(\mathbf{A}_i\) is the diagonal of marginal variances and \(\operatorname{Corr}(\mathbf{Y}_i)\) is a working correlation.

Model Selection and AIC

\[ \text{AIC} = -2\ell + 2p \]

Smaller AIC = better fit/complexity tradeoff. Use the LRT (deviance) for nested models; AIC for non-nested comparison.

Common Pitfalls

Further Reading

• Fitzmaurice, Laird & Ware (2011), Applied Longitudinal Analysis, Ch. 11 (course text)
• Dobson & Barnett (2018), An Introduction to Generalized Linear Models (GLM theory, IRLS)
• McCullagh & Nelder (1989), Generalized Linear Models (exponential family, links)
• Agresti (2015), Foundations of Linear and Generalized Linear Models (ordinal/multinomial)
• Hosmer, Lemeshow & Sturdivant (2013), Applied Logistic Regression (ROC, calibration, H-L)
• Zeileis, Kleiber & Jackman (2008), Regression Models for Count Data in R (ZIP, hurdle)

GLM Summary: Key Takeaways

Next Lecture: Marginal Models (Ch. 12); GEE estimation follows in Ch. 13.

Appendix A: IRLS Algorithm Details

Each iteration solves a weighted least-squares problem:

\[ \boldsymbol{\beta}^{(t+1)} = (\mathbf{X}^\top W^{(t)} \mathbf{X})^{-1} \mathbf{X}^\top W^{(t)} z^{(t)} \]

with working response and weights

\[ z^{(t)} = \eta^{(t)} + (y - \mu^{(t)}) \frac{d\eta}{d\mu}, \qquad W^{(t)} = \operatorname{diag}\!\left[\left(\tfrac{d\mu}{d\eta}\right)^2 \big/ v(\mu)\right]. \]

Here \(d\mu/d\eta\) is the slope of the inverse link; \(v(\mu)\) is the variance function.

Appendix B: Score Function and Fisher Information

For the canonical link (rows \(X_{ij}\), \(1\times p\)):

\[ U(\boldsymbol{\beta}) = \sum_i \sum_j X_{ij}^\top (Y_{ij} - \mu_{ij}), \qquad I(\boldsymbol{\beta}) = \sum_i \sum_j \frac{X_{ij}^\top X_{ij}}{\operatorname{Var}(Y_{ij})} . \]

Used for standard errors and likelihood inference.

Appendix C: Exponential Family Examples

Using FLW’s labeling: cumulant function \(a(\theta)\), data-only term \(b(y,\phi)\), scale \(\phi\).

Recall \(\mathbb{E}(Y)=a'(\theta)\) and \(\mathrm{Var}(Y)=\phi\,a''(\theta)\).