Regression Analysis and Generalized Linear Models

Regression analysis is the study of how the conditional mean (and sometimes other features) of a response $Y$ depends on covariates $\mathbf{x}$. The classical linear model assumes $E[Y\mid \mathbf{x}]={\mathbf{x}}^{\top }\mathbf{\beta }$ with homoscedastic Gaussian errors; the Generalized Linear Model relaxes both the distributional and the linearity-of-mean assumptions by positing that $Y$ follows an exponential-family distribution and that a monotone link function $g$ satisfies $g(\mu )=\eta ={\mathbf{x}}^{\top }\mathbf{\beta }$. This single framework unifies linear regression, logistic regression, Poisson regression, gamma regression, and more under one estimation algorithm and one inferential theory. A surprising fact: although these models look very different (continuous vs. binary vs. count data), their maximum likelihood estimation is identical in form — a weighted least squares iteration — because the exponential family endows them with a common geometric structure. GLMs are the workhorse of applied statistics, biostatistics, econometrics, and modern machine learning's "linear-then-link" layer.

The Core Insight

Linear regression is "least squares applied to a Gaussian." Replace the Gaussian by any exponential-family distribution and connect its mean to a linear predictor through a smooth invertible link. The exponential family's algebra guarantees that (i) the mean and variance are determined by a single natural parameter, (ii) the log-likelihood is concave in the natural parameter, and (iii) maximum likelihood estimation becomes Newton/Fisher-scoring that is exactly weighted least squares. Thus heterogeneous data types are tamed by a single, elegant, convex machine.

Formal Statement

The Linear Model

Given $n$ observations $(y_i,{\mathbf{x}}_i)$ with ${\mathbf{x}}_i\in {\mathbb{R}}^p$, the ordinary linear regression model is

equivalently in matrix form $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{\epsilon }$, where $\mathbf{X}\in {\mathbb{R}}^{n\times p}$ is the design matrix, $\mathbf{\beta }\in {\mathbb{R}}^p$ the coefficient vector, ${\sigma }^2$ the error variance, and $\mathbf{\epsilon }\sim \mathcal{N}(0,{\sigma }^2{\mathbf{I}}_n)$.

The Exponential Family

A distribution belongs to the (one-parameter) exponential dispersion family if its density/mass can be written

where $\theta$ is the natural (canonical) parameter, $\varphi >0$ the dispersion parameter, $b(\cdot )$ the cumulant function (a smooth convex function), $a(\varphi )$ typically $\varphi /w$ for known prior weight $w$, and $c(y,\varphi )$ a normalizing term independent of $\theta$.

Mean and Variance from the Cumulant Function

For any member of the exponential dispersion family [2],

where $V(\mu )=b^{\prime \prime }(\theta (\mu ))$ is the variance function, expressing the variance as a function of the mean.

Generalized Linear Model

A GLM consists of three components:

1. Random component: $Y_i$ independent with density [2] and means ${\mu }_i$.
2. Systematic component: a linear predictor ${\eta }_i={\mathbf{x}}_i^{\top }\mathbf{\beta }$.
3. Link function: a monotone, twice-differentiable $g$ with

The canonical link is $g=(b^{\prime }{)}^{-1}$, which makes ${\theta }_i={\eta }_i$ directly. For Gaussian, $g(\mu )=\mu$ (identity); for Bernoulli, $g(\mu )=\log \frac{\mu }{1-\mu }$ (logit); for Poisson, $g(\mu )=\log \mu$.

Score Equations and Fisher Information

The log-likelihood of a GLM yields score (gradient) equations

and the expected (Fisher) information matrix is

Why It Holds / Derivation

Maximum Likelihood is Iteratively Reweighted Least Squares

We derive Fisher scoring for [5] and show it equals weighted least squares.

1. Log-likelihood. From [2], the contribution of observation $i$ is

1. Score per observation. Differentiating with the chain rule and using ${\mu }_i=b^{\prime }({\theta }_i)$, $V({\mu }_i)=b^{\prime \prime }({\theta }_i)=\partial {\mu }_i/\partial {\theta }_i$,

1. Newton step with expected Hessian. Fisher scoring replaces the observed Hessian by $-\mathcal{I}$. The update is

1. Working response. Define the adjusted dependent variable

1. IRLS form. Substituting [6] and [10] into [9] gives

John Nelder & Robert Wedderburn | 1924–2010 / 1947–1975

