Regression Analysis and Generalized Linear Models

⚡ Key Takeaways

Regression analysis models the conditional distribution of a response variable given predictors; Generalized Linear Models (GLMs) extend ordinary linear regression to response distributions in the exponential family via a link function relating the mean to a linear predictor. Key result: maximum likelihood estimation in a GLM reduces to Iteratively Reweighted Least Squares (IRLS), which solves the score equations efficiently. Applies wherever responses are continuous, binary, count, or otherwise non-Gaussian. Foundational work by Nelder & Wedderburn (1972); roots in Legendre, Gauss, and Galton.

💡 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.

DefinitionThe Linear Model

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

1

\quad y_i = \mathbf{x}_i^\top\boldsymbol\beta + \varepsilon_i,\qquad \varepsilon_i \overset{\text{iid}}{\sim}\mathcal N(0,\sigma^2),

equivalently in matrix form $\mathbf{y} = \mathbf{X}\boldsymbol\beta + \boldsymbol\varepsilon$ , where $\mathbf{X}\in\mathbb{R}^{n\times p}$ is the design matrix, $\boldsymbol\beta\in\mathbb{R}^p$ the coefficient vector, $\sigma^2$ the error variance, and $\boldsymbol\varepsilon\sim\mathcal N(\mathbf 0,\sigma^2\mathbf I_n)$ . Derivation hint: stack the scalar equations row-wise; the independence of the $\varepsilon_i$ makes the covariance diagonal.

DefinitionThe Exponential Family

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

2

\quad f(y;\theta,\phi)=\exp\!\left\{\frac{y\theta-b(\theta)}{a(\phi)}+c(y,\phi)\right\},

where $\theta$ is the natural (canonical) parameter, $\phi>0$ the dispersion parameter, $b(\cdot)$ the cumulant function (a smooth convex function), $a(\phi)$ typically $\phi/w$ for known prior weight $w$ , and $c(y,\phi)$ a normalizing term independent of $\theta$ . Derivation hint: write Gaussian/Bernoulli/Poisson densities in exponential form and read off $\theta$ and $b(\theta)$ .

TheoremMean and Variance from the Cumulant Function

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

3

\quad \mu \equiv E[Y]=b'(\theta),\qquad \operatorname{Var}(Y)=a(\phi)\,b''(\theta)=a(\phi)\,V(\mu),

where $V(\mu)=b''(\theta(\mu))$ is the variance function, expressing the variance as a function of the mean. Derivation hint: differentiate the identity $\int f\,dy=1$ twice with respect to $\theta$ and use $E[\partial_\theta \log f]=0$ and $\operatorname{Var}(\partial_\theta\log f)=-E[\partial_\theta^2\log f]$ .

DefinitionGeneralized 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\boldsymbol\beta$ . 3. Link function: a monotone, twice-differentiable $g$ with

4

\quad g(\mu_i)=\eta_i=\mathbf{x}_i^\top\boldsymbol\beta,\qquad \mu_i=g^{-1}(\eta_i).

The canonical link is $g=( b')^{-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$ . Derivation hint: set $\theta=\eta$ and invert [3] $\mu=b'(\theta)$ to obtain $g=(b')^{-1}$ .

TheoremScore Equations and Fisher Information

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

5

\quad U(\boldsymbol\beta)=\sum_{i=1}^n \frac{(y_i-\mu_i)}{a(\phi)\,V(\mu_i)}\,\frac{\partial\mu_i}{\partial\eta_i}\,\mathbf{x}_i = \mathbf 0,

and the expected (Fisher) information matrix is

6

\quad \mathcal I(\boldsymbol\beta)=\mathbf X^\top\mathbf W\mathbf X,\qquad W_{ii}=\frac{1}{a(\phi)\,V(\mu_i)}\left(\frac{\partial\mu_i}{\partial\eta_i}\right)^2.

Derivation hint: apply the chain rule $\partial\ell/\partial\beta_j=\sum_i (\partial\ell_i/\partial\theta_i)(\partial\theta_i/\partial\mu_i)(\partial\mu_i/\partial\eta_i)(\partial\eta_i/\partial\beta_j)$ and use $\partial\ell_i/\partial\theta_i=(y_i-\mu_i)/a(\phi)$ .

ProofMaximum 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

7

\quad \ell_i = \frac{y_i\theta_i-b(\theta_i)}{a(\phi)}+c(y_i,\phi).

Derivation hint: take the logarithm of the exponential-family density. 2. Score per observation. Differentiating with the chain rule and using $\mu_i=b'(\theta_i)$ , $V(\mu_i)=b''(\theta_i)=\partial\mu_i/\partial\theta_i$ ,

8

\quad \frac{\partial\ell_i}{\partial\beta_j}=\frac{y_i-\mu_i}{a(\phi)V(\mu_i)}\,\frac{\partial\mu_i}{\partial\eta_i}\,x_{ij}.

Derivation hint: $\partial\theta_i/\partial\mu_i = 1/V(\mu_i)$ since $\mu=b'(\theta)$ . 3. Newton step with expected Hessian. Fisher scoring replaces the observed Hessian by $-\mathcal I$ . The update is

9

\quad \boldsymbol\beta^{(t+1)} = \boldsymbol\beta^{(t)} + \mathcal I(\boldsymbol\beta^{(t)})^{-1}U(\boldsymbol\beta^{(t)}).

Derivation hint: Newton–Raphson on the score, swapping observed for expected information. 4. Working response. Define the adjusted dependent variable

10

\quad z_i = \eta_i + (y_i-\mu_i)\,\frac{\partial\eta_i}{\partial\mu_i}.

Derivation hint: linearize $g(y_i)$ about $\mu_i$ : $g(y_i)\approx g(\mu_i)+(y_i-\mu_i)g'(\mu_i)$ . 5. IRLS form. Substituting [6] and [10] into [9] gives

11

\quad \boldsymbol\beta^{(t+1)}=\left(\mathbf X^\top\mathbf W\mathbf X\right)^{-1}\mathbf X^\top\mathbf W\,\mathbf z,

which is exactly the weighted least squares solution regressing the working response $\mathbf z$ on $\mathbf X$ with weights $\mathbf W$ . Iterating until convergence yields the MLE $\hat{\boldsymbol\beta}$ . Derivation hint: collect $U=\mathbf X^\top\mathbf W \mathbf D^{-1}(\mathbf y-\boldsymbol\mu)$ with $\mathbf D=\mathrm{diag}(\partial\mu_i/\partial\eta_i)$ , then absorb terms so $\mathcal I^{-1}U + \boldsymbol\beta = (\mathbf X^\top\mathbf W\mathbf X)^{-1}\mathbf X^\top\mathbf W\mathbf z$ .

∎

mathematician · 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.

🧠 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.

Figure

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

ExampleLogistic Regression: Probability of Defect

Given. An engineer tests $n=4$ components at temperatures $x=(10,20,30,40)$ °C with binary outcomes "defect" $y=(0,0,1,1)$ . We fit $\operatorname{logit}(\mu_i)=\beta_0+\beta_1 x_i$ , where $\mu_i=P(Y_i=1)$ . Variance function $V(\mu)=\mu(1-\mu)$ , canonical logit link, dispersion $\phi=1$ .

Steps. 1. Initialize $\boldsymbol\beta^{(0)}=(\,-2.5,\;0.1\,)^\top$ . Compute $\eta_i=\beta_0+\beta_1x_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}$ : $z=(-1.5-0.182/0.149,\; -0.5-0.378/0.235,\;0.5+0.378/0.235,\;1.5+0.182/0.149)=(-2.72,-2.11,2.11,2.72)$ . 5. Solve $\boldsymbol\beta^{(1)}=(\mathbf X^\top\mathbf W\mathbf X)^{-1}\mathbf X^\top\mathbf W\mathbf z$ . By the data's symmetry the estimate centers at $x=25$ ; one iteration moves toward $\hat\beta_1\approx0.13$ , $\hat\beta_0\approx-3.3$ , and IRLS converges within ~6 iterations.

Answer.

\boxed{\hat\beta_0\approx -3.3,\quad \hat\beta_1\approx 0.13\ \text{per °C};\quad \hat P(\text{defect}\mid 25^\circ)=0.5}

ExamplePoisson Regression: Insurance Claim Counts

Given. Policyholders' claim counts $y$ with exposure (years) $t$ . Model $\log\mu_i=\log t_i+\beta_0+\beta_1 a_i$ where $a_i$ is a risk indicator and $\log t_i$ is an offset. Data: $(a,t,y)=(0,2,1),(0,1,0),(1,3,4),(1,2,3)$ . Variance function $V(\mu)=\mu$ , log link, $\phi=1$ .

Steps. 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.

Answer.

\boxed{\hat\beta_0=-1.10,\quad \hat\beta_1=1.44,\quad \text{rate ratio}=4.20\ \text{(claims/year)}}

⚠️ Linearity is in the parameters, not the covariates

"Linear" in linear and generalized linear model refers to linearity in $\boldsymbol\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)]\neq g^{-1}(E[\eta])$ in general (Jensen's inequality).

CounterexamplePerfect 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.

connection

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.

🔬 Where this appears

⏳ History

- 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.

open-problem

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 $\phi$ 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.

Regression Analysis and Generalized Linear Models

Formal Statement

Why It Holds / Derivation

Building Intuition

Worked Examples

Common Misconceptions

Connections

Real-World Applications

Historical Development

Subtleties & Limits

See Also