Generalized Linear Models (GLM): A Comprehensive Technical Analysis

1. Introduction

1.1 The Problem

Ordinary least squares assumes a continuous response with roughly constant variance. Binary, count and skewed positive outcomes break those assumptions: the variance depends on the mean, and predictions can fall outside the possible range. Linear regression can still give useful average effects in some of these cases, but its predictions and its standard errors need care.

1.2 What This Paper Covers

• The three parts of a GLM: distribution, linear predictor and link
• Common families and links, and how to read their coefficients
• A reproducible comparison with linear regression
• Diagnostics, including overdispersion

2. Background

2.1 Structure of a GLM

A GLM (Nelder and Wedderburn, 1972) has a response from an exponential-family distribution with mean mu, a linear predictor eta = b₀ + b₁x₁ + ..., and a link function g with g(mu) = eta. The link keeps predictions in range: the logit keeps probabilities between 0 and 1, the log keeps means positive.

2.2 Common Families

3. Approach

3.1 Fitting

GLMs are fitted by maximum likelihood, usually through iteratively reweighted least squares, which statistical packages implement. Fit quality is reported through the deviance, a likelihood-based generalization of the residual sum of squares, and comparisons of nested models use differences in deviance.

3.2 Interpretation

On the link scale effects are additive. With a log link a coefficient multiplies the expected count or mean (a rate ratio), with a logit link it multiplies the odds (an odds ratio), and with an identity link it adds to the mean. For communication, translate to predicted probabilities or counts at realistic predictor values, because odds ratios are easily misread as risk ratios.

3.3 Checks

• Dispersion: for Poisson and binomial counts, the Pearson chi-square divided by residual degrees of freedom should be near 1. Much larger values mean overdispersion, which understates standard errors. For counts, use a negative binomial or quasi-Poisson model or robust standard errors; for binomial data, use a quasi-binomial model or robust standard errors, because a negative binomial model does not respect the binomial trial limit.
• Separation: in logistic regression, a predictor that perfectly separates outcomes makes coefficients diverge; check for it and consider penalized estimation.
• Residuals and influence: inspect deviance residuals and leverage for patterns and outliers.
• Held-out performance: compare predictions on data not used for fitting.

4. A Reproducible Comparison

The simulation draws 5,000 training and 5,000 test rows with two predictors. The first outcome is binary with a true logistic relationship; the second is a count with a true Poisson relationship using a log link. Each is fitted by linear regression and by the matching GLM, and judged on the test rows by log loss (binary) and mean Poisson deviance (counts), where lower is better. The last lines compute the dispersion check for Poisson counts and for counts with extra variation.

import numpy as np
import statsmodels.api as sm
from scipy import special

rng = np.random.default_rng(12)
n = 5000
def design(m):
    x = rng.normal(size=(m, 2))
    return x, sm.add_constant(x)

# 1. Binary outcome (converted or not), true model is logistic
x_tr, X_tr = design(n); x_te, X_te = design(n)
lin = lambda x: -1.5 + 1.2 * x[:, 0] - 0.8 * x[:, 1]
y_tr = rng.binomial(1, special.expit(lin(x_tr))); y_te = rng.binomial(1, special.expit(lin(x_te)))
ols = sm.OLS(y_tr, X_tr).fit(); logit = sm.Logit(y_tr, X_tr).fit(disp=0)
p_ols = ols.predict(X_te); p_logit = logit.predict(X_te)
logloss = lambda y, p: -np.mean(y * np.log(p) + (1 - y) * np.log(1 - p))
print("binary outcome, 5,000 test rows")
print(f"  linear regression: {np.mean((p_ols < 0) | (p_ols > 1)):.1%} of predicted probabilities fall outside [0, 1]; "
      f"log loss {logloss(y_te, np.clip(p_ols, 0.001, 0.999)):.3f} after clipping")
print(f"  logistic regression: all probabilities inside [0, 1]; log loss {logloss(y_te, p_logit):.3f}")

# 2. Count outcome (orders per customer), true model is Poisson with a log link
mu = lambda x: np.exp(0.2 + 0.6 * x[:, 0] - 0.4 * x[:, 1])
c_tr = rng.poisson(mu(x_tr)); c_te = rng.poisson(mu(x_te))
ols_c = sm.OLS(c_tr, X_tr).fit(); pois = sm.GLM(c_tr, X_tr, family=sm.families.Poisson()).fit()
pred_ols, pred_pois = ols_c.predict(X_te), pois.predict(X_te)
def poisson_deviance(y, m):
    m = np.maximum(m, 0.01)
    return np.mean(2 * (np.where(y > 0, y * np.log(np.where(y > 0, y, 1) / m), 0) - (y - m)))
print("count outcome, 5,000 test rows")
print(f"  linear regression: {np.mean(pred_ols < 0):.1%} of predicted counts are negative; mean deviance {poisson_deviance(c_te, pred_ols):.3f}")
print(f"  Poisson regression: mean deviance {poisson_deviance(c_te, pred_pois):.3f}")

# 3. A dispersion check: extra variation beyond Poisson shows up as Pearson chi-square / df well above 1
nb = rng.negative_binomial(2, 2 / (2 + mu(x_tr)))
for name, counts in (("Poisson counts", c_tr), ("overdispersed counts", nb)):
    fit = sm.GLM(counts, X_tr, family=sm.families.Poisson()).fit()
    print(f"  {name}: Pearson chi-square / df = {fit.pearson_chi2 / fit.df_resid:.2f}")

Output:

binary outcome, 5,000 test rows
  linear regression: 10.9% of predicted probabilities fall outside [0, 1]; log loss 0.458 after clipping
  logistic regression: all probabilities inside [0, 1]; log loss 0.443
count outcome, 5,000 test rows
  linear regression: 8.0% of predicted counts are negative; mean deviance 1.357
  Poisson regression: mean deviance 1.104
  Poisson counts: Pearson chi-square / df = 0.98
  overdispersed counts: Pearson chi-square / df = 1.83

The GLMs won on both measures. About 11% of the linear model's probabilities and 8% of its counts were impossible, so clipping or truncation was needed before scoring it, and its prediction error was higher even after clipping. These data were generated from exactly the GLM form, so the gap is the best case for the GLM; on real data the size of the gain depends on the conditional mean, the family and the link, and has to be measured on held-out observations. The dispersion check separated the two kinds of count data, 0.98 for Poisson counts and 1.83 for overdispersed ones.

5. Practical Guidance

6. Conclusion

Generalized linear models let regression fit outcomes that are binary, counts or skewed amounts while keeping predictions in range. In the simulation here, the matching GLM predicted held-out data better than linear regression and never produced impossible values. Choose the family from the outcome, check dispersion and separation, and compare against simpler models before relying on the result.