John Nelder, a British statistician at Rothamsted Experimental Station, and Robert Wedderburn, his young collaborator, published "Generalized Linear Models" in the Journal of the Royal Statistical Society in 1972. Their insight was to recognize that linear regression, logistic regression, log-linear models, and probit analysis — previously taught as separate techniques — were all instances of one unifying structure: an exponential-family response and a link function, fitted by a single algorithm. Wedderburn further introduced quasi-likelihood (1974), freeing estimation from full distributional assumptions by requiring only a mean–variance relationship. Tragically, Wedderburn died at age 28 from anaphylaxis; Nelder went on to develop the GLIM software that made GLMs computationally accessible and championed the framework for decades.

Building Intuition

Picture ordinary least squares as fitting a flat plane to a cloud of points, where every point pulls equally hard (constant variance). A GLM says: not all points should pull equally — points where the response is intrinsically more variable should pull less. The weights $W_{ii}$ encode "how much to trust this observation," and they depend on the current fitted mean. The link function is a change of scale: for counts you fit on a log scale (so the prediction can never go negative); for probabilities you fit on a logit scale (so predictions stay in $[0,1]$). IRLS then alternates: "given current beliefs about the means, recompute trust weights and a linearized target; refit by weighted least squares; repeat." It is least squares wearing different costumes for different data types.

GLM structure: linear predictor, link, and exponential-family response

Worked Examples

Logistic Regression: Probability of Defect

1. Initialize ${\mathbf{\beta }}^{(0)}=(-2.5,0.1{)}^{\top }$. Compute ${\eta }_i={\beta }_0+{\beta }_1{x}_i=(-1.5,-0.5,0.5,1.5)$.
2. Inverse link ${\mu }_i=1/(1+e^{-{\eta }_i})=(0.182,0.378,0.622,0.818)$.
3. Weights $W_{ii}={\mu }_i(1-{\mu }_i)=(0.149,0.235,0.235,0.149)$ (since for canonical link $\partial \mu /\partial \eta =V(\mu )$, so $W=V$).
4. Working response $z_i={\eta }_i+(y_i-{\mu }_i)/W_{ii}$:

1. Solve ${\mathbf{\beta }}^{(1)}=({\mathbf{X}}^{\top }\mathbf{WX}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{Wz}$. By the data's symmetry the estimate centers at $x=25$; one iteration moves toward ${\hat{\beta }}_1\approx 0.13$, ${\hat{\beta }}_0\approx -3.3$, and IRLS converges within ~6 iterations.

Poisson Regression: Insurance Claim Counts

1. Rate ${\lambda }_i={\mu }_i/t_i=\exp ({\beta }_0+{\beta }_1{a}_i)$. Total exposure-by-group: low-risk $t=3$, claims $=1$; high-risk $t=5$, claims $=7$.
2. For Poisson with log link and a single binary covariate, the MLE matches observed rates: ${\hat{\lambda }}_{\text{low}}=1/3=0.333$, ${\hat{\lambda }}_{\text{high}}=7/5=1.40$.
3. Recover coefficients: ${\hat{\beta }}_0=\log (0.333)=-1.099$; ${\hat{\beta }}_0+{\hat{\beta }}_1=\log (1.40)=0.336\Rightarrow {\hat{\beta }}_1=1.435$.
4. Rate ratio $=e^{{\hat{\beta }}_1}=e^{1.435}=4.20$: high-risk policyholders file claims at 4.2× the rate.

Common Misconceptions

Linearity is in the parameters, not the covariates

"Linear" in linear and generalized linear model refers to linearity in $\mathbf{\beta }$, not in $\mathbf{x}$. You may include $x^2$, $\log x$, splines, or interactions and remain within the (G)LM framework. Conversely, a model that is nonlinear in the parameters (e.g., $\mu ={\beta }_0{e}^{{\beta }_{1}x}$) is not a GLM even though it has an exponential shape — unless you can recast it through a valid link.

Do not confuse the link scale with the response scale

Coefficients are interpreted on the link scale. In Poisson log-link regression, ${\beta }_1$ is a log rate ratio; in logistic regression, ${\beta }_1$ is a log odds ratio — neither is a change in probability or in count. Predictions must be back-transformed via $g^{-1}$, and because $g^{-1}$ is nonlinear, $E[g^{-1}(\eta )]\ne g^{-1}(E[\eta ])$ in general (Jensen's inequality).

Perfect separation breaks logistic MLE

If a covariate perfectly separates the two classes — every $y_i=1$ has $x_i>c$ and every $y_i=0$ has $x_i<c$ — the likelihood is monotone increasing as ${\hat{\beta }}_1\to \infty$. The MLE does not exist (it diverges), even though IRLS appears to "converge" to huge coefficients with enormous standard errors. The fitted probabilities approach 0/1 and Wald inference becomes meaningless. Remedies: Firth's penalized likelihood, Bayesian priors, or exact logistic regression. This shows that the existence of a finite MLE is not guaranteed in GLMs, unlike in full-rank Gaussian OLS.

Connections

GLMs sit at a crossroads of statistics. They generalize the Ordinary Least Squares linear model and are themselves special cases of broader frameworks. The Exponential Family supplies their probabilistic backbone, and Maximum Likelihood Estimation supplies their inferential engine; the IRLS algorithm is just Newton–Raphson Method / Fisher scoring in disguise. When the dispersion or correlation structure is richer, GLMs extend to Generalized Estimating Equations (for correlated data) and Generalized Linear Mixed Models (with random effects). Wedderburn's Quasi-Likelihood drops the full distributional assumption, keeping only the mean–variance link. Regularized variants connect to Ridge Regression and Lasso Regression, and the logistic GLM is precisely the single-layer building block of Neural Networks. Inference relies on the Wald Test, Likelihood Ratio Test, and the Deviance goodness-of-fit measure.

Real-World Applications

1. Epidemiology — logistic regression: modeling disease presence (case/control) against exposures; odds ratios from the Framingham Heart Study quantified cardiovascular risk factors.
2. Actuarial science — Poisson and Gamma GLMs: insurance pricing decomposes premiums into claim frequency (Poisson/negative binomial) and claim severity (Gamma), as standardized in the Tweedie GLM used across the industry.
3. Ecology — count models: species abundance vs. habitat covariates with Poisson/negative-binomial regression, central to species distribution modeling.
4. Clinical trials — survival and binary endpoints: complementary log-log links connect GLMs to discrete-time hazard models; logistic GLMs analyze responder rates adjusting for covariates.
5. Marketing & credit scoring — logistic GLMs: the FICO-style probability-of-default and churn/click-through models are logistic regressions with engineered features.
6. Genomics — log-linear and negative binomial GLMs: RNA-seq differential expression (edgeR, DESeq2) fits per-gene negative-binomial GLMs to read counts.
7. Econometrics — gravity models of trade: Poisson pseudo-maximum-likelihood (Santos Silva & Tenreyro, 2006) estimates international trade flows robustly under heteroscedasticity.

Historical Development

• 1805–1809 — Legendre publishes the method of least squares; Gauss claims earlier use and ties it to the normal error law, founding linear regression.
• 1821–1823 — Gauss states the Gauss–Markov optimality of least squares (BLUE under uncorrelated, equal-variance errors).
• 1886 — Francis Galton coins "regression toward mediocrity" studying heredity of stature, giving regression its name.
• 1922 — R. A. Fisher formalizes maximum likelihood and the concept of statistical information.
• 1934–1944 — Bliss develops probit analysis for bioassay; Berkson (1944) introduces the logit alternative.
• 1972 — Nelder & Wedderburn unify these strands in "Generalized Linear Models," defining the link/random/systematic triad and IRLS.
• 1974 — Wedderburn introduces quasi-likelihood, requiring only a mean–variance relationship.
• 1983/1989 — McCullagh & Nelder's monograph Generalized Linear Models becomes the definitive reference; GLIM software popularizes the methods.
• 1986 — Liang & Zeger introduce Generalized Estimating Equations, extending GLMs to longitudinal/correlated data.

Subtleties & Limits

Several issues remain active or delicate. Existence and uniqueness of the MLE can fail under separation or sparse data; penalized and Bayesian regularization are still being refined for high-dimensional GLMs ($p\gg n$), where classical asymptotics break down and even the bias of the MLE behaves non-classically (Sur & Candès, 2019). Dispersion modeling — when $\varphi$ itself depends on covariates — pushes beyond standard GLMs into double GLMs and distributional regression (GAMLSS). Misspecified mean–variance links undermine standard errors; sandwich/robust variance estimators help but lose efficiency. Causal interpretation of GLM coefficients requires assumptions (no confounding, correct functional form) far stronger than the model itself supplies, a recurring source of error in applied work. Finally, inference after model selection (choosing covariates then testing them) invalidates naive p-values; post-selection inference for GLMs is an ongoing research frontier.

See Also