Generalized linear models · Probability & Statistics

Linear regression assumes the response is continuous and roughly Normal, with constant variance. Plenty of data isn't like that.

Binary: did the customer buy? Did the patient recover?
Counts: how many accidents this month?
Positive and skewed: insurance claim sizes.

Fitting a straight line to a 0/1 outcome produces predicted "probabilities" below 0 and above 1, and errors that are obviously not Normal or constant- variance. Generalized linear models fix this by keeping the linear predictor and changing two things around it.

The three components

A GLM has:

1. A random component — the distribution of Y, from the exponential family (Normal, Binomial, Poisson, Gamma, …).

2. A linear predictor — unchanged from ordinary regression:

$\eta = \mathbf{X}\boldsymbol\beta$

3. A link function g connecting the mean to the linear predictor:

$g(\mu) = \eta, \qquad \mu = E[Y \mid \mathbf x]$

The link is what makes it work. It maps the constrained mean — a probability in (0,1), a rate in (0,\infty) — onto the whole real line, where a linear predictor is free to roam.

Model	Distribution	Link	\mu lives in
Linear	Normal	identity: \mu	\mathbb{R}
Logistic	Binomial	logit: \log\frac{\mu}{1-\mu}	(0,1)
Poisson	Poisson	log: \log\mu	(0,\infty)
Gamma	Gamma	log or inverse	(0,\infty)

Ordinary linear regression is the special case: Normal distribution, identity link.

Note the variance is no longer a free constant — it's determined by the distribution. For a Binomial it's \mu(1-\mu); for a Poisson it's \mu. The heteroscedasticity is built into the model rather than being a violation of it.

Logistic regression

For a binary outcome, model the log-odds:

$\log\frac{p}{1-p} = \beta_0 + \beta_1x_1 + \cdots + \beta_px_p$

Invert to get the probability:

$p = \frac{1}{1 + e^{-\eta}} = \frac{e^\eta}{1 + e^\eta}$

This is the logistic (sigmoid) function, and it maps any real \eta into (0,1) — so predictions can never be impossible.

Interpreting coefficients

This is where people go wrong. \beta_j is not a change in probability.

$e^{\beta_j} = \text{odds ratio for a one-unit increase in } x_j$

So \beta_1 = 0.7 means e^{0.7} \approx 2.01 — the odds roughly double per unit.

Odds are not probability. Odds of 2 means p = 2/3, not 2. And the change in probability for a one-unit move depends on where you start: near p = 0.5 the effect is largest; near 0 or 1 it's tiny. A constant effect on the log-odds scale is a varying effect on the probability scale.

That's not a defect — it's exactly the S-shape you want, and it's why the log- odds scale is where the model is linear (as the delta method showed in §3).

Fitting

There's no closed form. The log-likelihood

$\ell(\boldsymbol\beta) = \sum_i \big[y_i\eta_i - \log(1 + e^{\eta_i})\big]$

is concave, so it has a unique maximum and is fitted by iteratively reweighted least squares (IRLS) or gradient methods — reliably, in practice.

Separation. If some predictor perfectly separates the classes, the MLE doesn't exist: the likelihood keeps increasing as the coefficient goes to infinity. Software reports huge coefficients and enormous standard errors. The fix is penalisation (§6, next lesson) or Firth's correction.

Poisson regression

For count data:

$\log\mu = \beta_0 + \beta_1x_1 + \cdots \implies \mu = e^{\beta_0 + \beta_1x_1 + \cdots}$

The log link guarantees \mu > 0, and it makes effects multiplicative: e^{\beta_j} is the rate ratio per unit of x_j.

Offsets. If units are observed for different exposures (person-years, area, time), model the rate rather than the count by adding \log(\text{exposure}) as a term with coefficient fixed at 1:

$\log\mu = \log(\text{exposure}) + \mathbf x^\top\boldsymbol\beta$

Overdispersion. The Poisson insists \operatorname{Var}(Y) = \mu. Real count data is usually more variable than that, which makes standard errors too small and p-values too optimistic. Check the ratio of Pearson \chi^2 to its degrees of freedom; if it's well above 1, switch to a negative binomial or quasi-Poisson model.

Inference

GLMs are fitted by maximum likelihood, so §4's machinery applies directly:

Standard errors from the inverse Fisher information.
Deviance = -2(\ell_{\text{model}} - \ell_{\text{saturated}}) plays the role of the residual sum of squares.
Nested models are compared by the likelihood ratio test with Wilks' \chi^2 (§5).
AIC = -2\ell + 2k for non-nested comparisons.

Worked example

Predict passing an exam from hours studied:

Hours	1	2	3	4	5	6	7	8
Passed	0	0	0	1	0	1	1	1

