{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "colab": {
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "id": "m14-l2-0",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesson 14.2 · Rosenbaum Bounds: How Much Hidden Bias Would Overturn a Matched Study?\n",
    "\n",
    "*Accelerated Foundation in Applied Causal Inference*\n",
    "\n",
    "**By the end of this lesson you can:**\n",
    "\n",
    "- Explain the sensitivity parameter Γ as a bound on hidden bias in treatment odds\n",
    "- Build matched pairs from the LaLonde data\n",
    "- Implement the Wilcoxon signed-rank Rosenbaum bound and find the Γ at which significance disappears\n",
    "- Interpret and report a sensitivity analysis in plain language\n",
    "\n",
    "> Read the full explanations on the course website, then run each cell here in order (**Runtime → Run all** works on a fresh runtime)."
   ]
  },
  {
   "id": "m14-l2-1",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The question sensitivity analysis answers"
   ]
  },
  {
   "id": "m14-l2-2",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "!pip install -q numpy pandas scipy statsmodels scikit-learn matplotlib causaldata\n",
    "import numpy as np, pandas as pd\n",
    "from scipy.stats import norm, rankdata\n",
    "import statsmodels.formula.api as smf\n",
    "from sklearn.neighbors import NearestNeighbors\n",
    "import matplotlib.pyplot as plt\n",
    "from causaldata import nsw_mixtape, cps_mixtape"
   ]
  },
  {
   "id": "m14-l2-3",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matched pairs from the LaLonde data"
   ]
  },
  {
   "id": "m14-l2-4",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "nsw = nsw_mixtape.load_pandas().data\n",
    "obs = pd.concat([nsw[nsw.treat == 1], cps_mixtape.load_pandas().data], ignore_index=True)\n",
    "obs[\"u74\"] = (obs.re74 == 0).astype(int); obs[\"u75\"] = (obs.re75 == 0).astype(int)\n",
    "ps = smf.logit(\"treat ~ age + I(age**2) + I(age**3) + educ + I(educ**2) + marr + nodegree + black + hisp\"\n",
    "               \" + re74 + re75 + u74 + u75 + educ:re74\", data=obs).fit(disp=0).predict(obs)\n",
    "obs[\"lps\"] = np.log(ps / (1 - ps))\n",
    "\n",
    "treated = obs[obs.treat == 1].sort_values(\"lps\", ascending=False)   # hardest-to-match first\n",
    "controls = obs[obs.treat == 0].copy()\n",
    "pairs = []\n",
    "for _, t in treated.iterrows():\n",
    "    nn = NearestNeighbors(n_neighbors=1).fit(controls[[\"lps\"]])\n",
    "    j = controls.index[nn.kneighbors([[t.lps]])[1][0][0]]\n",
    "    pairs.append((t.re78, controls.loc[j, \"re78\"]))\n",
    "    controls = controls.drop(j)\n",
    "pairs = pd.DataFrame(pairs, columns=[\"y_treated\", \"y_control\"])\n",
    "pairs[\"diff\"] = pairs.y_treated - pairs.y_control\n",
    "print(f\"{len(pairs)} pairs | mean difference ${pairs['diff'].mean():,.0f} | median ${pairs['diff'].median():,.0f}\")"
   ]
  },
  {
   "id": "m14-l2-5",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rosenbaum bounds step by step"
   ]
  },
  {
   "id": "m14-l2-6",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question it answers:** How strong would an unmeasured factor need to be to explain away the earnings gain from NSW training?"
   ]
  },
  {
   "id": "m14-l2-7",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "def rosenbaum_pvalue(diffs, gamma):\n",
    "    d = diffs[diffs != 0]\n",
    "    ranks = rankdata(np.abs(d))\n",
    "    t_plus = ranks[d > 0].sum()\n",
    "    n = len(d)\n",
    "    p_plus = gamma / (1 + gamma)\n",
    "    mean = p_plus * n * (n + 1) / 2\n",
    "    var = p_plus * (1 - p_plus) * n * (n + 1) * (2 * n + 1) / 6\n",
    "    return 1 - norm.cdf((t_plus - mean) / np.sqrt(var))          # one-sided, worst case\n",
    "\n",
    "gammas = np.round(np.arange(1.0, 2.01, 0.05), 2)\n",
    "pvals = [rosenbaum_pvalue(pairs[\"diff\"].to_numpy(), g) for g in gammas]\n",
    "table = pd.DataFrame({\"Gamma\": gammas, \"worst-case p\": pvals}).round(4)\n",
    "print(table.iloc[::2].to_string(index=False))\n",
    "crit = next((g for g, p in zip(gammas, pvals) if p > 0.05), None)\n",
    "print(f\"Significance at 5% is lost at Gamma ≈ {crit}\")\n",
    "\n",
    "plt.plot(gammas, pvals, marker=\"o\"); plt.axhline(0.05, color=\"red\", ls=\"--\")\n",
    "plt.xlabel(\"Γ (hidden bias)\"); plt.ylabel(\"Worst-case one-sided p-value\"); plt.title(\"Rosenbaum sensitivity analysis\"); plt.show()"
   ]
  },
  {
   "id": "m14-l2-8",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading Γ and reporting it"
   ]
  },
  {
   "id": "m14-l2-9",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✍️ Exercise\n",
    "\n",
    "Simulate a matched study with a large true effect (pair differences ~ Normal(3, 2), n = 300) and a small one\n",
    "(Normal(0.4, 2)). Find the critical Γ for each. What determines robustness to hidden bias: statistical\n",
    "significance, effect size, or both?"
   ]
  },
  {
   "id": "m14-l2-10",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Your answer here"
   ]
  },
  {
   "id": "m14-l2-11",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checklist and recap"
   ]
  }
 ]
}