{
"cells": [
{
"cell_type": "markdown",
"id": "ff51ab3d",
"metadata": {},
"source": [
"# Assessing Support and Overlap with OverRule\n",
"\n",
"## Motivation\n",
"\n",
"When performing backdoor adjustment, even if we correctly specify the causal graph and observe all the relevant variables $X$ for backdoor adjustment, we can only estimate causal effects for individuals with covariates $X$ that satisfy two conditions: \n",
"* Support: Simply put, we require that we observe similar individuals in our dataset. For discrete covariates $X$, this can be formalized as requiring that $$P(X) > 0$$\n",
"* Overlap: We require that there is some possibility of observing both the treatment and control for similar individuals. Formally, we require that $$1 > P(T = 1 \\mid X) > 0$$\n",
"\n",
"OverRule [1] is a method for learning Boolean Rule Sets that characterize the set of individuals who satisfy both conditions, and is demonstrated in this notebook on some simple examples to build intuition.\n",
"\n",
"## Acknowledgements\n",
"\n",
"Code for OverRule is adapted (with minimal modifications, but some simplifications) from https://github.com/clinicalml/overlap-code, under the MIT License.\n",
"\n",
"[1] Oberst, M., Johansson, F., Wei, D., Gao, T., Brat, G., Sontag, D., & Varshney, K. (2020). Characterization of\n",
"Overlap in Observational Studies. In S. Chiappa & R. Calandra (Eds.), Proceedings of the Twenty Third International\n",
"Conference on Artificial Intelligence and Statistics (Vol. 108, pp. 788–798). PMLR. https://arxiv.org/abs/1907.04138"
]
},
{
"cell_type": "markdown",
"id": "7d01b3fd",
"metadata": {},
"source": [
"# Table of contents \n",
"1. [Illustration on a simple 2D example](#2dexample)\n",
" 1. [Problem Data](#2d_problemdata)\n",
" 2. [Applying OverRule with default arguments](#2d_applyingOverRule)\n",
" 3. [Interpreting the output of OverRule](#2d_interpret_output)\n",
" 4. [Configuring OverRule](#2d_configuration)\n",
"2. [Illustration on Lalonde and PSID datasets](#lalonde\")"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "201dfed2",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"import dowhy.datasets\n",
"from dowhy import CausalModel\n",
"# Functional API\n",
"from dowhy.causal_refuters.assess_overlap import assess_support_and_overlap_overrule\n",
"\n",
"# Data classes to configure ruleset optimization\n",
"from dowhy.causal_refuters.assess_overlap_overrule import SupportConfig, OverlapConfig\n",
"\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"id": "c1b893ae",
"metadata": {},
"source": [
"# Illustration on a simple 2D example \n",
"\n",
"[Back to Table of Contents](#toc)"
]
},
{
"cell_type": "markdown",
"id": "fc102d2c",
"metadata": {},
"source": [
"## Problem Data \n",
"[Back to Table of Contents](#toc)\n",
"\n",
"In this example, we have a pair of binary covariates $X_1, X_2$, with simple violations of the above conditions:\n",
"* Support: There are no samples where $X_1 = X_2 = 1$, i.e., $P(X_1 = 1, X_2 = 1) = 0$ \n",
"* Overlap: Only individuals with $X_1 = 0, X_2 = 0$ have a chance of recieving both treatment $(T = 1)$ and control $(T = 0)$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "71a42a50",
"metadata": {},
"outputs": [],
"source": [
"test_data = pd.DataFrame(\n",
" np.array(\n",
" [\n",
" [0, 0, 1, 1],\n",
" [0, 0, 0, 0],\n",
" [0, 1, 1, 0],\n",
" [1, 0, 0, 0],\n",
" ]\n",
" * 50\n",
" ),\n",
" columns=[\"X1\", \"X2\", \"T\", \"Y\"],\n",
" )"
]
},
{
"cell_type": "markdown",
"id": "0c6e3f0b",
"metadata": {},
"source": [
"We can visualize these patterns as follows, where we **add some jitter** to the values to makes things easier to see"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "f9754054",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAGxCAYAAADxrrYqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMm0lEQVR4nO3deXhTZcI28PskTdK0NClQugBFllIQQXFAallksbQKoqgzg6MCgwijAirVUeoCoqPFHRVGXhh9HUWF0Rc3YKhQrAwKMhb4xIWllB1altKkpGmTNuf74zRp0mZrmzRpz/27rl7Qk3OS59DQO88uiKIogoiISKYUoS4AERFRKDEIiYhI1hiEREQkawxCIiKSNQYhERHJGoOQiIhkjUFIRESyxiAkIiJZYxASEZGsMQiJiEjWghqE27Ztw6RJk9C1a1cIgoDPP//c6/nr1q3D+PHj0aVLF+h0OqSnpyMvLy+YRSQiIpkLahCaTCZcddVVWL58uV/nb9u2DePHj8fGjRtRWFiIsWPHYtKkSdizZ08wi0lERDImtNai24Ig4LPPPsPkyZObdN0VV1yBKVOmYOHChX6db7PZcPr0acTExEAQhGaUlIiI2gNRFFFRUYGuXbtCofBc74toxTI1mc1mQ0VFBTp16uTxnOrqalRXVzu+P3XqFAYMGNAaxSMiojbgxIkT6N69u8fHwzoIX3nlFVy6dAl//OMfPZ6Tm5uLxYsXNzp+4sQJ6HS6YBaPiIjCmNFoRHJyMmJiYryeF7ZB+NFHH2Hx4sX44osvEB8f7/G8nJwcZGdnO76337hOp2MQEhGRz26ysAzCNWvW4N5778Unn3yCjIwMr+dqNBpoNJpWKhkREbU3YTeP8OOPP8aMGTPw8ccfY+LEiaEuDhERtXNBrRFeunQJRUVFju+PHDmCvXv3olOnTujRowdycnJw6tQpvP/++wCk5tDp06fjjTfeQFpaGkpKSgAAWq0Wer0+mEUlIiKZCur0iYKCAowdO7bR8enTp+O9997Dn//8Zxw9ehQFBQUAgDFjxuDbb7/1eL4/jEYj9Ho9DAYD+wiJKKzU1tbCarWGuhjthkqlglKp9Pi4v3nQavMIWwuDkIjCjSiKKCkpQXl5eaiL0u7ExsYiMTHR7YAYf/MgLAfLEBG1J/YQjI+PR1RUFBf7CABRFFFZWYmzZ88CAJKSkpr9XAxCIqIgqq2tdYRg586dQ12cdkWr1QIAzp49i/j4eK/NpN6E3ahRIqL2xN4nGBUVFeKStE/2f9eW9L0yCImIWgGbQ4MjEP+uDEIiIpI1BiEREckagzDcmcuBsmLg9B7pT3N5qEtERO2cIAhev5555pmgvK4oili4cCGSkpKg1WqRkZGBQ4cOBeW1nHHUaDgzngG+fBAo+rr+WN9MYNKbgK75Q4WJiLw5c+aM4+9r167FwoULceDAAcexDh06BOV1X3rpJbz55pv45z//iV69euHpp59GVlYWfv31V0RGRgblNQEGYfgylzcOQQA49DWw/iHg1pWANjYUJSOiEDCYrSivtMBgtkKvVSE2Sg29VhWU10pMTHT8Xa/XQxAEl2PBIIoili5diqeeegq33HILAOD9999HQkICPv/8c9xxxx1Be20GYbgylzUOQbuDeUBlGYOQSCZKjVVY8H8/4ZsD5xzHxvXrgtzbr0SCLng1paa67777sHr1aq/nXLp0ye3xI0eOoKSkxGXHIb1ej7S0NOzYsYNBKEtVBu+PV/t4nIjaBYPZ2igEAWDrgXPIWbcPr08ZHLSaYVM9++yzePTRR5t1rX2ThYSEBJfjCQkJjseChUEYriJ97Lbh63EiahfKKy2NQtBu6/6zuGiyhE0QxsfHe91IPVxx1Gi40naSBsa4k5olPU5E7Z7B7H3FlIqq8NnN4r777kOHDh28fnli74MsLS11OV5aWhr0/knWCMOVNlYaHbr+IalP0C41C7jpDfYPEsmEr9qeLkxqg0DLmkZ79eqFxMRE5OfnY/DgwQCk3SN++OEH3H///QEsZWMMwnCmS5JGh1aWSX2CGj0Q1YkhSCQjsVFqjOvXBVvdNI+O6x+P2Ch1CErlXkuaRgVBwMMPP4y//e1v6Nu3r2P6RNeuXTF58uTAFrQBBmG408Yy+IhkTK9VIff2K5Gzbh+27j/rOD6ufzxybxsUNv2DgfDYY4/BZDJh9uzZKC8vx8iRI7Fp06agziEEuDEvEVFQVVVV4ciRI+jVq1eLfqEbzFZcNFlQUWVFTKQKHaODN4+wLfH278uNeYmI2hG9VsXgCxKOGiUiIlljEBIRkawxCImISNYYhEREJGsMQiIikjWOGg1n5nJpF4oqg7S2qH1ZtYbHOM+QiKjZGIThytOmvNf9FfjgVsByqf4YN+olImo2No2GI2+b8n77IjB8ruux9Q9J1xARUZMxCN0xlwNlxcDpPdKfrR0y3jblLdoCJF/resy+US8RETUZg7Ah4xng/2YBb14NrBwj/blulnS8tfjalNda2fgYN+ologARBMHr1zPPPBOU1123bh0yMzPRuXNnCIKAvXv3BuV1GmIQOvPWJNmazY++Nt1VRzf9GiIiP505c8bxtXTpUuh0Opdjzd1qyReTyYSRI0fixRdfDMrze8LBMs68NUnamx9bOkLT3UjQhs9p35T3kJuypGQAx3e4HkvNAgSl1JTLkaRE7ZM/vzsCxHkjXL1eD0EQgr45LgBMnToVAHD06NGgv5YzBqEzX02SLW1+dB4Jqu4ADJ8H9B4DKCLq9hnsVL/tkqdNeUc9Ko0abXjs7eEcSUrUXnkaRR5m/8/vu+8+rF692us5ly5daqXS+I9B6MxX82JTmh8bfnpTRQMbH6sPwSkfADuWAQW59dc4v7E9bcoLAH/5j3RM3QE4sct1OgVQ35R768r6T4yt+GmSiALIV5eN8//zEGvJDvWhxCB05q1JcsAtUpiVFXsPE3O59PiGR4CizfXH+2YC194PFG+VaoI7lgFF+a7XNnxje9qU136srBj44gH39+LclNtGPk0SkRut0WUTIC3ZoT6UGITOPDVJDrgFyHwe+GKu9zAxngEOfwP88qn7kBNt0hzAHte61gSdNeWN7U9Tbhv6NElEbgS7yyaA2DTaXrhrklRHNw5BwDVMAClwhs9pHIJ2RVuA9LmAxeS9DP6+sf1pym1DnyaJyI1AdtkEGZtG25OGTZJlxb7DRIB0ztA/e39ua6X76Q/2wTM9rgVqa4BzBwFBkAbSaDt6aCL10pSbmiU9fvGI9/KE0adJInLDn//nYaKlTaNlZWU4fvw4Tp8+DQA4cOAAAGkUazBHrXIeoT/8aZqwn+Mu5JxFxQHHdwIp19cfsw+eObkLeP8W4J0MYPk1wKbHpSD7cp77Cf32ptzULNfjqVnATW9Ij7ehT5NE5IY//8/biS+//BJXX301Jk6cCAC44447cPXVV2PFihVBfV1BFEUxqK/QyoxGI/R6PQwGA3Q6XWCetKxYWmHGkwf3SPP4Lh6RgmXrC+5rkKlZwKS3AEulVIPcsgiIHyAdr7wg1f6O7wS+f8tpKkQWMDYHqLEAHXsCMQmNn9dc7n50qbkMMF8Eaq1AcYHr89rLwz5CoqCqqqrCkSNH0KtXL0RGRjb/idz9P+f/Xa//vv7mAZtG/eGraUKpAb562HVqBMQGo0azgBtfBiI09WGW9QLw1XzXgTMp10vXr50qhdahPGm06QeTPY/0bNiU626UaEqG6/O2w0+TRO2ap1Hk1GIMQn94m+A+8XXg3wvqQ8dySQqbkfOBzOcA42lAUAAndkqT3i8bLj2XSiuFYMOaY1E+AEEaXVqwRDpmX1vUPjjnxpel2p67KRyeRokWbQEUSmB2gVR75adJIiIADMImEIDLbwHS7peCSRUFVJwGbLXS3EBnlkuArQb4+knPcwUzn/e+w0S601ZLzv2OB/OkMnwwWfq+YS3R1yjRjMVAhNLvuyYiau8YhP4wl0sDVtwFTN8s19qbna+5giMf8f6a9lpg30xpgM0dH0mBeHwnUFtdf17D+YC+BvaUFQNr7uSEeiKiOkEdNbpt2zZMmjQJXbt2hSAI+Pzzz31eU1BQgN/97nfQaDRISUnBe++9F8wi+sdbLetQXuP9AQHfcwU1Hbw/ro6WQva6R4F3s6Twev8WaWRp5xSpL9LOeT9Cf3eu4Ia+REQAghyEJpMJV111FZYvX+7X+UeOHMHEiRMxduxY7N27Fw8//DDuvfde5OXl+b44mHzVskRb42O+plGIousUCmd9swBdd2DAzcAHt7mO9CzKBzbluO5SD9TPB7QP7HGn4c4V3NCXqNXYbG5+T1CLBeLfNahNozfeeCNuvPFGv89fsWIFevXqhVdffRUAcPnll2P79u14/fXXkZWV5ePqIPJVy4rtIQ2ccR5IYzzteaRp30xpAM0NS4C8p6RapeOxLGDSG4DVBHwxx/3rHfoauLbBGqP2Mnoa2JOSAaTPkQbyOOOEeqKgUqvVUCgUOH36NLp06QK1Wg1BEEJdrDZPFEVYLBacO3cOCoUCarW62c8VVn2EO3bsQEZGhsuxrKwsPPzwwx6vqa6uRnV1fZ+Z0WgMfMF8TZ+I7uK6LJu6gzSloke6tDpMo62UHgHeGS99P3wuMGq+NJJT3UEa2VlT1bRd6huuLuG8TJz5AlB9SRq1ap864YwT6omCSqFQoFevXjhz5oxjxRQKnKioKPTo0QMKRfMbOMMqCEtKSpCQ4DphPCEhAUajEWazGVqtttE1ubm5WLx4cXAL5q2WNepRwFoF6BIb7/Sg7iAF3chsabJ8pL5u26TbpOvtS6pVlgFRnaVa5CfTpbC6x0dzsL3p1dN8QPucI3MnYN2sNrE8E1F7pVar0aNHD9TU1KC2tjbUxWk3lEolIiIiWlzDDqsgbI6cnBxkZ2c7vjcajUhOTg78C6m0jadPnNgp7QXYc4TrwtvOcwoLlkhfqVnSBPovHvC8H2FKZv2k9+ICKWiLtjQuS2oWoOsGzNvjez6gc4gf/a4+fG21UpMuEbUKQRCgUqmgUqlCXRRqIKyCMDExEaWlpS7HSktLodPp3NYGAUCj0UCj0QS/cOYy33v/2Rfe9nTOqL9Kf/e0H2HR1wDqtmr6/i0pFBXKxk2rN73RtGkP9qZS+z6JnjYDJiKSobBadDs9PR35+a7hsHnzZqSnp4eoRE6asvC2J5oY6c8e13rfqin52voVajKfl2p+swukP29d2fzQ2vCo67JvAKdREJHsBbVGeOnSJRQVFTm+P3LkCPbu3YtOnTqhR48eyMnJwalTp/D+++8DkDZ1XLZsGR577DHcc8892Lp1K/71r39hw4YNwSymfwKxi4NKK9XAfM0xtA+EsVySaoSdevtXRm+c50I6b/lkMUn9jVUGLrlGRLIU1CD88ccfMXbsWMf39r686dOn47333sOZM2dw/Phxx+O9evXChg0bMH/+fLzxxhvo3r07/vGPf4R26oSdv3uCOZ/jHDi2WmnZtZteBwynvL+W80CYQA1mcWwT5aV/8mY2kRKR/HAbpiY9+Rn3C28799nZzzn6XX3gODeDXn4zMP5ZYONj7vsTUzKA7kOB03ua3hfojX0rqTE50uo07ppmuS0TEbUj/uYBg7Cp/NkTzFwOVBul6RCVF+rXCLXvB9j/JmDM48CWxa6jQlOzgBteBEQAUR52pW9JudfNkibVv3+L5/Pm7QE6B6AplogoxLgfYbD4syeY1Qysf6TBfoBO+wzuXw8M+4tU80ufWz8do1NvoONlvstgLpf6/KorgMhYwGaV/q6JARQqoKpc+rvzFk32aRRnf/X+3FxphohkhkEYaB73A2ywz6ClAvh+GTBcqB+0YjVL13sKWnM5YDoHlB8HlGpAqwc2zHdt5ky5XgrX/51Qv/ehvXlVlyQt3eas4cAZVbT3MhARtTMMwkDztlOF8z6DGp37QSue5vU13HV+TA5w8r9u5iI2CFznLZoAaUsn+4AeTwNnOLeQiGQkrOYRtgv+rBGakgHUVLufVN9wXp+5HLhQJDVpDp8jBaC6Q91cRDerzgD1cxGBxjtM2JtIU7M8T+zn3EIikhHWCAPN13zCqDhg9ONSM6inSfUH84BLpdLAmq/mu+9rtJq9v05NVf3fG/b76ZKASW9JtVdvmwdXlrGJlIjaPQZhoKmipQWz3Y0WTc0COvYEIjTA+UPen+dSKfDdG41rffamz/HPer8+JrH+7+7C2WqSplR4w4EzRCQDDMJAatiPB9TX4Arfk6ZGxNTtrhHd2ftzRcd5b/q88WXpud3VKlMypNGjgOdJ+VUG35sHc4smIpIB9hEGirfRoj+skILLefCJr53ka6rdP2ZnOicNvElx3b/RsQGv8aTnLZoAKeSO75TC1B1u0UREMsEgDBRvo0UP5jVeX9R50IqzvlnADbnSuqT2gTHuRMZIcxK7DwWmfg7c8ZH0Z/eh0vGOvb0v0K3tBJT+4iFMxwMTX2f/IBHJAptGA8Wf3Skasm+PdKkUMF8EIEr7EK4cK/UppmTUT8J33lm+bxag6gD0SJemSDSUmgV0iPe9T+GNLwGbHned2B/VGdAnA/puvu+ZiKgdYBAGSnN3p9DGSs2geU+4GRizBS5zAgEpHK+9T1q+LWMhoFC4rn2aMh6Y8Kp/tTldklQrtS8ZF3uZ741+iYjaGQZhoPi7O4U7VpOXgTGbgfGLgcQrpWXYTuwE/vsOkDgQKPkFuGYWkHZ//TJtJ3YCXz8pBZw/gebPknFERO0YgzBQ7H1+nnan8BY2vppVLx4F1twp/d0+GKbwPeCaexo3m9pd/wwDjojIDwzCQLL3+fnanaIhX82qnfsCs7ZKA2cEBSBEAOOeBlaOcR+CAOcAEhH5iUEYaM1pavTVrNohHtCmuh4vK/YcggDnABIR+YnTJ8KBp6kU3ppVvc1D5BxAIiK/cWPecOLPpr/OjGc890ly5wgikjluzNsWNbVZtbl9kkRE5MAgbOs4/YGIqEXYR0hERLLGICQiIlljEBIRkawxCImISNYYhEREJGsMQiIikjUGIRERyRqDkIiIZI1BSEREssYgJCIiWeMSa+2FuRwwl0mb/Ebqpd0nuPQaEZFPDML2wHgG+PJBoMhpP8O+mdLWTtyFgojIKzaNtnXm8sYhCEib/K5/SHqciIg8YhC2deayxiFodzBP2qKJiIg8YhC2dVUG749X+3iciEjmGIRtXaS+ZY8TEckcg7Ct03aSBsa4k5olPU5ERB4xCNs6baw0OjQ1y/V4ahZw0xucQkFE5AOnT7QHuiTg1pXSwJhqA6DRA1GcR0hE5A8GYXuhjWXwERE1A5tGiYhI1hiEREQka0EPwuXLl6Nnz56IjIxEWloadu3a5fX8pUuXol+/ftBqtUhOTsb8+fNRVVUV7GISEZFMBTUI165di+zsbCxatAi7d+/GVVddhaysLJw9e9bt+R999BEWLFiARYsW4bfffsM777yDtWvX4oknnghmMYmISMaCGoSvvfYaZs2ahRkzZmDAgAFYsWIFoqKi8O6777o9//vvv8eIESNw5513omfPnsjMzMSf/vQnn7VIIiKi5gpaEFosFhQWFiIjI6P+xRQKZGRkYMeOHW6vGT58OAoLCx3BV1xcjI0bN2LChAkeX6e6uhpGo9Hli4iIyF9Bmz5x/vx51NbWIiEhweV4QkIC9u/f7/aaO++8E+fPn8fIkSMhiiJqampw3333eW0azc3NxeLFiwNadiIiko+wGjVaUFCAF154AX//+9+xe/durFu3Dhs2bMBzzz3n8ZqcnBwYDAbH14kTJ1qxxERE1NYFrUYYFxcHpVKJ0tJSl+OlpaVITEx0e83TTz+NqVOn4t577wUADBo0CCaTCbNnz8aTTz4JhaJxbms0Gmg0msDfABERyULQaoRqtRpDhgxBfn6+45jNZkN+fj7S09PdXlNZWdko7JRKJQBAFMVgFZWIiGQsqEusZWdnY/r06Rg6dCiGDRuGpUuXwmQyYcaMGQCAadOmoVu3bsjNzQUATJo0Ca+99hquvvpqpKWloaioCE8//TQmTZrkCEQiIqJACmoQTpkyBefOncPChQtRUlKCwYMHY9OmTY4BNMePH3epAT711FMQBAFPPfUUTp06hS5dumDSpEl4/vnng1lMIiKSMUFsZ22ORqMRer0eBoMBOp0u1MUhIqIQ8TcPwmrUKBERUWtjEBIRkawxCImISNYYhEREJGsMQiIikjUGIRERyRqDkIiIZI1BSEREssYgJCIiWWMQEhGRrDEIiYhI1hiEREQkawxCIiKSNQYhERHJGoOQiIhkjUFIRESyxiAkIiJZYxASEZGsMQiJiEjWGIRERCRrDEIiIpI1BiEREckag5CIiGSNQUhERLLGICQiIlljEBIRkawxCImISNYYhEREJGsMQiIikjUGIRERyRqDkIiIZI1BSEREssYgJCIiWWMQEhGRrDEIiYhI1hiEREQkawxCIiKSNQYhERHJGoOQiIhkjUFIRESyxiAkIiJZC3oQLl++HD179kRkZCTS0tKwa9cur+eXl5djzpw5SEpKgkajQWpqKjZu3BjsYhIRkUxFBPPJ165di+zsbKxYsQJpaWlYunQpsrKycODAAcTHxzc632KxYPz48YiPj8enn36Kbt264dixY4iNjQ1mMYmISMYEURTFYD15WloarrnmGixbtgwAYLPZkJycjHnz5mHBggWNzl+xYgVefvll7N+/HyqVqlmvaTQaodfrYTAYoNPpWlR+IiJqu/zNg6A1jVosFhQWFiIjI6P+xRQKZGRkYMeOHW6v+fLLL5Geno45c+YgISEBAwcOxAsvvIDa2tpgFZOIiGQuaE2j58+fR21tLRISElyOJyQkYP/+/W6vKS4uxtatW3HXXXdh48aNKCoqwgMPPACr1YpFixa5vaa6uhrV1dWO741GY+BugoiI2r2wGjVqs9kQHx+PlStXYsiQIZgyZQqefPJJrFixwuM1ubm50Ov1jq/k5ORWLDEREbV1QQvCuLg4KJVKlJaWuhwvLS1FYmKi22uSkpKQmpoKpVLpOHb55ZejpKQEFovF7TU5OTkwGAyOrxMnTgTuJoiIqN0LWhCq1WoMGTIE+fn5jmM2mw35+flIT093e82IESNQVFQEm83mOHbw4EEkJSVBrVa7vUaj0UCn07l8ERER+SuoTaPZ2dlYtWoV/vnPf+K3337D/fffD5PJhBkzZgAApk2bhpycHMf5999/P8rKyvDQQw/h4MGD2LBhA1544QXMmTMnmMUkIiIZC+o8wilTpuDcuXNYuHAhSkpKMHjwYGzatMkxgOb48eNQKOqzODk5GXl5eZg/fz6uvPJKdOvWDQ899BAef/zxYBaTiIhkLKjzCEOB8wiJiAgIg3mEREREbQGDkIiIZI1BSEREssYgJCIiWWMQEhGRrDEIiYhI1hiEREQkawxCIiKSNQYhERHJGoOQiIhkjUFIRESyxiAkIiJZYxASEZGsMQiJiEjWGIRERCRrDEIiIpI1BiEREckag5CIiGSNQUhERLLGICQiIlljEBIRkawxCImISNYYhEREJGsMQiIikjUGIRERyRqDkIiIZI1BSEREssYgJCIiWWMQEhGRrDEIiYhI1hiEREQkawxCIiKSNQYhERHJGoOQiIhkjUFIRESyxiAkIiJZYxASEZGsMQiJiEjWGIRERCRrDEIiIpI1BiEREckag5CIiGQt6EG4fPly9OzZE5GRkUhLS8OuXbv8um7NmjUQBAGTJ08ObgGJiEjWghqEa9euRXZ2NhYtWoTdu3fjqquuQlZWFs6ePev1uqNHj+LRRx/FqFGjglk8IiKi4Abha6+9hlmzZmHGjBkYMGAAVqxYgaioKLz77rser6mtrcVdd92FxYsXo3fv3sEsHhERUfCC0GKxoLCwEBkZGfUvplAgIyMDO3bs8Hjds88+i/j4eMycOdOv16murobRaHT5IiIi8lfQgvD8+fOora1FQkKCy/GEhASUlJS4vWb79u145513sGrVKr9fJzc3F3q93vGVnJzconITEZG8hM2o0YqKCkydOhWrVq1CXFyc39fl5OTAYDA4vk6cOBHEUhIRUXsTEawnjouLg1KpRGlpqcvx0tJSJCYmNjr/8OHDOHr0KCZNmuQ4ZrPZpEJGRODAgQPo06dPo+s0Gg00Gk2AS09ERHIRtBqhWq3GkCFDkJ+f7zhms9mQn5+P9PT0Ruf3798f+/btw969ex1fN998M8aOHYu9e/eyyZOIiIIiaDVCAMjOzsb06dMxdOhQDBs2DEuXLoXJZMKMGTMAANOmTUO3bt2Qm5uLyMhIDBw40OX62NhYAGh0nIiIKFCCGoRTpkzBuXPnsHDhQpSUlGDw4MHYtGmTYwDN8ePHoVCETTclERHJkCCKohjqQgSS0WiEXq+HwWCATqcLdXGIiChE/M0DVseIiEjWGIRERCRrDEIiIpI1BiEREckag5CIiGSNQUhERLLGICQiIlljEBIRkawxCImISNYYhEREJGsMQiIikjUGIRERyRqDkIiIZI1BSEREssYgJCIiWWMQEhGRrDEIiYhI1hiEREQkawxCIiKSNQYhERHJGoOQiIhkjUFIRESyxiAkIiJZYxASEZGsMQiJiEjWGIRERCRrDEIiIpI1BiEREckag5CIiGSNQUhERLLGICQiIlljEBIRkawxCImISNYYhEREJGsMQiIikjUGIRERyRqDkIiIZI1BSEREssYgJCIiWWMQEhGRrDEIiYhI1hiEREQka0EPwuXLl6Nnz56IjIxEWloadu3a5fHcVatWYdSoUejYsSM6duyIjIwMr+cTERG1VFCDcO3atcjOzsaiRYuwe/duXHXVVcjKysLZs2fdnl9QUIA//elP+Oabb7Bjxw4kJycjMzMTp06dCmYxiYhIxgRRFMVgPXlaWhquueYaLFu2DABgs9mQnJyMefPmYcGCBT6vr62tRceOHbFs2TJMmzbNr9c0Go3Q6/UwGAzQ6XQtKj8REbVd/uZB0GqEFosFhYWFyMjIqH8xhQIZGRnYsWOHX89RWVkJq9WKTp06eTynuroaRqPR5YuIiMhfQQvC8+fPo7a2FgkJCS7HExISUFJS4tdzPP744+jatatLmDaUm5sLvV7v+EpOTm5RuYmISF7CdtTokiVLsGbNGnz22WeIjIz0eF5OTg4MBoPj68SJE61YSiIiausigvXEcXFxUCqVKC0tdTleWlqKxMREr9e+8sorWLJkCbZs2YIrr7zS67kajQYajabF5SUiInkKWo1QrVZjyJAhyM/Pdxyz2WzIz89Henq6x+teeuklPPfcc9i0aROGDh0arOIREREBCGKNEACys7Mxffp0DB06FMOGDcPSpUthMpkwY8YMAMC0adPQrVs35ObmAgBefPFFLFy4EB999BF69uzp6Evs0KEDOnToEMyiEhGRTAU1CKdMmYJz585h4cKFKCkpweDBg7Fp0ybHAJrjx49DoaivlL799tuwWCz4/e9/7/I8ixYtwjPPPBPMohIRkUwFdR5hKHAeIRERAWEwj5CIiKgtYBASEZGsBbWPkIiI5MlgtqK80gKD2Qq9VoXYKDX0WlWoi+UWg5CIiAKq1FiFBf/3E745cM5xbFy/Lsi9/Uok6DwvkBIqbBolIqKAMZitjUIQALYeOIecdftgMFtDVDLPGIRERBQw5ZWWRiFot3X/WVw0WVq5RL4xCImIKGB81fgqqlgjJCKidszXgBhdGA6YYRASEVHAxEapMa5fF7ePjesfj9gotcdrDWYrjl0w4aeT5Th2wdRq/YkcNUpERAGj16qQe/uVyFm3D1v3n3UcH9c/Hrm3DfJYYwzlSFMusUZERAFnMFtx0WRBRZUVMZEqdIz2PI/QYLbi4TV73A6yGdc/Hq9PGdysOYj+5gFrhEREMhTsCe96rcrv5/NnpGkwJ+MzCImIZKYpzZCtsUJMqEeaMgiJiGTE14R352bI1uq3C/VIU44aJSKSEX8nvLfmCjEtGWkaCKwREhHJiL/NkOWVFuw6Uob5GX0x5LJOMFtroFVFoPBYGVZuKw5ov11zR5oGCoMwjLSl1dqJqG3ytxmyosqKt+8egne2H8HrWw45Hh+d2gVv3z0EpurA9tsl6CLx+pTBfo80DSQGYZhoa6u1E1HbZG+G3OphqoK9GTJWq8bLeQfx7UHX8749eA4CgL9NHhjwsjVlpGkgsY8wDLTF1dqJqG2yN0OO6x/vcrxhM6S11tYoBO0KDp6DpdYW9LK2FtYIw4Dfc2jM5YC5DKgyAJF6QNsJ0Ma2almJqO3zpxmyorrG63OYfDzeljAIW5OHIPNV4zNVWwHjGeDLB4Gir+sf6JsJTHoT0CUFt9xE1O74aoYM9ZSG1sQgDDRPtTYvQRar9b4UXHKUtfG1AHDoa4jrH8L58W+hUtGBg2uIKGD87UtsD9hHGEjGM8D/zQLevBpYOUb6c90swHAS2PiY2yDD+oeQoK72Oocm2lbR+No6wsE8mA3nMPrlAsxfswelxqoA3xQRyZG/fYntAWuEgWIu91hrw4ZsoOvVwP4vG193MA8aS7nXOTTKS795fWmVtQKA+5UhiIiaKxRTGkIxjYxBGCjmMo+1NhzMA9Lu93xttQEJnXu7fcMBQI1K5/UHVavWASgD0DoL1BKRfLTmlIZQTSNj02igVBm8P26t9PxYpB6A9IbrGReNQd1j0TMuGlXWWjy8Zg/e/8mImj6Zbi+tTcmCENUJ0Wql41iwF6glIgq0UE4jY40wUOrCzBMxqjMEdw+kZkkDapyZy2GrMiDOUolV1wPmCC2sV70KYctTUO7/AlB3AIbPg9h7LGxQQGM1YPXd/XD36gMwWWrb1WguIpKHUG7FxCAMFG0niH0zIRxy0zyakgFExkJMzYJwMK/+eGoWcNMbrnMBjSWA4TgU3y4BivIBADEAxL6ZKB/3EkzDnkBCjAoReQsgFORCBSAOQKeUTKyb+jze+KGiXY3mIqLWFaqlHkO5FRODMFC0sRAnvAJhQzZQtKX+eEoGkD4Hwvs3o/aer6HMWgJUGwCNHohqMCHeXA4c3gr88qkjBO2EQ18jFgJKhj2PpE2PQzi82eVxRdHXSIGAlyf9HdGsERJRMzS3jy4Q4RnKeYsMwgCymS5C0X0okD5X6hNURQEndgJrpwKWSxAry4Duv/P8BOYyIDYZ6D4MGP4gYDEB6mjg+E7g+7cgHMpD/8znGoWgnbIoD9E1Rkh1RCJqi0JZI/N3n0JngRrgEsp5iwzCQDGXQ6mJAhKvBAQBKNkHfP8WYLnkOEXwtRxadQUQkwCc3AUU5NYfT7kemPKBFKjmiz6ew8egHSIKW6FcfL85fXTNDU93QrkVE4MwEOpWjRGcp084h5flEsTULCijGg+KcVmFJlIPbHi0UbOo9L0ADJ8LIdL7KjS1Gj1OXjBxKyeiNiYQodKS2mRT+ujsr2OpsQV0gEuotmJiELaUp4n0TuElnt4DoW5QjMFshaHSgnjhIjQb57uEp3jfdghF7ps9UbQFuO5RwHAaYkqma+jarx9wC05eEqCoOYouFiNqqnXYdVSJq/r2RLybT5Pc/5AofLR01GRLa5O+/u9HqSPw86ly6LVqLPziZ3xz4BxWTRvi9ZrmDHAJxVZMDMKW8jaRvmgLkPUChLT7AW2s4406qocK008/36ivTzCc9PFiAvDFAxBmbgb+/Vdpon4dccAtqLn+OXTf+AiUTs/bLSUT1T1fx7mKOHSJqf/PwP0PicJLS0ZNBqI26a2PbkxqF3z102molAJ2Hb2Ib+vO0aq8R4heq2oTH7g5ob6l/JlIX1cTtL9Rx/dU14eVugMwJgeY9gUQ28P7c9XWABnPADYrMPpxYM5/gbk/AvP2wpL5IhQNQhCQRpOqNmZjz8GjOF1uBsD9D4nCUUtGTfpTm/Tn9d2tLTomtQvuGdkLH+48hozLExwhCACFx8owOtX9OskTBiYiQqnAw2v2YPTLBbh52XdhuyYyg7ClfEyktz/u/EaNsBilx9QdpH7Ek7uA928BfvtK6lt0p/9NQHRn4OfPgLeGAKvGAcuvAfKeAFSRECymRiFoF3E4D/31NXhi3T6cMZhhCMB/GiIKLHuNzB1foyabWps0mK04dsGEn06W41jdmAKgvo+u4NEx+GruCGzJvg6De8Ti/tWFuPvay1BirHZ5npXbijFzZC+MaRCGEwYm4qmbBuDJz/a1iQ/cbBptKW0naTsldxPpnVaNcf6h16jrBrwMnwfsWFY/OOb7t6RghOAyF1HsmwXcsATC+myPO1hEXPe412KqrBUoOFiGUxfNjjVMPeESbUStryWjJptSm/TVLeLcR/fTyXIs3XIIADDksgaD/ersO1WO7MxUzB2XAp1WhSiVEkqFgFPl5pCtFNNUDMKW0sZKm+Ouf8ilz67hqjHOP/DNRy2Y1icTET2udZ0mYbkkjTIdPhdIn4saVQeUWLVI6KRHRMVpr4t6C+Of81rMmrqFuS9WWhAbJZ8NN4nakuaOmvR3Dl5T+xKd/2621uDX00aMTu2Cbw+eQ7RaibfvHoJ3th/By3kH61+vXxfMGdcXRh8fqMPpAzeDMBB0ScCtK4HKMo+rxji/UV/9thSjp7+A3rbjjdcftVwCCpYAAM5N2YRfKyPQdUc2hKF/9l6GWovUrNpw6gWAmj5ZKDheg4Xju2JU5wooqk/hv/f3xldFVXj121KYLLWOc9vbhptEbU1zRk3qtSosuf1K/OfQeSToImG21kCrikCpsQqj+sY5nq+pI1Odf29pVRFYua0YK6cOxaSrknBFVz1e3LQf3x5sHKo2EXj8xv5eyxxOH7gZhIGijXVdLq0B52aPH4ovwBqVADNsiPLylLGd4jAmToDiq6+B4XO8vrxgLpdWtHHTrHr02ueRJgJ9duZA+R+pH7ELgOkpmRgx9Xnc9kExTJbadrnhJpFciADW/3S6UZPnyL71K001tS/R+fdW4bEyjEyJgwgRX/10Bok6LQo8hGrBwXPIzkx11B4bCrcP3AzCVmRv9jCarSg1VmHN4UqpifRw4yZPW98sfLzvEm7sVoUkQFpmzUONDykZwLHtwPfLHM2qjiXeYpLQ2aaCbsujUDaYo6gs+hqpCgHfznsLl4QOrTJxlYgCz98mz1itCvMz+mLIZZ0ctcbCY2VYuU36MOxpusPrUwajvNKC23/XHU9+/jO+PXgOd6d5H+VeaqzGzJG9IEAKRrtw/MDNIPSm4cov2k5ea33+0GtVKK+0wFhlhUqhgC3zbxAr7oOgUDrWFK1JHoGqzFfx6lu/YPz0ZOlCLwNphGvvc6xgY29WtRPu+AixUXEeJ+oLB/MQl2VCXOeEFt0XEYWOv02eEUoFdh8vx+t1A2AAYHRqF7x99xB8tvukY7qDu4E0l3WOxrELJkcNz9ccQq1Kib988CNmXdcbM0f1QpW1Fh2j1OjROQrxdXOaw2WOYdCDcPny5Xj55ZdRUlKCq666Cm+99RaGDRvm8fxPPvkETz/9NI4ePYq+ffvixRdfxIQJE4JdzMbqlk1zGaDSN1MaGKNL8nqprx9utGjCmPgqCD8+B+Hb+oASUzJRM3s7Pt5nRF9DFEyW2vqBNYe/dhlIA2slxKjOsEUnQPk/I13WNHWhjoZQed77vXJ9UqI2zZ8mT4PZiic/29eoqfLbg+cgAHjhtkFepzssuX2QS9OpfQ6hp6bPEmMVTJZax6hTe03QHoLhtKhHUINw7dq1yM7OxooVK5CWloalS5ciKysLBw4cQHx8fKPzv//+e/zpT39Cbm4ubrrpJnz00UeYPHkydu/ejYEDBwazqK48LZtWN1UBt670WDP0+cM1nkHnE99AcLfVUtHXiMgTMHbki9BFnMeO6Z1Qq4lExYCXoBMEKIvy6mt8KRkQ0udAOLELYo/hbpdcQ0oGcHwH0CPd+/36mgtJRGHNV5NnlDoCZZeqPdYaCw6eg9lS67VWefxCJRJ0kYhWK2Gy1GLltmK8ffcQj02fkSolCh4d43b0ayAX6w4EQRRFMVhPnpaWhmuuuQbLli0DANhsNiQnJ2PevHlYsGBBo/OnTJkCk8mE9evXO45de+21GDx4MFasWOHXaxqNRuj1ehgMBuh03heo9qisGHjzas+Pz9sDdO7d6LDBbG3UrGA3rn883r69D9Rf/AXC8DnSBHoPxHvyILyb5fi+pt/NqBn/PDSVpyFUXqjf3ul76d9VvPsz4D+vQDjkNH2jbh9EqRY5DzhV6Hmuo5dgJ6Lwd7rcjJx1rrW90aldMHNkL3y08xiuH5CAHp2jMOV/dnp8jnX3D8dtb3/v8fFV04bg410nMGFgIh799CcAQLRaiVnX9caQyzpCFIHuHbXo3EHjNsScW8qi1RG4/rVvPb5WwaNj0DMu2p9b98rfPAjayjIWiwWFhYXIyMiofzGFAhkZGdixY4fba3bs2OFyPgBkZWV5PB8AqqurYTQaXb5azNeyaR6aEr210/9QfAHKqotSzc1i8vr0QuUFl+8jDnwJTeUZCKtvl7Z3EgRpu6c7PgSGz4Pwr7tRc+PLEO/JA+74CJj6OXDZcOD0Humc7sMgZr0AcfLb0mo2dg3mOhJR2+OtyfPd7UfwYEYqnvnyF8T6qGH5qoFpVRHYuv8srr6so2MZNnvT5/9+dxT9EmPQu0sHj3sWOi+1VnzeQ1dOndaeYxi0ptHz58+jtrYWCQmugzASEhKwf/9+t9eUlJS4Pb+kpMTj6+Tm5mLx4sUtL7AzP5dNa6hhO320WonZ1/XGkMs6QakAbOZj0gNqH590nB9Xd5DCLlIHzP4G2PRE470Kb10BmA0wR3VH5H9eg+LYd9LAmh3LgPxnAQACAFtKJqyzt6PswlkoI3WIiu2CaF1cUAYFEVHr8PYBvODgOcwa1RsmSy0sNTavfXox2givi27/eKwMAFBZXePXpH97DbC80gprrQ2Dk2Ox60gZTJZanwNtWnuOYZsfNZqTk4Ps7GzH90ajEcnJyS17Uj+XTWvIZUCM06oLr285hGi1Ejsf6Av1tC8AQQHckwcUFzTavNfRrwfUr0W6YxmgVAMdEqTBMkNnuOxcDwioyXwF17z5Mx65Lge/v16HmC2PQ2jQB6ko+hqKfwvYmJSDZ7cUY8LASrwxsRoRGx527WP0c1AQEYWer4EyldYaAMClqlrMHNkLSgEuYTeufzwW33wFzldUYfEtA4Evf3Fd4q1fPJ6edDmOnDdh9cw06LUqnCgzISZShR6do/3etd4+OvX+1YU+B9q09hzDoAVhXFwclEolSktLXY6XlpYiMTHR7TWJiYlNOh8ANBoNNBpNywvszM9l0xpyXoVh9nW98c72I46liNZN7Y3oLY8DztMYUjJcNu9FahYw7mlg9e3S4/a1SI//ANz4EvDvx1wH2Dht/quCFbOv641XtxVjfC81dB6mS0QczkNG+kI8C2BokqJxCAJ+DQoiovDQ8AO4vRXKPmimcwc1bhiQgB+OXsCqbcX49P7hmDGyF6qsNiR31GLfKQNuWLoNALBy6lDcODARM0b0dHn8pje3O1agGtOvC+4Z0QtT/mcn0np1cgwEtNcAa2wi/rb+10a1VPvo1FnX9fY50Kbd7EeoVqsxZMgQ5OfnY/LkyQCkwTL5+fmYO3eu22vS09ORn5+Phx9+2HFs8+bNSE/3MeoxGPxYNq0h51UYhlzWyTFX55HRCUjZ+QQUDXeHKNoCCErgnn9Lfx74N8SCJRBuXSGFo30t0jE5wL8XeN25/nzZBfx2RoXvHx6CGNNRr7emslYAAMb3VEP4j+f1S1FZxiAkCnP2D+A/HClzaYWyG9uvCxZNugIlhip076jF1v2leDnvIOZn9MXqnccctbL5GX2x8j/FLt87P25XcOAcIEqBtnTLIeSs24fnbx2IJ9ZJUy9Wz0zz2lQ7c1QvLLXU4v7VhZh1XW88ddPlMFtqW203eneC2jSanZ2N6dOnY+jQoRg2bBiWLl0Kk8mEGTNmAACmTZuGbt26ITdX6vN66KGHMHr0aLz66quYOHEi1qxZgx9//BErV64MZjE987FsmjuRKiUmDExEpKp+HNL4nmrH0maNHMoDxuYA790IWC5BACDWVEMYPrd+UE3DxbmdFW2RFuiGDvOvVaHDxvuhGO7+g4adfQFux3ZQnnB+IVHYs38A337ovKMVytk3B85BFH/B4B6x+H8nyjFzZC9Eq5UuH9YBOL631ypvGJiEAV11mDWqt8tUDKA+0ABpaoXRXIM703pg5sjeUHoYgml/3i4dNFg1bYhjLdQuMZEhX2UmqEE4ZcoUnDt3DgsXLkRJSQkGDx6MTZs2OQbEHD9+HApF/b/a8OHD8dFHH+Gpp57CE088gb59++Lzzz9v3TmELVReacGjn/6E1TPTHMd8Bk6tFfjTGqDaKE2AP74TtQNug2irlX5APkaZiqINFaIaqTtzpD0Juw/1ugD35iPSfoMRUbHey8X5hURtQoIuElf3iMUjn/w/t487amJbDkGsq83VNpg5Z7bWNBrbYOfcv2cPwypr/WL9x8tMmPV+IQDg0/sat+B5et6xdWuhtusgBIC5c+d6bAotKChodOwPf/gD/vCHPwS5VD60YBSlvePauTPYsf+gJ6JNWivUMVn+etiu+AP+72AVpqRkQvA1ylTXFX2UaihP1A2y+fFd4O5PASFCqnHWqemThcPX/g2vflCMGwYkQButk+YsVl5wHXxj76/0MCiIiMJLqbEKJy6avZ5jD66Cg+fw1E2XQym47n2jVUW4jG1w5ty/Z18pxnnkp/Pfvys632ggjKfn/SZEE+gbavOjRgPOzdJqYkomrBOXwhwZ7/OHZX/cuTN481ELZqRkul/9pf9NQPUl4PJbpLmBdYEUsXUxrvrdQqDnX4ELhzwvuN03E8L+DVCd2CWNKP1+mTSdomAJ0O1q4Nr7HcuxCbpk6KqUKJybApVSAcW/GwyUSckE/vItcGo30PM69g8StQH2VVpmjmy8yIcz57AyW2rRo3O0y3SJwmNlGJES51Jjc+bcHOo8ncL570D97z6lIGDrAWn0acNmWGfhsElv0CbUt0kellYTir6GYv187Pr1MEqNVV6fwt5xbarrDB7cIxb9L0uGeOOL0ihRZ/1vAjIWAz/8HXg7HVhzp7TizMldEIb+GX31NgirbwUulQA3vAT0zXK9PiVDCrrtr0uryiRfWz/SdP8GKQw/mAysuRPC6tuhqDiFxC3zEPnrp1D++9HGwVz0NfDvxwGbFdKmLkQU7uzzCO2tUO40DCtd3Z6Hubdf6Zgcv3Jbsc/XqrLWYkxqF9wzshdWbSvG2H71f7ez/+6bP74vPr0vHV/MGYFojdLr84Z6k17WCJ2ZyzzuAh9xOA/90hf6rMY7jxzduv+soxnhxzl9Edd9qOsWSbZaYNNjHkeDKjKfl5op858F/vMaMGMjcP1CaWWZ6gogQgMc3uq4TBRtEDwNrBk+D8K3S6TnHj7H5+AbTp8gahvs3TEepyT064I/j+iF+1dLfXjO8/TsW8PZJ8dHqb1HQs/O0Xjsxv64eMmCL+aOxKafz7j0G9qZLLWorrUhuVMUEnSROHbB+ziHUG/SyyB05mNpNZW1Alv3l/msxjd8c8VEqhCjqZaaHJ23SZr2hfvmTgAo2gIh6wVp6sT3b0nH1B3q5hLWb8OElPHAlA+BtXfBpk+GwnCi8a73gOvIUx+Db2Ct5PQJojbC/rvI5DQlwb7tUaRKibgOGvz+7e89br6tr6sdAlKoelpdZlz/eCgUAqKUCnTrLg2k233sYqMQtJ/bo1P9dkvqCEVYTaBviEHozMcoSfu0A3+q8c5vLkl040n6vtYcvVAE8eR/IUz5ACj5uW4u4RbXk4o2A4ICtgmvwaiIRYxecP9DdX4tf5d44/QJorDnvJCH87ZHQP2qMWtmX+vXPL2GLVrOz5N726BG2yN5O9d5z8Fnv/zF7Sa9Y1K74NlbrmjpP0GLMQideVlazXnaQbOr8XWT9GtNF2C4eB46Xaz3H4A6GkLRFogQIIxfDGx+2v15h/Igjn8OH/9kQIquFuPc7XrvHH6+dru3L/HG6RNEYc+f8EruFOX387lr0Wq4hZJ9F4lYrQpLbh+Eyupaj+uOllda8O9fSrHt0PlGtdXCYxdx/EIl/vGf4pDsQ2jHIHRWt7SauP4hCAfdTztocTVeG4uTlSqMfmc/Fo5X12+625BTIAlFmyGOftx9k2cdscqI5d9cwMqpQ1Ey+kV0VQhQOE2dsBlOQ7CPXLXvdi8oXUPfeesmTp8gajN8hVdTNW7Rknjbb9V52yTnsKy1SQPvGtZW7VZNGxKyfQjtGIQN6ZIg3LoSNaYLMF48j0tCNLYcteLVD4qR1rtzQNbBs3duv/ptKUZMfR59ICDisIe9BO00HeCNIlIH4AJmf/AjHhjbB78f/yaq085DZa2AVRWDghO1SL/2GvQS615r7VRg5HyI45+FYDwtLQR+Yqd0vOcIbs9E1MZ4Cq9A8Xcz3YZh6by4iDv2aR2hnEbBIHRHG4sIbSyU0T1gM1lwTX8rNgxODdg6eB000j+7yVKL2z4oxiPX5WBa5rOIuHi4ftNd+0LcdURlJAQvzZnC6R+x58ERWPOzATtPGnHjwCSM+0dh3QnSsOlotRKPXJeDjPSFSFBXo9SiQayqC3TdkuoHxgya4nNNVSJqv5xrc3qtCrFR0u89b9s92UMMQKOwdLfThH25tREpcbDU2rB6ZhoKj5XBVB2aaRQMQi8C8Qmr4ZtKq1bipxPljjeGyVKLZ7ecBoSumHb6PbfNpDV9snDaJCB59AIIHpozhbVToe5xLaYOvA13XT0Wl9SaRqO/7K+1LbULBveIxdIth/DV3O4Y1DGWwUdEXps+K6qsmJ/R12VnC+c1SCuqrBAENArLhtM6vC3jdtvvurfWrboQRFFsVzOnjUYj9Ho9DAYDdDofS5sFmbs31dh+XXDvyN4QIeIf/zniGEEVrVYi795UdP3PApe+PbFvFqpueA0f/FKFcd2APjE1ECrO1M9FPLFTWk3GXnuc+jnww9vArStRao1s1IFunwxrn/vz7V/H4LLOPkaRElG7ZzBb8fCaPW5rfRMGJiJnwuV46vOfXWp2o1O7YGbd75OND42CwWzFzcu+a3R9tFqJWdf1xsQrk2CtEfFS3n5pF4sGxvWPD2g/ob95wBphkHhqT5dWggeu6dURg3vEYsGE/jhRVolIlRJfHC6HruuTGJ220NG3F9s5AVXKGOws3oc3tlzArj93QvQHkz2/sNMcwITOvbHk9kE4fqESFystjlFa9hAMh/k7RBQevDV99kuMwcIvfva4BuniWwZ6/F3ivEfipaoaxESqcHVyLP5bt1u9s1D1EzIIg8Tbm+q/R8vw/I3JsFWWoUv1QURrtNh81IK/f1va4I1RJjVddq8fEaaynfH+wg3mAMbHREIUgb8XHHY7tDrUq74TUXjwttO9t7VCCw6ew9OTBjh+lzh3yTRlNwu7UCy3xiAMEk9vKvtu9Un5D0JRt5xbMoBpfTIxYurzuO2DYpc3hn3OoqO/0hznca6jpzmAgR5aTUTtj7ffB2ZrjddrK6trHM/hPKexKbtZ2IViuTUuuh0knt5Ujt3qG6xpGnH4a/TZ+RQeuS7Bccxt02XdXEekulmAO32O1F/oZg6gXqtCz7hoDOoei55x0QxBInJhX6HGnY7R3rtQnMPL/sG74NExmDAoye2yaoBUkxxyWUeXY6HqrmEQBomnN9X4nmpp81w3Ig7nIaNnXfOCt6bLuhVqMPe/wN3rpAEy3YdyDiARNVvD3SjsxvWPR7dYrceQdBde9g/eZmvjdUidOW/uG8ruGjaNBomnZY862C55uQpI0lSj4NExvpsutbHSV3S833MAPc0PIiICvHejeFvGzdtuPN70juuAr+aOCHl3DYMwiNy9qXRCqddrVNEd0bNTE6Yz2APRB2/zg0K1vh8RhR9P86ebM9bAeUHwhsb1j0dcjAZ6rfdVs1oDm0aDrGHfnDKqbmFvd4K0vqevpZG8jRYjIrJr6lgDb82t4TRqnTXC1mYf7OK8HRMghWCQ+vb8WRopXN6QRNS+tIVR6wzCULAPdqksk+b7afRBXd/TV40vFPN2iEg+gr0geEsxCEPFz769QPD1BgzFvB0ionDBPkIZ8DY/iMusEZHcMQhloK10WBMRhQKbRmWiLXRYExGFAoNQRsK9w5qIKBTYNEpERLLGICQiIllj02gwmMsBcxlQZZC2Q9IGb44gERG1DIMw0IxngC8fBJy3WeqbKa0mo0sKXbmIiMgtNo0Gkrm8cQgC0ia66x+SHiciorDCIAwkc1njELQ7mCctqUZERGGFQRhIVQbvj1f7eJyIiFodgzCQIvUte5yIiFodgzCQtK2/1yAREbUMgzCQ7HsNpma5Hg/iXoNERNQynD4RaK281yAREbUMgzAYWnGvQSIiahk2jRIRkayxRkhE1I4ZzFaUV1pgMFuh16oQG8Xt1xpiEBIRtVOlxios+L+f8M2Bc45j4/p1Qe7tVyJBFxnCkoWXoDWNlpWV4a677oJOp0NsbCxmzpyJS5cueT1/3rx56NevH7RaLXr06IEHH3wQBgMnoRMRNZXBbG0UggCw9cA55KzbB4PZGqKShZ+gBeFdd92FX375BZs3b8b69euxbds2zJ492+P5p0+fxunTp/HKK6/g559/xnvvvYdNmzZh5syZwSoiEVG7VV5paRSCdlv3n8VFk6WVSxS+BFEUxUA/6W+//YYBAwbgv//9L4YOHQoA2LRpEyZMmICTJ0+ia9eufj3PJ598grvvvhsmkwkREf614hqNRuj1ehgMBuh0umbfAxFRW/bTyXLcvOw7j49/NXcEBnWPbb0ChYC/eRCUGuGOHTsQGxvrCEEAyMjIgEKhwA8//OD389gL728IEhGRxNeAGB0HzDgEJQhLSkoQHx/vciwiIgKdOnVCSUmJX89x/vx5PPfcc16bUwGguroaRqPR5YuISO5io9QY16+L28fG9Y9HbJS6lUsUvpoUhAsWLIAgCF6/9u/f3+JCGY1GTJw4EQMGDMAzzzzj9dzc3Fzo9XrHV3Jycotfn4iordNrVci9/UqM6+9aKRnXPx65tw3iFAonTeojPHfuHC5cuOD1nN69e2P16tV45JFHcPHiRcfxmpoaREZG4pNPPsGtt97q8fqKigpkZWUhKioK69evR2Sk9yG+1dXVqK6udnxvNBqRnJzMPkIiIkijRy+aLKiosiImUoWO0fKZR+hvH2GTOt+6dOmCLl3cV7Wdpaeno7y8HIWFhRgyZAgAYOvWrbDZbEhLS/Na6KysLGg0Gnz55Zc+QxAANBoNNBqN/zdBRCQjeq1KNsHXXEHpI7z88stxww03YNasWdi1axe+++47zJ07F3fccYdjxOipU6fQv39/7Nq1C4AUgpmZmTCZTHjnnXdgNBpRUlKCkpIS1NbWBqOYREREwVtZ5sMPP8TcuXNx/fXXQ6FQ4Pbbb8ebb77peNxqteLAgQOorKwEAOzevdsxojQlJcXluY4cOYKePXsGq6hERCRjQZlHGEqcR0hERECI5xESERG1FQxCIiKSNQYhERHJGoOQiIhkjUFIRESyxiAkIiJZYxASEZGstbv9jezTIrkLBRGRvNlzwNd0+XYXhBUVFQDAXSiIiAiAlAt6vd7j4+1uZRmbzYbTp08jJiYGgiCEujgtZt9N48SJE+16pRzeZ/sih/uUwz0Cbfs+RVFERUUFunbtCoXCc09gu6sRKhQKdO/ePdTFCDidTtfm3oTNwftsX+Rwn3K4R6Dt3qe3mqAdB8sQEZGsMQiJiEjWGIRhTqPRYNGiRe1+82HeZ/sih/uUwz0C8rjPdjdYhoiIqClYIyQiIlljEBIRkawxCImISNYYhEREJGsMwjCwfPly9OzZE5GRkUhLS8OuXbu8nv/JJ5+gf//+iIyMxKBBg7Bx48ZWKmnLNOU+V61ahVGjRqFjx47o2LEjMjIyfP67hIum/jzt1qxZA0EQMHny5OAWMACaeo/l5eWYM2cOkpKSoNFokJqa2ibet029z6VLl6Jfv37QarVITk7G/PnzUVVV1UqlbZ5t27Zh0qRJ6Nq1KwRBwOeff+7zmoKCAvzud7+DRqNBSkoK3nvvvaCXM6hECqk1a9aIarVafPfdd8VffvlFnDVrlhgbGyuWlpa6Pf+7774TlUql+NJLL4m//vqr+NRTT4kqlUrct29fK5e8aZp6n3feeae4fPlycc+ePeJvv/0m/vnPfxb1er148uTJVi550zT1Pu2OHDkiduvWTRw1apR4yy23tE5hm6mp91hdXS0OHTpUnDBhgrh9+3bxyJEjYkFBgbh3795WLnnTNPU+P/zwQ1Gj0YgffviheOTIETEvL09MSkoS58+f38olb5qNGzeKTz75pLhu3ToRgPjZZ595Pb+4uFiMiooSs7OzxV9//VV86623RKVSKW7atKl1ChwEDMIQGzZsmDhnzhzH97W1tWLXrl3F3Nxct+f/8Y9/FCdOnOhyLC0tTfzLX/4S1HK2VFPvs6GamhoxJiZG/Oc//xmsIgZEc+6zpqZGHD58uPiPf/xDnD59etgHYVPv8e233xZ79+4tWiyW1ipiQDT1PufMmSOOGzfO5Vh2drY4YsSIoJYzkPwJwscee0y84oorXI5NmTJFzMrKCmLJgotNoyFksVhQWFiIjIwMxzGFQoGMjAzs2LHD7TU7duxwOR8AsrKyPJ4fDppznw1VVlbCarWiU6dOwSpmizX3Pp999lnEx8dj5syZrVHMFmnOPX755ZdIT0/HnDlzkJCQgIEDB+KFF15AbW1taxW7yZpzn8OHD0dhYaGj+bS4uBgbN27EhAkTWqXMraUt/g7ypd0tut2WnD9/HrW1tUhISHA5npCQgP3797u9pqSkxO35JSUlQStnSzXnPht6/PHH0bVr10b/AcNJc+5z+/bteOedd7B3795WKGHLNecei4uLsXXrVtx1113YuHEjioqK8MADD8BqtWLRokWtUewma8593nnnnTh//jxGjhwJURRRU1OD++67D0888URrFLnVePodZDQaYTabodVqQ1Sy5mONkMLekiVLsGbNGnz22WeIjIwMdXECpqKiAlOnTsWqVasQFxcX6uIEjc1mQ3x8PFauXIkhQ4ZgypQpePLJJ7FixYpQFy2gCgoK8MILL+Dvf/87du/ejXXr1mHDhg147rnnQl008oE1whCKi4uDUqlEaWmpy/HS0lIkJia6vSYxMbFJ54eD5tyn3SuvvIIlS5Zgy5YtuPLKK4NZzBZr6n0ePnwYR48exaRJkxzHbDYbACAiIgIHDhxAnz59glvoJmrOzzIpKQkqlQpKpdJx7PLLL0dJSQksFgvUanVQy9wczbnPp59+GlOnTsW9994LABg0aBBMJhNmz56NJ5980ut+eG2Jp99BOp2uTdYGAdYIQ0qtVmPIkCHIz893HLPZbMjPz0d6errba9LT013OB4DNmzd7PD8cNOc+AeCll17Cc889h02bNmHo0KGtUdQWaep99u/fH/v27cPevXsdXzfffDPGjh2LvXv3Ijk5uTWL75fm/CxHjBiBoqIiR8gDwMGDB5GUlBSWIQg07z4rKysbhZ09/MV2tKRzW/wd5FOoR+vI3Zo1a0SNRiO+99574q+//irOnj1bjI2NFUtKSkRRFMWpU6eKCxYscJz/3XffiREREeIrr7wi/vbbb+KiRYvazPSJptznkiVLRLVaLX766afimTNnHF8VFRWhugW/NPU+G2oLo0abeo/Hjx8XY2JixLlz54oHDhwQ169fL8bHx4t/+9vfQnULfmnqfS5atEiMiYkRP/74Y7G4uFj8+uuvxT59+oh//OMfQ3ULfqmoqBD37Nkj7tmzRwQgvvbaa+KePXvEY8eOiaIoigsWLBCnTp3qON8+feKvf/2r+Ntvv4nLly/n9Alqubfeekvs0aOHqFarxWHDhok7d+50PDZ69Ghx+vTpLuf/61//ElNTU0W1Wi1eccUV4oYNG1q5xM3TlPu87LLLRACNvhYtWtT6BW+ipv48nbWFIBTFpt/j999/L6alpYkajUbs3bu3+Pzzz4s1NTWtXOqma8p9Wq1W8ZlnnhH79OkjRkZGisnJyeIDDzwgXrx4sfUL3gTffPON2/9r9nubPn26OHr06EbXDB48WFSr1WLv3r3F//3f/231cgcSt2EiIiJZYx8hERHJGoOQiIhkjUFIRESyxiAkIiJZYxASEZGsMQiJiEjWGIRERCRrDEIiIpI1BiEREckag5CIiGSNQUhERLLGICQiIln7/xu2ijJtmn3vAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rng = np.random.default_rng(0)\n",
"jitter_data = test_data.copy()\n",
"jitter_data['X1'] = jitter_data['X1'] + rng.normal(0, 0.05, size=(200,)) \n",
"jitter_data['X2'] = jitter_data['X2'] + rng.normal(0, 0.05, size=(200,)) \n",
"\n",
"fig, ax = plt.subplots(1, 1, figsize=(5, 5))\n",
"for t in [0, 1]:\n",
" this_data = jitter_data.query('T == @t')\n",
" ax.scatter(this_data['X1'], this_data['X2'], label=f'T = {t}', edgecolor='w', linewidths=0.5)\n",
"\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "11abab7e",
"metadata": {},
"source": [
"Notice how the only region where overlap holds is when $X_1 = X_2 = 0$."
]
},
{
"cell_type": "markdown",
"id": "6a5d5c53",
"metadata": {},
"source": [
"## Applying OverRule with default arguments \n",
"[Back to Table of Contents](#toc)"
]
},
{
"cell_type": "markdown",
"id": "aa1618d9",
"metadata": {},
"source": [
"This method is currently implemented as a refutation method, though it can also be accessed through the functional API. We illustrate both methods below, before discussing their arguments.\n",
"\n",
"Note that OverRule returns rules in Disjunctive Normal Form (DNF), which can be read as an OR of ANDs, and are given when `print(refute)` is called."
]
},
{
"cell_type": "markdown",
"id": "2e165ecb",
"metadata": {},
"source": [
"### Example usage using CausalModel"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "cb218c53",
"metadata": {},
"outputs": [],
"source": [
"model = CausalModel(\n",
" data=test_data,\n",
" treatment=\"T\",\n",
" outcome=\"Y\",\n",
" common_causes=['X1', 'X2']\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "d1b9d144",
"metadata": {},
"outputs": [],
"source": [
"identified_estimand = model.identify_effect(proceed_when_unidentifiable=True)\n",
"\n",
"estimate = model.estimate_effect(identified_estimand, method_name=\"backdoor.linear_regression\")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "8b0e895d",
"metadata": {},
"outputs": [],
"source": [
"# We disucss configuration in more depth below.\n",
"support_config = SupportConfig(seed=0)\n",
"overlap_config = OverlapConfig()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "e6afb904",
"metadata": {},
"outputs": [],
"source": [
"refute = model.refute_estimate(\n",
" identified_estimand, \n",
" estimate, \n",
" method_name='assess_overlap', \n",
" support_config=support_config, \n",
" overlap_config=overlap_config,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "e11e8d3b",
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Rules cover 50.0% of all samples\n",
"Overall, 50.0% of samples meet the criteria for inclusion in the overlap set, \n",
"defined as (a) being covered by support rules and having propensity\n",
"score in (0.10, 0.90)\n",
"Rules capture 100.0% of samples which meet these criteria\n",
"\n",
"How to read rules: The following rules are given in Disjuntive Normal Form, \n",
"a series of AND clauses (e.g., X and Y and Z) joined by ORs. Hence, if a sample \n",
" satifies any of the clauses for the support rules, it is included in the support, \n",
"and likewise for the overlap rules.\n",
"\n",
"DETAILED RULES:\n",
"SUPPORT Rules: Found 2 rule(s), covering 100.0% of samples\n",
"\t Rule #0: (not X1)\n",
"\t\t [Covers 75.0% of samples]\n",
"\t OR Rule #1: (X1)\n",
"\t\t AND (not X2)\n",
"\t\t [Covers 25.0% of samples]\n",
"OVERLAP Rules: Found 1 rule(s), covering 50.0% of samples\n",
"\t Rule #0: (not X1)\n",
"\t\t AND (not X2)\n",
"\t\t [Covers 50.0% of samples]\n",
"\n"
]
}
],
"source": [
"print(refute)"
]
},
{
"cell_type": "markdown",
"id": "3dee0fdd",
"metadata": {},
"source": [
"### Example usage using the Functional API\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "91ff2447",
"metadata": {},
"source": [
"Note that the outcome $Y$ is not used as part of OverRule, so it is not required here."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "3ac5e83e",
"metadata": {},
"outputs": [],
"source": [
"refute = assess_support_and_overlap_overrule(\n",
" data=test_data, \n",
" backdoor_vars=['X1', 'X2'], \n",
" treatment_name='T', \n",
" support_config=support_config, \n",
" overlap_config=overlap_config\n",
")"
]
},
{
"cell_type": "markdown",
"id": "bbf343e9",
"metadata": {},
"source": [
"### Only fitting support or overlap rules\n",
"\n",
"You can also run OverRule to **only** learn either support or overlap rules, by passing in either `support_only=True` or `overlap_only=True`."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "176bc2a1",
"metadata": {},
"outputs": [],
"source": [
"refute = model.refute_estimate(\n",
" identified_estimand, \n",
" estimate, \n",
" method_name='assess_overlap', \n",
" support_only=True,\n",
" support_config=support_config,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "77d12a88",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Rules cover 100.0% of all samples\n",
"\n",
"How to read rules: The following rules are given in Disjuntive Normal Form, \n",
"a series of AND clauses (e.g., X and Y and Z) joined by ORs. Hence, if a sample \n",
" satifies any of the clauses for the support rules, it is included in the support, \n",
"and likewise for the overlap rules.\n",
"\n",
"DETAILED RULES:\n",
"SUPPORT Rules: Found 2 rule(s), covering 100.0% of samples\n",
"\t Rule #0: (not X1)\n",
"\t\t [Covers 75.0% of samples]\n",
"\t OR Rule #1: (X1)\n",
"\t\t AND (not X2)\n",
"\t\t [Covers 25.0% of samples]\n",
"No Overlap Rules Fitted (support_only=True).\n"
]
}
],
"source": [
"print(refute)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "add58e8f",
"metadata": {},
"outputs": [],
"source": [
"refute = model.refute_estimate(\n",
" identified_estimand, \n",
" estimate,\n",
" method_name='assess_overlap', \n",
" overlap_only=True, \n",
" overlap_config=overlap_config\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "f43178ac",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Rules cover 50.0% of all samples\n",
"Overall, 50.0% of samples meet the criteria for inclusion in the overlap set, \n",
"defined as (a) being covered by support rules and having propensity\n",
"score in (0.10, 0.90)\n",
"Rules capture 100.0% of samples which meet these criteria\n",
"\n",
"How to read rules: The following rules are given in Disjuntive Normal Form, \n",
"a series of AND clauses (e.g., X and Y and Z) joined by ORs. Hence, if a sample \n",
" satifies any of the clauses for the support rules, it is included in the support, \n",
"and likewise for the overlap rules.\n",
"\n",
"DETAILED RULES:\n",
"No Support Rules Fitted (overlap_only=True).\n",
"OVERLAP Rules: Found 1 rule(s), covering 50.0% of samples\n",
"\t Rule #0: (not X1)\n",
"\t\t AND (not X2)\n",
"\t\t [Covers 50.0% of samples]\n",
"\n"
]
}
],
"source": [
"print(refute)"
]
},
{
"cell_type": "markdown",
"id": "ac98cce5",
"metadata": {},
"source": [
"## Interpreting the output of OverRule \n",
"[Back to Table of Contents](#toc)\n",
"\n",
"We will use the following set of rules to illustrate how to interpret the output.\n",
"\n",
"```\n",
"SUPPORT Rules: Found 2 rule(s), covering 100.0% of samples\n",
"\t Rule #0: (not X1)\n",
" [Covers 75.0% of samples]\n",
"\t OR Rule #1: (X1)\n",
" AND (not X2)\n",
" [Covers 25.0% of samples]\n",
"OVERLAP Rules: Found 1 rule(s), covering 50.0% of samples\n",
"\t Rule #0: (not X2)\n",
" AND (not X1)\n",
" [Covers 50.0% of samples]\n",
"```\n",
"\n",
"For any sample, we can check if the sample is covered by the support by checking each rule (Rule 0 and Rule 1). The sample is in the support if ANY of these rules apply. Consider a point (X1 = 1, X2 = 1), which we know is not in the support of our original data. Since each binary input is treated as a boolean, we can check:\n",
"* Rule 0: (not X1) --> This rule is not satisfied by (X1 = 1, X2 = 1), because X1 = 1 is equivalent to X2 = True\n",
"* Rule 1: (X1) AND (not X2) --> This rule is not satisfied, because X2 = 1\n",
"\n",
"Similarly, we can check if a sample is covered by the overlap rules. Consider the point (X1 = 1, X2 = 0). This satisfies the support rules because it matches Rule 0. For the overlap rules, however, we can check:\n",
"* Rule 0: (not X2) AND (not X1) --> This rule is not satisfied, because X1 = 1.\n",
"\n",
"**Note**: When we later filter samples based on the rules learned by OverRule, we will filter out any samples that do not match BOTH a support rule, and an overlap rule."
]
},
{
"cell_type": "markdown",
"id": "cdc74f44",
"metadata": {},
"source": [
"## Configuring OverRule \n",
"[Back to Table of Contents](#toc)\n",
"\n",
"OverRule can be configured by passing in the following arguments\n",
"* `overlap_eps`: This should be a number $\\epsilon$ in (0, 0.5) which defines the overlap region as those points where $P(T = 1 \\mid X)$ is in $(\\epsilon, 1 - \\epsilon)$.\n",
"\n",
"* `cat_feats`: A list of any categorical features (automatically inferred if dtype is object).\n",
"\n",
"* `support_config`: A `SupportConfig` (see below), which controls optimization parameters\n",
"\n",
"* `overlap_config`: An `OverlapConfig` (see below), which controls optimization parameters"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "cc8734c5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SupportConfig(n_ref_multiplier=1.0, seed=None, alpha=0.98, lambda0=0.01, lambda1=0.001, K=20, D=20, B=10, iterMax=10, num_thresh=9, thresh_override=None, solver='ECOS', rounding='greedy_sweep')\n"
]
}
],
"source": [
"support_config = SupportConfig()\n",
"print(support_config)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "68fda6d1",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OverlapConfig(alpha=0.95, lambda0=0.001, lambda1=0.001, K=20, D=20, B=10, iterMax=10, num_thresh=9, thresh_override=None, solver='ECOS', rounding='greedy_sweep')\n"
]
}
],
"source": [
"overlap_config = OverlapConfig()\n",
"print(overlap_config)"
]
},
{
"cell_type": "markdown",
"id": "374ec6b5",
"metadata": {},
"source": [
"**Key Parameters**: For `SupportConfig` and `OverlapConfig`, the defaults are generally sensible, but some may be worth changing depending on the problem, and are detailed below."
]
},
{
"cell_type": "markdown",
"id": "bb2a9a2c",
"metadata": {},
"source": [
"* `alpha`: \n",
" + For learning support rules, we require that at least an `alpha` fraction of the original data is included, while minimizing the volume of the resulting support set. \n",
" + For learning overlap rules, at least an `alpha` fraction of the \"overlap points\" must be included, while excluding as many non-overlap points as possible. Here, the \"overlap points\" are defined as those where $P(T = 1 \\mid X)$ is within `[overlap_eps, 1 - overlap_eps]`, estimated using an XGBClassifier.\n",
"* `lambda0`, `lambda1`: These are regularization terms which control the complexity of the final rules. \n",
" + `lambda0` is a penalty added for every new rule (which can contain multiple literals). For instance, `(not X1) AND (not X2)` is a single rule. \n",
" + `lambda1` is a penalty added for every literal. For instance, `(not X1) AND (not X2)` contains two literals.\n",
"* `num_thresh`: OverRule turns all features into binary rules before learning rulesets. For continuous variables, it does so (by default) by converting into deciles. This default behavior can be changed via changing the number of thresholds, or specifying a dictionary (in `thresh_override`) that maps features to a `np.ndarray` of thresholds."
]
},
{
"cell_type": "markdown",
"id": "8861fb8c",
"metadata": {},
"source": [
"# Illustration on Lalonde and PSID datasets \n",
"[Back to Table of Contents](#toc)"
]
},
{
"cell_type": "markdown",
"id": "ce6eb962",
"metadata": {},
"source": [
"In this section, we show how OverRule can be applied in the course of doing another causal analysis. Here, we demonstrate both\n",
"* Using OverRule prior to conducting analysis\n",
"* Filtering samples using the learned rules\n",
"\n",
"Incorporating OverRule into a causal inference workflow looks like\n",
"* First, applying OverRule to learn rules, before any causal effect estimation is performed\n",
"* Second, filtering samples according to the learned rules\n",
"* Third, estimating causal effects on the filtered sample"
]
},
{
"cell_type": "markdown",
"id": "d2dff2a8",
"metadata": {},
"source": [
"For this illustration, we use two datasets to construct a somewhat artificial causal inference task:\n",
"* The original dataset from Lalonde (1986), which contains experimental samples, and which is loaded via `dowhy.datasets.lalonde_dataset()`. Because this is an *experimental* sample, we will see that **overlap holds for all the samples in this sample.**\n",
"* An observational dataset consisting entirely of non-treated samples, the PSID dataset, which was analyzed in Smith & Todd (2005) in conjunction with the original experimental dataset. This is loaded from `dowhy.datasets.psid_dataset()`\n",
"\n",
"By merging these datasets and then assessing overlap, we find that **in the combined dataset, overlap only holds for ~30% of samples.** Intuitively, many of the samples from the PSID are excluded, due to the restrictive inclusion criteria on the original experiment."
]
},
{
"cell_type": "code",
"execution_count": 48,
"id": "97cda5fa",
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"\n",
"import dowhy.datasets\n",
"from dowhy.causal_refuters.assess_overlap import assess_support_and_overlap_overrule\n",
"from dowhy.causal_refuters.assess_overlap_overrule import OverlapConfig, SupportConfig\n",
"\n",
"# Experimental sample from Lalonde\n",
"experiment_data = dowhy.datasets.lalonde_dataset()\n",
"\n",
"# Observational sample\n",
"observational_controls = dowhy.datasets.psid_dataset()\n",
"\n",
"experiment_data[\"exp_samp\"] = True\n",
"observational_controls[\"exp_samp\"] = False\n",
"\n",
"data = pd.concat([experiment_data, observational_controls], axis=0, ignore_index=True)\n",
"\n",
"support_config = SupportConfig(seed=0, lambda0=0.01, lambda1=0.01)\n",
"overlap_config = OverlapConfig(lambda0=0.01, lambda1=0.01)"
]
},
{
"cell_type": "markdown",
"id": "51a00dcc",
"metadata": {},
"source": [
"Before running any causal analysis, we first check for support and overlap. Note that the default setting for `alpha` is 0.98, which means that \n",
"* The support rules will try to cover at least 98% (`alpha`) of the data with a minimum volume set\n",
"* The overlap rules will try to cover at least 98% (`alpha`) of the samples which have propensity scores in `[overlap_eps, 1-overlap_eps]`"
]
},
{
"cell_type": "markdown",
"id": "d67ae8a7",
"metadata": {},
"source": [
"## Checking for Support/Overlap in the Experimental Sample"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "867c6981",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"All samples in the support region satisfy the overlap condition.\n",
"No Overlap Rules will be learned.\n",
"Rules cover 98.9% of all samples\n",
"\n",
"How to read rules: The following rules are given in Disjuntive Normal Form, \n",
"a series of AND clauses (e.g., X and Y and Z) joined by ORs. Hence, if a sample \n",
" satifies any of the clauses for the support rules, it is included in the support, \n",
"and likewise for the overlap rules.\n",
"\n",
"DETAILED RULES:\n",
"SUPPORT Rules: Found 4 rule(s), covering 98.9% of samples\n",
"\t Rule #0: ([re74 <= 0.000])\n",
"\t\t [Covers 73.3% of samples]\n",
"\t OR Rule #1: (black)\n",
"\t\t AND (not hisp)\n",
"\t\t AND ([age <= 35.000])\n",
"\t\t AND ([educ <= 12.000])\n",
"\t\t [Covers 71.5% of samples]\n",
"\t OR Rule #2: (not married)\n",
"\t\t AND ([re74 <= 7838.192])\n",
"\t\t AND ([re75 <= 5146.040])\n",
"\t\t [Covers 73.7% of samples]\n",
"\t OR Rule #3: (not black)\n",
"\t\t AND ([age <= 30.000])\n",
"\t\t AND ([educ <= 12.000])\n",
"\t\t AND ([educ > 8.000])\n",
"\t\t [Covers 10.6% of samples]\n",
"No Overlap Rules Fitted (support_only=True).\n"
]
}
],
"source": [
"refute_experiment = assess_support_and_overlap_overrule(\n",
" data=experiment_data,\n",
" backdoor_vars=\"nodegr+black+hisp+age+educ+married+re74+re75\".split(\"+\"),\n",
" treatment_name=\"treat\",\n",
" support_config=support_config,\n",
" overlap_config=overlap_config,\n",
")\n",
"\n",
"# Observe how everyone is in the overlap set, so we do not learn any overlap rules\n",
"print(refute_experiment)"
]
},
{
"cell_type": "markdown",
"id": "5faf90a0",
"metadata": {},
"source": [
"**Note**: As demonstrated above, if OverRule finds that all samples satisfy the overlap condition, then it will not bother to learn separate Overlap rules."
]
},
{
"cell_type": "markdown",
"id": "9af5d792",
"metadata": {},
"source": [
"## Assessing Support/Overlap on the combined sample"
]
},
{
"cell_type": "markdown",
"id": "060f1bdd",
"metadata": {},
"source": [
"Now, analyzing the combined dataset, we find a substantial lack of overlap"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "a7aa4ab7",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Rules cover 34.6% of all samples\n",
"Overall, 31.2% of samples meet the criteria for inclusion in the overlap set, \n",
"defined as (a) being covered by support rules and having propensity\n",
"score in (0.10, 0.90)\n",
"Rules capture 96.8% of samples which meet these criteria\n",
"\n",
"How to read rules: The following rules are given in Disjuntive Normal Form, \n",
"a series of AND clauses (e.g., X and Y and Z) joined by ORs. Hence, if a sample \n",
" satifies any of the clauses for the support rules, it is included in the support, \n",
"and likewise for the overlap rules.\n",
"\n",
"DETAILED RULES:\n",
"SUPPORT Rules: Found 2 rule(s), covering 98.1% of samples\n",
"\t Rule #0: ([re74 <= 33307.536])\n",
"\t\t AND ([re75 <= 32691.290])\n",
"\t\t [Covers 87.6% of samples]\n",
"\t OR Rule #1: (not nodegr)\n",
"\t\t AND (not hisp)\n",
"\t\t AND ([educ > 11.000])\n",
"\t\t [Covers 60.6% of samples]\n",
"OVERLAP Rules: Found 5 rule(s), covering 35.3% of samples\n",
"\t Rule #0: ([re74 <= 1282.723])\n",
"\t\t [Covers 20.0% of samples]\n",
"\t OR Rule #1: ([re75 <= 716.129])\n",
"\t\t [Covers 20.1% of samples]\n",
"\t OR Rule #2: (black)\n",
"\t\t AND (not married)\n",
"\t\t [Covers 15.2% of samples]\n",
"\t OR Rule #3: (not nodegr)\n",
"\t\t AND ([re74 <= 7837.067])\n",
"\t\t [Covers 12.2% of samples]\n",
"\t OR Rule #4: ([age <= 26.000])\n",
"\t\t AND ([educ > 11.000])\n",
"\t\t AND ([re75 <= 11637.097])\n",
"\t\t [Covers 7.7% of samples]\n",
"\n"
]
}
],
"source": [
"var_list = \"nodegr+black+hisp+age+educ+married+re74+re75\".split(\"+\")\n",
"refute = assess_support_and_overlap_overrule(\n",
" data=data,\n",
" backdoor_vars=var_list,\n",
" treatment_name=\"treat\",\n",
" support_config=support_config,\n",
" overlap_config=overlap_config,\n",
")\n",
"print(refute)"
]
},
{
"cell_type": "markdown",
"id": "e02a1e30",
"metadata": {},
"source": [
"## Filtering with OverRule before performing a causal analysis"
]
},
{
"cell_type": "markdown",
"id": "6058d9df",
"metadata": {},
"source": [
"OverRule can take any dataframe with the same columns as the original dataframe, and filter down to units that are included in both the support and overlap region."
]
},
{
"cell_type": "code",
"execution_count": 56,
"id": "6821c80e",
"metadata": {},
"outputs": [],
"source": [
"bool_mask = refute.predict_overlap_support(data)\n",
"\n",
"# This is a simple wrapper around refute.predict_overlap_support\n",
"filtered_df = refute.filter_dataframe(data).copy()\n",
"\n",
"assert np.array_equal(data[bool_mask], filtered_df)"
]
},
{
"cell_type": "code",
"execution_count": 47,
"id": "acf86d4d",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAKsAAAAPCAYAAACbZT/hAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHZUlEQVRoBeWa25EVOQxAe6kJYGAjADLgkcGQwbJEAGTAFl/wR0EGQAQ8MgAi4JEBbATAZDB7jq9l7G7f22742dpVlce2JEuyJMvuC9PZ2dn0b2wPHjy4MrcL3DHt0hz/X5r/H/dtTGl3aMeHYnk0VfDw4cPHefqV/jLtMbgvFcveIXyXIN7NDMf0zl3/NuO6HfQrEO7Sx9rgewdOOZ8ywrFwddft/mae+xVOvlfgF3rBadNfmfca/Tfn4ENHJu068OGPwL+oebfIy7yxx0P+Gdp3GFTZeDBm8D2NNfQXaLfBnVa4ifnmGLLmoI+UD8+abnNAnqfw0i3gFPz5kqxMPsLyiP61rPQ69CP9DdrBhM28Bj2CMTH+g/Vv6G/SkkzlduAVuF6ymEiCG1G/MrSvcTA4D0TRy3xibrIe04pexgZCZ9yQR2Cso2OPJbkzr3a5p4Snd27zEE+ZZ1Sevhz1z+i+tWE1ZvCkOKJfW59k2/Xp38yv0lJsM9+ojRP8+nPNR0O6kXOdpp97eXYCPhWYlKwovgNiHlyz2WCb8SXAjHvg+jvwm5yRIBH8+9AC16yF916DaCefoN9sUe0MunoN2Bxug3hHq/WamPOkNjjK0OnnaQHOraKxB/E6vnbmFnlb/LO6b43Jdo/E7DnsF+BPiZrXquMD4zq2W2xUzIiPRnW7n0WOgfNASEtxPOcEMCl61e09+BOYDdQhcO1pbomPNc73AnRPd7NmL/N+glVusck97J7Qz529mJAGPRzjjaBtz2gFoHvD1LqG5GUBm/1TFO8fjMbM/dSHLCRqUx3bYRvxw6iPRnWbZz1obs1UWeHS8U1w8srYpPS6SmXyrsN4A15XpilvSAZPbw9uwdNcOz2mFZybvIecN/Q+N+KAWPXmerXxSsXDtIHjPLP6equErIapmozKm5D1M/6pVHWHqzHLel38rSPBN65wjfZ2o42rPkJe+HNE9yK3WG8MH2lgwFElNHC9/kIPuQ+HTB2ZrknGi0MAzut/nkwLcfB5Nbnp32lWPt+sVoAEjF/bmHiCvzP2bWO1rZ8jTNNVsu9JYRWdWBtyDd4X5uJv0QyqMpuPNuij8ljaAmvX/LO270iEVnA78+r30IntxU+fCulG2Q1//F2xcdVHv6hbmyws8TGcDDvib2zkUCUZcc6EcANsIK7TDP4HWgPwaIhOjKrd0KuJOl/Cl+zK67zGvY6tVAkYW1FNfAPsAVGu/SqwTnu1p3ZK7PUa9IJn7GHwC3pRBUIRtJ68IE+ZftA/MI/se0vMtFedc9BWIfabJhtslH/NR5t0JwN2f4zfIobnKoZDwziFh3gMhg/3JzSrzguaX9pWvRr8mWpRbWsGx/CYlClR89wkNEmbigyP8uWz+kk3+UzquV7QC/AjweocX8kROE/13MaX8D4HHzwLYSAaeXMG1q76B56hfc9ld+YRMz82J+SWhGVsooZvm6IB7aCN0GP/Iz7apDvbafxO0FMKknjBZO29KRKRP3GC430T+NUeZZ4qHZJ+RnIBOKtfk2ziN4COvYQcNxTyDK5vX69uP4Diej6YWPBqh2uCX5EBTQAz0l8dDJRX4AJW5PX4F/5ZMP1ANPsGPRwz7DIGF2neQL7vrVj6Lz5qenuFnPx7yMbeusZHP6nb93BP9nSkQBr09joQAcQp6i7esaRNpSsFOfHuC9IHBp7ok0w7pj8oy4Xw+MHkm+uq8w6EXTreQBRgjZXSKvuZpm4d3gB0D43y66/7iXn4wgDvg3RQauI+ecEDfdU/8Gr30L4rO8MPoco+cMXP8oM3CQqAi2s28TEftXGTj0Z0F6N2A2/EYntNO5cncX3WNMdRWRclecboifLKD0fNyGlqkK/DY6UtDZxOMpnFhQOtXj1ZyR74vKqkm/wGogFwbtYkDfsLHZrOuExfKirjUq2hudee7pDROHJAnutG/CPf6r5lyvCrMdPv/goQ/hu1UfWbfJTtrbu57kTDFv1unoRN9ZrpKM98a0Wi1AxWNhOju7hilF5vPEg6XwjaIumR/T3TS/Iwfwa+fNwoIIOVMsnQptxMtCaBMq8bb/TBp5M8MHPZJnC8UX0e6I856At1FpmD8pSz6p+sbHXfmc9uKGbY6N78cf4i4xRHen2jL+uba9RGlqWn3KqPNuhWphD50n3mpMqKUAP1LQtPq/KG/mSSHskixdHOaJ7CGgy+V1gBeHSSTvGDKjmpENuBPLYa/KfB5m3LPP61q05qx+VNHALg9Zr356uSxIw9sTrYPST50YMrNoKzInu4yuFlrH1zXwzJY50w6p/RfU/YNBQzdGvnPPj6wT3Xz7ZRG9U95KMNumFNEHlwGoi6/83/5SJggIwGSMavtOu0R+DrDcnnW9Bgzt9AntQ6kXSS/wJRKhHzAuBNRnlcJ+iA9+Djy1xaVECvcx3uh1SzEeZWy/uZTpfABG70MtduZfbA26OuMhNzfaFPBPU3vvgJeUP+Qe7QvjUK3tGYxcGL/Sz8k+UN2SivgP6DPqp4HB7UnXndu4Ww+xPhPz4Vghy51fgMAAAAAElFTkSuQmCC\n",
"text/latex": [
"$\\displaystyle 0.345826235093697$"
],
"text/plain": [
"0.34582623509369675"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bool_mask.mean()"
]
},
{
"cell_type": "markdown",
"id": "2364ad4c",
"metadata": {},
"source": [
"## Performing causal inference after filtering\n",
"\n",
"Suppose that we were only given the combined sample, and wanted to perform causal effect estimation. Having used OverRule to characterize the Support/Overlap region, we can now filter the data and perform analysis on the reduced dataset."
]
},
{
"cell_type": "code",
"execution_count": 43,
"id": "1500a720",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Causal Estimate is -958.9002273474416\n"
]
}
],
"source": [
"# Note that we are using the filtered_df here\n",
"model = CausalModel(\n",
" data = filtered_df,\n",
" treatment='treat',\n",
" outcome='re78',\n",
" common_causes=\"nodegr+black+hisp+age+educ+married+re74+re75\".split(\"+\")\n",
")\n",
"identified_estimand = model.identify_effect(proceed_when_unidentifiable=True)\n",
"\n",
"estimate = model.estimate_effect(identified_estimand,\n",
" method_name=\"backdoor.propensity_score_weighting\",\n",
" target_units=\"ate\", \n",
" method_params={\"weighting_scheme\":\"ips_weight\"})\n",
"\n",
"print(\"Causal Estimate is \" + str(estimate.value))"
]
},
{
"cell_type": "markdown",
"id": "a7376c23",
"metadata": {},
"source": [
"**Disclaimer**: This is just an illustration of how to incorporate OverRule into a causal inference workflow. However, the task itself is somewhat artificial: Here, we have combined experimental data (from Lalonde) with observational data (from PSID), which can introduce bias into our estimates. We only do so here to demonstrate that OverRule can pick up on a lack of overlap between treatment and control in the combined sample."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.9.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}