Generalized Linear Models

Today

• Count response
• Poisson distribution
• Log-count regression
• Generalized linear models
• Poisson regression
• Maximum likelihood
• Rate ratios
• Inference
• Overdispersion

Distilled from the chapter’s six learning objectives, each compressed to a noun phrase: three components and special cases; linear regression’s failure and log-count regression; the score equations; coefficient interpretation and Wald inference; the drop-in-deviance test; Pearson residuals and overdispersion. The title slide carries the chapter title.

Timing: one 50-minute period. Optional slides, in the order they would go if time is short: Fitting, Conclusion. If time is still short, drop in this order: Observed hourly means, Variance function, Interaction test (keep the equations on Drop-in-deviance test).

Count Response

Bike rentals

ISLR2::Bikeshare: 8,645 hours of one year in Washington, D.C.

• Response variable: bikers, the number of bikes rented in the hour
• Explanatory variables:
  • hr: hour of the day, 0 to 23
  • workingday: working day or not
  • mnth: month
  • temp: temperature, rescaled to \([0, 1]\)
  • weathersit: weather

A count is a nonnegative integer, its mean must be positive and its spread typically grows with its mean, unlike the quantitative response of the regression lectures and the binary one of the logistic regression lectures. The question is how the hourly count depends on time of day, working day, month, temperature and weather. casual and registered are pieces of the response and no model uses them. 115 hours of the year have no record (365 days of 24 hours, against 8,645 records).

Hourly counts

• Distribution
• Hourly profile
• Mean and variance

Distribution: counts from 1 to 651, median 109, mean 143.8, strongly right-skewed. Hourly profile: working days peak twice, at hour 8 and hour 17, the commute; non-working days rise to a midday plateau; the two profiles differ in shape, not just in level, and the spread is larger where the mean is higher. Mean and variance: one point per hour-by-working-day group (48 of them) at its mean and variance on log scales; the variance climbs with the mean, on a line of slope 1.75 on these axes, and all 48 points lie above the dashed line where the variance equals the mean; month, temperature and weather vary within each group, so part of each group’s variance is variation they could explain.

The third tab is the same image as the chapter’s figure that opens the Generalized Linear Models section; say so aloud when the deck reaches that section.

Poisson distribution

\[P(Y = y) = \frac{\mu^{y} e^{-\mu}}{y!}, \qquad y = 0, 1, 2, \ldots\]

\(\mu > 0\) — mean

\[E[Y] = Var[Y] = \mu\]

The dashed line of the previous slide is the variance of the standard model for a count. The chapter writes mu for the Poisson mean (the symbol the GLMs below use for a mean) where ISLR2 writes lambda.

Poisson distribution, slider

viewof mu = Inputs.range([1, 40], {value: 1, step: 0.5, label: "Mean μ", width: 640})

{
  const pmf = transpose(poisson_pmf).filter(d => Math.abs(d.mu - mu) < 1e-9);
  const sd = Math.sqrt(mu);
  const xmax = Math.max(10, Math.ceil(mu + 4 * sd));
  const ymax = 1.15 * Math.max(...pmf.map(d => d.prob));
  const plot = Plot.plot({
    width: 900, height: 340, style: {fontSize: "16px"},
    x: {domain: [-0.5, xmax + 0.5], label: "Count y"},
    y: {domain: [0, ymax], label: "P(Y = y)"},
    marks: [
      Plot.rect([{x1: mu - sd, x2: mu + sd}],
                {x1: "x1", x2: "x2", y1: 0, y2: ymax,
                 fill: "#0072B2", fillOpacity: 0.12}),
      Plot.ruleX(pmf.filter(d => d.y <= xmax),
                 {x: "y", y1: 0, y2: "prob", strokeWidth: 2}),
      Plot.dot(pmf.filter(d => d.y <= xmax),
               {x: "y", y: "prob", r: 2.5, fill: "black"}),
      Plot.ruleX([mu], {stroke: "#0072B2", strokeDasharray: "6,4",
                        strokeWidth: 2})
    ]
  });
  return html`<p style="text-align: center"><b>mean = variance = ${mu.toFixed(1)}</b>,
    standard deviation = ${sd.toFixed(2)}, P(Y = 0) = ${Math.exp(-mu).toFixed(3)}</p>
    <div style="display: flex; justify-content: center">${plot}</div>`;
}

