{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "W3Vr9I8k5xYZ",
"toc-hr-collapsed": false
},
"source": [
"# Ch12 逆概率加权和边际结构模型\n",
"\n",
"IP Weighting and Marginal Structural Models "
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"toc-hr-collapsed": true
},
"source": [
"## 内容摘要"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "S2rrX3Lm5xYf"
},
"source": [
"本书第二部分围绕因果问题“戒烟对体重增加的平均因果效应是什么？”进行组织。在本章中，我们描述如何使用IP加权从观测数据中估计这种效应。尽管在第二章中介绍了IP加权，但我们仅将其描述为非参数方法。We now describe the use of models together with IP weighting，在额外的假设下，这些模型将使我们能够处理高维问题 with many covariates and nondichotomous treatments.\n",
"\n",
"为了估计戒烟对体重增加的影响，我们将使用NHEFS的真实数据，该缩写词代表国家卫生和营养检查调查数据I流行病学随访研究。NHEFS是由国家卫生统计中心和国家老龄研究所与美国公共卫生局其他机构共同发起的。\n",
"可以在 www.cdc.gov/nchs/nhanes/nhefs/ 上找到NHEFS的详细说明以及可公开获得的数据集和文档。\n",
"在本章及以后的章节中，我们将使用NHEFS数据的一部分，该数据可从本书的网站上获得。\n",
"\n",
"> 我们关注的是戒烟是否会引起体重增加？"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "kestfkDK5xY3"
},
"source": [
"### 12.1 The causal question\n",
"\n",
"> 戒烟是否会引起体重增加？\n",
"$\\theta = E[Y^{a=1}] - E[Y^{a=0}]$如何估计呢？\n",
"\n",
"Our goal is to estimate the average causal effect of smoking cessation (the treatment) A on weight gain (the outcome) Y ."
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "3NG_uc6t5xZ2"
},
"source": [
"### 12.2 Estimating IP weights via modeling\n",
"\n",
"必须回答的一个关键问题是： if $A$ is a continuous variable, then IPW 方法的推广是什么？更简单的如果 $A$ 是个离散变量，如何使用一组权重 $W$ 使得加权总体中，$L$ 的分布并不影响 $A$ 的分布。"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "S8ZdOh6h5xZ3"
},
"source": [
"IP weighting creates a pseudo-population in which the arrow from the covariates $L$ to the treatment $A$ is removed. 也就是说\n",
"\n",
"$L \\rightarrow A$ for population $S$, then we create a pseudo-population $S'$ in which the distribution of $A$ is not influenced by any information about $L$ keeping everything else unchange. 这里就导致一个悖论，在视角 $S$ 下，$L$ 是 $A$ 的原因，但是换一个视角 $S'$ 他就不再是其原因了。\n",
"\n",
"- 一种解释方式是 effect 是 cause 的一种测量。例如测量体重环境变化了，太空中重量就测出来是以个常数，或者说测不出来。 这种解释的问题在于新的环境意味着新的因果机制，pseudo-population 的 setting 变化只是视角的变化，而不是环境或者因果机制的变化。\n",
"- 另外一种解释方式是，只存在 individual-level causal effect 的因果关系不变性，在群体层次上的因果关系并存在 invariant 因果关系，群体层次上的因果关系是一种统计特性。具体来说具备不变性的是个体因果关系，而因果效应确实一个离开了特定总体就没有定义的量 which is a statistical quantity。这种解释的问题在于\n",
"\n",
"这里有一个问题 $A$ 是 $L$ 的原因, 那么 $A$ 一定对 $L$ 有因果效应。"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "mxjyJePD5xbe"
},
"source": [
"### 12.3 Stabilized IP weights\n",
"\n",
"Create the stabilized IP weights $SW^A = f(A) \\, / \\, f(A|L)$"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "rix4qsWk5xbi"
},
"source": [
"\"The effect estimate obtained in the pseudo-population created by weights $0.5 \\, / \\, f(A|L)$\n",
"is equal to that obtained in the pseudo-population created by weights $1 \\, / \\, f(A|L)$.\""
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "1dPCV3JQ5xc4",
"toc-hr-collapsed": false
},
"source": [
"### 12.4 Marginal structural models\n",
"\n",
"Models for the marginal mean of a counterfactual outcome are referred to as marginal structural mean models. 反事实结果的边际均值模型称为边际结构均值模型。"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"模型的各种变体，从离散的连续。\n",
"\n",
"Marginal structural mean models within a pseudo-population:\n",
"\n",
"$$E[Y^a] = \\beta_0 + \\beta_1 a$$ \n",
"$$\\Rightarrow E[Y^a] = \\beta_0 + \\beta_1 a + \\beta_2 a^2$$ \n",
"\n",
"\n",
"marginal structural logistic model\n",
"\n",
"$$logit Pr[Y^a = 1] = \\beta_0 + \\beta_1 a$$ "
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "qN_aHhL45xee"
},
"source": [
"### 12.5 Effect modification and marginal structural models\n",
"\n",
"Marginal structural models do not include covariates when the target parameter is the average causal effect in the population. However, one may include covariates–which may be non-confounders–in a marginal structural model to assess effect modification."
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "XZoJSbE95xef"
},
"source": [
"在第一部分中，我们讨论了 effect modification and confounding 是两个在逻辑上截然不同的概念。 但是，许多学生很难理解这种区别，因为经常使用相同的统计方法（分层（第4章）或回归（第15章）for confounder adjustment and detection of effect modification. 因此，使用边际结构模型来讲授这些概念可能会有一些优势， because then methods for confounder adjustment (IP weighting) are distinct from methods for detection of effect modification (adding treatment-covariate product terms to a marginal structural model).\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "gu14xRuZ5xev"
},
"source": [
"### 12.6 Censoring and missing data\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "true",
"colab_type": "text",
"id": "tAbhsGcA5xYg",
"toc-hr-collapsed": true
},
"source": [
"## Programs\n",
"\n",
"**Setup and imports**"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "1TMm0xOW5xYh"
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"import warnings\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"from collections import OrderedDict\n",
"import numpy as np\n",
"import pandas as pd\n",
"import statsmodels.api as sm\n",
"import scipy.stats\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "v7UvRwZ15xYq",
"outputId": "fe462c77-2858-472d-f4bb-e0e9171110f9"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING *** OLE2 inconsistency: SSCS size is 0 but SSAT size is non-zero\n"
]
},
{
"data": {
"text/plain": [
"(1629, 64)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Just a look at a couple basic details of the dataset\n",
"nhefs_all = pd.read_excel('NHEFS.xls')\n",
"nhefs_all.shape"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"data": {
"text/plain": [
"( Variable name Description\n",
" 0 active in your usual day, how active are you? in 1971...\n",
" 1 age age in 1971\n",
" 2 alcoholfreq how often do you drink? in 1971 0: almost ev...\n",
" 3 alcoholhowmuch when you drink, how much do you drink? in 1971\n",
" 4 alcoholpy have you had 1 drink past year? in 1971, 1:ev...,\n",
" ['in your usual day, how active are you? in 1971, 0:very active, 1:moderately active, 2:inactive',\n",
" 'age in 1971',\n",
" 'how often do you drink? in 1971 0: almost every day, 1: 2-3 times/week, 2: 1-4 times/month, 3: < 12 times/year, 4: no alcohol last year, 5: unknown',\n",
" 'when you drink, how much do you drink? in 1971',\n",
" 'have you had 1 drink past year? in 1971, 1:ever, 0:never; 2:missing'])"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"codebook = pd.read_excel('codebook.xls')\n",
"codebook.Description = codebook.Description.apply(lambda s: s.lower())\n",
"var_list = codebook.iloc[:, 1].tolist()\n",
"codebook.head(), var_list[:5]"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" qsmk \n",
" wt82_71 \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 0 \n",
" -10.093960 \n",
" \n",
" \n",
" 1 \n",
" 0 \n",
" 2.604970 \n",
" \n",
" \n",
" 2 \n",
" 0 \n",
" 9.414486 \n",
" \n",
" \n",
" 3 \n",
" 0 \n",
" 4.990117 \n",
" \n",
" \n",
" 4 \n",
" 0 \n",
" 4.989251 \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
" coef std err t P>|t| [0.025 0.975] \n",
" \n",
"\n",
" constant 1.9845 0.229 8.672 0.000 1.536 2.433 \n",
" \n",
"\n",
" qsmk 2.5406 0.451 5.632 0.000 1.656 3.425 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 7,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"ols = sm.OLS(nhefs.wt82_71, nhefs[['constant', 'qsmk']])\n",
"res = ols.fit()\n",
"res.summary().tables[1]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "qfvnavEW5xZo",
"outputId": "ed4a51c0-1521-4b33-d68e-e35fd4930de0"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" estimate 95% C.I.\n",
"difference 2.54 (1.7, 3.4)\n"
]
}
],
"source": [
"est = res.params.qsmk\n",
"conf_ints = res.conf_int(alpha=0.05, cols=None)\n",
"lo, hi = conf_ints[0]['qsmk'], conf_ints[1]['qsmk']\n",
"\n",
"print(' estimate 95% C.I.')\n",
"print('difference {:>6.2f} ({:>0.1f}, {:>0.1f})'.format(est, lo, hi))"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "GG8eIqMa5xZu"
},
"source": [
"Create Table 12.1 in the margin of pg 149."
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "7rIpPfIv5xZv"
},
"source": [
"As shown in Table 12.1, quitters and non-quitters also differed in their distribution of other variables such as sex, race, education, baseline weight, and intensity of smoking. If these variables are confounders, then they also need to be adjusted for in the analysis. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "XIEF4gyK5xZw",
"outputId": "c2fb3d20-9537-448b-b652-3a49805b40d0"
},
"outputs": [
{
"data": {
"text/html": [
" qsmk 1 0 \n",
" \n",
" Age, years \n",
" 46.2 \n",
" 42.8 \n",
" \n",
" \n",
" Men, % \n",
" 54.6 \n",
" 46.6 \n",
" \n",
" \n",
" White, % \n",
" 91.1 \n",
" 85.4 \n",
" \n",
" \n",
" University education, % \n",
" 15.4 \n",
" 9.9 \n",
" \n",
" \n",
" Weight, kg \n",
" 72.4 \n",
" 70.3 \n",
" \n",
" \n",
" Cigarettes/day \n",
" 18.6 \n",
" 21.2 \n",
" \n",
" \n",
" Years smoking \n",
" 26.0 \n",
" 24.1 \n",
" \n",
" \n",
" Little or no exercise, % \n",
" 40.7 \n",
" 37.9 \n",
" \n",
" \n",
" Inactive daily life, % \n",
" 11.2 \n",
" 8.9 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 9,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"summaries = OrderedDict((\n",
" ('age', 'mean'),\n",
" ('sex', lambda x: (100 * (x == 0)).mean()),\n",
" ('race', lambda x: (100 * (x == 0)).mean()),\n",
" ('university', lambda x: 100 * x.mean()),\n",
" ('wt71', 'mean'),\n",
" ('smokeintensity', 'mean'),\n",
" ('smokeyrs', 'mean'),\n",
" ('no_exercise', lambda x: 100 * x.mean()),\n",
" ('inactive', lambda x: 100 * x.mean())\n",
"))\n",
"\n",
"table = nhefs.groupby('qsmk').agg(summaries)\n",
"table.sort_index(ascending=False, inplace=True)\n",
"table = table.T\n",
"\n",
"table.index = [\n",
" 'Age, years',\n",
" 'Men, %',\n",
" 'White, %',\n",
" 'University education, %',\n",
" 'Weight, kg',\n",
" 'Cigarettes/day',\n",
" 'Years smoking',\n",
" 'Little or no exercise, %',\n",
" 'Inactive daily life, %'\n",
"]\n",
"\n",
"table.style.format(\"{:>0.1f}\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "gkJjAU5I5xZ5"
},
"source": [
"### Program 12.2"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "l4XJRrJQ5xZ9"
},
"source": [
"We're going to be modeling with squared terms and some categorical features. Here we'll explicitly add squared features and dummy features to the data. In later chapters we'll use Statsmodels' formula syntax."
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "qCIwhWGW5xZ_"
},
"source": [
"Squared features:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "FwAU3FhT5xaA"
},
"outputs": [],
"source": [
"for col in ['age', 'wt71', 'smokeintensity', 'smokeyrs']:\n",
" nhefs['{}^2'.format(col)] = nhefs[col] * nhefs[col]"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "_JF_pRQX5xaG"
},
"source": [
"Dummy features:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "VNnr5q5A5xaH"
},
"outputs": [],
"source": [
"edu_dummies = pd.get_dummies(nhefs.education, prefix='edu')\n",
"exercise_dummies = pd.get_dummies(nhefs.exercise, prefix='exercise')\n",
"active_dummies = pd.get_dummies(nhefs.active, prefix='active')\n",
"\n",
"nhefs = pd.concat(\n",
" [nhefs, edu_dummies, exercise_dummies, active_dummies],\n",
" axis=1\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "0X6_4oq75xaN"
},
"source": [
"We're going to be creating a lot of IP weights from logistic regressions so a function will help reduce the work. The following function creates the denominators of the IP weights."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "mX-jTBLe5xaP"
},
"outputs": [],
"source": [
"def logit_ip_f(y, X):\n",
" \"\"\"\n",
" Create the f(y|X) part of IP weights\n",
" from logistic regression\n",
" \n",
" Parameters\n",
" ----------\n",
" y : Pandas Series\n",
" X : Pandas DataFrame\n",
" \n",
" Returns\n",
" -------\n",
" Numpy array of IP weights\n",
" \n",
" \"\"\"\n",
" model = sm.Logit(y, X)\n",
" res = model.fit()\n",
" weights = np.zeros(X.shape[0])\n",
" weights[y == 1] = res.predict(X.loc[y == 1])\n",
" weights[y == 0] = (1 - res.predict(X.loc[y == 0]))\n",
" return weights"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "ossLEwFm5xaS",
"outputId": "e04f9791-7225-4c20-e9cd-a5f9530fbbb8"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.535408\n",
" Iterations 6\n"
]
}
],
"source": [
"X_ip = nhefs[[\n",
" 'constant',\n",
" 'sex', 'race', 'age', 'age^2', 'edu_2', 'edu_3', 'edu_4', 'edu_5',\n",
" 'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2', \n",
" 'exercise_1', 'exercise_2', 'active_1', 'active_2', 'wt71', 'wt71^2'\n",
"]]\n",
"\n",
"denoms = logit_ip_f(nhefs.qsmk, X_ip)\n",
"weights = 1 / denoms"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "79mVEfCr5xaX",
"outputId": "94d58591-b984-463d-c339-2f1693e0ed0f"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"IP weights\n",
" min: 1.05 expected: 1.05\n",
" max: 16.70 expected: 16.70\n",
" mean: 2.00 expected: 2.00\n"
]
}
],
"source": [
"print('IP weights')\n",
"print(' min: {:>5.2f} expected: 1.05'.format(weights.min()))\n",
"print(' max: {:>5.2f} expected: 16.70'.format(weights.max()))\n",
"print(' mean: {:>5.2f} expected: 2.00'.format(weights.mean()))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "5XRZb2DY5xae",
"outputId": "ea3c0912-592e-4a76-f4ea-3af448e12c34"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFlCAYAAAApldtwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAASuUlEQVR4nO3dfaykZ3nf8d9VNiSBNLHBCyW7Vhcah4SipFhb4gQ1anGSYkDYfwSJKA0raslSRQkJaYNppCK1UmXaKLyoFZWFHRsFQZDj1lZCXixDElUqbtZAeHNSrwi1Fzt4I4OTBqWJlat/nMftYfesd/fM8Zm9zn4+0mrmueeemfvx2fX3PM/MmVPdHQDg/PY31r0AAODMBBsABhBsABhAsAFgAMEGgAEEGwAG2LfuBTyZSy65pA8dOrTuZQDArrn33nv/pLv3nzx+Xgf70KFDOXr06LqXAQC7pqr+11bjTokDwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwADn9W/reiocuv7XdvTxvnjDq3b08QBgK46wAWAAwQaAAQQbAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGECwAWAAwQaAAQQbAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGECwAWAAwQaAAQQbAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAHOGOyqurmqHqmqz24ae1ZV3VVV9y+XFy/jVVXvqapjVfXpqrp8032OLPPvr6ojT83uAMDedDZH2LckecVJY9cnubu7L0ty97KdJFcluWz5c12S9yYbgU/y9iTfl+SlSd7+ROQBgDM7Y7C7+3eTPHrS8NVJbl2u35rkmk3j7+8NH09yUVU9L8k/TnJXdz/a3V9JcldO/SYAADiN7b6G/dzufjhJlsvnLOMHkjy4ad7xZex046eoquuq6mhVHT1x4sQ2lwcAe8tOv+msthjrJxk/dbD7xu4+3N2H9+/fv6OLA4CpthvsLy+nurNcPrKMH09y6aZ5B5M89CTjAMBZ2G6w70zyxDu9jyS5Y9P465d3i1+R5LHllPlvJvmRqrp4ebPZjyxjAMBZ2HemCVX1wST/MMklVXU8G+/2viHJh6vq2iQPJHntMv0jSV6Z5FiSryV5Q5J096NV9W+T/N4y799098lvZAMATuOMwe7uHzvNTVduMbeTvPE0j3NzkpvPaXUAQBKfdAYAIwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAOsFOyq+umq+lxVfbaqPlhV31RVz6+qe6rq/qr65ap6+jL3G5ftY8vth3ZiBwDgQrDtYFfVgSQ/meRwd784ydOSvC7JO5K8s7svS/KVJNcud7k2yVe6+zuSvHOZBwCchVVPie9L8s1VtS/JM5I8nOTlSW5bbr81yTXL9auX7Sy3X1lVteLzA8AFYdvB7u4vJfn5JA9kI9SPJbk3yVe7+/Fl2vEkB5brB5I8uNz38WX+s7f7/ABwIVnllPjF2Thqfn6Sb0/yzCRXbTG1n7jLk9y2+XGvq6qjVXX0xIkT210eAOwpq5wS/6Ekf9TdJ7r7r5LcnuQHkly0nCJPkoNJHlquH09yaZIst39bkkdPftDuvrG7D3f34f3796+wPADYO1YJ9gNJrqiqZyyvRV+Z5PNJPpbkR5c5R5LcsVy/c9nOcvtHu/uUI2wA4FSrvIZ9TzbePPaJJJ9ZHuvGJG9N8paqOpaN16hvWu5yU5JnL+NvSXL9CusGgAvKvjNPOb3ufnuSt580/IUkL91i7l8kee0qzwcAFyqfdAYAAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMsFKwq+qiqrqtqv6gqu6rqu+vqmdV1V1Vdf9yefEyt6rqPVV1rKo+XVWX78wuAMDet+oR9ruT/EZ3f1eS701yX5Lrk9zd3ZcluXvZTpKrkly2/LkuyXtXfG4AuGBsO9hV9a1JfjDJTUnS3X/Z3V9NcnWSW5dptya5Zrl+dZL394aPJ7moqp637ZUDwAVklSPsFyQ5keQXq+qTVfW+qnpmkud298NJslw+Z5l/IMmDm+5/fBn7OlV1XVUdraqjJ06cWGF5ALB3rBLsfUkuT/Le7n5Jkj/P/z/9vZXaYqxPGei+sbsPd/fh/fv3r7A8ANg7Vgn28STHu/ueZfu2bAT8y0+c6l4uH9k0/9JN9z+Y5KEVnh8ALhjbDnZ3/3GSB6vqhcvQlUk+n+TOJEeWsSNJ7liu35nk9cu7xa9I8tgTp84BgCe3b8X7vynJB6rq6Um+kOQN2fgm4MNVdW2SB5K8dpn7kSSvTHIsydeWuQDAWVgp2N39qSSHt7jpyi3mdpI3rvJ8AHCh8klnADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwwMrBrqqnVdUnq+pXl+3nV9U9VXV/Vf1yVT19Gf/GZfvYcvuhVZ8bAC4UO3GE/eYk923afkeSd3b3ZUm+kuTaZfzaJF/p7u9I8s5lHgBwFlYKdlUdTPKqJO9btivJy5Pctky5Nck1y/Wrl+0st1+5zAcAzmDVI+x3JfnZJH+9bD87yVe7+/Fl+3iSA8v1A0keTJLl9seW+QDAGWw72FX16iSPdPe9m4e3mNpncdvmx72uqo5W1dETJ05sd3kAsKescoT9siSvqaovJvlQNk6FvyvJRVW1b5lzMMlDy/XjSS5NkuX2b0vy6MkP2t03dvfh7j68f//+FZYHAHvHtoPd3W/r7oPdfSjJ65J8tLt/PMnHkvzoMu1IkjuW63cu21lu/2h3n3KEDQCc6qn4Oey3JnlLVR3LxmvUNy3jNyV59jL+liTXPwXPDQB70r4zTzmz7v7tJL+9XP9CkpduMecvkrx2J54PAC40PukMAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGECwAWAAwQaAAQQbAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGECwAWAAwQaAAQQbAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGECwAWAAwQaAAQQbAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGGDbwa6qS6vqY1V1X1V9rqrevIw/q6ruqqr7l8uLl/GqqvdU1bGq+nRVXb5TOwEAe90qR9iPJ/mZ7v7uJFckeWNVvSjJ9Unu7u7Lkty9bCfJVUkuW/5cl+S9Kzw3AFxQth3s7n64uz+xXP+zJPclOZDk6iS3LtNuTXLNcv3qJO/vDR9PclFVPW/bKweAC8iOvIZdVYeSvCTJPUme290PJxtRT/KcZdqBJA9uutvxZezkx7quqo5W1dETJ07sxPIAYLyVg11V35LkV5L8VHf/6ZNN3WKsTxnovrG7D3f34f3796+6PADYE1YKdlV9QzZi/YHuvn0Z/vITp7qXy0eW8eNJLt1094NJHlrl+QHgQrHKu8QryU1J7uvuX9h0051JjizXjyS5Y9P465d3i1+R5LEnTp0DAE9u3wr3fVmSn0jymar61DL2r5LckOTDVXVtkgeSvHa57SNJXpnkWJKvJXnDCs8NABeUbQe7u/9btn5dOkmu3GJ+J3njdp8PAC5kPukMAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGECwAWAAwQaAAQQbAAYQbAAYQLABYIBt/z5sNhy6/td29PG+eMOrdvTxANgbHGEDwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAPunsPOOT0wDYiiNsABhAsAFgAKfE97idPsWeOM0OsA6OsAFgAMEGgAEEGwAG8Bo258yPngHsPkfYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAzgg1NYOx/EAnBmjrABYABH2Ow5fqUosBc5wgaAAXb9CLuqXpHk3UmeluR93X3Dbq8BzpXX2YF129Uj7Kp6WpL/lOSqJC9K8mNV9aLdXAMATLTbR9gvTXKsu7+QJFX1oSRXJ/n8Lq8D1soRO3CudjvYB5I8uGn7eJLv2+U1wJ7zVLzRbqddiN9U+Mbs/DP5a7Lbwa4txvrrJlRdl+S6ZfN/V9UfPuWr2lmXJPmTdS9iB9mf89uY/al3nNW0MftzlnZ0f87yv+FTxddmC0/R1+RvbzW428E+nuTSTdsHkzy0eUJ335jkxt1c1E6qqqPdfXjd69gp9uf8Zn/Ob3tpf/bSviQz92e3f6zr95JcVlXPr6qnJ3ldkjt3eQ0AMM6uHmF39+NV9c+T/GY2fqzr5u7+3G6uAQAm2vWfw+7ujyT5yG4/7y4aezr/NOzP+c3+nN/20v7spX1JBu5PdfeZZwEAa+WjSQFgAMHeIVV1aVV9rKruq6rPVdWb172mVVXV06rqk1X1q+tey6qq6qKquq2q/mD5Gn3/ute0iqr66eXv2Wer6oNV9U3rXtO5qKqbq+qRqvrsprFnVdVdVXX/cnnxOtd4Lk6zP/9h+fv26ar6L1V10TrXeC622p9Nt/2LquqqumQda9uO0+1PVb2pqv5w+bf079e1vrMl2Dvn8SQ/093fneSKJG/cAx+7+uYk9617ETvk3Ul+o7u/K8n3ZvB+VdWBJD+Z5HB3vzgbb+B83XpXdc5uSfKKk8auT3J3d1+W5O5le4pbcur+3JXkxd39PUn+Z5K37faiVnBLTt2fVNWlSX44yQO7vaAV3ZKT9qeq/lE2Pmnze7r77yb5+TWs65wI9g7p7oe7+xPL9T/LRhAOrHdV21dVB5O8Ksn71r2WVVXVtyb5wSQ3JUl3/2V3f3W9q1rZviTfXFX7kjwjJ32ewfmuu383yaMnDV+d5Nbl+q1JrtnVRa1gq/3p7t/q7seXzY9n43MnRjjN1ydJ3pnkZ3PSB16d706zP/8syQ3d/X+WOY/s+sLOkWA/BarqUJKXJLlnvStZybuy8Q/zr9e9kB3wgiQnkvzicor/fVX1zHUvaru6+0vZOBp4IMnDSR7r7t9a76p2xHO7++Fk4xvgJM9Z83p20j9N8uvrXsQqquo1Sb7U3b+/7rXskO9M8g+q6p6q+p2q+vvrXtCZCPYOq6pvSfIrSX6qu/903evZjqp6dZJHuvveda9lh+xLcnmS93b3S5L8eWadbv06y2u7Vyd5fpJvT/LMqvon610Vp1NVP5eNl8w+sO61bFdVPSPJzyX51+teyw7al+TibLyE+S+TfLiqtvr47POGYO+gqvqGbMT6A919+7rXs4KXJXlNVX0xyYeSvLyqfmm9S1rJ8STHu/uJMx63ZSPgU/1Qkj/q7hPd/VdJbk/yA2te0074clU9L0mWy/P+FOWZVNWRJK9O8uM9+2do/042vkH8/eX/CweTfKKq/tZaV7Wa40lu7w3/IxtnE8/rN9IJ9g5ZvjO7Kcl93f0L617PKrr7bd19sLsPZePNTB/t7rFHcN39x0kerKoXLkNXZvavdH0gyRVV9Yzl792VGfwmuk3uTHJkuX4kyR1rXMvKquoVSd6a5DXd/bV1r2cV3f2Z7n5Odx9a/r9wPMnly7+tqf5rkpcnSVV9Z5Kn5zz/5SaCvXNeluQnsnE0+qnlzyvXvSj+nzcl+UBVfTrJ30vy79a8nm1bzhTcluQTST6TjX/Hoz61qao+mOS/J3lhVR2vqmuT3JDkh6vq/my8E/mGda7xXJxmf/5jkr+Z5K7l/wf/ea2LPAen2Z+xTrM/Nyd5wfKjXh9KcuR8Pwvik84AYABH2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAM8H8BQcwKzJLjJ5EAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light",
"tags": []
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(figsize=(8, 6))\n",
"\n",
"ax.hist(weights, bins=20);"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "cPQODKAj5xak"
},
"source": [
"Now, the main model"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "ItxzfYLM5xam"
},
"outputs": [],
"source": [
"y = nhefs.wt82_71\n",
"X = nhefs[['constant', 'qsmk']]"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "DxGUOjry5xas"
},
"source": [
"Weighted least squares gives the right coefficients, but the standard error is off."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "IP-ZN6dn5xat",
"outputId": "07a13ad4-8b01-4c94-c6ba-291e14656b20"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" coef std err t P>|t| [0.025 0.975] \n",
" \n",
"\n",
" constant 1.7800 0.288 6.175 0.000 1.215 2.345 \n",
" \n",
"\n",
" qsmk 3.4405 0.408 8.434 0.000 2.640 4.241 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 17,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"wls = sm.WLS(y, X, weights=weights) \n",
"res = wls.fit()\n",
"res.summary().tables[1]"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "UqvNDOJJ5xay"
},
"source": [
"GEE gives the right coefficients and better standard errors"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "vx6pV6_g5xaz",
"outputId": "66b51bd6-efb4-4b9c-faa0-f6b324b6f66f"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" coef std err z P>|z| [0.025 0.975] \n",
" \n",
"\n",
" constant 1.7800 0.225 7.920 0.000 1.340 2.220 \n",
" \n",
"\n",
" qsmk 3.4405 0.525 6.547 0.000 2.411 4.470 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 18,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"gee = sm.GEE(\n",
" nhefs.wt82_71,\n",
" nhefs[['constant', 'qsmk']],\n",
" groups=nhefs.seqn,\n",
" weights=weights\n",
")\n",
"res = gee.fit()\n",
"res.summary().tables[1]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "agP0AOsL5xa3",
"outputId": "8a12a01a-fc4a-4657-81f7-c998278f70ab"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" estimate 95% C.I.\n",
"theta_1 3.44 (2.4, 4.5)\n"
]
}
],
"source": [
"est = res.params.qsmk\n",
"conf_ints = res.conf_int(alpha=0.05, cols=None)\n",
"lo, hi = conf_ints[0]['qsmk'], conf_ints[1]['qsmk']\n",
"\n",
"print(' estimate 95% C.I.')\n",
"print('theta_1 {:>6.2f} ({:>0.1f}, {:>0.1f})'.format(est, lo, hi))"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "sNbP433q5xa6"
},
"source": [
"Here's a simple check that there is no association between `sex` and `qsmk`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "ELtHxkw35xa8",
"outputId": "c733a677-c6d3-41bc-c9d3-5f0eebec1277"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" qsmk \n",
" 0 \n",
" 1 \n",
" \n",
" \n",
" sex \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 763.607760 \n",
" 763.623497 \n",
" \n",
" \n",
" 1 \n",
" 801.748892 \n",
" 797.200691 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {
"needs_background": "light",
"tags": []
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(figsize=(8, 6))\n",
"\n",
"ax.axhline(0, c='gray')\n",
"ax.plot(crosstab.index, crosstab[0], label='non-quitters')\n",
"ax.plot(crosstab.index, crosstab[1], label='quitters')\n",
"ax.set_xlabel('age', fontsize=14)\n",
"ax.set_ylabel('count', fontsize=14)\n",
"ax.legend(fontsize=12);"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "zfrhqjAk5xbS"
},
"source": [
"We see that there are actually a few ages with zero counts"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "hzpYvIfS5xbW",
"outputId": "a1ce4453-0399-4f75-bd6f-cf67b6f30fd8"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" qsmk \n",
" 0 \n",
" 1 \n",
" \n",
" \n",
" age \n",
" \n",
" \n",
" \n",
" \n",
" 63 \n",
" 3 \n",
" 3 \n",
" \n",
" \n",
" 64 \n",
" 7 \n",
" 1 \n",
" \n",
" \n",
" 65 \n",
" 3 \n",
" 2 \n",
" \n",
" \n",
" 66 \n",
" 4 \n",
" 0 \n",
" \n",
" \n",
" 67 \n",
" 2 \n",
" 0 \n",
" \n",
" \n",
" 69 \n",
" 6 \n",
" 2 \n",
" \n",
" \n",
" 70 \n",
" 2 \n",
" 1 \n",
" \n",
" \n",
" 71 \n",
" 0 \n",
" 1 \n",
" \n",
" \n",
" 72 \n",
" 2 \n",
" 2 \n",
" \n",
" \n",
" 74 \n",
" 0 \n",
" 1 \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
" coef std err z P>|z| [0.025 0.975] \n",
" \n",
"\n",
" constant 1.7800 0.225 7.920 0.000 1.340 2.220 \n",
" \n",
"\n",
" qsmk 3.4405 0.525 6.547 0.000 2.411 4.470 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 25,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"gee = sm.GEE(\n",
" nhefs.wt82_71,\n",
" nhefs[['constant', 'qsmk']],\n",
" groups=nhefs.seqn,\n",
" weights=(0.5 * weights)\n",
")\n",
"res = gee.fit()\n",
"res.summary().tables[1]"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "IEUBU0aG5xbq"
},
"source": [
"\"Second, we need to estimate Pr[A=1] for the numerator of the weights. We can obtain a nonparametric estimate by the ratio 403/1566 or, equivalently, by fitting a saturated logistic model for Pr[A=1] with an intercept and no covariates.\" pg 154"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "oKO38NWD5xbt"
},
"outputs": [],
"source": [
"qsmk = (nhefs.qsmk == 1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "bgHLFLPs5xbw",
"outputId": "586a2e6d-6be3-4bfb-b784-1a0a7f5ecb80"
},
"outputs": [
{
"data": {
"text/plain": [
"0.25734355044699875"
]
},
"execution_count": 27,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"# option 1\n",
"qsmk_mean = qsmk.mean()\n",
"qsmk_mean"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "6KyEQU6M5xb1",
"outputId": "241340ee-5265-4add-825e-336c55da44e2"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.570260\n",
" Iterations 5\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
" coef std err z P>|z| [0.025 0.975] \n",
" \n",
"\n",
" constant -1.0598 0.058 -18.335 0.000 -1.173 -0.947 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 28,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"# option 2\n",
"lgt = sm.Logit(qsmk, nhefs.constant)\n",
"res = lgt.fit()\n",
"res.summary().tables[1]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "w3rHN59b5xb4"
},
"outputs": [],
"source": [
"lgt_pred = res.predict()"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "wtVq54Qb5xb8"
},
"source": [
"Check for equivalence"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "tKkKUP-g5xb9",
"outputId": "5b13d757-220c-41bf-98f8-ce1b8356b4e0"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"equivalent: True\n"
]
}
],
"source": [
"equivalent = np.all(np.isclose(lgt_pred, qsmk_mean))\n",
"print('equivalent: {}'.format(equivalent))"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "SQq-mxqX5xcD"
},
"source": [
"### Program 12.3"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "GvxxebIR5xcG"
},
"source": [
"Create stabilized IP weights. Shortcut: modify the IP weights already calculated."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "nesdlHrc5xcH"
},
"outputs": [],
"source": [
"s_weights = np.zeros(nhefs.shape[0])\n",
"s_weights[qsmk] = qsmk.mean() * weights[qsmk] # `qsmk` was defined a few cells ago\n",
"s_weights[~qsmk] = (1 - qsmk).mean() * weights[~qsmk]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "urqAaeDO5xcS",
"outputId": "c89dde16-69a9-4833-d44b-aabb5cf91b02"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Stabilized weights\n",
" min mean max\n",
"------------------\n",
"0.33 1.00 4.30\n"
]
}
],
"source": [
"print('Stabilized weights')\n",
"print(' min mean max')\n",
"print('------------------')\n",
"print('{:>04.2f} {:>04.2f} {:>04.2f}'.format(\n",
" s_weights.min(),\n",
" s_weights.mean(),\n",
" s_weights.max()\n",
"))"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "weZHsmvB5xcd"
},
"source": [
"Refit the model from the last section, using the new weights"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "Puz1UrKF5xcg",
"outputId": "398d879f-1a3e-4af3-c150-b9ac79f0041b"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" coef std err z P>|z| [0.025 0.975] \n",
" \n",
"\n",
" constant 1.7800 0.225 7.920 0.000 1.340 2.220 \n",
" \n",
"\n",
" qsmk 3.4405 0.525 6.547 0.000 2.411 4.470 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 33,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"gee = sm.GEE(\n",
" nhefs.wt82_71,\n",
" nhefs[['constant', 'qsmk']],\n",
" groups=nhefs.seqn,\n",
" weights=s_weights\n",
")\n",
"res = gee.fit()\n",
"res.summary().tables[1]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "JhW0TcYh5xck",
"outputId": "329116b8-f776-4e18-d2f0-f4c031eaeb63"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" estimate 95% C.I.\n",
"theta_1 3.44 (2.4, 4.5)\n"
]
}
],
"source": [
"est = res.params.qsmk\n",
"conf_ints = res.conf_int(alpha=0.05, cols=None)\n",
"lo, hi = conf_ints[0]['qsmk'], conf_ints[1]['qsmk']\n",
"\n",
"print(' estimate 95% C.I.')\n",
"print('theta_1 {:>6.2f} ({:>0.1f}, {:>0.1f})'.format(est, lo, hi))"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "sNF5CLEZ5xct"
},
"source": [
"The estimate is the same as in the previous section"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "Ybob4PS75xcx"
},
"source": [
"We can check again for no association between sex and qsmk in the the pseudo-population"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "GONEMTsI5xcy",
"outputId": "7bebc1f2-972f-494f-f9ed-c6a8ea090a4a"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" qsmk \n",
" 0 \n",
" 1 \n",
" \n",
" \n",
" sex \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 567.098228 \n",
" 196.513582 \n",
" \n",
" \n",
" 1 \n",
" 595.423986 \n",
" 205.154456 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {
"needs_background": "light",
"tags": []
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(figsize=(8, 6))\n",
"\n",
"A.hist(bins=30, ax=ax);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "HxX3tTUy5xds",
"outputId": "ff96fd89-2f3b-4a55-e65f-f84128dcad6e"
},
"outputs": [
{
"data": {
"text/plain": [
"(-2.057659208261618, 10.467830908151853)"
]
},
"execution_count": 43,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"A.mean(), A.std()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "YomIJH0W5xdw"
},
"outputs": [],
"source": [
"fA = scipy.stats.norm.pdf(A, A.mean(), A.std())"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "-1VHCB3y5xdy"
},
"source": [
"Then the stabilized IP weights are"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "PhhTqmii5xd1"
},
"outputs": [],
"source": [
"sw = fA / fAL"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "u-rvdWri5xd4",
"outputId": "8fda2d78-2ea1-4695-8029-ac86562dd5ab"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Stabilized weights\n",
" min mean max\n",
"------------------\n",
"0.19 1.00 5.10\n"
]
}
],
"source": [
"print('Stabilized weights')\n",
"print(' min mean max')\n",
"print('------------------')\n",
"print('{:>04.2f} {:>04.2f} {:>04.2f}'.format(\n",
" sw.min(),\n",
" sw.mean(),\n",
" sw.max()\n",
"))"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "GPJoQBsF5xd7"
},
"source": [
"Now fit the marginal structural model"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "YEbqpvtN5xd8"
},
"outputs": [],
"source": [
"y = intensity25.wt82_71\n",
"X = pd.DataFrame(OrderedDict((\n",
" ('constant', np.ones(y.shape[0])),\n",
" ('A', A),\n",
" ('A^2', A**2)\n",
")))\n",
"\n",
"model = sm.GEE(\n",
" y,\n",
" X,\n",
" groups=intensity25.seqn,\n",
" weights=sw\n",
")\n",
"res = model.fit()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "cSj00djn5xeC",
"outputId": "9352f6d4-ba8c-44f4-c0eb-131f85c3c354"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" coef std err z P>|z| [0.025 0.975] \n",
" \n",
"\n",
" constant 2.0045 0.295 6.792 0.000 1.426 2.583 \n",
" \n",
"\n",
" A -0.1090 0.032 -3.456 0.001 -0.171 -0.047 \n",
" \n",
"\n",
" A^2 0.0027 0.002 1.115 0.265 -0.002 0.007 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 48,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"res.summary().tables[1]"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "T9A2qnV75xeH"
},
"source": [
"To get the estimate and confidence interval for \"no change\", you can read off the values in the `constant` row above (because `A` and `A^2` will be zero).\n",
"\n",
"Getting Statmodels to calculate the estimate and confidence interval for when smoking increases by 20 cigarettes / day will take a couple extra steps. In Chapter 11, the regression result had a `get_prediction` method. The GEE result doesn't (yet?) have that _method_, so we'll use the hidden `get_prediction` _function_."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "0ZGG4gfc5xeI"
},
"outputs": [],
"source": [
"from statsmodels.regression._prediction import get_prediction"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "lCvR9WSD5xeL",
"outputId": "a6131514-35c0-4a03-be4c-9a177e5a1e22"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" mean_ci_lower \n",
" mean_ci_upper \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 2.0 \n",
" 1.4 \n",
" 2.6 \n",
" \n",
" \n",
" 1 \n",
" 0.9 \n",
" -1.7 \n",
" 3.5 \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" estimate \n",
" CI lower \n",
" CI upper \n",
" \n",
" \n",
" \n",
" \n",
" no change \n",
" 2.0 \n",
" 1.4 \n",
" 2.6 \n",
" \n",
" \n",
" +20 per day \n",
" 0.9 \n",
" -1.7 \n",
" 3.5 \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
" coef std err z P>|z| [0.025 0.975] \n",
" \n",
"\n",
" constant -1.4905 0.079 -18.881 0.000 -1.645 -1.336 \n",
" \n",
"\n",
" qsmk 0.0301 0.157 0.191 0.848 -0.278 0.338 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 52,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"model = sm.GEE(\n",
" nhefs.death,\n",
" nhefs[['constant', 'qsmk']],\n",
" groups=nhefs.seqn,\n",
" weights=s_weights,\n",
" family=sm.families.Binomial()\n",
")\n",
"res = model.fit()\n",
"res.summary().tables[1]"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "Os6FoWkv5xeb"
},
"source": [
"Odd ratio is $\\exp(\\hat{\\theta}_1)$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "X0p6ex9V5xeb",
"outputId": "e7bcf040-3f38-4a23-9268-abc3223d52a5"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" estimate 95% C.I.\n",
"odds ratio 1.03 (0.8, 1.4)\n"
]
}
],
"source": [
"est = np.exp(res.params.qsmk)\n",
"conf_ints = res.conf_int(alpha=0.05, cols=None)\n",
"lo = np.exp(conf_ints[0]['qsmk'])\n",
"hi = np.exp(conf_ints[1]['qsmk'])\n",
"\n",
"print(' estimate 95% C.I.')\n",
"print('odds ratio {:>6.2f} ({:>0.1f}, {:>0.1f})'.format(est, lo, hi))"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "0yompze_5xef"
},
"source": [
"### Program 12.6"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "1OMDoNfi5xeg"
},
"source": [
" Create the numerator of the IP weights. Reuse the basic `weights` for the denominator."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "8ecQsdoq5xeh",
"outputId": "1a1b0611-62df-4197-f2a1-2e6632f8efbd"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.567819\n",
" Iterations 5\n"
]
}
],
"source": [
"numer = logit_ip_f(nhefs.qsmk, nhefs[['constant', 'sex']])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "4B4afU4i5xej"
},
"outputs": [],
"source": [
"sw_AV = numer * weights"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "GAP0se4f5xem",
"outputId": "05d35726-5494-4603-819b-bb8c6b814c55"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Stabilized weights\n",
" min mean max\n",
"------------------\n",
"0.29 1.00 3.80\n"
]
}
],
"source": [
"print('Stabilized weights')\n",
"print(' min mean max')\n",
"print('------------------')\n",
"print('{:>04.2f} {:>04.2f} {:>04.2f}'.format(\n",
" sw_AV.min(),\n",
" sw_AV.mean(),\n",
" sw_AV.max()\n",
"))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "2HOn9xu65xeq",
"outputId": "e60ad35f-cd75-40b1-f735-5000cc5e582e"
},
"outputs": [
{
"data": {
"text/plain": [
"(1566, 83)"
]
},
"execution_count": 57,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"nhefs.shape"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "OoonBszO5xet",
"outputId": "45be7bdb-4f3a-4667-96c3-6ce4c9b9802b"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" coef std err z P>|z| [0.025 0.975] \n",
" \n",
"\n",
" constant 1.7844 0.310 5.752 0.000 1.176 2.393 \n",
" \n",
"\n",
" qsmk 3.5220 0.658 5.353 0.000 2.232 4.811 \n",
" \n",
"\n",
" sex -0.0087 0.449 -0.019 0.985 -0.890 0.872 \n",
" \n",
"\n",
" qsmk_and_female -0.1595 1.047 -0.152 0.879 -2.212 1.893 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 58,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"nhefs['qsmk_and_female'] = nhefs.qsmk * nhefs.sex\n",
"\n",
"model = sm.WLS(\n",
" nhefs.wt82_71,\n",
" nhefs[['constant', 'qsmk', 'sex', 'qsmk_and_female']],\n",
" weights=sw_AV\n",
")\n",
"res = model.fit(cov_type='cluster', cov_kwds={'groups': nhefs.seqn})\n",
"res.summary().tables[1]"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "SXgrOVwk5xew"
},
"source": [
"### Program 12.7"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "k0W7xmhC5xew"
},
"source": [
"We're going back to the original dataset"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "okk5Yj4F5xex",
"outputId": "7c498741-23f6-4be2-d417-244907e8b838"
},
"outputs": [
{
"data": {
"text/plain": [
"(1629, 64)"
]
},
"execution_count": 59,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"nhefs_all.shape"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "dtA_wcuV5xez"
},
"source": [
"We'll add features that were added to the reduced dataset that we've been using"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "z_b3LGhB5xez"
},
"source": [
"Add constant feature"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "AG64cc9c5xe4"
},
"outputs": [],
"source": [
"nhefs_all['constant'] = 1"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "3maWP6Ff5xe_"
},
"source": [
"Add dummy features"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "UcZAADB25xfA"
},
"outputs": [],
"source": [
"edu_dummies = pd.get_dummies(nhefs_all.education, prefix='edu')\n",
"exercise_dummies = pd.get_dummies(nhefs_all.exercise, prefix='exercise')\n",
"active_dummies = pd.get_dummies(nhefs_all.active, prefix='active')\n",
"\n",
"nhefs_all = pd.concat(\n",
" [nhefs_all, edu_dummies, exercise_dummies, active_dummies],\n",
" axis=1\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "SrToJKmO5xfD"
},
"source": [
"Add squared features"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "pbH7r8VZ5xfD"
},
"outputs": [],
"source": [
"for col in ['age', 'wt71', 'smokeintensity', 'smokeyrs']:\n",
" nhefs_all['{}^2'.format(col)] = nhefs_all[col] * nhefs_all[col]"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "i9skhn4h5xfF"
},
"source": [
"We'll also add a feature to track censored individuals"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "Se6bTv8b5xfG"
},
"outputs": [],
"source": [
"nhefs_all['censored'] = nhefs_all.wt82.isnull().astype('int')"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "uZ-MSFep5xfH"
},
"source": [
"Create the IP weights for treatment"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "F-5WskL95xfI"
},
"outputs": [],
"source": [
"X_ip = nhefs_all[[\n",
" 'constant', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
" 'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
" 'age', 'age^2', 'wt71', 'wt71^2',\n",
" 'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2'\n",
"]]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "2hanuMb45xfJ",
"outputId": "ae5093ce-9c83-48c6-c93e-882807c7cb2f",
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.542264\n",
" Iterations 6\n"
]
}
],
"source": [
"ip_denom = logit_ip_f(nhefs_all.qsmk, X_ip)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "KHYrzGOa5xfL",
"outputId": "d4740bfc-2718-4941-abcc-6ecd9fc967b5"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.575901\n",
" Iterations 5\n"
]
}
],
"source": [
"ip_numer = logit_ip_f(nhefs_all.qsmk, nhefs_all.constant)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "LFDS2Ji-5xfM"
},
"outputs": [],
"source": [
"sw_A = ip_numer / ip_denom"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "0D2Ke2Jl5xfO",
"outputId": "747db8e1-6255-4477-ab1b-7ab4d82bfb7b"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Stabilized weights\n",
" min mean max\n",
"------------------\n",
"0.33 1.00 4.21\n"
]
}
],
"source": [
"print('Stabilized weights')\n",
"print(' min mean max')\n",
"print('------------------')\n",
"print('{:>04.2f} {:>04.2f} {:>04.2f}'.format(\n",
" sw_A.min(),\n",
" sw_A.mean(),\n",
" sw_A.max()\n",
"))"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "9F-eAu7F5xfQ"
},
"source": [
"Now the IP weights for censoring"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "x-RYy7xJ5xfQ"
},
"outputs": [],
"source": [
"# same as previous, but with 'qsmk' added\n",
"\n",
"X_ip = nhefs_all[[\n",
" 'constant', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
" 'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
" 'age', 'age^2', 'wt71', 'wt71^2',\n",
" 'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2',\n",
" 'qsmk'\n",
"]]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "7UcdQP6-5xfT",
"outputId": "bab9068a-b066-4f57-e302-cdf84bb9af2d"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.142836\n",
" Iterations 8\n"
]
}
],
"source": [
"ip_denom = logit_ip_f(nhefs_all.censored, X_ip)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "8BGi1lpi5xfb",
"outputId": "13adef71-3a02-4dd3-c612-eee5280dce14"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.161989\n",
" Iterations 7\n"
]
}
],
"source": [
"ip_numer = logit_ip_f(\n",
" nhefs_all.censored,\n",
" nhefs_all[['constant', 'qsmk']]\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "Ep-9U1sK5xff"
},
"outputs": [],
"source": [
"sw_C = ip_numer / ip_denom\n",
"sw_C[nhefs_all.censored == 1] = 1"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "Yoq1Epg65xfk",
"outputId": "e1863baa-1b2e-46c4-b291-63d26cc88331"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Stabilized weights\n",
" min mean max\n",
"------------------\n",
"0.94 1.00 1.72\n"
]
}
],
"source": [
"print('Stabilized weights')\n",
"print(' min mean max')\n",
"print('------------------')\n",
"print('{:>04.2f} {:>04.2f} {:>04.2f}'.format(\n",
" sw_C.min(),\n",
" sw_C.mean(),\n",
" sw_C.max()\n",
"))"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "nsSj6SHV5xfr"
},
"source": [
"Now create the combined IP weights"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "T1pVFskw5xfs"
},
"outputs": [],
"source": [
"sw_AC = sw_A * sw_C"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "fECz2aYA5xft",
"outputId": "dd7799ac-b885-4812-ec42-ffb3b25b4cb6"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Stabilized weights\n",
" min mean max\n",
"------------------\n",
"0.35 1.00 4.09\n"
]
}
],
"source": [
"print('Stabilized weights')\n",
"print(' min mean max')\n",
"print('------------------')\n",
"print('{:>04.2f} {:>04.2f} {:>04.2f}'.format(\n",
" sw_AC.min(),\n",
" sw_AC.mean(),\n",
" sw_AC.max()\n",
"))"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab_type": "text",
"id": "AjMat9Ut5xfv"
},
"source": [
"Now model weight gain using the combined IP weights"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "0ze9Z_rf5xfv"
},
"outputs": [],
"source": [
"wls = sm.WLS(\n",
" nhefs.wt82_71,\n",
" nhefs[['constant', 'qsmk']],\n",
" weights=sw_AC[nhefs_all.censored == 0]\n",
") \n",
"res = wls.fit(cov_type='cluster', cov_kwds={'groups': nhefs.seqn})"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "NNVvYz6l5xfx",
"outputId": "9e540919-8f20-4db7-b5a9-3addfd8edb27",
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" coef std err z P>|z| [0.025 0.975] \n",
" \n",
"\n",
" constant 1.6620 0.233 7.136 0.000 1.206 2.118 \n",
" \n",
"\n",
" qsmk 3.4965 0.526 6.648 0.000 2.466 4.527 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 77,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"res.summary().tables[1]"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "true",
"colab_type": "text",
"id": "_fOBcWPg5xfz"
},
"source": [
"## 一个类组织本章代码\n",
"\n",
"我们把本章主要代码工作整理成类。"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "HPGd4e_b5xf0"
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"import warnings\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"from collections import OrderedDict\n",
"import numpy as np\n",
"import pandas as pd\n",
"import statsmodels.api as sm\n",
"import scipy.stats\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"Collapsed": "false",
"code_folding": [
0,
1,
6,
28,
49,
65,
83
],
"colab": {},
"colab_type": "code",
"id": "ZpHVrt6D5xf3",
"outputId": "b36c06cb-b4b0-4a51-a213-12ec58bb9573"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING *** OLE2 inconsistency: SSCS size is 0 but SSAT size is non-zero\n",
"Optimization terminated successfully.\n",
" Current function value: 0.567819\n",
" Iterations 5\n",
"Confounders for IP Weight is:\n",
" ['constant', 'sex', 'race', 'age', 'age^2', 'edu_2', 'edu_3', 'edu_4', 'edu_5', 'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2', 'exercise_1', 'exercise_2', 'active_1', 'active_2', 'wt71', 'wt71^2']\n",
"\n",
"Optimization terminated successfully.\n",
" Current function value: 0.535408\n",
" Iterations 6\n",
"Marignal Structural Models 的有关变量:\n",
" ['constant', 'qsmk', 'sex', 'qsmk_and_female']\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
" coef std err z P>|z| [0.025 0.975] \n",
" \n",
"\n",
" constant 1.7844 0.310 5.752 0.000 1.176 2.393 \n",
" \n",
"\n",
" qsmk 3.5220 0.658 5.353 0.000 2.232 4.811 \n",
" \n",
"\n",
" sex -0.0087 0.449 -0.019 0.985 -0.890 0.872 \n",
" \n",
"\n",
" qsmk_and_female -0.1595 1.047 -0.152 0.879 -2.212 1.893 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"class IPW(object):\n",
" def __init__(self):\n",
" self.help = \"This is a demonstration for IP Weighting!\"\n",
" self.data = self.get_data()\n",
" self.ipw_cols = None\n",
" \n",
" def get_data(self):\n",
" nhefs_all = pd.read_excel('NHEFS.xls')\n",
" restriction_cols = ['sex', 'age', 'race', 'wt82', 'ht', \n",
" 'school', 'alcoholpy', 'smokeintensity']\n",
" missing = nhefs_all[restriction_cols].isnull().any(axis=1)\n",
" nhefs = nhefs_all.loc[~missing]\n",
" nhefs['constant'] = 1\n",
" nhefs['university'] = (nhefs.education == 5).astype('int')\n",
" nhefs['inactive'] = (nhefs.active == 2).astype('int')\n",
" nhefs['no_exercise'] = (nhefs.exercise == 2).astype('int') \n",
" for col in ['age', 'wt71', 'smokeintensity', 'smokeyrs']:\n",
" nhefs['{}^2'.format(col)] = nhefs[col] * nhefs[col] \n",
" edu_dummies = pd.get_dummies(nhefs.education, prefix='edu')\n",
" exercise_dummies = pd.get_dummies(nhefs.exercise, prefix='exercise')\n",
" active_dummies = pd.get_dummies(nhefs.active, prefix='active')\n",
"\n",
" nhefs = pd.concat(\n",
" [nhefs, edu_dummies, exercise_dummies, active_dummies],\n",
" axis=1\n",
" )\n",
" return nhefs\n",
"\n",
" def _logit_ip_f(self, y, X):\n",
" model = sm.Logit(y, X)\n",
" res = model.fit()\n",
" weights = np.zeros(X.shape[0])\n",
" weights[y == 1] = res.predict(X.loc[y == 1])\n",
" weights[y == 0] = (1 - res.predict(X.loc[y == 0]))\n",
" return weights\n",
" \n",
" def treatment_probs(self, cols=None):\n",
" if cols is None:\n",
" cols = ['constant', 'sex', 'race', 'age', 'age^2', 'edu_2', 'edu_3', 'edu_4', 'edu_5',\n",
" 'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2', \n",
" 'exercise_1', 'exercise_2', 'active_1', 'active_2', 'wt71', 'wt71^2']\n",
" self.ipw_cols = cols\n",
" print('Confounders for IP Weight is:\\n', cols)\n",
" print()\n",
" \n",
" X_ip = self.data[cols] \n",
" \n",
" return self._logit_ip_f(self.data.qsmk, X_ip) \n",
" \n",
" def ipw(self):\n",
" nhefs = self.data\n",
" denoms = self.treatment_probs(self.ipw_cols)\n",
" weights = 1 / denoms\n",
" assert(weights.shape[0] == self.data.shape[0])\n",
" print('Min and Max of weights:', weights.min(), weights.max())\n",
" \n",
" gee = sm.GEE(\n",
" nhefs.wt82_71,\n",
" nhefs[['constant', 'qsmk']],\n",
" groups=nhefs.seqn,\n",
" weights=weights\n",
" )\n",
" res = gee.fit()\n",
" return res.summary().tables[1]\n",
" \n",
" def stabled_ipw(self):\n",
" nhefs = self.data\n",
" qsmk = (nhefs.qsmk == 1)\n",
" denoms = self.treatment_probs(self.ipw_cols)\n",
" weights = 1 / denoms\n",
" s_weights = np.zeros(nhefs.shape[0])\n",
" s_weights[qsmk] = qsmk.mean() * weights[qsmk] # `qsmk` was defined a few cells ago\n",
" s_weights[~qsmk] = (1 - qsmk).mean() * weights[~qsmk]\n",
" print(qsmk.mean())\n",
" gee = sm.GEE(\n",
" nhefs.wt82_71,\n",
" nhefs[['constant', 'qsmk']],\n",
" groups=nhefs.seqn,\n",
" weights=s_weights\n",
" )\n",
" res = gee.fit()\n",
" return res.summary().tables[1]\n",
"\n",
" def marginal_structural_models(self):\n",
" nhefs = self.data\n",
" numer = self._logit_ip_f(nhefs.qsmk, nhefs[['constant', 'sex']])\n",
" weights = 1 / self.treatment_probs()\n",
" sw_AV = numer * weights\n",
" nhefs['qsmk_and_female'] = nhefs.qsmk * nhefs.sex\n",
" model = sm.WLS(\n",
" nhefs.wt82_71,\n",
" nhefs[['constant', 'qsmk', 'sex', 'qsmk_and_female']],\n",
" weights=sw_AV\n",
" )\n",
" res = model.fit(cov_type='cluster', cov_kwds={'groups': nhefs.seqn})\n",
" print('Marignal Structural Models 的有关变量:\\n', ['constant', 'qsmk', 'sex', 'qsmk_and_female'] )\n",
" return res.summary().tables[1]\n",
"\n",
"g = IPW()\n",
"# g.data.shape\n",
"# g.ipw_cols = ['constant', 'sex', 'race', 'age']\n",
"# g.ipw()\n",
"# g.stabled_ipw()\n",
"g.marginal_structural_models()"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false",
"colab": {},
"colab_type": "code",
"id": "yzpPpKRv5xf5"
},
"source": [
"Table of Contents\n",
"\n",
"- IPW via modeling\n",
"- Stablized IPW\n",
"- Marginal structural models\n",
"- Effect modification and marginal structural models\n",
"- Censoring and missing data"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING *** OLE2 inconsistency: SSCS size is 0 but SSAT size is non-zero\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" qsmk \n",
" wt82_71 \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 0 \n",
" -10.093960 \n",
" \n",
" \n",
" 1 \n",
" 0 \n",
" 2.604970 \n",
" \n",
" \n",
" 2 \n",
" 0 \n",
" 9.414486 \n",
" \n",
" \n",
" 3 \n",
" 0 \n",
" 4.990117 \n",
" \n",
" \n",
" 4 \n",
" 0 \n",
" 4.989251 \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
" coef std err t P>|t| [0.025 0.975] \n",
" \n",
"\n",
" constant 1.9845 0.229 8.672 0.000 1.536 2.433 \n",
" \n",
"\n",
" qsmk 2.5406 0.451 5.632 0.000 1.656 3.425 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# 用线性回顾估计因果效应的\n",
"ols = sm.OLS(nhefs.wt82_71, nhefs[['constant', 'qsmk']])\n",
"res = ols.fit()\n",
"res.summary().tables[1]"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Confounders for IP Weight is:\n",
" ['constant', 'sex', 'race', 'age', 'age^2', 'edu_2', 'edu_3', 'edu_4', 'edu_5', 'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2', 'exercise_1', 'exercise_2', 'active_1', 'active_2', 'wt71', 'wt71^2']\n",
"\n",
"Optimization terminated successfully.\n",
" Current function value: 0.535408\n",
" Iterations 6\n",
"Min and Max of weights: 1.0537416280320455 16.700094350727504\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
" coef std err z P>|z| [0.025 0.975] \n",
" \n",
"\n",
" constant 1.7800 0.225 7.920 0.000 1.340 2.220 \n",
" \n",
"\n",
" qsmk 3.4405 0.525 6.547 0.000 2.411 4.470 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# IP Weight 方法\n",
"g.ipw()"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Confounders for IP Weight is:\n",
" ['constant', 'sex', 'race', 'age', 'age^2', 'edu_2', 'edu_3', 'edu_4', 'edu_5', 'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2', 'exercise_1', 'exercise_2', 'active_1', 'active_2', 'wt71', 'wt71^2']\n",
"\n",
"Optimization terminated successfully.\n",
" Current function value: 0.535408\n",
" Iterations 6\n",
"0.25734355044699875\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
" coef std err z P>|z| [0.025 0.975] \n",
" \n",
"\n",
" constant 1.7800 0.225 7.920 0.000 1.340 2.220 \n",
" \n",
"\n",
" qsmk 3.4405 0.525 6.547 0.000 2.411 4.470 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Stablized IPW\n",
"g.stabled_ipw()"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.567819\n",
" Iterations 5\n",
"Confounders for IP Weight is:\n",
" ['constant', 'sex', 'race', 'age', 'age^2', 'edu_2', 'edu_3', 'edu_4', 'edu_5', 'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2', 'exercise_1', 'exercise_2', 'active_1', 'active_2', 'wt71', 'wt71^2']\n",
"\n",
"Optimization terminated successfully.\n",
" Current function value: 0.535408\n",
" Iterations 6\n",
"Marignal Structural Models 的有关变量:\n",
" ['constant', 'qsmk', 'sex', 'qsmk_and_female']\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
" coef std err z P>|z| [0.025 0.975] \n",
" \n",
"\n",
" constant 1.7844 0.310 5.752 0.000 1.176 2.393 \n",
" \n",
"\n",
" qsmk 3.5220 0.658 5.353 0.000 2.232 4.811 \n",
" \n",
"\n",
" sex -0.0087 0.449 -0.019 0.985 -0.890 0.872 \n",
" \n",
"\n",
" qsmk_and_female -0.1595 1.047 -0.152 0.879 -2.212 1.893 \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# marginal_structural_models\n",
"g.marginal_structural_models()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"colab": {
"name": "chapter12-IPW.ipynb",
"provenance": [],
"toc_visible": true
},
"kernelspec": {
"display_name": "Python 3",
"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.7.4"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": false,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": true
}
},
"nbformat": 4,
"nbformat_minor": 4
}