df = pd.read_parquet("data/processed/motor_model.parquet")
bins = [17, 21, 25, 30, 40, 50, 65, 75, 120]
labels = ["17-20", "21-24", "25-29", "30-39",
          "40-49", "50-64", "65-74", "75+"]
df["age_band"] = pd.cut(df["driver_age"], bins,
                        right=False, labels=labels)
df["log_miles"] = np.log(df["annual_mileage"] / 12_000)

train = df[df["policy_year"] <= 2024]   # fit on these
hold = df[df["policy_year"] == 2025]    # keep in the vault

train["n_liab"].value_counts()
# 0: 62,293   1: 3,579   2: 213   3: 12

One row per vehicle per policy term. exposure is car-years earned: a policy canceled after three months earns about 0.25.

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

$$\mathrm{E}[Y] = \mu = b'(\theta) \qquad \mathrm{Var}(Y) = \phi\, V(\mu)$$

Each family is defined by its variance function V(μ): how the spread grows with the mean. φ is the dispersion.

$$g(\mu_i) = \eta_i = \mathbf{x}_i^\top\beta$$

• Random component: Y_i from an EDF distribution with mean μ_i
• Systematic component: the linear predictor η_i
• Link: g maps the mean’s range onto the whole real line

With the log link, μ = e^η is always positive and every coefficient is a multiplicative relativity:

$$\mu_i = e^{\beta_0}\,e^{\beta_1 x_{i1}}\cdots e^{\beta_p x_{ip}}$$

$$P(Y_i = y) = \frac{e^{-\mu_i}\mu_i^{\,y}}{y!} \qquad \ln\mu_i = \mathbf{x}_i^\top\beta$$

Log-likelihood and its gradient (the score):

$$\ell(\beta) = \sum_i \big(y_i\,\mathbf{x}_i^\top\beta - e^{\mathbf{x}_i^\top\beta} - \ln y_i!\big)$$

$$\frac{\partial\ell}{\partial\beta} = \sum_i (y_i - \mu_i)\,\mathbf{x}_i = \mathbf{X}^\top(y-\mu) = 0$$

Solved by iteratively reweighted least squares (IRLS): usually 5–8 iterations.

A vehicle insured for half a year should expect half the claims:

$$\mu_i = \text{Exposure}_i \times e^{\mathbf{x}_i^\top\beta} \iff \ln\mu_i = \ln(\text{Exposure}_i) + \mathbf{x}_i^\top\beta$$

ln(Exposure) enters the linear predictor with its coefficient fixed at 1. It’s an offset, not a feature.

import statsmodels.api as sm
import statsmodels.formula.api as smf

pois = smf.glm("n_liab ~ " + rhs, data=train,
               family=sm.families.Poisson(),
               offset=np.log(train["exposure"])).fit()

# identical: statsmodels takes the log for you
smf.glm(..., exposure=train["exposure"])

With the offset, e^β is a ratio of claims per car-year, the frequency, however long each policy was on the books.

rhs = ("C(vehicle_style, Treatment('Sedan'))"
       " + C(age_band, Treatment('40-49'))"
       " + urban + log_miles"
       " + C(channel, Treatment('Agent'))"
       " + holder_gender")
pois = smf.glm("n_liab ~ " + rhs, data=train,
               family=sm.families.Poisson(),
               offset=np.log(train["exposure"])).fit()

irr = np.exp(pois.params)
ci = np.exp(pois.conf_int())
pois.fittedvalues.sum()      # 4,041.00 = actual claims

Deviance 22,430.1 on 66,079 df · AIC 30,216.8

Likelihood-ratio test: drop channel and gender

lr = 2 * (full.llf - reduced.llf)      # 1.703
df_diff = reduced.df_resid - full.df_resid   # 3
scipy.stats.chi2.sf(lr, df_diff)       # p = 0.636
full.aic, reduced.aic                  # 30,216.8 → 30,212.5

No evidence that channel or gender affects frequency, and the truth agrees: neither was used to simulate claims. Drop them.

Age as a straight line?

Replace the bands with driver_age: one coefficient, IRR 0.925 per decade — “older is safer”, which is wrong for everyone over 70.

$$\text{NB2:}\quad \mathrm{Var}(Y) = \mu + \alpha\mu^2$$

Gamma-distributed rate multipliers with variance α. α = 0 is the Poisson.

# Diagnose: Pearson chi-square per degree of freedom
pois.pearson_chi2 / pois.df_resid      # 1.049

# Estimate alpha by maximum likelihood
nb = smf.negativebinomial("n_liab ~ " + rhs2, data=train,
        offset=np.log(train["exposure"])).fit()
nb.params["alpha"]                     # 0.649 (se 0.110)
2 * (nb.llf - pois.llf)                # LR = 50.4

• Poisson estimates stay consistent under overdispersion: the mean model is still right
• Its standard errors are too small, so its p-values are too optimistic
• Here the difference is 2–4%. On health visit counts (variance 3–6× the mean) it would be huge