• Mean \(\mu\) from 1 to 40
• Band: one standard deviation, \(\sqrt{\mu}\)

Work the slider live. At a mean near 1 most of the probability sits at 0, 1 and 2 and the distribution is strongly right-skewed; as mu grows it moves right, becomes nearly symmetric and widens, its standard deviation sqrt(mu) growing with the mean. Both axes rescale, so read the spread and heights from the axis labels. The busiest bike group has mean 395, so the Poisson standard deviation is 19.9, against an observed standard deviation of 164.4 (variance 27,031): the gap the mean-variance plot showed for every group.

The figure needs an internet connection. While the slider has focus reveal.js ignores the arrow keys, so click an empty part of the slide to go on.

Linear regression fails

• 895 of 8,645 hours have a negative fitted mean
• Residuals fan out to the right

The linear regression uses the 61 features of bike_formula (hour dummies, working day, their products, month dummies, temperature in degrees Celsius, weather dummies). Three failures. Negative fitted means: 895 hours, 10.4%, the smallest at -140.9 bikers, and a mean count cannot be negative. Variance growing with the mean: residual standard deviation 37.6 where the fitted value is below 50 and 90.5 where it is above 300, while every lm() test and interval assumes one common variance. Additive effects: a winter month, a cold hour or rain shifts every hour by the same number of riders, which at night is enough to push the fit below zero (460 of the 895 negative fitted values are in clear weather, and the median temperature of all 895 is 3.3 degrees Celsius against 15.5 for all hours).

Log-count linear regression

\[\log Y_i \mid x_i \stackrel{ind}{\sim} N\left(x_i\gamma, \sigma^2\right)\]

The regression lectures’ remedy for a positive response whose spread grows with its mean is least squares on log Y_i with the same features; gamma and sigma squared are this model’s coefficients and error variance. It needs every y_i > 0, and no hour has zero bikers, so log y_i is defined for all 8,645 hours. Its smallest back-transformed fitted value is 0.75 bikers, so no fitted value is negative, and the residuals no longer fan out to the right: they narrow as the mean grows (standard deviation 0.55 below 50 and 0.23 above 300 on the log scale), so one common sigma squared still does not describe the data.

Exponentiating a fitted value gives e to the x_i gamma-hat, which estimates e to the E[log Y_i | x_i] and not the mean count: by Jensen’s inequality it is at most the mean, and it is the median of Y_i given x_i under the normal model. Averaged within each of the 48 hour-by-working-day combinations the back-transformed fits fall below the observed mean in 45, with a median ratio of 0.93. If the contrast needs more time this is where to spend it; the Jensen derivation is in the chapter.

Generalized Linear Models

The chapter’s section opens with the mean-variance plot repeated; the deck does not repeat the image, since the third tab of Hourly counts already showed it. Say “the same plot, variance far above constant” aloud.

Three components

\[Y_i \mid x_i \stackrel{ind}{\sim} \text{a distribution with mean } \mu_i, \qquad g(\mu_i) = \eta_i, \qquad \eta_i = x_i\beta\]

• \(\mu_i = E[Y_i \mid x_i]\) — mean of the response
• \(\eta_i\) — linear predictor
• \(g\) — link function, strictly monotone

The Problems in Logistic Regression lecture wrote linear and logistic regression as one distributional statement each; both have the same three parts, a distribution, a linear predictor and a link, and a GLM keeps them while letting each vary. The distribution of Y_i given x_i is the random component; dummy variables, polynomials, step functions and interactions enter the linear predictor exactly as they did before. The link connects the two, mu_i = g^-1(eta_i), and g^-1 carries the real line onto the values a mean can take, (0, 1) for a probability and (0, infinity) for a count, so every beta gives a valid mean. This mu_i is the mean of observation i, unrelated to the class means mu_k of the discriminant analysis lectures; ISLR2 writes eta for the link itself, these notes write g. For the bike counts the linear regression is the GLM with a normal random component and the identity link, which does nothing to keep x_i beta-hat positive, so 895 fitted means fell below zero.

