{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# DoWhy example on the Lalonde dataset\n", "\n", "Thanks to [@mizuy](https://github.com/mizuy) for providing this example. Here we use the Lalonde dataset and apply IPW estimator to it. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "R[write to console]: Loading required package: MASS\n", "\n", "R[write to console]: ## \n", "## Matching (Version 4.9-7, Build Date: 2020-02-05)\n", "## See http://sekhon.berkeley.edu/matching for additional documentation.\n", "## Please cite software as:\n", "## Jasjeet S. Sekhon. 2011. ``Multivariate and Propensity Score Matching\n", "## Software with Automated Balance Optimization: The Matching package for R.''\n", "## Journal of Statistical Software, 42(7): 1-52. \n", "##\n", "\n", "\n" ] }, { "data": { "text/plain": [ "array(['Matching', 'MASS', 'tools', 'stats', 'graphics', 'grDevices',\n", " 'utils', 'datasets', 'methods', 'base'], dtype='\n", "WLS Regression Results\n", "\n", " Dep. Variable: re78 R-squared: 0.015\n", "\n", "\n", " Model: WLS Adj. R-squared: 0.013\n", "\n", "\n", " Method: Least Squares F-statistic: 6.743\n", "\n", "\n", " Date: Sun, 06 Jun 2021 Prob (F-statistic): 0.00973\n", "\n", "\n", " Time: 19:11:54 Log-Likelihood: -4544.7\n", "\n", "\n", " No. Observations: 445 AIC: 9093.\n", "\n", "\n", " Df Residuals: 443 BIC: 9102.\n", "\n", "\n", " Df Model: 1 \n", "\n", "\n", " Covariance Type: nonrobust \n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
Intercept 4555.0759 406.704 11.200 0.000 3755.767 5354.385
treat[T.True] 1639.8076 631.496 2.597 0.010 398.708 2880.907
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 303.262 Durbin-Watson: 2.085
Prob(Omnibus): 0.000 Jarque-Bera (JB): 4770.581
Skew: 2.709 Prob(JB): 0.00
Kurtosis: 18.097 Cond. No. 2.47


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified." ], "text/plain": [ "\n", "\"\"\"\n", " WLS Regression Results \n", "==============================================================================\n", "Dep. Variable: re78 R-squared: 0.015\n", "Model: WLS Adj. R-squared: 0.013\n", "Method: Least Squares F-statistic: 6.743\n", "Date: Sun, 06 Jun 2021 Prob (F-statistic): 0.00973\n", "Time: 19:11:54 Log-Likelihood: -4544.7\n", "No. Observations: 445 AIC: 9093.\n", "Df Residuals: 443 BIC: 9102.\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "=================================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "---------------------------------------------------------------------------------\n", "Intercept 4555.0759 406.704 11.200 0.000 3755.767 5354.385\n", "treat[T.True] 1639.8076 631.496 2.597 0.010 398.708 2880.907\n", "==============================================================================\n", "Omnibus: 303.262 Durbin-Watson: 2.085\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 4770.581\n", "Skew: 2.709 Prob(JB): 0.00\n", "Kurtosis: 18.097 Cond. No. 2.47\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model=CausalModel(\n", " data = lalonde,\n", " treatment='treat',\n", " outcome='re78',\n", " common_causes='nodegr+black+hisp+age+educ+married'.split('+'))\n", "identified_estimand = model.identify_effect(proceed_when_unidentifiable=True)\n", "estimate = model.estimate_effect(identified_estimand,\n", " method_name=\"backdoor.propensity_score_weighting\",\n", " target_units=\"ate\", \n", " method_params={\"weighting_scheme\":\"ips_weight\"})\n", "#print(estimate)\n", "print(\"Causal Estimate is \" + str(estimate.value))\n", "\n", "import statsmodels.formula.api as smf\n", "reg=smf.wls('re78~1+treat', data=lalonde, weights=lalonde.ips_stabilized_weight)\n", "res=reg.fit()\n", "res.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Interpret the estimate\n", "The plot below shows how the distribution of a confounder, \"married\" changes from the original data to the weighted data. In both datasets, we compare the distribution of \"married\" across treated and untreated units." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAHwCAYAAAC7apkrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyKklEQVR4nO3de5hcZZXv8e8iAYIkco2cmBCDGjXgwSAhIEYmXAYQVFABwTMQECeCUQcRFRlHooNHdFSO4BFkhCGoKAhyiIoXrgIiYNCIQGCIEiAhQAwQQLmFrPPH3p286VR3upOurr58P8/TT1ftveuttXdVr/717reqIjORJEmSVNmg1QVIkiRJfYkBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKgzYgBwR50TEv/XQWGMj4pmIGFJfvz4iPtgTY9fj/TwipvXUeN2439Mi4q8R8Uhv33cjEfG2iLh3HW87NSIWdrI+I+K16zj2WyPivvo5cPC6jNEXdPV5uz6Pw7pYn8emGSLiroiY2sVtF0TEPs2taOCxP3fpfgdNf17LbSMi/isinoiI29ZljFay73ZNX+y7/TIg1wfn2Yh4OiKejIibI+K4iFi5P5l5XGb+exfH6vRAZ+aDmTk8M1/qgdpnRsT32o3/9syctb5jd7OOscAngO0z83/05n13JDNvzMzXt7qOBr4AfLN+Dvy/VhfTbM18HHo6vDRDZu6Qmdev7zjrEwr6M/vz+rM/r2YK8I/AmMyc3OgxGgjsu32v7/bLgFx7Z2aOAF4FnA58Gjivp+8kIob29Jh9xFhgaWY+1tt33OiY9vHj/CrgrnW5YR/fL6lZ7M/rx/68yquABZn5t54YbAA/Z9TTMrPffQELgH3aLZsMrADeWF+/ADitvrw18FPgSeBx4EaqPw6+W9/mWeAZ4FPAOCCBY4EHgRuKZUPr8a4HvgTcBjwFXAFsWa+bCixsVC+wP/AC8GJ9f38sxvtgfXkD4LPAA8BjwIXAZvW6tjqm1bX9FfjXTo7TZvXtl9TjfbYef596n1fUdVzQ4LZTgYX1MXkMWAwcDBwA/Hd9HE9pd/x/Wx/jxcA3gY2K9QnMAO4D7i/G/zTwSP1YrHbsgFcCl9X13w98rFi3Sf0YPwHcDXyy/XFvtz8JfAz4S33c/gPYoFj/AWBePd4vgVfVy//M6s+Rjeu6ZtfHYD7wz8U4M4FLge/Vz40P1o/DefVxWQScBgzpoM6uHMeG+wEcDfymvs0y4B5g7+K219f1bFTX/j+Lda8A/g6MbPA4LABOAu6ox70YGFas/1Rd68P1+Am8tsG+fRF4CXiuPpbfLPbpuPq58STwf4FY22PTYPxZwCfqy6PrcWfU119T73PbsXoHMLe+v5uBHRv1F6rn2az6vufV+7rWYwNsyuo/Y89QPW8mA3Pq58ajwNdb3U97+gv7s/25+/35G8BD9eN1O/C2evmxVP3ipfpY3NrBY9Rhj2VVXzwDWEr9vGvw/LTv2ndX37dWN9N1+aJBA66XPwgcX1++gFUN+EvAOcCG9dfb2p4I7cdiVZO7sD7Ym9C4AS8C3lhvcxnwvbJxdVQvVYD6Xrv117OqAX+AKnS9GhgO/Bj4brva/rOu603A88CEDo7ThVS/HEbUt/1v4NiO6mx326nAcuBz9TH7Z6pGeFE93g71E3G7evudgd2AofV9zQNOaNdgrgK2rGtvG//LVKGzbdnCevsNqBrl56gay6upmtN+9frTqX6RbglsC9y5lv1J4Lp6+7H1sWg75gfVx3xCXf9ngZs7er5R/VL+FtUP5MT6uOxVPL4vUv2y2qDer8uBb1M9V15B9Yv7Qx3U2ZXj2NF+HF0f04/Xj9n7qJpHWzi4vtj2W8CXi3H/BfhJo+dGvf+3UTWaLeuajqvX7U/1C3QH4GVUfxg0bNTta2i3Tz8FNq/3aQmwf1cem3bjfKDYh/dT/XFzcbHuivryTlShYldgCFWgWQBs3ODn9XTg18AWwBiqhtzVY7PacayX/RY4sr48HNit1f20p7+wP9ufu9+f/wnYqq7vE1Q9ZVi97mjgpmLbRo9Rhz2WVX3xo/X4mzS4f/uufXfNfevt5tkTX3TcgG+h/oud1RvwF6gaUaO/rlYbi1VN7tUNlpUN+PRi/fZUf9UO6eDBKR/4mXTegK8BPlysez1V4Gr7wU2quVht628DDm+wX0PqmrYvln0IuL6jJ1G720+larBtf4WPqO9712Kb24GDO7j9CcDl7X4Y92o3/gus/hfxypqofogebDfmZ4D/qi//hfqHub4+fS37k+22/zBwTX3559S/mOrrG1D9Vf+qBo/ftlR/jY8otv8S9Vme+vG9oVi3DdUvyU2KZUcA13Xxud7oOHa0H0dTnU0ozwLcxqrGUD7PdqUKLG1BZA5wWKPnRr3//1Rc/wpwTn35fOBLxbrXsm6Nekpx/RLg5K48Nu3GeQ3VGYcNqALXh4rn0yzgxPry2cC/t7vtvcA/NHi8V/7Sr69/sBvHZrXjWC+7Afg8sHVXHv/++IX92f7czf7coL4ngDfVl4+mk4DMWnpsffsHu3rfnRwf++4g67v9eQ5yI6OpTue39x9Ufw39KiL+EhEnd2Gsh7qx/gGqvxy37lKVnXtlPV459lCqJtCmfFXz36n+Impv67qm9mON7kYtS3PVC1+erb8/Wqx/tu2+I+J1EfHTiHgkIp4C/jdrHo/2x3RJZj7XwX2/Cnhl/SKfJyPiSeAUVh2HV7LmY7A27bd/ZXFf3yju53EgaHysXgk8nplPtxur3La8n1dRPQ6Li/G/TXWWYw3rcBzL/QBYlHVH6GA9AJl5K9VzZ2pEvIGqwc5uVFOto+dc+8dhbT833R2/y49NZv4Z+BvVWf23UZ0deTgiXg/8A9UZibYxP9HuubUtDY4TXdu/rvw8tjkWeB1wT0T8LiLe0cm2A439eRX7cyEiToqIeRGxrB5rswb1daQrPbbT54t9d43x7bv07xfprSYidqF68G5qvy4zn87MT2Tmq4F3ASdGxN5tqzsYsqPlbbYtLo+lOovwV6onysuKuoZQzS/q6rgPUz2RyrGXs3rj64q/1jW1H2tRN8fpqrOp5l6Nz8yXUzXLaLdN+33v7Fg8BNyfmZsXXyMy84B6/WLWfAzWpv32Dxf39aF297VJZt7cYIyHgS0jYkS7scrjWu7XQ1RnN7Yuxn55Zu7QQY1dOY4d7QfA6IiITtaXZlH9a/NI4NJOfhl2ZjHVv8Aa1dbI2p7/7XXnsYGqGR9CNX9wUX19GtW/6uYWY36x3Zgvy8wfNBivu/tXWmNfM/O+zDyC6pf3l4FLI2LTbozZL9mf12B/rkXE26jmmB4GbJGZm1NNUWhfX0d1daXHru1xte+uzr7LAAjIEfHy+q+BH1L92+VPDbZ5R0S8tn4CL6P6F/mKevWjVPOnuuufImL7iHgZ1b8IL63/mv9vYFhEHBgRG1LN3dm4uN2jwLjyLY/a+QHw8YjYLiKGU/0le3FmLu9OcXUtlwBfjIgREfEq4ESquUrNMIJqAvwz9V/Gx6/neLcBT0fEpyNik4gYEhFvrH/RQrVvn4mILSJiDNX8srX5ZL39tlRzvy6ul59Tj7UDQERsFhGHNhogMx+ienHBlyJiWETsSPXXacPjmpmLgV8BX6ufqxtExGsi4h86qLErx7Gj/YCqAXwsIjas92ECcGUH9/U94N1UzfrCDrZZm0uAYyJiQv2zsLb3tu3uz1uXH5var4GPUP1LDap/LX6E6l+0bWfb/hM4LiJ2jcqm9c/riDWHW+15Nroeq6seBbaKiM3aFkTEP0XEyMxcQfVCFVjViwYc+3Nj9uc1altONQd2aER8Dnh5J9uv9hitQ4/tqAb77ir2Xfp3QP5JRDxN9VfJvwJfB47pYNvxwNVUr2j8LfCtzLyuXvcl4LNRnfI/qRv3/12qeXSPUL1Y62MAmbmMan7Sd6jOBvyN6tXAbX5Uf18aEb9vMO759dg3UL0y+Dm6Fv4a+Wh9/3+hOnNzUT1+M5xENUH/aaofhIs737xz9Q/VO6j+bXM/1RmX71D96w2q+UQP1Ot+RXXM1uYKqnl5c4GfUb/tVGZeTvVX5Q+j+vfancDbOxnnCKr5hg9TvTjk1My8upPtj6J6IcvdVHO1LgVGdbBtV45jw/2o3Ur1fP8r1auXD8nMpY3uqA77v6f6i/vGTurvUGb+HDiT6gUs86nmmUJ1RqeRbwCHRPWm/2d2YfzuPja/pvpl19aob6I6Y9h2ncycQ/Wipm9SPR7zqeYRNvIFqp/f+6l6yKWd7Fv72u+hClR/qfvLK6leXHNXRDxDdSwOz8xnOxunn7I/r539ufJL4BdUf7w8QHVMO5sy0Ogx6k6PbcS+u/r49l1WTRSX1A9ERFL9G3B+g3VHU70QY0o3xjsfeDgzP9tD9U2gaqYbd/esWn8QEcdTNdfunJ2S1I/Zd1urVX23P59BlrQeImIc8B7W8wMcIuLdEbFxRGxBddbhJwOlSUfEqKg+anyDqF508gmq/xpIUrfZd9eur/RdA7I0CEXEv1OdcfiPzLx/PYf7ENX7W/6Zav7o+s5v7Es2onpF/NPAtVT/Zv1WSyuS1C/Zd7usT/Rdp1hIkiRJBc8gS5IkSYWhrS5gfWy99dY5bty4VpchSd12++23/zUzR659y77PXiypv+qoF/frgDxu3DjmzJnT6jIkqdsioiuf/tgv2Isl9Vcd9WKnWEiSJEkFA7IkSZJUMCBLkiRJhX49B1lSz3rxxRdZuHAhzz33XKtLGTCGDRvGmDFj2HDDDVtdiqR+wD7cHN3txQZkSSstXLiQESNGMG7cOCKi1eX0e5nJ0qVLWbhwIdttt12ry5HUD9iHe9669GKnWEha6bnnnmOrrbayKfeQiGCrrbbyTJCkLrMP97x16cUGZEmrsSn3LI+npO6yb/S87h5TA7IkSZJUcA6ypA6NO/lnPTregtMPXPs2Cxbwjne8gzvvvHPlspkzZzJ8+HBOOumkhreZO3cuDz/8MAcccECP1Lmu402dOpWvfvWrTJo0qUfqkCT7cGv6sGeQJfV7c+fO5corr2y4bvny5T06niRpTQOtDxuQJfUbU6dO5dOf/jSTJ0/mda97HTfeeCMvvPACn/vc57j44ouZOHEiF198MTNnzuTII4/krW99K0ceeSRLlizhve99L7vssgu77LILv/nNbwC47bbbeMtb3sJOO+3E7rvvzr333ttwvL/97W984AMfYPLkyey0005cccUVADz77LMcfvjhTJgwgXe/+908++yzrTw8ktR0g6UPO8VCUr+yfPlybrvtNq688ko+//nPc/XVV/OFL3yBOXPm8M1vfhOo/hV49913c9NNN7HJJpvw/ve/n49//ONMmTKFBx98kP3224958+bxhje8gRtvvJGhQ4dy9dVXc8opp3DZZZetMd4pp5zCXnvtxfnnn8+TTz7J5MmT2Wefffj2t7/Ny172MubNm8cdd9zBm9/85lYeGknqFYOhDxuQJfUpHb3SuG35e97zHgB23nlnFixY0OE473rXu9hkk00AuPrqq7n77rtXrnvqqad45plnWLZsGdOmTeO+++4jInjxxRcbjvWrX/2K2bNn89WvfhWo3obpwQcf5IYbbuBjH/sYADvuuCM77rhj93ZWkvog+7ABWVIfs9VWW/HEE0+stuzxxx9f+ebuG2+8MQBDhgzpdF7bpptuuvLyihUruOWWWxg2bNhq23zkIx9hzz335PLLL2fBggVMnTq14ViZyWWXXcbrX//6ddklSepX7MPOQZbUxwwfPpxRo0Zx7bXXAlVT/sUvfsGUKVM6vM2IESN4+umnO1y/7777ctZZZ628PnfuXACWLVvG6NGjAbjgggs6HG+//fbjrLPOIjMB+MMf/gDAHnvswUUXXQTAnXfeyR133NGNPZWkvsk+7BlkSZ3oytsBNcOFF17IjBkzOPHEEwE49dRTec1rXtPh9nvuuSenn346EydO5DOf+cwa688880xmzJjBjjvuyPLly9ljjz0455xz+NSnPsW0adM47bTTOPDAAzsc79/+7d844YQT2HHHHVmxYgXbbbcdP/3pTzn++OM55phjmDBhAhMmTGDnnXfu+YMhaVCzD7emD0dbEu+PJk2alHPmzGl1GdKAMW/ePCZMmNDqMgacRsc1Im7PzAHxhsn2Yqnn2Iebpzu92CkWkiRJUsGALEmSJBUG5Rzknv7Yxr6oVXOWJKmrBnovtg9L/ZdnkCVJkqSCAVmSJEkqDMopFpIkNd3MzVpdQfPNXNbqCqSmMCBL6lhP/4Jfyy/TpUuXsvfeewPwyCOPMGTIEEaOHAnAbbfdxkYbbdTtu7z++uvZaKON2H333bt1u3HjxjFnzhy23nrrbt+nJPUY+3BL+rABWVKfsdVWW638dKWZM2cyfPhwTjrppJXrly9fztCh3Wtb119/PcOHD+92Y5akwcg+XHEOsqQ+7eijj+a4445j11135VOf+hR//vOf2X///dl5551529vexj333APAT37yE3bddVd22mkn9tlnHx599FEWLFjAOeecwxlnnMHEiRO58cYbWbJkCe9973vZZZdd2GWXXfjNb34DVGdN9t13X3bYYQc++MEP0p8/REmSetJg7MOeQZbU5y1cuJCbb76ZIUOGsPfee3POOecwfvx4br31Vj784Q9z7bXXMmXKFG655RYigu985zt85Stf4Wtf+xrHHXfcamdA3v/+9/Pxj3+cKVOm8OCDD7Lffvsxb948Pv/5zzNlyhQ+97nP8bOf/YzzzjuvxXstSX3HYOvDBmRJfd6hhx7KkCFDeOaZZ7j55ps59NBDV657/vnngap5v+9972Px4sW88MILbLfddg3Huvrqq7n77rtXXn/qqad45plnuOGGG/jxj38MwIEHHsgWW2zRxD2SpP5lsPVhA7KkPm/TTTcFYMWKFWy++eYr58eVPvrRj3LiiSfyrne9i+uvv56ZM2c2HGvFihXccsstDBs2rIkVS9LAMtj6sHOQJfUbL3/5y9luu+340Y9+BEBm8sc//hGAZcuWMXr0aABmzZq18jYjRozg6aefXnl933335ayzzlp5va3J77HHHlx00UUA/PznP+eJJ55o6r5IUn80WPqwZ5AldawPvsfp97//fY4//nhOO+00XnzxRQ4//HDe9KY3MXPmTA499FC22GIL9tprL+6//34A3vnOd3LIIYdwxRVXcNZZZ3HmmWcyY8YMdtxxR5YvX84ee+zBOeecw6mnnsoRRxzBDjvswO67787YsWNbvKeShH24RaI/v1J70qRJOWfOnG7fbtzJP2tCNX3LgtMPbHUJ6ofmzZvHhAkTWl3GgNPouEbE7Zk5qUUl9Sh7cWMLhr2/1SU0Xx8Mb/2dfbh5utOLnWIhSZIkFQzIkiRJUsGALGk1/XnaVV/k8ZTUXfaNntfdY2pAlrTSsGHDWLp0qc25h2QmS5cu7dNvZSSpb7EP97x16cW+i4WklcaMGcPChQtZsmRJq0sZMIYNG8aYMWNaXYakfsI+3Bzd7cUGZEkrbbjhhh1+8pH6p4gYBtwAbEzV8y/NzFMj4gLgH4C2tyE4OjPnRkQA3wAOAP5eL/9971cuDU724b7BgCxJA9vzwF6Z+UxEbAjcFBE/r9d9MjMvbbf924Hx9deuwNn1d0kaNJyDLEkDWFaeqa9uWH91NrnxIODC+na3AJtHxKhm1ylJfYkBWZIGuIgYEhFzgceAqzLz1nrVFyPijog4IyI2rpeNBh4qbr6wXtZ+zOkRMSci5jhXUtJAY0CWpAEuM1/KzInAGGByRLwR+AzwBmAXYEvg090c89zMnJSZk0aOHNnTJUtSSxmQJWmQyMwngeuA/TNzcT2N4nngv4DJ9WaLgG2Lm42pl0nSoGFAlqQBLCJGRsTm9eVNgH8E7mmbV1y/a8XBwJ31TWYDR0VlN2BZZi7u9cIlqYV8FwtJGthGAbMiYgjVSZFLMvOnEXFtRIwEApgLHFdvfyXVW7zNp3qbt2N6v2RJai0DsiQNYJl5B7BTg+V7dbB9AjOaXZck9WVOsZAkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSp0LSAHBHbRsR1EXF3RNwVEf9SL98yIq6KiPvq71vUyyMizoyI+RFxR0S8uVm1SZIkSR1p5hnk5cAnMnN7YDdgRkRsD5wMXJOZ44Fr6usAbwfG11/TgbObWJskSZLUUNMCcmYuzszf15efBuYBo4GDgFn1ZrOAg+vLBwEXZuUWYPOIGNWs+iRJkqRGemUOckSMA3YCbgW2yczF9apHgG3qy6OBh4qbLayXtR9rekTMiYg5S5YsaV7RkiRJGpSaHpAjYjhwGXBCZj5VrsvMBLI742XmuZk5KTMnjRw5sgcrlSRJkpockCNiQ6pw/P3M/HG9+NG2qRP198fq5YuAbYubj6mXSZIkSb2mme9iEcB5wLzM/HqxajYwrb48DbiiWH5U/W4WuwHLiqkYkiRJUq8Y2sSx3wocCfwpIubWy04BTgcuiYhjgQeAw+p1VwIHAPOBvwPHNLE2SZIkqaGmBeTMvAmIDlbv3WD7BGY0qx5JkiSpK/wkPUmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEkawCJiWETcFhF/jIi7IuLz9fLtIuLWiJgfERdHxEb18o3r6/Pr9eNaugOS1AIGZEka2J4H9srMNwETgf0jYjfgy8AZmfla4Ang2Hr7Y4En6uVn1NtJ0qBiQJakASwrz9RXN6y/EtgLuLRePgs4uL58UH2dev3eERG9U60k9Q0GZEka4CJiSETMBR4DrgL+DDyZmcvrTRYCo+vLo4GHAOr1y4CtGow5PSLmRMScJUuWNHkPJKl3GZAlaYDLzJcycyIwBpgMvKEHxjw3Mydl5qSRI0eu73CS1KcYkCVpkMjMJ4HrgLcAm0fE0HrVGGBRfXkRsC1AvX4zYGnvVipJrWVAlqQBLCJGRsTm9eVNgH8E5lEF5UPqzaYBV9SXZ9fXqddfm5nZawVLUh8wdO2bSJL6sVHArIgYQnVS5JLM/GlE3A38MCJOA/4AnFdvfx7w3YiYDzwOHN6KoiWplQzIkjSAZeYdwE4Nlv+Faj5y++XPAYf2QmmS1Gc5xUKSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpELTAnJEnB8Rj0XEncWymRGxKCLm1l8HFOs+ExHzI+LeiNivWXVJkiRJnWnmGeQLgP0bLD8jMyfWX1cCRMT2wOHADvVtvhURQ5pYmyQNChGxbURcFxF3R8RdEfEv9XJPWEhSB4Y2a+DMvCEixnVx84OAH2bm88D9ETEfmAz8tln1SdIgsRz4RGb+PiJGALdHxFX1ujMy86vlxu1OWLwSuDoiXpeZL/Vq1ZLUQk0LyJ34SEQcBcyhatpPAKOBW4ptFtbL1hAR04HpAGPHjm1yqf3YzM1aXUHzzVzW6gqkPi8zFwOL68tPR8Q8OuivNU9YSBr0evtFemcDrwEmUjXsr3V3gMw8NzMnZeakkSNH9nB5kjRw1f/V2wm4tV70kYi4o37NyBb1stHAQ8XNOjxhIUkDVa8G5Mx8NDNfyswVwH9SnZUAWARsW2w6pl4mSeoBETEcuAw4ITOfYj1PWETE9IiYExFzlixZ0tPlSlJL9WpAjohRxdV3A23vcDEbODwiNo6I7YDxwG29WZskDVQRsSFVOP5+Zv4Y1v+Ehf/NkzSQNW0OckT8AJgKbB0RC4FTgakRMRFIYAHwIYDMvCsiLgHupnpByQxfECJJ6y8iAjgPmJeZXy+Wj6rnJ8OaJywuioivU71IzxMWkgadZr6LxRENFp/XyfZfBL7YrHokaZB6K3Ak8KeImFsvOwU4whMWktRYK97FQpLUSzLzJiAarLqyk9t4wkLSoOZHTUuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVOhSQI6It3ZlmSSpOezDktR7unoG+awuLpMkNYd9WJJ6ydDOVkbEW4DdgZERcWKx6uXAkGYWJkla/z4cEdsCFwLbAAmcm5nfiIgtgYuBccAC4LDMfCIiAvgGcADwd+DozPx9z+2RJPV9azuDvBEwnCpIjyi+ngIOaW5pkiTWvw8vBz6RmdsDuwEzImJ74GTgmswcD1xTXwd4OzC+/poOnN1zuyJJ/UOnZ5Az89fAryPigsx8oJdqkiTV1rcPZ+ZiYHF9+emImAeMBg4CptabzQKuBz5dL78wMxO4JSI2j4hR9TiSNCh0GpALG0fEuVT/ilt5m8zcqxlFSZLWsN59OCLGATsBtwLbFKH3EaopGFCF54eKmy2sl60WkCNiOtUZZsaOHduN3ZCkvq+rAflHwDnAd4CXmleOJKkD69WHI2I4cBlwQmY+VU01rmRmRkR2Z7zMPBc4F2DSpEnduq0k9XVdDcjLM9N5aJLUOuvchyNiQ6pw/P3M/HG9+NG2qRMRMQp4rF6+CNi2uPmYepkkDRpdfZu3n0TEhyNiVERs2fbV1MokSaV16sP1u1KcB8zLzK8Xq2YD0+rL04AriuVHRWU3YJnzjyUNNl09g9zWRD9ZLEvg1T1bjiSpA+vah98KHAn8KSLm1stOAU4HLomIY4EHgMPqdVdSvcXbfKq3eTtmvSuXpH6mSwE5M7drdiGSpI6tax/OzJuA6GD13g22T2DGutyXJA0UXQrIEXFUo+WZeWHPliNJasQ+LEm9p6tTLHYpLg+jOuvwe6pPZ5IkNZ99WJJ6SVenWHy0vB4RmwM/bEZBkqQ12Yclqfd09V0s2vsb4LxkSWod+7AkNUlX5yD/hOrV0gBDgAnAJc0qSpK0OvuwJPWers5B/mpxeTnwQGYubEI9kqTG7MOS1Eu6NMUiM38N3AOMALYAXmhmUZKk1dmHJan3dCkgR8RhwG3AoVRvJn9rRBzSzMIkSavYhyWp93R1isW/Artk5mMAETESuBq4tFmFSZJWYx+WpF7S1Xex2KCtKdeWduO2kqT1Zx+WpF7S1TPIv4iIXwI/qK+/D7iyOSVJkhqwD0tSL+k0IEfEa4FtMvOTEfEeYEq96rfA95tdnCQNdvZhSep9azuD/H+AzwBk5o+BHwNExP+s172zibVJkuzDktTr1jZ/bZvM/FP7hfWycU2pSJJUsg9LUi9bW0DevJN1m/RgHZKkxjbvZJ19WJKaYG0BeU5E/HP7hRHxQeD25pQkSSrYhyWpl61tDvIJwOUR8b9Y1YgnARsB725iXZKkygnYhyWpV3UakDPzUWD3iNgTeGO9+GeZeW3TK5Mk2YclqQW69D7ImXkdcF2Ta5EkdcA+LEm9x09hkiRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKnQtIAcEedHxGMRcWexbMuIuCoi7qu/b1Evj4g4MyLmR8QdEfHmZtUlSZIkdaaZZ5AvAPZvt+xk4JrMHA9cU18HeDswvv6aDpzdxLokSZKkDjUtIGfmDcDj7RYfBMyqL88CDi6WX5iVW4DNI2JUs2qTJEmSOtLbc5C3yczF9eVHgG3qy6OBh4rtFtbL1hAR0yNiTkTMWbJkSfMqlSRJ0qDUshfpZWYCuQ63OzczJ2XmpJEjRzahMkmSJA1mvR2QH22bOlF/f6xevgjYtthuTL1MkiRJ6lW9HZBnA9Pqy9OAK4rlR9XvZrEbsKyYiiFJkiT1mqHNGjgifgBMBbaOiIXAqcDpwCURcSzwAHBYvfmVwAHAfODvwDHNqkuSJEnqTNMCcmYe0cGqvRtsm8CMZtUiSZIkdZWfpCdJkiQVDMiSJElSwYAsSZIkFQzIkiRJUsGALEmSJBUMyJIkSVLBgCxJkiQVDMiSJElSwYAsSZIkFQzIkiRJUsGALEmSJBUMyJI0gEXE+RHxWETcWSybGRGLImJu/XVAse4zETE/Iu6NiP1aU7UktZYBWZIGtguA/RssPyMzJ9ZfVwJExPbA4cAO9W2+FRFDeq1SSeojDMiSNIBl5g3A413c/CDgh5n5fGbeD8wHJjetOEnqowzIkjQ4fSQi7qinYGxRLxsNPFRss7BetoaImB4RcyJizpIlS5pdqyT1KgOyJA0+ZwOvASYCi4GvdXeAzDw3Mydl5qSRI0f2cHmS1FoGZEkaZDLz0cx8KTNXAP/JqmkUi4Bti03H1MskaVAxIEvSIBMRo4qr7wba3uFiNnB4RGwcEdsB44Hbers+SWq1oa0uQJLUPBHxA2AqsHVELAROBaZGxEQggQXAhwAy866IuAS4G1gOzMjMl1pQtiS1lAFZkgawzDyiweLzOtn+i8AXm1eRJPV9TrGQJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJGsAi4vyIeCwi7iyWbRkRV0XEffX3LerlERFnRsT8iLgjIt7cusolqXUMyJI0sF0A7N9u2cnANZk5Hrimvg7wdmB8/TUdOLuXapSkPsWALEkDWGbeADzebvFBwKz68izg4GL5hVm5Bdg8Ikb1SqGS1IcYkCVp8NkmMxfXlx8BtqkvjwYeKrZbWC9bQ0RMj4g5ETFnyZIlzatUklrAgCxJg1hmJpDrcLtzM3NSZk4aOXJkEyqTpNYxIEvS4PNo29SJ+vtj9fJFwLbFdmPqZZI0qBiQJWnwmQ1Mqy9PA64olh9Vv5vFbsCyYiqGJA0aQ1tdgCSpeSLiB8BUYOuIWAicCpwOXBIRxwIPAIfVm18JHADMB/4OHNPrBUtSH2BAlqQBLDOP6GDV3g22TWBGcyuSpL7PKRaSJElSwYAsSZIkFQzIkiRJUsGALEmSJBUMyJIkSVLBgCxJkiQVDMiSJElSwYAsSZIkFQzIkiRJUqEln6QXEQuAp4GXgOWZOSkitgQuBsYBC4DDMvOJVtQnSZKkwauVZ5D3zMyJmTmpvn4ycE1mjgeuqa9LkiRJvaovTbE4CJhVX54FHNy6UiRJkjRYtSogJ/CriLg9IqbXy7bJzMX15UeAbVpTmiRJkgazlsxBBqZk5qKIeAVwVUTcU67MzIyIbHTDOlBPBxg7dmzzK5UkSdKg0pIzyJm5qP7+GHA5MBl4NCJGAdTfH+vgtudm5qTMnDRy5MjeKlmSJEmDRK8H5IjYNCJGtF0G9gXuBGYD0+rNpgFX9HZtkiRJUiumWGwDXB4Rbfd/UWb+IiJ+B1wSEccCDwCHtaA2SZIkDXK9HpAz8y/AmxosXwrs3dv1SJIkSaW+9DZvkiRJUssZkCVJkqRCq97mTZIkSb1h5matrqD5Zi7r0eE8gyxJkiQVPIMsSZIGrXEn/6zVJTTdgmGtrqD/8QyyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQU/KERqkUHx5vSnH9jqEiRJ6jbPIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBT9JT1LzzNys1RU038xlra5AktTDPIMsSZIkFQzIkiRJUsGALEmSJBUMyJIkSVLBgCxJkiQVfBcLSRqkImIB8DTwErA8MydFxJbAxcA4YAFwWGY+0aoaJakVPIMsSYPbnpk5MTMn1ddPBq7JzPHANfV1SRpUDMiSpNJBwKz68izg4NaVIkmtYUCWpMErgV9FxO0RMb1etk1mLq4vPwJs05rSJKl1nIMsSYPXlMxcFBGvAK6KiHvKlZmZEZGNblgH6ukAY8eObX6lktSLPIMsSYNUZi6qvz8GXA5MBh6NiFEA9ffHOrjtuZk5KTMnjRw5srdKlqReYUCWpEEoIjaNiBFtl4F9gTuB2cC0erNpwBWtqVCSWscpFpI0OG0DXB4RUP0uuCgzfxERvwMuiYhjgQeAw1pYoyS1hAFZkgahzPwL8KYGy5cCe/d+RZLUdzjFQpIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpEKfC8gRsX9E3BsR8yPi5FbXI0mDjX1Y0mDXpwJyRAwB/i/wdmB74IiI2L61VUnS4GEflqQ+FpCBycD8zPxLZr4A/BA4qMU1SdJgYh+WNOgNbXUB7YwGHiquLwR2LTeIiOnA9PrqMxFxby/V1q8EbA38tdV1NNXno9UVaC18HnbqVT1ZRg9aax8Ge3FX+PxXX+FzsVMNe3FfC8hrlZnnAue2uo6+LiLmZOakVtehwc3n4cBlL147n//qK3wudl9fm2KxCNi2uD6mXiZJ6h32YUmDXl8LyL8DxkfEdhGxEXA4MLvFNUnSYGIfljTo9akpFpm5PCI+AvwSGAKcn5l3tbis/sp/faov8HnYz9iHe5TPf/UVPhe7KTKz1TVIkiRJfUZfm2IhSZIktZQBWZIkSSoYkPu5tX0kbERsHBEX1+tvjYhxLShTA1xEnB8Rj0XEnR2sj4g4s34e3hERb+7tGqVmsQ+rL7AP9ywDcj/WxY+EPRZ4IjNfC5wBfLl3q9QgcQGwfyfr3w6Mr7+mA2f3Qk1S09mH1YdcgH24xxiQ+7eufCTsQcCs+vKlwN4R4UcfqUdl5g3A451schBwYVZuATaPiFG9U53UVPZh9Qn24Z5lQO7fGn0k7OiOtsnM5cAyYKteqU5apSvPVak/sg+rv7APd4MBWZIkSSoYkPu3rnwk7MptImIosBmwtFeqk1bx44s1UNmH1V/Yh7vBgNy/deUjYWcD0+rLhwDXpp8Oo943GziqfhX1bsCyzFzc6qKkHmAfVn9hH+6GPvVR0+qejj4SNiK+AMzJzNnAecB3I2I+1eT9w1tXsQaqiPgBMBXYOiIWAqcCGwJk5jnAlcABwHzg78AxralU6ln2YfUV9uGe5UdNS5IkSQWnWEiSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDstRFEfHKiLi0m7e5ICIOaVZNkjTY2IvVGwzIUgP1p12tdj0zH85MG6wk9RJ7sVrFgKwBJSLGRcQ99dmC/46I70fEPhHxm4i4LyIm11+/jYg/RMTNEfH6+rZHR8TsiLgWuKbB9XERcWe97ZCI+I+I+F1E3BERH6qXR0R8MyLujYirgVe07GBIUovYi9Xf+Ul6GoheCxwKfIDqY2DfD0wB3gWcAhwFvK3+BKx9gP8NvLe+7ZuBHTPz8Yg4ut31ccV9HEv1MZ27RMTGwG8i4lfATsDrge2BbYC7gfObubOS1EfZi9VvGZA1EN2fmX8CiIi7gGsyMyPiT8A4YDNgVkSMB5L6ozhrV2Xm451cb7MvsGMxp20zYDywB/CDzHwJeLg+4yFJg5G9WP2WAVkD0fPF5RXF9RVUz/l/B67LzHfXZyKuL7b/W7ux2l9vE8BHM/OXqy2MOGAda5akgcZerH7LOcgajDYDFtWXj17HMX4JHB8RGwJExOsiYlPgBuB99by4UcCe61usJA1Q9mL1WQZkDUZfAb4UEX9g3f+L8h2qOW2/r18s8u16rMuB++p1FwK/Xf9yJWlAsherz4rMbHUNkiRJUp/hGWRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSp8P8BjfbB6AKQoTAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "estimate.interpret(method_name=\"confounder_distribution_interpreter\",var_type='discrete',\n", " var_name='married', fig_size = (10, 7), font_size = 12)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sanity check: compare to manual IPW estimate" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Causal Estimate is 1639.8075601428245\n" ] } ], "source": [ "df = model._data\n", "ps = df['ps']\n", "y = df['re78']\n", "z = df['treat']\n", "\n", "ey1 = z*y/ps / sum(z/ps)\n", "ey0 = (1-z)*y/(1-ps) / sum((1-z)/(1-ps))\n", "ate = ey1.sum()-ey0.sum()\n", "print(\"Causal Estimate is \" + str(ate))\n", "\n", "# correct -> Causal Estimate is 1634.9868359746906" ] } ], "metadata": { "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.8.5" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": false, "sideBar": true, "skip_h1_title": true, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }