{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-000",
   "metadata": {},
   "source": [
    "# Lesson 1 — Statistical Foundations of Predictive Modeling\n",
    "\n",
    "**Companion notebook** for the Lesson 1 slides. You'll fit OLS to health-claim severity,\n",
    "see why it struggles, move to a log-linear model, correct the retransformation bias,\n",
    "and use heteroscedasticity-robust standard errors.\n",
    "\n",
    "Data: `../data/health/encounters.parquet` (one row per encounter) and\n",
    "`../data/health/members.csv` (one row per member)."
   ]
  },
  {
   "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",
    "from statsmodels.stats.diagnostic import het_breuschpagan\n",
    "\n",
    "enc = pd.read_parquet(\"../data/health/encounters.parquet\")\n",
    "mem = pd.read_csv(\"../data/health/members.csv\")\n",
    "e = (enc[enc[\"service_date\"].dt.year == 2024]\n",
    "     .merge(mem, on=\"member_id\", validate=\"many_to_one\"))\n",
    "e[\"age10\"] = (e[\"age\"] - 40) / 10          # decades from age 40\n",
    "len(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-002",
   "metadata": {},
   "source": [
    "## 1. The shape of insurance data\n",
    "\n",
    "Compare the mean with the median, and look at how the spread grows with the mean\n",
    "across service categories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-003",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = e[\"allowed\"]\n",
    "print(f\"mean {a.mean():,.2f}  median {a.median():,.2f}  skewness {st.skew(a):.1f}\")\n",
    "print(f\"share of encounters above the mean: {(a > a.mean()).mean():.3f}\")\n",
    "\n",
    "by_cat = e.groupby(\"service_category\")[\"allowed\"].agg([\"count\", \"mean\", \"std\"])\n",
    "by_cat[\"cv\"] = by_cat[\"std\"] / by_cat[\"mean\"]\n",
    "by_cat.sort_values(\"mean\").round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-004",
   "metadata": {},
   "source": [
    "## 2. OLS on raw dollars\n",
    "\n",
    "`C(col, Treatment('level'))` creates dummy variables and names the base level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-005",
   "metadata": {},
   "outputs": [],
   "source": [
    "rhs = (\"C(service_category, Treatment('Office Visit')) + age10\"\n",
    "       \" + C(region, Treatment('Suburban')) + C(plan, Treatment('PPO'))\")\n",
    "ols = smf.ols(\"allowed ~ \" + rhs, data=e).fit()\n",
    "print(f\"R-squared {ols.rsquared:.3f}\")\n",
    "print(f\"negative fitted values: {(ols.fittedvalues < 0).sum():,}\")\n",
    "ols.params.round(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-006",
   "metadata": {},
   "source": [
    "Residuals vs fitted and a Q–Q plot tell you more than R². Look for a fan shape and\n",
    "for heavy tails."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-007",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "samp = e.sample(6_000, random_state=1).index\n",
    "fig, axes = plt.subplots(1, 2, figsize=(11, 3.8))\n",
    "axes[0].scatter(ols.fittedvalues[samp], ols.resid[samp], s=5, alpha=0.3)\n",
    "axes[0].axhline(0, color=\"orange\")\n",
    "axes[0].set(title=\"Residuals vs fitted\", xlabel=\"Fitted allowed\", ylabel=\"Residual\")\n",
    "sm.qqplot(ols.resid[samp] / ols.resid.std(), line=\"45\", ax=axes[1], markersize=3, alpha=0.4)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-008",
   "metadata": {},
   "source": [
    "## 3. The log-linear model\n",
    "\n",
    "Model `log(allowed)`: each coefficient becomes a multiplicative relativity `exp(beta)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-009",
   "metadata": {},
   "outputs": [],
   "source": [
    "e[\"log_allowed\"] = np.log(e[\"allowed\"])\n",
    "loglin = smf.ols(\"log_allowed ~ \" + rhs, data=e).fit()\n",
    "print(f\"R-squared {loglin.rsquared:.3f}, residual sd {np.sqrt(loglin.scale):.3f}\")\n",
    "np.exp(loglin.params).round(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-010",
   "metadata": {},
   "source": [
    "---\n",
    "## Your turn 1 — relativities with robust intervals\n",
    "\n",
    "Refit the log-linear model with `cov_type=\"HC3\"` and build a table of relativities\n",
    "(`exp` of the coefficients) with 95% confidence limits. Store the Metro relativity and\n",
    "the HMO relativity. Compare them with the truth used to simulate the data:\n",
    "Metro 1.15, HMO 0.88."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-011",
   "metadata": {},
   "outputs": [],
   "source": [
    "hc3 = metro = hmo = None\n",
    "...  # your code here\n",
    "check(\"Metro relativity\", None if metro is None else round(metro, 3), 1.163)\n",
    "check(\"HMO relativity\", None if hmo is None else round(hmo, 3), 0.874)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-012",
   "metadata": {},
   "source": [
    "## Your turn 2 — is the variance constant?\n",
    "\n",
    "1. Run the Breusch–Pagan test on the log-linear model’s residuals.\n",
    "2. Compute the residual standard deviation for each service category.\n",
    "\n",
    "Store the LM statistic (rounded to the nearest whole number) and the category with the\n",
    "largest residual standard deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-013",
   "metadata": {},
   "outputs": [],
   "source": [
    "bp_lm = widest = None\n",
    "...  # your code here\n",
    "check(\"Breusch-Pagan LM\", bp_lm, 21148)\n",
    "check(\"widest category\", widest, \"Pharmacy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-014",
   "metadata": {},
   "source": [
    "## Your turn 3 — how much do robust standard errors change?\n",
    "\n",
    "For the Specialist coefficient, compute the ratio HC3 SE ÷ ordinary SE. (Hint: `.bse`.)\n",
    "Robust standard errors are *correct*, not automatically bigger."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-015",
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio = None\n",
    "if hc3 is not None:\n",
    "    ...  # your code here\n",
    "check(\"HC3 / OLS standard error, Specialist\", None if ratio is None else round(ratio, 2), 0.43)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-016",
   "metadata": {},
   "source": [
    "## 4. Back to dollars: retransformation\n",
    "\n",
    "`exp(fitted)` estimates the **median**. Duan’s smearing factor, the mean of\n",
    "`exp(residual)`, turns it into a mean — but one global factor is wrong when the residual\n",
    "variance differs by category."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-017",
   "metadata": {},
   "outputs": [],
   "source": [
    "naive = np.exp(loglin.fittedvalues)\n",
    "smear = np.mean(np.exp(loglin.resid))\n",
    "print(f\"actual mean {e['allowed'].mean():,.1f}; naive {naive.mean():,.1f}; \"\n",
    "      f\"one smearing factor ({smear:.3f}) {(naive * smear).mean():,.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-018",
   "metadata": {},
   "source": [
    "## Your turn 4 — the lab: baseline severity for a reference member\n",
    "\n",
    "Compute a smearing factor **per service category**, then predict the median and the\n",
    "mean allowed amount for a reference member (age 40, Suburban, PPO) in every category.\n",
    "Store the predicted *mean* for Pharmacy and for Inpatient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-019",
   "metadata": {},
   "outputs": [],
   "source": [
    "ref = pd.DataFrame({\"service_category\": sorted(e[\"service_category\"].unique()),\n",
    "                    \"age10\": 0.0, \"region\": \"Suburban\", \"plan\": \"PPO\"})\n",
    "pharmacy_mean = inpatient_mean = None\n",
    "...  # your code here\n",
    "check(\"Pharmacy mean (reference member)\", None if pharmacy_mean is None else round(pharmacy_mean, 2), 161.83)\n",
    "check(\"Inpatient mean (reference member)\", None if inpatient_mean is None else round(inpatient_mean, 2), 36205.44)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-020",
   "metadata": {},
   "source": [
    "## 5. Preview: a log *link* instead of a log *transform*\n",
    "\n",
    "A Gamma GLM with a log link models the mean in dollars directly — no smearing needed.\n",
    "Lessons 2 and 3 build on this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-021",
   "metadata": {},
   "outputs": [],
   "source": [
    "glm = smf.glm(\"allowed ~ \" + rhs, data=e,\n",
    "              family=sm.families.Gamma(sm.families.links.Log())).fit()\n",
    "np.exp(glm.params).round(3)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}