Chapter18SensitivityAnalysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 18: Sensitivity Analysis for the Average Causal Effect with Unmeasured Confounding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "graphviz: 0.20.1\n",
      "scipy   : 1.13.0\n",
      "pandas  : 2.2.2\n",
      "numpy   : 1.26.4\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from utils import *\n",
    "\n",
    "np.random.seed(42)\n",
    "%load_ext autoreload\n",
    "%autoreload 1\n",
    "\n",
    "%load_ext watermark\n",
    "%watermark --iversions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression, LinearRegression\n",
    "\n",
    "\n",
    "def OS_est(\n",
    "    z,\n",
    "    y,\n",
    "    x,\n",
    "    omod=LinearRegression(),\n",
    "    pmod=LogisticRegression(C=1e4),\n",
    "    lb=0,\n",
    "    ub=1,\n",
    "    e1=1,\n",
    "    e0=1,\n",
    "):\n",
    "    pscore = pmod.fit(x, z).predict_proba(x)[:, 1]\n",
    "    pscore = np.clip(pscore, lb, ub)\n",
    "    # fitted potential outcomes\n",
    "    outcome1 = omod.fit(x[z == 1, :], y[z == 1]).predict(x)\n",
    "    outcome0 = omod.fit(x[z == 0, :], y[z == 0]).predict(x)\n",
    "\n",
    "    ## outcome regression estimator\n",
    "    ace_reg = (\n",
    "        np.mean(z * y)\n",
    "        + np.mean((1 - z) * outcome1 / e1)\n",
    "        - np.mean(z * outcome0 * e0)\n",
    "        - np.mean((1 - z) * y)\n",
    "    )\n",
    "    ## IPW estimators\n",
    "    w1 = pscore + (1 - pscore) / e1\n",
    "    w0 = pscore * e0 + (1 - pscore)\n",
    "    ace_ipw0 = np.mean(z * y * w1 / pscore) - np.mean((1 - z) * y * w0 / (1 - pscore))\n",
    "    ace_ipw = np.mean(z * y * w1 / pscore) / np.mean(z / pscore) - np.mean(\n",
    "        (1 - z) * y * w0 / (1 - pscore)\n",
    "    ) / np.mean((1 - z) / (1 - pscore))\n",
    "    ## doubly robust estimator\n",
    "    aug = outcome1 / pscore / e1 + outcome0 * e0 / (1 - pscore)\n",
    "    ace_dr = ace_ipw0 - np.mean((z - pscore) * aug)\n",
    "\n",
    "    return np.array([ace_reg, ace_ipw0, ace_ipw, ace_dr])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhanes_bmi = pd.read_csv(\"nhanes_bmi.csv\").iloc[:, 1:]\n",
    "z, y, x = (\n",
    "    nhanes_bmi.School_meal.values,\n",
    "    nhanes_bmi.BMI.values,\n",
    "    nhanes_bmi.iloc[:, 2:].values,\n",
    ")\n",
    "x = (x - x.mean(0)) / x.std(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "E1 = np.array([1 / 2, 1 / 1.7, 1 / 1.5, 1 / 1.3, 1, 1.3, 1.5, 1.7, 2])\n",
    "E0 = E1.copy()\n",
    "est = np.zeros((len(E1), len(E0)))\n",
    "\n",
    "for i in range(len(E1)):\n",
    "    for j in range(len(E0)):\n",
    "        est[i, j] = OS_est(z, y, x, e1=E1[i], e0=E0[j])[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0.500000</th>\n",
       "      <th>0.588235</th>\n",
       "      <th>0.666667</th>\n",
       "      <th>0.769231</th>\n",
       "      <th>1.000000</th>\n",
       "      <th>1.300000</th>\n",
       "      <th>1.500000</th>\n",
       "      <th>1.700000</th>\n",
       "      <th>2.000000</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.500000</th>\n",
       "      <td>14.62</td>\n",
       "      <td>13.62</td>\n",
       "      <td>12.73</td>\n",
       "      <td>11.57</td>\n",
       "      <td>8.96</td>\n",
       "      <td>5.57</td>\n",
       "      <td>3.30</td>\n",
       "      <td>1.04</td>\n",
       "      <td>-2.36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.588235</th>\n",
       "      <td>11.93</td>\n",
       "      <td>10.93</td>\n",
       "      <td>10.04</td>\n",
       "      <td>8.88</td>\n",
       "      <td>6.27</td>\n",
       "      <td>2.87</td>\n",
       "      <td>0.61</td>\n",
       "      <td>-1.66</td>\n",
       "      <td>-5.05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.666667</th>\n",
       "      <td>10.13</td>\n",
       "      <td>9.13</td>\n",
       "      <td>8.24</td>\n",
       "      <td>7.08</td>\n",
       "      <td>4.47</td>\n",
       "      <td>1.08</td>\n",
       "      <td>-1.19</td>\n",
       "      <td>-3.45</td>\n",
       "      <td>-6.85</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.769231</th>\n",
       "      <td>8.33</td>\n",
       "      <td>7.33</td>\n",
       "      <td>6.45</td>\n",
       "      <td>5.29</td>\n",
       "      <td>2.67</td>\n",
       "      <td>-0.72</td>\n",
       "      <td>-2.98</td>\n",
       "      <td>-5.25</td>\n",
       "      <td>-8.64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.000000</th>\n",
       "      <td>5.64</td>\n",
       "      <td>4.64</td>\n",
       "      <td>3.75</td>\n",
       "      <td>2.59</td>\n",
       "      <td>-0.02</td>\n",
       "      <td>-3.41</td>\n",
       "      <td>-5.68</td>\n",
       "      <td>-7.94</td>\n",
       "      <td>-11.34</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.300000</th>\n",
       "      <td>3.57</td>\n",
       "      <td>2.57</td>\n",
       "      <td>1.68</td>\n",
       "      <td>0.52</td>\n",
       "      <td>-2.09</td>\n",
       "      <td>-5.49</td>\n",
       "      <td>-7.75</td>\n",
       "      <td>-10.01</td>\n",
       "      <td>-13.41</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.500000</th>\n",
       "      <td>2.65</td>\n",
       "      <td>1.65</td>\n",
       "      <td>0.76</td>\n",
       "      <td>-0.40</td>\n",
       "      <td>-3.01</td>\n",
       "      <td>-6.41</td>\n",
       "      <td>-8.67</td>\n",
       "      <td>-10.93</td>\n",
       "      <td>-14.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.700000</th>\n",
       "      <td>1.94</td>\n",
       "      <td>0.94</td>\n",
       "      <td>0.06</td>\n",
       "      <td>-1.11</td>\n",
       "      <td>-3.72</td>\n",
       "      <td>-7.11</td>\n",
       "      <td>-9.38</td>\n",
       "      <td>-11.64</td>\n",
       "      <td>-15.03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2.000000</th>\n",
       "      <td>1.15</td>\n",
       "      <td>0.15</td>\n",
       "      <td>-0.74</td>\n",
       "      <td>-1.90</td>\n",
       "      <td>-4.51</td>\n",
       "      <td>-7.90</td>\n",
       "      <td>-10.17</td>\n",
       "      <td>-12.43</td>\n",
       "      <td>-15.83</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          0.500000  0.588235  0.666667  0.769231  1.000000  1.300000  \\\n",
       "0.500000     14.62     13.62     12.73     11.57      8.96      5.57   \n",
       "0.588235     11.93     10.93     10.04      8.88      6.27      2.87   \n",
       "0.666667     10.13      9.13      8.24      7.08      4.47      1.08   \n",
       "0.769231      8.33      7.33      6.45      5.29      2.67     -0.72   \n",
       "1.000000      5.64      4.64      3.75      2.59     -0.02     -3.41   \n",
       "1.300000      3.57      2.57      1.68      0.52     -2.09     -5.49   \n",
       "1.500000      2.65      1.65      0.76     -0.40     -3.01     -6.41   \n",
       "1.700000      1.94      0.94      0.06     -1.11     -3.72     -7.11   \n",
       "2.000000      1.15      0.15     -0.74     -1.90     -4.51     -7.90   \n",
       "\n",
       "          1.500000  1.700000  2.000000  \n",
       "0.500000      3.30      1.04     -2.36  \n",
       "0.588235      0.61     -1.66     -5.05  \n",
       "0.666667     -1.19     -3.45     -6.85  \n",
       "0.769231     -2.98     -5.25     -8.64  \n",
       "1.000000     -5.68     -7.94    -11.34  \n",
       "1.300000     -7.75    -10.01    -13.41  \n",
       "1.500000     -8.67    -10.93    -14.33  \n",
       "1.700000     -9.38    -11.64    -15.03  \n",
       "2.000000    -10.17    -12.43    -15.83  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(est, columns=E0, index=E1).round(2)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "metrics",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}