Nonlinear Regression, Robust Regression & Generalized Linear Models (Logistic & Poisson)

Throughout Units 1 through 7, every regression model—no matter how curved in the predictor space $\mathbf{x}$ (such as polynomial regression, cubic B-splines, or Fourier series)—shared one defining structural property: the conditional expectation function $\mathbb{E}[Y_i \mid \mathbf{x}_i] = \mathbf{x}_i^T \boldsymbol{\beta}$ had partial derivatives with respect to the unknown parameters $\beta_j$ that depended **only on the known covariates $\mathbf{x}_i$ and never on the unknown parameters $\boldsymbol{\beta}$ themselves**. In chemical kinetics, enzyme biochemistry, pharmacokinetics, population ecology, and solid-state physics, mechanistic differential equations frequently generate response trajectories of the general form: $$Y_i = f(\mathbf{x}_i, \boldsymbol{\theta}) + \varepsilon_i, \quad i = 1, 2, \dots, n$$ where $\boldsymbol{\theta} = (\theta_1, \theta_2, \dots, \theta_p)^T \in \Theta \subset \mathbb{R}^p$ is a $p \times 1$ parameter vector and $f(\mathbf{x}_i, \boldsymbol{\theta})$ is a known smooth nonlinear expectation function. A regression function $f(\mathbf{x}, \boldsymbol{\theta})$ is Nonlinear in the Parameters whenever at least one partial derivative $\frac{\partial f(\mathbf{x}_i, \boldsymbol{\theta})}{\partial \theta_j}$ is a function of the unknown parameter vector $\boldsymbol{\theta}$. Statisticians classify nonlinear mechanistic curves into two fundamental categories: 1. **Intrinsically Linear (Transformably Linear) Models**: The deterministic function $f(\mathbf{x}, \boldsymbol{\theta})$ can be transformed into an exact linear form $\mathbf{x}^T \boldsymbol{\beta}$ by a monotonic transformation of the response or covariates. For example, in a multiplicative exponential growth model with **multiplicative log-normal errors** $Y_i = \theta_1 e^{\theta_2 x_i} \cdot u_i$, taking natural logarithms yields $\ln(Y_i) = \ln(\theta_1) + \theta_2 x_i + \ln(u_i)$, which is an exact simple linear regression with $\beta_0 = \ln(\theta_1)$ and $\varepsilon_i = \ln(u_i)$. 2. **Intrinsically Nonlinear Models**: Either the functional form $f(\mathbf{x}, \boldsymbol{\theta})$ cannot be algebraically linearized by any transformation (such as a two-compartment pharmacokinetic sum of exponentials $f(x, \boldsymbol{\theta}) = \theta_1 e^{-\theta_2 x} + \theta_3 e^{-\theta_4 x}$ or a three-parameter asymptotic exponential $f(x, \boldsymbol{\theta}) = \theta_1 - \theta_2 e^{-\theta_3 x}$), **or** the stochastic disturbance $\varepsilon_i$ enters **additively** on the original scale ($Y_i = f(x_i, \boldsymbol{\theta}) + \varepsilon_i$ with $\text{Var}(\varepsilon_i) = \sigma^2$). ### Why Naive Algebraic Linearization Fails Under Additive Homoscedastic Errors Consider the Michaelis–Menten enzyme kinetics model $Y_i = \frac{\theta_1 x_i}{\theta_2 + x_i} + \varepsilon_i$ where $\varepsilon_i \overset{\text{iid}}{\sim} \mathcal{N}(0, \sigma^2)$ represents constant-variance instrument measurement noise on reaction velocity $Y_i$. Biochemists historically inverted both sides to form the **Lineweaver–Burk plot**, regressing $1/Y_i$ on $1/x_i$. By the Delta Method (Section 7.5), if $\text{Var}(Y_i) = \sigma^2$ is constant, then $\text{Var}(1/Y_i) \approx \frac{\sigma^2}{\mu_i^4}$! Thus, at low substrate concentrations $x_i \to 0$ (where $\mu_i$ is small and $1/x_i$ is huge), the transformed error variance explodes as $\mu_i^{-4}$, turning the least-accurate measurements at the lowest concentrations into extreme high-leverage points that corrupt the OLS slope! Fitting the model directly on its natural scale via **Nonlinear Least Squares (NLS)** preserves homoscedasticity and achieves full asymptotic efficiency.

