Generalized Linear Models: Zero to Hero Tutorial

This comprehensive tutorial takes you from the foundational concepts of Generalized Linear Models (GLMs) all the way through advanced model specification, estimation, diagnostics, and practical usage within the DataStatPro application. Whether you are a complete beginner or an experienced analyst, this guide is structured to build your understanding step by step.

Table of Contents

1. Prerequisites and Background Concepts
2. What are Generalized Linear Models?
3. The Mathematical Framework of GLMs
4. The Exponential Family of Distributions
5. Link Functions
6. GLM Distributions and Their Applications
7. Assumptions of GLMs
8. Parameter Estimation: Maximum Likelihood and IRLS
9. Model Fit and Evaluation
10. Hypothesis Testing and Inference
11. Model Diagnostics and Residuals
12. Model Selection and Variable Selection
13. Overdispersion and Underdispersion
14. Using the GLM Component
15. Computational and Formula Details
16. Worked Examples
17. Common Mistakes and How to Avoid Them
18. Troubleshooting
19. Quick Reference Cheat Sheet

1. Prerequisites and Background Concepts

Before diving into Generalized Linear Models, it is helpful to be familiar with the following foundational concepts. Do not worry if you are not — each concept is briefly explained here.

1.1 Probability Distributions

A probability distribution describes how the values of a random variable are distributed. Key distributions used in GLMs:

• Normal (Gaussian): Continuous, symmetric, unbounded. Characterised by mean $\mu$ and variance $\sigma^2$.
• Binomial: Discrete, counts successes in $n$ independent Bernoulli trials. Characterised by $n$ and probability $p$.
• Poisson: Discrete, counts of events in a fixed interval. Characterised by rate $\lambda$, with $E[Y] = \text{Var}(Y) = \lambda$.
• Gamma: Continuous, positive-valued, right-skewed. Characterised by shape $\alpha$ and rate $\beta$.
• Inverse Gaussian: Continuous, positive-valued, highly right-skewed. Models first-passage times.
• Negative Binomial: Discrete, counts with overdispersion relative to Poisson.

1.2 The Likelihood Function

The likelihood function $L(\boldsymbol{\theta}; \mathbf{y})$ measures how probable the observed data $\mathbf{y}$ are, given a parameter vector $\boldsymbol{\theta}$. For independent observations:

$L(\boldsymbol{\theta}; \mathbf{y}) = \prod_{i=1}^n f(y_i; \boldsymbol{\theta})$

The log-likelihood is more convenient to work with:

$\ell(\boldsymbol{\theta}; \mathbf{y}) = \sum_{i=1}^n \ln f(y_i; \boldsymbol{\theta})$

Maximum Likelihood Estimation (MLE) finds the parameter values that maximise $\ell(\boldsymbol{\theta}; \mathbf{y})$.

1.3 The Linear Predictor

A linear predictor $\eta$ is a weighted linear combination of predictor variables:

$\eta = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p = \mathbf{x}^T \boldsymbol{\beta}$

This is the core structure inherited from linear regression. In GLMs, $\eta$ is not the outcome itself but is transformed through a link function to relate to the mean of the response distribution.

1.4 Ordinary Linear Regression Recap

In ordinary linear regression (OLS):

$Y_i = \beta_0 + \beta_1 X_{i1} + \dots + \beta_p X_{ip} + \epsilon_i, \quad \epsilon_i \sim \mathcal{N}(0, \sigma^2)$

This model has three implicit components:

1. A distribution for the response: $Y_i \sim \mathcal{N}(\mu_i, \sigma^2)$.
2. A linear predictor: $\eta_i = \mathbf{x}_i^T \boldsymbol{\beta}$.
3. A link function connecting $\mu_i$ to $\eta_i$: $\eta_i = \mu_i$ (the identity link).

GLMs generalise this framework by allowing different distributions and link functions.

1.5 The Score Equations and the Information Matrix

The score vector is the gradient of the log-likelihood with respect to the parameters:

$\mathbf{s}(\boldsymbol{\beta}) = \frac{\partial \ell}{\partial \boldsymbol{\beta}}$

Setting $\mathbf{s}(\boldsymbol{\beta}) = \mathbf{0}$ gives the MLE. The Fisher information matrix is:

$\mathcal{I}(\boldsymbol{\beta}) = -E\left[\frac{\partial^2 \ell}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^T}\right]$

Its inverse $\mathcal{I}^{-1}(\boldsymbol{\beta})$ gives the asymptotic covariance matrix of the MLE, used to compute standard errors and confidence intervals.

2. What are Generalized Linear Models?

Generalized Linear Models (GLMs) are a unified class of regression models that extend ordinary linear regression to accommodate response variables with non-normal distributions. Introduced by Nelder and Wedderburn (1972), GLMs provide a coherent framework for modelling a wide variety of outcome types — counts, proportions, binary outcomes, continuous positive values, and more — using a single, elegant mathematical structure.

2.1 The Central Idea

Ordinary linear regression assumes the response $Y$ is normally distributed and that the mean $\mu$ equals the linear predictor directly: $\mu = \eta$. GLMs relax both restrictions:

1. The distribution of $Y$ can be any member of the exponential family (Normal, Binomial, Poisson, Gamma, Inverse Gaussian, etc.).
2. The link function $g(\mu) = \eta$ can be any monotone, differentiable function that maps the mean $\mu$ (constrained to its natural range) to the real line $(-\infty, +\infty)$.

This two-step generalisation unlocks an enormous range of practical modelling scenarios while preserving the interpretability of regression coefficients.

2.2 Real-World Applications

GLMs are among the most widely used statistical models in applied science, business, and policy:

• Insurance & Actuarial Science: Modelling claim counts (Poisson/Negative Binomial), claim severity (Gamma), and pure premiums (Tweedie).
• Public Health & Epidemiology: Modelling disease incidence rates (Poisson with offset), binary disease outcomes (Binomial/logistic), survival and time-to-event data.
• Ecology: Modelling species counts (Poisson, Negative Binomial), presence/absence (Binomial), and biomass (Gamma).
• Economics & Finance: Modelling discrete choices (Binomial), income (Gamma), financial durations (Inverse Gaussian).
• Marketing: Modelling purchase counts (Poisson), click-through rates (Binomial), and customer lifetime value (Gamma/Tweedie).
• Clinical Trials: Modelling adverse event counts (Poisson), binary treatment outcomes (Binomial), and length of stay (Gamma).
• Manufacturing & Quality Control: Modelling defect counts (Poisson), product lifetimes (Gamma/Inverse Gaussian).
• Social Sciences: Modelling ordered survey responses (Ordinal), multinomial choices (Multinomial), and event rates.

2.3 The Three Components of a GLM

Every GLM is fully specified by three components:

2.4 How GLMs Generalise Linear Regression

3. The Mathematical Framework of GLMs

3.1 The Three-Component Structure in Detail

A GLM specifies that the $i$-th response $Y_i$ has:

Random Component: $Y_i \sim f(y_i; \theta_i, \phi) \quad \text{(exponential family distribution)}$

With mean $E[Y_i] = \mu_i$ and variance $\text{Var}(Y_i) = \phi \cdot V(\mu_i)$, where:

• $\phi$ is the dispersion parameter (estimated or known).
• $V(\mu)$ is the variance function — a function of the mean that characterises the distribution.

Systematic Component: $\eta_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \dots + \beta_p X_{ip} = \mathbf{x}_i^T \boldsymbol{\beta}$

Link Function: $g(\mu_i) = \eta_i$

So the mean is related to the predictors through: $\mu_i = g^{-1}(\eta_i) = g^{-1}(\mathbf{x}_i^T \boldsymbol{\beta})$

Where $g^{-1}$ is the inverse link function (also called the mean function or response function).

3.2 The Variance Function $V(\mu)$

The variance function characterises how the variance of the response depends on its mean. Each exponential family distribution has a specific variance function:

3.3 The Dispersion Parameter $\phi$

The full variance of $Y_i$ is:

$\text{Var}(Y_i) = \frac{\phi \cdot V(\mu_i)}{w_i}$

Where $w_i$ is a known prior weight (e.g., $w_i = n_i$ for binomial proportions). The dispersion parameter $\phi$:

• Is known for Poisson ($\phi = 1$) and Binomial ($\phi = 1$).
• Is estimated for Normal, Gamma, and Inverse Gaussian.
• Can be estimated for Poisson and Binomial to account for overdispersion (quasi-GLM).

3.4 The Canonical Link

For each distribution, there is a canonical link function that arises naturally from the mathematical structure of the exponential family. Using the canonical link has desirable statistical properties (sufficient statistics, simpler score equations):

Non-canonical links can also be used and may be more interpretable in certain applications. The canonical link is the default in most GLM software but is not obligatory.

3.5 The Offset

An offset is a term added to the linear predictor with a fixed coefficient of 1:

$\eta_i = \mathbf{x}_i^T \boldsymbol{\beta} + \text{offset}_i$

Offsets are used when the response is a rate and observations have different exposure times or population sizes. The offset is known (not estimated) and is included on the linear predictor scale.

Example: Modelling disease incidence counts $Y_i$ for regions with different population sizes $N_i$. The rate per person is $\lambda_i = \mu_i / N_i$. Using a Poisson model with log link:

$\ln(\mu_i) = \mathbf{x}_i^T \boldsymbol{\beta} + \ln(N_i)$

Where $\ln(N_i)$ is the offset. The model then estimates the log rate $\ln(\lambda_i) = \mathbf{x}_i^T \boldsymbol{\beta}$.

4. The Exponential Family of Distributions

The exponential family is a broad class of distributions that share a common mathematical form, which is the foundation of the GLM framework.

4.1 The Exponential Family Form

A distribution belongs to the exponential family if its probability density (or mass) function can be written as:

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

Where:

• $\theta$ = natural (canonical) parameter — a function of the mean $\mu$.
• $\phi$ = dispersion parameter.
• $b(\theta)$ = cumulant function (log-normalising constant).
• $a(\phi)$ = dispersion function (typically $\phi/w$ for prior weight $w$).
• $c(y, \phi)$ = a function of the data and dispersion only (not $\theta$).

