PracticeVideosBlogFlashcards

6.1 Deviance, Likelihood Ratio Tests & Residual Diagnostics

Key Takeaways

  • Total deviance D(y, mu_hat) = 2*(l_sat - l_model) measures total discrepancy between the fitted GLM and the hypothetical saturated model that fits each observation perfectly.

  • In models with unknown dispersion (Gamma, Inverse Gaussian, Tweedie), scaled deviance D* = D / phi governs nested model comparison, and the dispersion parameter is estimated via the Pearson statistic phi_hat = chi^2 / (n - p).

  • The Likelihood Ratio Test (LRT) assesses nested models via the deviance difference Delta D / phi ~ chi^2_{Delta df}; for severity GLMs with estimated dispersion, this converts to an F-test.

  • The Bayesian Information Criterion (BIC) imposes a much heavier parameter penalty k*ln(n) than AIC (2k) in large P&C datasets (n > 100,000), preventing overfitting and selecting parsimonious rating structures.

  • Deviance residuals r_{D, i} normalize skewness better than Pearson residuals; residual plots against fitted values detect link misspecification and incorrect variance functions, while leverage and Cook's distance identify influential commercial shock losses.

Last updated: September 2026

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=∑(yi−μ^i)2RSS = \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 (MsatM_{\text{sat}}): A model with as many parameters as observations (p=np = n). It assigns one parameter per data point, achieving an exact fit: μ^i=yi\hat{\mu}_i = y_i for all ii. It represents the highest attainable log-likelihood for the dataset, denoted ℓ(y;y)\ell(\mathbf{y}; \mathbf{y}).
  2. The Current Model (MM): A parsimonious model with pp parameters (p≪np \ll n) yielding fitted values μ^i=g−1(xiTβ^)\hat{\mu}_i = g^{-1}(\mathbf{x}_i^T \hat{\boldsymbol{\beta}}) and log-likelihood ℓ(μ^;y)\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∗(y,μ^)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∗(y,μ^)=2[ℓ(y;y)−ℓ(μ^;y)]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 ii with prior weight ωi\omega_i is:

ℓi(μi;yi)=yiθi−b(θi)ϕ/ωi+c(yi,ϕ/ωi)\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∗(y,μ^)=2ϕ∑i=1nωi[yi(θ~i−θ^i)−(b(θ~i)−b(θ^i))]=D(y,μ^)ϕ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 θ~i\tilde{\theta}_i is the canonical parameter evaluated at yiy_i, and θ^i\hat{\theta}_i is evaluated at μ^i\hat{\mu}_i.

  • D(y,μ^)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(yi,μ^i)d(y_i, \hat{\mu}_i):