In the nonlinear regression model $Y_i = f(\mathbf{x}_i, \boldsymbol{\theta}) + \varepsilon_i$ with $\mathbb{E}[\boldsymbol{\varepsilon}] = \mathbf{0}_n$ and $\text{Cov}(\boldsymbol{\varepsilon}) = \sigma^2 \mathbf{I}_n$, the Nonlinear Least Squares (NLS) Estimator $\hat{\boldsymbol{\theta}}$ minimizes the residual sum of squares objective function $S(\boldsymbol{\theta})$:

where $\mathbf{f}(\boldsymbol{\theta}) = (f(\mathbf{x}_1, \boldsymbol{\theta}), \dots, f(\mathbf{x}_n, \boldsymbol{\theta}))^T \in \mathbb{R}^n$. Let $\mathbf{J}(\boldsymbol{\theta}) \in \mathbb{R}^{n \times p}$ denote the $n \times p$ local Jacobian (sensitivity) matrix of first partial derivatives of the expectation surface with respect to $\boldsymbol{\theta}$:

Setting the gradient of $S(\boldsymbol{\theta})$ equal to zero gives the $p$ Nonlinear Normal Equations:

Because both $\mathbf{J}(\hat{\boldsymbol{\theta}})$ and $\mathbf{f}(\hat{\boldsymbol{\theta}})$ are nonlinear functions of $\hat{\boldsymbol{\theta}}$, this system of $p$ nonlinear equations admits no closed-form analytical solution and must be solved iteratively.

1. The Gauss–Newton Taylor Linearization Algorithm

Given a current parameter iterate $\boldsymbol{\theta}^{(m)}$ at iteration $m$, we approximate the nonlinear expectation vector $\mathbf{f}(\boldsymbol{\theta})$ around $\boldsymbol{\theta}^{(m)}$ via a first-order multivariate Taylor series expansion:

Substituting this linear approximation into $\mathbf{Y} = \mathbf{f}(\boldsymbol{\theta}) + \boldsymbol{\varepsilon}$ and rearranging terms yields the Local Linearized Pseudo-Model at iteration $m$:

Notice that for fixed $\boldsymbol{\theta}^{(m)}$, estimating the parameter increment vector $\boldsymbol{\delta}^{(m)} = \boldsymbol{\theta}^{(m+1)} - \boldsymbol{\theta}^{(m)}$ is simply an Ordinary Least Squares regression of the current residual vector $\mathbf{e}^{(m)} = \mathbf{Y} - \mathbf{f}(\boldsymbol{\theta}^{(m)})$ on the current $n \times p$ Jacobian matrix $\mathbf{J}_m = \mathbf{J}(\boldsymbol{\theta}^{(m)})$! Solving the linearized normal equations $(\mathbf{J}_m^T \mathbf{J}_m) \hat{\boldsymbol{\delta}}^{(m)} = \mathbf{J}_m^T \mathbf{e}^{(m)}$ gives the Gauss–Newton Iteration:

To understand how Gauss–Newton relates to full Newton–Raphson optimization, observe that the exact $p \times p$ Hessian matrix of $S(\boldsymbol{\theta})$ is:

Because $\mathbb{E}[e_i(\boldsymbol{\theta})] \approx 0$ near the true parameter $\boldsymbol{\theta}$ (and the residuals are small when the model fits well), Gauss–Newton drops the second-derivative intrinsic curvature term $-2 \sum e_i \nabla^2 f_i$ and approximates the Hessian by the positive semi-definite matrix $2 \mathbf{J}_m^T \mathbf{J}_m$ (which is identical to Fisher's Method of Scoring under normal errors!).

2. Levenberg–Marquardt Trust-Region Damping

When the initial guess $\boldsymbol{\theta}^{(0)}$ is far from the optimum or when the Jacobian columns are nearly collinear, $\mathbf{J}_m^T \mathbf{J}_m$ can be nearly singular, causing the undamped Gauss–Newton step $\hat{\boldsymbol{\delta}}^{(m)}$ to overshoot wildly and increase $S(\boldsymbol{\theta})$. Kenneth Levenberg (1944) and Donald Marquardt (1963) resolved this instability by adding a positive ridge damping parameter $\lambda_m > 0$ (or diagonal scaling $\lambda_m \text{diag}(\mathbf{J}_m^T \mathbf{J}_m)$) to the approximate Hessian:

The Levenberg–Marquardt Algorithm interpolates smoothly between steepest-descent gradient descent (as $\lambda_m \to \infty$, $\hat{\boldsymbol{\delta}}^{(m)} \approx \frac{1}{\lambda_m} \mathbf{J}_m^T \mathbf{e}^{(m)} = -\frac{1}{2\lambda_m} \nabla S$, guaranteeing descent in a small trust region) and quadratic Gauss–Newton convergence (as $\lambda_m \to 0$ near the minimum).

3. Asymptotic Distribution and Inference at Convergence

In finite samples, the NLS estimator $\hat{\boldsymbol{\theta}}$ is slightly biased ($\mathbb{E}[\hat{\boldsymbol{\theta}}] \neq \boldsymbol{\theta}$), the residuals do not sum to zero in general ($\sum e_i \neq 0$ unless $\mathbf{1}_n \in \mathcal{C}(\mathbf{J}(\hat{\boldsymbol{\theta}}))$), and $\text{SST} \neq \text{SSR} + \text{SSE}$. However, as $n \to \infty$ under Jennrich's regularity conditions ($\frac{1}{n}\mathbf{J}(\boldsymbol{\theta})^T\mathbf{J}(\boldsymbol{\theta}) \to \mathbf{Q}(\boldsymbol{\theta}) \succ \mathbf{0}$), $\hat{\boldsymbol{\theta}}$ is consistent and asymptotically normal:

As demonstrated in Unit 4, Ordinary Least Squares is exquisitely sensitive to vertical outliers because its quadratic loss function $\rho_{\text{OLS}}(e_i) = \frac{1}{2} e_i^2$ has an **unbounded influence function** $\psi_{\text{OLS}}(e_i) = \rho_{\text{OLS}}'(e_i) = e_i$: a single gross outlier ($|e_i| \to \infty$) exerts infinite pull on the normal equations $\sum_{i=1}^n \psi(e_i) \mathbf{x}_i = \mathbf{0}_p$, resulting in a finite-sample breakdown point of $\varepsilon^* = \frac{1}{n} \to 0$. When disturbances follow heavy-tailed distributions (such as Laplace, Cauchy, or Student $t_3$) or contain unremovable contamination, **Robust Regression** replaces the squared loss $\frac{1}{2}e_i^2$ with a symmetric, less rapidly growing loss function $\rho(u)$. Let $\hat{\sigma}$ denote a robust, outlier-resistant scale estimator of the residuals—most commonly the normalized **Median Absolute Deviation (MAD)**: $$\hat{\sigma}_{\text{MAD}} = \frac{\text{median}_{1 \le i \le n} \left| e_i - \text{median}(e) \right|}{0.6745}$$ where the divisor $\Phi^{-1}(0.75) \approx 0.6745$ makes $\hat{\sigma}_{\text{MAD}}$ a consistent estimator of $\sigma$ when errors are $\mathcal{N}(0, \sigma^2)$ while maintaining a $50\%$ breakdown point! Defining the scaled residual $u_i = \frac{Y_i - \mathbf{x}_i^T \boldsymbol{\beta}}{\hat{\sigma}}$, Peter Huber (1964, 1973) formulated the class of **Regression M-Estimators** (Maximum-Likelihood-Type Estimators), which minimize: $$Q_M(\boldsymbol{\beta}) = \sum_{i=1}^n \rho\left( \frac{Y_i - \mathbf{x}_i^T \boldsymbol{\beta}}{\hat{\sigma}} \right)$$ Differentiating with respect to $\boldsymbol{\beta}$ and setting the gradient to zero yields the **$p$ Robust Estimating Equations** in terms of the **score (influence) function** $\psi(u) = \rho'(u)$: $$\sum_{i=1}^n \psi\left( \frac{Y_i - \mathbf{x}_i^T \hat{\boldsymbol{\beta}}}{\hat{\sigma}} \right) \mathbf{x}_i = \mathbf{0}_p$$ ### The Iteratively Reweighted Least Squares (IRLS) Algorithm for M-Estimation How do we solve the nonlinear system $\sum_{i=1}^n \psi(u_i) \mathbf{x}_i = \mathbf{0}_p$? Notice the algebraic trick of defining the **robust weight function** $w(u) = \frac{\psi(u)}{u}$ (with $w(0) = \lim_{u \to 0} \frac{\psi(u)}{u} = \psi'(0)$), so that $\psi(u_i) = w(u_i) u_i = w_i \left(\frac{Y_i - \mathbf{x}_i^T \hat{\boldsymbol{\beta}}}{\hat{\sigma}}\right)$. Substituting this into the estimating equations and multiplying by $\hat{\sigma}$ transforms the system into: $$\sum_{i=1}^n w_i \left( Y_i - \mathbf{x}_i^T \hat{\boldsymbol{\beta}} \right) \mathbf{x}_i = \mathbf{0}_p \iff \left( \mathbf{X}^T \mathbf{W}(\hat{\boldsymbol{\beta}}) \mathbf{X} \right) \hat{\boldsymbol{\beta}} = \mathbf{X}^T \mathbf{W}(\hat{\boldsymbol{\beta}}) \mathbf{Y}$$ where $\mathbf{W}(\hat{\boldsymbol{\beta}}) = \text{diag}(w_1, w_2, \dots, w_n)$! This immediately generates the **Iteratively Reweighted Least Squares (IRLS)** algorithm: 1. Start with an initial robust estimate $\hat{\boldsymbol{\beta}}^{(0)}$ (such as $L_1$ regression), compute residuals $e_i^{(m)}$ and $\hat{\sigma}_{\text{MAD}}$, and evaluate standardized residuals $u_i^{(m)} = e_i^{(m)} / \hat{\sigma}_{\text{MAD}}$. 2. Compute the observation weights $w_i^{(m)} = w(u_i^{(m)}) = \psi(u_i^{(m)}) / u_i^{(m)}$ and form $\mathbf{W}^{(m)} = \text{diag}(w_1^{(m)}, \dots, w_n^{(m)})$. 3. Update the parameter vector via Weighted Least Squares: $\hat{\boldsymbol{\beta}}^{(m+1)} = \left(\mathbf{X}^T \mathbf{W}^{(m)} \mathbf{X}\right)^{-1} \mathbf{X}^T \mathbf{W}^{(m)} \mathbf{Y}$, repeating until convergence! Note that while monotone M-estimators (Huber) and $L_1$ regression bound the influence of vertical outliers in $Y$, they remain vulnerable to **bad high-leverage points in $\mathbf{X}$** because the estimating equation multiplies $\psi(u_i)$ by the unbounded covariate vector $\mathbf{x}_i$. Achieving high breakdown ($50\%$) against simultaneous $(\mathbf{x}_i, Y_i)$ leverage contamination requires **Generalized M-estimators (Bounded-Influence GM / Mallows–Schweppe)**, Peter Rousseeuw's **Least Trimmed Squares (LTS)** (which minimizes the sum of the smallest $h \approx \lfloor n/2 \rfloor + \lfloor (p+1)/2 \rfloor$ squared residuals), or **MM-Estimation** (which initializes Tukey's biweight IRLS from a $50\%$-breakdown S-estimator).

When the response variable $Y_i$ is intrinsically discrete or bounded—such as a binary clinical outcome $Y_i \in \{0, 1\}$, a binomial proportion $Y_i / m_i \in [0, 1]$, a non-negative integer Poisson count $Y_i \in \{0, 1, 2, \dots\}$, or a strictly positive skewed Gamma survival time $Y_i > 0$—classical linear regression $Y_i = \mathbf{x}_i^T \boldsymbol{\beta} + \varepsilon_i$ breaks down fundamentally: an unbounded linear predictor $\mathbf{x}_i^T \boldsymbol{\beta} \in (-\infty, \infty)$ predicts impossible negative counts or probabilities outside $[0, 1]$, and the variance of $Y_i$ is functionally tied to its mean $\mu_i$. In a landmark 1972 paper, **John Nelder and Robert Wedderburn** unified linear, logistic, probit, Poisson, and Gamma regression under the **Generalized Linear Model (GLM)** framework, which consists of three structural components: ### 1. The Random Component: Exponential Dispersion Family Each independent response $Y_i$ follows a probability density or mass function in the **Exponential Dispersion Family** with canonical (natural) parameter $\theta_i$ and dispersion parameter $\phi > 0$: $$f(y_i; \theta_i, \phi) = \exp\left\{ \frac{y_i \theta_i - b(\theta_i)}{a(\phi)} + c(y_i, \phi) \right\}$$ where $a(\phi) = \phi / w_i$ ($w_i$ being a known prior weight) and $b(\theta_i)$ is a twice-differentiable convex **cumulant function**. Applying Bartlett's first and second log-likelihood identities ($\mathbb{E}[\frac{\partial \ell_i}{\partial \theta_i}] = 0$ and $\mathbb{E}[\frac{\partial^2 \ell_i}{\partial \theta_i^2}] + \mathbb{E}[(\frac{\partial \ell_i}{\partial \theta_i})^2] = 0$) immediately determines the mean $\mu_i$ and variance $\text{Var}(Y_i)$ in terms of derivatives of $b(\theta_i)$: $$\frac{\partial \ell_i}{\partial \theta_i} = \frac{Y_i - b'(\theta_i)}{a(\phi)} \implies \mu_i = \mathbb{E}[Y_i] = b'(\theta_i)$$ $$\frac{\partial^2 \ell_i}{\partial \theta_i^2} = -\frac{b''(\theta_i)}{a(\phi)} \implies \text{Var}(Y_i) = a(\phi) \, b''(\theta_i) = a(\phi) \, V(\mu_i)$$ where $V(\mu_i) = b''\left((b')^{-1}(\mu_i)\right) = \frac{d\mu_i}{d\theta_i}$ is the **GLM Variance Function**, capturing how the variance of $Y_i$ depends on its mean $\mu_i$! ### 2. The Systematic Component: Linear Predictor $\eta_i$ The $p$ covariates enter purely through the linear predictor $\eta_i = \mathbf{x}_i^T \boldsymbol{\beta}$ ($\boldsymbol{\eta} = \mathbf{X}\boldsymbol{\beta}$). ### 3. The Link Function $g(\mu_i) = \eta_i$ A smooth, strictly monotonic differentiable **link function** $g(\cdot)$ maps the restricted domain of the conditional mean $\mu_i = \mathbb{E}[Y_i \mid \mathbf{x}_i]$ onto the entire real line $(-\infty, \infty)$ of the linear predictor: $$\eta_i = g(\mu_i) = \mathbf{x}_i^T \boldsymbol{\beta} \iff \mu_i = g^{-1}(\mathbf{x}_i^T \boldsymbol{\beta})$$ Notice the fundamental distinction between transforming the data $\mathbb{E}[g(Y_i)]$ (Section 7.5) and transforming the population mean $g(\mathbb{E}[Y_i])$ in a GLM: a log-link GLM sets $\ln(\mathbb{E}[Y_i]) = \mathbf{x}_i^T\boldsymbol{\beta}$ without ever taking $\ln(Y_i)$, so observed zeros $Y_i = 0$ pose no problem whatsoever! The special link function $g(\mu_i) = \theta_i = (b')^{-1}(\mu_i)$ that equates the linear predictor directly to the natural parameter ($\theta_i = \eta_i = \mathbf{x}_i^T\boldsymbol{\beta}$) is called the **Canonical Link Function**. Under the canonical link, $\mathbf{X}^T \mathbf{Y}$ is a minimal sufficient statistic for $\boldsymbol{\beta}$, and the observed Hessian equals the expected Fisher Information matrix! ### 4. Maximum Likelihood Estimation via Fisher Scoring and IRLS By the chain rule $\frac{\partial \ell_i}{\partial \beta_j} = \frac{\partial \ell_i}{\partial \theta_i} \frac{d \theta_i}{d \mu_i} \frac{d \mu_i}{d \eta_i} \frac{\partial \eta_i}{\partial \beta_j}$, where $\frac{d\theta_i}{d\mu_i} = \frac{1}{V(\mu_i)}$, $\frac{d\mu_i}{d\eta_i} = \frac{1}{g'(\mu_i)}$, and $\frac{\partial \eta_i}{\partial \beta_j} = x_{ij}$, the **GLM Score Equations** and **Fisher Information Matrix** $\mathcal{I}(\boldsymbol{\beta})$ are: $$\mathbf{U}(\boldsymbol{\beta}) = \frac{1}{\phi} \sum_{i=1}^n \frac{Y_i - \mu_i}{V(\mu_i) \, g'(\mu_i)} \mathbf{x}_i = \mathbf{0}_p, \qquad \mathcal{I}(\boldsymbol{\beta}) = \frac{1}{\phi} \mathbf{X}^T \mathbf{W} \mathbf{X}$$ where $\mathbf{W} = \text{diag}(w_1, \dots, w_n)$ has **GLM working weights** $w_i = \frac{1}{V(\mu_i) [g'(\mu_i)]^2}$ (which simplify to $w_i = V(\mu_i)$ under the canonical link since $g'(\mu_i) = 1/V(\mu_i)$!).Defining the **working dependent variable (adjusted response)** $z_i = \eta_i + (Y_i - \mu_i) g'(\mu_i)$, each Fisher Scoring step $\boldsymbol{\beta}^{(m+1)} = \boldsymbol{\beta}^{(m)} + [\mathcal{I}(\boldsymbol{\beta}^{(m)})]^{-1} \mathbf{U}(\boldsymbol{\beta}^{(m)})$ reduces to an **Iteratively Reweighted Least Squares (IRLS)** update: $$\boldsymbol{\beta}^{(m+1)} = \left( \mathbf{X}^T \mathbf{W}^{(m)} \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{W}^{(m)} \mathbf{z}^{(m)}, \qquad \widehat{\text{Cov}}(\hat{\boldsymbol{\beta}}) = \hat{\phi} \left( \mathbf{X}^T \hat{\mathbf{W}} \mathbf{X} \right)^{-1}$$

When the response variable is a binary Bernoulli indicator $Y_i \in \{0, 1\}$ (or a grouped binomial count $Y_i \sim \text{Binomial}(m_i, \pi_i)$) representing the occurrence of an event such as disease presence, loan default, or customer churn, our target is the conditional success probability $\pi(\mathbf{x}_i) = \mathbb{P}(Y_i = 1 \mid \mathbf{x}_i) = \mathbb{E}[Y_i \mid \mathbf{x}_i] \in (0, 1)$.

Fitting a naive Linear Probability Model (LPM) $\pi(\mathbf{x}_i) = \mathbf{x}_i^T \boldsymbol{\beta}$ via Ordinary Least Squares suffers from three fatal statistical defects:

1. Unbounded Predictions: Because $\mathbf{x}_i^T \boldsymbol{\beta}$ is linear and unbounded, it inevitably predicts nonsensical probabilities $\hat{\pi}_i < 0$ or $\hat{\pi}_i > 1$ for extreme covariate values.

2. Inherent Heteroscedasticity: Since $Y_i \in \{0, 1\}$ is Bernoulli, its error term $\varepsilon_i = Y_i - \pi_i$ takes only two discrete values ($1 - \pi_i$ and $-\pi_i$) with variance $\text{Var}(\varepsilon_i) = \pi_i(1 - \pi_i)$, violating homoscedasticity and invalidating OLS standard errors.

3. Constant Marginal Effects: A one-unit increase in $x_j$ changes $\pi$ by $\beta_j$ identically whether $\pi = 0.50$ or $\pi = 0.99$.

1. The Logistic Sigmoid Function, Logit Link, and Odds Ratio Interpretation

Binary Logistic Regression resolves all three defects by modeling $\pi(\mathbf{x}_i)$ via the standard logistic cumulative distribution function (sigmoid curve) $\pi(\eta) = \frac{e^\eta}{1 + e^\eta} = \frac{1}{1 + e^{-\eta}} \in (0, 1)$, whose inverse is the canonical Logit (Log-Odds) Link Function:

Here, the ratio $\Omega(\mathbf{x}_i) = \frac{\pi(\mathbf{x}_i)}{1 - \pi(\mathbf{x}_i)} = \exp(\mathbf{x}_i^T \boldsymbol{\beta}) \in (0, \infty)$ is the Odds of Success at covariate vector $\mathbf{x}_i$. How do we interpret an individual slope coefficient $\beta_j$? Increase predictor $x_j$ by $1$ unit while holding all other covariates constant:

Thus, $\beta_j$ is the change in log-odds per unit increase in $x_j$, and $\exp(\beta_j)$ is the Multiplicative Odds Ratio ($\text{OR}_j$)! If $\beta_j > 0$, $e^{\beta_j} > 1$ (odds increase by $100(e^{\beta_j} - 1)\%$ per unit increase in $x_j$); if $\beta_j = 0$, $e^{\beta_j} = 1$ (odds unchanged); if $\beta_j < 0$, $0 < e^{\beta_j} < 1$ (odds decrease). Furthermore, on the probability scale, differentiating $\pi(\mathbf{x})$ with respect to $x_j$ gives the Marginal Probability Effect:

which achieves its steepest slope $\frac{1}{4}\beta_j$ at the median effective dose ($\pi = 0.50$, where $\mathbf{x}^T\boldsymbol{\beta} = 0$, so $\text{LD}_{50} = -\beta_0 / \beta_1$ in simple logistic regression) and flattens gracefully toward zero as $\pi \to 0$ or $\pi \to 1$.

2. Log-Likelihood, IRLS Working Equations, Wald Tests, and Deviance Goodness-of-Fit

For $n$ independent Bernoulli observations $Y_i \in \{0, 1\}$, the log-likelihood function is:

Setting the score vector to zero gives the Logistic Likelihood Equations $\mathbf{U}(\hat{\boldsymbol{\beta}}) = \mathbf{X}^T (\mathbf{Y} - \hat{\boldsymbol{\pi}}) = \mathbf{0}_p$. Specializing the GLM IRLS algorithm of Section 8.4 to the canonical logit link ($V(\pi_i) = \pi_i(1 - \pi_i)$, $g'(\pi_i) = \frac{1}{\pi_i(1 - \pi_i)}$, $\phi = 1$), the IRLS weight matrix and working response at iteration $m$ are:

At convergence, $\hat{\boldsymbol{\beta}} \overset{\text{approx}}{\sim} \mathcal{N}_p\left(\boldsymbol{\beta}, (\mathbf{X}^T \hat{\mathbf{W}} \mathbf{X})^{-1}\right)$. Hypothesis testing in logistic regression employs:

1. Wald Test: $z_j = \frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)} \overset{d}{\to} \mathcal{N}(0, 1)$ (or $W_j = z_j^2 \sim \chi_1^2$), and the $95\%$ CI for the Odds Ratio is $\exp\left(\hat{\beta}_j \pm 1.96 \, \text{SE}(\hat{\beta}_j)\right)$. *(Caution: the Hauck–Donner effect shows that for huge $|\hat{\beta}_j|$, $\hat{\pi}_i(1 - \hat{\pi}_i) \to 0$ inflates $\text{SE}(\hat{\beta}_j)$, making the Wald statistic lose power; the Likelihood Ratio / Deviance test is always preferred!).*