Fitting a logistic regression gives

$\log\frac{\hat p}{1-\hat p} = -5.77 + 1.28 \times \text{hours}$

Interpreting the slope: e^{1.28} \approx 3.60, so each extra hour of study multiplies the odds of passing by about 3.6.

Predicting at 5 hours:

$\eta = -5.77 + 1.28(5) = 0.64, \qquad \hat p = \frac{1}{1 + e^{-0.64}} \approx 0.655$

The probability effect is not constant. Going from 4 to 5 hours moves $\hat p$ from 0.345 to 0.655 — a 31 point jump. Going from 8 to 9 hours moves it from 0.989 to 0.997 — under 1 point. Same coefficient, wildly different practical effect, because the sigmoid flattens at the ends.

The 50% threshold sits where \eta = 0, i.e. at $-\beta_0/\beta_1 = 5.77/ 1.28 \approx 4.5$ hours. That's often the most interpretable single number a logistic model gives you.

Doing it in Python

Logistic regression fitted by maximising the likelihood directly:

import numpy as np
from scipy.optimize import minimize

hours = np.array([1, 2, 3, 4, 5, 6, 7, 8], dtype=float)
passed = np.array([0, 0, 0, 1, 0, 1, 1, 1], dtype=float)
X = np.column_stack([np.ones(len(hours)), hours])

def neg_loglik(beta):
    eta = X @ beta
    # log(1 + exp(eta)), computed stably for large |eta|
    return -np.sum(passed * eta - np.logaddexp(0, eta))

res = minimize(neg_loglik, x0=[0.0, 0.0])
b0, b1 = res.x
print(f"intercept {b0:.4f}   slope {b1:.4f}")
print(f"odds ratio per hour: e^{b1:.4f} = {np.exp(b1):.4f}")
print(f"50% threshold at {-b0/b1:.4f} hours")

sigmoid = lambda z: 1 / (1 + np.exp(-z))
print(f"\n{'hours':>7} {'p_hat':>8}")
for h in (1, 4, 5, 8, 9):
    print(f"{h:>7} {sigmoid(b0 + b1*h):>8.4f}")

Why the effect on probability isn't constant:

import numpy as np

b0, b1 = -5.7703, 1.2823
sigmoid = lambda z: 1 / (1 + np.exp(-z))

print(f"{'from -> to':>14} {'p change':>12} {'odds ratio':>12}")
for lo in (1, 4, 5, 8):
    p_lo, p_hi = sigmoid(b0 + b1*lo), sigmoid(b0 + b1*(lo+1))
    odds_lo = p_lo / (1 - p_lo)
    odds_hi = p_hi / (1 - p_hi)
    print(f"{f'{lo} -> {lo+1}':>14} {p_hi - p_lo:>12.4f} {odds_hi/odds_lo:>12.4f}")

print("\nThe odds ratio is CONSTANT (that's what the model assumes).")
print("The probability change is not -- it peaks near p = 0.5.")

Poisson regression, and detecting overdispersion:

import numpy as np
from scipy.optimize import minimize

rng = np.random.default_rng(0)
n = 300
x = rng.uniform(0, 3, n)
y = rng.poisson(np.exp(0.5 + 0.8 * x))          # genuinely Poisson
X = np.column_stack([np.ones(n), x])

def neg_loglik(beta):
    eta = X @ beta
    return -np.sum(y * eta - np.exp(eta))       # dropping the log(y!) constant

res = minimize(neg_loglik, x0=[0.0, 0.0])
print("fitted:", res.x.round(4), " (true 0.5, 0.8)")
print("rate ratio per unit x:", round(float(np.exp(res.x[1])), 4))

mu = np.exp(X @ res.x)
pearson = ((y - mu) ** 2 / mu).sum()
print(f"\nPearson chi2 / df = {pearson / (n - 2):.4f}   (near 1 = no overdispersion)")

# Now genuinely overdispersed data
y2 = rng.negative_binomial(2, 2 / (2 + np.exp(0.5 + 0.8 * x)))
res2 = minimize(lambda b: -np.sum(y2 * (X @ b) - np.exp(X @ b)), x0=[0.0, 0.0])
mu2 = np.exp(X @ res2.x)
print(f"overdispersed data: Pearson chi2 / df = "
      f"{((y2 - mu2)**2 / mu2).sum() / (n - 2):.4f}   <- well above 1")

Why a linear model on binary data is wrong:

import numpy as np

hours = np.array([1, 2, 3, 4, 5, 6, 7, 8], dtype=float)
passed = np.array([0, 0, 0, 1, 0, 1, 1, 1], dtype=float)
X = np.column_stack([np.ones(len(hours)), hours])

lin, *_ = np.linalg.lstsq(X, passed, rcond=None)
print("linear fit on 0/1 data:", lin.round(4))
print(f"\n{'hours':>7} {'linear pred':>13} {'valid?':>8}")
for h in (0, 1, 5, 10, 12):
    p = lin[0] + lin[1] * h
    print(f"{h:>7} {p:>13.4f} {'yes' if 0 <= p <= 1 else 'NO':>8}")

print("\nOutside the observed range it predicts probabilities below 0 and above 1.")
print("The logistic link makes that structurally impossible.")

Your turn

1. Logistic regression gives \beta_1 = 1.5. Interpret it.

2. Why not just fit a linear model to a 0/1 outcome?

3. In a Poisson regression, \beta_1 = 0.4. What does that mean?

Solutions

1. e^{1.5} \approx 4.48: each one-unit increase in x_1 multiplies the odds of the outcome by about 4.5, holding other predictors fixed.

Careful phrasing matters here:

✓ "The odds are 4.5 times higher."
✗ "The probability is 4.5 times higher." (False except when p is tiny.)
✗ "The probability increases by 1.5." (Wrong scale entirely.)

To say anything about probability you must pick a starting point. From p = 0.1 (odds 0.111), the new odds are 0.498, so p \approx 0.332 — a rise of 23 points. From p = 0.5, the new odds are 4.48, so p \approx 0.817 — a rise of 32 points. From p = 0.9, it goes to 0.976 — under 8 points.

2. Several reasons, in order of seriousness:

Impossible predictions. A line is unbounded, so it will predict $\hat p < 0$ and \hat p > 1 outside the middle of the data.
Heteroscedasticity by construction. For a binary outcome \operatorname{Var}(Y) = p(1-p), which varies with p. Constant variance is violated automatically, so the usual standard errors are wrong.
Non-Normal errors. With y \in \{0,1\} the residuals take only two values for a given \hat y — nothing like Normal, which matters for small-sample inference.
Wrong functional form. The true relationship is almost always S-shaped: saturating at both ends. A line can't represent that.

The linear probability model does get used, mainly in economics, because the coefficients are directly interpretable as probability changes and it's robust for average effects near the middle. But it needs robust standard errors, and its predictions must not be extrapolated.

3. With a log link, e^{0.4} \approx 1.49: a one-unit increase in x_1 multiplies the expected count by about 1.49 — a 49% increase.

Effects are multiplicative, not additive. Going from an expected 10 events to 14.9 is the same coefficient as going from 100 to 149.

If the model includes an offset for exposure, the same number is a rate ratio — 49% more events per unit of exposure — which is usually the more meaningful statement.

Check yourself in code

Fit a logistic regression to the exam data by maximum likelihood and report the odds ratio and the 50% threshold.

Print exactly this:

intercept -5.7703
slope 1.2823
odds ratio 3.6049
50% threshold 4.5

Round every value to 4 decimal places. Use scipy.optimize.minimize starting from [0.0, 0.0], and np.logaddexp(0, eta) for numerical stability.

import numpy as np
from scipy.optimize import minimize

hours = np.array([1, 2, 3, 4, 5, 6, 7, 8], dtype=float)
passed = np.array([0, 0, 0, 1, 0, 1, 1, 1], dtype=float)
X = np.column_stack([np.ones(len(hours)), hours])

def neg_loglik(beta):
    eta = X @ beta
    return -np.sum(passed * eta - np.logaddexp(0, eta))

res = minimize(neg_loglik, x0=[0.0, 0.0])
b0, b1 = res.x
print("intercept", round(b0, 4))

# Print the slope, the odds ratio e^slope, and the hours at which
# the predicted probability crosses 0.5.

import numpy as np
from scipy.optimize import minimize

hours = np.array([1, 2, 3, 4, 5, 6, 7, 8], dtype=float)
passed = np.array([0, 0, 0, 1, 0, 1, 1, 1], dtype=float)
X = np.column_stack([np.ones(len(hours)), hours])

def neg_loglik(beta):
    eta = X @ beta
    return -np.sum(passed * eta - np.logaddexp(0, eta))

res = minimize(neg_loglik, x0=[0.0, 0.0])
b0, b1 = res.x
print("intercept", round(b0, 4))
print("slope", round(b1, 4))
print("odds ratio", round(float(np.exp(b1)), 4))
print("50% threshold", round(-b0 / b1, 4))

A GLM keeps the linear predictor \mathbf X\boldsymbol\beta and wraps it in a distribution and a link function, so the mean stays inside its natural range. Logistic regression models log-odds — coefficients are odds ratios, never probability changes. Poisson regression models log-rates, with multiplicative effects and a variance assumption worth checking.

Next: what to do when you have too many predictors, and least squares overfits.

43. Generalized linear models