6.1 Deviance, Likelihood Ratio Tests & Residual Diagnostics

Deviance in Generalized Linear Models

In classical Ordinary Least Squares (OLS) regression, model goodness-of-fit is assessed using the Residual Sum of Squares ($RSS = \sum (y_i - \hat{\mu}_i)^2$). In Generalized Linear Models, where response distributions belong to the Exponential Dispersion Family, the general measure of discrepancy is deviance.

The Saturated Model vs. The Proposed Model

To understand deviance, consider two theoretical extremes of model complexity:

1. The Saturated Model ($M_{\text{sat}}$): A model with as many parameters as observations ($p = n$). It assigns one parameter per data point, achieving an exact fit: $\hat{\mu}_i = y_i$ for all $i$. It represents the highest attainable log-likelihood for the dataset, denoted $\ell(\mathbf{y}; \mathbf{y})$.
2. The Current Model ($M$): A parsimonious model with $p$ parameters ($p \ll n$) yielding fitted values $\hat{\mu}_i = g^{-1}(\mathbf{x}_i^T \hat{\boldsymbol{\beta}})$ and log-likelihood $\ell(\hat{\boldsymbol{\mu}}; \mathbf{y})$.

                               MODEL SPECTRUM

          Null Model                 Current Fitted Model             Saturated Model
         (Intercept Only)               (p parameters)                (n parameters)
       ────────────────────┼──────────────────────────────────────┼────────────────────►
           Worst Fit                    Parsimonious Fit                 Perfect Fit
          Lowest l(mu)                     l(mu_hat)                      Highest l(y)
                                └──────────────────┬───────────────────┘
                                                   ▼
                                        DEVIANCE D = 2 * [l(y) - l(mu_hat)]

Total and Scaled Deviance

The scaled deviance $D^*(\mathbf{y}, \hat{\boldsymbol{\mu}})$ is defined as twice the difference between the maximized log-likelihood of the saturated model and that of the current model:

$D^*(\mathbf{y}, \hat{\boldsymbol{\mu}}) = 2 \Big[ \ell(\mathbf{y}; \mathbf{y}) - \ell(\hat{\boldsymbol{\mu}}; \mathbf{y}) \Big]$

In the Exponential Dispersion Family, the log-likelihood for observation $i$ with prior weight $\omega_i$ is:

$\ell_i(\mu_i; y_i) = \frac{y_i \theta_i - b(\theta_i)}{\phi / \omega_i} + c(y_i, \phi / \omega_i)$

Substituting this into the scaled deviance expression:

$D^*(\mathbf{y}, \hat{\boldsymbol{\mu}}) = \frac{2}{\phi} \sum_{i=1}^n \omega_i \Big[ y_i (\tilde{\theta}_i - \hat{\theta}_i) - \big(b(\tilde{\theta}_i) - b(\hat{\theta}_i)\big) \Big] = \frac{D(\mathbf{y}, \hat{\boldsymbol{\mu}})}{\phi}$

Where $\tilde{\theta}_i$ is the canonical parameter evaluated at $y_i$, and $\hat{\theta}_i$ is evaluated at $\hat{\mu}_i$.

• $D(\mathbf{y}, \hat{\boldsymbol{\mu}})$ is the total (unscaled) deviance.
• $\phi$ is the dispersion parameter.
• The total deviance decomposes into a weighted sum of individual unit deviances $d(y_i, \hat{\mu}_i)$:

$D(\mathbf{y}, \hat{\boldsymbol{\mu}}) = \sum_{i=1}^n \omega_i d(y_i, \hat{\mu}_i)$

Unit Deviance Formulas Across Actuarial Distributions

Each Exponential Dispersion Family distribution has a unique unit deviance function:

[!NOTE] For the Poisson unit deviance, when $y = 0$, the term $y \ln(y / \mu)$ evaluates to $0$ by taking the continuous limit $\lim_{y \to 0^+} y \ln(y) = 0$, yielding $d(0, \mu) = 2\mu$.

---

Pearson Chi-Square Statistic and Dispersion Estimation

An alternative measure of model goodness-of-fit is the Pearson Chi-Square statistic $\chi^2$, which aggregates standardized squared residuals:

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

Estimating the Unknown Dispersion Parameter $\phi$

In Poisson and Binomial models, the dispersion parameter is fixed by definition ($\phi = 1$). However, in severity models (Gamma, Inverse Gaussian) and pure premium models (Tweedie), $\phi$ is an unknown nuisance parameter that must be estimated from data to calculate standard errors and conduct hypothesis tests.

Actuaries estimate $\phi$ using the Pearson estimator:

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

Where:

• $n$ is the number of observations in the training sample.
• $p$ is the number of estimated regression parameters (including the intercept).
• $n - p$ represents the residual degrees of freedom.

Alternatively, software can estimate dispersion using the deviance estimator $\hat{\phi}{\text{Deviance}} = \frac{D(\mathbf{y}, \hat{\boldsymbol{\mu}})}{n - p}$. However, casualty actuaries strongly favor the Pearson estimator $\hat{\phi}{\text{Pearson}}$ because $\mathbb{E}[\chi^2] = (n - p)\phi$ holds exactly under mild moment conditions, whereas the deviance estimator can be severely biased when distributions are highly skewed or sample sizes within rating cells are moderate.

---

Hypothesis Testing and Nested Model Comparison

In actuarial rate filing and pricing research, actuaries systematically evaluate whether adding rating variables (e.g., credit tier, vehicle telematics, roof age) significantly improves model performance.

The Likelihood Ratio Test (LRT) for Nested GLMs

Two models are nested if the smaller (reduced) model $M_0$ with $p_0$ parameters can be obtained by constraining $p_1 - p_0$ parameters of the larger (full) model $M_1$ to zero.

• Null Hypothesis $H_0$: The constrained parameters equal zero (the simpler model $M_0$ is sufficient).
• Alternative Hypothesis $H_1$: At least one constrained parameter is non-zero (the full model $M_1$ is required).

Under Wilks' theorem, the difference in scaled deviances follows an asymptotic Chi-Square distribution:

$\Delta D^* = D^*_0 - D^*_1 = \frac{D_0 - D_1}{\phi} \;\; \stackrel{a}{\sim} \;\; \chi^2_{\Delta df}, \quad \Delta df = p_1 - p_0$

The Actuarial $F$-Test When Dispersion is Estimated

In severity models (Gamma, Inverse Gaussian) where $\phi$ is unknown and must be estimated, using $\Delta D / \hat{\phi}$ against a $\chi^2$ critical value ignores the sampling uncertainty in $\hat{\phi}$. To maintain exact statistical size, actuaries apply the Analysis of Deviance $F$-test:

$F = \frac{(D_0 - D_1) / (p_1 - p_0)}{\hat{\phi}_{M_1}} = \frac{(D_0 - D_1) / \Delta df}{\chi^2_{M_1} / (n - p_1)} \;\; \sim \;\; F_{\Delta df, \; n - p_1}$

If the computed $F$-statistic exceeds the critical value $F_{1-\alpha}(\Delta df, n - p_1)$, the actuary rejects $H_0$ and concludes that the additional rating variables provide statistically significant explanatory power.

                             NESTED MODEL COMPARISON
                                        │
                ┌───────────────────────┴───────────────────────┐
                ▼                                               ▼
        Known Dispersion                                Estimated Dispersion
       (Poisson: phi = 1)                              (Gamma, Inverse Gaussian)
                │                                               │
                ▼                                               ▼
      Likelihood Ratio Test                             Actuarial F-Test
      Delta D ~ Chi-Square(Delta df)                   F = [(D0 - D1)/Delta df] / phi_hat
                                                       F ~ F(Delta df, n - p1)

Wald Tests and Score Tests

• Wald Test: Evaluates individual parameter significance via $W = \left(\frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)}\right)^2 \sim \chi^2_1$ (or $z = \hat{\beta}_j / \text{SE} \sim \mathcal{N}(0, 1)$). Computationally cheap because it requires fitting only the full model, but can suffer from numerical instability when coefficients are very large (the Hauck-Donner effect).
• Score Test (Lagrange Multiplier): Evaluates the slope of the log-likelihood function at the null hypothesis value. Requires fitting only the reduced model $M_0$, making it popular for screening thousands of candidate variables in big data pipelines.

---

Information Criteria: AIC vs. BIC in Massive P&C Datasets

When comparing non-nested models (e.g., comparing a model with Territory $\times$ Age against a model with Territory $\times$ Vehicle Use), deviance tests cannot be applied. Actuaries use penalized likelihood information criteria:

$\text{AIC} = -2\ln(L) + 2k$

$\text{BIC} = -2\ln(L) + k\ln(n)$

Where:

• $L$ is the maximized likelihood of the model.
• $k$ is the total number of estimated parameters (including intercept and dispersion).
• $n$ is the total number of observations (policies or claims).

Why BIC Dominates in Actuarial Rate Filings with Big Data

Consider the parameter penalty terms of AIC and BIC:

• AIC penalty per parameter: $2.0$ (constant, independent of sample size).
• BIC penalty per parameter: $\ln(n)$ (grows logarithmically with sample size).

The Overfitting Danger of AIC in Casualty Modeling

In an insurance portfolio with $2,000,000$ earned policies, $\ln(n) = \ln(2\times 10^6) \approx 14.5$. Under AIC, any variable that increases the log-likelihood by more than $1.0$ unit (corresponding to a trivial $p$-value of $p < 0.157$) will be retained in the model.

In massive insurance datasets, pure random noise or unrepeatable geographic micro-clusters easily achieve $p$-values of $0.01$ or $0.05$ simply due to vast sample size. Using AIC in large P&C datasets results in severe overparameterization, adding dozens of spurious rating levels that fail out-of-time validation and trigger regulatory objections. BIC enforces parsimony, penalizing unnecessary parameters $7\times$ more severely than AIC and producing stable, defensible rate indications.

---

GLM Residual Diagnostics in Actuarial Practice

Because GLM response variables are non-normal and heteroscedastic, classical OLS residual analysis cannot be applied directly. Actuaries rely on specialized GLM residual definitions.

Summary of Actuarial Residual Formulations

1. Raw Residual: $r_i = y_i - \hat{\mu}_i$.
   • Unstandardized and scale-dependent. Its variance is $\text{Var}(r_i) = \frac{\phi}{\omega_i} V(\mu_i)$, which increases directly with the mean. Raw residuals must never be plotted against fitted values to assess homoscedasticity.
2. Pearson Residual:

$r_{P, i} = \frac{y_i - \hat{\mu}_i}{\sqrt{V(\hat{\mu}_i) / \omega_i}}$

• Divides the raw error by the estimated standard deviation. Notice that $\sum_{i=1}^n r_{P, i}^2 = \chi^2$ (the Pearson Chi-Square statistic). However, for skewed distributions (Poisson, Gamma), Pearson residuals remain distinctly skewed.

1. Deviance Residual:

$r_{D, i} = \text{sign}(y_i - \hat{\mu}_i) \sqrt{\omega_i d(y_i, \hat{\mu}_i)}$

• Represents the signed square root of observation $i$'s unit deviance contribution. Notice that $\sum_{i=1}^n r_{D, i}^2 = D(\mathbf{y}, \hat{\boldsymbol{\mu}})$. Deviance residuals transform skewed errors toward approximate symmetry, behaving much more like standard normal errors than Pearson residuals.

1. Standardized (Studentized) Deviance Residual:

$r_{D, i}^{\text{std}} = \frac{r_{D, i}}{\sqrt{\hat{\phi}(1 - h_{ii})}}$

• Where $h_{ii}$ is the leverage of observation $i$ (the $i$-th diagonal element of the GLM hat matrix $\mathbf{H} = \mathbf{W}^{1/2} \mathbf{X} (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W}^{1/2}$). Standardizing adjusts for the fact that high-leverage points exert heavy pull on the regression surface, artificially depressing their unstandardized residuals.

---

Actuarial Diagnostic Plotting Workflows

Actuaries construct four core diagnostic plots to validate GLMs prior to filing rates with insurance regulators:

                               GLM DIAGNOSTIC WORKFLOW
                                          │
            ┌─────────────────────────────┼─────────────────────────────┐
            ▼                             ▼                             ▼
   Residuals vs. Fitted          Residuals vs. Predictors              Q-Q Plot
  Detects Link Misspecification    Detects Non-Linearity /        Checks Distributional
    & Wrong Variance Function          Missing Splines                 Tail Fit
                                          │
                                          ▼
                               Leverage vs. Cook's D
                             Identifies Distorting Shock
                                   Losses / Outliers

1. Standardized Deviance Residuals vs. Linear Predictor ($\hat{\eta}_i$)

• Ideal Pattern: A uniform, horizontal band centered at zero with constant dispersion roughly between $[-2, +2]$.
• Curved (Bowed) Pattern: Indicates link function misspecification (e.g., using an identity link when a log link is required) or omitting significant polynomial terms.
• Funnel / Fan Pattern: Indicates an incorrect variance function $V(\mu)$ (e.g., using a Poisson model with $V(\mu) = \mu$ when loss variance actually scales quadratically as Gamma $V(\mu) = \mu^2$).

2. Residuals vs. Individual Rating Variables

• Plotting residuals against individual continuous variables (such as driver age, vehicle age, or building square footage) reveals uncaptured nonlinearities.
• If residuals for young drivers ($< 25$) are systematically positive while mature drivers are negative, the linear predictor fails to capture the youthful risk curve, signaling the need for piece-wise splines or category re-binning.

3. Normal Q-Q Plot of Standardized Deviance Residuals

• Compares empirical quantiles of $r_{D, i}^{\text{std}}$ against standard normal quantiles $\mathcal{N}(0, 1)$.
• In severity modeling, Q-Q plots reveal right-tail inadequacies: if points drift steeply above the 45-degree reference line in the upper right quadrant, empirical claims exceed Gamma tail expectations, warranting an Inverse Gaussian or Pareto distribution.

4. Leverage ($h_{ii}$) and Cook's Distance ($D_i$)

In commercial property and casualty lines, a single $15,000,000 explosion loss or massive multi-vehicle catastrophe can distort an entire territory or class relativity. Actuaries isolate influential claims using Cook's Distance for GLMs:

$D_i = \frac{(r_{P, i}^{\text{std}})^2}{p} \left( \frac{h_{ii}}{1 - h_{ii}} \right)$

Where $r_{P, i}^{\text{std}}$ is the standardized Pearson residual and $h_{ii}$ is leverage.

• High Leverage ($h_{ii} > 2p / n$): Identifies policies with extreme or unusual combinations of predictor variables (e.g., a 16-year-old driver operating a $250,000 commercial dump truck).
• High Influence ($D_i > 4 / n$ or $D_i > 1.0$): Identifies observations that single-handedly alter coefficient estimates. In actuarial practice, high-influence claims are not silently deleted; they are scrutinized for data errors, capped at basic limits, or assigned to excess/catastrophe layers.

---

Actuarial Exam Traps & Pitfalls

[!WARNING] Trap 1: Using Total Deviance as an Absolute Goodness-of-Fit Test in Sparse Data Textbooks often state that $D(\mathbf{y}, \hat{\boldsymbol{\mu}}) \sim \chi^2_{n - p}$. However, this asymptotic result requires individual observation means $\mu_i$ to be large. In insurance frequency data where claim occurrence is rare (e.g., $\mu_i = 0.05$), the chi-square approximation for absolute deviance completely fails. However, the Likelihood Ratio Test for the difference between two nested models ($\Delta D \sim \chi^2_{\Delta df}$) remains highly accurate.

[!WARNING] Trap 2: Comparing AIC/BIC Across Different Distributions or Data Subsets Information criteria are valid only when comparing models fit to the exact same dataset with the exact same definition of the likelihood (including all normalizing constants). Comparing the AIC of a Gamma GLM against an OLS regression on $\ln(Y)$ is invalid because OLS evaluates the density of $\ln(Y)$, omitting the Jacobian transformation factor $1/y$.

[!WARNING] Trap 3: Relying on Raw Residuals to Evaluate Homoscedasticity Novice analysts frequently plot raw residuals $y_i - \hat{\mu}_i$ against $\hat{\mu}_i$ and claim the model is flawed because the scatter widens as $\hat{\mu}i$ increases. In insurance GLMs (Poisson, Gamma, Inverse Gaussian), variance is theoretically defined to increase with the mean! Homoscedasticity must be diagnosed using standardized deviance residuals $r{D, i}^{\text{std}}$.

[!WARNING] Trap 4: Misinterpreting Degrees of Freedom in Analysis of Deviance When running an $F$-test between nested models, candidates often divide $\Delta D$ by the total model degrees of freedom $p_1$ instead of the difference in parameters $\Delta df = p_1 - p_0$. The numerator in the $F$-test must strictly isolate the incremental change in degrees of freedom.