causal-inference-3-frontdoor-adjustment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\newcommand{P}[1]{\\mathrm{P}\\left( #1 \\right)}\n",
    "\\newcommand{Pc}[2]{\\mathrm{P}\\left( #1 \\mid #2 \\right)}\n",
    "\\newcommand{In}[2]{ #1 \\perp\\!\\!\\!\\perp #2}\n",
    "\\newcommand{Cin}[3]{ #1 \\perp\\!\\!\\!\\perp #2 \\, | \\, #3}\n",
    "\\newcommand{do}[1]{\\mathrm{do}\\left( #1 \\right)}\n",
    "\\newcommand{E}[1]{\\mathrm{E}\\left[ #1 \\right]}\n",
    "\\newcommand{Ec}[2]{\\mathrm{E}\\left[ #1 \\mid #2 \\right]}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a [previous post](http://www.degeneratestate.org/posts/2018/Jul/10/causal-inference-with-python-part-2-causal-graphical-models/) I looked at how we can describe causal systems mathematically using the framework of causal graphical model. These models give us the ability to be precise about how we can phrase and answer causal questions of the form:\n",
    "\n",
    "> \"What effect does changing $X$ have on $Y$?\"\n",
    "\n",
    "These causal graphical model show us exactly why causality is difficult: if there exist [\"backdoor paths\" - or confounding variables](https://en.wikipedia.org/wiki/Confounding), common causes for both $X$ and $Y$, then it is possible that any observed correlation between $X$ and $Y$ is due to these confounding paths, and not a direct causal relationship between $X$ and $Y$.\n",
    "\n",
    "In situations where we observe the confounding variables in a causal graphical model we can overcome this limitation by \"adjusting\" for the backdoor path (sometimes called \"covariate adjustment\", \"backdoor adjustment, \"controlling for variables\"). My [first post](http://www.degeneratestate.org/posts/2018/Mar/24/causal-inference-with-python-part-1-potential-outcomes/) contained a number of techniques for doing this. \n",
    "\n",
    "This leaves the question: if we do not observe enough variables to adjust for a backdoor, can we still make causal inferences? \n",
    "\n",
    "In some situations the answer to this question is yes. We will be examining one of these situations in this post. As in my previous post I will be using the [causalgraphicalmodels](https://github.com/ijmbarr/causalgraphicalmodels) python package for my examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we begin, we need to introduce some slightly new notation:\n",
    "\n",
    "# Unobserved Variables\n",
    "\n",
    "When we describe our system as a causal graphical model, the effect of unmeasured we would like to capture are those that cause [confounding](https://en.wikipedia.org/wiki/Confounding) relationships between variables because these change the [conditional independence](https://en.wikipedia.org/wiki/Conditional_independence) structure implied by the model. Unobserved variables which are parents of only one variable, or which act as mediators between variables can be absorbed into the non-parametric structure of the model without modification.\n",
    "\n",
    "To represent confounding relationships between two variables, we add a dashed edge between, with double headed arrows. We can create such a structure using the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first load some standard libraries\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"89pt\" height=\"228pt\"\n",
       " viewBox=\"0.00 0.00 89.00 227.60\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 223.6)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-223.6 85,-223.6 85,4 -4,4\"/>\n",
       "<!-- z -->\n",
       "<g id=\"node1\" class=\"node\"><title>z</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">z</text>\n",
       "</g>\n",
       "<!-- y -->\n",
       "<g id=\"node3\" class=\"node\"><title>y</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"54\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"54\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">y</text>\n",
       "</g>\n",
       "<!-- z&#45;&gt;y -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>z&#45;&gt;y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M33.3986,-72.411C36.5136,-64.3352 40.3337,-54.4312 43.8346,-45.3547\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"47.1265,-46.5458 47.4597,-35.9562 40.5955,-44.0267 47.1265,-46.5458\"/>\n",
       "</g>\n",
       "<!-- x -->\n",
       "<g id=\"node2\" class=\"node\"><title>x</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">x</text>\n",
       "</g>\n",
       "<!-- x&#45;&gt;z -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>x&#45;&gt;z</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M27,-143.697C27,-135.983 27,-126.712 27,-118.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"30.5001,-118.104 27,-108.104 23.5001,-118.104 30.5001,-118.104\"/>\n",
       "</g>\n",
       "<!-- Unobserved_0 -->\n",
       "<g id=\"node4\" class=\"node\"><title>Unobserved_0</title>\n",
       "<ellipse fill=\"black\" stroke=\"black\" cx=\"54\" cy=\"-217.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
       "</g>\n",
       "<!-- Unobserved_0&#45;&gt;x -->\n",
       "<g id=\"edge3\" class=\"edge\"><title>Unobserved_0&#45;&gt;x</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M53.5314,-215.866C51.8185,-212.453 45.6865,-200.235 39.7489,-188.403\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"42.7305,-186.541 35.1169,-179.174 36.4742,-189.681 42.7305,-186.541\"/>\n",
       "</g>\n",
       "<!-- Unobserved_0&#45;&gt;y -->\n",
       "<g id=\"edge4\" class=\"edge\"><title>Unobserved_0&#45;&gt;y</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M54.3041,-215.904C55.6687,-211.832 61.2373,-194.627 63,-180 68.7428,-132.345 67.0322,-119.83 63,-72 62.2834,-63.4991 60.9311,-54.3119 59.4884,-46.0122\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"62.9142,-45.2892 57.6538,-36.0925 56.0309,-46.5623 62.9142,-45.2892\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x7fef2427be80>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from causalgraphicalmodels import CausalGraphicalModel\n",
    "\n",
    "system_with_hidden_confounder = CausalGraphicalModel(\n",
    "    nodes=[\"x\", \"y\", \"z\"],\n",
    "    edges=[(\"x\", \"z\"), (\"z\", \"y\")],\n",
    "    latent_edges=[(\"x\", \"y\")]\n",
    ")\n",
    "\n",
    "system_with_hidden_confounder.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "It's important to realise that this notation is only a convenience: we could simply add another variable as a parent to each pair of confounded variables, and reason about the resulting graph. The reason we include them as hidden variables is that it allows us to automatically reason about whether or not it is possible to infer the results of interventions. \n",
    "\n",
    "Consider the graph above. The introduction of the hidden variable means that it is no longer possible to always use the backdoor adjustment to reason about about interventions - there is no valid backdoor adjustment for the graph above because the confounding variable is unobserved:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset()"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system_with_hidden_confounder.get_all_backdoor_adjustment_sets(\"x\", \"y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But all is not lost. The fact we observe variable $Z$ allows us to take another approach\n",
    "\n",
    "# Front Door Adjustment\n",
    "\n",
    " Let's begin by writing our target:\n",
    "\n",
    "$\\Pc{Y}{\\do{X}}$\n",
    "\n",
    "We can expand this by including variable $Z$ (as in my previous post, I'm assuming all variables are discrete - if they are continuous, replace the sum with a integral)\n",
    "\n",
    "$ = \\sum_{Z}\\Pc{Y}{\\do{X}, Z}\\Pc{Z}{\\do{X}}$\n",
    "\n",
    "Because there is only one directed path from $X$ to $Z$, their conditional distribution is the same as their interventional distribution and we can replace the $\\do{\\dots}$ in the second conditional with the observational value:\n",
    "\n",
    "$ = \\sum_{Z}\\Pc{Y}{\\do{X}, Z}\\Pc{Z}{X}$\n",
    "\n",
    "We can't make a similar adjustment to the first term, because there is a backdoor path between $X$ and $Y$. We can however replace the observed $Z$ with its intervential value $\\do{Z}$:\n",
    "\n",
    "$ = \\sum_{Z}\\Pc{Y}{\\do{X}, \\do{Z}}\\Pc{Z}{X}$\n",
    "\n",
    "Doing this, we can remove $\\do{X}$ from the first conditional, because when we intervene on $Z$, we \"screen off\" the intervention on $X$:\n",
    "\n",
    "$ = \\sum_{Z}\\Pc{Y}{\\do{Z}}\\Pc{Z}{X}$\n",
    "\n",
    "It might not look like we've achieved much yet by replacing one intervention with another, but $\\Pc{Y}{\\do{Z}}$ is something we know how to manipulate: we can use backdoor adjustment. Conditioning again on $X$ blocks the back door between $Y$ and $Z$, allowing us to write:\n",
    "\n",
    "$ = \\sum_{Z}\\Pc{Z}{X}\\left(\\sum_{X'}\\Pc{Y}{Z, X'}\\P{X'}\\right)$\n",
    "\n",
    "This is pretty amazing: we've managed to convert a causal statement into one about observations _despite there being unmeasured confounders_. There are still some strong assumptions here - we need to have observed mediation mechanism between $X$ and $Y$ that is somehow free from direct influence of the confounders, but it is not impossible to imagine situations where this might hold. I found it surprising that this kind of manipulation was possible.\n",
    "\n",
    "We can verify that in the causal graphical model above, $Z$ is a valid variable for front door adjustment using the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# {\"z\"} is a valid frontdoor adjustment set\n",
    "system_with_hidden_confounder.is_valid_frontdoor_adjustment_set(\"x\", \"y\", {\"z\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# but the empty set is not\n",
    "system_with_hidden_confounder.is_valid_frontdoor_adjustment_set(\"x\", \"y\", {})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can compute every valid frontdoor adjustment set using the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset({frozenset({'z'})})"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system_with_hidden_confounder.get_all_frontdoor_adjustment_sets(\"x\", \"y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The set of manipulations I carried out above to derived the frontdoor adjustment are not at all obvious. In deriving them I use a few imprecise arguments. The manipulation I carried out can be made precise, and are known as [do-calculus](http://ftp.cs.ucla.edu/pub/stat_ser/r402.pdf) - but that is a story for another blog post."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimating the ATE using the Front Door Adjustment\n",
    "\n",
    "In the previous section we derived the relationship \n",
    "\n",
    "$\\Pc{Y}{\\do{X}} = \\sum_{Z}\\Pc{Z}{X}\\left(\\sum_{X'}\\Pc{Y}{Z, X'}\\P{X'}\\right)$\n",
    "\n",
    "for valid frontdoor adjustment sets $Z$. \n",
    "\n",
    "The remaining question is: given observational data, how do we calculate it?\n",
    "\n",
    "To simplify things a bit, I am going to look at the simpler question:\n",
    "\n",
    "$\\Ec{Y}{\\do{X}} = \\sum_{Z}\\Pc{Z}{X}\\left(\\sum_{X'}\\Ec{Y}{Z, X'}\\P{X'}\\right)$\n",
    "\n",
    "$ = \\sum_{X', Y, Z} Y \\, \\Pc{Y}{Z, X'} \\Pc{Z}{X} \\P{X'}$\n",
    "\n",
    "Which is about the expectation of $Y$, under intervention $X$. \n",
    "\n",
    "We can reduce this to a single number, the average treatment effect (ATE):\n",
    "\n",
    "$ATE = \\Ec{Y}{\\do{X=1}} - \\Ec{Y}{\\do{X=0}}$\n",
    "\n",
    "The most direct way to estimate this is to estimate each separate conditional probability and combine the results.\n",
    "\n",
    "Let's start by defining a `StructuralCausalModel` we can sample from:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"89pt\" height=\"260pt\"\n",
       " viewBox=\"0.00 0.00 89.00 260.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 256)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-256 85,-256 85,4 -4,4\"/>\n",
       "<!-- z -->\n",
       "<g id=\"node1\" class=\"node\"><title>z</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">z</text>\n",
       "</g>\n",
       "<!-- y -->\n",
       "<g id=\"node4\" class=\"node\"><title>y</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"54\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"54\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">y</text>\n",
       "</g>\n",
       "<!-- z&#45;&gt;y -->\n",
       "<g id=\"edge4\" class=\"edge\"><title>z&#45;&gt;y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M33.3986,-72.411C36.5136,-64.3352 40.3337,-54.4312 43.8346,-45.3547\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"47.1265,-46.5458 47.4597,-35.9562 40.5955,-44.0267 47.1265,-46.5458\"/>\n",
       "</g>\n",
       "<!-- x -->\n",
       "<g id=\"node2\" class=\"node\"><title>x</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">x</text>\n",
       "</g>\n",
       "<!-- x&#45;&gt;z -->\n",
       "<g id=\"edge3\" class=\"edge\"><title>x&#45;&gt;z</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M27,-143.697C27,-135.983 27,-126.712 27,-118.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"30.5001,-118.104 27,-108.104 23.5001,-118.104 30.5001,-118.104\"/>\n",
       "</g>\n",
       "<!-- u -->\n",
       "<g id=\"node3\" class=\"node\"><title>u</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"54\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"54\" y=\"-230.3\" font-family=\"Times,serif\" font-size=\"14.00\">u</text>\n",
       "</g>\n",
       "<!-- u&#45;&gt;x -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>u&#45;&gt;x</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M47.6014,-216.411C44.4864,-208.335 40.6663,-198.431 37.1654,-189.355\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"40.4045,-188.027 33.5403,-179.956 33.8735,-190.546 40.4045,-188.027\"/>\n",
       "</g>\n",
       "<!-- u&#45;&gt;y -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>u&#45;&gt;y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M57.6538,-215.908C59.6758,-205.569 61.9808,-192.09 63,-180 67.0322,-132.17 67.0322,-119.83 63,-72 62.2834,-63.4991 60.9311,-54.3119 59.4884,-46.0122\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"62.9142,-45.2892 57.6538,-36.0925 56.0309,-46.5623 62.9142,-45.2892\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x7feeed4a3cc0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from causalgraphicalmodels import StructuralCausalModel\n",
    "from causalgraphicalmodels.csm import discrete_model, linear_model\n",
    "\n",
    "fd0 = StructuralCausalModel({\n",
    "    \"u\": lambda n_samples: np.random.binomial(n=1, p=0.5, size=n_samples),\n",
    "    \"x\": discrete_model([\"u\"], {(0,): (0,1), (1,): (1,0)}),\n",
    "    \"z\": discrete_model([\"x\"], {(0,): (0.95, 0.05), (1,): (0.05, 0.95)}),\n",
    "    \"y\": discrete_model([\"z\", \"u\"], {\n",
    "        (0, 0): (0.95, 0.05), \n",
    "        (0, 1): (0.95, 0.05),\n",
    "        (1, 0): (0.50, 0.50), \n",
    "        (1, 1): (0.95, 0.05),\n",
    "    }),\n",
    "})\n",
    "\n",
    "fd0.cgm.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This model relates four binary variables: $X$, $Z$, $Y$ and $U$. For the rest of this post I will assume that $U$ is unobserved. If it was we could adjust for it directly. Instead we have to use the frontdoor. \n",
    "\n",
    "Before that, we should have some idea what the actual ATE is. The following functions allow us to empirically estimate the ATE from a `StructuralCausalModel`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_conditional_expectation(df, x, y):\n",
    "    \"\"\"\n",
    "    Estimate E[Y|X=1] - E[Y|X=0]\n",
    "    \n",
    "    from a dataframe `df` of samples.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    df: pandas.DataFrame\n",
    "    \n",
    "    x: str\n",
    "    \n",
    "    y: str\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    estimates: tuple[float, float]\n",
    "        estiamted difference and standard deviation\n",
    "    \"\"\"\n",
    "    a = df.loc[lambda df: df[x] == 0]\n",
    "    b = df.loc[lambda df: df[x] == 1]\n",
    "    \n",
    "    n_a = a.shape[0]\n",
    "    n_b = b.shape[0]\n",
    "    \n",
    "    delta = b[y].mean() - a[y].mean()\n",
    "    delta_std = np.sqrt(a[y].var() / n_a + b[y].var() / n_b)\n",
    "    \n",
    "    return (delta, delta_std)\n",
    "\n",
    "\n",
    "def ab_test(scm, x, y, n=10000):\n",
    "    \"\"\"\n",
    "    Simulate an A/B on StructuralCausalModel scm\n",
    "    to estimate the quantity:\n",
    "    \n",
    "    E[Y|do(X) = 1] - E[Y|do(X) = 0]\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    scm: StructuralCausalmodel\n",
    "    \n",
    "    x: str\n",
    "    \n",
    "    y: str\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    estimate: tuple[estimates, estiamte_std]\n",
    "    \n",
    "    \"\"\"\n",
    "    scm_do = scm.do(x)\n",
    "    \n",
    "    n_a = int(n / 2)\n",
    "    n_b = n - n_a\n",
    "    set_variable = np.array([0]*n_a + [1]*n_b)\n",
    "    \n",
    "    samp = (\n",
    "        scm_do\n",
    "        .sample(\n",
    "            set_values={x: set_variable}, \n",
    "            n_samples=n)\n",
    "    )\n",
    "    \n",
    "    return estimate_conditional_expectation(samp, x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated ATE: 0.208\n"
     ]
    }
   ],
   "source": [
    "print(\"Estimated ATE: {:.3f}\".format(ab_test(fd0, \"x\", \"y\")[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now need a function to make the same calculation using our frontdoor adjustment formula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estiamte_ate_frontdoor_direct(df, x, y, zs):\n",
    "    \"\"\"\n",
    "    Estiamte the ATE of a system from a dataframe of samples `ds`\n",
    "    using frontdoor adjustment directly on ML estimates of probability.\n",
    "    \n",
    "    E[Y|do(X) = x'] = \\sum_{x,y,z} y P[y|x,z] P(z|x') P(X)\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    df: pandas.DataFrame\n",
    "    x: str\n",
    "    y: str\n",
    "    zs: list[str]\n",
    "    \n",
    "    Returns \n",
    "    -------\n",
    "    ATE: float\n",
    "    \"\"\"\n",
    "    zs_unique = [tuple(a) for a in np.unique(df[zs].values, axis=0)]\n",
    "    y_unique = np.unique(df[y].values, axis=0)    \n",
    "    \n",
    "    # P(X)\n",
    "    p_x = {\n",
    "        x_: np.mean(df[x] == x_)\n",
    "        for x_ in (0,1)\n",
    "    }\n",
    "    \n",
    "    # P(Z|X)\n",
    "    p_z_x = {\n",
    "        (x_, z_): np.mean(df\n",
    "            .loc[lambda df_: df_[x] == x_]\n",
    "            [zs]\n",
    "            .values == z_)\n",
    "        for x_ in (0,1)\n",
    "        for z_ in zs_unique\n",
    "    }\n",
    "    \n",
    "    # P(Y|X,Z)\n",
    "    p_y_xz = {\n",
    "        (x_, y_, z_): np.mean(df\n",
    "             .loc[lambda df_: df_[x] == x_]\n",
    "             .loc[lambda df_: (df_[zs].as_matrix() == z_).squeeze()]\n",
    "             [y]\n",
    "             .values == y_)\n",
    "        for x_ in (0,1)\n",
    "        for y_ in y_unique\n",
    "        for z_ in zs_unique\n",
    "    }\n",
    "    \n",
    "    \n",
    "    # ATE\n",
    "    E_y_do_x_0 = 0.0\n",
    "    E_y_do_x_1 = 0.0\n",
    "        \n",
    "    for y_ in y_unique:\n",
    "        for zs_ in  zs_unique:\n",
    "            for x_ in [0,1]:\n",
    "                E_y_do_x_0 += y_ * p_y_xz[(x_, y_, zs_)] * p_z_x[(0, zs_)] * p_x[x_]\n",
    "                E_y_do_x_1 += y_ * p_y_xz[(x_, y_, zs_)] * p_z_x[(1, zs_)] * p_x[x_]\n",
    "    \n",
    "    return E_y_do_x_1 - E_y_do_x_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.20228124029189454"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds = fd0.sample(10000)\n",
    "\n",
    "estiamte_ate_frontdoor_direct(ds, \"x\", \"y\", [\"z\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It works! \n",
    "\n",
    "But what happens when we have continuous values? Or sparse discrete values? In these cases directly estimating the conditional probabilities is impractical.\n",
    "\n",
    "There are a number of approaches you can take - which ones you want to use will depend on exactly which assumptions you are will to make about you data. Let's start by looking at what we can do when assume the data is linear. In doing so we will get a feeling for what the front door adjustment is doing.\n",
    "\n",
    "# Linear Relationship\n",
    "\n",
    "The linear approximation assumes that there are linear relationships between $X$, $Z$ and $Y$, or that:\n",
    "\n",
    "$Z = a_{z} + b_{z, x} X + \\epsilon_{z}$\n",
    "\n",
    "$Y = a_{y} + b_{y, z} Z + b_{y, u} U + \\epsilon_{y}$\n",
    "\n",
    "If these are the _true_ data generating equations then the causal effect of changing $X$ on $Y$ is given by the product of coefficients $ b_{z, x}  b_{y, z}$. Simple regression doesn't allow us to estimate this value because the unmeasured variable $U$ is a parent of both $X$ and $Y$ - and if $X$ and $U$ end up [collinear](https://en.wikipedia.org/wiki/Multicollinearity), our estimate of $b_{y, z}$ cannot be trusted.\n",
    "\n",
    "One way to side-step this issue is to perform a double regression:\n",
    "\n",
    "1. First we regress $Z$ on $X$. This not only gives us the coefficient $b_{z, x}$, but the way to estimate $\\Ec{Z}{X}$. If this estimate of $\\Ec{Z}{X}$ is accurate, we can use it to find the values of the noise in the first equation. Specifically: $ Z - \\Ec{Z}{X} = \\epsilon_{z}$\n",
    "\n",
    "2. This noise is only introduced at the node for $Z$, which means that by assumption it is independent of the confounding variable $U$. Regressing $\\epsilon_{z}$ on $Y$ then gives us an estimate of the coefficient $b_{y, z}$ where we don't have to worry about the collinearity problem.\n",
    "\n",
    "Together this procedure gives us the following estimate of "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "\n",
    "def estimate_ate_frontdoor_linear(df, x, y, z):\n",
    "    \"\"\"\n",
    "    Estiamte the ATE of a system from a dataframe of samples `df`\n",
    "    using frontdoor adjustment and assuming linear models.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    df: pandas.DataFrame\n",
    "    x: str\n",
    "    y: str\n",
    "    z: str\n",
    "    \n",
    "    Returns \n",
    "    -------\n",
    "    ATE: float\n",
    "    \"\"\"\n",
    "    x = df[x].values\n",
    "    y = df[y].values\n",
    "    z = df[z].values\n",
    "    \n",
    "    z_x_model = sm.OLS(z, sm.add_constant(x)).fit()\n",
    "    \n",
    "    z_bar = z_x_model.predict(sm.add_constant(x))\n",
    "    z_prime = z - z_bar\n",
    "\n",
    "    y_z_model = sm.OLS(y, sm.add_constant(z_prime)).fit()\n",
    "    \n",
    "    return y_z_model.params[1] * z_x_model.params[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how it does:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Different estimates of the ATE for a linear system\n",
      "--------------------------------------------------\n",
      "Naive estimator: 0.004\n",
      "AB test: 8.033\n",
      "Linear estimator: 8.012\n"
     ]
    }
   ],
   "source": [
    "# first define a linear model\n",
    "fd1 = StructuralCausalModel({\n",
    "    \"u\": lambda n_samples: np.random.binomial(n=1, p=0.5, size=n_samples),\n",
    "    \"x\": discrete_model([\"u\"], {(0,): (0,1), (1,): (1,0)}),\n",
    "    \"z\": linear_model([\"x\"], [8], -3),\n",
    "    \"y\": linear_model([\"z\", \"u\"], [1, 8], -2),\n",
    "})\n",
    "\n",
    "# sample from the model\n",
    "n = int(10e4)\n",
    "ds = fd1.sample(n)\n",
    "\n",
    "# evaluate the model\n",
    "print(\"Different estimates of the ATE for a linear system\")\n",
    "print(\"--------------------------------------------------\")\n",
    "print(\"Naive estimator: {:.3f}\".format(estimate_conditional_expectation(ds, \"x\", \"y\")[0]))\n",
    "print(\"AB test: {:.3f}\".format(ab_test(fd1, \"x\", \"y\", n=n)[0]))\n",
    "print(\"Linear estimator: {:.3f}\".format(estimate_ate_frontdoor_linear(ds, \"x\", \"y\", \"z\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks reasonably good.\n",
    "\n",
    "It is worth reflecting on what this process is doing: we are using the randomness introduced by the intermediate variable $Z$ to try and disentangle the confounded relationship between $X$ and $Y$. In some ways this is exactly what the more general frontdoor adjustment is doing: we combine our knowledge of the effect of an intervention at $X$ on $Z$ with the effect of an intervention on $Z$ using $X$ to remove it's bias."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Beyond linear models \n",
    "\n",
    "Linear models are powerful when the data is approximately linear, but often this is too restrictive an assumption. We may also want to introduce categoric variables which may not fit into a linear model framework. \n",
    "\n",
    "There are a number of approaches we can take. The one I am going to look at here is the semi-parametric estimator inspired by [this paper](https://arxiv.org/pdf/1210.4654.pdf). It requires us to make some assumptions about how we can model parts of our data, and leaves others to be estimated non-parametrically from the data. It is not only estimator that exists. \n",
    "\n",
    "The estimator is:\n",
    "\n",
    "$\\Ec{Y}{\\do{X=x}} = \\frac{1}{N}\\sum_{i} \\frac{g(z_{i}) I(x_{i} = x)}{\\hat{p}(x_{i})}$\n",
    "\n",
    "Where \n",
    "\n",
    "- $I(x_{i} = x)$ is the [indicator function](https://en.wikipedia.org/wiki/Indicator_function): one when $x_{i} = x$ and zero otherwise\n",
    "- $\\hat{p}(x_{i})$ is an estimate of the propensity of $x_{i}$ being observed\n",
    "- $g(z_{i})$ is an estimate of $\\Ec{\\Ec{y}{x,z}}{z}$\n",
    "- the sum runs over all the samples\n",
    "\n",
    "To understand how this estimator works, let's start expanding the probability distribution of our samples using Bayes theorem to get:\n",
    "\n",
    "$\\P{x,y,z} = \\Pc{y}{x,z}\\Pc{z}{x}\\P{x}$\n",
    "\n",
    "Averaging our estimator over our samples gives us an estimate of \n",
    "\n",
    "$\\E{\\frac{g(z_{i}) I(x_{i} = x')}{\\hat{p}(x_{i})}} = \\sum_{x,y,z} \\Pc{y}{x,z}\\Pc{z}{x}\\P{x} \\frac{g(z) I(x = x')}{\\hat{p}(x)}$\n",
    "\n",
    "If we assume that our model $\\hat{p}(x_{i})$ is accurate we can simplify this to \n",
    "\n",
    "$\\sum_{x,y,z} \\Pc{y}{x,z}\\Pc{z}{x} g(z) I(x = x')$\n",
    "\n",
    "The summing over $x$ and keeping only the non-zero values we get\n",
    "\n",
    "$\\sum_{z} \\Pc{z}{x'} g(z)$\n",
    "\n",
    "And replacing $g$ with the thing it is trying to estimate, we get\n",
    "\n",
    "$\\sum_{x,y,z} y \\Pc{y}{x,z}\\P{x}\\Pc{z}{x'}$\n",
    "\n",
    "Which, if you compare it with the definition of the frontdoor adjustment, is exactly the quantity we were trying to estimate. \n",
    "\n",
    "It is important to note that in demonstrating that this estimator matches the frontdoor adjustment we have ignored a number of import questions around asymptotic efficiency of the estimator. This is important and will affect the quality of your results, but they are also a question for another blog post.\n",
    "\n",
    "We have also not discussed how to calculate $g(z)$ or $\\hat{p}(x)$ - this is a much easier problem to address. For binary $x$, $\\hat{p}(x)$ can be estimated from the sample averages and $g(z)$ can be calculated using you favorite machine learning technique to estimate a conditional mean (although we must be careful to use stacking when computing the estimate to avoid training on the same data we are predicting on).\n",
    "\n",
    "Some example code is below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier\n",
    "from sklearn.model_selection import KFold\n",
    "\n",
    "def stack_predict_and_sum(x, z, y, model, n_splits=3):\n",
    "    \"\"\"\n",
    "    Estimate E[ E[Y|X,Z] |Z]\n",
    "    \n",
    "    Assumes that X is a binary value, and uses `model` to\n",
    "    predict the conditional mean. To prevent overfitting on \n",
    "    the data itself, we use stacking to predict on data that \n",
    "    wasn't used for fitting.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    x: numpy.array [Nx1] binary array\n",
    "    y: numpy.array [Nx1] target array\n",
    "    z: numpy.array [NxM] feature array\n",
    "    model: sklearn.estimator\n",
    "    \n",
    "    Returns \n",
    "    -------\n",
    "    numpy.array [Nx1]\n",
    "    \"\"\"\n",
    "    kf = KFold(n_splits=n_splits, shuffle=True)\n",
    "    output = np.zeros_like(y, dtype=np.float)\n",
    "    p_x = np.mean(x)\n",
    "    \n",
    "    for train_index, test_index in kf.split(x):\n",
    "        x_train = x[train_index]\n",
    "        x_test  = x[test_index]\n",
    "        \n",
    "        z_train = z[train_index]\n",
    "        z_test  = z[test_index]\n",
    "        \n",
    "        y_train = y[train_index]\n",
    "        y_test  = y[test_index]      \n",
    "        \n",
    "        # train model on train set\n",
    "        X = np.concatenate([x_train, z_train], axis=1)\n",
    "        model.fit(X, y_train.squeeze())\n",
    "        \n",
    "        # predict model on z from test set and x=0, x=1\n",
    "        X_pred_x_0 = np.concatenate([np.zeros_like(x_test), z_test], axis=1)\n",
    "        X_pred_x_1 = np.concatenate([np.ones_like(x_test), z_test], axis=1)\n",
    "        \n",
    "        y_pred_x_0 = model.predict(X_pred_x_0).reshape(-1,1)\n",
    "        y_pred_x_1 = model.predict(X_pred_x_1).reshape(-1,1)\n",
    "        \n",
    "        output[test_index] = (p_x * y_pred_x_1 + (1-p_x)*y_pred_x_0)\n",
    "       \n",
    "    return output\n",
    "\n",
    "\n",
    "def estimate_intervention_frontdoor_nonlinear(ds, x, y, z, model_y_zx):\n",
    "    \"\"\"\n",
    "    Estiamte the ATE of a system from a dataframe of samples `df`\n",
    "    using a semi-parametric frontdoor adjustment with a arbritary \n",
    "    model of the conditional expectation. \n",
    "    \n",
    "    Stacking is used to prevent overfitting on the data.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    df: pandas.DataFrame\n",
    "    x: str\n",
    "    y: str\n",
    "    z: List[str]\n",
    "    model_y_zx: sklearn.estimator \n",
    "    \n",
    "    Returns \n",
    "    -------\n",
    "    ATE: float\n",
    "    \"\"\"\n",
    "    x = ds[[x]].values\n",
    "    y = ds[[y]].values\n",
    "    z = ds[z].values\n",
    "\n",
    "    mask_x0 = ds.x.values == 0\n",
    "\n",
    "    # P(x)\n",
    "    p_x_0 = np.mean(mask_x0)\n",
    "    p_x_1 = 1 - p_x_0\n",
    "    p_x = np.where(mask_x0, p_x_0, p_x_1)\n",
    "\n",
    "    # E[ E[Y|X,Z] |Z]\n",
    "    y_z_mean = stack_predict_and_sum(x=x, y=y, z=z, model=model_y_zx).squeeze()\n",
    "\n",
    "    # sum over samples and estimate:\n",
    "    E_y_do_x0 = np.mean((y_z_mean * mask_x0) / (1 - p_x))\n",
    "    E_y_do_x1 = np.mean((y_z_mean * (~mask_x0)) / (p_x))\n",
    "    \n",
    "    return E_y_do_x1 - E_y_do_x0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how it performs. \n",
    "\n",
    "To evaluate the estimator we are going to use the following non-linear frontdoor-like system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Different estimates of the ATE\n",
      "------------------------------\n",
      "Naive estimator: -0.338\n",
      "AB test: -0.001\n",
      "Linear test: 0.272\n",
      "Non-linear: -0.012\n"
     ]
    }
   ],
   "source": [
    "# Define the model\n",
    "fd2 = StructuralCausalModel({\n",
    "    \"u\": lambda n_samples: np.random.binomial(n=1, p=0.5, size=n_samples),\n",
    "    \"x\": discrete_model([\"u\"], {(0,): (1.00,0.00), (1,): (0.00, 1.00)}),\n",
    "    \"z\": lambda x, n_samples: np.where(\n",
    "        x == 0,\n",
    "        np.random.triangular(0,0,1, size=n_samples),\n",
    "        np.random.triangular(0,1,1, size=n_samples),\n",
    "    ),\n",
    "    \"y\": lambda z, u, n_samples: np.where(\n",
    "        u == 0, \n",
    "        1-(4*(z - 0.5)**2),\n",
    "        (4*(z - 0.5)**2)\n",
    "    ),\n",
    "})\n",
    "\n",
    "# sample from the model\n",
    "n = int(10e3)\n",
    "ds = fd2.sample(n)\n",
    "\n",
    "# evaluate the model\n",
    "print(\"Different estimates of the ATE\")\n",
    "print(\"------------------------------\")\n",
    "print(\"Naive estimator: {:.3f}\".format(estimate_conditional_expectation(ds, \"x\", \"y\")[0]))\n",
    "print(\"AB test: {:.3f}\".format(ab_test(fd2, \"x\", \"y\", n=n)[0]))\n",
    "print(\"Linear test: {:.3f}\".format(estimate_ate_frontdoor_linear(ds, \"x\", \"y\", \"z\")))\n",
    "print(\"Non-linear: {:.3f}\".format(\n",
    "    estimate_intervention_frontdoor_nonlinear(\n",
    "        ds, \"x\", \"y\", [\"z\"], model_y_zx=RandomForestRegressor())))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The actual ATE is zero, which our estimate predicts correctly, but the linear ATE model predict incorrectly.\n",
    "\n",
    "This looks good - but it important to realise that there are limitations. Let's go back to the second `StructuralCausalModel` we created: the linear system and see how our estimator performs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Different estimates of the ATE for a linear system\n",
      "--------------------------------------------------\n",
      "Naive estimator: 0.050\n",
      "AB test: 8.092\n",
      "Linear estimator: 8.054\n",
      "Non-linear estimator: 1.424\n"
     ]
    }
   ],
   "source": [
    "# sample from the model\n",
    "n = int(10e3)\n",
    "ds = fd1.sample(n)\n",
    "\n",
    "print(\"Different estimates of the ATE for a linear system\")\n",
    "print(\"--------------------------------------------------\")\n",
    "print(\"Naive estimator: {:.3f}\".format(estimate_conditional_expectation(ds, \"x\", \"y\")[0]))\n",
    "print(\"AB test: {:.3f}\".format(ab_test(fd1, \"x\", \"y\", n=n)[0]))\n",
    "print(\"Linear estimator: {:.3f}\".format(estimate_ate_frontdoor_linear(ds, \"x\", \"y\", \"z\")))\n",
    "print(\"Non-linear estimator: {:.3f}\".format(\n",
    "    estimate_intervention_frontdoor_nonlinear(\n",
    "        ds, \"x\", \"y\", [\"z\"], model_y_zx=RandomForestRegressor())))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It doesn't do well. The reason for this is that we have violated one of the assumptions of the model we're using: good overlap of the covariates. We are using a [random forest](https://en.wikipedia.org/wiki/Random_forest) model to estimate the conditional expectation. This kind of model will do well when we have samples that have good overlap, but fail badly if they attempt to extrapolate beyond that range. \n",
    "\n",
    "Conversely, if a system is actually linear, a well fitted linear model will be able to extrapolate well beyond the range it was trained on. This is why the results of the linear ATE are so good in this case. \n",
    "\n",
    "If we plot the value of our $Z$ variable for the linear model, the lack of overlap becomes clear:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl43Fed5/v3txaptJT23XJseYm3\nxNmcQCCBQAKkQ0iavgxtGC6hm2eAO50JcB9mQg+56R6YuX2hL72kgblkCBnIENJNILShA4FsdAjE\niRMv8SZZ3iVZUkm29l117h+nylYUySpJVfVb6vt6Hj0lqX761dcl+VOnzu8sYoxBKaWUvwScLkAp\npVT6abgrpZQPabgrpZQPabgrpZQPabgrpZQPabgrpZQPabgrpZQPabgrpZQPabgrpZQPhZx64Kqq\nKrN69WqnHl4ppTzp1Vdf7THGVC90nGPhvnr1anbt2uXUwyullCeJyMlUjtNuGaWU8iENd6WU8iEN\nd6WU8iHH+tyVUmo5JicnaWtrY2xszOlSMiISidDY2Eg4HF7Sz2u4K6U8qa2tjWg0yurVqxERp8tJ\nK2MMvb29tLW10dTUtKRzaLeMUsqTxsbGqKys9F2wA4gIlZWVy3pXouGulPIsPwZ70nL/bRruSinl\nQ9rnrpTyhUd3nkrr+T76lkvSer6ZvvGNb/B3f/d3HD16lFgsRlVVVdofQ8M9l+16+MLn2/7EuTqU\nJ8wMz0wGXy54+9vfzu23385NN92UscfQbhmllFqiV155ha1btzI2Nsbw8DBbtmxh//79C/7cVVdd\nRabX1tKWu1JKLdG1117LHXfcwX333cfo6Cgf+9jHWLVqFVdeeeWcxz/66KNs3rw5K7VpuCul1DLc\nf//9XHvttUQiER544AGCwSB79uxxuiwNd6WUWo7e3l6GhoaYnJxkbGyMeDzOjTfeOOex2nJXSimP\n+PSnP81XvvIVjh8/zr333ss3vvENbbkrpVS6ODGC5/vf/z7hcJiPfvSjTE9P87a3vY1nn32Wd7/7\n3Rf9uQceeICvfe1rdHZ2snXrVm677Ta+853vpLU2DXellFqij3/843z84x8HIBgMsnPnzpR+7p57\n7uGee+7JZGk6FFIppfxIw10ppXxIw10ppXxIw10ppXxIw10ppXwopXAXkVtFpFlEWkXkixc57n8T\nESMi29JXolJKqcVacCikiASBbwLvAdqAV0RkhzHm4KzjosBngdTGAimlVDrNXOU0HTK4Uurx48fZ\nvn07vb29XHPNNTzyyCPk5eWl9TFSablfB7QaY44ZYyaAx4A75zjuK8BXAX/uVptLdj184UMplXb3\n3nsvn//852ltbaW8vJyHHnoo7Y+RSrivAE7P+Lot8b3zRORqYKUx5l8udiIR+ZSI7BKRXbFYbNHF\nqgwaH3K6AqU8ZylL/hpjePbZZ/nQhz4EwF133cVPf/rTtNe27BmqIhIA/gb4xELHGmMeBB4E2LZt\nm1nuY6s0iB2Gff8IP/8cNL0D3n2/0xUp5RlLWfK3pqaGsrIyQiEbv42NjbS3t6e9tlTCvR1YOePr\nxsT3kqLAZcDziQ1d64AdInKHMWZXugpVGTA1Bnt+CKE8ePvnYO8P4aFb4JLrYePtkFfkdIXKJV46\n1st3XzzO+OQ0H7pm5cI/kEMWu+RvT09PVupKJdxfAdaLSBM21LcDH03eaYzpB85vACgizwNf0GD3\ngJanYHwAtn0W3vNf4B3/EZ7/K3jpW9D5Olx/t27FpzjTP8qnH3kVgMnpOA/99hgfuW4lNSURhytz\nh8Uu+btp0yb6+vqYmpoiFArR1tbGihUr5jx+ORYMd2PMlIjcDTwFBIHvGmMOiMiXgV3GmB1pr0pl\n3vQUnN4JDVdC+Wr7vfxieN9/g3AB7Pw2vPIduOHz2oLPcV96Yj+T03E+8861TE7H+eZzrXznt8f5\nz7dtcro0V1jKkr/vete7ePzxx9m+fTvf+973uPPOucaoLE9Kfe7GmCeBJ2d9b87OWWPMTcsvS2Xc\n6Z0wOQJ1V7z5vpIVtpX+u2/Y1v1lf5T9+pQrNHcO8uzhbv7j+zZQXmiH6m1pKOWHO09xz83rKc53\n0cKyDryzXOqSv1/96lfZvn079913H1dddRWf/OQn016bi34zKquan4RAEGo22q9nD3ssb4L6K6Ht\nZdj4fgjlZ79G5biHXzxOJBzgo9ddwi/2dwJww7oqXm/v5yevtfHx61c7W6DDlrrk75o1a3j55Zcz\nWZouP5CTjIHD/wKV6yF0kX7T1W+3F13bX8tebco1BsYmeWJ3Ox+8qpHyogsTbFZWFLK+pphfJsJe\nuZOGey7qPw3njkPNAns5ljdBtB7adNJxLnpqfyfjU3E+vK3xTffdsrmWncfP0j8y6UBlKhUa7rko\n2RJPXkidjwjUXQ7nTsLkaMbLUu6yY28Hl1QUcuXKsjfdZwxMxw3Pt3Q7UNnMOvw7XWa5/zYN91zU\nsRsCYYg2LHxs5XrAwNmjGS9LuUdscJwXW3u444oGEvNX3qCxvIDi/BDPHHIu3CORCL29vb4MeGMM\nvb29RCJLH26qF1RzUcduqN0MwRR+/eWr7QtBz5GMl6Wc9ejOU+c/DwjEDdx2ef2cxwZEWFtdxO+P\n2XCd6wUg0xobG2lra8OvS5lEIhEaG9/cJZYqDfdcE49Dx57UhzcGw1DRBD0tma1LucrTh7pYUVbA\npvrovMc0VRWzt62f4z3DrKkuzmJ1VjgcpqmpKeuP6xUa7rnm3HEY74eGq8DEU/uZqvV2dM3IWSis\nyGx9ynETU3Geb45xbVMFP3z59LzHramyk9v+4ZlWrm2q4KNvuSRbJaoUaJ97runYbW8brkr9Z8pW\n2dszF591p/zhaGyIqbhhU13JRY+rLM4jmh/iWI+uKOpGGu65Zvcjtg/99CImUJQk+v3O7M1MTcpV\nWroGyQsGWF1VeNHjRITVVUWc6B3JUmVqMTTcc03fKShpsLNTU5VXCIWVtq9e+V5r9xBNVUWEAgvH\nw8qKQvpHJxkc0/HubqPhnkvi09DfBmVL6BstbdSWew44NzJB7/AE62pSu0C6srwAgLZzOg/CbTTc\nc0nPEZiegLIlrMdd2mgvxo72pb8u5RpHu23/+doUw72+tICAQNs57ZpxGw33XJK8mFq6lJZ74gVB\nW+++djQ2RHF+iNpoagvF5YUC1EQj2nJ3IQ33XNLxGgTzoLhm8T9bkthMoOtAemtSrnKyd4TVVUWL\nmpTUWF5A27lRX84U9TIN91zSsdt2r8gSfu35UXtRtftg+utSrtA3MkHf6CSrKy8+Sma2xvJCRien\nOXVWu2bcRMM9V0xP2q3zlnIxNalms91QW/nSyUQ4r6pc3M5bjYmLqnvb+tNek1o6Dfdc0X3Irs1e\nuozNjas3QqzZLgmofOdk7zB5oQB1i9wbtbYkQigg7D2tF9vdRMM9VyznYmrS6Dm7ofZv/zY9NSlX\nOdU7wsryAoKBxS0CFgwIDWUF7GvTcHcTDfdc0bEb8kuhqGrp54gmVggcPJOempRrjE9N0zUwTmP5\n4vrbkxrLC3i9vZ+p6RTXK1IZp+GeKzp2Q8OVdgOOpYrW2dtB3V7Nb1o6h5g2hoaygiX9fGN5AWOT\ncY506zozbqHhngumxu0QxsUsFjaXvCLIL9Fw96H9HfZi6IolhnvyReFAx0DaalLLo+GeC7r2Q3xy\n+eEOUFwLw85urabS7/X2fiLhAOWF4SX9fFVxPpFwgIMa7q6h4Z4LkhdTV1y9/HMVVcOQhrvfHGjv\np6G0YMk7KgVE2FAb5dAZDXe30HDPBR277QSk5QyDTCquhskRu3GH8oXJ6TiHOgeX3CWTtLmhhEOd\nAzpT1SU03HNB+27bJZOOfS6LEksX6J6qvnGka4iJqfiSL6YmbaovoW9kks6BsTRVppZDw93vJkYg\ndgga0tAlAxfCvbc1PedTjkteTE1HuAPaNeMSGu5+1/m63Su14cr0nK+wwq5N06std7/Y395PUV6Q\nyuK8ZZ1nY53dTPvQmcF0lKWWSTfI9rvOffY21pyeC6GBIBRWacvdR/a397OloZTAMrvtopEwKysK\nOKgtd1fQlrvfndkLBRUQKUvfOYtroEfD3Q+m44aDZwbYsuLim2GnanN9iXbLuISGu5/tehiOPW+X\nHEjHxdSkomo4e8xu26c87VhsiLHJOJc1lKblfJvqSzjeM8z/fPEEj+48lZZzqqXRcPez+DQMdlzY\naCNdimtgehz6T6f3vCrrDnfa/vHkxdDl2lRfgjHQpSNmHKfh7mdDXTbgSxvTe14dMeMbLV2DBAPC\nmurFreE+n82JF4kz/RruTtNw97OBNnubiZY7aL+7DzR3DrK6spBIOJiW8zWWFxDND3GmX/dUdZqG\nu5/1t0MgvLQ9Uy8mr9guIKYtd89r6RpkQ2IIYzqICBvrozqRyQV0KKSf9bdBScPS9ky9GBGoXKdj\n3T1uZGKKk70jrK0uTsvFz+Q5AiJ0DYzpMgQO05a7XxkDA+1QmuYumSQRO8xy18OZOb/KuNbuIQx2\nm7x0qi2JMDYZp390Mq3nVYuj4e5XfSftnqnp7m9PKqqx2+5NT2Tm/CrjmhMjZRa7Z+pCkufTETPO\nSincReRWEWkWkVYR+eIc939GRF4XkT0i8lsR2Zz+UtWinEnMTC1J80iZpOJqezscy8z5Vca1dA0S\nCggVy1x2YLbkO4HOgfG0nlctzoLhLiJB4JvAHwCbgY/MEd6PGmMuN8ZcCXwN+Ju0V6oWp/N1QKCk\nPjPnTw6HHNJw96rmriFqovnLXnZgtoK8IKUFYW25OyyVlvt1QKsx5pgxZgJ4DLhz5gHGmJnzjYsA\nvZLitM59dtekYHpbZeclN9oe0XD3qubOgbT3tyfVlUTo1LHujkol3FcAM6citiW+9wYi8mcichTb\ncr8nPeWpJes6YEfKZEooYodDareMJ/WNTNA1MJ6xcK8tiRAbHGdyOp6R86uFpe2CqjHmm8aYtcC9\nwH1zHSMinxKRXSKyKxbTUMiYsQG7NEA0Q10ySUVVMNyT2cdQGdHSNQSkf6RMUl1pPtPGcLxnOCPn\nVwtLJdzbgZn7szUmvjefx4A/nOsOY8yDxphtxpht1dXVqVepFifWbG+jdZl9nKJqbbl7VHOXHSlT\nW5KfkfMnXzSSa9eo7Esl3F8B1otIk4jkAduBHTMPEJH1M758P6CzW5wUO2RvM95yr4bxQftOQXlK\nS+cg0fwQpQXhjJy/OppPQGy/vnLGgjNUjTFTInI38BQQBL5rjDkgIl8GdhljdgB3i8gtwCRwDrgr\nk0WrBXQfhlCB3TUpk4oS777OHrV7tCrPaO4c5NK6KJLmkTJJoUCAquL882PpVfaltPyAMeZJ4MlZ\n37t/xuefTXNdajlih6D60vQvOzBbMtx7Ndy9xBjD6+39XL4iPWu4z6euNKLdMg7SGap+1H0Yqjdl\n/nGSwyHPHsv8Y6m06R4cZ3RyOmP97Ul1JRHazo0yOKbLEDhBw91vRvvsBh01GzP/WME8u31f79HM\nP5ZKm2RXSaZGyiQlz9/Spa13J2i4+80LX7e3/Rcb0JRGRdW2z115RjLcazIc7nU6YsZRGu5+M9hp\nbzM9DDKpqFpb7h7T3DVIcX6I4vzMrvhdVhimOD+kF1UdouHuN4OdtrukoDw7j1dUDaNnYeRsdh5P\nLVtL12DaV4Kci4hwaW2xttwdouHuN4NnoLgu8yNlks4Ph9SLql4QjxtaugYzfjE1aUNdCc2dg7px\nhwM03P1mqBNKstQlAxdGzGjXjCecPjfC2GQ84xdTkzbWRekfnaRLl//NOg13Pxk5a2eMFmd4ZupM\nhVX2XYJeVPWEw1kaKZOU3J/1sM5UzToNdz/pTi47kMWWezAEpSu15e4RLcmRMtHsdMtsTIS7XlTN\nPg13P4k5EO4AlWuhtzW7j6mWpLlrkJUVBeSHg1l5vLLCPGpLdBkCJ2i4+0msBYL5dmJRNlVtgJ4j\ndlNu5WotXYNsqI1m9TE31JXoiBkHaLj7SU8LFNdAhhaDmlf1pTA5DP1t2X1ctSgTU3GOxYZZn+Vw\n31gXpTU2xJRu3JFVGu5+0tNit9bLturEUgfJdeSVKx3vGWYqbs73g2fLhtooE1NxTvTqxh3ZpOHu\nF+ODMNDucLgfzv5jq5QlN+hYX5PtbpnkiBntmskmDXe/6Ensj1Jck/3HLqywk5k03F3tSNcgwYCw\nprooq4+7rqaYUEA42KHDIbMps4tLqOw5H+4OtNzBtt61W8bVmjsHWV1ZSCRLI2UAHt15CoCq4nwO\naLhnlbbc/aKnGSR4YcZotlVvsOGuI2Zcq6Vr8HwXSbY1lBWwv71flyHIIg13v+hpgYo1EHDgzdiu\nh2H0HIz3X1iVUrnK6MQ0J3vt0gPJ1nQ2NZRF6B2e0GUIskjD3S9iLVB1qXOPX5yYOKX97q50NDaE\nIXvLDsy2oqwAgP3t/Y48fi7ScPeD6Um7tku1g+GenBWr/e6udGH3pewsOzBbXWkEEbTfPYs03P3g\n3AmIT9mZok7JK4ZwobbcXaolMVKmssiZcM8PBWmqKmJ/h7bcs0XD3Q+SrWUnu2VEbOv9+G9sH7xy\nlZauQaqL8wkGsjx7eYbLGko5oN0yWaPh7gc9Lfa2ar2zdRTX2guqOiLCdVq6hhzrkkmanI7T0T/G\n//hX3dglGzTc/aCnBaL1EClxto7iOpgcgYkhZ+tQbzA4Nkl736hjF1OTGhIXVTv6Rx2tI1douPtB\nj8MjZZKSF1V1OKSrtHTZF1unw72+1D5+R9+Yo3XkCg13rzPG+WGQScnZsUMa7m5ypCu7uy/NpzAv\nRHlhmI4+bblng4a71w12wsSgnSHqtEgphAvsJt3KNZq7BikIBykrDDtdCg1lBRruWaLh7nU9yZEy\nDl9MhcSImXrtlnGZlq5BLq0tJpDtdf7n0FBWQO/wBINjk06X4nsa7l4XS46UcUHLHWy4D3ToiBkX\nae4cyvoGHfNpKLUXVXWFyMzTcPe6nhYI5UPzL9wxvjxaD1NjNuCV484OT9AzNJ71rfXm01Bm+/33\na7hnnIa71/U02wuZLnjLDdhwB+g+6GwdCoBvPW83Lm93ST93NBImGgnpZKYs0PXcvSrZSu/Y6+ya\nMrMlw73rAKx/j7O1qPOrMDo9UmamhtICXWMmC7Tl7mWTo3aZXac26JhLXqEdNdN9yOlKFNA9MEYk\nHKAk4p52XENZAUe6BxmdmHa6FF/TcPeyoW5768TWehcTrYfuA05XoYCugTFqoxHELd122H73uIHD\nndp6zyQNdy8bToR7kQvDPdYC01NOV5LTjDF0DYy7qksGLoyY0a6ZzNJw97LhGCDOba03n2g9TI/D\nWV0gykndg+OMTk47vmDYbGWFYUoLwhruGabh7mVD3VBY6czWehdToiNm3KAlsexAjcta7iLCloYS\nDura7hml4e5lwzEoqna6ijcrrgUJaLg77MLuS+4Kd4AtDSUc6hxkcjrudCm+peHuVcbYcC92YbgH\n8+xm3RrujmrpGqQoP0Rxvsve2QFbGkqZmIpzNKbLQ2dKSuEuIreKSLOItIrIF+e4//8UkYMisk9E\nnhGRVekvVb3BWD9MT7jvYmpSuBBOveSOWbM56kj3ELVRd/W3J21psHsPHGjXfvdMWTDcRSQIfBP4\nA2Az8BER2TzrsN3ANmPMVuBx4GvpLlTNMhyzt27slgF7UXW4x74Aqax6dOcpfvDSSQ6dGaDapeG+\n8/hZwkHhJ6+1OV2Kb6XScr8OaDXGHDPGTACPAXfOPMAY85wxZiTx5UtAY3rLVG+SDHe3jXFPitYB\n5sJYfJVVQ+NTjE3GXRvuARHqSiK068YdGZNKuK8ATs/4ui3xvfl8EvjFcopSKRiOQSBsZ4O6UXFi\nV6ahLmfryFGxQbvsgFvDHexM1TP9o8TjuoJoJqT1gqqIfAzYBvz1PPd/SkR2iciuWCyWzofOPcPd\ndny7uPSaeFG1rU3XdndEbCgR7sUuDvfSAsan4pw+N7LwwWrRUkmGdmDljK8bE997AxG5BfgScIcx\nZnyuExljHjTGbDPGbKuudmlfsVcMuXQYZFIwZF98tOXuiO7BcfKCAUoLnN99aT7JDbN1MlNmpBLu\nrwDrRaRJRPKA7cCOmQeIyFXAt7HBrp2smTY9BSM97g53sF0zuuWeI3oGx6mO5rtqTZnZakvyCQjs\n1+V/M2LBcDfGTAF3A08Bh4B/MsYcEJEvi8gdicP+GigGfiQie0RkxzynU+nQdxJM3L0XU5OidXbE\nzNScb+RUBsUS4e5moWCAmmhEW+4ZktLsBmPMk8CTs753/4zPb0lzXepieo/aW9e33GsBA72tULvF\n6WpyxsRUnL7RSapc3N+e1FAW4UBHP8YYV7/L8CKXXo1TF9Vrd9dx7QSmpGhixEzssLN15JiexMXU\nGpe33MH2u/cMTdA9qO/u0k3D3Yt6WyFcAHlFTldycUU1gECs2elKckq3B4ZBJtWfX/5X+93TTcPd\ni3pbE0MNXf42Nhi2I2Z0V6asig2OI0BlUZ7TpSyovtQuaqbLEKSfhrsX9R51f5dMUnGtttyzLDY0\nTkVRHqGg+/97R8JBKovy+OUBnQ+Rbu7/7as3mhyFgTb3X0xNitbB2aMwpWvMZEuPB0bKzFRbEjm/\nkbdKHw13r0nubuTGpX7nUlwH8SndlSlLpuOGniHvhXvv0Dhjk7phdjppuHvN+ZEyHgl3HTGTVW3n\nRpiKG1cvOzBbXWkEA7R269ru6aTh7jVeC/fi5IgZDfdsSG5+4a2Wu631cGLnKJUeGu5e03vUdnWE\n3Ld12pyCeVC+SsM9S452DwPeCvfKonxCAaG5U0fMpJOGu9f0tkLlOqerWJzqjTpiJkuOxoYoyg9R\nmOe+rfXmEwwINdF8mru0WyadNNy9prcVKtc6XcXiVG+AniN2wTOVUa3dQ57qb0+qLYloyz3NNNy9\nZOQsjPR6sOW+CeKTOmImC47GhjzVJZNUV2qHQ/aN6JDZdNFw95JkOHou3DfYW+13z6izwxOcG5n0\nZLjXlthrSHpRNX003L0kOVLGa+Fedam91X73jDo/Usaj3TIAzRruaaPh7iX7fwwIHP+N05UsTn4x\nlF6iLfcMS44T92LLvSQSorQgrC33NNJw95LhGBRWQsA7IyHOq96gLfcMa+kapCAcpKzQvVvrzUdE\n2FAXpaVLwz1dNNy9xO37pl5MzUboaYG4TjHPlNbuIdbVFBNw+2qh80huufeDl046XYovaLh7hTEw\n3O3dcK/eCNPjcO6E05X41pGuIdbXFjtdxpLVlkQYn4rTPzrpdCm+oOHuFYOdMD3hnQXDZqveaG+1\n3z0jBsYm6RwYY31N1OlSlqwmai+q6gqR6aHh7hVe2VpvPudHzGi4Z0LyYmpH36jDlSxdco2Z7sEx\nhyvxBw13rzgf7lXO1rFUkRIoWaEXVTOkNTF13wv7ps6nMC9ENBLSlnuaaLh7Rc8RCIShoNzpSpau\neqNuuZchR7oHCQWEcg9srXcxtdEIXQPack8HDXev6Gm2y+eKh39l1TpiJlNauuyyA14dKZNUU5JP\n9+AY8bhxuhTP8+CA6RzV05JYG92Ddj1sb0d6YGoM+k5CxRpna/KZ1u6h87M8vaw2GmFy2tDeN8rK\nikKny/E0DzcDc8jECPSdtptNe1m03t52HXC2Dp8ZGp+ivW/U0/3tSTWJi6o6mWn5NNy9oPcIYOwm\nHV5W0mC7lTr3O12Jrxzt9v7F1KTkcMgWXdt92TTcvaDniL31ardMUjAPKtZC5+tOV+IrR86Hu/e7\nZQrygpQWhLXlngYa7l4Qa7YtXq/OTp2p7jLo0nBPpyPdg+QFA54fKZNUE83XcE8DDXcv6GmBslUQ\n9N6CUG9Sdzn0nYKxfqcr8Y0jXUOsqS4iGPD2SJmk2pIIrd1DTOuImWXRcPeCnpYLG1543UCnvf3N\nXztbh48c7BhgU32J02WkTU00n/GpOKfPjjhdiqdpuLtdfNrOTk1O3/e6kgZ7O9DubB0+0Ts0TufA\nGJt9FO7JIZ3aNbM8Gu5ud+6EXTDML+EeKYVwEQx0OF2JLxw8YzeV3tLgn3BPjvpJXihWS6Ph7nY9\nLfbWL90yIlC6QlvuaXKww4a7n7pl8hMbjjx1oJNHd55yuhzP0nB3u2S4V613to50KmlILGE85XQl\nnnfwzAANpRHfjJRJqo1G6NYFxJZFw93tYi12mV8vLxg2W0kDxCcvrHSpluxAxwCbG0qdLiPtakry\niQ2N64iZZdBwdzs/jZRJKmm0t106U3U5hsanOBob4rIV/umSSaqNRpiOG84OTzhdimdpuLuZMXbC\nj5m+sPiWHxTXgAR1puoy7Wvrwxi4cmWZ06WkXXKNGV3+d+k03N1s8AxMjnp/TZnZAiGI1mq4L9Pe\n03Yi2BWNPgz3aAQBunRXpiXTcHez5OqJJSucrSMTSlfCmT323Ylakj2nz7G6stB3F1MB8kJ2OYWu\nfg33pUop3EXkVhFpFpFWEfniHPe/Q0ReE5EpEflQ+svMUcmWbUm9s3VkQulKGOm1SxGoJdl7up8r\nfNglk1RXEuGMhvuSLbhZh4gEgW8C7wHagFdEZIcx5uCMw04BnwC+kIkic1bXfjtKJuzDTQvKLrG3\nHbuhfJWztXjMoztP0T86SefAGFPTxrdjwetLIxw6M8DIxBSFebqv0GKl0nK/Dmg1xhwzxkwAjwF3\nzjzAGHPCGLMPiGegxtzVud+fXTIA0Qa7J2zHbqcr8aQTvcMArKr04Qt/Qn1pBAMc7tRlCJYilXBf\nAZye8XVb4nuLJiKfEpFdIrIrFost5RS5Y3LUbtKRXIvFb4IhqN2i4b5EJ3uHyQsGqC8tcLqUjKlL\n/NsOJZZYUIuT1QuqxpgHjTHbjDHbqqt9sDZ5JnUfAhO3LVy/argKOvZAXN/wLdaJnhEuqSj0zTK/\ncykvDJMfCmi4L1Eq4d4OrJzxdWPieyqTOvfZ21KfdssArLgaxvvh3HGnK/GU0YlpugbGWFXl3y4Z\nABGhrjRyfv0ctTiphPsrwHoRaRKRPGA7sCOzZSnaX4NIGRRWOV1J5jRcZW+1a2ZRTp4dxgCrK4uc\nLiXjVpQVcPDMAFPT+u5usRbblLhSAAAS70lEQVQMd2PMFHA38BRwCPgnY8wBEfmyiNwBICLXikgb\n8G+Ab4uIbm+/XO2v2Zat+PdtN9UbIRSx/1aVsmOxYUIB4ZIKf7fcARrLCxibjOvyv0uQUp+7MeZJ\nY8ylxpi1xpj/lvje/caYHYnPXzHGNBpjiowxlcaYLZks2vcmRqD7IDRc7XQlmbX7f0FxLbT8wulK\nPOVobIhLKgoJB/0/B7Gx3L6A7Wvrc7gS7/H/X4cXdSbWk1lxjdOVZF7ZJdDfZnecUgs6OzzBmf4x\n1lQXO11KVlQW5VESCbG3TffcXSwNdzd66b/b2+Ra7n5WttLuNJUL/9Y0eOlYLwDrqv3f3w72ourW\nxjJtuS+Bhrsb9Z2wF1Mj/lun+03KErNTT7/sbB0e8cKRGPmhACvK/d/fnnTFylIOnxlkdELf3S2G\nhrvbGANnj0FFk9OVZEdRDeQVw8kXna7E9eJxwzOHullfG/X1+PbZrl1dwVTc8Nqpc06X4ika7m7T\ndxLG+qFirdOVZIcIVK6D4y/oCpELONAxQPfgOBvrok6XklXXrConILAz0SWlUqPh7jYnf2dvcyXc\nwe4PO9hh37GoeT1zuAsRuLQ2t8I9Gglz2YpSXjp+1ulSPEXD3W1OvmhXgYzWOl1J9lSus7cnXnC2\nDpd79nA3V60sozg/91ZIfEtTBXtO9zE2qf3uqdJwd5uTv7P97ZJDv5qiGrvb1NHnnK7EtboHxtjX\n1s/Nm3LoRX+G69dWMjEV55UT2npPVQ4liAecO2m7JirXO11JdonAhluh9WmYGne6Gld6rrkbgHdv\nrHG4Emdcv6aKSDjA0we7nC7FMzTc3eToM/a2eqOzdThh4+0wMQTH/9XpSlzpmUPdNJRGcu5iKtjN\nSZ7Y3c4N66p5+lA3Ri+8p0TD3U1an4GSRjslP9c0vcMOiTz8c6crcZ2RiSleONLDzZtqET+vNbSA\n926upb1vlIO6BHBKNNzdYnrStlrXvdvfi4XNJ5QP698Dh35unwt13tOHuhmdnOb2rT7cS3cR3r2p\nhmBA+NneM06X4gka7m5x+mUYH4C1NztdiXO2/jGM9MCRXzldiav8bG8HtSX5HOke8u1+qamoKs7n\nXRtq+PFrbboEcAo03N3i0A4I5sO6W5yuxDnr3mO7pHb/L6crcY3+kUl+0xzj9q0NBHLxHd0sH97W\nSGxwnP/ys4M8uvNUTr/YLUTD3Q2MgUM/g3U3Q35urPY3p2AIrtgOLU/BYKfT1bjCE7vbmJiO88Gr\nfLwj1yK8a2MN1dF8fne0x+lSXE/D3Q3aX4OBdigog10PO12Ns66+y+4du+OenH8ujDE89sppLl9R\nymUrcmARuRSEgwH+3Y1NHI0Nc+rsiNPluJqGuxsc+AlIEGouc7oS51WuhZrNdqZujl9Y3dvWz+HO\nQdZUF2n3wwz/9i2rKAgHeT4x9l/NLffmMbvN1ATsfQxqL4O83FnG9U1mttLXvBNe+hZ05Pb2e4+9\nfIpwULiisczpUlxh5gvc29dV8fShLjr6Rh2syN205e60I0/ZESIr3+J0Je5RuR6i9XDsNzm7UuTQ\n+BQ79nawtbGMSDjodDmu87a1lUTCgfMzd9Wbabg77dXv2SCr3uB0Je4hAk3vtCtF5uhiYj/b28HI\nxDTXrip3uhRXioSDXL+mkgMdAxzs0ElNc9Fwd1LXQWj9NVzzJxDQ1tkbrLga8ooubDmYQ4wxPPzi\ncTbVl7CyIoe76hZww7pqIuEAf/PrZqdLcSUNdyft+A8QzLMhpt4omAer3g7Nv8i5dd5fONJDS9cQ\nn7yhKaeXG1hIQV6QG9fb9WZ0l6Y303B3Ss8Re8Hwkus13Oez6u0QCMHOB52uJKse+u1xqorz+cAV\nub3cQCretraSyqI8vv4rbb3PpuHulF/+uW2d5vKM1IVESqH+Cnj1u/D73Oieae0e5DctMa5cWcqP\nX213uhzXyw8F+T9uWsuLrb06sWkWDXcnHPxn29d+6fsgP/eWcF2UpnfaNd5Pv+R0JVnx0G9PEAoI\n1zVVOl2KZ3zsrauoK4nw/z7VrMsBz6Dhnm3nTsA//weovxJW3+h0Ne5XthIq1thRM3F/b7HWPTDG\nj19r46pLcnMrvaX6yWvtvGVNBa+d6uMvdhzQCV8JGu7ZNDUOP/qE/fzD37P9yWphTe+EkV44/C9O\nV5JRn//HPUxNx3nH+mqnS/GcbasqqCjK49cHu4hr6x3QcM+uX90HHbvhD78F5audrsY7ai+Dwip4\n4eu+ndQUGxzn5RNnuXJlGZXF+U6X4znBgHDzxhrO9I9xQMe9Axru2XPgCXj5QXjrn8FwLOcXxVqU\nQNCumHlmz4WtCH3mOy8cY2racNOlublHajpcsbKMmmg+Tx/qYjruz0bAYmi4Z0PvUXji01C2Ciqa\nnK7GmxqvhZIV8K9fd7qStOsdGueRl05yxcoyqqLaal+qgAi3bKolNjjOT3frSCMN90ybHIMf3WVX\nfbz6Lu1nX6pAyK6/c+p38NR9TleTVn/39BHGp+LctEH72pdrS0MJDWUR/vbpFsYm/X0BfiEa7plk\nDDz5Beh8Ha78t1BY4XRF3nbJW+0m2j7ahu9w5wA/2HmSj73lEmqiEafL8TwR4Q8uq6ft3Cjf/k1u\nzWyeTcM9k3Z+G3Y/Ajd+AWq3OF2N9wXz7B6zPc1w5Gmnq1m2yek4/+nxfZQWhPncLZc6XY5vrK0u\n5gNXNPCt51s5GhtyuhzHaLhnytFn4ak/h423w7u+5HQ1/rH6Riiqhl/ea9fC97B/eOYI+9r6+a9/\neDnlRXlOl+Mr/9f7N1GYF+Rzj+1hYio3N9PWcM+EM3vtePbiWhtGr33P6Yr8IxiCLR+E3lZ49itO\nV7NkO/Z28MCzrXzomkbev1XXkEm3mpIIf/VHW3m9vZ8vPfF6Ts5c1at76XbqJXj0w5BfAtd8AkI6\n+iHtajbDtj+F3z1g++E3vt/pihbl8VfbuPfH+7h2dTlbV5TqjMoMufWyOj5783r+/pkjRMJB/uID\nmwkFc6c9mzv/0kyLx+H334KHb4NA2AZ7oa4PkjHv+7+h4Sr4p7vsHAIP6Ogb5XOP7eYLP9pLU2UR\nt11Wn1Nhk02P7jzFoztPURPN51PvWMMjL53kjx98if3t/U6XljXack+HEy/C038Bba9AzRY7MiaX\n90PNhnAB/O8/hR98yHaB7f8xXPvvbEveJe+WjDGc7B3h5eNnefpQF881d2MMvHtjDTdtqCYU0GDP\nNBHhP9+2ic31JXz55we5/R9+y3VNFXzsrat47+ZaX29hKKn0RYnIrcDfA0HgO8aY/2fW/fnA94Fr\ngF7gj40xJy52zm3btpldu3YtsWwXiMftYla//Vs49pzthtl0O6y41m4Tp7JjegqOPWtnrk6N23dN\nNRvtUsG1l0Pd5XakUkF2Npnu7B/ja788zJHuIU70DjM4NgVAbUk+H9jawJ/c0MRvmmNZqUW90ejE\nNLtOnmXn8bOcHZ4gLxRgU12UP3vXOt5xabVngl5EXjXGbFvwuIXCXUSCQAvwHqANeAX4iDHm4Ixj\n/j2w1RjzGRHZDnzQGPPHFzuvJ8N9espeLG19Gvb8APpO2q6X1TfYjSWCOuLBMdMTUNIA7a/CmX32\n9zQyY33vslW2G2fF1dBwNTRcmZbllofHp3j5xFl+19rDC0d6ONw5CEBJJERTVRGrq4poqiyiOpqv\nuyq5RNwYjsWGeb29jwMdA4xMTFOcH+I9m2u57fJ63ra2kiIXr8qZznC/HvhLY8z7El//OYAx5q9m\nHPNU4pjfi0gI6ASqzUVOnpZwN8YuAxufhOlJiE8lbhOfJ5eIDYQSH8ELtxK0LeyJEZgYSnwMJz5m\nfD7aZ0P87HHo3GfvA6hcb2dM1m/VUHcjY2B8AAba7Ud/G/SdgtHkdmwClWvtAm5lq6DsEvviEC60\nO2PlFUEownQgzBQhJgjRNwaxMcOJs5Mc7B5jb8cQr57qI24gLxjgmlXl3LShmuHxaWpLNMy9YDpu\nOBYb4vX2fg50DDA6OU0oIGxtLGVrYxkrKwpZUVZAUX6QgnCQgrwg4cR1EiH5Jl0QgaAI4VCAcEAI\nBwOEgvY2GBCCIgQC6fl7SDXcU3l5WgGcnvF1G/CW+Y4xxkyJSD9QCaR/a5Tffwue+XIivCfTfvo5\n5ZdCUaV9q1+5zn7oJhvuJmJ3coqU2tE1SeND0H/KBv3AGYgdhpO/g8mROU8TTHzkA1FgJXA18EeJ\n+yfzw8QDYX72tseZiDYAEI2EM/bPUukVDAjra6Osr41y55WG4z3DHIsNMTA2yeOvtjE0PpW2xwoI\nhAI27P/iA5vZft0laTv3XLL63kNEPgV8KvHlkIhkeuPDKtLyAjOAfe3as/xTzS1NdWaF1jqnG5bz\nw/qcpp+r6/zIf4WPXPhysbWuSuWgVMK9HdtgSWpMfG+uY9oS3TKl2Aurb2CMeRDI2m7HIrIrlbcv\nTvNKnaC1ZoJX6gTv1OqVOiFztaYyFusVYL2INIlIHrAd2DHrmB3AXYnPPwQ8e7H+dqWUUpm1YMs9\n0Yd+N/AUtvvxu8aYAyLyZWCXMWYH8BDwiIi0AmexLwBKKaUcklKfuzHmSeDJWd+7f8bnY8C/SW9p\naZG1LqBl8kqdoLVmglfqBO/U6pU6IUO1pjSJSSmllLfo/GellPIhX4W7iPyjiOxJfJwQkTnHLibu\nez1xXNanyYrIX4pI+4xab5vnuFtFpFlEWkXki9muM1HDX4vIYRHZJyJPiMic8/idek4Xeo5EJD/x\nd9EqIjtFZHW2aptVx0oReU5EDorIARH57BzH3CQi/TP+Lu6f61xZqPWiv0uxHkg8p/tE5GqH6tww\n47naIyIDIvK5Wcc49pyKyHdFpFtE9s/4XoWI/FpEjiRuy+f52bsSxxwRkbvmOmZBxhhffgBfB+6f\n574TQJWDtf0l8IUFjgkCR4E1QB6wF9jsQK3vBUKJz78KfNUtz2kqzxHw74H/L/H5duAfHfqd1wNX\nJz6PYpf0mF3rTcDPnahvMb9L4DbgF9hJmm8Fdrqg5iB2ZvwqtzynwDuwc972z/je14AvJj7/4lz/\nn4AK4FjitjzxefliH99XLfcksfO+Pwz80OlaluE6oNUYc8wYMwE8BtyZ7SKMMb8yxiSn6b2Enefg\nFqk8R3cCyd1SHgduFgfWBTDGnDHGvJb4fBA4hJ3Z7UV3At831ktAmYg4vePIzcBRY8xJh+s4zxjz\nr9jRgzPN/Hv8HvCHc/zo+4BfG2POGmPOAb8Gbl3s4/sy3IEbgS5jzJF57jfAr0Tk1cSsWSfcnXhL\n+9153prNteyD02Hwp9gW21yceE5TeY7esDQGkFwawzGJrqGrgJ1z3H29iOwVkV+IiFMb7y70u3Tj\n3+Z25m/MueE5Tao1xpxJfN4J1M5xTFqeX/cufTYPEXkaqJvjri8ZY/458flHuHir/QZjTLuI1AC/\nFpHDiVfZrNQJ/HfgK9j/RF/BdiH9aToffzFSeU5F5EvAFPCDeU6T8efUD0SkGPgx8DljzMCsu1/D\ndisMJa7D/BRYn+0a8djvMjG58g7gz+e42y3P6ZsYY4yIZGy4oufC3Rhzy8XuF7v8wR9h15af7xzt\nidtuEXkC+/Y+rX+8C9WZJCL/A/j5HHelsuxDWqTwnH4CuB242SQ6Bec4R8af0zmkbWmMbBCRMDbY\nf2CM+cns+2eGvTHmSRH5lohUGWOyukZKCr/LrP1tpugPgNeMMV2z73DLczpDl4jUG2POJLqyuuc4\nph17rSCpEXh+sQ/kx26ZW4DDxpi2ue4UkSIRiSY/x14w3D/XsZkyq3/yg/M8firLPmSc2I1a/hNw\nhzFmzqUTHXxOPbM0RqKf/yHgkDHmb+Y5pi55PUBErsP+/8zqC1GKv8sdwMcTo2beCvTP6Gpwwrzv\n1N3wnM4y8+/xLuCf5zjmKeC9IlKe6LJ9b+J7i+PEVeRMfgD/E/jMrO81AE8mPl+DHVWxFziA7XrI\ndo2PAK8D+xK/7PrZdSa+vg07quKoE3UmamjlwpKYe7gw8sQVz+lczxHwZeyLEUAE+FHi3/EysMah\n5/EGbDfcvhnP5W3AZ5J/r8DdiedvL/bi9dscqHPO3+WsOgX4ZuI5fx3Y5sRzmqilCBvWpTO+54rn\nFPuCcwaYxPabfxJ7vecZ4AjwNFCROHYbdpe75M/+aeJvthX4k6U8vs5QVUopH/Jjt4xSSuU8DXel\nlPIhDXellPIhDXellPIhDXellPIhDXellPIhDXellPIhDXelEkTkMzPW/T4uIs85XZNSS6WTmJSa\nJbEGzLPA14wxP3O6HqWWQlvuSr3Z32PXoNFgV57luVUhlcqkxAqYq7BrkijlWdoto1SCiFyD3R3n\nRmN3wFHKs7RbRqkL7sbuW/lc4qLqd5wuSKml0pa7Ukr5kLbclVLKhzTclVLKhzTclVLKhzTclVLK\nhzTclVLKhzTclVLKhzTclVLKhzTclVLKh/5/9Ymr3EN+eyMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7feeed44be80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "z1 = ds.loc[lambda df: df.x == 1].z\n",
    "z0 = ds.loc[lambda df: df.x == 0].z\n",
    "\n",
    "sns.distplot(z1, label=\"x=1\")\n",
    "sns.distplot(z0, label=\"x=0\")\n",
    "\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to realise that this trade-off between ability-to-extrapolate and quality-of-the-model (note: I am not using these terms in a precise way) is nothing specific to front door adjustment: [we saw exactly the same thing in my first post which focused on techniques that applied backdoor adjustment](http://www.degeneratestate.org/posts/2018/Mar/24/causal-inference-with-python-part-1-potential-outcomes/#Covariate-Imbalance). This is not even specific to causal inference - any model building approach will require some trade-off at this level.\n",
    "\n",
    "The main difference between using models to make causal inferences and using them for predictions is that if our goal is just predictions we can use [cross-validation](https://en.wikipedia.org/wiki/Cross-validation_%28statistics%29) to assess the accuracy of our model. When doing causal inference we usually have to rely assumptions about the system - and any conclusions we draw from our models will only be as good as the assumptions we put in. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion\n",
    "\n",
    "In this post I've covered how we can use the front door adjustment approach to estimate the average treatment effect even in situations where we do not have access to the confounding variables of a system. While this is  powerful, it does require other assumptions about our ability to observe mediating variables. \n",
    "\n",
    "At a more general level, I have hoped to demonstrate that there is more to causal inference then just the back door adjustment, and that graphical models are an excellent way to express our knowledge of a system. It should also raise the question: are there more general causal inference approaches than just frontdoor and backdoor adjustment?\n",
    "\n",
    "In turns out there are. The system of rule for deriving such expressions is known as [do-calculus](http://ftp.cs.ucla.edu/pub/stat_ser/r402.pdf) and I plan to write my next post on this topic. \n",
    "\n",
    "A separate but more practical question is once we have these expressions about causal interventions, how do we _actually estimate_ there values? This is a fascinating topic in itself, and one I hope to tackle in another post. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Code\n",
    "\n",
    "You can find the notebook for this post on github [here](https://github.com/ijmbarr/notes-on-causal-inference)."
   ]
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}