{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-000",
   "metadata": {},
   "source": [
    "# Lesson 7 — Trees, Forests & Gradient Boosting\n",
    "\n",
    "**Companion notebook** for the Lesson 7 slides. You'll grow a Poisson regression tree,\n",
    "train XGBoost with a Poisson objective and an exposure offset, compare it fairly with\n",
    "GLMs, and explain it with SHAP values."
   ]
  },
  {
   "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 matplotlib.pyplot as plt\n",
    "import shap\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import xgboost as xgb\n",
    "from sklearn.metrics import mean_poisson_deviance\n",
    "from sklearn.model_selection import GroupKFold\n",
    "from sklearn.tree import DecisionTreeRegressor, export_text\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",
    "df[\"young_sports\"] = (df[\"driver_age\"] < 25) & (df[\"vehicle_style\"] == \"Sports\")\n",
    "train = df[df[\"policy_year\"] <= 2024].copy()\n",
    "hold = df[df[\"policy_year\"] == 2025].copy()\n",
    "\n",
    "FEATURES = [\"driver_age\", \"years_licensed\", \"n_drivers\", \"vehicle_age\", \"vehicle_value\",\n",
    "            \"annual_mileage\", \"urban\", \"coast_km\", \"latitude\", \"longitude\", \"vehicle_style\",\n",
    "            \"channel\", \"holder_gender\"]\n",
    "\n",
    "\n",
    "def design(d):\n",
    "    \"\"\"Numeric feature matrix for tree models: 0/1 flags and one-hot categories.\"\"\"\n",
    "    X = d[FEATURES].copy()\n",
    "    X[\"urban\"] = X[\"urban\"].astype(int)\n",
    "    X[\"male\"] = (X.pop(\"holder_gender\") == \"M\").astype(int)\n",
    "    X = pd.get_dummies(X, columns=[\"vehicle_style\", \"channel\"], dtype=int)\n",
    "    X.columns = [c.replace(\"vehicle_style_\", \"style_\") for c in X.columns]\n",
    "    return X\n",
    "\n",
    "\n",
    "Xtr = design(train)\n",
    "Xho = design(hold).reindex(columns=Xtr.columns, fill_value=0)\n",
    "ytr, etr = train[\"n_liab\"].to_numpy(), train[\"exposure\"].to_numpy()\n",
    "yho, eho = hold[\"n_liab\"].to_numpy(), hold[\"exposure\"].to_numpy()\n",
    "list(Xtr.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-002",
   "metadata": {},
   "source": [
    "## 1. A GLM with a hand-made interaction\n",
    "\n",
    "If you *know* young drivers in sports cars are extra risky, you can add the term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-003",
   "metadata": {},
   "outputs": [],
   "source": [
    "GLM = (\"C(vehicle_style, Treatment('Sedan')) + C(age_band, Treatment('40-49'))\"\n",
    "       \" + urban + log_miles\")\n",
    "glm = smf.glm(\"n_liab ~ \" + GLM, data=train, family=sm.families.Poisson(),\n",
    "              offset=np.log(train[\"exposure\"])).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-004",
   "metadata": {},
   "source": [
    "---\n",
    "## Your turn 1 — test the interaction\n",
    "\n",
    "Fit the GLM again with `+ young_sports`. Store the interaction's IRR and the AIC of the\n",
    "new model. (The truth used to simulate the data was 1.60.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-005",
   "metadata": {},
   "outputs": [],
   "source": [
    "young_sports_irr = aic_int = None\n",
    "...  # your code here\n",
    "check(\"young x sports IRR\", None if young_sports_irr is None else round(young_sports_irr, 3), 1.407)\n",
    "check(\"AIC with the interaction\", None if aic_int is None else round(aic_int, 1), 30203.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-006",
   "metadata": {},
   "source": [
    "## 2. One regression tree\n",
    "\n",
    "`criterion=\"poisson\"` splits to reduce Poisson deviance. Fit the claim **rate** with\n",
    "exposure as the sample weight: the tree equivalent of an offset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-007",
   "metadata": {},
   "outputs": [],
   "source": [
    "cart = DecisionTreeRegressor(criterion=\"poisson\", max_depth=3, min_samples_leaf=500, random_state=0)\n",
    "cart.fit(Xtr, ytr / etr, sample_weight=etr)\n",
    "print(export_text(cart, feature_names=list(Xtr.columns), decimals=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-008",
   "metadata": {},
   "source": [
    "## Your turn 2 — XGBoost with a Poisson objective\n",
    "\n",
    "Train on policy years 2021–2023 and early-stop on 2024:\n",
    "`base_margin=np.log(exposure)` is the offset. Use the parameters below with\n",
    "`num_boost_round=3000` and `early_stopping_rounds=100`. Store the number of trees\n",
    "(`best_iteration + 1`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-009",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {\"objective\": \"count:poisson\", \"eta\": 0.03, \"max_depth\": 3, \"min_child_weight\": 50,\n",
    "          \"subsample\": 0.8, \"colsample_bytree\": 0.8, \"lambda\": 1.0,\n",
    "          \"eval_metric\": \"poisson-nloglik\", \"seed\": 0}\n",
    "val = (train[\"policy_year\"] == 2024).to_numpy()\n",
    "n_rounds = None\n",
    "...  # your code here\n",
    "check(\"trees after early stopping\", n_rounds, 526)\n",
    "\n",
    "# Refit on all training years with the chosen number of trees\n",
    "final = None\n",
    "if n_rounds is not None:\n",
    "    final = xgb.train(params, xgb.DMatrix(Xtr, label=ytr, base_margin=np.log(etr)), num_boost_round=n_rounds)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-010",
   "metadata": {},
   "source": [
    "## Your turn 3 — a fair comparison\n",
    "\n",
    "Five-fold `GroupKFold(5, shuffle=True, random_state=0)` over the training years,\n",
    "grouped by `policyholder_id`. In each fold, fit and score (mean Poisson deviance on\n",
    "the left-out fold):\n",
    "- a GLM on raw features: `C(vehicle_style, Treatment('Sedan')) + driver_age + urban +\n",
    "  annual_mileage + channel + holder_gender`\n",
    "- the engineered GLM (`GLM` above)\n",
    "- XGBoost with `n_rounds` trees\n",
    "\n",
    "Store the mean deviance of each in a dict `cv`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-011",
   "metadata": {},
   "outputs": [],
   "source": [
    "NAIVE = \"C(vehicle_style, Treatment('Sedan')) + driver_age + urban + annual_mileage + channel + holder_gender\"\n",
    "cv = None\n",
    "if n_rounds is not None:\n",
    "    ...  # your code here\n",
    "check(\"CV deviance, raw GLM\", None if cv is None else round(cv[\"raw GLM\"], 5), 0.34327)\n",
    "check(\"CV deviance, GLM\", None if cv is None else round(cv[\"GLM\"], 5), 0.33994)\n",
    "check(\"CV deviance, XGBoost\", None if cv is None else round(cv[\"XGBoost\"], 5), 0.34082)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-012",
   "metadata": {},
   "source": [
    "## 3. SHAP values\n",
    "\n",
    "XGBoost computes exact TreeSHAP values with `pred_contribs=True`. With\n",
    "`base_margin = 0`, they explain the log of the **annual** claim rate: the last column\n",
    "is the base value, and each row sums to the log rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-013",
   "metadata": {},
   "outputs": [],
   "source": [
    "def shap_explain(model, X):\n",
    "    contrib = model.predict(xgb.DMatrix(X, base_margin=np.zeros(len(X))), pred_contribs=True)\n",
    "    return shap.Explanation(values=contrib[:, :-1], base_values=contrib[:, -1], data=X.to_numpy(),\n",
    "                            feature_names=list(X.columns))\n",
    "\n",
    "\n",
    "if final is not None:\n",
    "    samp = Xho.sample(3000, random_state=0)\n",
    "    sv = shap_explain(final, samp)\n",
    "    shap.plots.beeswarm(sv, max_display=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-014",
   "metadata": {},
   "source": [
    "## Your turn 4 — read the SHAP values\n",
    "\n",
    "1. The base rate: `exp` of the base value.\n",
    "2. The feature with the largest mean |SHAP value|.\n",
    "3. A dependence plot of `driver_age` (x) against its SHAP value (y), sports cars in a\n",
    "   different color. Do you see the U, and the young-sports-car effect?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-015",
   "metadata": {},
   "outputs": [],
   "source": [
    "base_rate = top_feature = None\n",
    "if final is not None:\n",
    "    ...  # your code here\n",
    "check(\"SHAP base rate\", None if base_rate is None else round(base_rate, 4), 0.0636)\n",
    "check(\"top feature\", top_feature, \"annual_mileage\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-016",
   "metadata": {},
   "source": [
    "## 4. One policy, explained\n",
    "\n",
    "A waterfall for a young, urban sports-car driver. `exp` of each contribution is a\n",
    "per-policy rating factor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-017",
   "metadata": {},
   "outputs": [],
   "source": [
    "if final is not None:\n",
    "    one = Xho[(Xho[\"driver_age\"] < 22) & (Xho[\"style_Sports\"] == 1) & (Xho[\"urban\"] == 1)].iloc[[0]]\n",
    "    sv1 = shap_explain(final, one)\n",
    "    shap.plots.waterfall(sv1[0], max_display=8)\n",
    "    print(\"annual frequency:\", float(np.exp(sv1.base_values[0] + sv1.values[0].sum())))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}