insurance-dispersion 0.1.2

Double GLM (DGLM) for joint modelling of mean and dispersion in insurance pricing

insurance-dispersion

Double GLM (DGLM) for joint modelling of mean and dispersion in non-life insurance pricing.

The problem

Standard GLMs assume a single scalar dispersion parameter phi shared across all observations. For a Gamma severity model, that means your fleet broker policy and your personal lines online policy are assumed to have identical volatility around the fitted mean. That assumption is almost always wrong.

Dispersion varies systematically with the same risk factors that drive the mean — and often with different factors entirely. Broker-sourced business tends to be more volatile (larger phi) because brokers aggregate heterogeneous risks. Fleet accounts have more predictable frequencies. High-limit policies show fat-tailed severity that the Gamma captures poorly with a flat dispersion assumption.

The Double GLM (Smyth 1989) solves this by adding a second regression model for phi:

Mean submodel:        g(mu_i)  = x_i^T beta    [standard GLM]
Dispersion submodel:  h(phi_i) = z_i^T alpha   [new: each obs gets its own phi]

Var[Y_i] = phi_i * V(mu_i)

This matters for:

• Risk-differentiated pricing: your pure premium estimate is mu_i, but quoting confidence depends on phi_i
• Reinsurance pricing: the tail risk on a policy is driven by both mu_i and phi_i
• Model validation: a well-specified mean model with poor dispersion fit still mispredicts volatility
• Credibility: low-phi risks can be priced more confidently than high-phi risks

Installation

pip install insurance-dispersion

Or from source:

git clone https://github.com/burning-cost/insurance-dispersion
cd insurance-dispersion
uv pip install -e .

Quick start

import numpy as np
import pandas as pd
from insurance_dispersion import DGLM
import insurance_dispersion.families as fam

# Synthetic claim severity data
rng = np.random.default_rng(42)
n = 500
df = pd.DataFrame({
    "vehicle_class":  rng.choice(["A", "B", "C"], size=n),
    "age_band":       rng.choice(["17-24", "25-35", "36-60"], size=n),
    "vehicle_value":  rng.uniform(5000, 40000, size=n),
    "channel":        rng.choice(["direct", "broker"], size=n),
    "limit_band":     rng.choice(["50k", "100k", "250k"], size=n),
    "earned_premium": rng.uniform(0.5, 1.0, size=n),
})
df["claim_amount"] = rng.gamma(shape=2.0, scale=1500.0, size=n)

# Fit a Gamma DGLM for claim severity
# Mean model: severity depends on vehicle class and age band
# Dispersion model: volatility depends on distribution channel and limit band
model = DGLM(
    formula="claim_amount ~ C(vehicle_class) + C(age_band) + log(vehicle_value)",
    dformula="~ C(channel) + C(limit_band)",
    family=fam.Gamma(),
    data=df,
    exposure="earned_premium",  # log-offset in mean only
    method="reml",              # REML correction (recommended)
)

result = model.fit()
print(result.summary())

Output:

Double GLM (DGLM) Results
============================================================
Family:      Gamma(link='log')
Method:      REML
Observations: 500
Converged:   True (after 8 iterations)
Log-lik:     -4182.3521
AIC:         8398.7042

Mean Submodel Coefficients:
------------------------------------------------------------
                            coef  exp_coef    se       z  p_value
Intercept               2.1543    8.6224  0.0321  67.12    0.0000
C(vehicle_class)[T.B]   0.1823    1.1999  0.0211   8.64    0.0000
...

Dispersion Submodel Coefficients:
------------------------------------------------------------
                          coef  exp_coef    se       z  p_value
Intercept             -0.8234    0.4390  0.0412 -19.99    0.0000
C(channel)[T.broker]   0.6112    1.8426  0.0518  11.80    0.0000
...

Factor tables

# Mean relativities: exp(beta) for each level vs. base
mean_rel = result.mean_relativities()
print(mean_rel[["exp_coef", "se", "p_value"]])

# Dispersion relativities: exp(alpha)
# Broker channel has 1.84x the dispersion of direct channel
disp_rel = result.dispersion_relativities()
print(disp_rel[["exp_coef", "se", "p_value"]])

Predictions

new_risk = pd.DataFrame({
    "vehicle_class": ["A", "B"],
    "age_band": ["25-35", "17-24"],
    "vehicle_value": [15000, 8000],
    "channel": ["direct", "broker"],
    "limit_band": ["100k", "50k"],
    "earned_premium": [1.0, 1.0],
})

# Expected severity
mu_pred = result.predict(new_risk, which="mean")

# Observation-level dispersion
phi_pred = result.predict(new_risk, which="dispersion")

# Predicted variance = phi_i * V(mu_i) = phi_i * mu_i^2 (Gamma)
var_pred = result.predict(new_risk, which="variance")

Overdispersion test

# Likelihood ratio test: constant phi vs. phi = f(channel, limit_band)
# Note: use method='ml' for formal testing — REML-based LRT is approximate
test = result.overdispersion_test()
print(f"LRT statistic: {test['statistic']:.2f}")
print(f"df: {test['df']}")
print(f"p-value: {test['p_value']:.4f}")
print(test["conclusion"])

Diagnostics

from insurance_dispersion import diagnostics

# Residuals
pearson_r = diagnostics.pearson_residuals(result)
deviance_r = diagnostics.deviance_residuals(result)
qr = diagnostics.quantile_residuals(result)  # ~ N(0,1) under true model

# QQ plot data
qq = diagnostics.qq_plot_data(result)
import matplotlib.pyplot as plt
plt.scatter(qq["theoretical"], qq["observed"], alpha=0.3, s=10)
plt.plot([-3, 3], [-3, 3], "r--")
plt.xlabel("N(0,1) quantiles")
plt.ylabel("Observed quantile residuals")

# Dispersion diagnostic
diag = diagnostics.dispersion_diagnostic(result)
plt.scatter(diag["fitted_phi"], diag["scaled_deviance"], alpha=0.2, s=8)
plt.axhline(1.0, color="red", linestyle="--")  # E[delta_i] = 1 under model
plt.xlabel("Fitted phi")
plt.ylabel("Scaled unit deviance")

Supported families

Family	Default link	Use case
Gamma()	log	Claim severity
InverseGaussian()	log	Heavy-tail severity
Tweedie(p=1.5)	log	Pure premium (compound Poisson-Gamma)
Gaussian()	identity	Reserve amounts, Gaussian responses
Poisson()	log	Claim frequency (extra-Poisson variation)
NegativeBinomial(alpha=1.0)	log	Overdispersed frequency

Algorithm

Alternating IRLS (Smyth 1989, Smyth & Verbyla 1999):

1. Initialise mu from intercept-only GLM, phi = 1
2. Mean step: IRLS for GLM(y ~ X, family, weights = prior_weights / phi_i)
3. Dispersion step: compute unit deviances d_i; fit Gamma GLM on delta_i = d_i with log link (E[d_i] = phi_i)
4. REML correction (method='reml'): subtract hat-matrix diagonal from delta_i before dispersion fit. Recommended when the mean model has many parameters.
5. Check convergence: relative change in Gamma deviance of the dispersion pseudo-response < epsilon
6. Repeat until convergence or maxit reached

Pure numpy/scipy. No ML frameworks, no statsmodels dependency.

Design choices

formulaic not patsy: patsy is unmaintained. formulaic has an active development community, cleaner model matrix schemas for prediction on new data, and better handling of interactions and transformations.

method='reml' default: the REML correction removes the contribution of estimating beta from the dispersion score. With even 10 mean parameters in a dataset of 500 observations this makes a material difference to the dispersion estimates. The correction is cheap (hat diagonal via QR) and almost always helps.

Exposure on mean only: log(exposure) enters as an offset in the mean linear predictor. Dispersion phi_i is per-unit-exposure: a 6-month policy has the same dispersion per claim as a 12-month policy with identical risk characteristics.

Log link for dispersion (only supported link): ensures phi_i > 0 always. The identity link option has been removed — it was accepted but not implemented in the fitting engine, silently defaulting to log link. Only log link is supported.

Databricks Notebook

A ready-to-run Databricks notebook benchmarking this library against standard approaches is available in burning-cost-examples.

Reference

• Smyth (1989): "Generalized Linear Models with Varying Dispersion", JRSS-B 51:47-60
• Smyth & Verbyla (1999): "Adjusted likelihood methods for modelling dispersion in GLMs", Environmetrics 10:695-709
• R dglm package: https://github.com/cran/dglm

Performance

Benchmarked against a constant-phi Tweedie GLM (statsmodels) on synthetic UK commercial property pure premium data: 25,000 policies, known DGP where phi varies 3–6x across distribution channels (direct vs broker SME vs broker large), temporal 70/30 train/test split. See notebooks/benchmark_dispersion.py for full methodology.

Metric	Tweedie GLM (const phi)	DGLM
Tweedie deviance (test)	1.842	1.847
Phi MAE vs true	0.312	0.089
Max channel A/E deviation	0.38	0.11
Variance ratio by channel	0.51–1.94	0.87–1.12
Overdispersion LRT p-value	not applicable	< 0.001
Fit time	1.2s	5.8s

The DGLM's Tweedie deviance on the test set is comparable to the constant-phi GLM (the mean submodel is essentially the same). The material improvement is in variance calibration: the constant-phi GLM assigns variance ratios of 0.51–1.94 across channels (substantially miscalibrated), while the DGLM achieves 0.87–1.12. Phi MAE drops from 0.312 to 0.089. This matters for reinsurance pricing and capital loading where the volatility estimate, not just the mean, drives the price.

The fit time penalty (3–6x slower) comes from the alternating IRLS. On 25k policies this is under 10 seconds; the cost is proportional to the number of iterations needed for convergence.

Gamma severity benchmark (unit test, benchmarks/benchmark.py)

3,000 Gamma severity observations, DGP with known channel-dependent dispersion (phi_direct=0.30, phi_broker=1.20). 2,400 train / 600 test. Converged in 8 iterations.

Metric	DGLM result	Target
Converged	Yes (8 iter)	—
Phi estimate (direct channel)	0.304	0.30
Phi estimate (broker channel)	1.446	1.20
90% PI coverage — all	88.8%	90.0%
90% PI coverage — direct	88.5%	90.0%
90% PI coverage — broker	89.4%	90.0%

Phi recovery is close for direct (0.304 vs true 0.30) and moderately close for broker (1.446 vs true 1.20). The broker estimate is upward-biased at n=957 training observations — expected given the fat-tailed Gamma DGP. Coverage across both channels is within 1.5pp of nominal.

Related Libraries

Library	What it does
insurance-distributional-glm	GAMLSS — the full RS algorithm for jointly modelling mean and all distributional parameters including shape
insurance-frequency-severity	Joint frequency-severity models with Sarmanov copula — extends dispersion modelling to the two-part structure