diff --git a/notebooks/motivating_example.ipynb b/notebooks/motivating_example.ipynb index 53b2a93..d7ec3f0 100644 --- a/notebooks/motivating_example.ipynb +++ b/notebooks/motivating_example.ipynb @@ -2,224 +2,259 @@ "cells": [ { "cell_type": "markdown", - "source": [ - "# Motivating example: Figure 4" - ], + "id": "abc410d51568d48c", "metadata": { - "collapsed": false - }, - "id": "abc410d51568d48c" - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "initial_id", - "metadata": { - "collapsed": true, - "ExecuteTime": { - "end_time": "2024-01-11T22:28:46.723639500Z", - "start_time": "2024-01-11T22:28:35.137759Z" + "collapsed": false, + "jupyter": { + "outputs_hidden": false } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Collecting git+https://github.com/y0-causal-inference/eliater.git@linear-regression\n", - " Cloning https://github.com/y0-causal-inference/eliater.git (to revision linear-regression) to c:\\users\\pnava\\appdata\\local\\temp\\pip-req-build-f3onojbm\n", - " Resolved https://github.com/y0-causal-inference/eliater.git to commit f666788b42cf32722a3bef39754a9bb19a375a92\n", - " Installing build dependencies: started\n", - " Installing build dependencies: finished with status 'done'\n", - " Getting requirements to build wheel: started\n", - " Getting requirements to build wheel: finished with status 'done'\n", - " Preparing metadata (pyproject.toml): started\n", - " Preparing metadata (pyproject.toml): finished with status 'done'\n", - "Requirement already satisfied: y0>=0.2.5 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from eliater==0.0.1.dev0) (0.2.5)\n", - "Requirement already satisfied: scipy in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from eliater==0.0.1.dev0) (1.11.2)\n", - "Requirement already satisfied: numpy in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from eliater==0.0.1.dev0) (1.25.2)\n", - "Requirement already satisfied: ananke-causal>=0.5.0 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from eliater==0.0.1.dev0) (0.5.0)\n", - "Requirement already satisfied: pgmpy>=0.1.24 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from eliater==0.0.1.dev0) (0.1.24)\n", - "Requirement already satisfied: matplotlib in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from eliater==0.0.1.dev0) (3.8.0)\n", - "Requirement already satisfied: pandas in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from eliater==0.0.1.dev0) (1.5.3)\n", - "Requirement already satisfied: seaborn in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from eliater==0.0.1.dev0) (0.13.0)\n", - "Requirement already satisfied: optimaladj>=0.0.4 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from eliater==0.0.1.dev0) (0.0.4)\n", - "Requirement already satisfied: jax<0.5.0,>=0.4.8 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from ananke-causal>=0.5.0->eliater==0.0.1.dev0) (0.4.21)\n", - "Requirement already satisfied: jaxlib<0.5.0,>=0.4.7 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from ananke-causal>=0.5.0->eliater==0.0.1.dev0) (0.4.21)\n", - "Requirement already satisfied: mystic<0.5.0,>=0.4.0 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from ananke-causal>=0.5.0->eliater==0.0.1.dev0) (0.4.1)\n", - "Requirement already satisfied: statsmodels<0.14.0,>=0.13.2 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from ananke-causal>=0.5.0->eliater==0.0.1.dev0) (0.13.5)\n", - "Requirement already satisfied: contourpy>=1.0.1 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from matplotlib->eliater==0.0.1.dev0) (1.1.1)\n", - "Requirement already satisfied: cycler>=0.10 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from matplotlib->eliater==0.0.1.dev0) (0.11.0)\n", - "Requirement already satisfied: fonttools>=4.22.0 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from matplotlib->eliater==0.0.1.dev0) (4.42.1)\n", - "Requirement already satisfied: kiwisolver>=1.0.1 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from matplotlib->eliater==0.0.1.dev0) (1.4.5)\n", - "Requirement already satisfied: packaging>=20.0 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from matplotlib->eliater==0.0.1.dev0) (23.1)\n", - "Requirement already satisfied: pillow>=6.2.0 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from matplotlib->eliater==0.0.1.dev0) (10.0.1)\n", - "Requirement already satisfied: pyparsing>=2.3.1 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from matplotlib->eliater==0.0.1.dev0) (3.1.1)\n", - "Requirement already satisfied: python-dateutil>=2.7 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from matplotlib->eliater==0.0.1.dev0) (2.8.2)\n", - "Requirement already satisfied: pytz>=2020.1 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from pandas->eliater==0.0.1.dev0) (2023.3.post1)\n", - "Requirement already satisfied: networkx in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from pgmpy>=0.1.24->eliater==0.0.1.dev0) (3.1)\n", - "Requirement already satisfied: scikit-learn in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from pgmpy>=0.1.24->eliater==0.0.1.dev0) (1.3.0)\n", - "Requirement already satisfied: torch in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from pgmpy>=0.1.24->eliater==0.0.1.dev0) (2.0.1)\n", - "Requirement already satisfied: tqdm in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from pgmpy>=0.1.24->eliater==0.0.1.dev0) (4.66.1)\n", - "Requirement already satisfied: joblib in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from pgmpy>=0.1.24->eliater==0.0.1.dev0) (1.3.2)\n", - "Requirement already satisfied: opt-einsum in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from pgmpy>=0.1.24->eliater==0.0.1.dev0) (3.3.0)\n", - "Requirement already satisfied: more-itertools in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from y0>=0.2.5->eliater==0.0.1.dev0) (10.1.0)\n", - "Requirement already satisfied: click in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from y0>=0.2.5->eliater==0.0.1.dev0) (8.1.7)\n", - "Requirement already satisfied: more-click in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from y0>=0.2.5->eliater==0.0.1.dev0) (0.1.2)\n", - "Requirement already satisfied: tabulate in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from y0>=0.2.5->eliater==0.0.1.dev0) (0.9.0)\n", - "Requirement already satisfied: ml-dtypes>=0.2.0 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from jax<0.5.0,>=0.4.8->ananke-causal>=0.5.0->eliater==0.0.1.dev0) (0.3.1)\n", - "Requirement already satisfied: dill>=0.3.7 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from mystic<0.5.0,>=0.4.0->ananke-causal>=0.5.0->eliater==0.0.1.dev0) (0.3.7)\n", - "Requirement already satisfied: klepto>=0.2.4 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from mystic<0.5.0,>=0.4.0->ananke-causal>=0.5.0->eliater==0.0.1.dev0) (0.2.4)\n", - "Requirement already satisfied: sympy>=0.6.7 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from mystic<0.5.0,>=0.4.0->ananke-causal>=0.5.0->eliater==0.0.1.dev0) (1.12)\n", - "Requirement already satisfied: mpmath>=0.19 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from mystic<0.5.0,>=0.4.0->ananke-causal>=0.5.0->eliater==0.0.1.dev0) (1.3.0)\n", - "Requirement already satisfied: six>=1.5 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from python-dateutil>=2.7->matplotlib->eliater==0.0.1.dev0) (1.16.0)\n", - "Requirement already satisfied: patsy>=0.5.2 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from statsmodels<0.14.0,>=0.13.2->ananke-causal>=0.5.0->eliater==0.0.1.dev0) (0.5.3)\n", - "Requirement already satisfied: colorama in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from click->y0>=0.2.5->eliater==0.0.1.dev0) (0.4.6)\n", - "Requirement already satisfied: threadpoolctl>=2.0.0 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from scikit-learn->pgmpy>=0.1.24->eliater==0.0.1.dev0) (3.2.0)\n", - "Requirement already satisfied: filelock in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from torch->pgmpy>=0.1.24->eliater==0.0.1.dev0) (3.12.3)\n", - "Requirement already satisfied: typing-extensions in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from torch->pgmpy>=0.1.24->eliater==0.0.1.dev0) (4.7.1)\n", - "Requirement already satisfied: jinja2 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from torch->pgmpy>=0.1.24->eliater==0.0.1.dev0) (3.1.2)\n", - "Requirement already satisfied: pox>=0.3.3 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from klepto>=0.2.4->mystic<0.5.0,>=0.4.0->ananke-causal>=0.5.0->eliater==0.0.1.dev0) (0.3.3)\n", - "Requirement already satisfied: MarkupSafe>=2.0 in c:\\users\\pnava\\pycharmprojects\\eliater\\venv\\lib\\site-packages (from jinja2->torch->pgmpy>=0.1.24->eliater==0.0.1.dev0) (2.1.3)\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - " Running command git clone --filter=blob:none --quiet https://github.com/y0-causal-inference/eliater.git 'C:\\Users\\pnava\\AppData\\Local\\Temp\\pip-req-build-f3onojbm'\n", - " Running command git checkout -b linear-regression --track origin/linear-regression\n", - " branch 'linear-regression' set up to track 'origin/linear-regression'.\n", - " Switched to a new branch 'linear-regression'\n", - "\n", - "[notice] A new release of pip is available: 23.3.1 -> 23.3.2\n", - "[notice] To update, run: python.exe -m pip install --upgrade pip\n" - ] - } - ], "source": [ - "!pip install git+https://github.com/y0-causal-inference/eliater.git@linear-regression" + "# Motivating example: Figure 4" ] }, { "cell_type": "markdown", - "source": [ - "Figure below is the motivating example in this paper: *Eliater: an open source software for causal query estimation from observational measurements of biomolecular networks*. This graph contains one mediator $M_1$ that connects the exposure $X$ to the outcome $Y$." - ], + "id": "bda3042434f63a26", "metadata": { - "collapsed": false + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } }, - "id": "bda3042434f63a26" + "source": [ + "Figure below is the motivating example in this paper: *Eliater: an open source software for causal query estimation from observational measurements of biomolecular networks*. This graph contains one mediator $M_1$ that connects the exposure $X$ to the outcome $Y$." + ] }, { "cell_type": "code", - "outputs": [], - "source": [ - "import numpy as np\n", - "\n", - "from eliater.examples.frontdoor_backdoor_discrete import (\n", - " single_mediator_with_multiple_confounders_nuisances_discrete_example,\n", - ")" - ], + "execution_count": 1, + "id": "d9f05326d469ff21", "metadata": { - "collapsed": false, "ExecuteTime": { "end_time": "2024-01-17T15:39:13.427336Z", "start_time": "2024-01-17T15:39:13.378250Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false } }, - "id": "d9f05326d469ff21", - "execution_count": 13 + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "from eliater.examples.frontdoor_backdoor_discrete import (\n", + " single_mediator_with_multiple_confounders_nuisances_discrete_example,\n", + ")\n", + "from eliater.network_validation import print_graph_falsifications\n", + "from y0.algorithm.identify import Identification\n", + "from y0.dsl import P, Variable, X, Y\n", + "from y0.algorithm.estimation import estimate_ace\n", + "from y0.algorithm.identify import Identification, identify_outcomes\n", + "from eliater.discover_latent_nodes import find_nuisance_variables, mark_nuisance_variables_as_latent\n", + "from eliater.discover_latent_nodes import remove_nuisance_variables\n", + "import seaborn as sns\n", + "import matplotlib.pyplot as plt" + ] }, { "cell_type": "code", "execution_count": 2, - "outputs": [], - "source": [ - "graph = single_mediator_with_multiple_confounders_nuisances_discrete_example.graph" - ], + "id": "5805a9bed9191161", "metadata": { - "collapsed": false, "ExecuteTime": { "end_time": "2024-01-17T15:25:59.267390Z", "start_time": "2024-01-17T15:25:59.243131Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false } }, - "id": "5805a9bed9191161" + "outputs": [], + "source": [ + "graph = single_mediator_with_multiple_confounders_nuisances_discrete_example.graph" + ] }, { "cell_type": "code", - "execution_count": 21, - "outputs": [], - "source": [ - "data = single_mediator_with_multiple_confounders_nuisances_discrete_example.generate_data(\n", - " num_samples=500, seed=500\n", - ")" - ], + "execution_count": 3, + "id": "24a3d8b93b3804f2", "metadata": { - "collapsed": false, "ExecuteTime": { "end_time": "2024-01-17T15:41:39.934008Z", "start_time": "2024-01-17T15:41:39.904426Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false } }, - "id": "24a3d8b93b3804f2" + "outputs": [], + "source": [ + "data = single_mediator_with_multiple_confounders_nuisances_discrete_example.generate_data(\n", + " num_samples=500, seed=500\n", + ")" + ] }, { "cell_type": "code", + "execution_count": 4, + "id": "cf358541c283c63c", + "metadata": { + "ExecuteTime": { + "end_time": "2024-01-17T15:41:41.133450Z", + "start_time": "2024-01-17T15:41:41.078923Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, "outputs": [ { "data": { - "text/plain": " X M1 Z1 Z2 Z3 R1 R2 R3 Y\n0 1 1 1 1 1 1 1 1 1\n1 1 1 1 0 0 1 1 1 1\n2 1 0 1 0 1 1 1 0 1\n3 1 1 1 1 1 1 1 1 0\n4 1 1 1 1 0 1 1 1 1", - "text/html": "
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n
XM1Z1Z2Z3R1R2R3Y
0111111111
1111001111
2101011101
3111111110
4111101111
\n
" + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
XM1Z1Z2Z3R1R2R3Y
0111111111
1111001111
2101011101
3111111110
4111101111
\n", + "
" + ], + "text/plain": [ + " X M1 Z1 Z2 Z3 R1 R2 R3 Y\n", + "0 1 1 1 1 1 1 1 1 1\n", + "1 1 1 1 0 0 1 1 1 1\n", + "2 1 0 1 0 1 1 1 0 1\n", + "3 1 1 1 1 1 1 1 1 0\n", + "4 1 1 1 1 0 1 1 1 1" + ] }, - "execution_count": 22, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.head()" - ], + ] + }, + { + "cell_type": "markdown", + "id": "7f8b8c77d0b42c3e", "metadata": { "collapsed": false, - "ExecuteTime": { - "end_time": "2024-01-17T15:41:41.133450Z", - "start_time": "2024-01-17T15:41:41.078923Z" + "jupyter": { + "outputs_hidden": false } }, - "id": "cf358541c283c63c", - "execution_count": 22 - }, - { - "cell_type": "markdown", "source": [ "## Step 1: Verify correctness of the network structure" - ], - "metadata": { - "collapsed": false - }, - "id": "7f8b8c77d0b42c3e" + ] }, { "cell_type": "code", - "outputs": [], - "source": [ - "from eliater.network_validation import print_graph_falsifications" - ], + "execution_count": 5, + "id": "4e3220a2c48cb9b4", "metadata": { - "collapsed": false, "ExecuteTime": { - "end_time": "2024-01-17T15:41:42.833594Z", - "start_time": "2024-01-17T15:41:42.813867Z" + "end_time": "2024-01-17T15:41:45.144807Z", + "start_time": "2024-01-17T15:41:43.311539Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false } }, - "id": "61f4ae2447e8e353", - "execution_count": 23 - }, - { - "cell_type": "code", - "execution_count": 24, "outputs": [ { "name": "stdout", @@ -230,156 +265,161 @@ "====== ======= ======= ========= ======== ===== ======= ===================\n", "left right given stats p dof p_adj p_adj_significant\n", "====== ======= ======= ========= ======== ===== ======= ===================\n", - "X Z2 Z1 0.34346 0.842206 2 1 False\n", - "X Z3 Z1 1.12103 0.570914 2 1 False\n", - "X Y M1|Z1 0.0346805 0.999851 4 1 False\n", - "M1 Z3 Z1 0.630044 0.729773 2 1 False\n", - "M1 Z2 Z1 1.51285 0.469341 2 1 False\n", + "R3 X M1|Z2 3.72468 0.444547 4 1 False\n", "Y Z1 X|Z2 3.842 0.427811 4 1 False\n", - "R3 Z3 R1|Y 1.20928 0.876569 4 1 False\n", - "R1 Z2 Z1 0.988663 0.609979 2 1 False\n", - "M1 Z1 X 1.44905 0.484554 2 1 False\n", - "R2 Z1 R1 0.151588 0.927007 2 1 False\n", - "R2 X R1 0.0114721 0.99428 2 1 False\n", - "R1 Z3 Z1 1.39424 0.498017 2 1 False\n", - "Y Z2 Z1|Z3 2.65857 0.616484 4 1 False\n", "M1 R2 R1 0.148438 0.928469 2 1 False\n", - "R1 R3 M1|R2 0.527356 0.970784 4 1 False\n", "R1 Z1 X 0.0530231 0.973837 2 1 False\n", - "R2 Z2 Z1 0.221862 0.895 2 1 False\n", + "R3 Z2 X|Z3 5.99107 0.112045 3 1 False\n", + "R3 Z3 R2|Y 0.149065 0.997357 4 1 False\n", "R2 Y R1 0.543749 0.76195 2 1 False\n", + "R1 Z2 X 3.57855 0.167081 2 1 False\n", + "M1 R3 R2|Y 0.901882 0.924291 4 1 False\n", "R1 X M1 0.227982 0.892266 2 1 False\n", - "R2 Z3 Z1 1.66956 0.43397 2 1 False\n", - "Z1 Z3 Z2 0.698232 0.705311 2 1 False\n", - "M1 R3 R1|Y 0.307719 0.98931 4 1 False\n", - "R3 Z1 R1|Y 1.29458 0.862294 4 1 False\n", - "R3 X M1|Z1 2.18017 0.702662 4 1 False\n", + "R1 Z3 X 1.41849 0.492015 2 1 False\n", + "X Z2 Z1 0.34346 0.842206 2 1 False\n", + "X Z3 Z2 1.73587 0.419817 2 1 False\n", + "R1 R3 R2|Y 0.818848 0.935903 4 1 False\n", "R1 Y M1 4.41965 0.10972 2 1 False\n", - "R3 Z2 Z1|Z3 2.85924 0.58165 4 1 False\n", + "R3 Z1 X|Z2 1.03831 0.791983 3 1 False\n", + "R2 Z2 X 0.529657 0.767338 2 1 False\n", + "R2 Z1 X 0.186598 0.910921 2 1 False\n", + "Z1 Z3 Z2 0.698232 0.705311 2 1 False\n", + "R2 Z3 X 4.48986 0.105935 2 1 False\n", + "M1 Z1 X 1.44905 0.484554 2 1 False\n", + "X Y M1|Z2 2.04254 0.727934 4 1 False\n", + "R2 X R1 0.0114721 0.99428 2 1 False\n", + "M1 Z3 X 0.617323 0.734429 2 1 False\n", + "Y Z2 X|Z3 1.1174 0.891501 4 1 False\n", + "M1 Z2 X 0 1 2 1 False\n", "====== ======= ======= ========= ======== ===== ======= ===================\n" ] } ], "source": [ "print_graph_falsifications(graph, data, method=\"chi-square\", verbose=True, significance_level=0.01)" - ], + ] + }, + { + "cell_type": "markdown", + "id": "8b40df88c10664e3", "metadata": { "collapsed": false, - "ExecuteTime": { - "end_time": "2024-01-17T15:41:45.144807Z", - "start_time": "2024-01-17T15:41:43.311539Z" + "jupyter": { + "outputs_hidden": false } }, - "id": "4e3220a2c48cb9b4" - }, - { - "cell_type": "markdown", "source": [ "All the d-separations implied by the network are validated by the data. No test failed. Hence, we can proceed to step 2." - ], - "metadata": { - "collapsed": false - }, - "id": "8b40df88c10664e3" + ] }, { "cell_type": "markdown", + "id": "4435aa5f3a8f55b9", + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, "source": [ "## Step 2: Check query identifiability\n", "\n", "The causal query of interest is the average treatment effect of $X$ on $Y$, defined as:\n", "$E[Y|do(X=1)] - E[Y|do(X=0)]$." - ], - "metadata": { - "collapsed": false - }, - "id": "4435aa5f3a8f55b9" + ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, + "id": "6a1d0da707ca1d2f", + "metadata": { + "ExecuteTime": { + "end_time": "2024-01-17T15:38:31.023301Z", + "start_time": "2024-01-17T15:38:30.975649Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, "outputs": [ { "data": { - "text/plain": "Identification(outcomes=\"{Y}, treatments=\"{X}\",conditions=\"set()\", graph=\"NxMixedGraph(directed=, undirected=)\", estimand=\"P(M1, R1, R2, R3, X, Y, Z1, Z2, Z3)\")" + "text/latex": [ + "$\\sum\\limits_{M_1, Z_1, Z_2, Z_3} P(M_1 | X, Z_1) P(Y | M_1, X, Z_1, Z_2, Z_3) P(Z_2 | Z_1) P(Z_3 | Z_1, Z_2) \\sum\\limits_{M_1, X, Y, Z_2, Z_3} \\sum\\limits_{R_1, R_2, R_3} P(M_1, R_1, R_2, R_3, X, Y, Z_1, Z_2, Z_3)$" + ], + "text/plain": [ + "Sum[M1, Z1, Z2, Z3](P(M1 | X, Z1) * P(Y | M1, X, Z1, Z2, Z3) * P(Z2 | Z1) * P(Z3 | Z1, Z2) * Sum[M1, X, Y, Z2, Z3](Sum[R1, R2, R3](P(M1, R1, R2, R3, X, Y, Z1, Z2, Z3))))" + ] }, - "execution_count": 5, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "from y0.algorithm.identify import Identification\n", - "from y0.dsl import P, Variable\n", - "\n", - "id_in = Identification.from_expression(\n", - " query=P(Variable(\"Y\") @ Variable(\"X\")),\n", - " graph=graph,\n", - ")\n", - "id_in" - ], + "identify_outcomes(graph=graph, treatments=X, outcomes=Y)" + ] + }, + { + "cell_type": "markdown", + "id": "5e27342ad3164b42", "metadata": { "collapsed": false, - "ExecuteTime": { - "end_time": "2024-01-17T15:38:31.023301Z", - "start_time": "2024-01-17T15:38:30.975649Z" + "jupyter": { + "outputs_hidden": false } }, - "id": "6a1d0da707ca1d2f" - }, - { - "cell_type": "markdown", "source": [ "The query is identifiable. Hence we can proceed to step 3." - ], - "metadata": { - "collapsed": false - }, - "id": "5e27342ad3164b42" + ] }, { "cell_type": "markdown", - "source": [ - "## Step 3: Find nuisance variables and mark them as latent" - ], + "id": "2ed3073afbc7030a", "metadata": { - "collapsed": false + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } }, - "id": "2ed3073afbc7030a" + "source": [ + "## Step 3: Find nuisance variables and mark them as latent" + ] }, { - "cell_type": "code", - "outputs": [], - "source": [ - "from eliater.discover_latent_nodes import find_nuisance_variables, mark_nuisance_variables_as_latent" - ], + "cell_type": "markdown", + "id": "b7661cedf357ad01", "metadata": { "collapsed": false, - "ExecuteTime": { - "end_time": "2024-01-17T15:38:39.681358Z", - "start_time": "2024-01-17T15:38:39.658772Z" + "jupyter": { + "outputs_hidden": false } }, - "id": "2f30339a39dca602", - "execution_count": 6 - }, - { - "cell_type": "markdown", "source": [ "This function finds the nuisance variables for the input graph." - ], - "metadata": { - "collapsed": false - }, - "id": "b7661cedf357ad01" + ] }, { "cell_type": "code", "execution_count": 7, + "id": "c094ba6186dfecf6", + "metadata": { + "ExecuteTime": { + "end_time": "2024-01-17T15:38:41.126011Z", + "start_time": "2024-01-17T15:38:41.054670Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, "outputs": [ { "data": { - "text/plain": "{R1, R2, R3}" + "text/plain": [ + "{R1, R2, R3}" + ] }, "execution_count": 7, "metadata": {}, @@ -387,222 +427,246 @@ } ], "source": [ - "nuisance_variables = find_nuisance_variables(\n", - " graph, treatments=Variable(\"X\"), outcomes=Variable(\"Y\")\n", - ")\n", + "nuisance_variables = find_nuisance_variables(graph, treatments=X, outcomes=Y)\n", "nuisance_variables" - ], + ] + }, + { + "cell_type": "markdown", + "id": "79a36765bb1640f6", "metadata": { "collapsed": false, - "ExecuteTime": { - "end_time": "2024-01-17T15:38:41.126011Z", - "start_time": "2024-01-17T15:38:41.054670Z" + "jupyter": { + "outputs_hidden": false } }, - "id": "c094ba6186dfecf6" - }, - { - "cell_type": "markdown", "source": [ "The nuisance variables are $R_1$, $R_2$, and $R_3$." - ], - "metadata": { - "collapsed": false - }, - "id": "79a36765bb1640f6" + ] }, { "cell_type": "markdown", - "source": [ - "## Step 4: Simplify the network" - ], + "id": "b86a9300fa3d3881", "metadata": { - "collapsed": false + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } }, - "id": "b86a9300fa3d3881" - }, - { - "cell_type": "markdown", "source": [ - "The following function finds the nuisance variable (step 3), marks them as latent and then applies Evan's simplification rules to remove the nuisance variables. As a result, running the 'find_nuisance_variables' and 'mark_nuisance_variables_as_latent' functions in step 3 is not necessary to get the value of step 4. However, we called them to illustrate the results. The new graph obtained in step 4 does not contain the nuisance variables. " - ], - "metadata": { - "collapsed": false - }, - "id": "92e386966d986cec" + "## Step 4: Simplify the network" + ] }, { - "cell_type": "code", - "outputs": [], - "source": [ - "from eliater.discover_latent_nodes import remove_nuisance_variables" - ], + "cell_type": "markdown", + "id": "92e386966d986cec", "metadata": { "collapsed": false, - "ExecuteTime": { - "end_time": "2024-01-17T15:38:43.659886Z", - "start_time": "2024-01-17T15:38:43.614147Z" + "jupyter": { + "outputs_hidden": false } }, - "id": "18dce64dfc819025", - "execution_count": 8 + "source": [ + "The following function finds the nuisance variable (step 3), marks them as latent and then applies Evan's simplification rules to remove the nuisance variables. As a result, running the 'find_nuisance_variables' and 'mark_nuisance_variables_as_latent' functions in step 3 is not necessary to get the value of step 4. However, we called them to illustrate the results. The new graph obtained in step 4 does not contain the nuisance variables. " + ] }, { "cell_type": "code", - "execution_count": 9, - "outputs": [], - "source": [ - "new_graph = remove_nuisance_variables(graph, treatments=Variable(\"X\"), outcomes=Variable(\"Y\"))" - ], + "execution_count": 8, + "id": "e1d0f93c0d993112", "metadata": { - "collapsed": false, "ExecuteTime": { "end_time": "2024-01-17T15:38:44.092780Z", "start_time": "2024-01-17T15:38:44.067590Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false } }, - "id": "e1d0f93c0d993112" + "outputs": [], + "source": [ + "new_graph = remove_nuisance_variables(graph, treatments=X, outcomes=Y)" + ] }, { "cell_type": "markdown", - "source": [ - "## Step 5: Estimate the query" - ], + "id": "f1b8adaf922f3096", "metadata": { - "collapsed": false + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } }, - "id": "f1b8adaf922f3096" + "source": [ + "## Step 5: Estimate the query" + ] }, { "cell_type": "code", - "execution_count": 10, - "outputs": [], - "source": [ - "from y0.algorithm.estimation import estimate_ace" - ], + "execution_count": 9, + "id": "55a147fe3a3ee98d", "metadata": { - "collapsed": false, "ExecuteTime": { - "end_time": "2024-01-17T15:38:45.639668Z", - "start_time": "2024-01-17T15:38:45.605103Z" + "end_time": "2024-01-17T15:41:01.399162Z", + "start_time": "2024-01-17T15:41:01.075146Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false } }, - "id": "9544bf2cd9459cc" - }, - { - "cell_type": "code", - "execution_count": 19, "outputs": [ { "data": { - "text/plain": "0.20915697893053198" + "text/plain": [ + "0.20915697893041166" + ] }, - "execution_count": 19, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ATE_value = estimate_ace(\n", - " graph=new_graph, treatments=Variable(\"X\"), outcomes=Variable(\"Y\"), data=data\n", + " graph=new_graph, treatments=X, outcomes=Y, data=data\n", ")\n", "ATE_value" - ], + ] + }, + { + "cell_type": "markdown", + "id": "8b1fbdecffb16328", "metadata": { "collapsed": false, - "ExecuteTime": { - "end_time": "2024-01-17T15:41:01.399162Z", - "start_time": "2024-01-17T15:41:01.075146Z" + "jupyter": { + "outputs_hidden": false } }, - "id": "55a147fe3a3ee98d" - }, - { - "cell_type": "markdown", "source": [ "The ATE amounts to 0.21 meaning that the average effect that $X$ has on $Y$ is negative." - ], - "metadata": { - "collapsed": false - }, - "id": "8b1fbdecffb16328" + ] }, { "cell_type": "markdown", + "id": "e362c47381bc281f", + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, "source": [ - "#Evaluation Criterion\n", + "# Evaluation Criterion\n", + "\n", "As we used synthetic data set, we were able to generate two interventional data sets where in\n", "one X was set to 1, and the other one X is set to 0. The ATE was calculated by subtracting the average value of Y obtained from each interventional data,\n", "resulting in the ground truth ATE=0.01. The ATE indicates that increase in X can increase Y levels." - ], - "metadata": { - "collapsed": false - }, - "id": "e362c47381bc281f" + ] }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 23, + "id": "7aa1b23565adae07", + "metadata": { + "ExecuteTime": { + "end_time": "2024-01-17T15:41:11.036266Z", + "start_time": "2024-01-17T15:41:10.992234Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "0.010000000000000009\n" + "The true ATE is 0.01\n" ] } ], "source": [ "# get interventional data where EGFR is set to 1\n", "intv_data_X_1 = single_mediator_with_multiple_confounders_nuisances_discrete_example.generate_data(\n", - " num_samples=500, seed=500, treatments={Variable(\"X\"): 1}\n", + " num_samples=500, seed=500, treatments={X: 1}\n", ")\n", "\n", "# get interventional data where EGFR is set to 0\n", "intv_data_X_0 = single_mediator_with_multiple_confounders_nuisances_discrete_example.generate_data(\n", - " num_samples=500, seed=500, treatments={Variable(\"X\"): 0}\n", + " num_samples=500, seed=500, treatments={X: 0}\n", ")\n", "\n", - "# get the true value of ATE\n", - "print(np.mean(intv_data_X_1[\"Y\"]) - np.mean(intv_data_X_0[\"Y\"]))" - ], - "metadata": { - "collapsed": false, - "ExecuteTime": { - "end_time": "2024-01-17T15:41:11.036266Z", - "start_time": "2024-01-17T15:41:10.992234Z" - } - }, - "id": "7aa1b23565adae07" + "true_ate = intv_data_X_1.mean()[\"Y\"] - intv_data_X_0.mean()[\"Y\"] \n", + "print(f\"The true ATE is {true_ate:.04}\")" + ] + }, + { + "cell_type": "markdown", + "id": "7ae64b2d-be03-4053-a43b-40ec8ebc39fe", + "metadata": {}, + "source": [ + "If you don't set the random seed and take a sampling over many generated datasets, you can plot the ATE. TODO: how to determine the convidence in the ATE?" + ] }, { "cell_type": "code", - "execution_count": null, - "outputs": [], - "source": [], - "metadata": { - "collapsed": false - }, - "id": "b53fbb13d2711b89" + "execution_count": 39, + "id": "e73d76fc-a451-47f5-a92a-005a3ba7d371", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAEECAYAAAAoOz76AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIl0lEQVR4nO3dd1gUV/s38O+C9K5UC0VQwd4RYxfFxF5iiwbsvSbRmDwGNSZ2JVGjMU9EY+MXjRqT2BDFWLD33kBsqKAUBWl73j98mYdxF1gQWGC/n+vaS/fMmTP3nJ2ZvZmdM6MQQggQERER6QA9bQdAREREVFyY+BAREZHOYOJDREREOoOJDxEREekMJj5ERESkM5j4EBERkc5g4kNEREQ6g4kPERER6QwmPkRERKQzmPgQUZFTKBSYNWtWobUXHh4OhUKB8PDwQmuTiHQDEx+iUuLy5cvo06cPXFxcYGxsjEqVKqFDhw5Yvny5tkMrFX766ScoFAp4e3vLyl1dXaFQKPJ8rVu3DgByrTN69GgtrBkR5Uc5bQdARHk7fvw42rZtC2dnZ4wYMQKOjo548OABTpw4gR9++AETJkzQdogl3qZNm+Dq6opTp07hzp078PDwAAAEBQXh1atXUr3du3djy5YtWLZsGWxtbaXy5s2bS//v0KEDPv30U5VlVK9evQjXgIgKAxMfolLgu+++g5WVFU6fPg1ra2vZtGfPnmknqFIkMjISx48fx/bt2zFq1Chs2rQJgYGBAIAePXrI6sbExGDLli3o0aMHXF1d1bZXvXp1DBo0qIijJqKiwJ+6iEqBu3fvolatWipJDwDY29vL3gcHB6Ndu3awt7eHkZERatasiVWrVqnM5+rqii5duiA8PByNGzeGiYkJ6tSpI103s337dtSpUwfGxsZo1KgRzp8/L5s/ICAA5ubmuHfvHvz8/GBmZoaKFStizpw5EELkuU6PHj3C0KFD4eDgACMjI9SqVQtr165Vqffw4UP06NEDZmZmsLe3x5QpU5Camppn+9lt2rQJNjY26Ny5M/r06YNNmzbla34iKjt4xoeoFHBxcUFERASuXLmC2rVr51p31apVqFWrFrp164Zy5crhr7/+wtixY6FUKjFu3DhZ3Tt37mDgwIEYNWoUBg0ahMWLF6Nr165YvXo1vvrqK4wdOxYAMG/ePPTt2xc3b96Ent7//l7KzMxEp06d0KxZMyxcuBB79+5FYGAgMjIyMGfOnBxjfPr0KZo1awaFQoHx48fDzs4Oe/bswbBhw5CYmIjJkycDAFJSUtC+fXtER0dj4sSJqFixIjZs2ICDBw/mq/82bdqEXr16wdDQEAMGDMCqVatw+vRpNGnSJF/tZHnz5g1iY2NVyi0tLWFoaFigNomomAgiKvH2798v9PX1hb6+vvDx8RHTpk0T+/btE2lpaSp1k5OTVcr8/PxE1apVZWUuLi4CgDh+/LhUtm/fPgFAmJiYiPv370vlP//8swAgDh06JJX5+/sLAGLChAlSmVKpFJ07dxaGhobi+fPnUjkAERgYKL0fNmyYcHJyErGxsbKY+vfvL6ysrKR1CAoKEgDE77//LtV5/fq18PDwUIknJ2fOnBEARGhoqBRj5cqVxaRJk9TWX7RokQAgIiMj1U4HkONry5YtecZDRNrFn7qISoEOHTogIiIC3bp1w8WLF7Fw4UL4+fmhUqVK2LVrl6yuiYmJ9P+EhATExsaidevWuHfvHhISEmR1a9asCR8fH+l91oindu3awdnZWaX83r17KrGNHz9e+n/WGZy0tDQcOHBA7boIIfDHH3+ga9euEEIgNjZWevn5+SEhIQHnzp0D8PZCYycnJ/Tp00ea39TUFCNHjsy9w7LZtGkTHBwc0LZtWynGfv36ISQkBJmZmRq3k1337t0RGhqq8spaBhGVXPypi6iUaNKkCbZv3460tDRcvHgRO3bswLJly9CnTx9cuHABNWvWBAAcO3YMgYGBiIiIQHJysqyNhIQEWFlZSe+zJzcApGlVqlRRW/7y5UtZuZ6eHqpWrSoryxrZFBUVpXY9nj9/jvj4eKxZswZr1qxRWyfrgu379+/Dw8MDCoVCNr1GjRpq53tXZmYmQkJC0LZtW0RGRkrl3t7eWLJkCcLCwtCxY0eN2squcuXK8PX1zfd8RKR9THyIShlDQ0M0adIETZo0QfXq1TFkyBBs3boVgYGBuHv3Ltq3bw9PT08sXboUVapUgaGhIXbv3o1ly5ZBqVTK2tLX11e7jJzKhQYXLeclK4ZBgwbB399fbZ26deu+93IA4ODBg3jy5AlCQkIQEhKiMn3Tpk0FSnyIqPRi4kNUijVu3BgA8OTJEwDAX3/9hdTUVOzatUt2NufQoUNFsnylUol79+7J7l9z69YtAMhxKLidnR0sLCyQmZmZ51kTFxcXXLlyBUII2VmfmzdvahTfpk2bYG9vj5UrV6pM2759O3bs2IHVq1fLfh4korKN1/gQlQKHDh1Se7Zl9+7dAP7300/WmZrsdRMSEhAcHFxksa1YsUL6vxACK1asgIGBAdq3b6+2vr6+Pnr37o0//vgDV65cUZn+/Plz6f8fffQRHj9+jG3btkllycnJOf5Ell1KSgq2b9+OLl26oE+fPiqv8ePHIykpSeUaKSIq23jGh6gUmDBhApKTk9GzZ094enoiLS0Nx48fx//93//B1dUVQ4YMAQB07NgRhoaG6Nq1K0aNGoVXr17hl19+gb29vXRWqDAZGxtj79698Pf3h7e3N/bs2YN//vkHX331Fezs7HKcb/78+Th06BC8vb0xYsQI1KxZEy9evMC5c+dw4MABvHjxAgAwYsQIrFixAp9++inOnj0LJycnbNiwAaampnnGtmvXLiQlJaFbt25qpzdr1gx2dnbYtGkT+vXrl6/1vnXrFjZu3KhS7uDggA4dOuSrLSIqZlocUUZEGtqzZ48YOnSo8PT0FObm5sLQ0FB4eHiICRMmiKdPn8rq7tq1S9StW1cYGxsLV1dXsWDBArF27VqVIdouLi6ic+fOKssCIMaNGycri4yMFADEokWLpDJ/f39hZmYm7t69Kzp27ChMTU2Fg4ODCAwMFJmZmSptZh/OLoQQT58+FePGjRNVqlQRBgYGwtHRUbRv316sWbNGVu/+/fuiW7duwtTUVNja2opJkyaJvXv35jmcvWvXrsLY2Fi8fv06xzoBAQHCwMBANqz+fYazt27dOsdlEVHJoBCiEK5WJCKdExAQgG3btsmec0VEVNLxGh8iIiLSGUx8iIiISGcw8SEiIiKdwWt8iIiISGfwjA8RERHpDCY+REREpDPylfjMmjULCoUCsbGxedZ1dXVFQEBAQeMqdAEBATA3Ny/25bq6uqJLly7FvtySQqFQYNasWdoOQyYqKgoKhQLr1q0r1HYDAgJUHtNQEte/qGUdJ4pKeHg4FAoFwsPDC7VdXfysSps2bdqgTZs20vui2pdLg6LaD7Irin3i3c9QG3jGh6gM2r17N7/E1WC/lDxZyUtOr/nz5xfJco8fP45Zs2YhPj6+SNrXZdeuXcOsWbMQFRWl7VDU4iMriIpQSkoKypUr/t1s9+7dWLlyJb/k35Fbv2jrs6K3BgwYgI8++kilvEGDBjnO4+LigpSUFBgYGOR7ecePH8fs2bMREBAAa2vrfM9PObt27Rpmz56NNm3aqJwF379/v3aCyoZ7eRnw+vVrmJmZaTuMIieEwJs3b0rVk7SNjY3zrKMrn19Jp8lnVVokJydr9DyzkqRhw4YYNGhQvuZRKBQl7nMrjX1fnAwNDbUdQsF+6oqNjUXfvn1haWmJChUqYNKkSXjz5k2u8+T0u/+6deugUChUTont2bMHLVu2hJmZGSwsLNC5c2dcvXpVVic9PR03btzI18MX7927Bz8/P5iZmaFixYqYM2eOylOvFy9ejObNm6NChQowMTFBo0aNZE+Hzm7jxo1o2rQpTE1NYWNjg1atWuWZ0a5fvx7lypXDF198IZXFxcVh8ODBsLS0hLW1Nfz9/XHx4kWV36+zrlW6e/cuPvroI1hYWOCTTz4B8PYL9LPPPkOVKlVgZGSEGjVqYPHixbL1y+038Xd/z836zO7cuSP9VWRlZYUhQ4YgOTlZNm9qaiqmTJkCOzs7WFhYoFu3bnj48GGu/ZCXrOuj9u3bh8aNG8PExAQ///wzACA+Ph6TJ0+W1tXDwwMLFiyAUqmUtREfH4+AgABYWVlJ/aru1HZ+t6WdO3eidu3aMDY2Ru3atbFjxw619XLq02vXrmHgwIGwsbFBixYtpOkbN25Eo0aNYGJigvLly6N///548OCBSrsnT57ERx99BBsbG5iZmaFu3br44YcfALzdRlauXCktP+uVlz179qB169awsLCApaUlmjRpgs2bN0vTjxw5go8//hjOzs4wMjJClSpVMGXKFKSkpGjUZ3ntKzldT6DJ9YKaxJZXv6hb/vnz5/Hhhx/C0tIS5ubmaN++PU6cOCGrk3UMO3bsGKZOnQo7OzuYmZmhZ8+esifN59f9+/fRrVs3mJmZwd7eHlOmTMG+fftUruto06YNateujbNnz6JVq1YwNTXFV199BeDtfhkYGAgPDw+pX6ZNm4bU1FSV5Wmy7WUt69q1a2jbti1MTU1RqVIlLFy4UKW96Oho3Lhxo8Drrwl1x7NLly4hICAAVatWhbGxMRwdHTF06FDExcVJdWbNmiUdf93c3KRtIfv3UH76Q13f5+TGjRvo27cv7OzsYGJigho1auDrr7+Wpt+/fx9jx45FjRo1YGJiggoVKuDjjz/W+Gej3I4NWTGru8ZG3TWK79IktnXr1uHjjz8GALRt21bq26xtVt3ynz17hmHDhsHBwQHGxsaoV68e1q9fL6uT9VkvXrwYa9asgbu7O4yMjNCkSROcPn1ao77JUqAzPn379oWrqyvmzZuHEydO4Mcff8TLly/x22+/FaQ5FRs2bIC/vz/8/PywYMECJCcnY9WqVWjRogXOnz8vfTiPHj2Cl5cX/P39Nbq4LTMzE506dUKzZs2wcOFC7N27F4GBgcjIyMCcOXOkej/88AO6deuGTz75BGlpaQgJCcHHH3+Mv//+G507d5bqzZ49G7NmzULz5s0xZ84cGBoa4uTJkzh48CA6duyoNoY1a9Zg9OjR+OqrrzB37lwAgFKpRNeuXXHq1CmMGTMGnp6e+PPPP+Hv76+2jYyMDPj5+aFFixZYvHgxTE1NIYRAt27dcOjQIQwbNgz169fHvn378MUXX+DRo0dYtmyZhr2vqm/fvnBzc8O8efNw7tw5/Pe//4W9vT0WLFgg1Rk+fDg2btyIgQMHonnz5jh48KCsrwrq5s2bGDBgAEaNGoURI0agRo0aSE5ORuvWrfHo0SOMGjUKzs7OOH78OGbMmIEnT54gKCgIwNszRN27d8fRo0cxevRoeHl5YceOHWr7NT/b0v79+9G7d2/UrFkT8+bNQ1xcHIYMGYLKlStrvF4ff/wxqlWrhu+//15KTL/77jvMnDkTffv2xfDhw/H8+XMsX74crVq1wvnz56XT8aGhoejSpQucnJwwadIkODo64vr16/j7778xadIkjBo1Co8fP0ZoaCg2bNigUTzr1q3D0KFDUatWLcyYMQPW1tY4f/489u7di4EDBwIAtm7diuTkZIwZMwYVKlTAqVOnsHz5cjx8+BBbt27Ntf2C7Cv5oUls+e2Xq1evomXLlrC0tMS0adNgYGCAn3/+GW3atMHhw4fh7e0tqz9hwgTY2NggMDAQUVFRCAoKwvjx4/F///d/+V6f169fo127dnjy5In0GW/evBmHDh1SWz8uLg4ffvgh+vfvj0GDBsHBwQFKpRLdunXD0aNHMXLkSHh5eeHy5ctYtmwZbt26hZ07d0rza7rtAcDLly/RqVMn9OrVC3379sW2bdswffp01KlTBx9++KFU79NPP8Xhw4dV/rDMSXJystqBM9bW1vn6CTI0NBT37t3DkCFD4OjoiKtXr2LNmjW4evUqTpw4AYVCgV69euHWrVvYsmULli1bBltbWwCAnZ1dvvtDXd/n5NKlS2jZsiUMDAwwcuRIuLq64u7du/jrr7/w3XffAQBOnz6N48ePo3///qhcuTKioqKwatUqtGnTBteuXcv1bFJex4b3pUlsrVq1wsSJE/Hjjz/iq6++gpeXFwBI/74rJSUFbdq0wZ07dzB+/Hi4ublh69atCAgIQHx8vErcmzdvRlJSEkaNGgWFQoGFCxeiV69euHfvnuY/eebniaaBgYECgOjWrZusfOzYsQKAuHjxolTm4uIi/P39VeZ9V3BwsOxJyElJScLa2lqMGDFCVi8mJkZYWVnJyrOeGJ19OTnx9/cXAMSECROkMqVSKTp37iwMDQ3F8+fPpfLk5GTZvGlpaaJ27dqiXbt2Utnt27eFnp6e6Nmzp8qTqJVKpawfsp6A/cMPPwiFQiG+/fZbWf0//vhDABBBQUFSWWZmpmjXrp0AIIKDg1XW48svv5S1sXPnTgFAzJ07V1bep08foVAoxJ07d4QQ/+uz7G1mwTtP0M76zIYOHSqr17NnT1GhQgXp/YULFwQAMXbsWFm9gQMHqn0qt6ZcXFwEALF3715Z+bfffivMzMzErVu3ZOVffvml0NfXF9HR0UKI//XJwoULpToZGRmiZcuWKn2Qn22pfv36wsnJScTHx0tl+/fvFwCEi4uLrG5OfTpgwABZvaioKKGvry++++47Wfnly5dFuXLlpPKMjAzh5uYmXFxcxMuXL2V1s29348aNU7u/qRMfHy8sLCyEt7e3SElJybHNd/cLIYSYN2+eUCgU4v79+yrrmEXTfSWnbeXdY8mhQ4dUnsyuaWy59cu7y+/Ro4cwNDQUd+/elcoeP34sLCwsRKtWraSyrGOYr6+vbH2mTJki9PX1ZduJppYsWSIAiJ07d0plKSkpwtPTU2XdW7duLQCI1atXy9rYsGGD0NPTE0eOHJGVr169WgAQx44dE0Jovu1lX9Zvv/0mlaWmpgpHR0fRu3dv2fxZdfOSte/l9IqIiJC12bp1a5V5s+/L6raFLVu2CADi33//lcoWLVok++7JUpD+eLfvc9KqVSthYWEh2yaFyHs/i4iIUOn3d/cDTY8N7/ZhFn9//zyPX5rGtnXrVpXtNKflBwUFCQBi48aNUllaWprw8fER5ubmIjExUQjxv8+6QoUK4sWLF1LdP//8UwAQf/31l8qyclKgn7rGjRsnez9hwgQAby8cfF+hoaGIj4/HgAEDEBsbK7309fXh7e0t+4vH1dUVQoh8DWUcP3689H+FQoHx48cjLS0NBw4ckMqzX0Py8uVLJCQkoGXLljh37pxUvnPnTiiVSnzzzTfQ05N3o7qfFRYuXIhJkyZhwYIF+M9//iObtnfvXhgYGGDEiBFSmZ6enko/ZzdmzBjZ+927d0NfXx8TJ06UlX/22WcQQmDPnj05tpWX0aNHy963bNkScXFxSExMlJYNQGXZkydPLvAys7i5ucHPz09WtnXrVrRs2RI2NjaybcTX1xeZmZn4999/pbjKlSsn6yt9fX1pe81O023pyZMnuHDhAvz9/WFlZSWVd+jQATVr1tR4vd7t0+3bt0OpVKJv376ydXJ0dES1atWk7f78+fOIjIzE5MmTVS7ILOgQ8tDQUCQlJeHLL79UuV4ie5vZ94vXr18jNjYWzZs3hxAC58+fz7H9/O4rBVHQ2HKSmZmJ/fv3o0ePHqhatapU7uTkhIEDB+Lo0aPS9p9l5MiRsvVp2bIlMjMzcf/+/Xwvf+/evahUqRK6desmlRkbG8uOEdkZGRlhyJAhsrKtW7fCy8sLnp6esm2qXbt2ACBtU5pue1nMzc1l1+IYGhqiadOmuHfvnqxeeHi4xmd7gLf9FxoaqvLKz34FyLeFN2/eIDY2Fs2aNQMA2TE8J/ntD3V9r87z58/x77//YujQoXB2dpZNy2k/S09PR1xcHDw8PGBtbZ1r/EVxbHhXQWPLze7du+Ho6IgBAwZIZQYGBpg4cSJevXqFw4cPy+r369cPNjY20vuWLVsCgMr2l5sC/dRVrVo12Xt3d3fo6ekVytC127dvA4C0c77L0tKywG3r6enJDmIAUL16dQCQxf73339j7ty5uHDhguy38Owbz927d6Gnp6fRTnn48GH8888/mD59uuy6niz379+Hk5OTyilMDw8Pte2VK1dO5WeV+/fvo2LFirCwsJCVZ51eLMjBN8u7O2nWRvfy5UtYWlri/v370NPTg7u7u6xejRo1CrzMLG5ubiplt2/fxqVLl6TT0u969uwZgP/167v3b3qfuLL68d19IKtdTXf+d9fr9u3bEEKobReAdAr37t27AIDatWtrHHNeNG0zOjoa33zzDXbt2oWXL1/KpiUkJOTavqb7SkEVNLacPH/+HMnJyWq3FS8vLyiVSjx48AC1atWSynPbT/Lr/v37cHd3V/nCyumYUKlSJZWLRm/fvo3r16/nuZ9ouu1lqVy5skpcNjY2uHTpUs4rpIFq1arB19f3vdoAgBcvXmD27NkICQmR1jGLJttCfvtDXd+rk/XFnNd+lpKSgnnz5iE4OBiPHj2SJY957WeatP8+Chpbbu7fv49q1aqp/FGU03dXYexnhTKqS5NsMqc6mZmZsvdZF6du2LABjo6OKvWLerjpkSNH0K1bN7Rq1Qo//fQTnJycYGBggODgYNmFnvlRq1YtxMfHY8OGDRg1apTaL/P8MDIyUtlINKXp55Cdvr6+2vL8/DVXUOpGcCmVSnTo0AHTpk1TO09WMluSvbteSqUSCoUCe/bsUdvf2rj5ZnaZmZno0KEDXrx4genTp8PT0xNmZmZ49OgRAgICVC4qL+xll9TYsiuJ+0mdOnWwdOlStfNUqVJFqpefbU+b66mJvn374vjx4/jiiy9Qv359mJubQ6lUolOnThptC/ntj8IeZTphwgQEBwdj8uTJ8PHxgZWVFRQKBfr3718o27JCoVD7WeW1nxVHbJoojO2vQFnE7du3ZV/ed+7cgVKpzPWK8KysLD4+XnYa7t1sLuusgb29faFk/9kplUrcu3dP9sV469YtAJBi/+OPP2BsbIx9+/bByMhIqhccHKwSp1KpxLVr11C/fv1cl2tra4tt27ahRYsWaN++PY4ePYqKFStK011cXHDo0CGVYZB37tzReN1cXFxw4MABJCUlyc76ZI2qcHFxASD/HLJ7nzNCLi4uUCqVuHv3ruwv5Js3bxa4zdy4u7vj1atXeW4fLi4uCAsLw6tXr2QHq/eJK6sfs85MZvc+7bq7u0MIATc3t1wTt6z948qVK7muf35ObWdvM6czCpcvX8atW7ewfv16fPrpp1J5aGioRu1rsq/Y2NiobJdpaWl5jrTLT2ya9oudnR1MTU3VfqY3btyAnp6elDgUBRcXF1y7dg1CCFnM+TkmuLu74+LFi2jfvn2u663ptlcavHz5EmFhYZg9eza++eYbqVzd/ppTnxRVf2T92nDlypVc623btg3+/v5YsmSJVPbmzZs8b7So6bHBxsZG7c9CmnwHaBpbfo4/Li4uuHTpEpRKpewP+ne/uwpTgU4bZA0JzbJ8+XIAkF3R/66sDyXr+gvg7W/x7w5Z8/Pzg6WlJb7//nukp6ertJN9eGhBhrOvWLFC+r8QAitWrICBgQHat28P4G02qVAoZNlvVFSUbAQEAPTo0QN6enqYM2eOSqarLvOsXLkyDhw4gJSUFHTo0EE2tNLPzw/p6en45ZdfpDKlUqnSz7n56KOPkJmZKVs/AFi2bBkUCoX02VhaWsLW1lb2OQDATz/9pPGy3pXV9o8//igrzxpdVdj69u2LiIgI7Nu3T2VafHw8MjIyALztk4yMDKxatUqanpmZKW2v2Wm6LTk5OaF+/fpYv3697NRuaGgorl27VtBVQq9evaCvr4/Zs2erbD9CCGl7adiwIdzc3BAUFKRysMk+X9Z9gTS5K23Hjh1hYWGBefPmqdyWIqvNrL+ysi9DCCEbJpsTTfcVd3d3le1yzZo1ef4lmp/YNO0XfX19dOzYEX/++afsZ/CnT59i8+bNaNGixXv97J4XPz8/PHr0CLt27ZLK3rx5IztG5KVv37549OiR2nlSUlLw+vVrAJpve/lVHMPZ36VuWwDUH4ty2haKqj/s7OzQqlUrrF27FtHR0SrtZl+Hd5e7fPnyPPcDTY8N7u7uuHHjhuy79OLFizh27Fie66BpbPk5/nz00UeIiYmRjX7MyMjA8uXLYW5ujtatW+fZRn4V6IxPZGQkunXrhk6dOiEiIkIaxlyvXr0c5+nYsSOcnZ0xbNgwfPHFF9DX18fatWthZ2cn2wgsLS2xatUqDB48GA0bNkT//v2lOv/88w8++OAD6cs9v8PZjY2NsXfvXvj7+8Pb2xt79uzBP//8g6+++kr6Hbxz585YunQpOnXqhIEDB+LZs2dYuXIlPDw8ZL9he3h44Ouvv8a3336Lli1bolevXjAyMsLp06dRsWJFzJs3T2X5Hh4e2L9/P9q0aQM/Pz8cPHgQlpaW6NGjB5o2bYrPPvsMd+7cgaenJ3bt2oUXL14A0Cx77tq1K9q2bYuvv/4aUVFRqFevHvbv348///wTkydPll1/M3z4cMyfPx/Dhw9H48aN8e+//0pnvgqifv36GDBgAH766SckJCSgefPmCAsLy/GvU4VCgdatWxf4GTNffPEFdu3ahS5duiAgIACNGjXC69evcfnyZWzbtg1RUVGwtbVF165d8cEHH+DLL79EVFQUatasie3bt6v9LTo/29K8efPQuXNntGjRAkOHDsWLFy+wfPly1KpVC69evSrQOrm7u2Pu3LmYMWMGoqKi0KNHD1hYWCAyMhI7duzAyJEj8fnnn0NPTw+rVq1C165dUb9+fQwZMgROTk64ceMGrl69KiWDjRo1AvD2gnM/Pz/o6+ujf//+apdtaWmJZcuWYfjw4WjSpIl0f6GLFy8iOTkZ69evh6enJ9zd3fH555/j0aNHsLS0xB9//KHR7+qa7ivDhw/H6NGj0bt3b3To0AEXL17Evn37pKHGOclPbPnpl7lz5yI0NBQtWrTA2LFjUa5cOfz8889ITU1Ve98aTaxbtw5DhgxBcHBwrvcmGjVqFFasWIEBAwZg0qRJcHJywqZNm6SLzzU5JgwePBi///47Ro8ejUOHDuGDDz5AZmYmbty4gd9//126P5am215+5Xc4+7lz57Bx40aVcnd3d/j4+GjUhqWlJVq1aoWFCxciPT0dlSpVwv79+xEZGalSN2tb+Prrr9G/f38YGBiga9euRdYfwNs/Dlu0aIGGDRti5MiRcHNzQ1RUFP755x9cuHABANClSxds2LABVlZWqFmzJiIiInDgwAFUqFAh17Y1PTYMHToUS5cuhZ+fH4YNG4Znz55h9erVqFWrlsoF++/SNLb69etDX18fCxYsQEJCAoyMjNCuXTvY29urtDly5Ej8/PPPCAgIwNmzZ+Hq6opt27bh2LFjCAoKUrlutVBoPP5L/G+Y6rVr10SfPn2EhYWFsLGxEePHj1cZBvvuEFQhhDh79qzw9vYWhoaGwtnZWSxdulRlOHuWQ4cOCT8/P2FlZSWMjY2Fu7u7CAgIEGfOnJHq5Hc4u5mZmbh7967o2LGjMDU1FQ4ODiIwMFBliO2vv/4qqlWrJoyMjISnp6cIDg7OcTj+2rVrRYMGDYSRkZGwsbERrVu3FqGhobJ+yBrOnuXkyZPSkNis4YHPnz8XAwcOFBYWFsLKykoEBASIY8eOCQAiJCREZT3USUpKElOmTBEVK1YUBgYGolq1amLRokWyoYxCvB2SOGzYMGFlZSUsLCxE3759xbNnz3Icep19qL8QqrcgEOLtUNuJEyeKChUqCDMzM9G1a1fx4MEDlTaTkpIEANG/f3+165Cdur7L3s6MGTOEh4eHMDQ0FLa2tqJ58+Zi8eLFIi0tTaoXFxcnBg8eLCwtLYWVlZUYPHiwOH/+/HsNZxfi7S0IvLy8hJGRkahZs6bYvn27RsNBc+rT7O22aNFCmJmZCTMzM+Hp6SnGjRsnbt68Kat39OhR0aFDB2FhYSHMzMxE3bp1xfLly6XpGRkZYsKECcLOzk4oFAqNhhXv2rVLNG/eXJiYmAhLS0vRtGlTsWXLFmn6tWvXhK+vrzA3Nxe2trZixIgR4uLFiyp9WdB9JTMzU0yfPl3Y2toKU1NT4efnJ+7cuaPRcHZNY8utX979rIQQ4ty5c8LPz0+Ym5sLU1NT0bZtW3H8+HFZnaz94fTp07JydXEuX75c7S0a1Ll3757o3LmzMDExEXZ2duKzzz6Tbn1x4sQJqV7r1q1FrVq11LaRlpYmFixYIGrVqiX1e6NGjcTs2bNFQkKCrK4m215Oy1K37RfWcPbsn70mw9kfPnwoevbsKaytrYWVlZX4+OOPxePHj9V+vt9++62oVKmS0NPTUzmmvU9/5ObKlStSfMbGxqJGjRpi5syZ0vSXL1+KIUOGCFtbW2Fubi78/PzEjRs3NNoPhMj72CCEEBs3bhRVq1YVhoaGon79+mLfvn0aHb80jU0IIX755RdRtWpVoa+vL4tT3XD6p0+fSu0aGhqKOnXqqNxyJeuzXrRokUqfqvtsc6P4/zNRCbRz50707NkTR48exQcffKDtcArF7t270aVLF1y8eBF16tTRdjhExapv376IiorCqVOnCjR/UFAQpkyZgocPH6JSpUqFHB2RbmDiU0KkpKTIRgdkZmaiY8eOOHPmDGJiYkrV86lyk3Un6YKOkCMqrYQQcHBwwMaNGzW6W/W7x4Q3b96gQYMGyMzMfK+fpol0HR9SWkJMmDABKSkp8PHxQWpqKrZv347jx4/j+++/LzNJDwAsWrRI2yEQaYVCoVC5r0xuevXqBWdnZ9SvXx8JCQnYuHEjbty4gU2bNhVhlERlH8/4lBCbN2/GkiVLcOfOHbx58wYeHh4YM2aM7E7TRKQ7goKC8N///hdRUVHIzMxEzZo1MW3aNPTr10/boRGVakx8iIiISGcU7Pa/RERERKUQEx8iIiLSGby4WQuUSiUeP34MCwuLQntqLhERlS5CCCQlJaFixYoFfv4i5R8THy14/PhxkT7nh4iISo8HDx6gcuXK2g5DZzDx0YKsW3A/ePCgSJ/3ozNiY4Fsj+QAANy9C+TxqAMiovyKfR0L9x/lx5u7E+/C1iz/x5vExERUqVKlaB7LQDli4qMFWT9vWVpaMvEpDKmpqmUWFgD7logKWap+KmAsL7OwtIClWcGPN7zkoXjxR0UiIiLSGUx8iIiISGcw8SEiIiKdwWt8qFSKjo5GbGwsAKDcy5eoq+V4iIiodGDiQ6VOdHQ0PD29kJKSDACwBfD8nToPHz5EZTu7Yo+NiIhKNiY+VOrExsYiJSUZ3kMDYenkCqvXiUDQJFmdFy9egHfFICKidzHxoVLL0skV5Z1rwCLppbZDISKiUoIXNxMREZHOYOJDREREOoOJDxEREekMJj5ERESkM5j4EBERkc5g4kNEREQ6g4kPERER6QwmPkRERKQzmPgQERGRzmDiQ0RERDqDiQ8RERHpDCY+REREpDOY+BAREZHOYOJDREREOoOJDxEREekMJj5ERESkM8ppOwAq+6KjoxEbG6t2mq2tLZydnYs5IiIi0lVMfKhIRUdHw9PTCykpyWqnm5iY4saN60x+iIioWDDxoSIVGxuLlJRkeA8NhKWTq2xa4pMonFw7G7GxsUx8iIioWDDxoWJh6eSK8s41tB0GERHpOF7cTERERDqDiQ8RERHpDCY+73j06BEGDRqEChUqwMTEBHXq1MGZM2ek6UIIfPPNN3BycoKJiQl8fX1x+/ZtLUZMREREmmLik83Lly/xwQcfwMDAAHv27MG1a9ewZMkS2NjYSHUWLlyIH3/8EatXr8bJkydhZmYGPz8/vHnzRouRExERkSZ4cXM2CxYsQJUqVRAcHCyVubm5Sf8XQiAoKAj/+c9/0L17dwDAb7/9BgcHB+zcuRP9+/cv9piJiIhIczzjk82uXbvQuHFjfPzxx7C3t0eDBg3wyy+/SNMjIyMRExMDX19fqczKygre3t6IiIjIsd3U1FQkJibKXkRERFT8mPhkc+/ePaxatQrVqlXDvn37MGbMGEycOBHr168HAMTExAAAHBwcZPM5ODhI09SZN28erKyspFeVKlWKbiWIiIgoR0x8slEqlWjYsCG+//57NGjQACNHjsSIESOwevXq92p3xowZSEhIkF4PHjwopIiJiIgoP5j4ZOPk5ISaNWvKyry8vBAdHQ0AcHR0BAA8ffpUVufp06fSNHWMjIxgaWkpexEREVHx48XN2XzwwQe4efOmrOzWrVtwcXEB8PZCZ0dHR4SFhaF+/foAgMTERJw8eRJjxowp7nDLvJwebnr9+nUtRENERGUBE59spkyZgubNm+P7779H3759cerUKaxZswZr1qwBACgUCkyePBlz585FtWrV4ObmhpkzZ6JixYro0aOHdoMvY/J6uCkApKemFWNERERUFjDxyaZJkybYsWMHZsyYgTlz5sDNzQ1BQUH45JNPpDrTpk3D69evMXLkSMTHx6NFixbYu3cvjI2NtRh52ZPbw02fXI7AlV1rkJGRke92czqLBAC2trZ8WCoRURnHxOcdXbp0QZcuXXKcrlAoMGfOHMyZM6cYo9Jd6h5umvgkqkBt5XUWycTEFDduXGfyQ0RUhjHxIZ2R21mkxCdROLl2NmJjY5n4EBGVYUx8SOeoO4tERES6gcPZiYiISGcw8SEiIiKdwcSHiIiIdAYTHyIiItIZTHyIiIhIZzDxISIiIp3BxIeIiIh0BhMfIiIi0hlMfIiIiEhnMPEhIiIincHEh4iIiHQGEx8iIiLSGUx8iIiISGcw8SEiIiKdUWYSn6pVqyIuLk6lPD4+HlWrVtVCRERERFTSlJnEJyoqCpmZmSrlqampePTokRYiIiIiopKmnLYDeF+7du2S/r9v3z5YWVlJ7zMzMxEWFgZXV1ctREaaun79ukZlRERE76vUJz49evQAACgUCvj7+8umGRgYwNXVFUuWLNFCZJSXlIQ4AAoMGjQoxzrpqWnFFxAREZV5pT7xUSqVAAA3NzecPn0atra2Wo6INJWenARAoP7A6bBz85RNe3I5Ald2rUFGRoZ2giMiojKp1Cc+WSIjI7UdAhWQub0zyjvXkJUlPonSTjBERFSmlZnEBwDCwsIQFhaGZ8+eSWeCsqxdu1ZLUREREVFJUWYSn9mzZ2POnDlo3LgxnJycoFAotB0SERERlTBlJvFZvXo11q1bh8GDB2s7FCIiIiqhysx9fNLS0tC8eXNth0FEREQlWJk54zN8+HBs3rwZM2fO1HYoVALcunULGTY2sjLeG4iIiMpM4vPmzRusWbMGBw4cQN26dWFgYCCbvnTpUi1FRtowZuxYxOYwjfcGIiLSXWUm8bl06RLq168PALhy5YpsWkEvdJ4/fz5mzJiBSZMmISgoCMDbBOuzzz5DSEgIUlNT4efnh59++gkODg7vEz4Vsjq9JqCcV0NZGe8NREREZSbxOXToUKG2d/r0afz888+oW7eurHzKlCn4559/sHXrVlhZWWH8+PHo1asXjh07VqjLp/djZlsJZrw3EBERvaPMXNxcmF69eoVPPvkEv/zyC2yyXSeSkJCAX3/9FUuXLkW7du3QqFEjBAcH4/jx4zhx4oQWIyYiIiJNlJkzPm3bts31J62DBw9q3Na4cePQuXNn+Pr6Yu7cuVL52bNnkZ6eDl9fX6nM09MTzs7OiIiIQLNmzQoWPBERERWLMpP4ZF3fkyU9PR0XLlzAlStXVB5empuQkBCcO3cOp0+fVpkWExMDQ0NDWFtby8odHBwQExOTY5upqalITU2V3icmJmocDxERERWeMpP4LFu2TG35rFmz8OrVK43aePDgASZNmoTQ0FAYGxsXWmzz5s3D7NmzC609IiIiKpgyf43PoEGDNH5O19mzZ/Hs2TM0bNgQ5cqVQ7ly5XD48GH8+OOPKFeuHBwcHJCWlob4+HjZfE+fPoWjo2OO7c6YMQMJCQnS68GDB++zSkRERFRAZeaMT04iIiI0PnvTvn17XL58WVY2ZMgQeHp6Yvr06ahSpQoMDAwQFhaG3r17AwBu3ryJ6Oho+Pj45NiukZERjIyMCr4SREREVCjKTOLTq1cv2XshBJ48eYIzZ85ofDdnCwsL1K5dW1ZmZmaGChUqSOXDhg3D1KlTUb58eVhaWmLChAnw8fHhhc1ERESlQJlJfKysrGTv9fT0UKNGDcyZMwcdO3YstOUsW7YMenp66N27t+wGhkRERFTylZnEJzg4uEjaDQ8Pl703NjbGypUrsXLlyiJZHhERERWdMpP4ZDl79qz0MMpatWqhQYMGWo6IiIiISooyk/g8e/YM/fv3R3h4uHSfnfj4eLRt2xYhISGws7PTboBERESkdWVmOPuECROQlJSEq1ev4sWLF3jx4gWuXLmCxMRETJw4UdvhERERUQlQZs747N27FwcOHICXl5dUVrNmTaxcubJQL24mIiKi0qvMnPFRKpUwMDBQKTcwMIBSqdRCRERERFTSlJnEp127dpg0aRIeP34slT169AhTpkxB+/bttRgZERERlRRlJvFZsWIFEhMT4erqCnd3d7i7u8PNzQ2JiYlYvny5tsMjIiKiEqDMXONTpUoVnDt3DgcOHMCNGzcAAF5eXvD19dVyZERERFRSlPozPgcPHkTNmjWRmJgIhUKBDh06YMKECZgwYQKaNGmCWrVq4ciRI9oOk4iIiEqAUp/4BAUFYcSIEbC0tFSZZmVlhVGjRmHp0qVaiIyIiIhKmlKf+Fy8eBGdOnXKcXrHjh1x9uzZYoyIiIiISqpSn/g8ffpU7TD2LOXKlcPz58+LMSIiIiIqqUr9xc2VKlXClStX4OHhoXb6pUuX4OTkVMxRUVkTHR2N2NhYtdNsbW3h7OxczBEREVFBlPrE56OPPsLMmTPRqVMnGBsby6alpKQgMDAQXbp00VJ0VBZER0fD09MLKSnJaqebmJjixo3rTH6IiEqBUp/4/Oc//8H27dtRvXp1jB8/HjVq1AAA3LhxAytXrkRmZia+/vprLUdJpVlsbCxSUpLhPTQQlk6usmmJT6Jwcu1sxMbGMvEhIioFSn3i4+DggOPHj2PMmDGYMWMGhBAAAIVCAT8/P6xcuRIODg5ajpLKAksnV5R3rqHtMIiI6D2U+sQHAFxcXLB79268fPkSd+7cgRAC1apVg42NjbZDIyIiohKkTCQ+WWxsbNCkSRNth0FEREQlVKkfzk5ERESkKSY+REREpDPK1E9dVLR4LxsiIirtmPiQRngvGyIiKguY+JBGNLmXzZEjR+Dl5SWbdv369WKMkoiIKHdMfChf1N3LJiUhDoACgwYNynG+9NS0Io6MiIgob0x86L2lJycBEKg/cDrs3Dxl055cjsCVXWuQkZGhneDySd0ZKp61IiIqO5j4UKExt3dWORuU+CRKO8HkE89aERHpBiY+RChbZ62IiChnTHyIsinNZ62IiChvvIEhERER6QwmPtnMmzcPTZo0gYWFBezt7dGjRw/cvHlTVufNmzcYN24cKlSoAHNzc/Tu3RtPnz7VUsRERESUH0x8sjl8+DDGjRuHEydOIDQ0FOnp6ejYsSNev34t1ZkyZQr++usvbN26FYcPH8bjx4/Rq1cvLUZNREREmuI1Ptns3btX9n7dunWwt7fH2bNn0apVKyQkJODXX3/F5s2b0a5dOwBAcHAwvLy8cOLECTRr1kwbYRMREZGGeMYnFwkJCQCA8uXLAwDOnj2L9PR0+Pr6SnU8PT3h7OyMiIiIHNtJTU1FYmKi7EVERETFj4lPDpRKJSZPnowPPvgAtWvXBgDExMTA0NAQ1tbWsroODg6IiYnJsa158+bByspKelWpUqUoQyciIqIcMPHJwbhx43DlyhWEhIS8d1szZsxAQkKC9Hrw4EEhREhERET5xWt81Bg/fjz+/vtv/Pvvv6hcubJU7ujoiLS0NMTHx8vO+jx9+hSOjo45tmdkZAQjI6OiDJmIiIg0wDM+2QghMH78eOzYsQMHDx6Em5ubbHqjRo1gYGCAsLAwqezmzZuIjo6Gj49PcYdLRERE+cQzPtmMGzcOmzdvxp9//gkLCwvpuh0rKyuYmJjAysoKw4YNw9SpU1G+fHlYWlpiwoQJ8PHx4YguIiKiUoCJTzarVq0CALRp00ZWHhwcjICAAADAsmXLoKenh969eyM1NRV+fn746aefijlSIiIiKggmPtkIIfKsY2xsjJUrV2LlypXFEBEREREVJl7jQ0RERDqDiQ8RERHpDCY+REREpDOY+BAREZHOYOJDREREOoOJDxEREekMJj5ERESkM5j4EBERkc5g4kNEREQ6g3duJpno6GjExsaqlF+/fl0L0RARERUuJj4kiY6OhqenF1JSknOsk56aVowRERERFS4mPiSJjY1FSkoyvIcGwtLJVTbtyeUIXNm1BhkZGdoJjoiIqBAw8SEVlk6uKO9cQ1aW+CRKO8EQEREVIiY+RIUgp2ugUlNTYWRkpHaara0tnJ2dizIsIiJ6BxMfoveQkhAHQIFBgwapr6BQAEKonWRiYoobN64z+SEiKkZMfIjeQ3pyEgCB+gOnw87NUzYt67ooddMSn0Th5NrZiI2NZeJDRFSMmPgQFQJze+ccr4tSN42IiLSDiU8ZldP9eABeW0JERLqLiU8ZlNf9eHhtCRER6SomPmVQbvfjybq25MiRI/Dy8pJN492ZiYiorGPiU4apux9PnqOQwLszExFR2cXER8doMgqJd2cmIqKyiomPjsptFBIREVFZpaftAIiIiIiKCxMfIiIi0hn8qasUy+lePRydVfrxPkxEREWDiU8plde9egCOziqteB8mIqKiw8SnlMrtXj0cnVV6qDs7d/369Tzvw8RnfBERFQwTnwJauXIlFi1ahJiYGNSrVw/Lly9H06ZNiz0Odffq4eiskk+T+ymZlK/IZ3wRERUyJj4F8H//93+YOnUqVq9eDW9vbwQFBcHPzw83b96Evb29tsOjUoD3UyIi0g6O6iqApUuXYsSIERgyZAhq1qyJ1atXw9TUFGvXrtV2aFTKZN1PKfvLzNZJ22EREZVZPOOTT2lpaTh79ixmzJghlenp6cHX1xcRERFq50lNTUVqaqr0PiEhAQCQmJhY4DhevXoFAHhx/yYyUlNk0xKf3H+7nEe3YVBOUeanvXmdhHd7MuFJJF7ftixRcRbKtJhoAG8///fZfoioYJJeJwFv3ilLTIJRplG+28rah4UQhREaaUgh2OP58vjxY1SqVAnHjx+Hj4+PVD5t2jQcPnwYJ0+eVJln1qxZmD17dnGGSUREpcSDBw9QuXJlbYehM3jGpxjMmDEDU6dOld7Hx8fDxcUF0dHRsLKy0mJkJU9iYiKqVKmCBw8ewNLSMu8ZdAz7J2fsm5yxb3Knrf4RQiApKQkVK1YstmUSE598s7W1hb6+Pp4+fSorf/r0KRwdHdXOY2RkBCMj1dOgVlZWPAjlwNLSkn2TC/ZPztg3OWPf5E4b/cM/fosfL27OJ0NDQzRq1AhhYWFSmVKpRFhYmOynLyIiIip5eManAKZOnQp/f380btwYTZs2RVBQEF6/fo0hQ4ZoOzQiIiLKBROfAujXrx+eP3+Ob775BjExMahfvz727t0LBwcHjeY3MjJCYGCg2p+/dB37Jnfsn5yxb3LGvskd+0e3cFQXERER6Qxe40NEREQ6g4kPERER6QwmPkRERKQzmPgQERGRzmDiU0RevHiBTz75BJaWlrC2tsawYcOk52vl5M2bNxg3bhwqVKgAc3Nz9O7dW+VGiQqFQuUVEhJSlKvy3lauXAlXV1cYGxvD29sbp06dyrX+1q1b4enpCWNjY9SpUwe7d++WTRdC4JtvvoGTkxNMTEzg6+uL27dvF+UqFJnC7puAgACV7aNTp05FuQpFJj99c/XqVfTu3Ruurq5QKBQICgp67zZLusLun1mzZqlsO56enkW4BkUnP33zyy+/oGXLlrCxsYGNjQ18fX1V6pelYw4BEFQkOnXqJOrVqydOnDghjhw5Ijw8PMSAAQNynWf06NGiSpUqIiwsTJw5c0Y0a9ZMNG/eXFYHgAgODhZPnjyRXikpKUW5Ku8lJCREGBoairVr14qrV6+KESNGCGtra/H06VO19Y8dOyb09fXFwoULxbVr18R//vMfYWBgIC5fvizVmT9/vrCyshI7d+4UFy9eFN26dRNubm4luh/UKYq+8ff3F506dZJtHy9evCiuVSo0+e2bU6dOic8//1xs2bJFODo6imXLlr13myVZUfRPYGCgqFWrlmzbef78eRGvSeHLb98MHDhQrFy5Upw/f15cv35dBAQECCsrK/Hw4UOpTlk55tBbTHyKwLVr1wQAcfr0aalsz549QqFQiEePHqmdJz4+XhgYGIitW7dKZdevXxcAREREhFQGQOzYsaPIYi9sTZs2FePGjZPeZ2ZmiooVK4p58+aprd+3b1/RuXNnWZm3t7cYNWqUEEIIpVIpHB0dxaJFi6Tp8fHxwsjISGzZsqUI1qDoFHbfCPE28enevXuRxFuc8ts32bm4uKj9Yn+fNkuaouifwMBAUa9evUKMUjve93POyMgQFhYWYv369UKIsnXMobf4U1cRiIiIgLW1NRo3biyV+fr6Qk9PT+3T2wHg7NmzSE9Ph6+vr1Tm6ekJZ2dnREREyOqOGzcOtra2aNq0KdauXQtRQm/FlJaWhrNnz8rWSU9PD76+virrlCUiIkJWHwD8/Pyk+pGRkYiJiZHVsbKygre3d45tlkRF0TdZwsPDYW9vjxo1amDMmDGIi4sr/BUoQgXpG220qS1FuS63b99GxYoVUbVqVXzyySeIjo5+33CLVWH0TXJyMtLT01G+fHkAZeeYQ//DxKcIxMTEwN7eXlZWrlw5lC9fHjExMTnOY2hoCGtra1m5g4ODbJ45c+bg999/R2hoKHr37o2xY8di+fLlhb4OhSE2NhaZmZkqd7R+d52yi4mJybV+1r/5abMkKoq+AYBOnTrht99+Q1hYGBYsWIDDhw/jww8/RGZmZuGvRBEpSN9oo01tKap18fb2xrp167B3716sWrUKkZGRaNmyJZKSkt435GJTGH0zffp0VKxYUUp0ysoxh/6Hj6zIhy+//BILFizItc7169eLNIaZM2dK/2/QoAFev36NRYsWYeLEiUW6XCod+vfvL/2/Tp06qFu3Ltzd3REeHo727dtrMTIq6T788EPp/3Xr1oW3tzdcXFzw+++/Y9iwYVqMrPjMnz8fISEhCA8Ph7GxsbbDoSLCMz758Nlnn+H69eu5vqpWrQpHR0c8e/ZMNm9GRgZevHgBR0dHtW07OjoiLS0N8fHxsvKnT5/mOA/w9q+0hw8fIjU19b3Xr7DZ2tpCX19fZWRabuvk6OiYa/2sf/PTZklUFH2jTtWqVWFra4s7d+68f9DFpCB9o402taW41sXa2hrVq1fXmW1n8eLFmD9/Pvbv34+6detK5WXlmEP/w8QnH+zs7ODp6Znry9DQED4+PoiPj8fZs2eleQ8ePAilUglvb2+1bTdq1AgGBgYICwuTym7evIno6Gj4+PjkGNOFCxdgY2NTIh+uZ2hoiEaNGsnWSalUIiwsLMd18vHxkdUHgNDQUKm+m5sbHB0dZXUSExNx8uTJXPuppCmKvlHn4cOHiIuLg5OTU+EEXgwK0jfaaFNbimtdXr16hbt37+rEtrNw4UJ8++232Lt3r+zaTKDsHHMoG21fXV1WderUSTRo0ECcPHlSHD16VFSrVk02nP3hw4eiRo0a4uTJk1LZ6NGjhbOzszh48KA4c+aM8PHxET4+PtL0Xbt2iV9++UVcvnxZ3L59W/z000/C1NRUfPPNN8W6bvkREhIijIyMxLp168S1a9fEyJEjhbW1tYiJiRFCCDF48GDx5ZdfSvWPHTsmypUrJxYvXiyuX78uAgMD1Q5nt7a2Fn/++ae4dOmS6N69e6kcWlrYfZOUlCQ+//xzERERISIjI8WBAwdEw4YNRbVq1cSbN2+0so4Fld++SU1NFefPnxfnz58XTk5O4vPPPxfnz58Xt2/f1rjN0qQo+uezzz4T4eHhIjIyUhw7dkz4+voKW1tb8ezZs2Jfv/eR376ZP3++MDQ0FNu2bZMN5U9KSpLVKQvHHHqLiU8RiYuLEwMGDBDm5ubC0tJSDBkyRLYjRUZGCgDi0KFDUllKSooYO3assLGxEaampqJnz57iyZMn0vQ9e/aI+vXrC3Nzc2FmZibq1asnVq9eLTIzM4tz1fJt+fLlwtnZWRgaGoqmTZuKEydOSNNat24t/P39ZfV///13Ub16dWFoaChq1aol/vnnH9l0pVIpZs6cKRwcHISRkZFo3769uHnzZnGsSqErzL5JTk4WHTt2FHZ2dsLAwEC4uLiIESNGlMovdiHy1zdZ+9O7r9atW2vcZmlT2P3Tr18/4eTkJAwNDUWlSpVEv379xJ07d4pxjQpPfvrGxcVFbd8EBgZKdcrSMYeEUAhRQsdCExERERUyXuNDREREOoOJDxEREekMJj5ERESkM5j4EBERkc5g4kNEREQ6g4kPERER6QwmPkRERKQzmPgQUbFQKBTYuXPne7UREBCAHj16FEo8RKSbmPgQlRHPnz/HmDFj4OzsDCMjIzg6OsLPzw/Hjh3Tdmha4enpCSMjI8TExAAAwsPDoVAocn2Fh4dj3bp1aqfxad1EZUM5bQdARIWjd+/eSEtLw/r161G1alU8ffoUYWFhiIuL03Zoxe7o0aNISUlBnz59sH79ekyfPh3NmzfHkydPpDqTJk1CYmIigoODpbLy5csjKioKlpaWuHnzpqxNhUJRbPETUdHhGR+iMiA+Ph5HjhzBggUL0LZtW7i4uKBp06aYMWMGunXrJtVbunQp6tSpAzMzM1SpUgVjx47Fq1evpOnr1q2DtbU1/v77b9SoUQOmpqbo06cPkpOTsX79eri6usLGxgYTJ05EZmamNJ+rqyu+/fZbDBgwAGZmZqhUqRJWrlyZa8wPHjxA3759YW1tjfLly6N79+6IioqSpmdmZmLq1KmwtrZGhQoVMG3aNGj6hJ1ff/0VAwcOxODBg7F27VoAb5/c7ejoKL1MTEykM2NZL0NDQwBvk5zs5Y6OjnBwcNBo2URUsjHxISoDzM3NYW5ujp07dyI1NTXHenp6evjxxx9x9epVrF+/HgcPHsS0adNkdZKTk/Hjjz8iJCQEe/fuRXh4OHr27Indu3dj9+7d2LBhA37++Wds27ZNNt+iRYtQr149nD9/Hl9++SUmTZqE0NBQtXGkp6fDz88PFhYWOHLkCI4dOwZzc3N06tQJaWlpAIAlS5Zg3bp1WLt2LY4ePYoXL15gx44defZFUlIStm7dikGDBqFDhw5ISEjAkSNH8pyPiHSElh+SSkSFZNu2bcLGxkYYGxuL5s2bixkzZoiLFy/mOs/WrVtFhQoVpPfBwcECgOyp3KNGjRKmpqYiKSlJKvPz8xOjRo2S3ru4uIhOnTrJ2u7Xr5/48MMPpfcAxI4dO4QQQmzYsEHUqFFDKJVKaXpqaqowMTER+/btE0II4eTkJBYuXChNT09PF5UrVxbdu3fPdZ3WrFkj6tevL72fNGmSylPuhRDC399fbVtZfWBmZiZ7vbt+RFQ68YwPURnRu3dvPH78GLt27UKnTp0QHh6Ohg0bYt26dVKdAwcOoH379qhUqRIsLCwwePBgxMXFITk5WapjamoKd3d36b2DgwNcXV1hbm4uK3v27Jls+T4+Pirvr1+/rjbWixcv4s6dO7CwsJDOVpUvXx5v3rzB3bt3kZCQgCdPnsDb21uap1y5cmjcuHGe/bB27VoMGjRIej9o0CBs3boVSUlJec6bxcLCAhcuXJC9/vvf/2o8PxGVXLy4magMMTY2RocOHdChQwfMnDkTw4cPR2BgIAICAhAVFYUuXbpgzJgx+O6771C+fHkcPXoUw4YNQ1paGkxNTQEABgYGsjYVCoXaMqVSWeA4X716hUaNGmHTpk0q0+zs7Arc7rVr13DixAmcOnUK06dPl8ozMzMREhKCESNGaNSOnp4ePDw8ChwHEZVcPONDVIbVrFkTr1+/BgCcPXsWSqUSS5YsQbNmzVC9enU8fvy40JZ14sQJlfdeXl5q6zZs2BC3b9+Gvb09PDw8ZC8rKytYWVnByckJJ0+elObJyMjA2bNnc43h119/RatWrXDx4kXZ2ZqpU6fi119/ff+VJKJSj4kPURkQFxeHdu3aYePGjbh06RIiIyOxdetWLFy4EN27dwcAeHh4ID09HcuXL8e9e/ewYcMGrF69utBiOHbsGBYuXIhbt25h5cqV2Lp1KyZNmqS27ieffAJbW1t0794dR44cQWRkJMLDwzFx4kQ8fPgQwNvh5vPnz8fOnTtx48YNjB07FvHx8TkuPz09HRs2bMCAAQNQu3Zt2Wv48OE4efIkrl69qtG6CCEQExOj8nqfs1xEVDIw8SEqA8zNzeHt7Y1ly5ahVatWqF27NmbOnIkRI0ZgxYoVAIB69eph6dKlWLBgAWrXro1NmzZh3rx5hRbDZ599hjNnzqBBgwaYO3culi5dCj8/P7V1TU1N8e+//8LZ2Rm9evWCl5cXhg0bhjdv3sDS0lJqb/DgwfD394ePjw8sLCzQs2fPHJe/a9cuxMXFqa3j5eUFLy8vjc/6JCYmwsnJSeX17nVNRFT6KITQ8MYYREQ5cHV1xeTJkzF58mRth0JElCue8SEiIiKdwcSHiIiIdAZ/6iIiIiKdwTM+REREpDOY+BAREZHOYOJDREREOoOJDxEREekMJj5ERESkM5j4EBERkc5g4kNEREQ6g4kPERER6QwmPkRERKQz/h/1kif2Ln+mpwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "background_ates = []\n", + "for _ in range(500):\n", + " intv_data_X_1 = single_mediator_with_multiple_confounders_nuisances_discrete_example.generate_data(\n", + " num_samples=500, treatments={X: 1}\n", + " )\n", + " intv_data_X_0 = single_mediator_with_multiple_confounders_nuisances_discrete_example.generate_data(\n", + " num_samples=500, treatments={X: 0}\n", + " )\n", + " background_ate = intv_data_X_1.mean()[\"Y\"] - intv_data_X_0.mean()[\"Y\"]\n", + " background_ates.append(background_ate)\n", + "\n", + "fix, ax = plt.subplots(1,1, figsize=(5, 2))\n", + "sns.histplot(background_ates, ax=ax)\n", + "ax.axvline(true_ate, linewidth=3, color=\"red\")\n", + "ax.axvline(ATE_value, linewidth=3, color=\"green\")\n", + "ax.set_xlabel(\"Sampled ATE\")\n", + "ax.set_title(\"Sampled ATE\\nblue: background, red: direct calculation, green: Eliater calculation\")\n", + "plt.show()" + ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.6" + "pygments_lexer": "ipython3", + "version": "3.11.6" } }, "nbformat": 4,