{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# NHEFS Dataset\n", "**NHANS (National Health and Nutrition Examionation Survey) Epidemiologic Followup Study**\n", "\n", "The effect of smoking on health measurements was a decades-long debate. \n", "As Judea Pearl puts it, it was \"probably the highest-profile medical question of the century\".\n", "The debate has divided not only the public across the US, but also physicians, epidemiologists and statisticians alike. \n", "Many known academic figures had took stand in this debate, researchers like \n", "[Abraham Lilienfeld](https://en.wikipedia.org/wiki/Abraham_Lilienfeld), \n", "[Austin Bradford Hill](https://en.wikipedia.org/wiki/Austin_Bradford_Hill), \n", "[Richard Doll](https://en.wikipedia.org/wiki/Richard_Doll), and even \n", "[Ronald Fisher](https://en.wikipedia.org/wiki/Ronald_Fisher).\n", "\n", "In fact, since a randomized control trial on the issue was always regraded unethical, discovering the detrimintal effects of smoking is one of the most known examples of deriving a causal effect from observational studies.\n", "If you believe smoking is not good for your health, you must also believe in causal inference.\n", "\n", "In this example we will analyze the effect of smoking cessation on weight gain from real-world data.\n", "People have noticed that that smoking cessation was associated with positive weight gaining.\n", "Formaly, the hypothesis made was that people who quit smoking gain more weight than they would have if they would've kept smoking.\n", "Actually, this association was used by pro-smoking advocators as the adverse effects of higher BMI (due to weught gain) have masked the positive effects that quitters benefitted from.\n", "\n", "To quantify the effect of smoking cessation on weight gain, we will use data from the NHEFS. The NHANS (National Health and Nutrition Examionation Survey) Epidemiologic Followup Study (NHEFS) is a national longitudinal study that was jointly initiated by the American National Center for Health Statistics and the American National Institute on Aging. \n", "The study was designed to follow up on smokers, whom some of which quit, and record their weight (among other measurements) for a period of 11 years: from 1971 to 1982. \n", "We will follow an analysis suggested by Hernán and Robins in their [Causal Inference book](https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/) to answer the question how much smoking cessation contributes to weight gain in the population." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Data\n", "First, let's download the dataset from the [Causal Inference book's](https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/) webpage." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1629, 64)\n" ] }, { "data": { "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", " \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", "
seqnqsmkdeathyrdthmodthdadthsbpdbpsexage...birthcontrolpregnanciescholesterolhightax82price71price82tax71tax82price71_82tax71_82
023300NaNNaNNaN175.096.0042...2NaN197.00.02.1835941.7399901.1022950.4619750.4437870.640381
123500NaNNaNNaN123.080.0036...2NaN301.00.02.3466801.7973631.3649900.5718990.5493160.792969
224400NaNNaNNaN115.075.0156...02.0157.00.01.5695801.5134280.5512700.2309880.0561980.320251
32450185.02.014.0148.078.0068...2NaN174.00.01.5065921.4519040.5249020.2199710.0547940.304993
425200NaNNaNNaN118.077.0040...2NaN216.00.02.3466801.7973631.3649900.5718990.5493160.792969
\n", "

5 rows × 64 columns

\n", "
" ], "text/plain": [ " seqn qsmk death yrdth modth dadth sbp dbp sex age ... \\\n", "0 233 0 0 NaN NaN NaN 175.0 96.0 0 42 ... \n", "1 235 0 0 NaN NaN NaN 123.0 80.0 0 36 ... \n", "2 244 0 0 NaN NaN NaN 115.0 75.0 1 56 ... \n", "3 245 0 1 85.0 2.0 14.0 148.0 78.0 0 68 ... \n", "4 252 0 0 NaN NaN NaN 118.0 77.0 0 40 ... \n", "\n", " birthcontrol pregnancies cholesterol hightax82 price71 price82 \\\n", "0 2 NaN 197.0 0.0 2.183594 1.739990 \n", "1 2 NaN 301.0 0.0 2.346680 1.797363 \n", "2 0 2.0 157.0 0.0 1.569580 1.513428 \n", "3 2 NaN 174.0 0.0 1.506592 1.451904 \n", "4 2 NaN 216.0 0.0 2.346680 1.797363 \n", "\n", " tax71 tax82 price71_82 tax71_82 \n", "0 1.102295 0.461975 0.443787 0.640381 \n", "1 1.364990 0.571899 0.549316 0.792969 \n", "2 0.551270 0.230988 0.056198 0.320251 \n", "3 0.524902 0.219971 0.054794 0.304993 \n", "4 1.364990 0.571899 0.549316 0.792969 \n", "\n", "[5 rows x 64 columns]" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "nhefs = pd.read_csv(\"https://cdn1.sph.harvard.edu/wp-content/uploads/sites/1268/1268/20/nhefs.csv\")\n", "print(nhefs.shape)\n", "nhefs.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see the dataset has 1629 inidividuals in it. \n", "\n", "### Restriction criteria\n", "Since the study focused on smokers, some also elderly, not all of the participants have survived the entire 11 years of follow-up period.\n", "Although not the best practice in causal analysis, we will, for the sake of simplicity, discard from the dataset individuals whom we don't have their weight measured at the end period." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1566, 64)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nhefs = nhefs.dropna(subset=[\"wt82\"]) # weight meausred at 1982 - the end of the followup period\n", "nhefs.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Variables selection\n", "Second, for our analysis to be accurate, we must identify the variables that act as confounders - meaning, the variables that affect both the outcome and the treatment assignment. \n", "For example, *age* of participants can affect both the outcome (weight), as people, on average, increase their weight as they grow older (this has several physiological reasons, such as change in metabolism). However, age can also effect the treatment (smoking cessation) as we can hypothesis that older people will be less likely to change their habits. \n", "Another example can be the *education level* of participants. More educated people might more prone to quit smoking (as they may be more exposed to the debate and understand the risks depicted in the arguments), but they also might be more conscious to the benefits of healthier food which in turn reduce the weight they gain over the years." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1566, 9)\n" ] }, { "data": { "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", "
activeageeducationexerciseracesexsmokeintensitysmokeyrswt71
00421210302979.04
10362000202458.63
20562211202656.81
3168121035359.42
41402100201987.09
\n", "
" ], "text/plain": [ " active age education exercise race sex smokeintensity smokeyrs \\\n", "0 0 42 1 2 1 0 30 29 \n", "1 0 36 2 0 0 0 20 24 \n", "2 0 56 2 2 1 1 20 26 \n", "3 1 68 1 2 1 0 3 53 \n", "4 1 40 2 1 0 0 20 19 \n", "\n", " wt71 \n", "0 79.04 \n", "1 58.63 \n", "2 56.81 \n", "3 59.42 \n", "4 87.09 " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "confounders = [\"active\", # Measure of daily activity on 1971.\n", " \"age\", # Age in 1971.\n", " \"education\", # Amount of education in 1971: from 8th-grade to college level education.\n", " \"exercise\", # Measure of recreational excercise.\n", " \"race\", # Caucasian or not.\n", " \"sex\", # Female or male.\n", " \"smokeintensity\", # Number of Cigarettes smoked per day in 1971. \n", " \"smokeyrs\", # Years of smoking.\n", " \"wt71\"] # Weight in Kilograms in 1971.\n", "X = nhefs[confounders]\n", "print(X.shape)\n", "X.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We continue to define our input matrix for the regression analysis: \n", "We square our continuous variables and encode the categorical ones using a one-hot encoding scheme." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "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", " \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", "
ageracesexsmokeintensitysmokeyrswt71active_1active_2education_2education_3education_4education_5exercise_1exercise_2age^2wt71^2smokeintensity^2smokeyrs^2
04210302979.040000000117646247.3216900841
13600202458.630010000012963437.4769400576
25611202656.810010000131363227.3761400676
3681035359.421000000146243530.736492809
44000201987.091010001016007584.6681400361
\n", "
" ], "text/plain": [ " age race sex smokeintensity smokeyrs wt71 active_1 active_2 \\\n", "0 42 1 0 30 29 79.04 0 0 \n", "1 36 0 0 20 24 58.63 0 0 \n", "2 56 1 1 20 26 56.81 0 0 \n", "3 68 1 0 3 53 59.42 1 0 \n", "4 40 0 0 20 19 87.09 1 0 \n", "\n", " education_2 education_3 education_4 education_5 exercise_1 exercise_2 \\\n", "0 0 0 0 0 0 1 \n", "1 1 0 0 0 0 0 \n", "2 1 0 0 0 0 1 \n", "3 0 0 0 0 0 1 \n", "4 1 0 0 0 1 0 \n", "\n", " age^2 wt71^2 smokeintensity^2 smokeyrs^2 \n", "0 1764 6247.3216 900 841 \n", "1 1296 3437.4769 400 576 \n", "2 3136 3227.3761 400 676 \n", "3 4624 3530.7364 9 2809 \n", "4 1600 7584.6681 400 361 " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X = pd.get_dummies(X, columns=[\"active\", \"education\", \"exercise\"], drop_first=True)\n", "X = X.join(X[['age', 'wt71', 'smokeintensity', 'smokeyrs']]**2, rsuffix=\"^2\")\n", "X.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lastly, we extract our treatment variable `qsmk` (quit smoking), and our outcome variable `wt82_71` (the difference in participants' weight between 1982 and 1971)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "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", "
qsmkwt82_71
00-10.093960
102.604970
209.414486
304.990117
404.989251
\n", "
" ], "text/plain": [ " qsmk wt82_71\n", "0 0 -10.093960\n", "1 0 2.604970\n", "2 0 9.414486\n", "3 0 4.990117\n", "4 0 4.989251" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = nhefs[\"qsmk\"]\n", "y = nhefs[\"wt82_71\"]\n", "pd.concat([a, y], axis=\"columns\").head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Model\n", "After defining the confounders and further building our design matrix `X`, we can continue to define the causal model.\n", "\n", "In this case we will use an Inverse Treatment Probability Weighting (IPTW, or IPW) causal model. \n", "Briefly, this model will model the probability of a participants' to quit smoking and use it to emulate two equal-sized populations: one of quitters and another of persisters. \n", "In this synthetic population, we could use the actual change in weight to estimate what would have happen if everyone were to quit or everyone were to persist smoking.\n", "\n", "But first, we will need to use a machine learning model to estimate the propensity score $\\Pr[A=1|X]$ - the probability of each participant to quit smoking. \n", "Following the design matrix we prepared above, and given the binary nature of our treatment, we will choose a logistic regression for this task." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from sklearn.linear_model import LogisticRegression\n", "\n", "learner = LogisticRegression(penalty='none', # No regularization, new in scikit-learn 0.21.*\n", " solver='lbfgs',\n", " max_iter=500) # Increaed to achieve convergence with 'lbfgs' solver" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once we defined a learner, we can simply plug it into the causal model" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "from causallib.estimation import IPW\n", "\n", "ipw = IPW(learner)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimating Causal Effect\n", "Once we defined the causal model (yes, that all it took), we can move on to estimate the effect of smoking cessation on weight gain.\n", "\n", "First, we will fit our causal model. \n", "Second, we'll predict the potential outcomes: what would be the weight difference if everyone were to quit smoking or everyone were to persist smoking. \n", "Third, we will use the two potential outcomes to estimate the effect: the difference of the two potential outcomes." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "ipw.fit(X, a)\n", "outcomes = ipw.estimate_population_outcome(X, a, y)\n", "effect = ipw.estimate_effect(outcomes[1], outcomes[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Checking the potential outcomes, we can see that the average difference in weight (across 11 years) if everyone would have quit smoking (`1`) is 5.28kg (11.64lbs) while the averge weight difference if everyone would have persist smoking (`0`, not quit) is 1.77kg (3.9lbs)." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 1.765489\n", "1 5.283029\n", "dtype: float64" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "outcomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Therefore we can conclude that the average additive effect (`diff`) of smoking cessation on weight gain is 3.52kg (7.76lbs)." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "diff 3.51754\n", "dtype: float64" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "effect" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Indeed, we can conclude that smoking cessation accounts for a weight gain of 3.5kg on average." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Unadjusted estimation\n", "To compare, we can ask what conclusion would have we achive if we did not control for confounding. \n", "What would've been the result if we haven't adjusted for treatment assignment biases and treated the data as if it came from a randomized control trial?" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "diff 2.540581\n", "dtype: float64" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from causallib.estimation import MarginalOutcomeEstimator\n", "\n", "moe = MarginalOutcomeEstimator(None).fit(X, a, y)\n", "outcomes = moe.estimate_population_outcome(X, a, y)\n", "moe.estimate_effect(outcomes[1], outcomes[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that when failing to adjust for confounding factors, we underestimate the effect of smoking cessation on weight gain by 1kg (2.2lbs).\n", "\n", "To see why, we can examine the how balanced were the treatment groups before and after IP weighting. \n", "We can check that using a _Love plot_, a graphical way to present distribution distance between the treated and untreated. It plots the standardized-mean-difference between the treatment groups for each covariate marginally." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAFzCAYAAACzRzkmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABXy0lEQVR4nO3de5xVZb3H8c9XHGUQhQQLBQ2ow6iADgIaIGCCggkKiFhhHbTyWJnWSY5SeU2PKB6vWaSmWI0lKVCgIgIKhJqgDFcZvEAKigI6XOQ2DL/zx/Ns2HPfw+y5MPN7v17zYu+11n7Ws9aM/uZZa83zlZnhnHPOufIdUtsdcM455w4GXjCdc865FHjBdM4551LgBdM555xLgRdM55xzLgVeMJ1zzrkUHFrbHXDVp2XLlta2bdva7oZzzh1U3njjjY1mdkzx5V4w67G2bduycOHC2u6Gc84dVCT9u7TlfknWOeecS4EXTOeccy4FXjCdc865FHjBdM4551LgBdM555xLgRdM55xzLgVeMJ1zzrkUeMF0zjnnUuAF0znnnEuBF0xX46YsWkevsbNpd/2z9Bo7mymL1tV2l5xz9cXW9XD/qbD147Q37QXT1agpi9YxZtJS1uXvwIB1+TsYM2mpF03nXHrMuQvy34c5d6a9aZ9LthZJmgIcDzQG7jezhyV9D7gOyAcWA7vM7CpJxwDjgRPix39qZvNrvtdVM+6FPHYUFAIw8LCV+5a/NHUl+Yub07FjR7p3705BQQE5OTklPp+dnU12djbbt29n4sSJJdZ369aNTp06sXnzZiZPnlxifY8ePcjKymLjxo1MmzatxPo+ffrQvn171q9fz/Tp00us79evH8cffzwffPABs2bNKrF+4MCBtGrVivfee4+5c+eWWD9o0CBatmxJXl4er776aon1Q4cOpVmzZixbtqzUeYBHjBhBkyZNyM3NJTc3t8T6kSNHkpGRwYIFC1i+fHmJ9aNGjQLglVdeYdWqVUXWZWRkMHLkSADmzJnD6tWri6xv0qQJI0aMAGDmzJmsXbu2yPqjjjqKYcOGATB9+nTWr19fZH2LFi0YPHgwAFOnTmXTpk1F1rdq1YqBAwcCMGnSJLZs2VJkfZs2bejfvz8AEydOZPv27UXWt2vXjr59+wKQk5NDQUFBkfUdOnSgZ8+eAEyYMIHi/GevHvzsbV3P9DfXMND2Qm4O9L0OjvxSib4cKB9h1q7Lzawr0A24WlJr4Abga0Av4MSkbe8H7jWz7sBFwKOlNSjpCkkLJS3csGFD9fb+AHyYv6PU5bv3FNZwT5xz9c6cu8AsvLa9aR9lyhKNuxon6WZgaHzbFrgDOMnM/jOuvxroEEeYnwAfJn38GCDLzLaV1X63bt2srqWV9Bo7m3WlFM3WzTOZf/3ZtdAj51y9kLh3uWfn/mWHNoZrllR6lCnpDTPrVny5jzBriaSzgP5ADzM7FVgErCznI4cAXzOz7PjVurxiWVeNHpBFZkajIssyMxoxekBWLfXIOVcvzLkrjCqTpXmU6QWz9jQDPjOz7ZJOJFyGPQLoK+kLkg4lXHpNmAH8JPFGUnZNdjZdhnRpzR3DOtO6eSYijCzvGNaZIV1a13bXnHMHs7znoHB30WWFu8PyNPFLsrVE0uHAFMKl2DygOXAz0AEYDXxKGHGuNbNfSmoJPAScRHhYa66ZXVnePuriJVnnnKvryrok60/J1hIz2wWcV3y5pIXxadlDgcmEooqZbQQuqdFOOuec28cvydY9N0vKBZYBq4kF0znnXO3yEWYdY2bX1nYfnHPOleQjTOeccy4FXjCdc865FHjBdM4551LgBdM555xLgRdM55xzLgVeMJ1zzrkUeMGsQySdLOkTSdPjxAWJ5cdLeknSCknLJV1Tm/10zrmGyAtmHSHpOGAiIb1kOfBw0uo9wM/N7GTCnLM/lnRyzfeycqYsWkevsbNpd/2z9Bo720OinXPVL5FasvXjtDftBbOKJE2R9EYc+V0Rl31P0ipJr0t6RNJv4vJjJD0jaUH86hWXHwU8BVxhZvPN7OfABkm3ApjZR2b2Zny9FXgLqNOzlU9ZtI4xk5ayLn8HBqzL38GYSUu9aDrnqtecuyD//bRnYYJPvl5lko42s08lZQILgAHAfOA0YCswG1gcMy2fBH5rZv+UdALwgpmdVMn9tQXmAp3MbEt529bm5OvJuZcDD9ufWnb4oY3ockJzT7evD+n2wPTp01m/fn2R9S1atGDw4MEATJ06lU2bNhVZ36pVKwYOHAjApEmT2LKl6I9xmzZt6N+/PwATJ05k+/btRda3a9eOvn37ApCTk0NBQUGR9R06dKBnz54ATJgwgeL8Z68e/+wdlcngFVeHTMwDzMIEz8OsTldLWgy8BhwPfAeYY2afmlkB8LekbfsDv4lzxf4DOEpS01R3FLd9BvhpWcVS0hWSFkpauGHDhgM7ojT4sJSQaIDdewpruCfOuQbjoyX7MzHTnIUJPsKskhgCfRtwbsy1fBm4DxhqZv8Zt7ka6BBHmBuBNma2s/QWy91XBjCNMCq9J5XP1JURZrLWzTOZf/3ZtdAj51y9lrh3uSfpf68HOMr0EWb1qJEQaEkC/gC8lWqxrG2jB2SRmdGoyLLMjEaMHpBVSz1yztVrc+7aP7pMSPMo0wtm1UwHDpX0FjCWcFl2HfC/wOuEe5lrgM1x+6uBbpKWSFoBlBsAnaQX4VLv2ZJy49c30ncY6TekS2vuGNaZ1s0zEWFkecewzgzpUqefVXLOHazynoPC3UWXFe4Oy9PEL8lWA0lNzWxbUgj0Y2ZW8umAalabl2Sdc+5g5Zdka5aHQDvnXD3jAdLVwEOgnXOu/vERpnPOOZcCL5jOOedcCrxgOueccynwgumcc86lwAumc845lwIvmM4551wKvGA655xzKWgQBVPSWZJK5uyUvf2tkvqn0GbPqveuRLuPJsKhJf0i3e3XBA+Odq4GVGNQsitdgyiYlWVmN5rZzAo2OwtIe8E0s++b2Yr49qArmB4c7VwNqcagZFe6OjnTj6QjgIlAG6AR8GvgTuAvwHnAHuAK4A7gq8A4MxsfUz3uitsYcJuZPVWs7e7Aw8BwoDlwD9AU2AiMMrOPJE0AppnZ05LWAE8Ag4EM4GJgJ2Hi9EJJlxISSFYC44ET4q5+ambzJd0cl7WP/95nZg+Udoxm9lSMCLs29i8zTrG3HHgX+NTM7ovHcTvwiZndf6DnuTqMeyGPHQUh8zI5OPqlqSvJX+zB0fU6vNeDo2vuZ++lmbB2G6NsL+TmQN/rDigo2VVOXR1hDgQ+NLNTzawTIRUE4H0zywbmARMIReVrwC1x/TAgGziVENY8TtKxiUbjJdTxwIXA+8CDwHAz6wo8BtxeRn82mtlpwO+Aa81sTWznXjPLNrN5wP3xfXdCpNejSZ8/ERgAnA7cFLMtyzpGAMzsemBHbH9k7N9343EcAnwT+HPxjtZ2gLQHRztXA/I/2P+6GoKSXenqZFqJpA6E7MinCCO9eXGk18vM1km6HOhhZj+I278PnALcBCw1s8fi8j8BfwO2EPIkdxDCnj+U1Al4BXgv7rYR8JGZnVvKCDOx3zOA282sfxw5bjOzu+O+PgE+TDqMY4AswmixwMxuj9u9BZwDNCl+jHH9y4SivFDSNjNrmnReXgT+B/gS8H0zG17eeayNtBIPjnaumqUxKNmV7qBKKzGzVcBpwFLgNkk3xlW74r97k14n3ld0efkjwqXULvG9gOVxBJdtZp3N7NwyPpvYV2E5+zkE+FpSe63NbFuxz+9ro5xjLM+jwCjgMsKIs87x4GjnqlkNBCW70tXJginpOGC7mf0ZGEcoLKmYB1wiqZGkY4A+hCBngHzgfOAOSWcBecAxknrEfWZI6liJbm4Fjkx6P4NwLzNxDNnlfTjFYyyIl28TJhMu5XYHXqhEX2uMB0c7V81qICjZla5OPvQDdCbcf9wLFAA/BJ5O4XOTgR7AYsJDP/9jZuslnQhgZh9LGgQ8D1xOuAf6gKRmhHNxH+EBm1RMBZ6WdCGhUF4NPCRpSWxrLuHBoMocY3EPA0skvWlmI81st6SXgHwzq7M3BYd0ae0F0rnq8vOVFW/jqkWdvIfpShcf9nkTuNjM3q5o+9q4h+mccwe7g+oepispTmbwDjArlWLpnHMuverqJVlXTJzMoH1t98M55xoqH2E655xzKfCC6ZxzzqXAC6ZzzjmXAi+YzjnnXAq8YDrnnHMp8ILpnHPOpcALZhkqGzqdpn2eLOkTSdMlHZq0/HhJL0laIWm5pGtqsl/OOee8YNYoSY3KWXccIR9zKGF6voeTVu8Bfm5mJxPizH4cJzKocVMWraPX2Nm0u/5Zeo2d7cHQrn5KJIJs/bi2e+LqkHpTMCUdIelZSYslLZN0iaQ1ku6QlBszIk+T9IKkdyVdGT8nSePiZ5ZKuqSUtrtLWiTpK5K6Spoj6Y3Y1rFx+ZtJ2/9H4n3sw53x/cWSro4jxSWS/hq3OYoQ83WFmc03s58DGyTdCmBmH5nZm/H1VuAtoMYna52yaB1jJi1lXf4ODFiXv4Mxk5Z60XT1z5y7IP99TwBxRdSnmX4SgcznA8QJ1e8khk5LupcQOt0LaAwsI4RAJ4dOtwQWSJqbaDSGTj9ICJ3+CPgTcKGZbYjF9XYzu1zSZknZZpZLiN96PKlvm2IANZI+BNqZ2S5JzQHMbAvQO/lgzOy60g5SUltCRNm/DuAcVcm4F/LYURDmfB942P4JoF+aupL8xc3rV6L93Lkl1g8aNIiWLVuSl5fHq6++WmL90KFDadasGcuWLaO0OXxHjBhBkyZNyM3NJTc3t8T6kSNHkpGRwYIFC1i+vGQGwKhRowB45ZVXWLVqVZF1GRkZjBw5EoA5c+awevXqIuubNGnCiBEjAJg5cyZr164tsv6oo45i2LBhAEyfPp3169cXWd+iRQsGDx4MwNSpU9m0aVOR9a1atWLgwIEATJo0iS1bthRZ36ZNG/r37w/AxIkT2b59e5H17dq1o2/fvgDk5ORQUFBQZH2HDh3o2bMnABMmTKC4tP7sPf0UrN3GKNsLuTnQ9zrPmXRAPRphEnIlz4mjud5mtjku/0fS+n+Z2VYz2wAkCtaZwF/MrNDMPgbmEOKzAE4iXBodbGbvEwKhOwEvSsoFfgW0ids+ClwWL7teAjyZ1Lenkl4vAXIkXUq41JoySU2BZ4CfxiJb2jZXxNH0wg0bNlSm+Qp9WEowNMDuPXU2OMW5ysv/YP9rz5l0SepVWomko4FvAD8AZhEivLqZ2UZJo+Lrq+K2a4BuwC+BpWb2WFz+J+BvwBbgNsJo9CYze1ZSZ+BhM+tRyr4bE4rhaGCkmY1I3o+ZbYzvGxFyOgcD5wGdzazCwhlzMacBL5jZPamcj3SnlfQaO5t1pRTN1s0zmX/92Wnbj3O1JnHvcs/O/csObQzXLPFRZgNS79NKajt02sx2EkKdf0fRy7HJfTwEON7MXgKuA5oBTVM4NgF/AN5KtVhWh9EDssjMKPrcUmZGI0YPyKqlHjmXZnPuCqPKZD7KdFG9KZiEQObX46XSmwijw1RMJowMFwOziaHTiZXxMu0g4CHCvcPhwJ2SFgO5QM+ktnKAvcCMMvbVCPizpKXAIuABM8tPoY+9gO8AZ8cHmHIlfSPF40ubIV1ac8ewzrRunokII8s7hnX2sGhXf+Q9B4W7iy4r3B2WuwavXl2SrW2SrgWamdkNtd0X8ABp55w7EGVdkq1PT8nWKkmTga8AfjPPOefqIS+YaWJmQ2u7D84556pPfbqH6ZxzzlUbL5jOOedcCrxgOueccynwgumcc86lwAumc845lwIvmM4551wKvGA655xzKfCC6VyKPDy7nvPQaFcBL5i1JAZX+/k/SHh4dgPgodGuAj7TTw2K4c8vEMKfuxImi+8MZAJPm9lNcbvuwP3AEcAuoB+wHRgLnAUcDjxkZr+v4UNosDw8u56HZ5/TM4RFe2i0K4cXzJr3H8B/mtlrko42s09jRuYsSacAKwmB05eY2QJJRwE7gO8Bm82su6TDgfmSZphZkf87SLoCuALghBNOqMnjqtc8PLueS471SsR5Daq1JD1XR3laSQ2KI8yXzKxdfH8lobgdChwL/ARYDow3s17FPvs0cAphpAkhS/O/zKysKDFPK0kjD8+uxzw02hVT7wOkDyKfA0hqB1wL9DOzU4BngcblfE7AT8wsO361K69YuvTy8Ox6zEOjXYq8YNaeowjFc7OkLwHnxeV5wLHxPiaSjpR0KOHe5w8lZcTlHSQdUQv9bpA8PLse89BolyK/h1lLzGyxpEWEe5YfAPPj8t2SLgEelJRJuH/ZH3gUaAu8KUnABmBILXS9wRrSpbUXyPro5ysr3sY5/B5mveb3MJ1zrvL8HqZzzjlXBV4wnXPOuRR4wXTOOedS4AXTOeecS4EXTOeccy4FXjCdc865FHjBdM4551LgBdM555xLgRdM55xzLgUNrmBKGiXpN2luc4ikk5Pe3yqpfxrbP0vSZkm58evGdLVdmimL1tFr7GzaXf8svcbO9pBkd3BKpJBs/bi2e+LqiQZXMKvJEGBfwTSzG81sZpr3MS8pqeTWNLe9z5RF6xgzaSnr8ndgwLr8HYyZtNSLpjv4zLkL8t/31BGXNvVu8nVJlwJXA4cB/wJ+BHwXGAPkA4uBXXHbCcA0M3s6vt9mZk3j6+uAS4G9wPNmdr2kHxDyKw8D3gG+A2QDFwB9Jf0KuAi4IdGupH7A3YRzvQD4oZntkrQGeAIYDGQAF5tZrc8CPe6FPHYUhFDkgYft785LU1eSv7g5HTt2pHv37hQUFJCTk1Pi89nZ2WRnZ7N9+3YmTpxYYn23bt3o1KkTmzdvZvLkySXW9+jRg6ysLDZu3Mi0adNKrO/Tpw/t27dn/fr1TJ8+vcT6fv36cfzxx/PBBx8wa9asEusHDhxIq1ateO+995g7d26J9YMGDaJly5bk5eXx6quvllg/dOhQmjVrxrJlyyhtnt4RI0bQpEkTcnNzyc3NLbF+5MiRZGRksGDBApYvX15i/ahRowB45ZVXWLVqVZF1GRkZjBw5EoA5c+awenWR7HCaNGnCiBEjAJg5cyZr164tsv6oo45i2LBhAEyfPp3169cXWd+iRQsGDx4MwNSpU9m0aVOR9a1atWLgwIEATJo0iS1bthRZ36ZNG/r3DxdWJk6cyPbt24usb9euHX379gUgJyeHgoKCIus7dOhAz549AZgwYQLFVepn7y9/hrXbGGV7ITcH+l7n2ZauyurVCFPSScAlQC8zywYKCUXvFqAXcCZJI8Fy2jkPuBA4w8xOBe6KqyaZWfe47C3ge2b2CvAPYHQc/b2b1E5jYAJwiZl1JhTNHybtaqOZnQb8jpCNWZ4ekhZLel5Sx3L6foWkhZIWbtiwoaJDLeHDUkKSAXbvKax0W87VmvwP9r/2bEuXJvUqrUTSVcAvgE/iokQ81hIz+27c5mqgg5ldVdYIU9L/ASvN7JFi7fcFbgOaA02BF8zsylLamQBMA94GHjSzPnF5P+DHZjYsjjB7mdk6SWcAt5tZqfc9JR0F7DWzbZK+AdxvZv9R0fk4kLSSXmNns66Uotm6eSbzrz+7Um05VysS9y737Ny/7NDGcM0SH2W6lDSUtBIBTyTd68sCbi5n+z3EcyDpEMKl1vJMAK6Ko8VbgMZV7O+u+G8h5VweN7MtZrYtvn4OyJDUsor7LtXoAVlkZjQqsiwzoxGjB2RVx+6cS785d4VRZTIfZbo0qG8FcxYwXNIXASQdDSwi3F9sISkDuDhp+zVA1/j6AsK9RIAXgcskNUlqB+BI4KPYzsikdrbGdcXlAW0lfTW+/w4wp7IHJalVDI1G0umE79um8j91YIZ0ac0dwzrTunkmIows7xjW2YOT3cEj7zko3F10WeHusNy5KqhXD/2Y2Yr44M2MOGIsAH5MGGW+SnjoJzfpI48Af5e0GJgOfB7bmS4pG1goaTfwHOFS7w2EB4k2xH8TRfKvwCPxcu/wpP7slHQZ8DdJiYd+xh/AoQ0HfihpD+ES8zetGq+lD+nS2gukO3j9vNafnXP1VL26h+mKOpB7mM4519A1lHuYzjnnXLWoV5dkD3bx8u01xRbPN7Mf10Z/nHPO7ecFsw4xs8eBx2u7H84550ryS7LOOedcCrxgOueccynwgumcc86lwAumc845lwIvmM4551wKvGC6GucB1S6tPCja1RAvmK5GeUC1SzsPinY1xP8OsxZJOgKYCLQBGgG/JgRT30OID9sIjAK2A68DF5hZnqS/ALOLx48dDDyg2gOqIY0B1YW7PSja1RgfYdaugcCHZnaqmXUiTAD/IDDczLoCjxFyMjcDVwETJH0T+EJZxbKqAdLVzQOqXVp5ULSrQT75ei2S1AGYATxFCJz+DHgFeC9u0gj4yMzOjds/DFwEnGpma0u2WFRdnHzdA6pd2nhQtKsmPvl6HWRmq4DTgKXAbYRiuDwpALtzUrE8BDiJcHn2C7XV56rygGqXNh4U7WqYF8xaJOk4YLuZ/RkYB5wBHCOpR1yfIalj3PxnwFvAt4HHY4j1QccDql3aeFC0q2F+SbYWSRpAKJR7CWHXPwT2AA8AzQgPZd0HzAWmAKeb2VZJ9wBbzeym8tqvi5dknXOurivrkqw/JVuLzOwF4IVSVvUpZdlJSZ/772rrlHPOuVL5JVnnnHMuBV4wnXPOuRR4wXTOOedS4AXTOeecS4EXTOeccy4FXjCdc865FHjBdM4551LgBdM555xLgRfMaiDpLEkls6Occ84dtLxgHiQkNap4q5ozZdE6eo2dTbvrn6XX2NkeAO1qRyKxZOvHtd0T1wB4wSQEOUt6VtJiScskXSJpjaQ7JOXGfMnTJL0g6V1JV8bPSdK4+Jmlki4ppe3ukhZJ+oqkrpLmSHojtnVsXP5m0vb/kXgf+3BnfH+xpKslrZC0RNJfa+wEFTNl0TrGTFrKuvwdGLAufwdjJi31oulq3py7IP99TyhxNcLnkg0SQc7nA0hqBtwJvG9m2ZLuBSYAvYDGwDJgPDAMyAZOBVoCCyTNTTQqqSchEPpC4CPgT8CFZrYhFtfbzexySZslZZtZLnAZ8HhS3zaZ2WmxvQ+Bdma2S1LzajkTKRj3Qh47CkLg88DDVu5b/tLUleQvbk7Hjh3p3r07BQUF5OTklPh8dnY22dnZbN++nYkTJ5ZY361bNzp16sTmzZuZPHlyifU9evQgKyuLjRs3Mm1aySvfffr0oX379qxfv57p06eXWN+vXz+OP/54PvjgA2bNmlVi/cCBA2nVqhXvvfcec+fOLbF+0KBBtGzZkry8PF599dUS64cOHUqzZs1YtmwZpU1+P2LECJo0aUJubi65ubkl1o8cOZKMjAwWLFjA8uXLS6wfNWoUAK+88gqrVq0qsi4jI4ORI0cCMGfOHFavXl1kfZMmTRgxYgQAM2fOZO3aorGqRx11FMOGDQNg+vTprF+/vsj6Fi1aMHjwYACmTp3Kpk2biqxv1aoVAwcOBGDSpEls2bKlyPo2bdrQv39/ACZOnMj27duLrG/Xrh19+/YFICcnh4KCgiLrO3ToQM+ePQGY8IeHYe02RtleyM2Bvtd5DqarVj7CDJYC58TRXG8z2xyX/yNp/b/MbKuZbQASBetM4C9mVmhmHwNzgO7xMycBDwODzex9IAvoBLwoKRf4FdAmbvsocFm87HoJ8GRS355Ker0EyJF0KSHVpARJV8QR8cINGzYc0MmoyIelBEAD7N5TWC37c65U+R/sf+05mK4GeLxXJOlo4BvAD4BZwOVANzPbKGlUfH1V3HYN0A34JbDUzB6Ly/8E/A3YQgiEbgzcZGbPSuoMPGxmPUrZd2NCMRwNjDSzEcn7MbON8X0jQpLJYOA8oLOZlVo4ofrivXqNnc26Uopm6+aZzL/+7LTvz7kSEvcu9+zcv+zQxnDNEh9luiorK94rpRGmpIslHRlf/0rSJEmnpbuTtaWUIOdUj20ecImkRpKOIRSz1+O6fOB84A5JZwF5lBEObWY7CTFfv6Po5djkPh4CHG9mLwHXEfIym1buSNNj9IAsMjOKPoOUmdGI0QOyaqM7riGac1cYVSbzUaarZqlekr0hBhefCfQH/kD4n3t90Rl4PV4qvYkwOkzFZMLIcDEwG/gfM9t30ydeph0EPAR0AYYDd0paDOQCPZPayiEESc8oY1+NgD9LWgosAh4ws/wU+5lWQ7q05o5hnWndPBMRRpZ3DOvMkC6ta6M7riHKew4KdxddVrg7LHeumqR0SVbSIjPrIukOwiXIJxPLqr+LDYOka4FmZnZDutqsrkuyzjlXn5V1STbVp2TXSfo9cA5hhHQ4/sBQ2kiaDHwF8BuAzjlXR6VaMEcQ/vTibjPLl3Qs4QEVlwZmNrS2++Ccc658KY0SzWw78Anhzygg/EnD29XVKeecc66uSfUp2ZsIT2aOiYsygD9XV6ecc865uibV+5BDgQuAzwHM7EPgyOrqlHPOOVfXpFowd1t4nNYgzL1afV1yzjnn6p5UC+bE+JRsc0k/AGYSpnNzzjnnGoSUnpI1s7slnUOY8i0LuNHMXqzWnjnnnHN1SEoFU9KdZnYd8GIpy5xzzrl6L9VLsueUsuy8dHakIZL0i6TXWTF7M/G1RdJP47qLJS2XtFdSidknapqHR7syeaCzq8fKLZiSfhjnLs2KocWJr9WEOVRd1ewrmGaWZ2bZZpYNdAW2E+aqhZC/OQwoGc5Ywzw82pXLA51dPVbRJdkngeeBO4Drk5ZvNbNPq61X9YSk0cAuM3sghlCfamZnSzob+DGQGSd8X25mI5M+2g9418z+DWBmb8X2avYASuHh0R4eDWWERx99JAOX5ITUEA90dvVQuSNMM9tsZmvM7Fvxf947CH9a0lTSCTXSw4PbPKB3fN2NcN4y4rIZwI44qhxZ7HPfBP5yIDus7gBpD492ZVr3xv7ILY/acvVQqmklg4F7gOMIU+R9GXjLzDpWb/cObrE45gHZwCRgOfBX4NfA1cDrZta02GcOAz4EOsZ4sOR1LwPXmllKESTVkVbi4dGuVB7o7OqRKgVIE/IhvwasMrN2hEuGr6Wxf/WSmRUAq4FRwCuEEefXga8Cb5XxsfOAN4sXy7rCw6NdqTzQ2TUAqRbMAjPbBBwi6RAze4lwidFVbB5wLeGBnXnAlcCiOHNSQRyFJvsWB3g5tiZ4eLQrlQc6uwYg1XivfElNCf/Tz5H0CXFeWVehecAvgVfN7HNJO+MygIeBJZLeNLORccrBc4D/Sm5A0lDgQeAY4FlJuWY2oOYOoaghXVp7gXRF/Xxlxds4d5BL9R7mEcBOQMBIoBmQE0edro6qjnuYzjlX35V1DzPVqfGSR5NPpK1Xzjnn3EGi3IIp6Z9mdqakrcSkksQqwMzsqGrtnXPOOVdHlFswzezM+K9nXzrnnGvQKnxKVlIjSX5H3znnXINWYcE0s0Igz2f2cc4515Cl+mclXwCWS3qdpD8nMbMLqqVXzjnnXB2TasG8oVp74ZxzztVxqf5ZyZzq7ohzzjlXl6U0NZ6kr0laIGmbpN2SCiVtqe7OFevDWZJK5jmVvf2tkvqn0GbPqveuRLuPSjo5vv5FRdsnfa5pTBp5T9JxxdblSMqTtEzSY6VMqVeneeh0PeRh0a6BSXUu2d8Q5jh9G8gEvg88VF2dSgczu9HMZlaw2VlA2gummX3fzFbEtykVTEmHAhOBPwGjgb9LSv471xzgRKAz+78HBwUPna6nPCzaNTCp3sPEzN6R1Cg+Nfu4pEXAmLK2j9PpTQTaAI0IkVZ3EiYWPw/YA1xBCKf+KjDOzMYrpCTfFbcx4DYze6pY290J87AOB5oToseaAhuBUWb2kaQJwDQze1rSGsIMRYOBDOBiwlR/VwKFki4FfgKsBMYDiSeCf2pm8yXdHJe1j//eF0OhSxyjmT2ViOGK/dsXEg28C3xqZvfF47gd+MTM7gd+DzxvZg/GdYXAXyVdaGYFZrZvFuv48FWbss59XeOh0/UwdLpwN6zdxigPi3YNSKoFc3vMacyVdBfwERWPTgcCH5rZ+QCSmhEK5vtmli3pXmAC0AtoDCwjFKthhPzIU4GWwAJJ+/4vFC+hPghcGPvxJ+BCM9sg6RLgduDyUvqz0cxOk/QjQqbk9yWNB7aZ2d2x7SeBe83sn/HPaF4AToqfP5EQzXUk4c9sflfGMe5jZtdLusrMsuP6toRczPskHUIIij49bvu9Yp+dAkwpfhDxUux3gGtKOUYkXUH4RYQTTqgbfwnkodP1UP4H+18nYrwG3VN7/XGuBqQ6+fqXgY+Bw4CfESZf/62ZvVPOZzoAM4CnCCO9eXGk18vM1km6HOhhZj+I278PnALcBCw1s8fi8j8BfwO2AH8AdgDnmtmHkjoRcibfi7ttBHxkZueWMsJM7PcM4HYz6x9HjskF8xNCeHPCMUAWYbRYYGa3x+3eIqSKNCl+jHH9y8SgZ0nbkkOiJb0I/A/wJeD7Zja8wm9A0fP6CPC5mf20om3ryuTrHjpdz3hYtKvnqhog3ZUwd+wWM7vFzP67vGJJ2HgVcBqwFLhN0o1x1a74796k14n3FY14PyJcSu0S3wtYbmbZ8auzmZ1bxmcT+yosZz+HAF9Laq+1mW0r9vl9bZRzjOV5lBAofRnwWArb7yPpJkIR/+/KfK62eeh0PeNh0a6BSrVgDgZWSfqTpEHxAZVyxac8t5vZn4FxhMKSinnAJXFKvmOAPsDrcV0+cD5wh6SzgDzgGEk94j4zJHVMcT8AWwmXWBNmEO5lJo4hu7wPp3iMxUOiJxMu5XYnXPJNiaTvAwOAb5kV/79V3eah0/WMh0W7BirVv8O8LP5P/zzC07IPSXrRzMp7UrMzME7SXqAA+CHwdAq7mwz0ABYTHvr5HzNbL+nE2JePJQ0CnifcqxwOPBDvHx4K3Ed4wCYVU4GnJV1IKJRXx2NbEtuaS3gwqDLHWFyRkGgz2y3pJSA/PkCVqvHAv4FXw3NRTDKzWyvx+VrlodP1iIdFuwYqpXuY+zYORXMg4XJiHzNrWV0dq6/iwz5vAheb2dvVua+6cg/TOecOJlW6hynpvPgQzdvARYT7cK3S2sMGIE5m8A4wq7qLpXPOufRK9c9Kvkt4EvS/zGxXRRu70sXJDNrXdj+cc85VXqr3ML8l6UvAOfH+2etm9km19sw555yrQ1K9JHsx4UnVi4ERwL8kVervB51zzrmDWaqXZH8FdE+MKuOfe8wktadenXPOuYNeqn+HeUixS7CbKvFZ55xz7qCX6ghzuqQXCBOnA1wC+F8pO+ecazDKLZiSvgp8ycxGSxoGnBlXvUqIm3LOOecahIouq95HmPQcM5sU55D9b8JsPPdVb9cOHskh0ZKyJOUmfW2R9NO47mJJyyXtlVTij2IlHSvpHUlvSjoyaXkTSc9KWhk/P7ZGDsw559w+FRXML5nZ0uIL47K21dKjg9O+gmlmeYnJ2wmT1m8n/IIBIcJsGGHKvSJigZwCXEfI7ny62By0d5vZiYSJ53tJOq8ajmOfKYvW0WvsbNpd/yy9xs72sGdXVCKxZOvHtd0T52pMRQWzeTnrMtPYjzpN0mhJV8fX90qaHV+fLekZYki0pOKXqfsB75rZvwHM7C0zyyul/QzC/eE7zeyZGCj9D+CR+LntZvZSfL2bMLVetQVIT1m0jjGTlrIufwcGrMvfwZhJS71ouv3m3AX573tCiWtQKnroZ6GkH5jZI8kLY3LGG9XXrTpnHvBz4AGgG3B4LHK9CQknAxIh0cV8k/0PSpXJzAqAQcWWPVTatpKaE9Jj7k+9+5Uz7oU8dhSEeeEHHrZ/ou2Xpq4kf3FzOnbsSPfu3SkoKCAnp+St7OzsbLKzs9m+fTsTJ04ssb5bt2506tSJzZs3M3ny5BLre/ToQVZWFhs3bmTatGkl1vfp04f27duzfv16pk+fXmJ9v379OP744/nggw+YNWtWifUDBw6kVatWvPfee8ydW2Kwz6BBg2jZsiV5eXm8+uqrJdYPHTqUZs2asWzZMkqbq3fEiBE0adKE3NxccnNzS6wfOXIkGRkZLFiwgOXLS+YEjBo1CoBXXnmFVatWFVmXkZHByJEjAZgzZw6rV68usr5JkyaMGDECgJkzZ7J27doi64866iiGDRsGwPTp01m/fn2R9S1atGDw4MEATJ06lU2bNhVZ36pVKwb2yobcnBDplZsDfa/zHEzXIFRUMH8KTJY0kv0FshshSHpoNfarrnkD6CrpKEIu5puE89CbkHBSgqTDgAuAMenqRIxV+wvwgJm9V8Y2VwBXAJxwwgkHtJ8PSwl7Bti9pzLhKq7eSs7DTORgDrqndvvkXA1IKa1E0teBTvHtcjObXa29qoMkzQL+DrQElgAdCIWpHbDVzJoW2/5C4MelBVpLehm41swqFSUi6TFgm5mVWqSLO9C0kl5jZ7OulKLZunkm868/u9LtuXokce9yz879yw5tDNcs8VGmqzeqlFZiZi+Z2YPxq8EVy2gecC3hgZ15hJzMRRZ+4ygeEg0hN7TCy7GpknQb0Iww6q9WowdkkZnRqMiyzIxGjB6QVd27dnVd8ugyITHKdK6e89l6UjcPOBZ41cw+BnbGZbA/JDoHQNIRwDnApOQGJA2VtJYQkP1snAyiQpLaAL8ETgbejA8YlRfeXSVDurTmjmGdad08ExFGlncM6+wB0A7ynoPC3UWXFe4Oy52r5yoVIO0OLh4g7ZxzlVelS7LOOedcQ+cF0znnnEuBF0znnHMuBV4wnXPOuRR4wXTOOedS4AXTOeecS4EXTOeccy4FXjCdc865FHjBdM4551LgBdM1WA02JNvDn507IA2uYEoaJek3aW5ziKSTk97fKql/GtsfKWmJpKWSXpF0arrabqgadEi2hz87d0AqysN0qRkCTANWAJjZjWlufzXQ18w+k3QeYbL3M9K8jwalwYZkF+6GtdsY5eHPzlVavRthSrpU0usx0eP3khpJukzSKkmvA72Stp0gaXjS+21Jr6+LI7rFksbGZT+QtCAue0ZSE0k9CUHR4+I+v5LcrqR+khbFth6TdHhcvkbSLZLejOtOLOuYzOwVM/ssvn0NaFPO8V8haaGkhRs2bDigc9gQNNiQ7PwP9r/2WC7nKqVepZVIOgm4CxhmZgWSfgv8C/g10BXYDLxEyLG8StIEYJqZPR0/v83MmsZR3A1AfzPbLuloM/tUUgsz2xS3vQ342MweLKWdCYQR5zTgbaCfma2S9EfgTTO7T9Ia4P/i538EnGZmFUZ2SboWODGVbT2tpGwNMiTbw5+dS0lDSSvpRyiMCyTlxvc/A142sw1mtht4KoV2+gOPm9l2ADP7NC7vJGmepKXASKBjBe1kAavNbFV8/wTQJ2l9Ii/zDaBtRZ2S9HXge8B1KRyDK0eDDMn28GfnqqS+FUwBT5hZdvzKAm4uZ/s9xHMg6RDgsAranwBcZWadgVuAxlXs7674byEV3E+WdArwKHBhYpTrDlyDDMn28GfnqqS+PfQzC/i7pHvN7BNJRwOLgPsltQC2ABcDi+P2awgj0omE+5AZcfmLwI2ScpIvyQJHAh9JyiCMMBOPVG6N64rLA9pK+qqZvQN8B5hT2YOSdAJhNPqdpNGqq6IhXVrX7wJZ3M9XVryNc65M9WqEaWYrgF8BMyQtIRS+YwmjzFeB+cBbSR95BOgraTHQA/g8tjMd+AewMF7avTZufwPhnuh8IPn/Pn8FRseHe76S1J+dwGXA3+Jl3L3A+AM4tBuBFsBv44NFfmPSOedqWL166McV5Q/9OOdc5TWUh36cc865alHf7mEe1CRdBlxTbPF8M/txbfTHOefcfl4w6xAzexx4vLb74ZxzriS/JOucc86lwAumc845lwIvmM4551wKvGA655xzKfCC6ZxzzqXAC2aaSDorRn0l3l8p6btpbH+cpJUxSHqypObparu6TFm0jl5jZ9Pu+mfpNXZ2wwhnTiSCbP24tnvinEszL5jpcxawr2Ca2Xgz+2Ma238R6GRmpwCrgDFpbDvtpixax5hJS1mXvwMD1uXvYMykpfW/aM65C/Lf9wQQ5+oh/zvMCkiaAhxPSCa538weljQQ+F+gEbCRELl1JVAo6VLgJ4RosW2ETMw/mtnpsb22wFQz6yypK3AP0DS2M8rMPiqtH2Y2I+nta8Dw0rarK8a9kMeOghDGPPCw/dPuvjR1JfmLm9OxY0e6d+9OQUEBOTk5JT6fnZ1NdnY227dvZ+LEiSXWd+vWjU6dOrF582YmT55cYn2PHj3Iyspi48aNTJs2rcT6Pn360L59e9avX8/06dNLrO/Xrx/HH388H3zwAbNmzSqxfuDAgbRq1Yr33nuPuXPnhoWFu2HtNkbZXsjNgb7Xec6kc/WIjzArdrmZdQW6AVdL+hJh0vaLzOxU4GIzW0OYVP3eGCs2L/FhM1sJHCapXVx0CfBUTDx5EBge238MuD3VPgHPl7ZC0hWSFkpauGHDhkofbLp8WEo4M8DuPYU13JMalP/B/teeM+lcveOTr1dA0s3A0Pi2LXA3cKKZjSxlu21mdnfx95J+Aew1s7GS3iQUzcOBV4D3YhONgI/M7NwK+vNLQvEeZhV882pz8vVeY2ezrpSi2bp5JvOvP7sWelTNEvcu9+zcv+zQxnDNEh9lOneQ8cnXD4Cks4D+QI84mlwE5B5AU08BIyR1AMzM3iaEXS9PCrvunEKxHAUMAkZWVCxr2+gBWWRmNCqyLDOjEaMHZNVSj6rZnLvCqDKZjzKdq1e8YJavGfBZDJE+Efga4V5mn8Ql1hhSDWWHSGNm7wKFhDzNp+LiPOAYST1iOxmSOpbVkXjf9H+AC8xse5WPrJoN6dKaO4Z1pnXzTEQYWd4xrHP9DWzOey7cw0xWuDssd87VC35JthySDgemEC7F5gHNCWHUmYSHfg4BPjGzc+Lo8WlCSPS+h36SLtFeC4wD2sV7nkjKBh4gFOZDgfvM7JEy+vIO4TLuprjoNTO7srz+ex6mc85VXlmXZL1g1mNeMJ1zrvL8HqZzzjlXBf53mHWMpIeAXsUW3x+zMp1zztUSL5h1jJn9uLb74JxzriS/JOucc86lwAumc845lwIvmM4551wKvGA655xzKfCC6ZxzzqXAC6ZzNRj63CBDtZ2rJ7xgViNJz0lqnqa2xklaKWmJpMnpatdRY6HPDTZU27l6wqfGqyJJh5rZnhrYz7nAbDPbI+lOADO7rrzP+NR4KYijywl7BoMOgTbdoNFhAGkPuV70fj67Yh7o9N0n7tuu3kaeOXeQ8qnxkki6VNLrknIl/V7SGXHk1ljSEZKWS+oUXz8Wt10k6cL4+VGS/iFpNjBLUlNJj0taGtu5KG63RlLL2M6zkhZLWibpkri+q6Q5kt6Q9IKkY8vqs5nNSCrMrwFtyji2OhEgfdAoHsuVHAKdZmWFZ5cVtu2cq1sa3AhT0knAXYQA5gJJvyUUoA6E6K5MYK2Z3SHpf4EVZvbneAn0daALcDFwG3CKmX0aR3yHm9lP4z6+YGafSVpDCHvuCww0sx/E9c2A7cAc4EIz2xCL6AAzuzyFY5gKPGVmfy5vOx9hVqCGQ58bXKi2cwcpH2Hu1w/oCiyQlBvftwduBc4hFLi74rbnAtfH7V4mFNQT4roXzezT+Lo/8FBiB2b2WbF9LgXOkXSnpN5mthnIAjoBL8b2f0UZo8Zkkn4J7AFKXiN0lVPDoc8NLlTbuXqmIc4lK+AJMxtTZGG4HNoUyCAUxs/jtheZWV6xbc+I61NiZqsknQZ8A7hN0ixgMrDczHqk3HFpFDAI6GcN7dJAdSgv9HnQPWnfXSI8e9wLeXyYv4PjmmcyekBW/Q3Vdq6eaYiXZE8G/g70MrNPJB0NHAk8CPwVaAcca2ZXxUuyRwE/MTOT1MXMFsXC1c3MroptjgUal3NJ9jDgUzPbKWkQ8H1gBLAC+I6ZvSopA+hgZsvL6PdA4B6gr5mldHPSL8k651zllXVJtsGNMM1shaRfATMkHQIUEApogZk9KakR8Iqks4FfA/cBS+K2qwkjvOJuAx6StAwoBG4BJiWt7wyMk7Q37u+HZrZb0nDggXhP89C4r1ILJvAb4HDCJVyA18zsygM9D8455yqnwY0wGxIfYTrnXOX5Qz/OOedcFTS4S7J1naSHgF7FFt9vZo/XRn+cc84FXjDrGDP7cW33wTnnXEl+SdY555xLgRdM55xzLgVeMJ1zzrkUeMF0zjnnUuAF0znnnEtBgyuYMZrrN2luc0icci/x/lZJ/dPY/oUxNiw3Rnedma62yzNl0Tp6jZ1Nu+ufpdfY2Qdf0HEijWTrx7XdE+dcPdDgCmY1GQLsK5hmdqOZzUxj+7OAU80sG7gceDSNbZdqyqJ1jJm0lHX5OzBgXf4OxkxaenAVzTl3Qf771ZY+4pxrWOrd32FKuhS4mjDh+b+AHwHfBcYA+cBiYFfcdgIwzcyeju+3mVnT+Po64FJgL/C8mV0v6QfAFbHtd4DvANnABUDfOEftRcANiXYl9QPuJpzrBYR5ZHfFidmfAAYTElIuNrOVpR2TmW1LensEUO3zGY57IY8dBSHweOBh+7v10tSV5C9uTseOHenevTsFBQXk5JRMGsvOziY7O5vt27czceLEEuu7detGp06d2Lx5M5MnTy6xvkePHmRlZbFx40amTZtWYn2fPn1o374969evZ/r06SXW9+uRzfG5OSGuKzcH+l5XLRmXzrmGo16NMGM49CWEJJJswkTolxImQ+8FnEnSSLCcds4DLgTOMLNT2Z+POcnMusdlbwHfM7NXgH8Ao80s28zeTWqnMTABuMTMOhOK5g+TdrXRzE4DfgdcW0GfhkpaCTxLGGWWtd0V8bLtwg0bUgo1KdWHpQQdA+zeU3jAbdaoRTn7sy6rMePSOddw1KvJ1yVdBfwC+CQuygR2AEvM7Ltxm6sJMVpXlTXClPR/wEoze6RY+30JySTNCdmZL5jZlaW0MwGYBrwNPGhmfeLyfsCPzWxYHGH2MrN1MV/zdjOr8L6npD7AjalsW5XJ13uNnc26Uopm6+aZzL/+7ANqs8Yk7l3u2bl/2aGN4ZolPsp0zlWooUy+ngiHzo5fWcDN5Wy/h3gOYnzXYRW0PwG4Ko4WbyEETVfFrvhvISleHjezuUB7SS2ruO9yjR6QRWZGoyLLMjMaMXpAVnXuNj3m3LV/dJngo0znXBXVt4I5Cxgu6YsAMRx6EeH+YosY0nxx0vZrgK7x9QWEe4kALwKXSWqS1A6EoOmPYjsjk9rZGtcVlwe0lfTV+P47wJzKHpSkryqGYEo6jZCLuamy7VTGkC6tuWNYZ1o3z0SEkeUdwzozpEvr6txteuQ9B4W7iy4r3B2WO+fcAapXD/2UEQ79Y8Io81XCQz+5SR95BPi7pMXAdODz2M50SdnAQkm7gecIl3pvIDxItCH+myiSfwUeiZd7hyf1Z6eky4C/SUo89DP+AA7tIuC7kgoIl5gvsRq4lj6kS+uDo0AW9/NSn51yzrkqqVf3MF1RHiDtnHOV11DuYTrnnHPVol5dkj3Yxcu31xRbPN8zMp1zrvZ5waxDzOxx4PHa7odzzrmS/JKsc845lwIvmM4551wKvGA655xzKfCC6ZxzzqXAC6ZzzjmXAi+YrlLqVKi0B0Q752qQF8xqJOk5Sc3T1NbFkpZL2iupxAwUNaHOhUp7QLRzrgb532FWkaRDzWxPaevM7Btp3NUyYBjw+zS2WSl1KlT675Ng7TZGeUC0c66GNMgRpqRLJb0uKVfS7yWdIWmJpMaSjogjuU7x9WNx20WSLoyfHyXpH5JmA7MkNZX0uKSlsZ2L4nZrJLWM7TwrabGkZZIuieu7Spoj6Q1JL0g6tqw+m9lbZpaXwrGlJUC6NHUqVDr/g/2vPbrLOVcDGtzk65JOAu4ChplZgaTfAq8BHQj5lpnAWjO7Q9L/AivM7M/x0urrQBdCRNhtwClm9qmkO4HDzeyncR9fMLPPYkh0N6AvMNDMfhDXNwO2E6K+LjSzDbGIDjCzyyvo/8vAtWZW4azq6Z58vc6ESntAtHOuGvnk6/v1I2RgLpCUG9+3B24FziEUuLvitucC18ftXiYU1BPiuhfN7NP4uj/wUGIHZvZZsX0uBc6RdKek3ma2GcgCOgEvxvZ/BbRJ21FWgzoTKu0B0c65WtAQ72EKeMLMxhRZGC6HNiWESDcmZGMKuKj4pVBJZ8T1KTGzVTH4+RvAbZJmAZOB5WbWoyoHU5MS2ZjjXsjjw/wdHNc8k9EDsmo+M7O8gOhB99RsX5xzDUZDLJizCKHR95rZJ5KOJgRBP0gIiG4H3AlcBbwA/ETST8zMJHUxs0WltPkiIaj6p7D/kmxipaTjgE/jpd184PvAWOAYST3M7FVJGUAHM1teTcedFnUiVNoDop1ztaDBFUwzWyHpV8AMSYcABcDfgQIze1JSI+AVSWcDvwbuA5bEbVcDg0pp9jbgIUnLgELgFmBS0vrOwDhJe+P+fmhmuyUNBx6I9zQPjfsqtWBKGkoo6scAz0rKNbMBVTkXzjnnUtfgHvppSNL90I9zzjUE/tCPc845VwUN7pJsXSfpIaBXscX3x3Bp55xztcQLZh1jZj+u7T4455wryS/JOueccynwgumcc86lwC/JOudcLSsoKGDt2rXs3Lmz4o1d2jRu3Jg2bdqQkZGR0vZeMJ1zrpatXbuWI488krZt2yKptrvTIJgZmzZtYu3atbRr1y6lz/glWeecq2U7d+6kRYsWXixrkCRatGhRqVF9gyuYMZrrN2luc4ikk5Pe3yqpfzr3EdvtLmlPnCGoWkxZtI5eY2fT7vpn6TV2ds2EQyfSR7Z+XP37cq6O8mJZ8yp7zhtcwawmQ4B9BdPMbjSzmencQZyy705gRjrbTTZl0TrGTFrKuvwdGLAufwdjJi2t/qI55y7If9/TRpw7yHz/+99nxYoV5W4zatQonn766RLL16xZw5NPPlnpfZbVXk2od/cwJV0KXA0cBvwL+BHwXWAMkA8sBnbFbScA08zs6fh+m5k1ja+vAy4F9gLPm9n1kn4AXBHbfgf4DpANXAD0jXPUXkSYxH2amT0tqR9wN+FcLyDMI7srZmU+AQwmJKRcbGblzSr+E+AZoHvVzlDZxr2Qx46CEAY98LD9XXlp6kryFzenY8eOdO/enYKCAnJyckp8Pjs7m+zsbLZv387EiRNLrO/WrRudOnVi8+bNTJ48OSws3A1rtzHK9kJuDvS9zjMtnTtIPProowf82UTB/Pa3v53GHlWvejXCjOHQlwC9zCybMBH6pYTJ0HsBZ5I0EiynnfOAC4EzzOxU9udjTjKz7nHZW8D3zOwV4B/AaDPLNrN3k9ppDEwALjGzzoSi+cOkXW00s9OA3wHXltOf1sDQuF1Ffb9C0kJJCzds2FDR5kV8WEo4NMDuPYWVaqdS8j/Y/9ozLZ1LSbpvnYwbN44HHngAgJ/97GecfXYIhJ89ezYjR45kxowZ9OjRg9NOO42LL76Ybdu2AXDWWWeRmK/6D3/4Ax06dOD000/nBz/4AVddddW+9ufOnUvPnj1p3779vtHh9ddfz7x588jOzubee++lsLCQ0aNH0717d0455RR+//vfA+HhnKuuuoqsrCz69+/PJ598UqVjrYr6NsJMDocGyAR6Ai+b2QYASU8BHSpopz/wuJltB0gKiu4k6TagOSE784UK2skCVpvZqvj+CUIM2H3xfSLR5A1gWDnt3AdcZ2Z7K7rmbmYPAw9DmHy9gv4VcVzzTNbFojl994n7lrdunsnvRp29731GRgajRo0qs50mTZqUu75Zs2ZhfeLepcWb7oW7fZTpXAUSt04SV4MSt06AA47e6927N//3f//H1VdfzcKFC9m1axcFBQXMmzePU045hdtuu42ZM2dyxBFHcOedd3LPPfdw44037vv8hx9+yK9//WvefPNNjjzySM4++2xOPfXUfes/+ugj/vnPf7Jy5UouuOAChg8fztixY7n77ruZNm0aAA8//DDNmjVjwYIF7Nq1i169enHuueeyaNEi8vLyWLFiBR9//DEnn3wyl19++YGeviqpVyNM9odDZ8evLODmcrbfQzwHMb7rsAranwBcFUeLtxCCpqtiV/y3kPJ/eekG/DVexh0O/FbSkCruu4TRA7LIzGhUZFlmRiNGD8hK966COXeFUWUyH2U6V67kWycJOwoKGfdCXhmfqFjXrl1544032LJlC4cffjg9evRg4cKFzJs3j8zMTFasWEGvXr3Izs7miSee4N///neRz7/++uv07duXo48+moyMDC6++OIi64cMGcIhhxzCySefzMcfl/5w34wZM/jjH/9IdnY2Z5xxBps2beLtt99m7ty5fOtb36JRo0Ycd9xx+0a/taG+jTBLC4deBNwvqQWwBbiYcB8TYA1hRDqRcB8y8derLwI3Ssoxs+2Sjo6jzCOBj2LY80ggcR1ka1xXXB7QVtJXzSxxz3NOZQ/KzPb9kVDSfdcplW2nIonfTse9kMeH+Ts4rnkmowdkVV9gdN5zYVSZrHB3WD7onurZp3MHubJunZS1PBUZGRm0a9eOCRMm0LNnT0455RReeukl3nnnHdq1a8c555zDX/7ylwNu//DDD9/3uqxISTPjwQcfZMCAojG/zz333AHvN93q1QjTzFYAiXDoJYTCdyxhlPkqMJ9w7zHhEcLDOouBHsDnsZ3phPuSCyXlsv/+4g2EB4nmA8kP6PwVGC1pkaSvJPVnJ3AZ8DdJSwkPEI1P4yGn3ZAurZl//dmsHns+868/u/qKJcDPV8LNm0t+/by8Z5+ca9iOa55ZqeWp6t27N3fffTd9+vShd+/ejB8/ni5duvC1r32N+fPn88477wDw+eefs2rVqiKf7d69O3PmzOGzzz5jz549PPPMMxXu78gjj2Tr1q373g8YMIDf/e53FBQUALBq1So+//xz+vTpw1NPPUVhYSEfffQRL730UpWOsyrq2wgTM3sKeKrY4teAEvFYZvYx8LWkRdclrRsLjC22/e8o5cEbM5tP0YeJRiWtmwV0KeUzbZNeLwTOKnk0JZnZqAo3cs7VW6MHZBW5hwnpuXXSu3dvbr/9dnr06MERRxxB48aN6d27N8cccwwTJkzgW9/6Frt2hbtIt912Gx067H8UpHXr1vziF7/g9NNP5+ijj+bEE0+kWbNm5e7vlFNOoVGjRpx66qmMGjWKa665hjVr1nDaaadhZhxzzDFMmTKFoUOHMnv2bE4++WROOOEEevToUaXjrAqVNTx2B79u3bpZ4gk251zd9dZbb3HSSSelvP2URetq7tZJirZt20bTpk3Zs2cPQ4cO5fLLL2fo0KG12qdUlHbuJb1hZt2Kb1vvRpgHM0mXAdcUWzzfMzKdc8mGdGld6wWyuJtvvpmZM2eyc+dOzj33XIYMGVLbXUo7L5h1iJk9TimXjp1zrq67++67a7sL1a5ePfTjnHPOVRcvmM4551wKvGA655xzKfCC6ZxzzqXAC6Zzzrka07Nnzwq3adu2LRs3biyx/OWXX+aVV16p9D7Laq+yvGCmiaSzJPVMen+lpO+msf2LJS2XtFdSib8POqhVY4B0rQRiO+fKdCAFL+FAC2a6eMFMn7MIySgAmNl4M/tjGttfRkg0mZvGNuuGagqQrrVAbOdqQhp/0VyzZg2dOnXa9/7uu+/m5ptv5qyzzuK6667j9NNPp0OHDsybNw+A888/nyVLlgDQpUsXbr31VgBuvPFGHnnkESBEhiWium666aZ9bTdt2hSAvXv38qMf/YgTTzyRc845h2984xtFgqEffPBBTjvtNDp37szKlStZs2YN48eP59577yU7O5t58+axYcMGLrroIrp370737t2ZP38+AJs2beLcc8+lY8eOfP/73y9z/trK8r/DrICkKcDxhGSS+83sYUkDgf8FGgEbge8BVwKFMcD6J4SosW3ANOCPZnZ6bK8tMNXMOkvqCtxDiArbCIwys49K64eZvRU/X01HWku2rofcHCbYRfDGNvj4YWgUQmOqGlj9j3WZ7ChoyhHazecW2kykOtS1P/p2rtKSf9GsxrCCPXv28Prrr/Pcc89xyy23MHPmTHr37s28efP48pe/zKGHHrqvUM2bN4/x48czY8YM3n77bV5//XXMjAsuuIC5c+fSp0+ffe1OmjSJNWvWsGLFCj755BNOOumkIrFdLVu25M033+S3v/0td999N48++ihXXnklTZs25dprw/Te3/72t/nZz37GmWeeyfvvv8+AAQN46623uOWWWzjzzDO58cYbefbZZ/nDH/6QlnPhBbNil5vZp5IyCTmbfydM2t7HzFYnkkwkjQe2mdndAJL6AZjZSkmHSWpnZqsJAddPxcSTB4ELzWyDpEuA24EqBb1JugK4AuCEE06oSlM1o3jEV/4H0OIrZW9fCZ9t313q8qqkOjhXJ8RfNLG91Z4hO2xYiOrt2rUra9asAcK8sw888ADt2rXj/PPP58UXX2T79u2sXr2arKwsHnnkEWbMmEGXLmEa7W3btvH2228XKZj//Oc/ufjiiznkkENo1aoVX//618vc76RJkyjNzJkzWbFixb73W7ZsYdu2bcydO3ffZ84//3y+8IUvpOVceMGs2NWSEhMiHk8oRnNj8UsOly7PREKhHBv/vYQQLt0JeDGOGhsBpY4uK6MqAdI1LvEffeFuRvE3MGBHY7hoSZH/+A80sPqRsbMhf8e+0WVCVVMdnKt1yb9oJjJkqzDKPPTQQ9m7d/8vrjt37tz3OhHN1ahRI/bs2QOEdJKFCxfSvn17zjnnHDZu3MgjjzxC165dQ5fMGDNmDP/1X/91wH0qbb/F7d27l9dee43GjasaTZwav4dZDklnAf2BHmZ2KiFbM/cAmnoKGCGpA2Bm9jYh7Hp5Uth1ZzM7Nz09P0hUc4B0jQdiO1cTkn7RBMK/uTlVupf5pS99iU8++YRNmzaxa9cupk2bVu72hx12GMcffzx/+9vf6NGjR5FoMAhRXY899hjbtm0DYN26dXzyySdF2ujVqxfPPPMMe/fu5eOPP+bll1+usJ/FI8HOPfdcHnzwwX3vc3NzAejTpw9PPvkkAM8//zyfffZZhW2nwgtm+ZoBn8UQ6RMJUWCNgT6S2gHEkGooO0QaM3sXKCTkaSaix/KAYyT1iO1kSOpYbUdSF5UXIJ0GQ7q05o5hnWndPBMBrZtncsewzn7/0h3cquEXzYyMDG688UZOP/10zjnnHE488cQKP9O7d2+++MUvkpmZSe/evVm7di29e/cGQiH79re/TY8ePejcuTPDhw8vUugALrroItq0acPJJ5/MpZdeymmnnVZhJNjgwYOZPHnyvod+HnjgARYuXMgpp5zCySefzPjxIW74pptuYu7cuXTs2JFJkyal7faUx3uVQ9LhwBSgLaHANSeEUWcSHvo5BPjEzM6Jo8enCSHR+x76SbqneS0wDmhnZmvismzgAUJhPhS4z8weKaMvQwn3PI8B8oFcMxtQ2rYJHu/l3MGhUvFe/3cibC3l7s2Rxx504euJSLBNmzZx+umnM3/+fFq1alWjffB4rzQxs13AeWWsfr7YtquAU5IWzSu2/m7g7mLLcoE+pMDMJgOTU9nWOVePHWRFsTyDBg0iPz+f3bt3c8MNN9R4sawsL5jOOedqRSr3LesSL5h1jKSHgF7FFt8fszKdc87VEi+YdYyZ/bi2++Ccq3lmVv8mJqnjKvsMjz8l65xztaxx48Zs2rQpbVO4uYqZGZs2barU33D6CNM552pZmzZtWLt2LRs2bKjtrjQojRs3pk2bNilv7wXTOedqWUZGBu3atavtbrgK+CVZ55xzLgVeMJ1zzrkUeMF0zjnnUuBT49VjkjYA/z6Aj7Yk5HPWNd6vyqurffN+VU5d7RfU3b5VpV9fNrNjii/0gulKkLSwtHkUa5v3q/Lqat+8X5VTV/sFdbdv1dEvvyTrnHPOpcALpnPOOZcCL5iuNA/XdgfK4P2qvLraN+9X5dTVfkHd7Vva++X3MJ1zzrkU+AjTOeecS4EXzHpO0kBJeZLekXR9KesPl/RUXP8vSW2T1o2Jy/MkDUi1zersl6RzJL0haWn89+ykz7wc28yNX1+s4b61lbQjaf/jkz7TNfb5HUkP6ABiKarQr5FJfcqVtFdSdlxX5XOWQr/6SHpT0h5Jw4ut+09Jb8ev/0xaXhPnq9R+ScqW9Kqk5ZKWSLokad0ESauTzld2ZftVlb7FdYVJ+/9H0vJ28fv+Tvw5OKym+iXp68V+xnZKGhLX1dQ5+29JK+L3bJakLyetS8/PmZn5Vz39AhoB7wLtgcOAxcDJxbb5ETA+vv4m8FR8fXLc/nCgXWynUSptVnO/ugDHxdedgHVJn3kZ6FaL56wtsKyMdl8HvgYIeB44r6b6VWybzsC76TpnKfarLXAK8EdgeNLyo4H34r9fiK+/UIPnq6x+dQD+I74+DvgIaB7fT0jetqbPWVy3rYx2JwLfjK/HAz+syX4V+75+CjSp4XP29aR9/pD9/12m7efMR5j12+nAO2b2npntBv4KXFhsmwuBJ+Lrp4F+8besC4G/mtkuM1sNvBPbS6XNauuXmS0ysw/j8uVApqTDK7n/aulbWQ1KOhY4ysxes/Bf6R+BIbXUr2/Fz6ZLhf0yszVmtgTYW+yzA4AXzexTM/sMeBEYWFPnq6x+mdkqM3s7vv4Q+AQo8UfsVVCVc1aq+H0+m/B9h/BzMKSW+jUceN7Mtldy/1Xt20tJ+3wNSMSQpO3nzAtm/dYa+CDp/dq4rNRtzGwPsBloUc5nU2mzOvuV7CLgTTPblbTs8XjZ54YDuYyXhr61k7RI0hxJvZO2X1tBm9Xdr4RLgL8UW1aVc1aVn4fyfsZq4nxVSNLphBHNu0mLb4+X/e49wF/Wqtq3xpIWSnotcdmT8H3Oj9/3A2kzHf1K+CYlf8Zq+px9jzBiLO+zlf4584LpDkqSOgJ3Av+VtHikmXUGesev79Rwtz4CTjCzLsB/A09KOqqG+1AmSWcA281sWdLi2j5ndVYcgfwJuMzMEiOqMcCJQHfCJb7raqFrX7Ywg823gfskfaUW+lCqeM46Ay8kLa7RcybpUqAbMC7dbXvBrN/WAccnvW8Tl5W6jaRDgWbApnI+m0qb1dkvJLUBJgPfNbN9v/mb2br471bgScJlnMo64L7Fy9ebYh/eIIxKOsTtk1Nqa/ycRSV+80/DOavKz0N5P2M1cb7KFH/ReRb4pZm9llhuZh9ZsAt4nOr7GStT0vfsPcI96C6E73Pz+H2vdJvp6Fc0AphsZgVJ/a2xcyapP/BL4IKkK0/p+zmryo1Y/6rbX4SA8PcID+0kbpR3LLbNjyn6oMjE+LojRR/6eY9w473CNqu5X83j9sNKabNlfJ1BuJdzZQ2fs2OARvF1+/gf39HxffGHC75RU/2K7w+J/WmfznNWmZ8Hij38QRhtrCY8iPGF+LrGzlc5/ToMmAX8tJRtj43/CrgPGFsdP2Pl9O0LwOHxdUvgbeLDL8DfKPrQz49qql9Jy18Dvl4b54zwi8O7xAe2quPnrFKd9q+D7wv4BrAq/iD9Mi67lfAbGEDj+B/aO/GHJ/l/qL+Mn8sj6emx0tqsqX4BvwI+B3KTvr4IHAG8ASwhPAx0P7F41WDfLor7zgXeBAYntdkNWBbb/A1x0pAa/F6eBbxWrL20nLMU+tWdcH/oc8JIaHnSZy+P/X2HcOmzJs9Xqf0CLgUKiv2MZcd1s4GlsW9/BppW089YWX3rGfe/OP77vaQ228fv+zvx5+DwGv5etiX8UnZIsTZr6pzNBD5O+p79I90/Zz7Tj3POOZcCv4fpnHPOpcALpnPOOZcCL5jOOedcCrxgOueccynwgumcc86lwAumc2WQNESSSToxadlZkqaloe0JxVMoStnmLEk9K9luE0k5MYFhmaR/SmoqqbmkH1Wt10X201bSsoq3LLeNUs9BXL5d0pFJy+6L34uWVdlnBf05S9LmOLVhnqS5kgYlrb9S0nfj6xPjdIKLJH1F0tWS3pKUU139c7XPC6ZzZfsW8M/4b204i/B3d5VxDfCxmXU2s06EOTULCBM+pK1gVlbSDDSpeoc4ubakQwgTi1d21pkDMc/MuphZFnA18BtJ/QDMbLyZ/TFuNwR4Om77LuHcnmNmI1PZyQGcD1cHeMF0rhSSmgJnEgrON4utPkrSs3EUMl7SIZIaxZHRsji6+1lsJztOkr1E0mRJXyhlX2sSIydJ3RQyKtsCVwI/iyOZ3pKOkfSMpAXxq1cpXT+WpMJiZnkWpggbC3wltjUujjpnKWQbLpWUKE5t40jpEYU8yBmSMuO6rpIWS1pMmFWIpM/Mi229mRgVxxHbPIXMxhUKfhPP20zChBNl+SthongIvzjMBxITiyPpUkmvx+P5vaRGcfnvFCYmXy7plmLn+Jak4z2RCphZLuEP46+Kbdws6VpJ3wB+CvxQ0ksKuaftgecl/UzSEZIei/1blHRuR0n6h6TZwKwKtpskabpCfuNdSccxMB7DYkmz4rJS23HV4EBmXPAv/6rvX8BI4A/x9StA1/j6LGAn4X+QjQhRQcOBroQIocTnm8d/lwB94+tbgfvi6wnEqcWANeyfoq4b8HJ8fTNwbVKbTwJnxtcnAG+V0u9sQhzVq8Bt7M91bEtSVidhqrGj4uuWhBGd4nZ72D+zzUTg0qRj6RNfj0u0BzQBGsfX/wEsTDpXnwPt4vth8Xw1IuRM5lP69GoT4jl9jTCV2SNA38R5Ak4CpgIZcfvfEuYVhv1TnjUizLN6StI5/kl8/SPg0VL2exYwrZTz+Vbx70cp35vk7+H/Jp2z5oTZaY4ARhFmyTk6he3eI8wF3Bj4N2Eu1GMIqRvtih1rqe3U9n9D9fHLLws4V7pvEaaKgzDa+RZhGjmA1y1MfI2kvxBGorOA9pIeJEzaPUNSM0LhnBM/9wRhyrID1R84WfsTuI6S1NTMtiUWmFmupPbAuXH7BZJ6ADuKtSXgfyX1IWQbtga+FNettjC6Ih5zW0nN47HMjcv/BJwXX2cQLl1mA4WECecTXreQpwrQB/iLmRUCH8aRVnkmEUb3Z1A0laYf4ReUBfFcZBJ+SQAYIekKwi8ExxKC0JcktZc4pmEV7DvhQCLizgUukHRtfN+Y8AsOxFzGFLabZWabASStAL5M+OVhbuJ8ptDOWwfQd1cOL5jOFSPpaMI9s86SjDBaMUmj4ybF55M0M/tM0qmEsNorCakNP0txl3vYf3ukcTnbHQJ8zcx2ltdYLKCTgEmS9hLm4Hym2GYjCSOWrmZWIGlN0r6T80ULCQWpPD8jzOF5auxjcv8+r+Cz5XmKUNyeMLO9Sb8oKC4bk7yxpHbAtUD3+P2YQNHzmTiuQlL/f18XKl94BFxkZnnF+ncGRc9HedsV/x6U199S23Hp5/cwnStpOPAnM/uymbU1s+MJCQeJQOjTJbVTeBjlEuCf8R7kIWb2DGGC+NPiCOEz7Q+S/g4wh5LWEEZMECZwT9gKHJn0fgbwk8SbOKIrQlIvxfukkg4jjLD+XUpbzYBPYrH8OmEEUyYzywfyJZ0ZFyU/3NIM+MhCZuR3CL9glGYucInC/d5jga9XsM9/EwIAflts1SxguKQvxuM8WtKXgaMIBWmzpC+xfwR8QCSdAtwAPFTJj74A/ESxwkvqUsXtEl4D+sRfDBK/2B1IO+4AecF0rqRvEfI2kz3D/qdlFxCSDd4iFNLJhEuaL0vKJSQyJEY//wmMk7SEcD/s1lL2dwtwv6SFhNFEwlRgaHywpTfhqc1uCg8QrSCMZIv7CjBH0lJgEbAQeMZCTud8hYeSxgE5sa2lwHeBlRWfFi4DHorHmHyp8rfAfyo8DHQiZY8qJxPiqFYAfyTcZy2Xmf3ekjJP47IVhF9KZsTz+iIhQmpxPOaVhPu981M4puJ6xwdn8giF8mozm1XJNn5NuEy9RNLy+L4q2wFgZhuAKwhXDhYTRuCVbscdOE8rcc4551LgI0znnHMuBV4wnXPOuRR4wXTOOedS4AXTOeecS4EXTOeccy4FXjCdc865FHjBdM4551LgBdM555xLwf8DX0f4lRmA2dUAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "from causallib.evaluation import evaluate\n", "import matplotlib.pyplot as plt\n", "\n", "evaluation_results = evaluate(ipw, X, a, y)\n", "\n", "f, ax = plt.subplots(figsize=(6, 6))\n", "fig = evaluation_results.plot_covariate_balance(kind=\"love\", phase=\"train\", ax=ax)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that prior to applying IP-weighting (orange triangles), the absolute mean-difference across treatment groups was up to 0.2 standard-deviations. \n", "This can easily bias a simple marginal comparison (i.e., simply check the outcome of treated against the outcome of untreated), since we see the groups are not marginally similar. \n", "For example, one could argue that difference in outcome is contributed to the fact that the two treatment groups differ on their average age. \n", "Since the treatment groups have different distribution of characteristics, it should be accounted for in a causal inference analysis, for-example: the importance-sampling-like scheme applied by the IPW model. \n", "After balancing (blue dots), we see that the marginal difference between the groups' covariate distribution is much smaller on average, suggesting that the IPW model had successfully balanced the dataset." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.8.10 ('causal')", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" }, "vscode": { "interpreter": { "hash": "69f88c614e1e3c90b1e7d1bcdbafb917f9042edfca4477b05464240f7cf41bf5" } } }, "nbformat": 4, "nbformat_minor": 4 }