2. Residual Deviance $D$ and Likelihood Ratio Drop-in-Deviance Test: The Residual Deviance is twice the log-likelihood ratio between the saturated model ($\tilde{\pi}_i = Y_i / m_i$) and the fitted model:

To test a nested Reduced model $M_R$ ($p_R$ parameters) against a Full model $M_F$ ($p_F$ parameters), the Drop-in-Deviance (Likelihood Ratio) Statistic is:

For ungrouped binary data ($m_i = 1$), overall goodness of fit is evaluated via the Hosmer–Lemeshow Test, which partitions observations into $G = 10$ deciles of predicted risk $\hat{\pi}_i$ and compares observed vs. expected event counts via a $\chi_{G-2}^2 = \chi_8^2$ statistic.

When the response variable $Y_i \in \{0, 1, 2, \dots\}$ represents a non-negative integer count of discrete events occurring within a fixed spatial or temporal window—such as hospital emergency admissions per day, highway traffic accidents per intersection, or defective microchips per silicon wafer—the benchmark generalized linear model is Poisson Log-Linear Regression.

1. The Poisson Log-Linear Model, Rate Ratios, and Exposure Offsets

Assume $Y_i \sim \text{Poisson}(\mu_i)$ independently with probability mass function $f(y_i; \mu_i) = \frac{e^{-\mu_i} \mu_i^{y_i}}{y_i!}$ for $y_i = 0, 1, 2, \dots$, where the mean and variance satisfy the equidispersion property $\mathbb{E}[Y_i] = \text{Var}(Y_i) = \mu_i > 0$. Using the canonical logarithmic link function $g(\mu_i) = \ln(\mu_i) = \eta_i$ guarantees that the predicted mean count $\mu_i = \exp(\mathbf{x}_i^T \boldsymbol{\beta}) > 0$ is strictly positive for all $\boldsymbol{\beta} \in \mathbb{R}^p$:

Because the covariates act multiplicatively on the expected count $\mu_i$, exponentiating an individual regression slope $\beta_j$ yields the Incidence Rate Ratio ($\text{IRR}_j$) (or Relative Risk):

That is, a one-unit increase in $x_j$ multiplies the expected count by $e^{\beta_j}$ (a $100(e^{\beta_j} - 1)\%$ percentage change).

In many epidemiological and actuarial studies, the counts $Y_i$ are observed over unequal exposure denominators $t_i > 0$ (such as person-years at risk, population size, or kilometers of highway), and we wish to model the underlying event rate $\lambda_i = \frac{\mu_i}{t_i}$. Taking logarithms of $\frac{\mu_i}{t_i} = \exp(\mathbf{x}_i^T \boldsymbol{\beta})$ gives:

Here, $\ln(t_i)$ enters the linear predictor as an Offset Variable—a known regressor whose slope coefficient is constrained to equal exactly $1.0$ without being estimated!

2. Score Equations, IRLS Weights, and Poisson Deviance

The Poisson log-likelihood is $\ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[ Y_i (\mathbf{x}_i^T\boldsymbol{\beta}) - \exp(\mathbf{x}_i^T\boldsymbol{\beta}) - \ln(Y_i!) \right]$, yielding the score equations $\mathbf{U}(\hat{\boldsymbol{\beta}}) = \mathbf{X}^T(\mathbf{Y} - \hat{\boldsymbol{\mu}}) = \mathbf{0}_p$, the IRLS weight matrix $\mathbf{W} = \text{diag}(\mu_1, \dots, \mu_n)$, working response $z_i = \eta_i + \frac{Y_i - \mu_i}{\mu_i}$, and asymptotic covariance $\widehat{\text{Cov}}(\hat{\boldsymbol{\beta}}) = (\mathbf{X}^T \hat{\mathbf{W}} \mathbf{X})^{-1}$. The Poisson Residual Deviance and Pearson Chi-Square Statistic are:

(with $0 \ln(0) \equiv 0$). Notice that whenever an intercept is included, the first score equation $\sum_{i=1}^n (Y_i - \hat{\mu}_i) = 0$ forces the sum of fitted counts to equal the sum of observed counts ($\sum \hat{\mu}_i = \sum Y_i$), simplifying the deviance to $D = 2 \sum_{i=1}^n Y_i \ln(Y_i / \hat{\mu}_i)$.

3. Overdispersion Diagnostics and Quasi-Poisson vs Negative Binomial Remedies

In real-world count data, unobserved individual heterogeneity (frailty) or event clustering almost always causes the empirical variance to exceed the conditional mean ($\text{Var}(Y_i) > \mu_i$), a phenomenon called Overdispersion. When a well-specified Poisson model holds ($\phi = 1$), both $D / (n - p)$ and $X^2 / (n - p)$ should hover near $1.0$. If the estimated dispersion parameter $\hat{\phi}$:

exceeds unity substantially, standard Poisson standard errors are severely biased downward (producing spuriously large Wald $z$-statistics). Two standard remedies are used:

1. Quasi-Poisson GLM (Wedderburn Quasi-Likelihood): Postulates $\mathbb{E}[Y_i] = \mu_i$ and $\text{Var}(Y_i) = \phi \mu_i$ ($\phi > 1$). Because $\phi$ cancels out of the score equation $\frac{1}{\phi}\mathbf{X}^T(\mathbf{Y} - \boldsymbol{\mu}) = \mathbf{0}_p$, the point estimates $\hat{\boldsymbol{\beta}}$ are identical to ordinary Poisson MLE, while all standard errors are multiplied by the overdispersion inflation factor $\sqrt{\hat{\phi}}$, and drop-in-deviance tests use $F = \frac{(D_R - D_F)/q}{\hat{\phi}} \sim F_{q, n-p}$.

2. Negative Binomial Regression (Poisson–Gamma Mixture): Assumes $Y_i \mid u_i \sim \text{Poisson}(\mu_i u_i)$ where the unobserved frailty $u_i \sim \text{Gamma}(\alpha^{-1}, \alpha^{-1})$ has mean $1$ and variance $\alpha > 0$. Integrating out $u_i$ yields a Negative Binomial marginal distribution with mean $\mathbb{E}[Y_i] = \mu_i$ and quadratic variance function $\text{Var}(Y_i) = \mu_i + \alpha \mu_i^2$, allowing full likelihood inference and testing $H_0: \alpha = 0$.