{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-000",
   "metadata": {},
   "source": [
    "# Lesson 2 — GLMs for Claim Frequency\n",
    "\n",
    "**Companion notebook** for the Lesson 2 slides. You'll build a Poisson frequency model\n",
    "for motor liability claims with an exposure offset, select factors with likelihood-ratio\n",
    "tests, diagnose overdispersion, and read incidence rate ratios (IRRs).\n",
    "\n",
    "Data: `../data/processed/motor_model.parquet`, one row per vehicle per policy term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-001",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check(label, got, expected, tol=1e-6):\n",
    "    \"\"\"Compare your answer with the expected value; prints ✅ or ❌ (never raises).\"\"\"\n",
    "    import math\n",
    "    try:\n",
    "        if isinstance(expected, float) or isinstance(got, float):\n",
    "            ok = math.isclose(float(got), float(expected), rel_tol=tol, abs_tol=tol)\n",
    "        elif isinstance(expected, (list, tuple)) and expected and isinstance(expected[0], str):\n",
    "            ok = list(got) == list(expected)\n",
    "        elif hasattr(expected, \"__len__\") and not isinstance(expected, str):\n",
    "            import numpy as np\n",
    "            ok = np.allclose(np.asarray(got, dtype=float), np.asarray(expected, dtype=float), rtol=tol, atol=tol)\n",
    "        else:\n",
    "            ok = got == expected\n",
    "    except Exception:\n",
    "        ok = False\n",
    "    print((\"✅ \" if ok else \"❌ \") + f\"{label}: got {got!r}, expected {expected!r}\")\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as st\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "df = pd.read_parquet(\"../data/processed/motor_model.parquet\")\n",
    "bins = [17, 21, 25, 30, 40, 50, 65, 75, 120]\n",
    "labels = [\"17-20\", \"21-24\", \"25-29\", \"30-39\", \"40-49\", \"50-64\", \"65-74\", \"75+\"]\n",
    "df[\"age_band\"] = pd.cut(df[\"driver_age\"], bins, right=False, labels=labels)\n",
    "df[\"log_miles\"] = np.log(df[\"annual_mileage\"] / 12_000)\n",
    "\n",
    "train = df[df[\"policy_year\"] <= 2024].copy()   # fit on these\n",
    "hold = df[df[\"policy_year\"] == 2025].copy()    # the vault: untouched until Lesson 7\n",
    "print(f\"training: {len(train):,} vehicle-terms, {train['exposure'].sum():,.1f} car-years, \"\n",
    "      f\"{train['n_liab'].sum():,} claims\")\n",
    "train[\"n_liab\"].value_counts().sort_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-002",
   "metadata": {},
   "source": [
    "---\n",
    "## Your turn 1 — one-way frequency by age band\n",
    "\n",
    "Frequency = claims ÷ car-years, computed as a **ratio of sums** within each band.\n",
    "Store the result as a Series `freq_by_age` indexed by band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-003",
   "metadata": {},
   "outputs": [],
   "source": [
    "freq_by_age = None\n",
    "...  # your code here\n",
    "check(\"frequency 17-20\", None if freq_by_age is None else round(freq_by_age[\"17-20\"], 4), 0.1162)\n",
    "check(\"frequency 50-64\", None if freq_by_age is None else round(freq_by_age[\"50-64\"], 4), 0.0520)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-004",
   "metadata": {},
   "source": [
    "## 1. The Poisson GLM with an exposure offset\n",
    "\n",
    "`offset=np.log(exposure)` fixes the coefficient of log-exposure at 1, so the model\n",
    "describes claims **per car-year**. `exposure=` does the same thing for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-005",
   "metadata": {},
   "outputs": [],
   "source": [
    "rhs_full = (\"C(vehicle_style, Treatment('Sedan')) + C(age_band, Treatment('40-49'))\"\n",
    "            \" + urban + log_miles + C(channel, Treatment('Agent')) + holder_gender\")\n",
    "full = smf.glm(\"n_liab ~ \" + rhs_full, data=train, family=sm.families.Poisson(),\n",
    "               offset=np.log(train[\"exposure\"])).fit()\n",
    "same = smf.glm(\"n_liab ~ \" + rhs_full, data=train, family=sm.families.Poisson(),\n",
    "               exposure=train[\"exposure\"]).fit()\n",
    "print(\"offset vs exposure=, largest difference:\", float(np.abs(full.params - same.params).max()))\n",
    "print(f\"deviance {full.deviance:,.1f} on {full.df_resid:,.0f} df; AIC {full.aic:,.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-006",
   "metadata": {},
   "source": [
    "## Your turn 2 — IRRs and the balance property\n",
    "\n",
    "1. Build a table of IRRs (`exp` of the coefficients) with 95% limits and p-values.\n",
    "2. Check the balance property: the sum of fitted claims equals the actual claims.\n",
    "\n",
    "Store the urban IRR and the total fitted claims."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-007",
   "metadata": {},
   "outputs": [],
   "source": [
    "urban_irr = fitted_total = None\n",
    "...  # your code here\n",
    "check(\"urban IRR\", None if urban_irr is None else round(urban_irr, 3), 1.366)\n",
    "check(\"fitted claims\", None if fitted_total is None else round(fitted_total, 2), 4041.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-008",
   "metadata": {},
   "source": [
    "## Your turn 3 — do channel and gender belong?\n",
    "\n",
    "Fit the reduced model without channel and gender, then compute the likelihood-ratio\n",
    "statistic `2 × (llf_full − llf_reduced)` and its p-value on the difference in degrees of\n",
    "freedom. Keep the reduced model as `reduced`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-009",
   "metadata": {},
   "outputs": [],
   "source": [
    "rhs = (\"C(vehicle_style, Treatment('Sedan')) + C(age_band, Treatment('40-49'))\"\n",
    "       \" + urban + log_miles\")\n",
    "reduced = lr = p_value = None\n",
    "...  # your code here\n",
    "check(\"LR statistic\", None if lr is None else round(lr, 3), 1.703)\n",
    "check(\"p-value\", None if p_value is None else round(p_value, 3), 0.636)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-010",
   "metadata": {},
   "source": [
    "## 2. Age as a straight line misses the U\n",
    "\n",
    "A GLM is linear in the predictor. With `driver_age` as one number, the model can only\n",
    "say \"older is safer\" (or riskier) everywhere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-011",
   "metadata": {},
   "outputs": [],
   "source": [
    "linear = smf.glm(\"n_liab ~ C(vehicle_style, Treatment('Sedan')) + driver_age + urban + log_miles\",\n",
    "                 data=train, family=sm.families.Poisson(), offset=np.log(train[\"exposure\"])).fit()\n",
    "print(f\"IRR per decade {np.exp(10 * linear.params['driver_age']):.4f}; AIC linear {linear.aic:,.1f}\")\n",
    "if reduced is not None:\n",
    "    print(f\"AIC with bands {reduced.aic:,.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-012",
   "metadata": {},
   "source": [
    "## 3. Overdispersion\n",
    "\n",
    "Pearson χ² per degree of freedom near 1 is consistent with Poisson. The negative\n",
    "binomial (NB2) adds α: Var(Y) = μ + αμ². Health visit counts are a clearer example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-013",
   "metadata": {},
   "outputs": [],
   "source": [
    "enc = pd.read_parquet(\"../data/health/encounters.parquet\")\n",
    "er = (enc[(enc[\"service_date\"].dt.year == 2024) & (enc[\"service_category\"] == \"Emergency\")]\n",
    "      .groupby(\"member_id\").size()\n",
    "      .reindex(pd.read_csv(\"../data/health/members.csv\")[\"member_id\"], fill_value=0))\n",
    "print(f\"ER visits per member: mean {er.mean():.3f}, variance {er.var():.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-014",
   "metadata": {},
   "source": [
    "## Your turn 4 — a negative binomial frequency model\n",
    "\n",
    "1. Compute Pearson χ²/df for the reduced Poisson model.\n",
    "2. Fit `smf.negativebinomial` with the same formula and offset; store α.\n",
    "3. Compare the IRRs and standard errors with the Poisson model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-015",
   "metadata": {},
   "outputs": [],
   "source": [
    "dispersion = alpha = None\n",
    "if reduced is not None:\n",
    "    ...  # your code here\n",
    "check(\"Pearson chi2 / df\", None if dispersion is None else round(dispersion, 3), 1.049)\n",
    "check(\"NB alpha\", None if alpha is None else round(alpha, 3), 0.649)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-016",
   "metadata": {},
   "source": [
    "## Your turn 5 — read the model like a rating manual\n",
    "\n",
    "Using the reduced Poisson model:\n",
    "1. What is the mileage **elasticity** (the `log_miles` coefficient)?\n",
    "2. What multiple of the base rate does a 19-year-old in an urban sports car driving\n",
    "   15,000 miles a year get? Multiply the relevant IRRs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-017",
   "metadata": {},
   "outputs": [],
   "source": [
    "elasticity = multiple = None\n",
    "if reduced is not None:\n",
    "    ...  # your code here\n",
    "check(\"mileage elasticity\", None if elasticity is None else round(elasticity, 3), 0.488)\n",
    "check(\"multiple of the base rate\", None if multiple is None else round(multiple, 2), 5.42)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-018",
   "metadata": {},
   "source": [
    "## 4. Compare with the truth\n",
    "\n",
    "The book was simulated with Sports 1.40, Minivan 0.85, urban 1.30 and a mileage\n",
    "elasticity of 0.5. Sports comes out far higher. Why? The data contain an extra 60% for\n",
    "young drivers **in sports cars**, and the model has no term for it. Lesson 7 finds it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-019",
   "metadata": {},
   "outputs": [],
   "source": [
    "if reduced is not None:\n",
    "    display(np.exp(reduced.params).round(3))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}