Special cases

Model	Response	Distribution of \(Y_i\)	Link \(g(\mu_i)\)	Mean \(\mu_i\)	Variance	glm() family
Linear regression	Quantitative	\(N(\mu_i, \sigma^2)\)	identity, \(\mu_i\)	\(\eta_i\)	\(\sigma^2\)	gaussian
Logistic regression	Binary (\(0/1\))	\(\text{Bernoulli}(\mu_i)\)	logit, \(\log\left[\mu_i/(1 - \mu_i)\right]\)	\(e^{\eta_i}/\left(1 + e^{\eta_i}\right)\)	\(\mu_i(1 - \mu_i)\)	binomial
Poisson regression	Count	\(\text{Poisson}(\mu_i)\)	log, \(\log \mu_i\)	\(e^{\eta_i}\)	\(\mu_i\)	poisson

For each model the table gives the response type, the distribution of Y_i, the link, the mean as a function of eta_i, the variance and the glm() family that fits it (family = gaussian gives the same estimates as lm()). In logistic regression the mean mu_i = P(Y_i = 1 | x_i) is the probability the earlier lectures wrote p(x_i), so the credit card default lectures are this table’s logistic regression row.

Variance function

\[Var[Y_i \mid x_i] = \phi\, V(\mu_i)\]

• \(V\) — variance function
• \(\phi\) — dispersion parameter

Distribution	\(V(\mu)\)	\(\phi\)
Normal	\(1\)	\(\sigma^2\)
Bernoulli	\(\mu(1 - \mu)\)	\(1\)
Poisson	\(\mu\)	\(1\)

Third to drop if time is short (the table on Special cases already has a Variance column). Each distribution fixes how the variance depends on the mean; linear regression’s constant variance is the special case V(mu) = 1, the assumption the fanning residuals contradicted; the other two models require the variance to change with the mean, as the bike counts’ variance did. The Bernoulli and the Poisson have no variance parameter separate from the mean, so phi = 1, which is the assumption the overdispersion slides test.

Transforming the mean

\[\text{log-count linear regression: } E[\log Y_i \mid x_i] = x_i\gamma, \qquad \text{Poisson regression: } \log E[Y_i \mid x_i] = x_i\beta\]

A GLM transforms the mean, not the response

A log link never takes \(\log y_i\): a count of 0 is fine

Both models put the linear predictor on a log scale, of different quantities. The log-count model takes the logarithm of y_i, so it fails for a count of zero (bikers has none, but counts in general can); the log link takes the logarithm of the mean, never of y_i.

Poisson Regression

Observed hourly means

Rain scales the profile: 0.54 times clear at 3 a.m., 0.63 times at 5 p.m.

Second to drop if time is short. Each weather row repeats the hourly profile of the clear row at a lower level; on working days the rain/snow mean is 0.54 times the clear mean at 3 a.m. and 0.63 times at 5 p.m., while the differences are 2.2 and 161.1 riders, so rain is a multiplicative effect, not a constant number of riders at every hour.

Log-linear model

\[Y_i \mid x_i \stackrel{ind}{\sim} \text{Poisson}(\mu_i), \qquad \log \mu_i = x_i\beta = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}\]

Poisson random component, log link

\[\mu_i = e^{x_i\beta} = e^{\beta_0} \times e^{\beta_1 x_{i1}} \times \cdots \times e^{\beta_p x_{ip}}\]

Features add on the log scale, multiply on the scale of the mean

Poisson regression is the GLM with a Poisson random component and a log link; the mean is positive for every beta and every row x_i. A weather dummy multiplies every hour’s mean by the same factor e to the beta_j, the scaling of the rain/snow row in the previous figure.

Maximum likelihood estimation

\[\ell(\beta) = \sum_{i=1}^n \left[y_i \log \mu_i - \mu_i - \log(y_i!)\right] = \sum_{i=1}^n \left[y_i\, x_i\beta - e^{x_i\beta} - \log(y_i!)\right]\]