Key relationships derived from $b(\theta)$:

$\mu = E[Y] = b'(\theta) \quad \text{(first derivative of } b\text{)}$

$\text{Var}(Y) = a(\phi) \cdot b''(\theta) = \phi \cdot V(\mu) \quad \text{(second derivative of } b\text{)}$

This elegant structure means that all moments of the distribution follow automatically from the cumulant function $b(\theta)$.

4.2 Major Exponential Family Distributions for GLMs

4.2.1 Normal Distribution

$f(y; \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{(y-\mu)^2}{2\sigma^2}\right\}$

• Natural parameter: $\theta = \mu$
• Cumulant function: $b(\theta) = \theta^2/2$
• Variance function: $V(\mu) = 1$
• Dispersion: $\phi = \sigma^2$
• Support: $(-\infty, +\infty)$

4.2.2 Binomial Distribution

$f(y; n, p) = \binom{n}{y} p^y (1-p)^{n-y}, \quad y \in \{0, 1, \dots, n\}$

• Natural parameter: $\theta = \ln(p/(1-p))$ (log odds)
• Cumulant function: $b(\theta) = n \ln(1 + e^\theta)$
• Variance function: $V(\mu) = \mu(1-\mu)$ (where $\mu = p$)
• Dispersion: $\phi = 1$ (known)
• Support: $\{0, 1, \dots, n\}$

4.2.3 Poisson Distribution

$f(y; \lambda) = \frac{e^{-\lambda} \lambda^y}{y!}, \quad y \in \{0, 1, 2, \dots\}$

• Natural parameter: $\theta = \ln(\lambda)$
• Cumulant function: $b(\theta) = e^\theta$
• Variance function: $V(\mu) = \mu$
• Dispersion: $\phi = 1$ (known)
• Support: $\{0, 1, 2, \dots\}$

4.2.4 Gamma Distribution

$f(y; \alpha, \beta) = \frac{y^{\alpha-1} e^{-y\beta} \beta^\alpha}{\Gamma(\alpha)}, \quad y > 0$

• Natural parameter: $\theta = -1/\mu$ (inverse of mean)
• Cumulant function: $b(\theta) = -\ln(-\theta)$
• Variance function: $V(\mu) = \mu^2$
• Dispersion: $\phi = 1/\alpha$ (inverse of shape)
• Support: $(0, +\infty)$

4.2.5 Inverse Gaussian Distribution

$f(y; \mu, \lambda) = \sqrt{\frac{\lambda}{2\pi y^3}} \exp\left\{-\frac{\lambda(y-\mu)^2}{2\mu^2 y}\right\}, \quad y > 0$

• Natural parameter: $\theta = -1/(2\mu^2)$
• Cumulant function: $b(\theta) = -\sqrt{-2\theta}$
• Variance function: $V(\mu) = \mu^3$
• Dispersion: $\phi = 1/\lambda$
• Support: $(0, +\infty)$

4.3 The Negative Binomial Distribution

While the Negative Binomial is not a member of the exponential family in its most general form, it can be treated as a quasi-exponential family or as a Poisson-Gamma mixture:

$f(y; \mu, k) = \frac{\Gamma(y+k)}{\Gamma(k) y!} \left(\frac{k}{k+\mu}\right)^k \left(\frac{\mu}{k+\mu}\right)^y, \quad y \in \{0, 1, 2, \dots\}$

• Mean: $E[Y] = \mu$
• Variance: $\text{Var}(Y) = \mu + \mu^2/k$
• Overdispersion parameter: $k > 0$ (smaller $k$ → more overdispersion; as $k \to \infty$, Negative Binomial → Poisson)
• Support: $\{0, 1, 2, \dots\}$

4.4 The Tweedie Distribution

The Tweedie distribution is a special case of the exponential dispersion model with power variance function $V(\mu) = \mu^p$:

The Tweedie distribution with $1 < p < 2$ is particularly valuable in insurance (pure premium modelling) and ecology (biomass data) because it naturally accommodates data with a mass at zero and a continuous positive distribution for non-zero values.

5. Link Functions

The link function $g(\mu) = \eta$ is the bridge between the mean of the response distribution and the linear predictor. Choosing an appropriate link function is a key modelling decision.

5.1 Requirements for a Valid Link Function

A valid link function must be:

1. Monotone: Strictly increasing or decreasing.
2. Differentiable: $g'(\mu)$ must exist and be non-zero.
3. Range-compatible: The range of $g^{-1}(\eta)$ must match the support of the mean $\mu$.

5.2 Commonly Used Link Functions

5.3 Logit Link (Binomial GLM)

$g(\mu) = \text{logit}(\mu) = \ln\left(\frac{\mu}{1-\mu}\right)$

• Maps probabilities $(0,1)$ to the real line $(-\infty, +\infty)$.
• Produces odds ratio interpretations: $e^{\beta_j}$ is the multiplicative change in odds per unit increase in $X_j$.
• Symmetric around $\mu = 0.5$.

5.4 Probit Link (Binomial GLM)

$g(\mu) = \Phi^{-1}(\mu)$

Where $\Phi^{-1}$ is the quantile function of the standard normal distribution.

• Maps probabilities $(0,1)$ to the real line via the normal distribution.
• Produces probit (z-score) interpretations.
• Very similar to logit but with lighter tails; differs mainly for extreme probabilities.

5.5 Complementary Log-Log Link (Binomial GLM)

$g(\mu) = \ln(-\ln(1-\mu))$

• Asymmetric: Approaches 1 faster from below than it approaches 0.
• Appropriate when the probability approaches 1 quickly but approaches 0 slowly — common in survival/hazard models.
• Produces a proportional hazards interpretation: $e^{\beta_j}$ is the multiplicative change in the hazard.

5.6 Log Link (Poisson, Negative Binomial, Gamma GLM)

$g(\mu) = \ln(\mu)$

• Ensures $\mu > 0$ (positivity constraint satisfied automatically).
• Produces multiplicative interpretations: $e^{\beta_j}$ is the multiplicative change in the mean per unit increase in $X_j$.
• Most commonly used link for count and continuous positive data.

5.7 Inverse Link (Gamma GLM)

$g(\mu) = \frac{1}{\mu}$

• The canonical link for the Gamma distribution.
• Less commonly used than the log link because the inverse parameterisation is less interpretable.
• Coefficients represent changes in the reciprocal of the mean.

5.8 Choosing the Link Function

6. GLM Distributions and Their Applications

6.1 Binomial GLM (Logistic, Probit, Cloglog Regression)

Use when: Response is a binary outcome (0/1, yes/no, success/failure) or a proportion $y/n$ where both $y$ (successes) and $n$ (trials) are known.

Model:

$Y_i \sim \text{Binomial}(n_i, \mu_i), \quad g(\mu_i) = \eta_i = \mathbf{x}_i^T \boldsymbol{\beta}$

Default link: Logit.

Interpretation (logit link): $e^{\beta_j}$ is the odds ratio — the multiplicative change in the odds of success for a one-unit increase in $X_j$.

Special cases:

• $n_i = 1$ for all $i$: Binary logistic regression.
• $n_i > 1$: Grouped binomial (proportions) regression.

Applications: Disease diagnosis, credit default, customer churn, election outcomes, clinical trial response rates.

6.2 Poisson GLM (Poisson Regression)

Use when: Response is a count of events that could in principle be any non-negative integer, arising from a process with a constant rate.

Model:

$Y_i \sim \text{Poisson}(\mu_i), \quad \ln(\mu_i) = \mathbf{x}_i^T \boldsymbol{\beta} + \text{offset}_i$

Default link: Log.

Interpretation (log link): $e^{\beta_j}$ is the rate ratio (or incidence rate ratio) — the multiplicative change in the expected count for a one-unit increase in $X_j$.

Key assumption: $E[Y_i] = \text{Var}(Y_i) = \mu_i$ (equidispersion). Violations lead to overdispersion (Section 13).

Applications: Number of accidents, hospital admissions, species counts, web page visits, insurance claims frequency.

6.3 Negative Binomial GLM

Use when: Response is a count variable with overdispersion (variance exceeds the mean) — the most common departure from Poisson assumptions.

Model:

$Y_i \sim \text{NegBin}(\mu_i, k), \quad \ln(\mu_i) = \mathbf{x}_i^T \boldsymbol{\beta}$

$\text{Var}(Y_i) = \mu_i + \frac{\mu_i^2}{k}$

Interpretation: Same as Poisson (log link, rate ratios), but with an additional overdispersion parameter $k$ estimated from the data.

Applications: Same as Poisson but when the Poisson assumption of equidispersion is violated — common in ecology (species abundance), healthcare (hospitalisation counts), and insurance.

6.4 Gamma GLM

Use when: Response is continuous and strictly positive, with variance that increases proportionally to the square of the mean (coefficient of variation is roughly constant).

Model:

$Y_i \sim \text{Gamma}(\mu_i, \phi), \quad g(\mu_i) = \eta_i$

Common links: Log (most interpretable), inverse (canonical), identity.

Interpretation (log link): $e^{\beta_j}$ is the multiplicative change in the mean response per unit increase in $X_j$.

Applications: Insurance claim severity (cost per claim), income, hospital costs, reaction times, survival times (without censoring), environmental concentrations.

6.5 Inverse Gaussian GLM

Use when: Response is continuous, strictly positive, and highly right-skewed, with variance increasing as the cube of the mean — more extreme than Gamma.

Model:

$Y_i \sim \text{InverseGaussian}(\mu_i, \phi), \quad g(\mu_i) = \eta_i$

Common links: Inverse squared (canonical), log, inverse.

Applications: First-passage times, repair times, extreme claim sizes, some types of survival data.

6.6 Gaussian GLM (Standard Linear Regression)

Use when: Response is continuous, unbounded, approximately normally distributed, with constant variance.

Model:

$Y_i \sim \mathcal{N}(\mu_i, \sigma^2), \quad \mu_i = \mathbf{x}_i^T \boldsymbol{\beta}$

This is identical to OLS with the identity link. Including it in the GLM framework confirms that linear regression is a special case of GLMs.

Applications: All classical linear regression applications.

6.7 Tweedie GLM

Use when: Response contains exact zeros mixed with positive continuous values (a "zero-inflated" continuous distribution), or when the appropriate power variance is uncertain.

Model:

$Y_i \sim \text{Tweedie}(\mu_i, \phi, p), \quad \ln(\mu_i) = \mathbf{x}_i^T \boldsymbol{\beta}$

The power $p \in (1, 2)$ is estimated from the data or set by domain knowledge.

Applications: Insurance pure premium (frequency × severity), rainfall amounts, ecological biomass, fisheries catch data.

6.8 Quasi-GLMs

When the distributional assumption is uncertain or violated, quasi-GLMs relax the full distributional assumption and specify only the mean and variance function:

$E[Y_i] = \mu_i, \quad \text{Var}(Y_i) = \phi \cdot V(\mu_i)$

The dispersion parameter $\phi$ is estimated from the data (not fixed at 1), providing valid inference even when the count data are overdispersed or underdispersed.

Common quasi-GLMs:

• Quasi-Poisson: Poisson mean function with estimated $\phi$.
• Quasi-Binomial: Binomial mean function with estimated $\phi$.

⚠️ Quasi-GLMs do not have a full likelihood, so AIC/BIC cannot be computed. Use deviance and F-tests for model comparison instead.

6.9 Summary of GLM Distributions

7. Assumptions of GLMs

7.1 Correct Distributional Family

The chosen distribution must be appropriate for the type of response variable. For example:

• Using Gaussian for count data ignores the non-negativity and discreteness.
• Using Poisson for overdispersed counts ignores the excess variance, leading to underestimated standard errors.

How to check: Inspect the response variable's distribution (histogram, range), consider the data-generating process, and verify using residual diagnostics and goodness-of-fit tests.

7.2 Correct Link Function

The link function must be appropriate for the chosen distribution and the expected relationship between predictors and the mean.

How to check: Inspect residual plots; compare alternative link functions using AIC; use added-variable plots for the link function.

7.3 Linearity on the Link Scale

GLMs assume a linear relationship between the predictors $\mathbf{x}_i$ and the transformed mean $g(\mu_i)$:

$g(\mu_i) = \beta_0 + \beta_1 X_{i1} + \dots + \beta_p X_{ip}$

This means the relationship between each $X_j$ and $g(\mu)$ must be linear, even if the relationship between $X_j$ and $\mu$ itself is non-linear.

How to check: Partial residual plots (component-plus-residual plots); LOESS-smoothed plots of residuals vs. each predictor.

7.4 Independence of Observations

Observations must be independent of each other. Clustered, longitudinal, or spatial data may have within-group correlations that violate this assumption.

How to check: Consider the study design. For clustered data, use Generalised Estimating Equations (GEE) or mixed models (GLMM) instead.

7.5 Correct Specification of the Variance Function

The variance function must correctly describe how variability changes with the mean. Misspecification leads to:

• Incorrect standard errors (too small if variance is underestimated).
• Invalid hypothesis tests and confidence intervals.

How to check: Residual vs. fitted value plots; scale-location plots; Pearson $\chi^2$ / deviance tests for dispersion.

7.6 No Perfect Multicollinearity

As in OLS, perfect multicollinearity (one predictor is a perfect linear function of others) prevents estimation. Near-multicollinearity inflates standard errors.

How to check: Variance Inflation Factor (VIF); condition number of the design matrix.

7.7 No Complete Separation (for Binomial GLMs)

For logistic regression and other binomial GLMs, complete separation (a predictor or combination perfectly predicts the outcome) causes the MLE to diverge to $\pm\infty$.

How to check: Warning messages from the fitting algorithm; extremely large coefficient estimates with huge standard errors.

7.8 Sufficient Sample Size

GLM inference is based on asymptotic (large-sample) theory. The adequacy of asymptotic approximations depends on:

• Total sample size $n$.
• For Binomial: Expected counts $n_i \mu_i \geq 5$ in each cell.
• For Poisson: Expected counts $\mu_i \geq 1$ (preferably $\geq 5$) in most cells.

Small expected counts reduce the reliability of likelihood ratio tests, Wald tests, and residual diagnostics.

8. Parameter Estimation: Maximum Likelihood and IRLS

8.1 The Log-Likelihood for GLMs

For independent observations, the log-likelihood is:

$\ell(\boldsymbol{\beta}; \mathbf{y}) = \sum_{i=1}^n \ell_i(\boldsymbol{\beta}; y_i) = \sum_{i=1}^n \frac{y_i \theta_i - b(\theta_i)}{a(\phi)} + c(y_i, \phi)$

Where $\theta_i = \theta(\mu_i)$ and $\mu_i = g^{-1}(\mathbf{x}_i^T \boldsymbol{\beta})$ depend on the regression coefficients $\boldsymbol{\beta}$ through the link function.

8.2 The Score Equations

Setting the gradient of the log-likelihood to zero gives the score equations (MLE conditions):

$\frac{\partial \ell}{\partial \beta_j} = \sum_{i=1}^n \frac{(y_i - \mu_i)}{a(\phi) V(\mu_i)} \frac{\partial \mu_i}{\partial \eta_i} x_{ij} = 0, \quad j = 0, 1, \dots, p$

In matrix form:

$\mathbf{X}^T \mathbf{W}^{1/2} \mathbf{D} (\mathbf{y} - \boldsymbol{\mu}) = \mathbf{0}$

Where $\mathbf{D} = \text{diag}(\partial \mu_i / \partial \eta_i)$ and $\mathbf{W} = \text{diag}(w_i / (V(\mu_i) (g'(\mu_i))^2))$.

These equations are generally non-linear in $\boldsymbol{\beta}$ and require iterative solution.

8.3 Iteratively Reweighted Least Squares (IRLS)

The standard algorithm for fitting GLMs is Iteratively Reweighted Least Squares (IRLS), a Newton-Raphson optimisation applied to the log-likelihood.

At each iteration $t$:

Step 1: Compute the adjusted dependent variable (working response):

$z_i^{(t)} = \hat{\eta}_i^{(t)} + (y_i - \hat{\mu}_i^{(t)}) \frac{d\eta_i}{d\mu_i}\Bigg|_{\hat{\mu}_i^{(t)}}= \hat{\eta}_i^{(t)} + (y_i - \hat{\mu}_i^{(t)}) g'(\hat{\mu}_i^{(t)})$

Step 2: Compute the working weights:

$w_i^{(t)} = \frac{w_i}{V(\hat{\mu}_i^{(t)}) \left[g'(\hat{\mu}_i^{(t)})\right]^2}$

Step 3: Solve the weighted least squares problem:

$\boldsymbol{\beta}^{(t+1)} = \left(\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X}\right)^{-1} \mathbf{X}^T \mathbf{W}^{(t)} \mathbf{z}^{(t)}$

Convergence: Repeat until $\|\boldsymbol{\beta}^{(t+1)} - \boldsymbol{\beta}^{(t)}\| < \epsilon$ (e.g., $\epsilon = 10^{-8}$) or the change in deviance is negligible.

Starting values: Typically $\hat{\mu}_i^{(0)} = y_i + \delta$ (small constant to avoid boundary issues) or the overall mean $\bar{y}$.

8.4 The Fisher Information Matrix and Standard Errors

At convergence, the observed Fisher information matrix is:

$\mathcal{I}(\hat{\boldsymbol{\beta}}) = \mathbf{X}^T \hat{\mathbf{W}} \mathbf{X}$

Where $\hat{\mathbf{W}}$ is the weight matrix evaluated at $\hat{\boldsymbol{\beta}}$. The asymptotic covariance matrix of $\hat{\boldsymbol{\beta}}$ is:

$\text{Cov}(\hat{\boldsymbol{\beta}}) = \phi \left(\mathbf{X}^T \hat{\mathbf{W}} \mathbf{X}\right)^{-1}$

The standard error of $\hat{\beta}_j$:

$SE(\hat{\beta}_j) = \sqrt{\hat{\phi} \left[\left(\mathbf{X}^T \hat{\mathbf{W}} \mathbf{X}\right)^{-1}\right]_{jj}}$

For known-dispersion models (Poisson, Binomial with $\phi = 1$):

$SE(\hat{\beta}_j) = \sqrt{\left[\left(\mathbf{X}^T \hat{\mathbf{W}} \mathbf{X}\right)^{-1}\right]_{jj}}$

8.5 Estimating the Dispersion Parameter

For distributions with estimated dispersion (Normal, Gamma, Inverse Gaussian), $\phi$ is estimated after obtaining $\hat{\boldsymbol{\beta}}$:

Method of Moments (Pearson $\chi^2$):

$\hat{\phi}_{Pearson} = \frac{1}{n - p - 1} \sum_{i=1}^n \frac{(y_i - \hat{\mu}_i)^2}{V(\hat{\mu}_i)}$

Maximum Likelihood / Deviance Estimator:

$\hat{\phi}_{deviance} = \frac{D(\mathbf{y}, \hat{\boldsymbol{\mu}})}{n - p - 1}$

Where $D(\mathbf{y}, \hat{\boldsymbol{\mu}})$ is the residual deviance (see Section 9).

💡 The Pearson $\chi^2$ estimator of $\phi$ is generally preferred for its robustness. For Poisson and Binomial, $\phi = 1$ is known; if the Pearson estimator gives $\hat{\phi} \gg 1$, overdispersion is present (Section 13).

9. Model Fit and Evaluation

9.1 The Deviance

The deviance $D(\mathbf{y}, \hat{\boldsymbol{\mu}})$ is the primary goodness-of-fit measure for GLMs. It is defined as twice the log-likelihood difference between the saturated model (perfect fit, one parameter per observation) and the fitted model:

$D(\mathbf{y}, \hat{\boldsymbol{\mu}}) = 2\left[\ell(\text{saturated}) - \ell(\hat{\boldsymbol{\mu}})\right] = 2\sum_{i=1}^n d_i$

Where the deviance contribution of observation $i$ is:

$d_i = y_i \ln\left(\frac{y_i}{\hat{\mu}_i}\right) - (y_i - \hat{\mu}_i) \quad \text{(for Poisson)}$

Or more generally, twice the contribution of observation $i$ to the log-likelihood difference.

Deviance for each distribution:

9.2 Null Deviance and Residual Deviance

The null deviance $D_0$ is the deviance of the null model (intercept only):

$D_0 = D(\mathbf{y}, \hat{\mu}_0) \quad \text{where } \hat{\mu}_0 = \bar{y}$

The residual deviance $D_r$ is the deviance of the fitted model:

$D_r = D(\mathbf{y}, \hat{\boldsymbol{\mu}})$

The difference $D_0 - D_r$ measures how much the predictors have reduced the unexplained deviance — analogous to the regression sum of squares in linear regression.

For known-dispersion models (Poisson, Binomial), the residual deviance approximately follows $\chi^2_{n-p-1}$ when the model is correct. A residual deviance much larger than $n - p - 1$ suggests poor fit or overdispersion.

9.3 The Pearson $\chi^2$ Statistic

The Pearson $\chi^2$ statistic provides an alternative goodness-of-fit measure:

$X^2 = \sum_{i=1}^n \frac{(y_i - \hat{\mu}_i)^2}{V(\hat{\mu}_i)/w_i}$

Under the correct model with known dispersion, $X^2 \approx \chi^2_{n-p-1}$ for large samples. The Pearson dispersion estimate is $\hat{\phi} = X^2 / (n - p - 1)$.

9.4 Pseudo R² Measures

Since ordinary $R^2$ is not directly meaningful for non-Gaussian GLMs, several pseudo R² measures have been developed:

McFadden's Pseudo R²:

$R^2_{McFadden} = 1 - \frac{\ell(\hat{\boldsymbol{\mu}})}{\ell(\hat{\mu}_0)} = 1 - \frac{D_r + 2\ell_{sat}}{D_0 + 2\ell_{sat}}$

Simplified using deviances:

$R^2_{McFadden} = 1 - \frac{D_r}{D_0}$

Cox-Snell Pseudo R²:

$R^2_{CS} = 1 - \left(\frac{L_0}{L_{\hat{\boldsymbol{\mu}}}}\right)^{2/n}$

Nagelkerke Pseudo R² (scaled to reach maximum of 1):

$R^2_{N} = \frac{R^2_{CS}}{1 - L_0^{2/n}}$

Deviance R² (common in GLM literature):

$R^2_D = 1 - \frac{D_r}{D_0} = R^2_{McFadden}$

Interpretation of McFadden's $R^2$ for GLMs:

9.5 AIC and BIC

Akaike Information Criterion:

$AIC = -2\ell(\hat{\boldsymbol{\mu}}) + 2(p+1)$

Bayesian Information Criterion:

$BIC = -2\ell(\hat{\boldsymbol{\mu}}) + (p+1)\ln(n)$

Where $p + 1$ is the number of estimated regression parameters (including the intercept). For models where $\phi$ is estimated, include it as an additional parameter.

Lower AIC/BIC indicates a better model (adjusted for complexity). AIC favours predictive accuracy; BIC imposes a stronger penalty for complexity and prefers parsimonious models.

⚠️ AIC and BIC require a proper likelihood. They cannot be computed for quasi-GLMs, which use a pseudo-likelihood. For quasi-models, use the F-test for model comparison.

10. Hypothesis Testing and Inference

10.1 Wald Tests for Individual Coefficients

For each coefficient $\beta_j$, the Wald test tests $H_0: \beta_j = 0$:

$z_j = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \sim \mathcal{N}(0, 1) \quad \text{(approximately, for large } n\text{)}$

Two-sided p-value:

$p\text{-value} = 2 \times P(Z > |z_j|) = 2 \times (1 - \Phi(|z_j|))$

A $(1-\alpha) \times 100\%$ Wald confidence interval for $\beta_j$:

$\hat{\beta}_j \pm z_{\alpha/2} \times SE(\hat{\beta}_j)$

For the effect on the original response scale, exponentiate:

• Log link: $e^{\hat{\beta}_j}$ is the rate ratio or mean ratio (Poisson/Gamma).
• Logit link: $e^{\hat{\beta}_j}$ is the odds ratio (Binomial).

Confidence interval on the original scale: $\left[e^{\hat{\beta}_j - z_{\alpha/2} SE(\hat{\beta}_j)}, \; e^{\hat{\beta}_j + z_{\alpha/2} SE(\hat{\beta}_j)}\right]$.

10.2 Likelihood Ratio Test (LRT)

The likelihood ratio test compares two nested models: a smaller (restricted) model $M_0$ and a larger (full) model $M_1$:

$\Lambda = 2\left[\ell(M_1) - \ell(M_0)\right] = D(M_0) - D(M_1)$

Under $H_0$ (the restrictions hold), $\Lambda \sim \chi^2_{df}$ where $df$ is the difference in the number of parameters between $M_1$ and $M_0$.

For testing a single coefficient ($H_0: \beta_j = 0$), $df = 1$. For testing a group of $q$ coefficients jointly, $df = q$.

💡 The LRT is generally preferred over the Wald test for GLMs because it is more accurate in small samples and avoids the Wald test's known deficiencies (e.g., the Hauck-Donner effect, where Wald $z$-values can decrease for very large effects).

10.3 Score Test (Rao Test)

The score test evaluates whether the gradient of the log-likelihood (the score) is significantly different from zero at the restricted (null) parameter values:

$S = \mathbf{s}(\hat{\boldsymbol{\beta}}_0)^T \mathcal{I}^{-1}(\hat{\boldsymbol{\beta}}_0) \mathbf{s}(\hat{\boldsymbol{\beta}}_0) \sim \chi^2_{df}$

The score test only requires fitting the null model (not the full model), making it computationally convenient when the null model is much simpler.

10.4 Analysis of Deviance Table

The analysis of deviance is the GLM analogue of the ANOVA table in linear regression. It sequentially adds predictors and reports the reduction in deviance:

For overdispersed models (quasi-GLMs), use an F-test instead of the $\chi^2$ test, dividing the deviance change by the estimated dispersion $\hat{\phi}$:

$F = \frac{(D_0 - D_r)/p}{\hat{\phi}} \sim F_{p, n-p-1}$

10.5 Confidence Intervals for the Mean Response

A confidence interval for the mean response $\mu_{new}$ at a new predictor vector $\mathbf{x}_{new}$ is constructed on the linear predictor scale (where the asymptotic normality applies) and back-transformed:

Linear predictor and its SE:

$\hat{\eta}_{new} = \mathbf{x}_{new}^T \hat{\boldsymbol{\beta}}, \quad SE(\hat{\eta}_{new}) = \sqrt{\hat{\phi} \cdot \mathbf{x}_{new}^T (\mathbf{X}^T \hat{\mathbf{W}} \mathbf{X})^{-1} \mathbf{x}_{new}}$

Confidence interval on the $\eta$ scale:

$\hat{\eta}_{new} \pm z_{\alpha/2} \times SE(\hat{\eta}_{new})$

Back-transform to the $\mu$ scale using $g^{-1}$:

$\left[g^{-1}\left(\hat{\eta}_{new} - z_{\alpha/2} \times SE(\hat{\eta}_{new})\right), \; g^{-1}\left(\hat{\eta}_{new} + z_{\alpha/2} \times SE(\hat{\eta}_{new})\right)\right]$

💡 Constructing confidence intervals on the link scale and back-transforming (rather than constructing them directly on the $\mu$ scale) ensures the bounds respect the natural constraints of $\mu$ (e.g., positivity for Poisson/Gamma, $[0,1]$ for Binomial).

10.6 Profile Likelihood Confidence Intervals

Profile likelihood confidence intervals are more accurate than Wald intervals, especially in small samples or when the likelihood is asymmetric:

$CI_{profile} = \left\{\beta_j : 2[\ell(\hat{\boldsymbol{\beta}}) - \ell(\hat{\boldsymbol{\beta}}_{(\beta_j)})] \leq \chi^2_{1, \alpha}\right\}$

Where $\hat{\boldsymbol{\beta}}_{(\beta_j)}$ is the MLE with $\beta_j$ fixed at a test value. The DataStatPro application computes both Wald and profile likelihood CIs.

11. Model Diagnostics and Residuals

11.1 Types of GLM Residuals

Unlike linear regression, which has a single natural residual $y_i - \hat{\mu}_i$, GLMs have several types of residuals, each useful for different diagnostic purposes.

11.1.1 Raw (Response) Residuals

$r_i^{raw} = y_i - \hat{\mu}_i$

Simple but not standardised — larger values of $\mu_i$ tend to produce larger raw residuals even if the fit is equally good.

11.1.2 Pearson Residuals

$r_i^P = \frac{y_i - \hat{\mu}_i}{\sqrt{V(\hat{\mu}_i)/w_i}}$

Standardised by the expected standard deviation under the model. Pearson residuals should be approximately $\mathcal{N}(0, 1)$ for large samples if the model is correct.

Standardised Pearson residuals (adjusted for leverage):

$r_i^{PS} = \frac{r_i^P}{\sqrt{\hat{\phi}(1 - h_{ii})}}$

Where $h_{ii}$ is the leverage (hat matrix diagonal). Values $|r_i^{PS}| > 2$ warrant investigation.

11.1.3 Deviance Residuals

$r_i^D = \text{sign}(y_i - \hat{\mu}_i) \sqrt{d_i}$

Where $d_i$ is the deviance contribution of observation $i$ (see Section 9.1). The sum of squared deviance residuals equals the total deviance: $\sum_i (r_i^D)^2 = D_r$.

Deviance residuals are generally preferred for normality assessments because they are closer to normally distributed than Pearson residuals in many GLMs.

Standardised deviance residuals:

$r_i^{DS} = \frac{r_i^D}{\sqrt{\hat{\phi}(1 - h_{ii})}}$

11.1.4 Anscombe Residuals

Anscombe residuals are constructed using a variance-stabilising transformation $A(\mu)$ chosen so that residuals are approximately normally distributed:

$r_i^A = \frac{A(y_i) - A(\hat{\mu}_i)}{A'(\hat{\mu}_i)\sqrt{V(\hat{\mu}_i)/w_i}}$

The Anscombe transformation for each distribution: