diff --git a/docs/source/glossary.rst b/docs/source/glossary.rst index c6aedbb2..6f31530d 100644 --- a/docs/source/glossary.rst +++ b/docs/source/glossary.rst @@ -53,6 +53,9 @@ Glossary One-group posttest-only design A design where a single group is exposed to a treatment and assessed on an outcome measure. There is no pretest measure or comparison group. + Parallel trends assumption + An assumption made in difference in differences designs that the trends (over time) in the outcome variable would have been the same between the treatment and control groups in the absence of the treatment. + Panel data Time series data collected on multiple units where the same units are observed at each time point. diff --git a/docs/source/index.rst b/docs/source/index.rst index 4e3cc9c8..3e312037 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -141,6 +141,7 @@ Documentation outline :caption: Knowledge Base design_notation.md + quasi_dags.ipynb glossary.rst .. toctree:: diff --git a/docs/source/quasi_dags.ipynb b/docs/source/quasi_dags.ipynb new file mode 100644 index 00000000..dff18e53 --- /dev/null +++ b/docs/source/quasi_dags.ipynb @@ -0,0 +1,461 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Causal DAGS for Quasi-Experiments" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This page provides an overview of causal Directed Acyclic Graphs (DAG's) for some of the most common quasi-experiments. It takes inspiration from a paper by {cite:t}`steiner2017graphical`, and the books by {cite:t}`cunningham2021causal` and {cite:t}`huntington2021effect`, and readers are encouraged to consult these sources for more details." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "tags": [ + "remove-input" + ] + }, + "outputs": [], + "source": [ + "import daft\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "tags": [ + "remove-input" + ] + }, + "outputs": [], + "source": [ + "ff = \"times new roman\"\n", + "plt.rcParams[\"font.family\"] = ff\n", + "\n", + "GRID_UNIT = 2.0\n", + "DPI = 200\n", + "NODE_EC = \"none\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Before we take a look at randomized controlled trials (RCTs) and quasi-experiments, let's first consider the concept of confounding. Confounding occurs when a variable (or variables) causally influence both the treatment and the outcome and is very common in observational studies. This can lead to biased estimates of the treatment effect (the causal effect of $Z \\rightarrow Y$). The following causal DAG illustrates the concept of confounding. Note that the confounder is written as a vector because there may be multiple confounding variables, $\\mathbf{X}=x_1, x_2,x_3$." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "tags": [ + "remove-input" + ] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAATMAAAEMCAYAAACodFEmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAB7CAAAewgFu0HU+AAAQAElEQVR4nO3de2jV9R/H8ddxaqV4g5kXnJVGOZfzsnlJXIlZm5akKaywpHQ6/8myzMgmJBGZVkoUkbJRIKWpZavpzMVW2siQ8oLTLlNzapaZyJzT3b6/P/p58Ox6pmfne877+3zAwHPOd/Ku1tPP+5zt6HMcxxEARLl2bg8AAKFAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMAJhAzACYQMwAmEDMcM0uXbqkQYMGyefzBXw8/fTTzX5ecXGxYmJiGnxeQUFBmCaHRT7HcRy3h0D0Ki4uVkpKiurq6vz3+Xw+ffvtt0pJSWlw/eXLlzVs2DAdPnw44P6MjAytXbu2zeeFXZzMcF3Gjh2rBQsWBNznOI7mzJmjysrKBtcvW7asQcj69eunN998s03nhH2czHDdLl68qMTERJWWlgbcv2jRIq1cudJ/+6efftLo0aNVU1MTcF1eXp4mT54clllhFzFDSBQVFWnChAm6+sspJiZGxcXFGjVqlGpqapScnKx9+/YFfN6sWbP00UcfhXtcGMSaiZAYP3685s+fH3BfbW2tZs+eraqqKi1fvrxByHr37q3Vq1eHcUpYxskMIXPhwgXddddd+uOPPwLunzlzpjZu3KiqqqqA+z/77DNNmzYtnCPCMGKGkNqxY4ceeOCBFq9LT0/X+vXrwzARvIKYIeQyMjKUnZ3d5OM9e/bUwYMH1bNnzzBOBeuIGULu/PnzSkhI0MmTJxt9fP369UpPTw/zVLCOFwAQct26ddOECRMafaxz585BraFAaxEzhNx3332ndevWNfpYRUWFFi5cGOaJ4AWsmQipyspKJSYm6vfff2/2um3btiktLS1MU8ELOJkhpLKyshqErH379g2umzdvnsrLy8M1FjyAmCFkdu/e3eCbYH0+n7788kvFx8cH3F9WVqYXXnghjNPBOmKGkKiqqtLs2bMD3j1DkjIzM5WWlqacnBy1axf45bZmzRoVFhaGc0wYRswQEsuWLVNJSUnAfXFxcVqxYoUkacyYMQ2e+HccRxkZGbp48WLY5oRdvACA6/bzzz/7f5j8alu3btWkSZP8tysrKzV06FD99ttvAdc988wz/Iwmrhsxw3WpqanRyJEjtXfv3oD7m3o3jF27dumee+4JeHeNdu3aaefOnRo7dmxbjwvDWDNxXV5//fUGIevVq5dWrVrV6PXjxo1r8LbadXV1mjNnji5dutRWY8IDOJnhmpWUlGj48OEN3g1j06ZNmj59epOfV1FRocTERB05ciTg/hdffFHLly9vk1lhHzEDYAJrJgATiBkAE4gZABOIGQATiBkAE4gZABOIGQATiBkAE4gZABOIGQATiBkAE4gZABOIGQATiBkAE4gZABOIGQATiBkAE4gZABOIGQATiBkAE4gZABOIGQATiBkAE4gZJEmnT59WZWWl22O02tGjR8Vf/QqJmOH/+vTpo06dOunjjz92e5SgOI6jadOmacCAAZo6darb4yACEDNo3rx5/l/PnDlT//77r4vTBOett97Sli1bJEm5ubnav3+/uwPBdT6HM7qnlZaW6vbbb/ffvvfee1VUVOTeQK3g8/kCbtfV1TW4D97Byczjrg6ZJBUWFro0SeuVlJQE3Gbd9DZi5mFXr5eStHv37qg62cTHxys9Pd1/m3XT21gzPSqa18v6WDchcTLzrGheL+tj3YREzDwp2tfL+lg3IbFmeo6l9bI+1k1v42TmMZbWy/pYN72NmHmItfWyPtZNb2PN9AjL62V9rJvexMnMIyyvl/WxbnoTMfMA6+tlfayb3sSaaZyX1sv6WDe9hZOZcV5aL+tj3fQWYmaY19bL+lg3vYU10ygvr5f1sW56Ayczo7y8XtbHuukNxMwgr6+X9bFuegNrpjGsl01j3bSNk5kxrJdNY920jZgZwnrZPNZN21gzjWC9DB7rpk2czIxgvQwe66ZNxMwA1svWYd20iTUzyrFeXjvWTVs4mUU51strx7ppCzGLYqyX14d10xbWzCjFehk6rJs2cDKLUqyXocO6aQMxi0Ksl6HFumkDa2YEO3PmjMrKyjRixAj/fayXbae5ddNxHH399ddKTU11YzQEgZNZBMvOztaoUaO0dOlSVVVVSWK9bEtNrZunTp3SlClTlJaWxoktkjmISDU1Nc6tt97qSHIkOUOGDHFSUlL8tyU5u3fvdntMc9LT0wP+HWdlZTndu3f3354/f77bI6IJrJkRKi8vTw899FCTj7Netp3mnn/s3LmzTp06pa5du4ZxIgSDNTNCvf/++80+/vbbb4dpEm9xHEcvv/xyk49XVFRo3bp1YZwIweJkFoGOHj2qgQMHqrn/NDExMVqyZImysrLUsWPHME5n16lTp5SZmamvvvqq2esSEhJ04MABXkGOMJzMItCaNWuaDZkk1dbW6tVXX1VycrJOnz4dpsnsKigoUEJCQoshk6SDBw9q165dYZgKrUHMIszly5eVnZ0d9PWpqanq1atXG07kDSNHjtSgQYOCvr6lpwEQfsQswmzevFlnzpwJ6tpFixZpxYoVrDsh0K1bN+Xn52vMmDFBXb9p0yb99ddfbTwVWoOYRZhg/8QnZKHXmqBVV1crJycnDFMhWLwAEEEOHDigxMTEFq8jZG3r/PnzSktL0w8//NDsdf3799eRI0cUExMTpsnQHE5mESSYUxkha3vBntCOHz+ubdu2hWkqtISTWYQoLy9X3759deHChSavIWThFcwJbdKkSdq6dWsYp0JTOJlFiHXr1hGyCBPMCS0/P19Hjx4N41RoCjGLAI7jNLtiEjL3tBQ0x3H0wQcfhHkqNIY1MwLs2rVLKSkpjT5GyCJDcytnbGysTpw4oRtuuMGFyXAFJ7MI0NSpjJBFjuZOaP/88482bdrkwlS4Giczl/3999/q16+fqqurA+4nZJGpqRPa2LFj9f3337s0FSROZq7LyckhZFGkqRNacXExb9zoMmLmotra2gZPHhOyyNdU0Ph5TXcRMxf9+uuvOnnypP82IYsejQWtsLCwxXc7QdvhOTOX7d+/X0899ZQmTJhAyKLQlefQ7rzzTq1atUo9evRweyTPImYRoKamRjExMYQsStXU1Kh9+/Zuj+F5xAyACTxnBsAEYgbABGIGwARiBsAEYgbABGIGwARiBsAEYgbABGIGwARiBsAEYgbABGIGwARiBsAEYgbABGIGwARiBsAEYlbPlbeuvtaP1atXu/2PgChy7tw5denSRT6fT3FxcaqpqWnxc2prazVp0iT/19wnn3wShkkjHzGrZ8+ePdf1+UOGDAnRJPCCHj16aO7cuZKkEydOBPWXCT///PPKz8+XJGVlZemxxx5r0xmjBW+bXc+RI0d08eLFoK4tLy9Xenq6ysrKJEkjRozQzp071alTp7YcEcaUlZVp4MCBqq6u1ujRoxv8BcNXW7t2rebNmydJmj59ujZu3MjfHXGFg2tSWVnpjB8/3pHkSHLi4+OdM2fOuD0WotSsWbP8X0vFxcWNXlNYWOh06NDBkeQMHz7cqaioCPOUkY018xpUV1drxowZKioqkiTddtttKigoUGxsrLuDIWotXrzYf8Jq7HnX0tJSzZgxQ9XV1erdu7dyc3PZAOohZq1UV1enJ554Qnl5eZKkvn37qqCgQH379nV5MkSzhIQETZ48WZK0efNmHT9+3P/Y+fPnNWXKFJ09e1Y33nijvvjiC/Xr18+tUSMWMWulzMxMbdiwQZIUGxurHTt2aMCAAS5PBQsWL14s6b9XK999913/rx999FEdOnRIkpSTk6NRo0a5NmNEc3vPjSbPPfec/3mNrl27Onv27HF7JBgzZswYR5LTvXt358KFC86CBQv8X3NLly51e7yIxquZQVq2bJleeeUVSVKnTp20fft2jRs3zt2hPOjZZ59VbGyskpKSlJSUpJtvvtntkULq888/1yOPPCJJuu+++/TNN99I4pXLYBCzIKxevVoLFy6UJHXs2FG5ublKTU11eSpvGjRokH755Rf/7bi4OH/YLATOcRzFx8cH/DPyLT/BIWYtyMnJUUZGhhzHUUxMjDZs2KDp06e7PZZn1Y9ZY6I9cNnZ2crIyJAk9enTRz/++CNP+AehvdsDRLJPP/1Uc+fOleM48vl8ys7OJmRRoKysTGVlZdqyZYv/vmgK3MCBA/2/zszMJGTBcu3ZugiXl5fn/wZFSc4777zj9kjNeumll/yz8hHcR1xcnDN16lRnx44dbv/nC7Bq1Sr/jFu2bHF7nKhBzBpRVFTk3HTTTf4vqNdee83tkVrkdhii+SPSXiV88skn/bMdO3bM7XGiBt9nVs+ePXs0ZcoUVVZWSvrve3+WLFni8lRoS+3aRdb/Bnv37pX03w+h33LLLe4OE0V4zuwqBw8eVFpamsrLyyVJ8+fP1xtvvOHyVMHJz8/3v4xv2cqVK6/79/D5fLrjjjuUnJyspKQkPfjggyGYLDSqq6tVUlIiSRo6dKjL00QXYvZ/paWluv/++3X27FlJ0syZM/Xee++5PFXwUlNTPfHtIrm5uS2+mnm1+uFKSkrS8OHD1aVLlzac8todOnRIVVVVkqRhw4a5O0yUIWaSTp48qYkTJ+rPP/+UJD388MP68MMPI279QPOiLVyNubJiSsSstTwfs3PnzmnixIk6duyYpP++j2np0qU6fPhw0L9HXFycunXr1kYTojEWwtWYffv2+X9NzFrH8zHbvn17QLgOHz6s5OTkVv0excXFuvvuu0M9GhqRlZWl/v37mwhXY66czDp27KjBgwe7O0yU8XzMDhw4cF2fHxMTw5+gYfT444+7PUKbunIyGzx4sDp06ODyNNGFH2cCYALPcAMwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEwgZgBMIGYATCBmAEw4X95pcLi6dHvaAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pgm = daft.PGM(dpi=DPI, grid_unit=GRID_UNIT, node_ec=NODE_EC)\n", + "\n", + "pgm.add_node(\"z\", \"$Z$\", 1, 0)\n", + "pgm.add_node(\"x\", \"$\\mathbf{X}$\", 1.5, 0.75)\n", + "pgm.add_node(\"y\", \"$Y$\", 2, 0)\n", + "\n", + "pgm.add_edge(\"z\", \"y\")\n", + "pgm.add_edge(\"x\", \"y\")\n", + "pgm.add_edge(\"x\", \"z\")\n", + "\n", + "pgm.render();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "One way to tell that our estimate of the causal relationship $Z \\rightarrow Y$ may be biased is the presence of a backdoor path, $Z \\leftarrow \\mathbf{X} \\rightarrow Y$. This path type is known as a \"fork\". Because $\\mathbf{X}$ is a common cause of $Z$ and $Y$, any observed statistical relation between $Z$ and $Y$ may be due to the confounding effect of $\\mathbf{X}$. \n", + "\n", + "Backdoor paths are problematic because they introduce _statistical associations_ between variables that do not reflect the true causal relationships, potentially leading to biased causal estimates. For example, if we ran a regression of the form `y ~ z`, and observe a main effect of $Z$ on $Y$, we have no way of knowing if this represents a true causal impact of $Z$ on $Y$, or if it is due to the confounding effect of $\\mathbf{X}$. \n", + "\n", + "One approach is to \"close the backdoor path\" by conditioning on the confounding variables. Practically, this could involve including confounders $\\mathbf{X}$ as a covariate in a regression model such as: `y ~ z + x₁ + x₂ + x₃`. Without explaining why, the coefficient for the main effect of $Z$ would now be an unbiased estimate of the _causal_ effect of $Z \\rightarrow Y$.\n", + "\n", + "However, unless we are very sure that we have accurate measures of _all_ confounding variables (maybe there is an $x_4$ that we don't know about or couldn't measure), it is still possible that our estimate of the causal effect is biased." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This leads us to Randomized Controlled Trials (RCTs) which are considered the gold standard for estimating causal effects. One reason for this is that we (as experimenters) intervene in the system by assigning units to treatment by {term}`random assignment`. Because of this intervention, any causal influence of the confounders upon the treatment $\\mathbf{X} \\rightarrow Z$ is broken - treamtent is now soley determined by the randomisation process, $R \\rightarrow T$. The following causal DAG illustrates the structure of an RCT." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "tags": [ + "remove-input" + ] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdEAAAEMCAYAAACbY4xqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAB7CAAAewgFu0HU+AAASY0lEQVR4nO3daWxUZR+G8XvaQkqxNBjQglQUEkJZCrVIkRRCAKWFICoIGgQUyuIHVEAwKiTywYhgAiEaXllliYJFLGBZ0ipFSCOGsNi0rAWkRSQGCUKn0O28HwiTblOmD505s1y/ZJLpzBnzN9P24nl6ZsZhWZYlAADQZGF2DwAAQKAiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAAsqdO3fUvXt3ORyOWpfZs2c3+ri8vDyFh4fXe1xOTo6PJkcwcliWZdk9BAA0RV5engYNGqTq6mrXbQ6HQwcPHtSgQYPqHX/37l317dtXp0+frnV7enq61qxZ4/V5EbxYiQIIOAMHDtQ777xT6zbLsjRt2jSVlZXVO37x4sX1AtqpUyd98cUXXp0TwY+VKICA5HQ6lZCQoKKiolq3v//++1q2bJnr62PHjik5OVmVlZW1jsvKytLIkSN9MiuCFxEFELByc3M1dOhQ1fw1Fh4erry8PPXv31+VlZXq16+fTp48WetxkydP1saNG309LoIQ27kAAtaQIUM0a9asWrdVVVVp6tSpKi8v15IlS+oFNDY2VitWrPDhlAhmrEQBBLTbt2+rV69e+vPPP2vdPnHiRGVkZKi8vLzW7Tt27NDLL7/syxERxIgogICXnZ2tF1544YHHTZgwQVu3bvXBRAgVRBRAUEhPT9e6devc3t++fXsVFBSoffv2PpwKwY6IAggKN2/eVM+ePXXlypUG79+6dasmTJjg46kQ7DixCEBQiImJ0dChQxu8r3Xr1h5t9wJNRUQBBIVff/1VW7ZsafC+0tJSzZkzx8cTIRSwnQsg4JWVlSkhIUHnz59v9Li9e/cqNTXVR1MhFLASBRDwFi5cWC+gERER9Y6bMWOGbt265auxEAKIKICAduTIkXpvnuBwOLR7927Fx8fXur24uFjz58/34XQIdkQUQMAqLy/X1KlTa32aiyTNnDlTqampWr9+vcLCav+aW716tQ4cOODLMRHEiCiAgLV48WIVFhbWui0uLk5Lly6VJA0YMKDeCUWWZSk9PV1Op9NncyJ4cWIRgIB0/Phx15vM17Rnzx6lpaW5vi4rK1OfPn107ty5Wse9++67vIcuHhoRBRBwKisr9eyzz+rEiRO1bnf36SyHDx/W4MGDa33aS1hYmA4dOqSBAwd6e1wEMbZzAQSczz77rF5AH3/8cS1fvrzB41NSUjR79uxat1VXV2vatGm6c+eOt8ZECGAlCiCgFBYWKjExsd6ns2zfvl1jx451+7jS0lIlJCTowoULtW7/4IMPtGTJEq/MiuBHRAEAMMR2LgAAhogoAACGiCgAAIaIKAAAhogoAACGiCgAAIaIKAAAhogoAACGiCgAAIaIKAAAhogoAACGiCgAAIaIKAAAhogoAACGiCgAAIaIKAAAhogoAACGiCgAAIaIKAAAhogoAACGiCgAAIaIKAAAhogoAACGiCgANIFlWbp48aLdY8BPEFEA8NCJEycUFhamLl26yOl02j0O/IDDsizL7iEAIBA4HA7X9fDwcFVWVto4DfwBK1EA8NDq1atd16uqqrRp0yYbp4E/YCUKAE1QczUqSaWlpYqKirJpGtiNlSgANMF///1X6+s2bdrYNAn8AREFgCaIjo5mWxcubOcCgAG2dSGxEgUAI2zrQiKiAGCEbV1IbOcCwENhWze0sRIFgIfAtm5oI6IA8BDY1g1tbOcCQDNgWzc0sRIFgGbAtm5oIqIA0AzY1g1NbOcCQDNiWze0sBIFgGbEtm5oIaIA0IzY1g0tbOcCgBewrRsaWIkCgBewrRsaiCgAeAHbuqGB7VwA8CK2dYMbK1EA8CK2dYMbEQUAL2JbN7ixnQsAPsC2bnBiJQoAPsC2bnAiogDgA2zrBie2cwHAh9jWDS6sRAHAh9jWDS5EFAB8iG3d4MJ2LgDYgG3d4MBKFABswLZucCCiAGADT7Z1S0pKVFBQ4OvR0ARs5wKAjRra1m3VqpU2bNigOXPmKDU1Vdu2bbNpOjwIEQUAG926daveVm5aWpr27t0rSYqIiFBxcbFiY2PtGA8PwHYuANio7rauJFdAJamyslJr16719VjwECtRALBZSUmJ4uLi3N4fFxenCxcuKCIiwodTwROsRAHAJpZlacOGDerVq1ejxxUXFysrK8tHU6EpiCgA2KCiokIvvviipk6dqps3bz7w+FWrVvlgKjQVEQUAG7Ro0ULjx4+vd3auO/v371dRUZGXp0JTEVEAsMmkSZO0ceNGj0P6v//9z8sToak4sQgAbLZ582ZNmTJFD/p1/Oijj6qkpEStWrXy0WR4EFaiAGAzT1ek//77rzIyMnw0FTxBRAHAD3gaUk4w8i9s5wKAH/Fka/fYsWNKTEz04VRwh5UoAPgRT1akrEb9BytRAPBDja1Io6Ki9NdffykmJsaGyVATK1EA8EONrUidTme9j02DPViJAoAfc7cijY+PV0FBgcevMYV3sBIFAD/mbkV66tQpHTx40KapcB8RBQA/5y6knGBkP7ZzASBA1N3a5QO77cdKFAACRN0VaWVlpQ4dOmTzVKGNiAJAALkf0s6dOysrK0uvvvqq3SOFNLZzASAAVVZWKiIiwu4xQh4RBQDAENu5AAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqIAABgiogAAGCKiAAAYIqKNWLVqlRwOR4OX1q1bq3v37po1a5ZOnTpl96gIcEuXLnX7vebJZcWKFXb/LyCA3LhxQ9HR0XI4HIqLi1NlZeUDH1NVVaW0tDTX99x3333ng0n9HxFtxIkTJ9ze53Q6debMGX399ddKTEzUtm3bfDcYgs7Ro0cf6vG9e/dupkkQCtq2bavp06dLkkpKSrR9+/YHPmbevHnat2+fJGnhwoV6/fXXvTpjoHBYlmXZPYS/Sk5O1u+//66YmBgdPnzYdXt5ebmKioq0YsUK5eXlSZJatWqlc+fO6YknnrBrXASwCxcuyOl0enTsrVu3NGHCBBUXF0uSnnnmGR06dEhRUVHeHBFBpri4WF27dlVFRYWSk5P122+/uT12zZo1mjFjhiRp7NixysjIkMPh8NWo/s1Cg6qqqqyoqChLkpWSkuL2mAEDBliSLEnWsmXLfDwlQk1ZWZk1ZMgQ1/dcfHy89c8//9g9FgLU5MmTXd9LeXl5DR5z4MABq0WLFpYkKzEx0SotLfXxlP6N7Vw3zpw541oZJCQkNHhMWFiY3n77bdfXBQUFPpkNoamiokLjxo1Tbm6uJOnpp59WTk6O2rVrZ+9gCFgLFixwrSgb+rt6UVGRxo0bp4qKCsXGxmrXrl3seNRBRN2o+ffQxv7e1LlzZ9d1T/44D5iorq7WpEmTlJWVJUnq2LGjcnJy1LFjR5snQyDr2bOnRo4cKUn64YcfdPnyZdd9N2/e1OjRo3X9+nVFRkZq586d6tSpk12j+i0i6kbNiLpbiUrStWvXXNeffvppb46EEDZz5kzXyWvt2rVTdna2unTpYvNUCAYLFiyQdO/s2y+//NJ1/bXXXnO98mD9+vXq37+/bTP6MyLqRs2I9urVy+1xmZmZrutjxozx4kQIVfPmzdPatWslSW3atNG+ffvUo0cPm6dCsBg8eLAGDBgg6d4JRKWlpZo7d67rTNxFixZxJm4jODvXjdjYWF27dk1PPfWULl682OAxmZmZGjt2rKqrqzVu3DhlZGT4eMrQ895776ldu3ZKSkpSUlKSHnvsMbtH8qrFixfrk08+kSRFRUVp//79SklJsXeoEJOTk6PMzEzX91yPHj0UERFh91jN6scff9Qrr7wiSRo2bJh+/vlnSZyJ6wki2oC///5bHTp0kCSNHj1au3btct139+5dnT17Vhs2bNDKlStVVVWllJQU7dmzR9HR0XaNHDK6d++uM2fOuL6Oi4tz/XILtrCuWLFCc+bMkSS1bNlSu3bt0ogRI2yeKvRs2rRJU6ZMcX0dGRmpPn36qF+/fkETVsuyFB8fX+tni5dOeSZwn3UvOn78uOv67t273f4rLCkpSVOnTtWMGTMC+gcokBUXF6u4uLjWtnowhHX9+vWaO3euJCk8PFzffvstAfUTd+7c0ZEjR3TkyBHXbYEeVofDofnz5ys9PV2S1KFDB+3cuZOAeiAwnmEfa+ydimq6ffu20tLSAuYHJVQEeli///57TZ8+XZZlyeFwaN26dRo7dqzdY6ERwRDWrl27uq7PnDmTM3E9ZeeLVP3V+PHjXS9Azs3NtfLz8638/HzryJEj1ubNm63ExETX/YMGDbJ7XMuyLOvDDz90zcTFs0tcXJz10ksvWdnZ2XY/fS5ZWVmuF7ZLslauXGn3SG45nU7bn8NAu0RGRlrJycnW7Nmz7X766lm+fLlrzszMTLvHCRhEtAHdunWzJFnt2rVr8P6ysjKrZ8+erm+4o0eP+njC+uz+5RDIl0WLFtn99FmWZVm5ublWq1atXHN9+umndo/UqBs3btj+3AXqJSIiwu6nr54333zTNd+lS5fsHidg8BKXOpxOp86fPy9JSkxMbPCYyMhILVy40PX1li1bfDIbvCMszP4fg6NHj2r06NEqKyuTdO+1ex999JHNUyGU3P8zVtu2bWu9iQwa55+b8zY6efKkqqurJUl9+/Z1e9yYMWP0yCOP6Pbt29qxY4eWL1/uowkbtm/fPtdp6cFs2bJlD/3fcDgc6tatm+tvVaNGjWqGycwVFBQoNTVVt27dkiTNmjVLn3/+ua0zeaJ169aaP3++3WN4XXZ2tsfnSTQmJibG9bfRfv36PfxgzaiiokKFhYWSpD59+tg8TWAhonXU/GFxtxKV7n1qy/Dhw5WZmanLly/rjz/+aPSdjbxtxIgRIXH25q5du2qdhv8gdYOZlJSkxMREv3k5UlFRkZ5//nldv35dkjRx4kR99dVXNk/lmRYtWmjp0qV2j+F1dV/i4omawbwfzS5duvjt6y1PnTql8vJySY0vHlAfEa2jZkQf9M00atQo1xmgu3fvtjWi8P9g1nXlyhUNHz5cV69elXRvd+Obb77xi+1leC7QgtmQpvzeQ21EtI7730xRUVHq1q1bo8eOHDlSDodDlmXpp59+0scff+yDCSEFXjDrunHjhoYPH65Lly5JuvcmEosWLdLp06c9/m/ExcUpJibGSxOiIcEQzIacPHnSdZ2INg3vWFRDdXW1oqOj5XQ6H/ghtfclJSXp2LFjCgsL09WrV/32tYfBYsuWLXryyScDKpgN2bp160O/H2leXp6ee+65ZpoI7uTn56uwsDBogtmQYcOG6ZdfflHLli11+/ZttWjRwu6RAgYr0RrOnj3r+gxRT/81NmrUKB07dkzV1dXKysrSW2+95cUJ8cYbb9g9QrPIz89/qMeHh4ezYvCR3r17N/pxiMHg/kq0R48eBLSJWIkCAGCIMxgAADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADBERAEAMEREAQAwREQBADD0f1uMy4xvJGA6AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pgm = daft.PGM(dpi=DPI, grid_unit=GRID_UNIT, node_ec=NODE_EC)\n", + "\n", + "pgm.add_node(\"r\", \"$R$\", 0, 0)\n", + "pgm.add_node(\"z\", \"$Z$\", 1, 0)\n", + "pgm.add_node(\"x\", \"$\\mathbf{X}$\", 1.5, 0.75)\n", + "pgm.add_node(\"y\", \"$Y$\", 2, 0)\n", + "\n", + "pgm.add_edge(\"r\", \"z\")\n", + "pgm.add_edge(\"z\", \"y\")\n", + "pgm.add_edge(\"x\", \"y\")\n", + "\n", + "pgm.render();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The new variable $R$ represents the random assignment of units to the treatment group. This means that the treatment effect $Z \\rightarrow Y$ can be estimated without bias." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Instrumental Variables\n", + "\n", + "In quasi-experiments, we cannot randomly assign subjects to treatment groups. So confounders $\\mathbf{X}$ will still influence treatment assignment. In the instrumental variable (IV) approach, the causal effect of $Z \\rightarrow Y$ is identifiable if we have an IV that causally influences the treatment $Z$ but not the outcome $Y$." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "tags": [ + "remove-input" + ] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdEAAAEMCAYAAACbY4xqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAB7CAAAewgFu0HU+AAAVMElEQVR4nO3deWwU9f/H8de0RWsVEBRFpRqKolAoVrEgEaN4FFAjChENCkZBPIJGEa/QKBo8UCMqoilQFK2iKJ5osQgIpoIh0qrlUFGwiLdISg96zfcPfp0f23P3w+7O7szzkTTZnZ3dvMzWvvi897O7lm3btgAAQMgS3A4AAEC8okQBADBEiQIAYIgSBQDAECUKAIAhShQAAEOUKAAAhihRAAAMUaIAABiiRAEAMESJAgBgiBIFAMAQJQoAgCFKFAAAQ5QoAACGKFEAAAxRogAAGKJEAQAwRIkCAGCIEgUAwBAlCgCAIUoUAABDlCgAAIYoUQAADFGiAAAYokQBADBEiQIAYIgSBQDAECUKAIAhShQAAEOUKAAAhihRAAAMUaIAABiiRAEAMESJAgBgiBIFAMAQJQoAgCFKFAAAQ5QoAACGKFEAAAxRogAAGKJEAQAwRIkCAGCIEgUAwBAlCgCAIUoUAABDlCgAAIYoUQAADFGiAAAYokQBADBEiQIAYIgSBQDAECUKAIAhShQAAEOUKAAAhihRAAAMUaIAABiiRAEAMESJAgBgiBIFAMAQJQogrlRXV+u0006TZVkBP1OmTGnzfkVFRUpMTGx2vxUrVkQpObzIsm3bdjsEAISiqKhIQ4cOVUNDg3PMsix9/vnnGjp0aLPz9+3bp9NPP11btmwJOD5x4kTNmzcv4nnhXaxEAcSdIUOG6Pbbbw84Ztu2brzxRlVVVTU7f8aMGc0KtEePHnrqqacimhPex0oUQFyqrKxURkaGtm3bFnD87rvv1pNPPulc//rrrzVo0CDV1dUFnLds2TKNHDkyKlnhXZQogLi1evVqDRs2TAf+GUtMTFRRUZGysrJUV1engQMHqqSkJOB+48eP1yuvvBLtuPAgxrkA4tZ5552nm2++OeBYfX29brjhBtXU1Ojxxx9vVqDdu3fX7Nmzo5gSXsZKFEBc27t3r/r166cdO3YEHB83bpyWLFmimpqagONLly7VFVdcEc2I8DBKFEDcKyws1MUXX9zueWPHjtXixYujkAh+QYkC8ISJEydqwYIFrd7erVs3lZaWqlu3blFMBa+jRAF4wp49e5Senq5ff/21xdsXL16ssWPHRjkVvI6NRQA8oXPnzho2bFiLtx1++OFBjXuBUFGiADxhzZo1eu2111q8raKiQnfeeWeUE8EPGOcCiHtVVVXKyMjQjz/+2OZ5n3zyiYYPHx6lVPADVqIA4t706dObFWhSUlKz82666SaVl5dHKxZ8gBIFENfWr1/f7MMTLMvShx9+qD59+gQcLysr07Rp06KYDl5HiQKIWzU1NbrhhhsCvs1FkiZPnqzhw4crLy9PCQmBf+Zyc3O1atWqaMaEh1GiAOLWjBkztGnTpoBjqampmjVrliRp8ODBzTYU2batiRMnqrKyMmo54V1sLAIQlzZu3Oh8yPyBPv74Y40YMcK5XlVVpQEDBuiHH34IOO+OO+7gM3Rx0ChRAHGnrq5OZ511loqLiwOOt/btLF988YXOPffcgG97SUhI0Nq1azVkyJBIx4WHMc4FEHcee+yxZgV67LHH6plnnmnx/HPOOUdTpkwJONbQ0KAbb7xR1dXVkYoJH2AlCiCubNq0SZmZmc2+neXtt9/W6NGjW71fRUWFMjIy9NNPPwUcv/fee/X4449HJCu8jxIFAMAQ41wAAAxRogAAGKJEAQAwRIkCAGCIEgUAwBAlCgCAIUoUAABDlCgAAIYoUQAADFGiAAAYokQBADBEiQIAYIgSBQDAECUKAIAhShQAAEOUKAAAhihRAAAMUaIAABiiRAEAMESJAgBgiBIFAMAQJQoAgCFKFAAAQ5QoANf8/PPPbkcImW3bcZkbkUGJAnBFfn6+0tLS1LlzZ1VXV7sdJyjFxcVKSEhQWlqaKisr3Y6DGGDZtm27HQKAv1RVVSklJcW5fsopp+j77793MVFwLMtyLicmJqqurs7FNIgFrEQBRN2RRx4ZcH3Dhg3uBAlRbm6uc7m+vl6LFi1yMQ1iAStRAFGVn5+va6+91rk+d+5c3XLLLS4mCs2Bq1FJqqioCFhVw18oUQBR03SMK+3fqBNPysvL1alTJ+c6Y11/Y5wLIGqajnH37NnjTpCD0LFjR8a6cLASBRAV8T7GbYqxLiRKFEAUeGGM2xRjXUiMcwFEgRfGuE0x1oXEShRAhHltjNsUY11/o0QBRIwXx7hNMdb1N8a5ACLGi2Pcphjr+hsrUQAR4fUxblOMdf2JEgUQdn4Y4zbFWNefGOcCCDs/jHGbYqzrT6xEAYSV38a4TTHW9RdKFEDY+HGM2xRjXX9hnAsgbPw4xm2Ksa6/sBIFEBZ+H+M2xVjXHyhRAAeNMW5zjHX9gXEugIPGGLc5xrr+wEoUwEFhjNs2xrreRokCMMYYt32Mdb2NcS4AY4xx28dY19tYiQIwwhg3NIx1vYkSBRAyxrihY6zrTYxzAYSsS5cuAdcZ47aPsa43sRIFEBLGuAeHsa63UKIAgsYY9+Ax1vUWxrkAgsYY9+Ax1vUWVqIAgsIYN7wY63oDJQqgXYxxw4+xrjcwzgUQYOXKlaqpqQk4xhg3/IIZ6+7cuVOlpaXRjoYQUKIAHHv37tWoUaOUlZWl4uJiSfvHuPv27XPOmTt3bsAKCuYmTZoUcH3ChAmqrKyUbdvKy8tTenq6Hn74YZfSIRiMcwE4cnNzNXnyZElSUlKS7r33Xs2cOTPgHP5khFfTsa4kjRgxQp988omk/c9DWVmZunfv7kY8tIOVKABJ+8tx7ty5zvW6urpmBcoYN/yajnUlOQUq7X8e5s+fH+1YCBIlCkCStG7dOpWUlLR6e0ZGhg477LAoJvKPESNGtHl7bm4um45iFCUKQJICVqEt+eabb5SVldVm0SI0tm1r4cKF6tevX5vnlZWVadmyZVFKhVDwmigA/f333zrhhBOa7cptSVJSkp599lndeuutUUjmXbW1tbryyiv10UcfBXV+dna2CgoKIpwKoWIlCkALFy4MqkAlqXv37rr44osjnMj7OnTooKuuuqrZhy60Zvny5dq2bVuEUyFUlCjgcw0NDXrppZeCOrdHjx5atWqVTj755Ain8ofrrrtOr7zyStBFGuzzhOhhnAv4XEFBQbsbWyQKNJJeffVVTZgwod23D3Xt2lU7d+5kg1cMYSUK+Fx7G4okCjTSgl2R/vvvv1qyZEmUUiEYrEQBH9uxY4fS0tLU0NDQ6jkUaPQEsyIdPHiwvvzyyyimQltYiQI+lpubS4HGkGBWpOvWrdPGjRujmAptoUQBn6qpqWnzk3AoUHcEU6QvvvhiFBOhLZQo4FPvvvuu/vzzzxZvo0Dd1V6R5ufn8xGMMYISBXyqtQ1FFGhsaKtIKysrm31tGtzBxiLAh0pLS1v8qDkKNPa0ttmoT58+Ki0tDfo9pogMVqKAD7X0mhoFGptaW5Fu3rxZn3/+uUup0IgSBXxm7969zUaBFGhsa61I2WDkPkoU8JnXX39d5eXlznUKND60VKRLly7V77//7mIqUKKAz6xcudK5TIHGl6ZFWldXp7Vr17qcyt/YWAT4TG1trWbOnKlFixbp008/pUDj0KuvvqqcnBzNnTtXI0eOdDuOr1GigE/V1dUpKSnJ7RgwxPMXGyhRAAAM8ZooAACGKFEAAAxRogAAGKJEAQAwRIkCAGCIEgUAwBAlCgCAIUoUAABDlCgAAIYoUQAADFGiAAAYokQBADBEiQIAYIgSBQDAECUKAIAhShQAAEOUqKQXX3xRlmXJsizl5+c7x88//3xZlqUuXbqE/Jjjx493HvPtt98OZ1x40KxZs5zfF5Of2bNnu/2fgDiye/dudezYUZZlKTU1VXV1de3ep76+XiNGjHB+5954440oJI19lKik4uJi5/Lpp5/uXO7Xr58k6b///tOuXbuCfrySkhKnjIcMGaIxY8aEJSe8a8OGDQd1//79+4cpCfygS5cumjRpkiRp586dQf1Df+rUqSooKJAkTZ8+Xddcc01EM8YLy7Zt2+0Qbhs0aJC++uorJScna+/evUpMTJQkvfTSS7rlllskSStWrNAFF1wQ1ONlZ2fr008/lSR9+eWXGjx4cGSCwzN++uknVVZWBnVueXm5xo4dq7KyMknSGWecobVr1yolJSWSEeExZWVl6tWrl2prazVo0CCtW7eu1XPnzZunm266SZI0evRoLVmyRJZlRStqbLN9rr6+3k5JSbEl2QMHDgy4bc2aNbYkW5L93HPPBfV4hYWFzn3Gjh0bicjwsaqqKvu8885zfsf69Olj//XXX27HQpwaP36887tUVFTU4jmrVq2yO3ToYEuyMzMz7YqKiiinjG2+H+du3brVWQEcOMqVpPT0dOfypk2b2n0s27Z1zz33SJIOPfRQPfbYY+ELCt+rra3VmDFjtHr1aklSz549tWLFCh199NHuBkPcuueee5wVZUuvq2/btk1jxoxRbW2tunfvrg8++ICJRxO+L9EDXw/NzMwMuK1r16467rjjJEmbN29u97Hy8/O1ceNGSdLtt9+unj17hi8ofK2hoUHXXXedli1bJkk6/vjjtWLFCh1//PEuJ0M8S09P18iRIyVJ77zzjn755Rfntj179uiyyy7TP//8o+TkZL3//vvq0aOHW1FjFiXayqaiRo2r0fZWovv27VNOTo4k6aijjtIDDzwQtozA5MmT9eabb0qSjj76aBUWFiotLc3lVPCCxulZfX295syZ41y++uqrncVDXl6esrKyXMsYyyjR/ytRy7KUkZHR7PbGEv3rr7/0999/t/o4c+bM0fbt2yVJDz74oI488shwR4VPTZ06VfPnz5ckderUSQUFBerbt6/LqeAV5557rrP5cd68eaqoqNBdd93l7MTNyclhJ25b3H5R1m3HHnusLck+5ZRTWrw9NzfXeeF9zZo1LZ6ze/duu2vXrrYku3fv3nZNTU0kI/vaHXfcYT/yyCP2xx9/bP/xxx9ux4m4hx56yPn9S0lJsdeuXet2JN8pLCy0b7vtNjsvL88uKSmxa2tr3Y4UdkuXLnV+zy644ALn8ujRo+2Ghga348W0JLfKOxb8/vvv+uOPPyS1PMqV/v+9otL+ke7QoUObnfPoo4/q33//lSQ98cQT6tChQ/jDQpJUUFCgrVu3OtdTU1N15plnBvwcc8wxLiYMn9mzZ+uhhx6SJB1yyCFaunSpzjnnHHdD+dCuXbv0wgsvONeTk5M1YMAADRw40Pmd69u3r5KS4vfP6ahRo3Tqqadq69at+uyzzyTtf+vUokWLeCtLO+L3WQ+Dxk1AUvNNRY0O3KHb0uaisrIyPf/885L2j0VGjRoV3pBoU1lZmcrKyvTee+85x7xQrHl5ebrrrrskSYmJiXr99deVnZ3tcipIUnV1tdavX6/169c7x+K9WC3L0rRp0zRx4kRJ0nHHHaf333+fnbhBiI9nOELa21Qk7X8NqkePHtq5c2eLm4umT5+u6upqWZalp59+OkJJEYp4L9a33npLkyZNkm3bsixLCxYs0OjRo92OhTZ4oVh79erlXJ48eTI7cYPl9jzZTVdddZUz+9+1a1er52VnZ9uS7BNOOCHgeElJiZ2QkGBLsseNGxfpuG26//77nf8WfoL7SU1NtUeNGmUXFha6+twdaNmyZc4b2xXCh3y4obKy0vXnMN5+kpOT7UGDBtlTpkxx++lr5plnnnFyvvfee27HiRu+LtHevXvbkuxjjjmmzfOmTp3q/HLt2bPHOd5YrsnJyfaOHTsiHbdNbv9xiOefnJwcV5+7RqtXr7YPO+wwJ9fMmTPdjtSm3bt3u/7cxetPUlKS209fM9dff72Tb/v27W7HiRu+fYtLZWWlfvzxR0mtj3IbtfS66MqVK7V8+XJJ0p133qkTTzwxMkERcQkJ7v9vsGHDBl122WWqqqqStP+9e7zXGNHU+PJWly5ddNJJJ7kbJo7E5nA+CkpKStTQ0CCp9U1FjZp+/F9WVpamTZsmSerWrZvuu+++yAUNUkFBgbOrzsuefPLJg34My7LUu3dv57WqSy65JAzJzJWWlmr48OEqLy+XJN1888164oknXM0UjMMPP9z5/8DLCgsLA/ZPmOrcubPz2ujAgQMPPlgY1dbWOns+BgwY4HKa+OLbEg1mU1Gjvn37yrIs2batzZs364033tDXX38tSZoxY4Y6deoUwaTByc7O9sXuzQ8++CDgLS7taVqYZ555pjIzM9WxY8cIpgzetm3bdNFFF+mff/6RJI0bNy7g7RSxrEOHDpo1a5bbMSJu0aJFmjBhQkj3ObAwG0szLS0tZt8usnnzZtXU1Ehq/+8hAlGiav+X5ogjjtBJJ52k7du3q7i42PnuvT59+jjfyQf3xXphNvXrr7/qwgsv1G+//SZJuvzyy/Xyyy/HxHgZwYu3wmxJKH8PEcj3JZqSkqLevXu3e356erq2b9+uwsJC59isWbNidru618VbYTa1e/duXXjhhc5HRZ522mnKycnRli1bgn6M1NRUde7cOUIJ0RIvFGZLSkpKnMuUaGh82QANDQ367rvvJEn9+/cP6l/+6enpzjdoSNKwYcN06aWXRiwjWjZ9+nSdeOKJcVWYLVm+fHlAYW7ZsiXk18mKiop09tlnhzsamsjMzNTixYs9U5gtaVxUHHLIIXwuc4h8WaLff/99q98h2poDNxclJCToqaeeikQ0tOPaa691O0JYfPvttwd1/8TERFYMUdK/f3/179/f7RgR1bgS7du3Lx9bGiLLtm3b7RAAAMQjdjAAAGCIEgUAwBAlCgCAIUoUAABDlCgAAIYoUQAADFGiAAAYokQBADBEiQIAYIgSBQDAECUKAIAhShQAAEOUKAAAhihRAAAMUaIAABiiRAEAMESJAgBgiBIFAMAQJQoAgCFKFAAAQ5QoAACGKFEAAAxRogAAGKJEAQAwRIkCAGCIEgUAwBAlCgCAIUoUAABDlCgAAIYoUQAADFGiAAAYokQBADBEiQIAYIgSBQDAECUKAIAhShQAAEOUKAAAhihRAAAMUaIAABiiRAEAMESJAgBgiBIFAMAQJQoAgCFKFAAAQ5QoAACGKFEAAAxRogAAGKJEAQAwRIkCAGCIEgUAwBAlCgCAIUoUAABDlCgAAIYoUQAADFGiAAAYokQBADBEiQIAYIgSBQDAECUKAIAhShQAAEOUKAAAhihRAAAMUaIAABiiRAEAMESJAgBgiBIFAMAQJQoAgKH/AcCXPLig+oBWAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pgm = daft.PGM(dpi=DPI, grid_unit=GRID_UNIT, node_ec=NODE_EC)\n", + "\n", + "pgm.add_node(\"iv\", \"$IV$\", 0, 0)\n", + "pgm.add_node(\"z\", \"$Z$\", 1, 0)\n", + "pgm.add_node(\"y\", \"$Y$\", 2, 0)\n", + "pgm.add_node(\"x\", \"$\\mathbf{X}$\", 1.5, 0.75)\n", + "pgm.add_edge(\"iv\", \"z\")\n", + "pgm.add_edge(\"x\", \"z\")\n", + "pgm.add_edge(\"x\", \"y\")\n", + "pgm.add_edge(\"z\", \"y\")\n", + "\n", + "pgm.render();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's try to get some intuition of why having the $IV$ helps:\n", + "* The presence of $\\mathbf{X}$ is a confounder because it influences both $Z$ and $Y$.\n", + "* But the $IV$ helps overcome this confounding because it is not influenced by $\\mathbf{X}$.\n", + "* Any association between the $IV$ and $Y$ must be through the treatment $Z$.\n", + "* This means that the $IV$ can be used to estimate the causal effect of $Z \\rightarrow Y$, without being confounded by $\\mathbf{X}$. Informally, the $IV$ causes some variation in the treatment $Z$ that is not due to $\\mathbf{X}$, and this variation can be used to estimate the causal effect of $Z \\rightarrow Y$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Readers are referred to {cite:t}`steiner2017graphical,cunningham2021causal` or {cite:t}`huntington2021effect` for a more in-depth discussion of the IV approach from the causal DAG perspective." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Interrupted Time Series\n", + "\n", + "A causal DAG for interrupted time series quasi-experiment is given in Chapter 17 of {cite:t}`huntington2021effect`, though they are labelled as [Event Studies](https://theeffectbook.net/ch-EventStudies.html). These kinds of studies are suited to situations where an intervention is made at a given point in time at which we move from untreated to treated. Typically, we consider situations where there are a 'decent' number of observations over time. Here's the causal DAG - note that $\\text{time}$ represents all the things changing over time such as the time index as well as time-varying predictor variables." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "tags": [ + "remove-input" + ] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdEAAAEzCAYAAABuTzkSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAB7CAAAewgFu0HU+AAAhKklEQVR4nO3deXRU9d3H8U8IkIUACViMQTYXkKWCAuKGSUGFChxQQKB4RDQgiKiIwnlkK1QqxyoFxIK1AoeWTbBQVESUTUEQREREK8hOZE9igSSEJN/njzS3CUkguWQy2/t1zhyTe+/M/G6IvLm/uXMnxMxMAACg1Cp4ewAAAPgrIgoAgEtEFAAAl4goAAAuEVEAAFwiogAAuEREAQBwiYgCAOASEQUAwCUiCgCAS0QUAACXiCgAAC4RUQAAXCKiAAC4REQBAHCJiAIA4BIRBQDAJSIKAIBLRBQAAJeIKAAALhFRAABcIqIAALhERAEAcImIAgDgEhEFAMAlIgoAgEtEFAAAl4goAAAuEVEAAFwiogAAuEREAQBwiYgCAOASEQUAwCUiCgCAS0QUAACXiCgAAC4RUQAAXCKiAAC4REQBAHCJiAIA4BIRBQDAJSIKAIBLRBQAAJeIKAAALhFRAABcIqIAALhERAEAcImIAgDgEhEFAMAlIgoAgEtEFAAAl4goAAAuEVEAAFwiogAAuEREAQBwiYgCAOASEQUAwCUiCgCAS0QUAACXiCgAAC4RUQAAXCKigIetWbNGNWvW1IMPPigz8/ZwAJQhIgp42NKlS5WcnKxly5YpOTnZ28MBUIYqensAQKAbPHiwtm/frrZt26pmzZreHg6AMhRizC8BZWL37t1auHChxo4d6+2hACgnTOcCZWTSpEnKycnx9jAAlCMiCpSBlStXas6cOd4eBoByRkSBK7R69Wr17NmTM2+BIEREgSswefJk9erVS+np6ZJyp3Sjo6MVHR2t2bNnKycnRytXrtSDDz6oG2+8sdD9jx07plGjRqlGjRo6cOCAJOmdd95RkyZNVKVKFbVt21bbt293tt+2bZs6d+6satWqKTY2VmPGjCl2CjktLU0TJ05UixYtVL16dUVHR6tDhw5av3592f8ggGBlAK5YfHy8SbJx48Y5y7Zs2WIdOnSwihUrmiSrV6+es+6XX36xxMRECwsLM0kmyfbt22cDBw60KlWqWGxsrLO8Vq1advLkSXv//fctPDzcrr32WgsPD3fWT5w4sdB4jhw5Ys2aNbOXX37ZfvnlF8vMzLS3337bwsLCrEKFCjZ79mzP/1CAIMCRKOAhrVu31sqVKzVs2LBC66pUqaLx48frtddec5ZNmDBBzZs31+nTp3X06FEtW7ZMISEhOnHihF544QX95S9/0ddff63Dhw8rJSVFnTt3liRNmzatwGObmXr27KmOHTtq1KhRqlatmipVqqTExESNHDlSOTk5Gjx4sH7++WfP/gCAIEBEAQ9r0KBBoWWhoaGKi4tTQkKCs6x///566qmnFBYWJknq2rWrbr/9dknSoUOH9P7776tx48aSpPDwcL3wwguSpOPHjyslJcV5nA8++ECbNm1SYmJioee94447JEkZGRlavHhx2ewgEMS42ALgYeHh4cWui4yMdL6uW7duofX169fXpk2bVLduXYWGhhZYFxcX53x99uxZxcTESMq9QpIktWnTptDj5eTkOJHet29fKfYCQFGIKOBhISEhxa6rUOHSk0GVK1cudl2lSpWcr7Ozs52vd+7cKSn3pKVLBRzAlWM6FwgwqampkqSkpCTvDgQIAkQUCDB5U8Tr1q275Hb5j14BuENEgQDTqFEjSbln7WZlZRW5zaFDh/Tyyy+X57CAgEREgTKQ99pmRkaGl0cidenSRZL07bffaujQoYWupJSTk6NnnnlG9957rzeGBwQUIgqUgapVq0qSNmzYoOzsbP3000+aMGGCpNwrB0lyrmqUX/5lRQU470gyMzOz0Lr807H5H6d3795q2rSpJGnmzJmKj4/Xe++9px07dmj58uW6//77deHCBd11112l3k8ABRFRoAzkvZ9z48aNiouLU9u2bdWvXz+dP3/eecvJiRMnCrxOaWZatGiR8/2SJUsKHDUeOXLE2X7Tpk0FThQyMy1cuND5fvHixc7l/ypVqqSlS5c6b5n5/PPP1aNHD7Vo0UJdu3bVkSNHuFg+UFa8er0kIECkpaVZ3759LTIy0m677TbbunWrrV271rnkX/5b+/btLSsrq8Al//JuYWFhdvDgQRs4cGChdSEhIdavXz87ePCgVa5cucj75nfq1Cl7/vnnrX79+la5cmWrU6eOPffcc5acnOylnxIQePhQbgAAXGI6FwAAl4goAAAuEVEAAFwiogAAuEREAQBwiYgCAOASEQUAwCUiCgCAS0QUAACXiCgAAC4RUQAAXCKiAAC4REQBAHCJiAIA4BIRBQDAJSIKAIBLRBQAAJeIKAAALhFRoISaNm2qVatWeXsYHpWUlKSwsDCdP3/e20MB/EKImZm3BwH4upiYGKWmpkqSPv74Y91///3eHZAHJCUl6dprr3W+z8jIUFhYmBdHBPg+jkSBy8gfUEkaNmyY9wbjQbNmzSrwfXh4OEekwGVwJApcwsUBjYmJUXJysvcG5GGDBw/WzJkzCyzjiBQoHkeiQDGCLaCSNGPGDA0aNKjAMo5IgeIRUaAIwRjQPIQUKDkiClwkmAOah5ACJUNEgXwI6P8QUuDyiCjwXwS0MEIKXBoRBURAL4WQAsUjogh6BPTyCClQNCKKoEZAS46QAoURUQQtAlp6hBQoiIgiKBFQ9wgp8D9EFEGHgF45QgrkIqIIKgS07BBSgIgiiBDQskdIEeyIKIICAfUcQopgRkQR8Aio5xFSBCsiioBGQMsPIUUwIqIIWAS0/BFSBBsiioBEQL2HkCKYEFEEHALqfYQUwYKIIqAQUN9BSBEMiCgCBgH1PYQUgY6IIiAQUN9FSBHIiCj8HgH1fYQUgYqIwq8RUP9BSBGIiCj8FgH1P4QUgYaIwi8RUP9FSBFIiCj8DgH1f4QUgYKIwq8Q0MBBSBEIiCj8BgENPIQU/o6Iwi8Q0MBFSOHPiCh8HgENfIQU/oqIwqcR0OBBSOGPiCh8FgENPoQU/oaIwicR0OBFSOFPiCh8DgEFIYW/IKLwKQQUeQgp/AERhc8goLgYIYWvI6LwCQQUxSGk8GVEFF5HQHE5hBS+iojCqwgoSoqQwhcRUXgNAUVpEVL4GiIKryCgcIuQwpcQUZQ7AoorRUjhK4goyhUBRVkhpPAFRBTlhoCirBFSeBsRRbkgoPAUQgpvIqLwOAIKTyOk8BYiCo8ioCgvhBTeQEThMQQU5Y2QorwRUXgEAYW3EFKUJyKKMkdA4W2EFOWFiKJMEVD4CkKK8kBEUWYIKHwNIYWnEVGUCQIKX0VI4UlEFFeMgMLXEVJ4ChHFFSGg8BeEFJ5AROEaAYW/IaQoa0QUrhBQ+CtCirJERFFqBBT+jpCirBBRlAoBRaAgpCgLRBQlRkARaAgprhQRRYkQUAQqQoorQURxWQQUgY6Qwi0iiksioAgWhBRuEFEUi4Ai2BBSlBYRRZEIKIIVIUVpEFEUQkAR7AgpSoqIogACCuQipCgJIgoHAQUKIqS4HCIKSQQUKA4hxaUQURBQ4DIIKYpDRIMcAQVKhpCiKEQ0iBFQoHQIKS5GRIMUAQXcIaTIj4gGIQIKXBlCijxENMgQUKBsEFJIRDSoEFCgbBFSENEgQUABzyCkwY2IBgECCngWIQ1eRDTAEVCgfBDS4EREA8C7776rNWvWFFpOQIHyVdKQ7t27V6+99lp5Dg0eEmJm5u1BwL1jx46pSZMmysjI0AcffKB27dpJIqCANw0ePFgzZ84ssCwjI0NhYWHau3evEhISlJSUpA0bNujOO+/00ihRFjgS9WNmpkGDBiklJUXp6enq3Lmz1qxZQ0ABLyvuiPT7779XQkKCjhw5IjNT//79lZ6e7qVRoixwJOrH5s2bp0ceeeSS2xBQwHuKOiK92PDhw5na9WNE1E/lTeOmpKQUuw0BBbzvciENCQlhWtePMZ3rh/JP417KkiVLymlEAIrzwgsvXHI907r+jYj6ofnz5+tf//rXZbfLe40UgHfknUR0Obt379aYMWM8PyCUOaZz/UxJpnHzi4iIKHDWLoDykRfQI0eOlGh7pnX9E0eifqSk07j5RURE6PTp0x4cFYCinDp1ShERESXenmld/0RE/UhJp3HzdOvWTbt27VLPnj09OCoARWnTpo2++eYbDRs2TCEhISW6D9O6/ofpXD9RmmncGjVqaPr06erdu3eJ/+cF4DkbNmzQ448/rj179lx2W6Z1/QtHon6gNNO4eUefffr0IaCAj7j77rtLfFTKtK5/IaJ+oCTTuDVq1ND8+fP1z3/+U7GxseU0MgAlFRkZqcmTJ+uzzz7TjTfeeMltmdb1H0zn+riSTON269ZNM2bMIJ6An0hLS9Po0aM1ZcoUFfdXMNO6/oEjUR92uWlcjj4B/1SSo1Kmdf0DEfVhl5rG5bVPwP9d7rVSpnV9H9O5Pqq4aVzOvAUCU3Fn8DKt69s4EvVBxU3jcvQJBK7ijkqZ1vVtRNQHnT17Vjt37nS+57VPIDgU91rp0aNHtX//fi+ODMVhOtdHHTp0SL/5zW908803c+YtEITyzuCdNWuWPvzwQ911113eHhKKQEQvsmrVKk2ZMkVff/21srOz1bZtW02cOFGNGzcu97GcOXNGUVFRTN0CQezMmTOqWrWqt4eBYhDRfCZOnKjVq1dr2bJlysnJ0YMPPqh169apSpUq2rZtmxo1auTtIQIAfAivif7X2rVrNXr0aA0bNkzVqlVTdHS0Fi1apLi4OJ07d07btm3z9hABAD6GiP7X66+/LkmqW7eus6xWrVrauHGjZs+erV69emnUqFE6cOCAl0YIAPA1TOcq9xTyKlWqKD09XXv27NENN9xQaJvU1FRdd911+vrrr1W/fv3yHyQAwOdwJCopJSXFeQ9WaGhokds89dRTpfowbABA4COiyj37Lc/FZ8JmZ2frmWee0YIFC8p7WAAAHxdwET1//rymTp2qVq1aKTY2VlFRUWratKnGjh2rjIyMAtuOHz9e0dHR+vWvf+0su/nmmxUdHa3o6Gh9+OGHuu+++zRnzpwi158/f95ZbmaaN2+e4uPjVatWLUVGRqpVq1b629/+VmiMP/zwg5599llFR0frwIED2r17t+655x5Vq1ZN48aNK/sfCgDAMyyAZGRk2O23326VK1e2lStXmpnZnj17rHnz5ibJHn744SLvt3//fpNkkmz//v2lXp+ZmWk9evSwnj172qFDh8zMbP369Va3bl2TZI8//riZmZ08edIeeughq1ixovN4mzdvtkaNGjnLKlSoYJmZmWXzAwEAeFRARXT27Nkmye68884Cy5cuXeoEKi0trdD9rjSiI0eOtNtvv92ysrIKLF+zZo1zv/fee89Z/t577znLO3XqZD/++KPt2rXLunTpYiNHjnS38wCAclexnA98Pers2bOSct+akt9NN90kScrJyVFycrJq165dZs958uRJTZ06VVOmTCl0UtIdd9zhfP33v/9dDz30kCSpadOmzvKHHnpIDRs2lCQtX768zMYFAPC8gIroE088obCwMN1///3OspycHG3fvt35Pisrq0yf86OPPlJGRoaGDx+ukSNHFlofFhYmSQXeX1q5cmXn63bt2pXpeAAA5SegIhoREaEBAwZIkk6cOKEZM2ZoxYoVat68ubONlfHbYvM+bWXJkiXq2LFjie7DtXABIDAE3Nm56enpGjFihO6++241bNhQGzdu1EsvveSx50tNTZUkJSUleew5AAC+KaAievLkSbVu3VorVqzQ5s2b1adPH1Ws6NmD7cjISEnSunXrLrlddna2R8cB//bqq68qJCTE9W3KlCne3gX4kZSUFFWtWlUhISGqU6dOiV7mys7O1m9/+1vnd473zucKqIgOHDhQu3bt0qhRo1SjRo1yec68T3ZZsmSJDh48WOQ2OTk5evLJJ8tlPPBPX3311RXdP/97nYHLiYmJcV76OnLkiJYsWXLZ+wwfPlwrV66UJI0ePVp9+vTx6Bj9RUC9Jrpq1SpJxV+6T8oN2sXyHyUW9ZpphQr/+7fGxRdseOCBBxQaGqqMjAw9/PDDWrVqlapXr15gmz//+c9q0KBByXYCQWnSpEkaO3ZsibY9c+aMevXqpcOHD0uSbr311gJnggMlMWzYME2fPl0XLlzQlClT1Lt372K3ffvttzV16lRJUvfu3TVhwoTyGqbv8/Z7bMpS7dq1TZK1bNnSTp8+bWZmmzdvtvj4eOd9mRs2bLDNmzfbxo0bnftt3brVWb9ly5ZCj5ucnOysnzlzppmZvfbaa7Z9+3YzMxswYICzvn79+jZz5kz7+uuvbe3atTZo0CCLi4uz1NRU5/H27NnjbL9t2zYP/kQQaNLT0y0hIcH5/WncuLGdPHnS28OCn3r00Ued36UvvviiyG3Wrl1rlSpVMkl2yy232Llz58p5lL4toCL6/PPPO78Q4eHhds0111iDBg1s3bp1zvLq1atb+/btnasCpaenF/hFeuyxxywjI6PQYzds2NC5YENcXJx1797dWXfu3Dlr37698xj5b1FRUQWCnZWVZWPHjnXW9+rVy86ePev5Hw78XmZmpnXq1Mn53WnQoIElJSV5e1jwY999952FhIQUe0W3n376yWrWrGmSLDY21g4fPuyFUfq2gIpoWlqaDRkyxK666iqrUaOGDRw40FJSUszM7LHHHrOoqCjr37+/E6033njD+RdW/lvFihWtdu3aBR5769at1qxZM6tWrZolJibamTNnCqy/cOGCTZ061Zo3b25hYWFWo0YN69Gjh33//ffONqdPny5wyb+8W4UKFZzLFAJFyc7Otl69ejm/M3FxcbZ3715vDwsBIO8fZqGhoXbw4EFneWpqqjVu3Ng5KPnyyy+9OErfxeeJAn5gwIABzocZXHXVVVq/fr2aNGni5VEhEHz22WeKj4+XJL344ot69dVXlZ2drc6dOzsnEs2fP58TiYpBROFXnnvuOV111VVq2bKlWrZsWegSj4Fo+PDhmjx5siSpWrVqWrNmjVq2bOnlUQWPTz/9VMuWLXN+55o0aeLxt86VtzvuuEObN29WdHS0jhw5opdeeknTpk2TJI0ZM4YTiS6BiMKv3HTTTfrxxx+d7+vUqeP85RaIYR0/frx+//vfS8p9T/LHH3+su+++27uDCjJz585Vv379nO/Dw8PVvHlztWrVKmDCunTpUufa3u3bt9fq1asl5Z6Ju3jxYq6ydglEFH7l4ogWJVDCOmXKFA0bNkxS7vWWly9frg4dOnh5VMHn4ogWxd/DamZq3Lhxgf+3br31Vn3++efOBWVQNCIKv1KSiBbF38I6a9YsJSYmyswUGhqqRYsWqXv37t4eVlAqSUSL4m9hfeedd5SYmChJuuaaa7RlyxZde+21Xh6VH/DO+Uwoa//3f/9X5FtsuBV/q1OnjnXr1s0++eQTb//xFbBo0SKrUKGCSbKQkBCbM2eOt4dUpLS0NK//GfrbLTw83Nq0aWNDhw719h9fIWvXrnXG+fvf/97bw/EbAXXZv2D2yiuveHsIfufw4cNatmyZPvvsM28PxbFixQo98sgjzpW1pk6d6uooqDycP3/e20PwOxkZGfryyy81Y8YMbw+lkG+++cb5ukWLFl4bh78hogh6+S/r6E3r169Xjx49dOHCBUnSxIkTNXToUC+PCsFix44dztdEtOR8c3IepbZy5UrnjLpA9qc//emKHyMkJEQNGzZ0Xqvq1KlTGYzsynz11Vfq0qWL0tPTJUkjRozw6Ef4lYUqVaroxRdf9PYwPO6TTz4pcJTmVvXq1Z3XRlu1anXlAytjefsYExOjevXqeXcwfoQTi+BXSnti0cXBbNmypW655RZVrVrVg6MsnV27dik+Pl6nT5+WJA0aNMgnp/uClZsTi/IHMy+a1113nc++VeTChQuKiopSZmamEhIStHbtWm8PyW9wJIqA4Q/BvNjevXt13333OQHt27ev3nzzTS+PCqXhb8Esyg8//KDMzExJTOWWFhGFX/LHYF4sKSlJ9957r44ePSpJ6tq1q+bMmeMzr9GisEAIZlE4qcg9Igq/Mnr0aNWtW9fvgnmxlJQU3XvvvTpw4ICk3GnqMWPG6N///neJH6NOnTqFPrsWZe+WW27RwoULAyaYReGkIvfK7TXR3bt3a+HChSX+4GFvGDVqlAYMGKD69et7eyhlLjU1VaNHj9b06dO9PRRIWrhw4RVf0PuLL77gw7hRJtq3b681a9aocuXKOnv2rCpVquTtIfmNcps3mjRpkvPeN1+Umpoa0CdzTJs2TWfPnvX2MPBfO3fuvKL7h4aGcsSAMpN3JNqkSRMCWkrlEtGVK1dqzpw55fFUrj311FNKSUnx9jA8YseOHVyMwcdMnDhRlvt5vq5uWVlZioiI8PZuIECcOnVKZqbt27d7eyh+x+MRXb16tXr27ClffSdNdna2nnnmGS1YsMDbQ/GIHTt2qGPHjsrIyPD2UAAg4Hg0opMnT1avXr2cN5BPmjRJ0dHRio6O1uzZs3Xo0CGNHj1asbGxWrdunY4dO6YuXbqoatWqzsW38yQnJ2vkyJFq0qSJqlSpopo1a6p79+4FXhDPb/ny5WrXrp3q1auniIgI1a9fX4mJiTp8+LCzzblz53TfffcVOEq++eabnTHmXdZs3759Gjp0qKpUqSIpN7yvvvqqrrvuOlWtWlUPPPCA9u/f7zzGp59+qoSEBEVFRalevXp64403iv0ZlWa/1q5dqy5duqhdu3aSpKNHj+rxxx9XzZo1VatWLY0YMaLAlPmCBQvUsWNH5+0T8+fPd/Zt/PjxxY4JAFBC5XGB3vj4eJNk48aNMzOzCxcu2COPPGKRkZHOBY9XrFhhrVu3toiICGfZ3r17zczsu+++sxtuuMFmzpxp586ds3Pnztkrr7xiISEhFh4ebh9//HGB5xs2bJhJsgkTJlhWVpalpKTY7373O5Nk9evXt3PnzhXYfv/+/c5z7t+/31l+5MgR69mzp4WGhjrrz58/b506dbLq1atbzZo1neWNGjWy8+fP21tvvWUVK1a0OnXqWMWKFZ318+bNK/RzKel+rV271jp06OA8Vnx8vO3cudPi4uLs6quvtqpVqzrr/vCHPxR6nn79+pkk69ev3xX8KQIALuaViObZtm2b85d/QkKCbdiwwQ4fPmy9e/e2xx57zLKzs+3s2bN2/fXX2/Tp0ws97qOPPmqS7Oqrr7a0tDQzM9u3b5/zmJmZmc62ycnJFhISYpJs5cqVBR6nuIimpaXZiRMn7MUXX3TWDxkyxBYsWGBZWVlmZjZt2jRn3cCBA+3hhx+2Q4cOmZlZSkqKtW7d2iRZmzZtCjxnafYrb9+6detmkuymm26yBx54wL744gszy/1HSZcuXUySxcbGFno8IgoAnuHVd3U3bdrU+bpt27a66667dO2112rBggWaPXu2KlSooFmzZungwYNFXnYr7/T+48ePa9WqVZJyp2glKTo6usBZZjExMc7nR546dapE44uIiNCvfvUr3XPPPc6yESNGqHfv3goNDZUkPf3004qLi5OU+wkNixYtUp06dZwxDBkyRFLuFUHyK81+5Z1A0rBhQ0m508mLFi1ytqtYsaKefvppSdKxY8eUmppaov0DAFwZr15soXLlys7Xea/zXWzp0qXKzs4u8sNhs7OzFRYWJin38mmS1KxZMy1evFi1a9cusO2OHTuUlZUlSc5/Syr/J7vXrVu3wLqQkBDVq1dPP//8sxo0aFDovnmBvfjtJaXdL0nOsri4OEVFRRW4T164JenMmTOKjo4uya4BAK6AVyNakit/7Ny5U7Vq1dKxY8dK/Lg9evSQlHtk+I9//ENz587VDTfc4Ky3Up4pfLnLsOX/x8DF8o6GL36PrJv9utQ48r/dITs7u8SPCQBwz+cv0pmamqrTp0+X+i0a77zzjpo0aaLjx4/r/fff16xZswodvXmT2/0CAPgOn49oZGSksrKytHHjxktul3f0lZOToz59+ujZZ5/VkiVLNGrUKJ+8vmhp9wsA4Ht8PqKNGjWSJL3++uvFbrNp0ybNmjVLkvTGG29o4cKF6tu3r2699dZyGaMbpd0vAIDvKZeI5r2W52bqskuXLpKkjz76qMhL16Wnp2v48OHq0KGDJDln6eadPVuUi1+fzP9aY3lNr5Z2v67Elfz8AQDFK5eI5n1k1YYNG5Sdna2ffvpJEyZMKHCWbHJycpH3HTJkiGJjYyVJL730krp27aoPP/xQ33zzjd599121bdtWt9xyi3PWbExMjCRp3rx5+vbbbyVJhw4d0hNPPKGkpCRJuVf6OX78uObPn19gfJK0fv16SblHiHmfsZd3xSWp6BDl7Ufeh9rml386Nv/jlHa/JOk///mPJDlXUipO/ufJv39btmxRenq6Tpw4oeHDh1/yMQAAJVAeb0b94x//6FyQoFatWhYbG2sHDhywt99+21l+zz332KlTp4q8/6ZNmywmJsbZNv/trrvuci5GYGa2fPlyZ12FChWsTp06FhUVZfPnz7cWLVqYJAsLC7Mbb7zRuSiCmVnDhg2d+8TFxVn37t3NzCwrK8v69+/vPObcuXMLjO3bb7+1qKgok2StW7e2lJQUZ11mZqY9+eSTzn1nz57ter9SU1OtWbNmJsmioqLs+++/L/BYr7/+unPfCRMmWE5OjrNu/vz5zroaNWpYTEyMbd26tWR/eACAYpVLRNPS0qxv374WGRlpt912m23durXIeISEhNhbb71V5GMcOnTIEhMTLS4uzipXrmzXX3+9jRs3rkBo8rz55ptWv359i4yMtI4dO9quXbvMzOyvf/2rRUVFWUJCgu3Zs6fAfbZu3WrNmjWzatWqWWJiop05c8b27dtnlSpVKjTOiIgIM7MCl+LLvw9jx461zz//vMBl//JuDRs2LPV+zZ0717naUv5bQkKCJSUlWXh4eKF1YWFhlpGRYWZm2dnZ9txzz1nVqlWtadOmha7YBABwp9w+lBsAgEDj82fnAgDgq4goAAAuEVEAAFwiogAAuEREAQBwiYgCAOASEQUAwCUiCgCAS0QUAACXiCgAAC4RUQAAXCKiAAC4REQBAHCJiAIA4BIRBQDAJSIKAIBLRBQAAJeIKAAALhFRAABcIqIAALhERAEAcImIAgDgEhEFAMAlIgoAgEtEFAAAl4goAAAuEVEAAFwiogAAuEREAQBwiYgCAOASEQUAwCUiCgCAS0QUAACXiCgAAC4RUQAAXCKiAAC4REQBAHCJiAIA4BIRBQDAJSIKAIBLRBQAAJeIKAAALhFRAABcIqIAALhERAEAcImIAgDgEhEFAMAlIgoAgEtEFAAAl4goAAAuEVEAAFz6f5x+AzkljiY5AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pgm = daft.PGM(dpi=DPI, grid_unit=GRID_UNIT, node_ec=NODE_EC)\n", + "\n", + "pgm.add_node(\"a\", \"after\\ntreatment\", -1, 0)\n", + "pgm.add_node(\"z\", \"$Z$\", 0, 0)\n", + "pgm.add_node(\"y\", \"$Y$\", 1, 0)\n", + "pgm.add_node(\"t\", \"time\", 0, 1)\n", + "\n", + "pgm.add_edge(\"a\", \"z\")\n", + "pgm.add_edge(\"t\", \"a\")\n", + "pgm.add_edge(\"t\", \"y\")\n", + "pgm.add_edge(\"z\", \"y\")\n", + "\n", + "pgm.render();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What we want to understand is the causal effect of the treatment upon the outcome, $Z \\rightarrow Y$. But we have a back door path between $Z$ and $Y$ which will make this hard, $Z \\leftarrow \\text{after treatment} \\leftarrow \\text{time} \\rightarrow Y$.\n", + "\n", + ":::{note}\n", + "Below is an attempt to explain one way that we can deal with this. Though it is a bit of a brain-twister and can take some time to get your head around. Thanks to Nick Huntington-Klein for some clarification in [this twitter thread](https://twitter.com/inferencelab/status/1783882438063661374).\n", + ":::\n", + "\n", + "One approach we can use is:\n", + "1. We want to close the backdoor path, and one way to do this is to split the dataset into two parts: pre-treatment and post-treatment. By fitting a model only to the pre-treatment data, we have removed any variation in $\\text{after treatment}$ (all values are $0$), so there is now no variation in $Z$ caused by $\\text{time}$. This is one way to close a backdoor path, and means that a model fitted to this data (e.g. $Y_{\\text{pre}} \\sim f(\\text{time}_{\\text{pre}})$) will not be biased by the backdoor path.\n", + "2. However, our goal is to estimate the causal effects of the treatment $Z \\rightarrow Y$, but we have just removed any variation in $Z$ and it does not appear in the aforementioned model, $Y_{\\text{pre}} \\sim f(\\text{time}_{\\text{pre}})$, so our work is not done. One way to deal with this is to use the model to predict what would have happened in the post-treatment era if no treatment had been given. If we make the assumption that nothing would have changed in the absence of treatment, then this will be an unbiased estimate of the counterfactual. By comparing the counterfactual with the observed post-treatment data, we can estimate the treatment effect $Z \\rightarrow Y$. By focussing only on the post-treatment data we are looking at empirical outcomes $Y_\\text{post}$ which are affected by treatment $Z = 1$, but have closed the back door because all $\\text{after treatment} = 1$. The final comparison (subtraction) between the counterfactual estimate and the observed post-treatment data gives us the estimated treatment effect." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Difference in Differences\n", + "\n", + "Difference in Difference studies involve comparing the change in outcomes over time between a treatment and control group. The causal DAG for this is given in Chapter 18 of {cite:t}`huntington2021effect`:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "tags": [ + "remove-input" + ] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAATMAAAEzCAYAAABdWOReAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAB7CAAAewgFu0HU+AAAbmklEQVR4nO3deXCU9R3H8c8mJCEhkN1AMaYkIExBkgCBcBjRhiPIIFBgoAUrCM5EpDPFIgiKyiVQQKpDWyx1OEwLgigKhOEo9yEqhwMaEWQ4AiEcUhIoOSDHfvsH7mM22c1ukt19nv3t5zWzM2Gf3fBNeHjzXHkwiYiAiMjPBek9ABGRJzBmRKQExoyIlMCYEZESGDMiUgJjRkRKYMyISAmMGREpgTEjIiUwZkSkBMaMiJTAmBGREhgzIlICY0ZESmDMiEgJjBkRKYExIyIlMGZEpATGjIiUwJgRkRIYMyJSAmNGREpgzIhICYwZESmBMSMiJTBmRKQExoyIlMCYEZESGDMiUgJjRkRKYMyISAmMGREpgTEjIiUwZkSkBMaMiJTAmBGREhgzIlICY0ZESmDMiEgJjBkRKYExIyIlMGZEpATGjIiUwJgRkRIYMyJSAmNGREpgzIhICYwZESmBMSMiJTBmRKQExoyIlMCYEZESGDMiUgJjRkRKYMyISAmMGREpgTEjIiUwZkSkBMaMiJTAmBGREhgzIlJCwMVs7969aNq0KYYNGwYR0XscIvKQgIvZxo0bkZ+fj02bNiE/P1/vcYjIQxroPYCv/eEPf8CJEyfw5JNPomnTpnqPQ0QeYhJF97XOnj2Ljz76CDNnztR7FCLyAWV3MxcuXAir1ar3GETkI0rGbMeOHcjMzNR7DCLyIeVitmfPHvz2t7/lmUqiAKNUzN59912MHDkSJSUlAB7saprNZpjNZnzwwQewWq3YsWMHhg0bhl/96lfV3n/9+nW88cYbiI6ORk5ODgBg5cqVSEhIQKNGjfDkk0/ixIkT2uu//vprDBo0CE2aNEFMTAxmzJjhdNe2uLgY8+fPR3JyMqKiomA2m9G/f38cOHDA898I8pri4mLMmTMHCQkJaNy4MSwWC37zm9/g2LFj1V578+ZNLFq0CK1bt0ZmZibu3r2L5557Dk2aNMHgwYNx//597bUVFRXIzMzEr3/9a8TGxqJx48bo0qULFi9ebPc6ALBYLDCZTNqj8l7IokWL7Jb16tXL7r3nz5/HtGnT0KxZM+Tk5KCoqAgvv/wyHn74YURGRqJfv34Ovxa/IApKS0sTADJr1iztuaNHj0r//v2lQYMGAkBatmypLbtz545kZGRIWFiYABAAcuHCBRk/frw0atRIYmJitOebN28uN2/elC1btkjDhg2lRYsW0rBhQ235/Pnzq81z5coVSUpKknnz5smdO3ektLRUli9fLmFhYRIUFCQffPCB978pVG8FBQXSqVMniYyMlK1bt0pFRYUcP35cWrVqJSaTScLDwyUqKkosFoskJydLdHS0tl6sWrVKBgwYIJGRkdpze/bsEZEH61///v2lbdu2cvToURERycvLk5EjRwoASUpKkhs3btjNsnXrVu3zVF1/8vPzpVu3bgJA0tLStOfGjBljt45///330qNHD2natKk0b95ce75hw4ayf/9+r38/PS1gYmYzderUajErLy+XvLw8+fvf/679gY4bN07ee+89uXfvnoiIbNq0SUwmkwCQsWPHyoABA+T7778XEZGSkhIZNGiQAJCHHnrI7vezWq2Smpoqr7zySrVZZs6cqa08eXl5nvsGkFc8//zzAkBef/11u+c3btwoAKRBgwZy7tw57flr165JSEiIAJDHH39cPvnkE8nPz5cXX3xRhg4dKoWFhSIiMnr0aAkKCpLs7Gy7z1tRUSE9e/YUANKjRw8pLy+3Wx4bG+swZiIi06dPt4uZzZYtW7R1fNiwYbJy5UqxWq0iIrJ27VoJDQ0VABIfHy/379+v67dKFwEXs3/84x/VYmaTnZ2t/UEfOHCg2vLU1FQBIL179662Yu3fv197b35+vvZ8VlaWAJAzZ85U+3zbt2/X3rNkyZLaf6HkM/fv39e2wDdu3Gi3zGq1SlRUlMMtc9tW/ZgxYxx+3kOHDgkA6d69u8Plhw8fttu6q6xly5ZOYzZr1iyHMTt37pz2+VavXl3tfXPnztWWr1u3zuFMRqXUMTN3NGzY0OmyiIgI7eP4+Phqy1u1aqUtCw4OtlsWGxurfVxYWKh9vHHjRgBAjx49tON3tsfvfvc7hIWFISwsDBcuXKjT10O+kZ+fj3v37gEATCaT3TKTyYRHHnkEAJCbm2u3LDQ0FADQp08fh5935cqVAIBOnTo5XP74449r692///3vug1fSeX19oknnqi2/I9//CPCwsIAAPv376/37+dLAfcTAFVXxMqCgmpuu23FdCQkJET7uKKiQvs4OzsbwIOTCzWFlIytefPmaNSoEYqKinDx4sVqy+Wns+cPPfSQ3fM1rW8AcPDgQQBAVFSU09d069YNOTk5OHnyZC2nrj2z2Yzk5GQcOXKkWpiNLuC2zHzt9u3bAIC8vDx9B6F6CQoKwgsvvAAA+Ne//mV36c+dO3dw7tw5BAcHY+TIkbX6vLb1ori42OlrbFv9Vc9qektcXBwA+N3lTYyZl9l2XV1tslfemiNjWrhwIQYNGoSTJ09i8uTJKCgowI0bN5CRkYHi4mL85S9/Qfv27Wv1ORs1agQADrf2bMLDwwH8HBlvs+2BVN3KNDrGzMvatWsHAPjb3/6G8vJyh6+5fPky5s2b58uxqA7CwsLw1ltvISEhAUePHkWbNm3QsWNHFBcXY+fOnZg0aVKtP2fXrl0BAF999ZXT9ePu3bsAgN69e9s936DBg6NENf3YXk1bV86W3bhxAwCQmprq9L1GpGTMbMe+bAds9TR48GAAwLfffouJEydWW4GsViteeuklpKen6zEe1UJ2djb69OmDzZs34/Dhw8jPz8eNGzewdevWOv/5jR07FgBQUFCALVu2OHzN6dOnAQDjxo2ze75x48YAHhyPrerSpUsAgNLSUqe/d1FRUbXnSktL8fXXX6NRo0YYPny46y/AQJSMme0P+fPPP0dFRQXOnTuHt956C8DPxyZsPyVQWeXnHIXQ9i+noxWk8m5i5c8zatQoJCYmAgD++c9/Ii0tDZ9++im++eYbZGVl4amnnkJZWRl69uxZ66+TfGvChAkoKirCsWPH8O233+LMmTM4e/Yszp8/j7y8PIfHtGzrjLN7540aNUo7q/jGG29UW+9yc3Nx6NAhjB07Fo899pjdskcffRTAg7OctmNvV69exZgxY7QTT+fPn8f9+/cdbr05OvSxcuVK3L59G9OnT/e/W2TpeV2It/z5z3+2u2I/JiZGcnJy5N69e5Kenq4t27dvn/Yeq9WqXZsDQObOnatdTCgikpubK7/85S8FgLRq1UquXLli99558+Zp750zZ45UVFRoy8+ePSvx8fHa8sqPdu3ayY8//uiT7wvVT5cuXRz+GdoeoaGhMnjwYLl8+bKIiOzYsUOCgoIEgLRv314uXbrk8PP++OOPkpiYKACkX79+2oW33333nXTu3FkGDhwod+/erfa+Xbt2aRdyh4SESHx8vISFhUlmZqbduty6dWt55513RETk4sWL2vPR0dHy4Ycfyv3796W8vFzWrl0rERERMnLkSLv1118oGbPi4mJ59tlnJSIiQrp37y7Hjh2Tffv2aT/KVPnRt29fKS8vt/sxD9sjLCxMLl26JOPHj6+2zGQyydixY+XSpUvaVdNV31vZf//7X5k8ebK0atVKQkNDJS4uTiZNmmR3gS0Z27Vr1yQ+Pl46duwoMTExEhERocWq8iMpKUmSk5MdBq/qTw/YFBcXy7x58yQpKUnCw8OlRYsWkpaWJqtXr64xLKtXr5Y2bdpIeHi4PPHEE9o/0LNmzZLExERZtmyZ9pMGIvYx2759uwwcOFCioqLEbDZLt27dZMWKFX4ZMhERZW/OSORpb7/9NgoKCrBgwYJqy8rLy1FQUIADBw5gzJgxOHr0KDp06KDDlDXLycnRLvC9ePGidkGuCgLuolmiuti2bRsWLFiAq1evOlzeoEED/OIXv8CIESOwYMEClxdgk+fxO07kgtVqxYQJExASEuIyUsePH0dpaSkSEhJ8NB3ZMGZELvzvf//D1atXcfPmTaSlpWH79u3Vzjrevn0b77//PoYMGYIVK1a4/DEmvVQ+4+rojL4/Y8yIXDCbzVi6dCnCwsJw5MgRPP3004iMjERcXBzatm2L2NhYREdHY8GCBdi8eTN69Oih98gOVVRUYO3atdqv16xZo9RPnvAEAJGbLl68iGXLlmHnzp3a9VtNmzZFcnIyhg4dinHjxml3nDCaPXv2YMCAASgrK7N7PiQkBAcPHqx2DZs/YsyISAnczSQiJTBmRKQExoyIlMCYEZESGDMiUgJjRkRKYMyISAmMGREpgTEjIiUwZkSkBMaMiJTAmBGREhgzAxs9ejQWLVqk9xj0k8TEROzcuVPvMcgJ3jbboH7/+99j3bp12q9fffVVHachi8WC27dvo3///vjPf/6Dp556Su+RqAreAsiARKTa7ZkXLlzIoOnEFjKbhIQEnDp1Sr+ByCHGzKBu374Ni8Vi9xyD5ntVQ2axWJz+h76kLx4zMyiz2YyCggK751577TUeQ/Mhhsy/MGYGxqDphyHzP4yZwTFovseQ+SfGzA8waL7DkPkvxsxPMGjex5D5N8bMjzBo3sOQ+T/GzM8waJ7HkKmBMfNDDJrnMGTq4EWzbsjJycGGDRsAAKmpqejZs6fOEz3AC2vrx+ghW7NmDa5fvw6TyYQpU6boPY7hMWZumDZtGhYvXqz92kjfMgatboweMgAwmUzax0Za54yKu5l+jructecPIaPaY8wUwKC5jyFTF2OmCAbNNYZMbYyZQhg05xgy9TFmimHQqmPIAgNjpiAG7WcMWeBgzBTFoDFkgYYxU1ggB40hCzyMmeICMWgMWWBizAJAIAWNIQtcjFmACISgMWSBjTELICoHjSEjxizAqBg0howAxiwgqRQ0hoxsGLMApULQGDKqjDELYP4cNIaMqmLMApw/Bo0hI0cYM/KroDFk5AxjRgD8I2gMGdWEMSONkYPGkJErjBnZMWLQGDJyB2NG1RgpaAwZuYsxI4eMEDSGjGqDMSOn9AwaQ0a1xZhRjfQIGkNGdcGYkUu+DBpDRnXFmJFbfBE0hozqgzEjt3kzaAwZ1RdjRrXijaAxZOQJjBnVmieDxpCRpzBmVCeeCBpDRp7EmFGd1SdoDBl5GmNG9VKXoDFk5A2MGdVbbYLGkJG3MGbkEe4EjSEjb2LMyGNqChpDRt7GmJFHOQsaQ0bexpiRxzkKmg1DRt7CmJFXmM1mREVFVXv+1Vdf1WEaCgSMGXmFxWLBnTt3qj2v9y24SV2MGXlc1YP9ZrPZbjmDRt7AmJFHOTprWVBQoPstuEl9jBl5TE2XXxjh/xQgtTFm5BHuXEfGoJE3MWZUb7W5IJZBI29hzKhe6nJlP4NG3sCYUZ3V50eUGDTyNMaM6sQTP2vJoJEnMWZUa578oXEGjTyFMaNa8cbdLxg08gTGjNzmzdv4MGhUX4wZucUX9yNj0Kg+GDNyyZc3VmTQqK4YM6qRHneIZdCoLhgzckrPW10zaFRbjBk5ZIR79jNoVBsN9B7AiPbu3YusrCzt1ytWrLBbPmnSJO3jpKQkZGRk+Go0nzBCyGxsQbNYLNpzr732GgC17lorIpg6dSrKy8sdLq+8zgHAtGnTEBsb64PJ/IhQNVlZWQLArcfUqVP1HtejzGaz3ddnsVj0HklERAoKCqp97xcuXKj3WB7VoUMHt9a54OBgKS4u1ntcw+FupgMpKSluv7Zr165enMS3jLRFVlUg7HK6u94lJiYiPDzcy9P4H8bMgdjYWMTExLj12tqEz8iMHDIb1YPm7rqkyjrnaYyZE+6sMGazGa1bt/bBNN7lDyGzUTlojFn9MGZOuLPCdOnSBSaTyQfTeI8/hcxG1aB16tQJQUGu/0oyZo4xZk64cyzM31cqfwyZjYpBi4iIQGJiYo2vCQ4ORqdOnXw0kX9hzJxwJ1T+fPDfn0Nmo2LQXK13PPjvHGPmhDsnAfx1y0yFkNmoFjRX65S/rnO+wJjVoKYVx18P/qsUMhuVgsaY1R1jVoOaVhx/PPivYshsVAmaq5MAjJlzjFkNajom5m8rlcohs1EhaDWdBODB/5oxZjWoKVj+dPA/EEJmo0LQnK13PPhfM8asBjWdBPCXLbNACpmNvwfN2brlL+ucXhgzFxytQEY7+P/xxx9j79691Z4PxJDZ1CZoCxcuxKVLl3w1mkuMWR3p/ZPuRjdz5sxqdy3o06eP3mNprl27JhaLRcLDw2XPnj3a80a9+4WvubrbxqxZswSApKeni9Vq1XHSnxUVFUlQUFC1ub/88ku9RzM0xswFR7cDMsptf6xWqwwZMkSbyxY0hsyes6DZQmZ7vP/++3qPqql6OyDe9sc1xsyFvLy8an8R1q9fr/dYIiKyZs0al/e+CvSQ2TgKWtVHZGSk5OTk6D2qiIiMGzfObraOHTvqPZLh8ZiZC45OAhjh2MX169cxceLEGl8TSMfIXHF0DK2qwsJCZGRkQER8NJVzVdcxI6xzRseYuaHyimSEg/8iggkTJrj8y7lhwwYfTeQfzGYzpkyZUuNrdu/ejeXLl/toIucYs9pjzNxQeUUywpX/a9euxebNm12+btCgQQ7Pcgaq2bNn45133nH5uilTpuh+drPqTwIwZq4xZm7o1asXunfvjpSUFAwaNEjXWdzZvbQpKSlh0H4ye/ZszJkzx63XGmF3MyIiAsOHD0dKSgq6d+/OK//dwP+dyQ29e/fGkSNH9B7D7d3LysLDw3Hr1i0vTmV85eXlKCoqQlBQEKxWq1vvse1ujh8/3svTOffxxx/r9nv7I5MY4WgnueXDDz/E6NGj3X790KFDsWzZMrf/PwPVffnll3j++efxww8/uPX6yMhIfPfdd2jZsqWXJyNP4G6mn6jN7mV0dDTWrl2Lzz77jCGrJDU1FSdOnMArr7zi1u2pjbC7Se5jzPxAbXYvhw4dilOnTuGZZ57R/USFEYWHh2Px4sX4/PPP0a5dO5evN8rZTXKNu5l+wJ3dy+joaCxduhSjRo1ixNxUUlKCmTNn4t13363xWBp3N/0DY2Zw169fR0JCQo1bZTw2Vj/uHEtLT0/Hzp07+Q+FgXE308Bc7V7y2JhnuHMsjbubxsctMwOrafeSW2PeUdNWGnc3jY1bZgbl7Owlt8a8q6atNJ7dNDbGzICc7V7yTKVv1HTGk7ubxsWYGVBhYSGys7O1X3NrTB/OttIOHz6s41TkDI+ZGdTly5fRu3dvdOzYkcfGDMB2LC0pKQnr1q1DSEiI3iNRFYyZgd29exeRkZHcpTSIkpISNGjQgCEzKMaMiJTAY2ZEpATGjIiUwJgRkRIYMyJSAmNGREpgzIhICYwZESmBMSMiJTBmRKQExoyIlMCYEZESGDMiUgJjVsXbb78Nk8lU58eSJUv0/hLIjxQUFKBx48YwmUyIi4tDeXm5y/dUVFRgwIAB2jq3bt06H0xqfIxZFcePH6/X+zt06OChSSgQWCwWvPDCCwCAK1euYMOGDS7fM2XKFOzYsQMA8Oabb+KZZ57x6oz+grcAquLChQsoLi5267V3797FyJEjkZubCwDo0qULDh06hIiICG+OSIrJzc1FmzZtUFZWhh49euCrr75y+trly5dj/PjxAIDhw4fjk08+4f3ubITqpKSkRHr16iUABIC0b99ebt68qfdY5Keee+45bV364osvHL5m3759EhISIgCkc+fOUlRU5OMpjY27mXVQVlaGESNGYP/+/QCARx55BLt370azZs30HYz81rRp07QtLEfHXc+fP48RI0agrKwMMTExyMrK4h5AFYxZLVmtVowZMwZbt24FAMTGxmL37t2IjY3VeTLyZ4mJiXj66acBAJ9++ikuX76sLbtz5w4GDx6MW7duoWHDhti8eTNatGih16iGxZjV0osvvoj169cDAJo1a4Zdu3ahdevWOk9FKpg2bRqAB2crly5dqn08atQonD59GgCwatUqdO/eXbcZDU3v/Vx/MnnyZO24RpMmTeT48eN6jxRw/vSnP8ncuXNl27ZtcuPGDb3H8bjHHntMAIjZbJbCwkJ56aWXtHVuxowZeo9naIyZm2bPnq2tVBEREXLo0CG9RwpI7dq10/4cAEhcXJwMHTpUmcB99tln2tfWt29f7ePhw4eL1WrVezxD46UZbliyZAlefvllAEBoaCiysrLQv39/nacKTI8++ih++OGHGl8TFxeHlJQUu0fz5s19NGH9iAjat29v9zXykh/3MGYurFq1ChkZGRARBAcHY/369Rg+fLjeYwUsd2LmiD8FbuXKlcjIyAAAPPzwwzh69CgP+LtDz81Co1u/fr0EBQUJADGZTJKZman3SE5Nnz7dbveLD9cP2y7qrl279P7js7Nv3z5txtmzZ+s9jt/g2Uwntm3bhtGjR8NqtQIA/vrXv2Ls2LE6T+XcggUL9B7B7+Tm5mLTpk04ePCg3qPYOXnypPZxcnKybnP4G8bMgQMHDmgXKALA/PnzMXHiRJ2nIm8JCjLWX4NvvvlG+5gxc18DvQcwmuPHj2Pw4MEoKSkB8ODan9dff13nqVzbsWMH9uzZo/cYXrd48eJ6fw6TyYS2bduia9euSElJwcCBAz0wmefYtswsFgtatmyp7zB+hCcAKjl16hTS0tJw69YtAMCECROwbNkynaeiymp7AqBquFJSUtC5c2c0btzYi1PWXVlZGSIjI1FaWopevXph3759eo/kN7hl9pPz58+jX79+WsieffZZvPfeezpPRbXhb+Fy5PTp0ygtLQXAXczaYswA5OXlIT09HdeuXQMADBkyBJmZmYY7lkI/UyFcjvDgf90FfMwKCgqQnp6OnJwcAA92Y2bMmIEzZ864/Tni4uIQFRXlpQmpsjfffBPx8fFKhMsRHvyvu4A/ZvbRRx/V+06dX3zxBVJTUz00EQWyvn37Yu/evQgNDUVhYSFCQkL0HslvBPx+VHZ2dr3eHxwczH9ByWNsW2YJCQkMWS0F/JYZEakh4LfMiEgNjBkRKYExIyIlMGZEpATGjIiUwJgRkRIYMyJSAmNGREpgzIhICYwZESmBMSMiJTBmRKQExoyIlMCYEZESGDMiUgJjRkRKYMyISAmMGREpgTEjIiUwZkSkBMaMiJTAmBGREhgzIlICY0ZESmDMiEgJjBkRKYExIyIlMGZEpATGjIiUwJgRkRIYMyJSAmNGREpgzIhICYwZESmBMSMiJTBmRKQExoyIlMCYEZESGDMiUgJjRkRKYMyISAmMGREpgTEjIiUwZkSkBMaMiJTAmBGREhgzIlICY0ZESmDMiEgJjBkRKYExIyIlMGZEpATGjIiUwJgRkRIYMyJSAmNGREpgzIhICYwZESmBMSMiJTBmRKQExoyIlMCYEZESGDMiUgJjRkRKYMyISAmMGREp4f/x46BQsXZs+AAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pgm = daft.PGM(dpi=DPI, grid_unit=GRID_UNIT, node_ec=NODE_EC)\n", + "\n", + "pgm.add_node(\"z\", \"$Z$\", 0, 0)\n", + "pgm.add_node(\"y\", \"$Y$\", 1, 0)\n", + "pgm.add_node(\"t\", \"time\", 0, 1)\n", + "pgm.add_node(\"g\", \"group\", 1, 1)\n", + "pgm.add_edge(\"t\", \"z\")\n", + "pgm.add_edge(\"t\", \"y\")\n", + "pgm.add_edge(\"g\", \"z\")\n", + "pgm.add_edge(\"g\", \"y\")\n", + "pgm.add_edge(\"z\", \"y\")\n", + "pgm.render();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + ":::{note}\n", + "For our explanation below, we will assume we are dealing with the simplest case of a two-group, two-time period design, the so called \"classical\" 2$\\times$2 difference-in-differences design. \n", + ":::\n", + "\n", + "Our goal is to estimate the causal effect of the treatment on the outcome, $Z \\rightarrow Y$, but now we have _two_ backdoor paths:\n", + "1. $Z \\leftarrow \\text{time} \\rightarrow Y$\n", + "2. $Z \\leftarrow \\text{group} \\rightarrow Y$\n", + "\n", + "From a regression point of view, both $time$ and $group$ are binary variables. In this situation, treatment is given to the treatment group ($\\text{group}=1$) at time $\\text{time}=1$.\n", + "\n", + "The causal effect of the treatment upon the outcome is typically estimated by fitting a regression model of the form `y ~ time + group + time:group`. The interaction term `time:group` captures the causal effect of $Z \\rightarrow Y$. \n", + "\n", + "We can note that this interaction term $\\text{time} \\times \\text{group}$ encodes the values of $Z$, which as we said above, is equal to 1 for only the treatment group at time 1. So another way to think about the inclusion of an interaction effect is that we are simply conditioning on all the observed data ($Z$, $\\text{time}$, $\\text{group}$, $Y$) to estimate the causal effect of $Z \\rightarrow Y$.\n", + "\n", + ":::{warning}\n", + "Achieving an unbiased estimate is strongly dependent upon the {term}`parallel trends assumption`. That is, we assume that the treatment and control groups would have followed the same trajectory (over time) in the absence of treatment. This is a strong assumption and should be carefully considered when interpreting the results of a difference-in-differences study. In the case of the classic 2$\\times$2 design we cannot assess the validity of this assumption empirically, so it is important to consider the plausibility of this assumption in the context of the particular example. \n", + ":::" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Synthetic Control\n", + "\n", + ":::{warning}\n", + "While many texts cover the synthetic control method, they typically do not provide a causal DAG-based treatment. So this section is pending - we hope to update it soon.\n", + ":::" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Regression Discontinuity\n", + " \n", + "The regression discontinuity design is similar to the interrupted time series design, but rather than the the treatment being at a specific point in time, treatment is based on a cutoff value $\\lambda$ along some running variable $RV$. This running variable could be a test score, age, spatial location, etc. The running variable may also influence the outcome $RV \\rightarrow Y$. The running variable may also be associated with a set of variables $\\mathbf{X}$ that influence the outcome, $RV - - - - \\mathbf{X} \\rightarrow Y$." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "tags": [ + "remove-input" + ] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAAFFCAYAAABVK2F7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAB7CAAAewgFu0HU+AAA+BElEQVR4nO3dd3RT5f8H8Hc6aOmggyK0tpSChbLLnjIEFFGGgKCgIBsFEUQQEBkHB4oooIAKSFEQRRyALBkiUBAEGS2bIntji9Dd5PP7g1/ut2mSNk2Tpjd5v87JOem9z715bnLz6Ts397nRiIiAiIiIiEo8N0d3gIiIiIgsw+BGREREpBIMbkREREQqweBGREREpBIMbkREREQqweBGREREpBIMbkREREQqweBGREREpBIMbkREREQqweBGREREpBIMbkREREQqweBGREREpBIMbkREREQqweBGREREpBIMbkREREQqweBGREREpBIMbkREREQqweBGREREpBIMbkREREQqweBGVAzS09PRrFkzhIeH4/jx447ujlNLT0/Hl19+iZo1a2LatGmO7o7Lu379OmJiYlC9enXcuHGj2B535MiRCAgIwMqVKwu9bGJiIsLDw9G8eXOkp6fboXekd/36dUyZMgUPPfQQduzY4ejuqAKDWyFs3boVr776KsLDw6HRaIxunp6eKFOmDKKiotCqVSuMHj0aGzduRE5OjqO7Tg52/Phx/Pnnn7hy5Qo2btzo6O44rTlz5iA6OhrDhg1jQC6Cb775Bt7e3ibrnEajQbly5SxeV3x8PE6dOoWTJ09i9+7dduy1oWXLluG///7Dt99+W+hlN2/ejCtXrmDv3r3cj+xEq9ViyJAhqFKlCmbMmIFbt245ukvqIVRoFy9eFI1GIwDkkUcekXXr1sm5c+fk33//lTNnzsjq1aulZ8+eSpsqVarI5s2bHd1tciCtVisDBw6Udu3aydWrVx3dHaeVnp4u9+7dk6CgIAEgU6dOdXSXVEur1covv/wibm5uAkAAyNixY+X27duFWs/9+/elW7du0q1bN7l//76demts3rx50rhxY/njjz8KvezVq1elffv2MmDAANFqtXboHYk8eL/Gx8cr+9fvv//u6C6pAoOblR566CEBIK1btzbbZvfu3Uo7Nzc3Wbx4sU370KdPH5uuj4qOr0nJ0KhRIwY3G6lZs6YAEB8fH6cKMcnJyTJixAhHd8PlpaamMrgVEr8qtVLp0qULbNOiRQvs3r0bAQEB0Ol0GDZsGPbs2WOTx9+5c6fN1kW2wdek5PD29nZ0F5xGcHAwAKBs2bJwc3Oefxnz5s3D/fv3Hd0Nl8f3auE5z7uwhIqOjsann34K4MF3+sOHD4eIFGmd9+7dw5AhQ4q8HrIdviYlizMFDEfTP5fO9JweOXIE77//vqO7QXCu/aq48BkrBn379kXlypUBAAkJCUU6Of3OnTvo1KkTTp8+bavuURHxNSFSjyNHjqBjx47IyMhwdFeIrMLgVgzc3Nzw7LPPKn+vXr3aYH5mZibmzp2Lhg0bokKFCvDz80PNmjUxZcoUg+Jy9OhRtGzZEgcPHgQAXLx4EYGBgQgMDESHDh2UdlqtFl9//TWaN2+O8PBw+Pj4IDo6GqNHj0ZycrJV2xAfH48uXbqgXLly8PHxQUxMDN59911kZmaaXUZEsGLFCrRu3RoPPfQQfHx80LBhQyxevNio7eXLlzFhwgQEBwfj/PnzEBHMnz8f1apVg6+vL9q3b4+kpCSzj7Vx40Z07NgRoaGh8Pb2Rq1atfDhhx8iOzvboN3FixcxefJkVKhQATt27MD169fRuXNn+Pv7Y/DgwcoRM1u/JklJSZgwYQLKly9vNOQ9NTUVn332GaKjoxEXFwcA2LJlC1q0aAFfX1/Url0bmzZtMrvtFy5cwJAhQxAREQFfX1+Eh4dj2LBhuHbtmtllCpKYmIjnn38eFSpUgI+PD6KiojBu3DjcvXvXZPu//voLAwcOhI+PDwBg3759qFevHgIDA/Hll18q7W7cuIExY8agZs2aCAoKQnBwMFq0aKFsd24pKSmYM2cOqlWrpsxfsGABqlWrhtKlSyM2NhbLly+3aHt2796Ntm3bwtfXF1WrVsV3331XuCeErHbkyBGMGDECAQEBOH/+vMG89PR0fP7556hataryGu/atQtt2rRR6syqVauU9jdu3MCwYcNQvnx5lClTBr179zZZ05KTkzFv3jyjS8KsXLkSHTt2xJ07dwAA3377rfJ+nT59utIuPj4e/fr1M3lKTFFrlTW1tCC//vor2rVrh7JlyyrvjQULFkCn0xm11Wq1WLNmDTp27IjHHnsMABAXF4eHH34YUVFR2Ldvn9L26NGj6NWrF6Kjo+Hn54cKFSqgS5cuJkcHJyUlYfz48QgJCcH58+eRmpqKMWPGIDQ0FH5+fujQoQP++usvi7bnq6++Qq1ateDj44PmzZvj0KFDVj4zTsqB59epWmRkZIGDE3Jbs2aNcgJmzZo1lekZGRnStGlTKVWqlGzatElERM6cOSN169YVANKrVy+jdS1dulQASGRkpNE8nU4n3bt3FwCydOlS0el0cvXqVWnXrp0AkCZNmohOpyvUti5evFjc3Nyke/fucufOHbl//75MmDBBAIinp6f4+/tLQECANGvWTFkmKytLevbsKc8++6xcvHhRRET++OMPqVixogCQgQMHiohISkqKjBw5Unx9fZXn5/Tp09KjRw/x8fGR0NBQZXq1atUkOzvbqH+vvfaatG3bVo4fPy4iIocPH5Y6deoIAHn88cclOztbsrOz5YUXXhAfHx9lfRs2bJBGjRpJ6dKllWlJSUk2fU20Wq0MGjRIypUrZ/IE3Pfff18qVaqkzFu6dKl89NFHUqpUKYmIiFBG9JUqVUrOnj1r9LgHDx6UoKAgqV69uhw/flyys7Pl+++/F29vb3F3dxdfX18JCAiQkJAQuXHjhkWv94YNG8Tb21tatmwply5dkoyMDPn0008FgHh4eIifn58EBARIRESEnDx5Up544gllBDUAOX78uISHhyvTIiIilOewXLlyUrlyZTlz5ozodDrZtm2bBAQECABZsGCB0odJkyZJhQoVDJ6XUaNGiY+Pj7IP6W8zZ8402obWrVsrgxNWrlwpXl5eUrFiRfHy8lIGC+3Zs8ei58PV6Z9LU/UmPykpKfLcc8+Jv7+/8lr9888/yvz333/f4H2xdOlSWbZsmXh7exvsP25ubvL7779LUlKSRERESNmyZZVRwwCkQ4cOBo87e/ZsqV69ujLf1ACV/v37CwDp37+/wfRz585J165dxdvbW1k+9/YUtVZZU0sLMnXqVAEgI0aMkPv370tycrK8+OKLSt0oU6aMBAYGSu/evWXNmjUSGxur9LN169ayfPlyCQwMVKa9+OKLIiLy008/ibu7u3Tu3FlSUlIkMzNT5s6dq/T1yJEjIiLy77//yosvvqi8t/Q1oEmTJlK2bFllgB4A8fb2lh07dhhtQ+7a+Oqrr4q3t7dERkYq+0C5cuXk3r17Fj8nzo7BzUqFDW4HDhxQdk5/f39luv4ffvPmzQ3a//zzz0rRSktLM5iXX3D7/fffBYCEhYUZTD906JDy+CdPnrRsI+XBP1tPT0/x8fGR5ORkZbpOp1MKwIABA4yWe/PNN6Vp06aSk5NjMH379u1KP3788UfRarWSkZEhv/zyizK9c+fO8umnn0pGRoaIiHz//ffKPH2Q0lu4cKFERUUZvamTkpKUN/3s2bOV6QcPHlTW1aZNG9m9e7dcunRJnnvuOXnppZdEq9Xa/DURETl27JjJ4JaWliZ37twRd3d3ASCPPfaYjBgxQm7evCkiIqdPn5YyZcoIAJkwYYLBOrOzsyU6OloAyG+//WYw77XXXlMuV5P3NchPcnKyBAcHK/+UcuvatasAkLZt2xotN3v2bGX7+vbtKzdv3pTdu3dL+/btZc6cOSLyv3+WkyZNMtnXpk2bGky/ffu2lCpVSgDI008/La+88oryvB86dEgiIiKUMJmQkGCwrD5stGzZUnr16iWXLl0SkQeXedAHZf0/KMqftcFNb/369SaD2+XLl+Xw4cPKvG7dusnrr7+uXG7k5MmTShhq166ddOjQQTZu3Cg6nU50Op28++67yrInTpwweMycnBwl+BcmuOnNnz/fKLgVtVZZW0vz88cffwgAqVixokFQTEtLU7Z/+vTpRst17txZAEj16tWlf//+kpOTI0uWLJFHH31UqU/6/3F5a4v+A2zeerRu3Tpl25955hlZsmSJcpDg22+/Vd7LFStWlMzMTINl9cu1a9dOpk6dqlw2ZteuXUptXLJkSaGeG2fG4Galwga3M2fOKDunh4eHMl1/JKNbt24G7U+cOKG0v3z5ssG8/EKC/s0TGxtrMD09PV1Z3+7duy3bSBGZPn26AJC6desazfv4448FgHh5eRlcJuDmzZvi7e0tn3/+udEyufuRe5uPHz+uTDd1zbuwsDABIJ988okyLTMzUypUqGBUQPTKly9v9FxkZGQoj/P222+bXM7Wr4nIg0JqKrjp6T+VTpw40Wje888/LwCka9euBtP1RRuAwT8CEZG///5bmRcfH2+yT6YsW7ZMAEhAQIDRvJ9++klZ55UrVwzmbdq0SZm3bds2k+vu0aOHAFCCnN7ChQuVkJlXeHi4AJBBgwYZzdu5c6fymMOGDTOYpw8bvXr1MjrCPGXKFAEgDRo0MNlPMlTU4Jb7vZ07uOmVLVvWbMCYOHGiAJDg4GC5deuWwbzMzEzlKM9PP/1ktGzjxo2tDm4bNmwwCm6mtsfSWiViXS0tyIABA0zWBhGRUaNGCQCJjo42mqc/yhcUFGT0vOrpX5fDhw8bTO/du7cAkMGDBxtMP3v2rPK8fPPNN0brmzFjhjJ/5cqVBvP00+fPn2+0XKtWrQR4cA1BeoDnuBWTe/fuKfeDgoKU+4MGDcKXX36JOXPmKNN0Op3Bd/qF+eWFJ598EnFxcfjmm28MpuvPwSrs+q5evQoA0Gg0RvP0Ay4yMzMNrnq9ceNGZGRkYOzYscr5I/pbhQoV4OXlBS8vL4PzXby8vJT7VatWNXqs8PBwAIbP459//onr169j7ty5Ro8TGBiIu3fvwsvLCxcuXFCWKVWqlHJff35HXrZ+TYCCh7zrt9/UtkdERAAw3Hbgf68NYPz66F8bALh06ZLF/bTk9Ta1Tkue148//hhffvklhgwZokxLTU3FyZMnAZh+Tt3d3QEALVu2NJr36KOPokmTJgBg9qdyqlevbrQtYWFhAID//vvP5DJkWwXt+/rzIitWrGg0r1KlSgAAf39/hISEGMwrVaoUypYtCwAmL+tRlMtM5LesNbUKsK6WFsSSdZp6/+vfr3Xq1DF6XvV+/fVXrFq1CnXr1jV4PP25s3nfr/r3KmD6/Tpy5EjluTP3fq1Ro4bRNL5fjXk4ugOu4vr168r9yMhI5X7p0qWVf2Q3b97EwoULsWHDBoM3ixTiEhPu7u7o378/gAc7+pIlS7Bq1So0atTIqvVFRUUBgHISbu4CoV+Pl5cXAgMDlekJCQkAHgzC6Nixo0WPU9CQcP1Jwlqt1uhxPv74YwwfPtyixzFV4Ew9li1fE0seN7/tN7XtwP9eGwD4559/EBsba7J/5cuXt7if+nWmpKQgJSXF4HXNb52WPK8VK1ZUntfExER8+umnOHbsmPLzSYV9TgGgTZs22LdvX6HCqf6fR96BK2QfRdn3c38gMMXT0xOA8XvDksfNT37LWlOrAOtqaUH06/znn3+M5unXaer9b8lz07RpUzRt2hQAsGnTJixcuBBarRZpaWkG67dUYGAgYmNj+X61AR5xKyZ79+5V7rdq1cpgXnp6OsaPH4+WLVuiatWqiI+Px6RJk6x+LK1Wiw8++ACxsbHw8vLC9u3bMW/ePKvW9eKLL8LX1xcpKSlYs2aNwbwDBw4AALp3727wKTQlJQUAcOXKFes2wEL2fBxbvyb20KRJE9SrVw8AjEZl6l+byMhItGjRwuJ1dunSBaGhofmus1mzZsqRkMK6dOkSnnnmGQwfPhz9+vXD7t270bVrV6vWBfzvaKQ1oc+aZYisZU0tLciQIUPg7u6OI0eOGI281K+zT58+Vvd5165dqF+/Pr777jvMmTMHv/76K2rWrGn1+vh+tQ0Gt2Ly888/K/eff/555f6tW7fQqFEjbNiwAX/++Seef/55eHhYfyA0PT0d7dq1w9y5c7F9+3a88sorFv3KgzlhYWFYtWoVAgICMGrUKMTHxyMrKwvr16/HvHnzULVqVcydO9dgGf1XH+YOh+uZ+pRcGPZ6HFu/Jva0evVqVK1aFZ9//jm+/vprZGRkICEhAaNHj4afnx9WrFihHJWwROnSpbFmzRqEhoZi+vTpWL9+PbKyshAfH49p06ahfPnyJi/dYYn4+HjUqVMHvr6+2LlzZ6ECpTn6IzKFOapI5AjW1NKC1K9fH1988QW8vLzQv39/JCYmIjMzE19//TW+++47NG/eHJMnT7aqv59++inatGmD4cOHIy4uzuAIv7X4frUNBrdisG7dOiQmJgIAnnjiCTRs2FCZN3ToUBw7dgxvvfWW8tMyRTFp0iT88ccfePXVV60+KpJXx44d0blzZ1SvXh19+vRBuXLlMH78eIwcORL79+9XvurSq1atGoAHoSL3+WW56X8CrCj0j7N7926Daw/llfucKkvY+jWxp8qVK2Ps2LGoU6cOPvzwQ5QvXx5PP/00GjVqhIMHD1oVjho1aoQBAwagTp06GD16NMqVK4eXXnoJzzzzDP7++2+T5/UUJDMzEz179kRKSgpmzZpV6Kulm/u0fePGDQAPjgJSybB//3789ttvju5GiVTYWmqJvn37olWrVoiIiMDjjz+OChUqYM6cOXjnnXewfft25QNuYRw6dAijR4/GI488gqFDhxZ6eb5f7atkHkZwIvoLRgIwuhgpAKXA5T6xM6+8F1HU/9MzdeVva9ZXkEGDBsHNzQ2bN2+2qH2nTp3g7u6OjIwM9OrVC7/99hsCAgIM2nzyySdF/gTXqlUrlClTBv/99x/69u2LHTt2KCcG6/3444/5Phem2Po1sae4uDjMnj0bR48eLdRXLPmZPn064uPjsWPHjiKdJ5Tb0aNHlfM8zT2v+e2XqampJqfrfxtWf14nOd4nn3yCJUuWOLob+XLU+7WwtbQgOTk56NKlC5o1a2ZwAeGi2rp1K3Q6ndX/R0y9X7OysnDw4EH4+vqiR48eNumnq+IRNytZ8oZPSkpCq1atcO3aNZQrVw4bN240GjmlH2H64Ycf4t9//wXw4KrzL730ktLm2rVr2Ldvn/JPyt/fH8CDE+dPnjyJ7OxsvPzyy8jMzFTWt2DBAly8eBEAcOLECXTv3t1gfWfOnMGvv/5a4DZs2bIFcXFxSE5Oxp49e3D8+HGcOnUKZ86cwfnz55UrkOdWqVIlDBw4EMCDT9+xsbH44osvcOjQIezYsQMvv/wyPv74Y4wcOVJZJveIofyuIJ6enq7c9/Pzw5tvvgngwXNdr149zJ49GwcOHMDu3bsxceJEDB06FG+99ZayTO6RUPrnOy9bvyYAlBN6825D3u23dNsB4O7duxgxYgR0Oh127tyJxMREnDp1CqdPn8a5c+dw/fr1Qn9NfPr0aUyfPh2ZmZnYuXMnjh07pqzzn3/+wa1bt0x+mi7oec09kvqtt96CVquFVqvFqlWrMGPGDGW5zMxMLFq0yGiUoKmvwxMTE/Hbb7+hQ4cOeOKJJwzm6f9xZGVlmd1WU68DGdOfFG5Jzfvmm2+QlZVlcJSnoH1fP83U+vX7lbnXUb9/m1qv/nFNzdO/X/fv34/09HTcvHkTY8eOtajP1tQqwLpaWpClS5diy5YtuHLlCvbv348TJ04o67xw4YJyHnBe+ue1oBp44sQJLFu2DMCDkbszZ85UTvu5du0asrOz8dlnnxktb+r9umTJEqSkpGDixInKaGDAMOTx/WqhYr78iFO4ePGicnHXWrVqyaFDh5QLMGZmZsrBgwfljTfeEG9vb9FoNNKtWzdJSkoyua7XX39duYaNt7e3hIaGSlRUlOzYsUOZHhAQIO3atZOsrCzl8T08PASAlC5dWsqWLSuffvqpiIjMmzdPWc7T01PCw8MlJCREtm3bplwd28/PT+rVqycpKSkFbuvatWuV9Zm7RUZGysKFCw2WS01NVX6tIe/Nz8/P4NpiOp1O3n77bWX+O++8Y3DtrbNnzyp9r1u3rty9e1eZp9Vq5YUXXjD5OJ6enrJ69WqDx1m0aJEyv1WrVsrFPu35moiILFmyRFl20KBBBhfF3bp1qzKvU6dOyr4kInL//n1p0qSJABAfHx/l1yFERG7cuGHwqw+mboGBgTJq1CijCwabc/jwYeXXGszdypcvL9OnT1e2ISMjQ7meFAB5/fXXjS6wqdPppGHDhkqboKAgCQwMlPbt2xs8NyEhIQYX6NVfL9HLy0smTpyoXK9u37598sgjj0jNmjWNfhHi4MGDyvPSsGFDuXPnjjJP/4se+sfbuXOnRc+LK9NfSw+A7Nu3z2i+TqeT8+fPy+TJk8Xd3V2WL19uME9/3TwAMm3aNIP39rZt25R5Tz75pMG+n5aWJh07dhTgwUWv875W27dvVy7O+uSTT0pqaqoyLzExUdkH6tSpY1AzRB5cEFb/uMHBwRIUFCR//fWXiDy4eO9LL72kzF+6dKnB9lhbq6ytpfnJXe/N3WJiYuTHH39Ulrl9+7byyzJubm6yatUqo2sdXrlyRfz8/JR1hIWFSenSpWXs2LHy6quvGtSC9evXi4jIP//8Y/CcrlixQjIzMyUnJ0e+/fZb8fHxkd69extcp06n0ynXcQQeXI8xd+24ffu2REVFCQB5+OGH5dq1axY/N86Mwa0QduzYIWPHjjX6yR39Tf/zQlFRUdKuXTuZMWOGHD16NN91pqWlyYgRIyQkJESCg4Nl6NChyj+nl156Sfz8/GTAgAHKlaT14uLiJDQ0VEJDQ+Wjjz5Spmu1WpkyZYqEhoaKv7+/wVXj3377bfH19ZWuXbta/PNHIiLvvPOOhIWFSeXKlSUwMFC5Anbe23fffWewXHZ2tsydO1fq1q0rXl5eEhwcLD179jQIHzk5OQY/laK/lSpVSk6fPq1c4Tv3zd3d3eAq2jqdTr7++mtp2rSp+Pj4iL+/v3Ts2FH27t1r0J/cP5Ojv2k0Gvniiy/s+po88sgjJkPluXPnTIZb/fbNmjVL+ceU+5b7oqF//PGHBAYGSkxMjISEhIiXl5fJ4DV8+HCLX+/ly5dLSEiIREdHS3BwsMnXB3jwU1N///23yT56eHgYvM4iIhcuXJBOnTqJn5+fREZGygcffCBarVbu3bsnsbGxEhISIrNmzTJYRh/c5syZI8OGDZPy5ctLmTJlJDo6WiZPniz//fefQXv9T/3kfY2nTJkiK1asUMJ17pupi5e6uitXrsgvv/yiXGw196106dISEBAgAQEBUqZMGYPXv1SpUsoHQnPvbS8vL8nJyTG573t4eMiKFStk+fLlJvfjSpUqiYhITEyMyffNrl27ZNCgQUbz3NzcZNmyZcr2abVaGT16tPj7+0vNmjWVXzk4d+6cyfpWpUoVm9Qqa2upOTqdToYMGSKVKlWSyMhICQgIEE9PT5Pvgb179yo/WZX3lveXdkQe/L+rV6+eeHt7S2xsrKxdu1ZEHvwKUPny5aVatWqyceNGpX3u4LZx40Z56qmnJCAgQAIDA6VRo0ayePFio4sL63/9JO8+sGXLFhk+fLjBT+nptyP3L+G4KgY3yte1a9ekUaNGRlfKF3lQNFJTU+XEiRPSpUsX6d69uwN66LpycnKkS5cuBsUzt4yMDLl48aJMnTpVgoODLVrn/fv3pWXLlsrvEOam0+kkPT1dkpKSZOjQoVK/fv0i9d8S+uCW+6gHkRrZo5aeOHFCmjZtavQBRkSUD0WHDx+Wpk2byuuvv17kbchP7uBm6hcyyHZ4jhuZpdVq0aNHDzzzzDPK1atz02g08PHxQUxMDAYPHlzokYJUNJMmTUJGRobZixx7eXkhIiICY8eOtfi1GThwIGrUqIE6deoYzdNoNPD29kblypUxcuRIvt5EFrJHLb137x46d+6MUaNGKefs5ebm5gY/Pz/UrVsXffr04fvVifCVJLOWLl2KPXv2FDhaUUTw1VdfoWfPnsXUMzp9+jQ++ugji0aSLl682KLXZuvWrVi1apVN10lE9qmlH3zwAc6ePVvgOrOysrB8+XKO5HQiDG5k1okTJwAAEydOxPTp041+VkVE8Pfff6NLly5wd3dH7969HdFNl3Tq1CnodDqsW7cOL7zwAg4ePGg0PP/ChQt488038cUXX+Ddd98tcJ361/uzzz7D6NGjcfz4caMRpCdPnsTAgQNx4MABjB492mbbY45+1B5HlJGa2aOW6tc5cOBAzJkzR/kNUT2tVotdu3ahbdu2aNmypfLzVfaSe4Qt36925tAvaqlEu3HjhjRu3Njg5NCAgACpUqWKREVFia+vr2g0Ghk1apRkZ2c7ursuJScnx2AkJ/BgBGylSpXkkUcekeDgYAEgbdu2tXggSmpqqnTq1MlowE1UVJRUqVJFypQpIwCkd+/eRgMz7GHz5s3Kycnt2rWTe/fu2f0xiezBHrX01KlTUrVqVYN1BgcHS5UqVaRSpUri5eUlnp6e8t5779l56x7Uo9yjhydNmmQwcp5sSyPCHwAj83Q6HVatWoVvv/0WBw4cwO3bt+Ht7Y2IiAi0bdsWw4YNQ+3atR3dTZe1Z88efPnll9i9ezcuX74MNzc3PPTQQ2jevDn69u2Lp556qtDr3LhxI+Li4vDnn3/i+vXr8PDwQFhYGB599FEMGjTIJj9VlR/99Qjzfmp3c3PDzJkzMW7cOLs+PpE92KOWZmVlIS4uDj/88AOOHDmC5ORk+Pn5ISoqCu3bt8crr7xis1/QMWfbtm148sknjX4E3tPTEzt37rT7kT5XxOBGREREpBI8x42IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIiIiFSCwY2IiIhIJRjciIioxBERlClTBm5ubihfvjx69eqFCxcuOLpbRA7nUsFt4cKF0Gg0Jm++vr6IiYnB8OHDceLECbPraNu2LTQaDYKCggr9+P369VMeb/Xq1UXZFFK5jIwMxMTEGO2Hr776ar7L7dmzB+7u7kbLbd26tZh6TrbCepS/pKQk3Lt3DyKCmzdv4ocffkCnTp0c3S0qJNY6OxAXMnToUAFQ4M3Ly0u+++47k+sYOXKk0u7KlSsWP/bhw4fFzc1NAEjz5s1ttUmkYvHx8co+ob9pNBrZuXOnyfYZGRkSExNjtL8OHjy4mHtOtsB6lL/79+9LQkKCbNq0SaKiopTtPHDggKO7RoXEWmdbLhXcGjduLAAkICBAEhISlNvBgwdl1apV0rx5c2UHKV26tFy+fNloHQsXLlTabN261eLHfvzxx5Xl9u7da8vNIhUbPXq0UXGKjo6WtLQ0o7YTJ040ahseHi4pKSkO6DkVFeuR5ZYvX670d/HixY7uDlmBtc52XCa4abVa8fHxEQDSsmVLs22aNm2q7CizZs0yarNz505l/rx58yx67C1btijL9O7du0jbQc4lNTVVqlSpYlSk3njjDYN2Bw8eFA8PD6N269evd1DPqShYjwrn3LlzSp9Hjx7t6O6QFVjrbMdlznE7deoU0tLSAAB16tQx2cbNzQ0vv/yy8vexY8eM2tSsWVO5f/z48QIfV0Qwfvx4AICXlxfef//9QvWbnJuPjw8WL14MjUZjMP2TTz7B/v37AQA5OTkYOHAgcnJyDNr069eP5/yoFOtR4VSqVAn+/v4AgMTERAf3hqzBWmc7LhPcDh8+rNyvXbu22XaRkZHK/bw7DwAEBwcjNDQUAPI9aVhvxYoVOHToEABg1KhRiIqKsrTL5CLatGmD4cOHG0zTarUYOHAgsrKyMHPmTBw5csRgfoUKFTBnzpxi7CXZEutR4Wg0GkRHRwNgcFMz1jrbcMngZu4TLgDcuHFDuW+uqOk/5Rb0CTczMxNvv/02AKBs2bKYNGmSpd0lF/Phhx8a/JMGHhxhGThwIGbMmGHUfsGCBVaNJKSSwZnrkVarxcSJE3Hz5k2brfPgwYPKc3b9+nXcuXPHZuum4sVaV3QuGdxq1apltt0vv/yi3O/atavJNvpCeevWLdy+fdvsuj777DOcP38eADB16lQEBgZa3F9yLX5+fli0aJHR9BUrViArK8tgWu/evfHMM88UV9fIDpy5Hr300kuYOXMm2rRpYxA8raXVajF06FDodDplmqmvjUkdWOtswNEn2RWX8uXLCwCpVKmS2TY///yzMmS5Z8+eZtt9+eWXygmT5oYzJycnS3BwsACQqlWrSlZWVpG3gZzfoEGD8r00RLly5eTmzZuO7iYVkTPXo+XLlyv9rl69uly7dq1I6/vkk0+M3gfz58+3UW/JUVjrrOcSwe3atWvKztC5c2eDeRkZGXL06FEZM2aMuLu7K6O8/vvvP7Pr27Nnj7K+zz//3GSbcePGKW1+/vlnW24OObGUlBR5+OGHzRYzc9fzIvVwhXqUO7zFxMRYHd4uXbokfn5+gv+/3px+G15++WUb95iKG2ud9VwiuG3YsCHfZK+/NWjQQObPny/Z2dn5ru/u3bvKMq+99prR/IsXL4q3t7cAkFatWtlpqwzdvHlT7ty5Y/aTdFZWlty5c0du3rwpN2/elNu3b0tmZqbJttnZ2YVq+++//yptb926JRkZGUVum5OTI8nJyUrbmzdvSnp6epHbarVaSUlJMWhb0rz44osm909fX1/5999/Hd09KiJXqEciD8KbPnxWq1ZNrl69Wuh1dOvWTQCIv7+/XLp0SQIDAwWAPProo3boMRU31jrreBT8Zar65T6fJD/379/Hk08+CQ+P/J+WMmXKIDw8HJcvXzZ5QvDkyZORkZEBjUaD2bNnW9PlQluwYAEAoH379mjRooXR/EOHDmHjxo0G01q3bo02bdoYtU1MTMSaNWsMpjVr1gyPP/64UdtTp04Z/VxOw4YN8dRTTxm1PXfuHFauXGkwrW7duujWrZtR20uXLmHZsmUG02rUqIFnn33WqO3169exePFig2nR0dHo06ePUdvbt29j4cKFBtOmTp1q1M5Rdu7cieXLl5ucl5qaijFjxiAuLq54O0U2pZZ6FBcXhwEDBljcPj+nTp3C448/joSEBIuXWbt2rXKO33vvvYfw8HDUrl0bu3bt4shSJ8BaZz2XGJyQu1Du2LEDCQkJSEhIwL59+/DNN9+gXr16AB4Ul/79+1u0TnMjuY4eParsjH369EHDhg1tsAXkCtLT0zFo0CCIiNk2y5Ytw6ZNm4qxV2RrrlqPkpOTLW6bmpqq/JZlkyZN8MorrwD436VTkpOTcfXqVdt3kooFa13RuNQRt5CQELRu3dpgXuPGjdGzZ080bNgQx44dw65du3Dw4EE0aNAg33XWqlULmzdvxpUrV/Dff/+hTJkyAIDx48dDp9PB29sb7733nl22h5zT5MmTcfbsWYNpHh4eRtfvGjp0KI4dO6ZckJTURS316JlnnkHTpk0LtUxu8fHxGDZsGLRaLUJCQrBhwwaLl50yZQouXrwIT09PLFq0CG5uD44x5L7mXWJiIsLCwqzuHzkOa13ROP0Rt7S0NGUH0X+Szcvb2xuTJ09W/jZ3+Da33Fcs11/4cvv27di8eTMAYMyYMahYsaLV/SbXsm/fPqOLTGo0Gqxbtw7Vq1c3mH7p0iWMGzeuGHtHtqKmehQQEICYmBirbnfu3MHo0aOV0LZt27Z8r1eX25EjRzB37lwAwBtvvGEQ1nKvg1+XqhNrnQ04+iQ7e8s94mrcuHFm26WlpSmjlypWrFjgevft26es96uvvhKdTif169dXhjHfvXvXlptRIA5OUO/ghMzMTKlRo4bRCbrDhw8XEZG9e/cqI/T0N41GI9u3b3dwz6mwXKEexcfHi7+/v/LYR48etXhZrVYrjRs3FgBSpUoVox8gzz0QY8CAAbbuOtkZa51tOH1wW7BggbIDfPvtt/m21Y9gAiBHjhzJt+29e/dEo9EoBXjFihXKsgsWLLDlJpCTmzRpklEhi4iIMLgExNixY43aVK5cWVJTUx3YcyosV6hHTzzxhBLaEhISCrXsZ599pvR7y5YtJttERkYKAGnUqJEtukvFiLXONpw+uA0dOlR58Y8fP55v20WLFilt33nnnQLXXalSJQEgHTp0kKioKAEeXHCyoOH7RHp///23eHh4GBWqDRs2GLRLS0uT6Ohoo3amLv9AJZcr1KO7d+9K165dJTExsVDLXb16VQICAgSA9OvXz2y7p59+WrlkhE6nK2p3qZiw1tmO0wc3/WF3Hx8fycnJybftlStXlE+tTZs2LXDdTz31lNHOtW7dOlt1nZxcdna2xMbGGu1D5v5p7dq1S9k/9Tc3NzeJj48v5p6TtViPzHv22WcFgISEhMitW7fMtst91CYpKakYe0jWYq2zLacenKDT6ZQTWGvXrg13d/d824eFhSknDO/fv7/AH0nOfUIwADz22GN4+umni9BjciXvv/++0TW9ypcvj08++cRk+5YtWyqXSNDT6XQYNGgQMjIy7NVNshHWI/M2btyIH374AQAwe/ZshISEmG2bd2QplXysdbbl1MHt9OnTSEtLAwDExsZatIz+wrE6nQ7r16/Pt23uQunm5oaPPvrIuo6Syzl+/Djeeecdo+nz589HcHCw2eXee+89VK5c2WDayZMnMW3aNFt3kWyM9ci09PR0jBgxAgDQrl079OvXL9/2DG7qwlpnexqRfK6AR0REREQlhlMfcSMiIiJyJgxuRERERCrB4EZERESkEgxuRERERCrB4EZERESkEgxuRERERCrB4EZERESkEgxuRERERCrB4EZERESkEgxuRERERCrB4EZERESkEgxuRERERCrB4EZERESkEgxuRERERCrB4EZERESkEgxuRERERCrB4EZERESkEgxuRERERCrB4FaCvfDCC/jggw8c3Q0iIiIqITwc3QEyrU+fPli5cqXy95tvvunA3hAREVFJoBERcXQnyJCIwM3N8GDozJkzGd6IiIhcHL8qLYE0Gg2Sk5MNpk2YMIFfmxIREbk4BrcSKjAwkOGNiIiIDDC4lWAMb0RERJQbg1sJx/BGREREegxuKsDwRkRERACDm2owvBERERGDm4owvBEREbk2BjeVYXgjIiJyXQxuKsTwRkRE5Jr4ywkWOH/+PFavXg0AaNasGVq0aOHgHj2QkpKCoKAgg2n8hQXnsXz5cly/fh0ajQZjx451dHfIBZw9exa//PILAODRRx9FkyZNHNshcglLly7FnTt34Onpiddee83R3SnxGNwsMH78eMyaNUv5uyQ9ZQxvzkuj0Sj3S9I+R85r4MCBWLp0KQAgICAAKSkpju0QuQTWusLhV6Uqx69NiYiIXAeDmxNgeCMiInINDG5OguGNiIjI+TG4ORGGNyIiIufG4OZkGN6IiIicF4ObE2J4IyIick4Mbk6K4Y2IiMj5MLg5MYY3IiIi58Lg5uQY3oiIiJwHg5sLYHgjIiJyDgxuLoLhjYiISP0Y3FwIwxsREZG6Mbi5GIY3IiIi9WJwc0EMb0REROrE4OaiGN6IiIjUh8HNhTG8ERERqQuDm4tjeCMiIlIPBjdieCMiIlIJBjcCwPBGRESkBgxupGB4IyIiKtkY3MgAwxsREVHJxeBGRhjeiIiISiYGNzKJ4Y2IiKjkYXAjsxjeiIiIShYGN8oXwxsREVHJweBGBWJ4IyIiKhkY3MgiDG9ERESOx+BGFmN4IyIiciwGNyoUhjciIiLHYXCjQmN4IyIicgwGN7IKwxsREVHxY3AjqzG8ERERFS8GNyoShjciIqLiw+BGRcbwRkREVDwY3MgmGN6IiIjsj8GNbIbhjYiIyL4Y3MimGN6IiIjsh8GNbI7hjYiIyD4Y3MguGN6IiIhsj8GN7IbhjYiIyLYY3MiuGN6IiIhsh8GN7I7hjYiIyDYY3KhYMLwREREVHYMbFRuGNyIioqJhcKNixfBGRERkPQY3KnYMb0RERNZhcCOHYHgjIiIqPAY3chiGNyIiosJhcCOHYngjIiKyHIMbORzDGxERkWUY3KhEYHgjIiIqGIMblRgMb0RERPljcKMSheGNiIjIPAY3KnEY3oiIiExjcKMSieGNiIjIGIMblVgMb0RERIY8HN2Bkmj79u1Yu3at8vfixYsN5o8ePVq5X6tWLQwePLi4uuZy9OEtKChImTZhwgQAwJtvvumobtmciGDcuHHIyckxOT/3PgcA48ePR1hYWDH0jJzZpk2bsGnTJuXvpUuXKvfv3r1rsN/Vq1cP/fv3L87ukRPKycnBuHHjICIm5+etdZMnT0ZISEgx9Ew9NGLu2XNh69atQ5cuXSxqO27cOHz44Yd27hGlpKQYhDcAmDlzplOFtzp16iAhIaHAdu7u7rh37x5Kly5dDL0iZ/b999/jueees6jt1KlTMW3aNPt2iFxClSpVcO7cuQLbeXt74969e/Dw4DGm3PhVqQkNGjSwuG3Dhg3t2BPSc4WvTS3d72rWrMnQRjbBWkeOYOm+FBsby9BmAoObCWFhYahQoYJFbQtT+KhonD28WbovcZ8jW6lSpQoCAgIsasv9jmyFta5oGNzMsGSHCQwMROXKlYuhN6TnzOGNxYyKm0ajQf369QtsFxoaitDQ0GLoEbkC1rqiYXAzw5Idpn79+tBoNMXQG8rNWcNb3bp14eZW8FuSxYxsyZL9ifsc2ZIlHxYA7nfmMLiZYcl38NypHMcZw5uPjw9q1qyZbxt3d3fUrVu3mHpEroC1jopbUFBQgd9WeXt7o0aNGsXUI3VhcDPDkkLFk3UdyxnDW0H7HQcmkK2x1pEjFLTfcWCCeQxuZlgyQIGfQh3P2cJbQfsU9zmyNUsGKHC/I1tjrbMeg1s+8ttxODCh5HCm8MZiRsWtoAEKHJhA9lDQUVzWOvMY3PKR347DgQkli7OEt4IGKLCYkT3kt19xnyN7KGiAAvc78xjc8pHfJwLuVCWPM4S3/AYocGAC2QtrHRW3/AYocGBC/hjc8pFfweLJuiWTM4Q3c/sdByaQvbDWkSOY2+84MCF/DG75yG+AAj+FllxqD2/m9i3uc2Qv+Q1Q4H5H9sJaZx0GtwKY2oE4MKHkK0x4mzlzJi5cuFBcXSsQixkVN3MDFDgwgeyJtc46DG4FMLUDcWCCOlgS3qZNm4aJEydi8ODBEJHi7qJJ5gYosJiRPZnav7jPkT0xuFmHXyIXwNT5Hdyp1EMf3oKCgpRpEyZMAACkp6dj+vTpAICtW7di0aJFGDp0qEP6mZt+gEJCQoIyjQMTyN5Y66i46QconDt3TpnGgQkF4xG3ApgqXDxZV13MHXnThza9sWPHlpivTPPudxyYQPbGWkeOkHe/48CEgjG4FcDUAAV+ClUfU+Etr/v375eYr0zz7mPc58jeTA1Q4H5H9sZaV3gMbhbIvSNxYIJ6BQYGYuzYsfm20X9l6mgsZlTc8g5Q4MAEKg6sdYXH4GaB3DsSByao17Rp0zB79uwC25WEr0zzDlBgMaPikHs/4z5HxSHvaGbudwXjF8kWaNOmDTZt2gStVounn37a0d0hK0ybNs3onDZz9F+Z/vbbbw4L6T4+PujRowfOnTvHgQlUbB577DHs3LkTWq0WnTp1cnR3yAUEBwejc+fOuHr1KkqVKsWBCRbQSEk4oYfIjnJycjBx4kR8/PHH0Ol0Fi/3xRdflIhRpkRERHoMbuQy9u7diwEDBuDUqVMWtffz80NiYiIiIyPt3DMiIiLL8Bw3chnNmjXDoUOH8MYbb5i8wG1eJWmUKREREcDgRi6mdOnSmDVrFnbv3o1q1aoV2L6kjDIlIiIC+FUpubD09HRMmTKlwHPf+JUpERGVFAxu5PIsOfetffv2Dh1lSkREBPCrUiKLzn3jV6ZERFQS8IgbUS75HX3jV6ZERORoPOJGlEt+R984ypSIiByNwY0oj/xGnvIrUyIiciQGNyIzzB19i4+Pd2CviIjIlfEcNyIL6M99q1WrFlauXAlPT09Hd4mIiFwQgxuRhdLT0+Hh4cHQRkREDsPgRkRERKQSPMeNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsEtjw8//BAajcbq25w5cxy9CaQiycnJ8Pf3h0ajQUREBHJycgpcRqvV4sknn1T2uZUrVxZDT8nZsNZRcWKtsx0GtzwOHDhQpOVr165to56QKwgKCsKQIUMAAJcvX8bq1asLXGbs2LHYtGkTAGDy5Ml4/vnn7dpHck6sdVScWOtsRyMi4uhOlCTnzp1DWlqaRW3v3buH3r1749KlSwCA+vXrY9euXfDx8bFnF8nJXLp0CVWqVEF2djaaNGmCP//802zbRYsWYejQoQCAHj164IcffoBGoymurpITYa2j4sZaZyNCVklPT5c2bdoIAAEg1atXl1u3bjm6W6RS/fr1U/alPXv2mGzz+++/i6enpwCQevXqSWpqajH3klwRax3ZEmtd0fGrUitkZ2ejZ8+e2LFjBwAgKioKW7duRUhIiGM7Rqo1fvx45dOkqXOHkpKS0LNnT2RnZ6NChQpYu3Ytj3aQ3bHWka2x1hUdg1sh6XQ6vPjii1i/fj0AICwsDFu3bkVYWJiDe0ZqVrNmTXTq1AkA8OOPP+LixYvKvLt376Jz5864c+cOvL29sWbNGoSHhzuqq+QiWOvIHljrio7BrZCGDRuG77//HgAQEhKCLVu2oHLlyg7uFTmD8ePHA3gwkuqzzz5T7j/33HM4ceIEAOCrr75C48aNHdZHch2sdWQvrHVF5OjvatXk9ddfV76bL1OmjBw4cMDRXXI5r732msyYMUM2bNggN27ccHR3bK5p06YCQAIDA+X+/fsyatQoZZ97++23Hd09chGsdY43fPhweffdd2XTpk1OeU4ha531OKrUQtOnT8e0adMAAD4+Pti8eTNatmzp2E65oJiYGJw6dUr5OyIiAg0aNDC4PfTQQw7sYdH8/PPP6N69OwCgXbt22LZtGwCOqqLiw1pXMkRERODy5cvK35GRkUa1Ts3nGrLWWY/BzQJz5szBmDFjAAClSpXC2rVr8cQTTzi4V64pb3AzRc1hTkRQvXp1g23kpReouLDWlRx5g5spag5zrHXWY3ArwFdffYXBgwdDRODu7o7vv/8ePXr0cHS3XJYlwc0UNYW5JUuWYPDgwQCA0NBQ7N+/nyfokt2x1pUslgQ3U9QU5ljrrOSo72jV4Pvvvxc3NzcBIBqNRuLi4hzdJbMmTpyonB/Am2W3iIgI6datm2zZssXRL5+B33//XenjtGnTHN0dcgFqqnWvvfaaw2uH2m6RkZHSvXt32bFjh6NfPgOsddbhqFIzNmzYgBdeeAE6nQ4AMHfuXPTv39/BvTLv/fffd3QXVOfSpUv45ZdfsHPnTkd3xcDhw4eV+7GxsQ7rB7kGtdW6uXPnOroLqnPhwgX89NNP+f5SgSOw1lmHwc2EP/74Q7kAIAC8++67ePXVVx3cK7IXN7eS9TY4cuSIcp/FjOyJtc61lLQT/lnrrOPh6A6UNAcOHEDnzp2Rnp4O4MH1ZiZNmuTgXhVs06ZNyqgcZzZr1qwir0Oj0aBq1apo2LAhGjRogKeeesoGPbMd/afQoKAgREZGOrYz5LTUWuvWrl2LXbt2ObobdmerWhcTE6PUus6dO9ugZ7bDWmcdDk7I5dixY2jdujXu3LkDABg+fDgWLlzo4F5RboUdnJA3pDVo0AD16tWDv7+/HXtpvezsbPj5+SErKwtt2rTB77//7ugukRNirSv5Cjs4IW9Ia9CgAWJjY+Hn52fHXlqPtc56POL2/5KSktChQwelkPXt2xfz5893cK+oMNQW0kw5ceIEsrKyAPCrA7IP1jr1U1tIM4W1znoMbgCuXLmC9u3b49q1awCArl27Ii4ursSd+0T/4wwhzRSerEv2xFqnPs4Q0kxhrbOeywe35ORktG/fHufPnwfw4Ku4t99+GydPnrR4HREREQgICLBTDym3yZMno2LFik4R0kzhybpkL6x16jJjxgw88sgjThHSTGGts57Ln+P23Xff4fnnny/SOvbs2YNmzZrZqEfkytq1a4ft27ejVKlSuH//Pjw9PR3dJXISrHVUkrDWWc/lj48nJCQUaXl3d3d+WiCb0X8KrVGjBgsZ2RRrHZUkrHXWc/kjbkRERERq4fJH3IiIiIjUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCUY3IiIiIhUgsGNiIiISCX+D40c5WtkujfRAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pgm = daft.PGM(dpi=DPI, grid_unit=GRID_UNIT, node_ec=NODE_EC)\n", + "\n", + "# data generating graph\n", + "pgm.add_node(\"a\", \"$RV$\", 0, 1)\n", + "pgm.add_node(\"z\", \"$Z$\", 0, 0)\n", + "pgm.add_node(\"x\", \"$\\mathbf{X}$\", 1, 1)\n", + "pgm.add_node(\"y\", \"$Y$\", 1, 0)\n", + "pgm.add_edge(\"a\", \"z\")\n", + "pgm.add_edge(\"a\", \"y\")\n", + "pgm.add_edge(\n", + " \"a\",\n", + " \"x\",\n", + " plot_params={\"ec\": \"grey\", \"lw\": 1.5, \"ls\": \":\", \"head_length\": 0, \"head_width\": 0},\n", + ")\n", + "pgm.add_edge(\"z\", \"y\")\n", + "pgm.add_edge(\"x\", \"y\")\n", + "pgm.add_text(0, 1.3, \"Data generating graph\")\n", + "\n", + "# limiting graph\n", + "x_offset = 2\n", + "pgm.add_node(\"a2\", r\"$RV \\rightarrow \\lambda$\", 0 + x_offset, 1)\n", + "pgm.add_node(\"z2\", \"$Z$\", 0 + x_offset, 0)\n", + "pgm.add_node(\"x2\", \"$\\mathbf{X}$\", 1 + x_offset, 1)\n", + "pgm.add_node(\"y2\", \"$Y$\", 1 + x_offset, 0)\n", + "pgm.add_edge(\"a2\", \"z2\")\n", + "pgm.add_edge(\"z2\", \"y2\")\n", + "pgm.add_edge(\"x2\", \"y2\")\n", + "pgm.add_text(x_offset, 1.3, \"Limiting graph\")\n", + "\n", + "pgm.render();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see from the data generating graph (left) that the $RV$ is a confounding variable as it influences both the treatment $Z$ and the outcome $Y$. \n", + "\n", + "If we tried to identify the causal effect of $Z \\rightarrow Y$ by conditioning on the running variable ($RV=rv$), we would eliminate any variation in $Z$ or $Y$ caused by $RV$. And because $Z$ is constant for any given value of $RV$, then the $Z \\rightarrow Y$ path would disappear and we could not estimate the causal effect.\n", + "\n", + "Identification of the causal effect of $Z \\rightarrow Y$ is done with a limiting graph (right). The $RV$ node is replaced by a subset of the data where $RV$ is close to the cutoff value $\\lambda$, hence the name \"limiting graph\" and the symbol $RV \\rightarrow \\lambda$.\n", + "\n", + "In the limit, this eliminates variation in the running variable and so breaks the $RV \\rightarrow Y$ path. The causal effect of $Z \\rightarrow Y$ can be estimated by comparing the outcomes of units just above and just below the cutoff value $\\lambda$.\n", + "\n", + "Readers are referred to {cite:t}`steiner2017graphical` and [Chapter 6](https://mixtape.scunning.com/06-regression_discontinuity) of {cite:t}`cunningham2021causal` who discuss limiting graphs in more detail. Chapter 20 of {cite:t}`huntington2021effect` also covers regression discontinuity designs, but presents simplified (and non-kosher, in his own words) causal DAG." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## References\n", + ":::{bibliography}\n", + ":filter: docname in docnames\n", + ":::" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "CausalPy", + "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.8" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/source/references.bib b/docs/source/references.bib index 08034e5a..93acf8a8 100644 --- a/docs/source/references.bib +++ b/docs/source/references.bib @@ -76,3 +76,28 @@ @book{shadish_cook_cambell_2002 year={2002}, publisher={Houghton Mifflin Boston, MA} } + +@article{steiner2017graphical, + title={Graphical models for quasi-experimental designs}, + author={Steiner, Peter M and Kim, Yongnam and Hall, Courtney E and Su, Dan}, + journal={Sociological methods \& research}, + volume={46}, + number={2}, + pages={155--188}, + year={2017}, + publisher={SAGE Publications Sage CA: Los Angeles, CA} +} + +@book{cunningham2021causal, + title={Causal inference: The mixtape}, + author={Cunningham, Scott}, + year={2021}, + publisher={Yale university press} +} + +@book{huntington2021effect, + title={The effect: An introduction to research design and causality}, + author={Huntington-Klein, Nick}, + year={2021}, + publisher={Chapman and Hall/CRC} +} diff --git a/pyproject.toml b/pyproject.toml index abd63362..eb8d61df 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -52,6 +52,7 @@ dependencies = [ dev = ["pathlib", "pre-commit", "twine", "interrogate"] docs = [ "ipykernel", + "daft", "linkify-it-py", "myst-nb", "pathlib",