D(y,μ^)=∑i=1nωid(yi,μ^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:

DistributionVariance Function V(μ)V(\mu)Dispersion ϕ\phiUnit Deviance d(y,μ)d(y, \mu)
Normal11σ2\sigma^2 (estimated)(y−μ)2(y - \mu)^2
Poissonμ\mu11 (fixed)2[yln⁡(yμ)−(y−μ)]2 \left[ y \ln\left(\frac{y}{\mu}\right) - (y - \mu) \right]
Gammaμ2\mu^2ϕ\phi (estimated)2[−ln⁡(yμ)+y−μμ]2 \left[ -\ln\left(\frac{y}{\mu}\right) + \frac{y - \mu}{\mu} \right]
Inverse Gaussianμ3\mu^3ϕ\phi (estimated)(y−μ)2μ2y\frac{(y - \mu)^2}{\mu^2 y}
Tweedie (1<p<21 < p < 2)μp\mu^pϕ\phi (estimated)2[y2−p(1−p)(2−p)−yμ1−p1−p+μ2−p2−p]2 \left[ \frac{y^{2-p}}{(1-p)(2-p)} - \frac{y\mu^{1-p}}{1-p} + \frac{\mu^{2-p}}{2-p} \right]

Note

For the Poisson unit deviance, when y=0y = 0, the term yln⁡(y/μ)y \ln(y / \mu) evaluates to 00 by taking the continuous limit lim⁡y→0+yln⁡(y)=0\lim_{y \to 0^+} y \ln(y) = 0, yielding d(0,μ)=2μ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 χ2\chi^2, which aggregates standardized squared residuals:

χ2=∑i=1nωi(yi−μ^i)2V(μ^i)\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 (ϕ=1\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:

ϕ^Pearson=χ2n−p=1n−p∑i=1nωi(yi−μ^i)2V(μ^i)\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:

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

Alternatively, software can estimate dispersion using the deviance estimator ϕ^Deviance=D(y,μ^)n−p\hat{\phi}_{\text{Deviance}} = \frac{D(\mathbf{y}, \hat{\boldsymbol{\mu}})}{n - p}. However, casualty actuaries strongly favor the Pearson estimator ϕ^Pearson\hat{\phi}_{\text{Pearson}} because E[χ2]=(n−p)ϕ\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 M0M_0 with p0p_0 parameters can be obtained by constraining p1−p0p_1 - p_0 parameters of the larger (full) model M1M_1 to zero.

  • Null Hypothesis H0H_0: The constrained parameters equal zero (the simpler model M0M_0 is sufficient).
  • Alternative Hypothesis H1H_1: At least one constrained parameter is non-zero (the full model M1M_1 is required).

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

ΔD∗=D0∗−D1∗=D0−D1ϕ    ∼a    χΔdf2,Δdf=p1−p0\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 FF-Test When Dispersion is Estimated

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

F=(D0−D1)/(p1−p0)ϕ^M1=(D0−D1)/ΔdfχM12/(n−p1)    ∼    FΔdf,  n−p1F = \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 FF-statistic exceeds the critical value F1−α(Δdf,n−p1)F_{1-\alpha}(\Delta df, n - p_1), the actuary rejects H0H_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=(β^jSE(β^j))2∼χ12W = \left(\frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)}\right)^2 \sim \chi^2_1 (or z=β^j/SE∼N(0,1)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 M0M_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:

AIC=−2ln⁡(L)+2k\text{AIC} = -2\ln(L) + 2k BIC=−2ln⁡(L)+kln⁡(n)\text{BIC} = -2\ln(L) + k\ln(n)

Where:

  • LL is the maximized likelihood of the model.
  • kk is the total number of estimated parameters (including intercept and dispersion).
  • nn 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.02.0 (constant, independent of sample size).
  • BIC penalty per parameter: ln⁡(n)\ln(n) (grows logarithmically with sample size).
Dataset Size (nn)Insurance ContextAIC Penalty (2k2k)BIC Penalty (kln⁡(n)k \ln(n))Penalty Ratio (BIC / AIC)
1,0001,000Small commercial niche line2.02.0ln⁡(1000)≈6.91\ln(1000) \approx 6.913.5×3.5\times
50,00050,000Regional personal lines book2.02.0ln⁡(50000)≈10.82\ln(50000) \approx 10.825.4×5.4\times
500,000500,000State-wide personal auto2.02.0ln⁡(500000)≈13.12\ln(500000) \approx 13.126.6×6.6\times
5,000,0005,000,000National multi-state book2.02.0ln⁡(5000000)≈15.42\ln(5000000) \approx 15.427.7×7.7\times

The Overfitting Danger of AIC in Casualty Modeling

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

In massive insurance datasets, pure random noise or unrepeatable geographic micro-clusters easily achieve pp-values of 0.010.01 or 0.050.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×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: ri=yi−μ^ir_i = y_i - \hat{\mu}_i.
    • Unstandardized and scale-dependent. Its variance is Var(ri)=ϕωiV(μi)\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:
rP,i=yi−μ^iV(μ^i)/ωir_{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 ∑i=1nrP,i2=χ2\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:
rD,i=sign(yi−μ^i)ωid(yi,μ^i)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 ii's unit deviance contribution. Notice that ∑i=1nrD,i2=D(y,μ^)\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:
rD,istd=rD,iϕ^(1−hii)r_{D, i}^{\text{std}} = \frac{r_{D, i}}{\sqrt{\hat{\phi}(1 - h_{ii})}}
  • Where hiih_{ii} is the leverage of observation ii (the ii-th diagonal element of the GLM hat matrix H=W1/2X(XTWX)−1XTW1/2\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 (η^i\hat{\eta}_i)

  • Ideal Pattern: A uniform, horizontal band centered at zero with constant dispersion roughly between [−2,+2][-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(μ)V(\mu) (e.g., using a Poisson model with V(μ)=μV(\mu) = \mu when loss variance actually scales quadratically as Gamma V(μ)=μ2V(\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< 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 rD,istdr_{D, i}^{\text{std}} against standard normal quantiles N(0,1)\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 (hiih_{ii}) and Cook's Distance (DiD_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:

Di=(rP,istd)2p(hii1−hii)D_i = \frac{(r_{P, i}^{\text{std}})^2}{p} \left( \frac{h_{ii}}{1 - h_{ii}} \right)

Where rP,istdr_{P, i}^{\text{std}} is the standardized Pearson residual and hiih_{ii} is leverage.

  • High Leverage (hii>2p/nh_{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 (Di>4/nD_i > 4 / n or Di>1.0D_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(y,μ^)∼χn−p2D(\mathbf{y}, \hat{\boldsymbol{\mu}}) \sim \chi^2_{n - p}. However, this asymptotic result requires individual observation means μi\mu_i to be large. In insurance frequency data where claim occurrence is rare (e.g., μi=0.05\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 (ΔD∼χΔdf2\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)\ln(Y) is invalid because OLS evaluates the density of ln⁡(Y)\ln(Y), omitting the Jacobian transformation factor 1/y1/y.

Warning

Trap 3: Relying on Raw Residuals to Evaluate Homoscedasticity Novice analysts frequently plot raw residuals yi−μ^iy_i - \hat{\mu}_i against μ^i\hat{\mu}_i and claim the model is flawed because the scatter widens as μ^i\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 rD,istdr_{D, i}^{\text{std}}.

Warning

Trap 4: Misinterpreting Degrees of Freedom in Analysis of Deviance When running an FF-test between nested models, candidates often divide ΔD\Delta D by the total model degrees of freedom p1p_1 instead of the difference in parameters Δdf=p1−p0\Delta df = p_1 - p_0. The numerator in the FF-test must strictly isolate the incremental change in degrees of freedom.

Test Your Knowledge

An actuary compares two nested Gamma severity GLMs: Model 0 (baseline rating variables, p0 = 8 parameters, Deviance = 14,250) and Model 1 (adds 4 credit score tier indicators, p1 = 12 parameters, Deviance = 14,190). The models are estimated on n = 5,012 claims, and Model 1 has a Pearson dispersion estimate phi_hat = 0.75. What is the test statistic for evaluating whether credit score significantly improves severity prediction?

A

Likelihood Ratio Chi-Square statistic = (14,250 - 14,190) = 60.00 with 4 degrees of freedom.

B

Likelihood Ratio Chi-Square statistic = 60.00 / 0.75 = 80.00 with 12 degrees of freedom.

C

F-statistic = [(14,250 - 14,190) / 4] / 0.75 = 20.00, evaluated against an F-distribution with (4, 5,000) degrees of freedom.

D

F-statistic = [(14,250 - 14,190) / 12] / 0.75 = 6.67, evaluated against an F-distribution with (12, 5,000) degrees of freedom.

Test Your Knowledge

A casualty insurer builds a commercial property frequency GLM on a database of 1,200,000 earned policy years. Why does the actuarial team select the Bayesian Information Criterion (BIC) rather than Akaike Information Criterion (AIC) to finalize the rating structure for a regulatory rate filing?

A

BIC mathematically guarantees that all predicted claim frequencies remain strictly positive, whereas AIC can generate negative frequency relativities.

B

AIC cannot be computed for non-canonical link functions, whereas BIC is valid for all link functions.

C

BIC minimizes the mean squared prediction error on holdout data, whereas AIC maximizes the Gini coefficient.

D

Because ln(n) = ln(1,200,000) ≈ 14.0, BIC penalizes each additional parameter seven times more heavily than AIC (penalty of 2.0), preventing the inclusion of noise variables that appear artificially significant due to massive sample size.

Test Your Knowledge

During GLM diagnostic review of a commercial liability severity model, the actuary plots standardized deviance residuals against the linear predictor and observes a distinct bowed (U-shaped) curve. What modeling defect does this diagnostic plot indicate?

A

The link function is misspecified (e.g., using an identity or reciprocal link instead of a log link) or the model is missing quadratic/nonlinear predictor terms.

B

The model has an incorrect dispersion parameter phi that should be re-estimated using the deviance estimator.

C

The dataset contains extreme claim frequency overdispersion that requires a Negative Binomial distribution.

D

The observations were improperly weighted by earned vehicle exposure instead of policy count.

Sections you finish are checked off in the contents.