\[\sum_{i=1}^n x_{ij}\left(y_i - \mu_i\right) = 0, \qquad j = 0, 1, \ldots, p\]

• Same form as the logistic regression score equations
• No closed form: \(\mu_i = e^{x_i\beta}\) is nonlinear in \(\beta\)

Solved by iteratively reweighted least squares

The lecture’s derivation, the slowest slide after the Hourly counts tabset. The log-likelihood is the sum of the logs of the Poisson mass functions of independent observations; log(y_i!) does not involve beta and so does not affect the maximizer. Differentiate with respect to beta_j using the derivative of e to the x_i beta, x_ij times e to the x_i beta, and set to zero, with x_i0 = 1. They have exactly the form of the logistic regression score equations with mu_i in place of p(x_i); at the solution the residuals y_i - mu_i-hat, where mu_i-hat = e to the x_i beta-hat is the fitted mean, sum to zero and are orthogonal to every feature. Because mu_i is nonlinear in beta the equations in general have no closed-form solution: with only the hour-by-working-day dummies the solution is the cell means, but with a quantitative feature such as temperature it has none. glm() solves them with the algorithm it used for logistic regression, each iteration a weighted least squares problem with weights recomputed from the current estimate.

Fitting

bike_formula <- bikers ~ hr * working + mnth + temp_c + weather
bike_glm <- glm(bike_formula, data = bike, family = poisson)
printCoefmat(summary(bike_glm)$coefficients[
  c("(Intercept)", "workingYes", "temp_c", "weathercloudy/misty",
    "weatherrain/snow"), ])

                       Estimate  Std. Error  z value  Pr(>|z|)    
(Intercept)          3.36291632  0.01208439  278.286 < 2.2e-16 ***
workingYes          -0.91915334  0.01608340  -57.149 < 2.2e-16 ***
temp_c               0.01605829  0.00024332   65.998 < 2.2e-16 ***
weathercloudy/misty -0.06709450  0.00218683  -30.681 < 2.2e-16 ***
weatherrain/snow    -0.58063152  0.00406292 -142.910 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Optional: the score equations on Maximum likelihood estimation already state the consequence. Same formula as lm(), with family = poisson for the Poisson random component and the log link; the model has 62 coefficients and glm() reports convergence after 4 iterations. The intercept’s score equation makes the fitted means sum to the observed counts, 1,243,103 each; a dummy variable’s equation makes the fitted means average to the observed mean where that dummy is 1, so both models reproduce the observed mean of each of the 48 hour-by-working-day combinations (largest gap 6.9 x 10^-8 bikers). The printout is five of the 62 rows of summary(bike_glm)$coefficients: the estimate, standard error, Wald z statistic and p-value of the intercept and the working day, temperature and weather coefficients; the rate ratios slide interprets them.

Fitted means

Pooled over weather, linear and Poisson regression agree, since both reproduce each hour-by-working-day mean; they differ only in how they split each combination across the weather rows. In the rain/snow row linear regression subtracts its constant weather coefficient from every hour and its lines fall as low as -83.7 bikers, where the observed means never fall below 2.6; the Poisson lines scale the profile down instead and follow the points closely, overnight included. The log-count lines also scale the profile and stay positive (smallest overnight 2.1 bikers in the rain/snow row) but lie below the observed mean at 122 of the 144 points, against 68 for the Poisson lines, and total 1,202,378 bikers against the 1,243,103 observed: consistent with a fit that estimates e to the E[log Y_i | x_i] rather than the mean count.

Coefficient interpretation

\(e^{\beta_0}\) — mean count when every feature is zero

\[\log \mu(x_j + 1) - \log \mu(x_j) = \beta_j, \qquad \frac{\mu(x_j + 1)}{\mu(x_j)} = e^{\beta_j}\]

\(e^{\beta_j}\) — rate ratio: factor multiplying the mean count per one-unit increase in feature \(j\), other features fixed

The intercept is read first. Write mu(x) for the mean count at a row of feature values x = (1, x_1, …, x_p); at x = (1, 0, …, 0), log mu(x) = beta_0, so e to the beta_0 is the mean count at the reference levels of the dummy variables. Increasing feature j by one unit with every other feature fixed, every term of the linear predictor except the jth cancels. e to the beta_j is the odds ratio of logistic regression with the mean count in place of the odds; e to the beta_j = 1 means no association; c units multiply the mean by e to the c beta_j; for a dummy variable e to the beta_j is the mean at the indicated level divided by the mean at the reference level; with an interaction x_j x_l the rate ratio for x_j depends on x_l.

Rate ratios

Coefficient	Estimate	\(e^{\widehat\beta_j}\)
Intercept	3.3629	28.8733
Working day (at hour 0)	-0.9192	0.3989
Temperature (per degree C)	0.0161	1.0162
Cloudy/misty vs. clear	-0.0671	0.9351
Rain/snow vs. clear	-0.5806	0.5595

Intercept: 28.9 bikers at hour 0, non-working day, January, clear, 0 °C (linear regression: −10.9)

Rain/snow: ×0.560, a 44% reduction at every hour

5 °C: ×1.084, an increase of 8.4%

The intercept is read first: e to the beta_0-hat = 28.9 is a realistic winter night, unlike the linear regression’s -10.9 bikers for the same conditions. Because hour and working day interact, the working-day coefficient compares the two kinds of day at the reference hour only: at midnight a working day’s estimated mean is 0.40 times a non-working day’s, and at 8 a.m. the ratio adds the hr8:workingYes coefficient to the exponent, e to the (-0.919 + 2.442) = 4.59. Each additional degree Celsius multiplies the mean by 1.0162, holding hour, working day, month and weather fixed; rain or snow multiplies the mean by 0.560 relative to clear weather, the scaling in the rain/snow row of the fitted-means figure, and cloudy or misty weather by 0.935. The Fitting slide printed the same estimates on the log scale.

Inference

\[z = \frac{\widehat\beta_j - \beta_j^0}{\text{SE}\left(\widehat\beta_j\right)}\]

\(\approx N(0, 1)\) under \(H_0: \beta_j = \beta_j^0\)

\[\widehat\beta_j \pm z_{1-\alpha/2}\,\text{SE}\left(\widehat\beta_j\right)\]

Exponentiate the endpoints: interval for the rate ratio \(e^{\beta_j}\)

• \(|z|\) from 30.7 to 142.9, so every 95% interval excludes 1
• Working day at hour 0: (0.386, 0.412)
• 5 °C: (1.081, 1.086)
• Cloudy/misty: (0.931, 0.939)
• Rain/snow: (0.555, 0.564)

Inference repeats the logistic regression lectures: in large samples beta_j-hat is approximately normal, SE(beta_j-hat) comes from the curvature of the log-likelihood at beta-hat, and the Wald interval is what confint.default() returns; exponentiating its endpoints gives an interval for the rate ratio. For the four rate ratios of the table |z| ranges from 30.7 to 142.9. If the last click is too crowded, keep only the |z| range and the rain/snow interval.

Drop-in-deviance test

\[D = 2\left[\ell_{\text{sat}} - \ell\left(\widehat\beta\right)\right]\]

• \(D\) — deviance, not an indicator
• \(\ell_{\text{sat}}\) — saturated model, \(\mu_i = y_i\)

\[G^2 = D_r - D_f = -2\left[\ell\left(\widehat\beta_r\right) - \ell\left(\widehat\beta_f\right)\right]\]

Under \(H_0\): the \(k\) dropped coefficients are zero, \(G^2 \approx \chi^2_k\)

anova(reduced, full, test = "Chisq")

For a count response glm() reports the residual deviance. The saturated model gives each observation its own mean mu_i = y_i, the largest log-likelihood any mean model can reach, and does not involve beta, so for a reduced model nested in a full one G squared = D_r - D_f, the statistic of the multiple logistic regression lecture (for a binary response l_sat = 0, which is why that lecture wrote D = -2 l(beta-hat)). A group of coefficients, here every product of an hour dummy with the working-day indicator, is tested together; no single z-test asks whether working and non-working days share one hourly profile up to a constant factor.

Interaction test

bike_additive <- glm(bikers ~ hr + working + mnth + temp_c + weather,
                     data = bike, family = poisson)
interaction_test <- anova(bike_additive, bike_glm, test = "Chisq")
interaction_test

Analysis of Deviance Table

Model 1: bikers ~ hr + working + mnth + temp_c + weather
Model 2: bikers ~ hr * working + mnth + temp_c + weather
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1      8606     228046                          
2      8583      99377 23   128669 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\(G^2 = 128{,}669\) on \(k = 23\) degrees of freedom, against a \(\chi^2_{23}\) mean of 23

Fourth to drop if time is short (keep Drop-in-deviance test). The full model has the 23 hour-by-working-day products, the additive one omits them; the deviance drops by 128,669 on 23 degrees of freedom, against a chi-squared 23 mean of 23: the data are inconsistent with one hourly profile, up to a constant factor, shared by both kinds of day. This test assumes Var[Y_i | x_i] = mu_i, which the next slides question.

Overdispersion

Pearson residuals

Overdispersion: \(Var[Y_i \mid x_i] > \mu_i\)

\[\frac{y_i - \widehat\mu_i}{\sqrt{\widehat\mu_i}}\]

Pearson residual: standardized by the Poisson standard deviation

\[\widehat\phi = \frac{1}{n - (p+1)}\sum_{i=1}^n \frac{\left(y_i - \widehat\mu_i\right)^2}{\widehat\mu_i}\]

• \(\widehat\phi\) — estimates the dispersion parameter, fixed at 1 by the Poisson model
• Well above 1: overdispersion

The Poisson model has no variance parameter, so every standard error and deviance test above is computed under Var[Y_i | x_i] = mu_i. If that holds then E[(Y_i - mu_i)^2 / mu_i | x_i] = 1, so the average squared Pearson residual should be close to 1. The divisor n - (p + 1) accounts for the p + 1 estimated coefficients, as in sigma-hat squared = RSS / [n - (p + 1)].

Variance check

• All 20 bins above the dashed line
• \(\widehat\phi = 11.3\) (additive model: 25.6)

The figure bins the hours into 20 groups by fitted mean; all 20 bins have an average squared residual above the dashed line. The average squared Pearson residual is phi-hat = 11.3: the hourly counts have about 11 times the variance around their fitted means that a Poisson distribution allows. The additive model’s phi-hat is 25.6, so the hour-by-working-day features explained part of the excess variation, but phi-hat remains far above 1.

Effects of overdispersion

• \(\widehat\beta\), fitted means, rate ratios: still sensible if the mean model is right
• Standard errors too small: \(z\) too large, intervals too narrow
• \(G^2\) too large against \(\chi^2_k\): \(G^2/\widehat\phi \approx \chi^2_k\)

\[G^2 / \widehat\phi = 11{,}380\]

495 times its 23 degrees of freedom: the interaction test’s conclusion survives

The score equations involve y_i only through y_i - mu_i, never its variance, so if the mean model is right beta-hat and the fitted means stay sensible; the standard errors, computed as if the variance were mu_i, are too small. If Var[Y_i | x_i] = phi mu_i then G squared is approximately phi times chi-squared k, so G squared over phi-hat is approximately chi-squared k, and for the bike counts G squared over phi-hat is 11,380, still 495 times its 23 degrees of freedom. Overdispersion has two sources, as in the problems in logistic regression lecture: explanatory variables missing from the model that vary among observations with the same x_i (plausibly, conditions shared by the hours of one day), and events within one count that are not independent, such as riders who rent together; adding features addresses only the first.

Conclusion

• Generalized linear model
• Poisson regression
• Rate ratio
• Log-count linear regression
• Score equations
• Wald tests
• Drop-in-deviance test
• Overdispersion

Optional; the lecture may end on Effects of overdispersion with the next lecture’s preview said aloud: the resampling unit begins on Wednesday 2026-10-14 with cross-validation, which estimates test error by repeatedly holding out part of the data and uses those estimates to choose among models and flexibility settings (ISLR2 Chapter 5).