Generalised linear models

Why this matters

Ordinary regression assumes a continuous response with constant, roughly normal noise. Many responses are nothing like that. Did the assay pass or fail? Did the customer churn? How many defects are on this board, how many colonies on this plate, how many calls arrived in this hour? Forcing these into OLS gives predictions outside the possible range (probabilities above 1, negative counts) and wrong standard errors. Generalised linear models (GLMs) keep the familiar linear predictorThe weighted sum b0 + b1*x1 + ... + bp*xp that a GLM shares with ordinary regression. but connect it to the response in a way that respects what the response is.

Three ingredients

Every GLM has:

1. A distribution for y from the exponential family: normal, binomial (yes/no or successes out of n), Poisson (counts), gamma (positive, skewed amounts), and others.
2. A linear predictor, η = b0 + b1*x1 + ... + bp*xp, exactly as in multiple regression.
3. A link functionThe function g that connects the mean of the response to the linear predictor, such as logit for yes/no data or log for counts. g that connects the mean to the linear predictor: g(mean of y) = η.

The link guarantees the mean stays in range: a logistic model can never predict a probability of 1.3, and a Poisson model can never predict −2 defects.

GLMs are fitted by maximum likelihood, usually with an algorithm called iteratively reweighted least squares (IRLS), which is Newton's method in disguise: repeatedly fit a weighted regression, update the weights, and stop when the coefficients settle. You will write it in the activity.

Logistic regression: reading the coefficients

For a yes/no outcome, the model is linear in the log-odds:

log( p / (1 - p) ) = b0 + b1*x1 + ...

So exp(b1) is an odds ratioexp(b) in a logistic model: the factor by which the odds of yes multiply for a one-unit increase in the predictor.: the multiplicative change in the odds of "yes" for a one-unit increase in x1, holding the rest fixed. If a logistic model for assay failure has b_humidity = 0.12 per percentage point of humidity, then exp(0.12) ≈ 1.13: each extra point multiplies the odds of failure by about 1.13 (13% higher odds). Odds ratios are not risk ratios; when the outcome is rare they are close, but when it is common they are not.

Poisson regression: reading the coefficients

For counts, the model is linear in the log of the mean:

log(μ) = b0 + b1*x1 + ...      so      μ = exp(b0) * exp(b1)^x1 * ...

exp(b1) is a rate ratioexp(b) in a Poisson model: the factor by which the expected count multiplies for a one-unit increase in the predictor.. In the activity's solder example the slope is about 0.8 per unit of temperature offsetA term with its coefficient fixed at 1, usually log(exposure), so a count model describes a rate per unit of exposure., and exp(0.8) ≈ 2.2: moving one unit hotter roughly doubles the expected defects per board.

When counts come from different amounts of exposure (boards of different area, observation windows of different length), add log(exposure) as an offset, a term with its coefficient fixed at 1, so the model describes a rate per unit of exposure.

Worked example with statsmodels

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

# Logistic: pass/fail by humidity and operator
logit = smf.glm("failed ~ humidity + C(operator)", data=assays,
                family=sm.families.Binomial()).fit()
print(np.exp(logit.params))          # odds ratios
print(np.exp(logit.conf_int()))      # their 95% intervals

# Poisson with an exposure offset: defects per square cm of board
pois = smf.glm("defects ~ offset_temp + C(line)", data=boards,
               family=sm.families.Poisson(),
               offset=np.log(boards["area_cm2"])).fit()
print(np.exp(pois.params))           # rate ratios
print(pois.pearson_chi2 / pois.df_resid)   # dispersion check

Exponentiate both the coefficients and the interval ends; never exponentiate a standard error.

Overdispersion

The Poisson distribution has variance equal to its mean. Real counts are often more variable (some boards come from a bad reel, some plates were contaminated). The quick check is the dispersion ratioPearson chi-square divided by residual degrees of freedom; well above 1 signals overdispersion and too-small standard errors.: Pearson chi-square divided by residual degrees of freedom. Near 1 is fine. Well above 1 (say 1.5 or more) means the Poisson standard errors are too small and p-values too optimistic. Remedies:

• Quasi-Poisson (scale the standard errors by the dispersion), or
• Negative binomial regression (sm.families.NegativeBinomial()), which adds a parameter for extra variance, or
• a mixed model with a random effect for the clustering source (module 3).

Logistic models on grouped data (successes out of n) can be overdispersed in the same way.

Model checking for GLMs

Ordinary residuals are not very useful because the variance changes with the mean. Use deviance or Pearson residuals instead (model.resid_deviance, model.resid_pearson) and plot them against fitted values and each predictor, as in the previous lesson. Compare nested models with a likelihood-ratio test (difference in deviance, chi-square with the difference in parameters) and non-nested ones with AIC.

Common mistakes

• Reading logistic coefficients as changes in probability. They are changes in log-odds. To talk about probabilities, compute predicted probabilities at meaningful predictor values.
• Ignoring exposure in count models, so a line that simply ran longer looks worse.
• Skipping the dispersion check and reporting tiny p-values from an overdispersed Poisson fit.
• Perfect separation in logistic regression: if a predictor splits yes and no perfectly, the estimates run off to infinity. Software warns about it; consider penalised (Firth) regression or collect more varied data.

Quick recap

• A GLM = a distribution + a linear predictor + a link function.
• Logistic (binomial, logit): exp(b) is an odds ratio. Poisson (log): exp(b) is a rate ratio.
• Use offsets for exposure; check overdispersion; use deviance residuals.

Try it

In the Do tab you will turn a worked logistic IRLS loop into a Poisson regression, fit it to 500 simulated circuit boards and compute the rate ratio.