{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# LaLonde Dataset\n", "\n", "Economists have long-hypothesized that training programs could improve the labor market prospects of participants.\n", "In an attempt to test (or demonstrate) this, the National Supported Work (NSW) Demonstration was initiated using combined private and federal funding. This program was implmented during 1975 to 1979 in 15 locations across the US.\n", "The program provided 6-18 month training for individuals who had faced economic and social problems (such as women receiving Aid to Families with Dependent Children, former drug addicts, ex-convicts, and former juvenile delinquents).\n", "\n", "Participants were randomly assigned into experimental group (Support Work Programs) and control groups.\n", "However, due to the long duration of the study, participants joining the program at the beginning had different charachteristics than people joining later. \n", "Therefore, his covariate shift should be adjusted for in order to estimate the true causal effect of the job-program on future employment.\n", "\n", "\n", "This dataset had become a common benchamrk for causal analysis over the years. \n", "Original analysis of the study was done by [Robert LaLonde](https://en.wikipedia.org/wiki/Robert_LaLonde) and published in his 1986 [Evaluating the Econometric Evaluations of Training Programs with Experimental Data](http://people.hbs.edu/nashraf/LaLonde_1986.pdf) paper. \n", "However, we will follow a later propensity-based analysis made by Dehejia and Wahba in their 1999 [Causal Effects in Non-Experimental Studies: Reevaluating the Evaluation of Training Programs](https://users.nber.org/~rdehejia/papers/dehejia_wahba_jasa.pdf)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Data\n", "First, let's download the dataset from [Rajeev Dehejia's webpage](https://users.nber.org/~rdehejia/nswdata2.html)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(445, 10)\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", "
trainingageeducationblackhispanicmarriedno_degreere74re75re78
2840.028.012.01.00.01.00.00.0000.00000.0000
3580.019.011.01.00.00.01.01626.6230.00000.0000
1171.020.012.01.00.00.00.00.000377.56861652.6370
4120.025.011.01.00.00.01.015209.9903072.7260284.6584
701.027.09.01.00.00.01.00.0000.00000.0000
\n", "
" ], "text/plain": [ " training age education black hispanic married no_degree \\\n", "284 0.0 28.0 12.0 1.0 0.0 1.0 0.0 \n", "358 0.0 19.0 11.0 1.0 0.0 0.0 1.0 \n", "117 1.0 20.0 12.0 1.0 0.0 0.0 0.0 \n", "412 0.0 25.0 11.0 1.0 0.0 0.0 1.0 \n", "70 1.0 27.0 9.0 1.0 0.0 0.0 1.0 \n", "\n", " re74 re75 re78 \n", "284 0.000 0.0000 0.0000 \n", "358 1626.623 0.0000 0.0000 \n", "117 0.000 377.5686 1652.6370 \n", "412 15209.990 3072.7260 284.6584 \n", "70 0.000 0.0000 0.0000 " ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "\n", "columns = [\"training\", # Treatment assignment indicator\n", " \"age\", # Age of participant\n", " \"education\", # Years of education\n", " \"black\", # Indicate whether individual is black\n", " \"hispanic\", # Indicate whether individual is hispanic\n", " \"married\", # Indicate whether individual is married\n", " \"no_degree\", # Indicate if individual has no high-school diploma\n", " \"re74\", # Real earnings in 1974, prior to study participation\n", " \"re75\", # Real earnings in 1975, prior to study participation\n", " \"re78\"] # Real earnings in 1978, after study end\n", "\n", "treated = pd.read_csv(\"http://www.nber.org/~rdehejia/data/nswre74_treated.txt\", \n", " delim_whitespace=True, header=None, names=columns)\n", "control = pd.read_csv(\"http://www.nber.org/~rdehejia/data/nswre74_control.txt\",\n", " delim_whitespace=True, header=None, names=columns)\n", "lalonde = pd.concat([treated, control], ignore_index=True)\n", "lalonde = lalonde.sample(frac=1.0, random_state=42) # Shuffle\n", "\n", "print(lalonde.shape)\n", "lalonde.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see the dataset has 445 inidividuals in it, as described on Dehejia's page. \n", "\n", "### Design matrix\n", "#### Earning indications\n", "Following the analysis performed by Gelman et al. on their [arm](https://cran.r-project.org/web/packages/arm/index.html) R library, we will create two indicator variables indicating no earnings in 1974 and 1975." ] }, { "cell_type": "code", "execution_count": 2, "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", "
trainingageeducationblackhispanicmarriedno_degreere74re75re78re74=0re75=0
2840.028.012.01.00.01.00.00.0000.00000.000011
3580.019.011.01.00.00.01.01626.6230.00000.000001
1171.020.012.01.00.00.00.00.000377.56861652.637010
4120.025.011.01.00.00.01.015209.9903072.7260284.658400
701.027.09.01.00.00.01.00.0000.00000.000011
\n", "
" ], "text/plain": [ " training age education black hispanic married no_degree \\\n", "284 0.0 28.0 12.0 1.0 0.0 1.0 0.0 \n", "358 0.0 19.0 11.0 1.0 0.0 0.0 1.0 \n", "117 1.0 20.0 12.0 1.0 0.0 0.0 0.0 \n", "412 0.0 25.0 11.0 1.0 0.0 0.0 1.0 \n", "70 1.0 27.0 9.0 1.0 0.0 0.0 1.0 \n", "\n", " re74 re75 re78 re74=0 re75=0 \n", "284 0.000 0.0000 0.0000 1 1 \n", "358 1626.623 0.0000 0.0000 0 1 \n", "117 0.000 377.5686 1652.6370 1 0 \n", "412 15209.990 3072.7260 284.6584 0 0 \n", "70 0.000 0.0000 0.0000 1 1 " ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lalonde = lalonde.join((lalonde[[\"re74\", \"re75\"]] == 0).astype(int), rsuffix=(\"=0\"))\n", "lalonde.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Factorizing education\n", "Since years of schooling are not to be taken by their numerical values, we will factorize it into indicator variables." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(445, 24)\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", "
trainingageblackhispanicmarriedno_degreere74re75re78re74=0...education_7.0education_8.0education_9.0education_10.0education_11.0education_12.0education_13.0education_14.0education_15.0education_16.0
2840.028.01.00.01.00.00.0000.00000.00001...0000010000
3580.019.01.00.00.01.01626.6230.00000.00000...0000100000
1171.020.01.00.00.00.00.000377.56861652.63701...0000010000
4120.025.01.00.00.01.015209.9903072.7260284.65840...0000100000
701.027.01.00.00.01.00.0000.00000.00001...0010000000
\n", "

5 rows × 24 columns

\n", "
" ], "text/plain": [ " training age black hispanic married no_degree re74 \\\n", "284 0.0 28.0 1.0 0.0 1.0 0.0 0.000 \n", "358 0.0 19.0 1.0 0.0 0.0 1.0 1626.623 \n", "117 1.0 20.0 1.0 0.0 0.0 0.0 0.000 \n", "412 0.0 25.0 1.0 0.0 0.0 1.0 15209.990 \n", "70 1.0 27.0 1.0 0.0 0.0 1.0 0.000 \n", "\n", " re75 re78 re74=0 ... education_7.0 education_8.0 \\\n", "284 0.0000 0.0000 1 ... 0 0 \n", "358 0.0000 0.0000 0 ... 0 0 \n", "117 377.5686 1652.6370 1 ... 0 0 \n", "412 3072.7260 284.6584 0 ... 0 0 \n", "70 0.0000 0.0000 1 ... 0 0 \n", "\n", " education_9.0 education_10.0 education_11.0 education_12.0 \\\n", "284 0 0 0 1 \n", "358 0 0 1 0 \n", "117 0 0 0 1 \n", "412 0 0 1 0 \n", "70 1 0 0 0 \n", "\n", " education_13.0 education_14.0 education_15.0 education_16.0 \n", "284 0 0 0 0 \n", "358 0 0 0 0 \n", "117 0 0 0 0 \n", "412 0 0 0 0 \n", "70 0 0 0 0 \n", "\n", "[5 rows x 24 columns]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lalonde = pd.get_dummies(lalonde, columns=[\"education\"], drop_first=True)\n", "print(lalonde.shape)\n", "lalonde.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Variables selection\n", "Lastly, we extract the covariates, treatment and outcome variables" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "((445, 22), (445,), (445,))" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = lalonde.pop(\"training\")\n", "y = lalonde.pop(\"re78\")\n", "X = lalonde\n", "X.shape, a.shape, y.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Model\n", "After defining the design matrix `X`, we can continue to define the causal model.\n", "\n", "In the spirit of Dehejia and Wahba's propensity-based analysis, \n", "we will use an Inverse Treatment Probability Weighting (IPTW, or IPW) causal model. \n", "Briefly, this model will model the probability of participants to be assigned to job-training program and use it to emulate two equal-sized populations: one of those being assigned to the program and another of ones who don't.\n", "In this synthetic population, we could use the actual earnings in 1978 to estimate what would have happen if everyone were to join the program or everyone were to not be part of it at all.\n", "\n", "Before we define the causal model itself, we will need to use a machine learning model to estimate the propensity score $\\Pr[A=1|X]$ - the probability of each participant to be assigned to job training. \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": 5, "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": 6, "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 job-training on yearly earnings.\n", "\n", "First, we will fit our causal model. \n", "Second, we'll predict the potential outcomes: what would be the earnings in 1978 if everyone were to join the job-training or everyone were to not join the training program.\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": 7, "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 earnings at the end of the study (1978) if everyone would have join the program (`1`) is 6141\\$ while the averge earnings if everyone would have not join the program (`0`) is 4598\\$." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0 4598.414070\n", "1.0 6141.051856\n", "dtype: float64" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "outcomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Therefore we can conclude that the average additive effect (`diff`) of job training on income is 1543\\$." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "diff 1542.637786\n", "dtype: float64" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "effect" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Indeed, we can conclude that job-training accounts for an income gain of 1543\\$ 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": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "diff 1794.342404\n", "dtype: float64" ] }, "execution_count": 10, "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 overestimate the effect of job training on income by 252\\$. \n", "This seemingly negligible bias accounts for \\~16\\% of our estimation, and probably due to the study being randomized in nature.\n", "\n", "To see why the differences between estimation methods, 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": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAGeCAYAAAAHVAHJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xm8lGX9//HXW0GBUCioSKTAUkwxj4qlUmhiiaWG5JqmmGWWS8s3k/q2aMs3MzNtTzTRX5aiggspLrjgLqBHxAX3ApUUDQRRNj+/P+5rYDhn5pw5hzln5sy8n4/HeTBzL9d9zaBc57qXz1sRgZmZWb3aqNIdMDMzqyQPhGZmVtc8EJqZWV3zQGhmZnXNA6GZmdU1D4RmZlbXPBCalZmk2yV9uYPa/r6kCzqi7Woi6UhJN3VQ289L2qcj2rauyQOh1a30D+Kbkpbl/fy+0v3KkbSXpAX5yyLi/yKiQwbZcpE0UdLPNqSNiLg0Ij5drj6ZtaRbpTtgVmEHRMQtle5EPZHULSJWV7ofZjmeEZo1IWlTSYslDctb9u40e3yPpHdKmirpFUn/Ta+3LNLW6ZL+lvd+sKSQ1C29P1bS45KWSnpW0lfT8ncANwBb5M1WtyjQ3oGSHk39vV3Sh/PWPS/pO5LmSFoi6XJJPYr0cyNJP5D0L0kvS7pEUp8mfT5G0r8lLZL0v0XaOR44Evhu6vN1eX05TdIc4A1J3SSNl/RM+uyPSToor51xku7Kex+STpD0VPrO/yBJLXznV6bPu1TSg5J2bLJZQ6HvpbW/29SvZ1O7z0k6Mm/dl9Lf5X8l3SjpA4X6Z9XHA6FZExGxApgMHJG3+FDgjoh4mez/m4uADwDvB94E2ntK9WVgf2Bz4FjgN5J2jog3gP2AFyOid/p5MX9HSdsA/wC+CbwbuB64TtImTfo9GhgCfAQYV6Qf49LPJ4GtgN4FPtPHgaHAKOBH+YNuTkScD1wKnJX6fEDe6iOAzwJ904zwGeATQB/gDOBvkt5XpH+QfU+7Ajumz7VvC9t+DrgCeBfwd+BqSd3z1hf7Xor+3aZfTn4L7BcRmwF7AI1p3Rjg+8BYsr+LO8n+bqwL8EBo9e7qNJvK/XwlLf876w+EX0jLiIhXI+KqiFgeEUuBnwN7tufgEfHPiHgmMncAN5ENDqU4DPhnRNwcEauAs4GeZP9A5/w2Il6MiNeA64CGIm0dCZwTEc9GxDLge8DhuZlrckZEvBkRDwMPkw1IbfHbiJgfEW8CRMQVqW9vR8TlwFPAR1vY/8yIWBwR/wZua+GzAMyOiCvT93IO0APYrUlfmn0vJfzdvg0Mk9QzIl6KiEfT8q8Cv4iIx9Mg/39ks07PCrsAD4RW78ZERN+8nwlp+a1AT0kfS/+YNQBTACT1kvSXdBrxdWAG0FfSxm09uKT9JN0n6TVJi4HPAP1L3H0L4F+5NxHxNjAfGJi3zcK818vJZnqttpVedwPe2462ipmf/0bS0ZIac7+EAMNo+bO35fhrj5W+lwVkn7HFtlr6u02z9MOAE4CXJP1T0rapjQ8A5+V9ltcAsf7fhVUpD4RmBaR/PCeRzQq/AExNMwSA/yE7RfixiNgcGJmWF7pm9QbQK+/9gNwLSZsCV5HN5N4bEX3JTm/m2mktGuZFsn+Ac+0JGAS80Nrna60tstOCq4H/tKOtYv1euzz9cjEBOAnolz77XAp/h+0xKO9YGwFbkn3G1rT4dxsRN0bEp4D3AU+kzwDZwPvVJr9U9YyIe8rzcawjeSA0K+7vZDOAI9PrnM3Irh0tlvQu4McttNEIjJT0/nTzyffy1m0CbAq8AqyWtB+Q/8jAf4B+uZtWCpgEfFbSqHT963+AFUB7/vH9B/AtSUMk9SY7tXd5O+/u/A/ZdcaWvINsYHwFspuGyGaE5bKLpLHp1O43yb6X+0rYr+jfraT3Krs56R2pvWXAmrT6z8D3JG2ftu0j6ZDyfRzrSB4Ird5dp/WfI5ySWxER95PN6LYgu4Mz51yya3GLyP5xnVas8Yi4GbgcmAPMBqbmrVsKnEI2oP2XbOZ5bd76J8gGqGfTKbf8U3tExDzgKOB3qS8HkD0OsrKtXwLwV+D/kZ0KfA54Czi5He0AXAhsl/p8daENIuIx4NfAvWQD5w7A3e08XiHXkP0S81/gi8DYdL2wNS393W5E9svGi2SnPvcEvg4QEVOAXwKXpVOqc8ludrIuQA7mNbNaIul04EMRcVSl+2Jdg2eEZmZW1zwQmplZXfOpUTMzq2ueEZqZWV3zQGhmZnXNA6GZmdU1D4RmZlbXPBCamVld80BoZmZ1zQOhmZnVNQ+EZmZW1zwQmplZXfNAaGZmdc0DoZmZ1TUPhGZmVtc8EJqZWV3zQGhmZnXNA6GZmdU1D4RmZlbXPBCamVld80BoZmZ1zQOhmZnVNQ+EZmZW1zwQmplZXfNAaGZmdc0DoZmZ1TUPhGZmVtc8EJqZWV3zQGhmZnXNA6GZmdU1D4RmZlbXPBCamVld61bpDljr+vfvH4MHD650N8zMupTZs2cvioh3t7adB8IuYPDgwcyaNavS3TAz61Ik/auU7Xxq1MzM6poHQjMzq2seCM3MrK55IDQzs7rmgdDMzOqaB0IzM6trHgjNzKyueSA0M7O65oHQzMzqmgdCMzOrax4IzcysrnkgTCQNljS30v3oSFc/9AIjzryVIeP/yYgzb+Xqh16odJfMzApbuhDO2xGW/qfDD+WBsANI2rjSfWjq6ode4HuTH+GFxW8SwAuL3+R7kx/xYGhm1emOs2Dxv+GOX3b4obps+oSkwcANwF3AHsALwOeAocCfgV7AM8CXIuK/RdrYBfgrsDy1k1u+MXAmsBewKfCHiPiLpI2A3wN7As+R/SLx14i4UtLzqa1PA7+XNBP4A/Du1P5XIuIJSe9O/Xt/Otw3I+LuDf9GWvarG+fx5qo1AIze5Im1y2+77gkWP9yX7bffnl133ZVVq1Zx6aWXNtu/oaGBhoYGli9fzqRJk5qtHz58OMOGDWPJkiVMmTKl2frdd9+doUOHsmjRIqZOndps/ciRI9lqq61YuHAh06ZNa7Z+1KhRDBo0iPnz5zN9+vRm60ePHs2AAQN49tlnmTFjRrP1+++/P/3792fevHnce++9zdYfdNBB9OnTh7lz5xZM+jj00EPp1asXjY2NNDY2Nlt/5JFH0r17d2bOnMmjjz7abP24ceMAuOeee3jyySfXW9e9e3eOPPJIAO644w6ee+659db36tWLQw89FIBbbrmFBQsWrLd+8803Z+zYsQBMmzaNhQsXrre+X79+HHDAAQBcd911vPrqq+utHzBgAKNHjwZg8uTJvP766+ut33LLLdlnn30AmDRpEsuXL19v/ZAhQ9hzzz0BuPTSS1m1atV667fZZhv22GMPACZOnEhT/m/P/+1Bk//21qyEBcsYF29D46Ww52mw2XubfbZy6eozwq3JBqntgcXA54FLgNMi4iPAI8CPW9j/IuCUiNi9yfLjgCURsSuwK/AVSUOAscBgYAfgy0DT/d6KiI9HxGXA+cDJEbEL8B3gj2mb84DfpLY/D1xQqGOSjpc0S9KsV155pbXvoVUvLn6z4PKVq9dscNtmZmW1eP661/F2h88KFREdeoCOkmaEN0fE1un9aUAP4LiIeH9a9kHgiojYucD+fYBH8rb9CPD3iBgm6UrgI2QzOYA+wFeBzwAPR8RFaZ/JaZ/cjHDPiPiXpN7AK8C8vENuGhEflvQy8GLe8ncD20bE0mKfdfjw4bGheYQjzryVFwoMhgP79uTu8XtvUNtmZmWTuza4+q11y7r1gG/MafOsUNLsiBje2nZdfUa4Iu/1GqBvG/YVUOy3AJHN5hrSz5CIuCktb8kb6c+NgMV5+zdExIfz1u2et3xgS4NguZy671B6dl//0mXP7htz6r5DO/rQZmalu+OsbBaYr4NnhV19IGxqCfBfSZ9I778I3FFow4hYDCyR9PG06Mi81TcCX5PUHUDSNpLeQXYd8fOSNpL0XrJriIXafh14TtIhaX9J2jGtvgk4KbetpIa2f8y2G7PTQH4xdgcG9u2JyGaCvxi7A2N2GtgZhzczK82867NrhPnWrMyWd5Aue7NMC44B/iypF/AscGwL2x4L/FXScrLBL+cCsmuBD0oS2WnOMcBVwChgLvAkcD/Z4FvIkcCfJP0A6A5cBjwMnAL8QdIcsu9/BnBC2z9m243ZaaAHPjOrbv/zROvblFmXvUZYKZJ6R8QySf2AB4AREbGwtf02RDmuEZqZ1ZtSrxHW4oywo02V1BfYBPhpRw+CZmbWsepiIJT0B2BEk8Xn5e7+bIuI2KssnTIzs6pQFwNhRJxY6T6YmVl1qrW7Rs3MzNrEA6GZmdW1ig6EksZJ+n2Z2xwjabu89z+RtE8Z299W0r2SVkj6TpN1oyXNk/S0pPFF9t9U0uVpm/tThRyrIKdyWE3qxPSGrq4WZ4RjgLUDYUT8KCJuKWP7r5E9C3h2/sJUqPsPwH7p+EfkD8h5jgP+GxEfAn4DdHxpdSvKqRxWszoxvaGr69CbZSQdRTZobEL28PnXgaOB7wEvkT2UviJtOxGYGhFXpvfLIqJ3ev1dsioxbwM3RMR4SV8Bjk9tP53WNwAHAnumB9k/D/ww166kUWQDWDdgJvC1iFiR6oReDBxA9vD7IRFR8KnOiHgZeFnSZ5us+ijwdEQ8m/p8GVkaxmNNtvsccHp6fSVZUoXCD3RWhFM5nIwANZjKsWYlw1+4n2GdlN7Q1XXYjFDSh4HDyB44byCrBXoUcAbZowyfIm/m1kI7+5HN8j4WETsCZ6VVkyNi17TscbJi2/cA1wKnpjqez+S10wOYCBwWETuQDYZfyzvUolSc+09kaRFtNRDIK5nOgrSs6HYRsZqsMk2/phuVO33CCnMqh9WkxfPJ5g10SnpDV9dhlWUknQR8H3g5LeoJvAnMiYij0zanANtExEnFZoSSfg08ERETmrS/J/AzskLbvYEbI+KEAu1MBKYCTwG/i4iRafko4MSIGJtmhCMi4gVJHwN+HhEtXleUdDqwLCLOTu8PAfaNiC+n918EPhoRJzfZ79G03YL0/pm03fq/juZxZZmO41QOqzllTG/o6qohfULAxXkpC0PJTgkWG3lX5/qT6ntuktdOoX0mAiel2d0ZZBFMrfWnJbkkizW075TxAmBQ3vstWT9uqdl2krqRRTy91o7jWRk4lcNqTgXSG7q6jhwIpwMHS3oPgKR3AQ8Be0nql5IdDsnb/nlgl/T6c2TX6iBLa/hSKqKdawdgM+Cl1E5+csTStK6pJ4DBkj6U3hdNpminmcDWkoZI2gQ4nOw0bVPXkhUGBzgYuNXXByvHqRxWcyqQ3tDVddjNMhHxWLph5SZJGwGrgBPJZoX3kt0s8yCQ+3V8AnCNpAfIBtE3UjvTUlTRLEkrgevJTrn+kOwGnH+RJdHnBr/LgAnptOvBef15S9KxwBVpJjYT+HNbP5ekAcAsYHPgbUnfBLaLiNfT6eAb02f6a0Q8mvb5CTArIq4FLgT+n6SnyWaCh7e1D1ZeTuWwmlKB9IauzukTXYCvEZqZtV01XCM0MzOrenVRdLs90mnUbzRZfLcLeJuZ1RYPhEWkiKY2xzSZmVnX4lOjZmZW1zwQmplZXXP6RNvbL5g+IWmQpNskPS7pUUlNry/mtpOk36b0iTmSdi5X36z9nEBRY5y8YG1QizPCiqRPkFXG+Z+I+DCwG3BikfSJ/YCt08/xZLVNrYKcQFGDnLxgbeD0iTKlT0TES+kzERFLJT1OVmC7UPrEJamazH2S+kp6X9rfKsAJFDWWQLFmJSxYxjgnL1iJnD6xzoamT+T3eTCwE9ng31RJKRVOn+g8TqCoMYvz/vdyjU0rgdMnypQ+kbe8N1kN059HxOQC+/0T+EVE3JXeTwe+GxGzix3LlWU6lhMoaoiTFyxPNVSWqbf0CVIB8KuASwsNgkmpKRXWSZxAUUOcvGDt4PSJMkmD94XA4xFxTgubXgscne4e3Q1Y4uuDleUEihri5AVrB6dPtFGx9AngI2SD6yOScncMfD8irpd0QurDn1P/P0N2g89y4Ni29sHKzwkUNcLJC9YOTp/oAnyN0Mys7arhGqGZmVnVc9HtIpw+YWZWHzwQFuH0CTOz+uBTo2ZmVtc8ELZA0p2SGtPPi5KuTsv3krQkb92P2tDmuyTdLOmp9Oc7O+4TmJlZa+p+IEzP8xX8HiLiE7mCAGSPfOQ/JH9nXrGAn7ThkOOB6RGxNdljIuPb3fkuwKkOtkGcImGdoC4HQkmDU1zSH8meZfxiilZ6UNIVqUxa/vabAXsDV5fh8J8jK/BN+nNMGdqsSk51sA3mFAnrBPV8s8xQsofZf0Q209snIt6QdBrwbSB/lncQ2Szu9bxlu0t6mKw82nci4lHITqdSuLLNd1Ic1HtzlWQi4qVc5Z1a5FSHGkt1yLP55pszduxYAKZNm8bChQvXW9+vXz8OOOAAAK677jpeffXV9dYPGDCA0aNHAzB58mRef/319dZvueWW7POxYVl6hFMkrIPV80D4r4i4T9L+ZJVh7s6qpLEJ2WnQfEcAF+S9fxD4QEQsk/QZspni1pCdTi1H5yQdTxYzxfvf//5yNNnpnOpgGyS/bmiuXuj+LVUvNGufuqwsk2KSpkbEMEkHAF+IiCOKbNuPLDdxYES8VWSb54HhEbGotRmhpHnAXmk2+D7g9lSQvKiuWlnGqQ7Wbk6RsDJwZZnS3QeMyBXjltRL0jZ56w8hGzTX/h8paUAqso2kj5J9j6/C+jfYNPm5Je1+LXBMen0McE2HfroKcqqDtZtTJKwT1f1AGBGvAOOAf0iaQzYwbpu3yeHAP5rsdjAwN10j/C1weJQ+tT4T+JSkp8jCic/cgO5XNac6WLs5RcI6UV2eGu1quuqpUTOzSvKpUTMzsxJ4IDQzs7rmgdDMzOqaB0IzM6trHgjNzKyu1cVAmGqLzi2w/CeS9umkPlwgabvOOJaZmZWunkusERElxyeV4Vhf7qxjmbXV1Q+9wK9unMeLi99ki749OXXfoZV93nPpQvjrvvClm1xJxjpcXcwIk40lTZD0qKSbJPWUNFHSwQCSzpT0mKQ5ks5OyyZK+nPKJXwy1SXNzTDvTGkVD0raIy3fS9Ltkq6U9ISkS/Mq0NwuaXh6PTrt97Ck5tWczTpRVaaEOHXCOlE9zQi3Bo6IiK9ImgR8PrdC0rvIEia2jYiQ1Ddvv8HAnsAHgdtSKbaXgU9FxFuStiarPJN7aHMnYHuyVIq7gRHAXXnHejcwARgZEc+lY5tVTNWlhKxZCQuWMc6pE9ZJ6mlG+FxE5LJsZpMNcDmvA28BF0gaCyzPWzcpIt6OiKeAZ8nKr3UHJkh6BLiCLL0i54GIWBARbwONTY4DsBswIyKeA4iI1wp1VtLxkmZJmvXKK6+0/dOalajqUkIWz1/32vVFrRPU04xwRd7rNUDP3JuIWJ2KZ48iqy16ElkQL0DTGnQBfAv4D7Aj2S8T+akUTY/T9DtWgTabiYjzgfMhK7HW2vZm7bVF355rU0KmrVxXZndg3578ady6lJDu3buvzTgspFevXi2u79OnT4vr+/fvz7jPj85SJ3I17tes9KzQOlw9zQiLSon0fSLieuCbQEPe6kMkbSTpg8BWwDygD/BSmvV9Edi4aZstuBfYU9KQdGyfGrWKqqqUEKdOWAV4IMxsBkxN6RN3kM34cualZTcAJ6Q4pj8Cx0i6D9gGeKPUA6W0i+OBySm94vLyfASz9qmqlBCnTlgFOH2iBZImkmURXlnJfjh9wsys7Zw+YWZmVoJ6ulmmzSJiXKX7YGZmHcszQjMzq2seCM3MrK55IDQzs7rmgdDMzOpaVQ6EksZJ+n2Z2xyTH4PUURFMknaVtCZXzLvA+l0kPSLpaUm/zRXlNqu0qx96gRFn3sqQ8f9kxJm3bljR7aULswoxS/9Tvg6adZCqHAg7yBjyaoJGxI8i4pZyHkDSxsAvgRtb2OxPZA/Ub51+RpezD2btUfYECqdHWBdSkccnJB0FnAJsAtwPfB04Gvge8BLwJKlmZ9OH2iUti4je6fV3yUqcvQ3cEBHjJX2FbKDZBHg6rW8ADiQrbfYDsuSJH+balTQKOJvs+5gJfC0iVkh6HrgYOICs0PYhEbGuPH9zJwNXAbsW+dzvAzaPiHvT+0vIBugbSv7yzDpAWRMo/vE3p0dYl9LpM0JJHwYOA0ZERANZYeqjgDPIIos+xfppDsXa2Y9sEPlYROwInJVWTY6IXdOyx4HjIuIe4Frg1IhoiIhn8trpAUwEDouIHcgGw6/lHWpRROxMNpP7Tgv9GUgW5fTnFro9EFiQ935BWlaoPadPWKcpawKF0yOsi6nEjHAUsAswM10e6wnsAdye6nAi6XKyGp4t2Qe4KCKWw3pxRsMk/QzoC/Sm5dOUAEPJIpqeTO8vBk4Ezk3vJ6c/ZwNjW2jnXOC0iFjTwmW/QisK1rhz+oR1prIlUKx5nXFvnu/0COtSKnGNUMDFaWbWEBFDgdMpHk20mtTPdGPJJnntFNpnInBSmt2dAfQooT8tycUqFYpUyjccuCydTj0Y+KOkMU22WQBsmfd+S7IAX7OKKlsChdMjrAuqxEA4HThY0ntgbQzRQ8BekvpJ6g4ckrf982QzSIDPkV2rA7gJ+JKkXnntQJYk8VJq58i8dpamdU09AQxOyfOQXVO8o60fKiKGRMTgiBgMXAl8PSKubrLNS8BSSbulQf1o4Jq2Hsus3MqWQOH0COuCOv3UaEQ8lm5YuUnSRsAqslORp5Nl9b0EPMi6jL8JwDWSHiAbRN9I7UyT1ADMkrQSuB74PtlNMPcD/wIeYd3gdxlZqvwpZDO2XH/eknQscIWk3M0yLV3nazNJjel6KGTXHyeSnRK+Ad8oY1VizE4DNzx66X9aupfMrDo5hqkLcAyTmVnbOYbJzMysBI5haqN0GvUbTRbfHREnVqI/Zma2YTwQtlFEXARcVOl+mJlZefjUqJmZ1TUPhGZmVtcqOhB2xZQJSdtKulfSCknNSq5J2ljSQ5KmFtl/U0mXp/SJ+yUNLlffrLLKmt5QrZwqYTWoFmeEHZ0y8RpZwfCzi6z/BlmN02KOA/4bER8CfkOWVmFdXNnTG6qVUyWsBnXozTK1mDIRES8DL0v6bIHPuyXwWeDnwLeLfC2fIyseAFkFmt9LUviBzi6trOkNkyY1Wz98+HCGDRvGkiVLmDJlSrP1u+++O0OHDmXRokVMndr8ZMTIkSPZaqutWLhwIdOmTWu2ftSoUQwaNIj58+czffr0ZutHjx7NgHfAsw/dwVZOlbAa02EzwlpNmWjFucB3yQbsYgYC8wEiYjWwBOjXdCOnT3QtZU1vqFb5dURdP9RqSEfOCGs1ZaIgSfsDL0fEbEl7tbRpgWXNZoNOn+haypbe0KtXi+v79OnT4vr+/fu3uH7AgAEtrh80aFDh9UsXQuOlbPW2UyWs9nTkNcJaTZkoZgRwYDrNehmwt6S/FdhuATAIINU27UN23dG6sLKlN1Qrp0pYDevIgbAmUyaKiYjvRcSWKX3icODWiDiqwKbXAsek1wen7Tzj6+LKlt5QrZwqYTWsw06N1mrKhKQBwCxgc+BtSd8EtouI11vY5yfArIi4FrgQ+H+SniabCR7e1j5YdSpLekO1cqqE1TCnT3QBTp8wM2s7p0+YmZmVwEW3i3DKhJlZffBAWIRTJszM6oNPjZqZWV3zQGhmZnXN6RNtb79g+oSkHpIekPSwpEclnVFkf6dPtFFdpDpUmlMlrI7V4oywUukTK4C9U+3TBmC0pN0K7O/0iTaom1SHSnOqhNUxp0+UKX0iVYdZlt52Tz+FHtJ0+kQb1EWqw4ABPPvss8yYMaPZ+v3335/+/fszb9487r333mbrDzroIPr06cPcuXMp9KzpoYceSq9evWhsbKSxsbHZ+iOPPJLub73KzAcb2dWpElannD6xzganT6RQ3kbgZeDmiLi/wGZOn2iDukh1qLQ7zmLt72yuH2p1qMMqy0g6iawU2stpUU/gTWBORBydtjkF2CYiTio2I5T0a+CJiJjQpP09gfXSJyLihALtTASmAk8Bv4uIkWn5KODEiBibZoQjIuIFSR8Dfh4RLV5XlHQ6sCwimgX0SuoLTAFOjoi5TdY9CuwbEQvS+2eAj0bEq8WOVc+VZUaceevaVId8A/v25O7xexfYw9okd21w9VvrlnXrAd+Y41mhdXnVUFmm3tIn1oqIxcDtwOgCq50+0QY1n+pQaU6VMHP6RPs/3vokvTvNBJHUkyxHsdB1RqdPtEHNpzpUmlMlzJw+0dbPVSx9AngfcLGkjcl+wZgUEVPTPk6f2AA1nepQaU6VMHP6RFdQz9cIzczaqxquEZqZmVU9F90uwukTZmb1wQNhEU6fMDOrDz41amZmdc0DoZmZ1bWqHAi7aCrFqZIa089cSWvynnnM325ISp14KqVQbFKoPausqk+8cFqEWdlU5UDYQTo0lSIifpWrokNWVPyOiChUMeaXwG8iYmvgv2RpFFZFukTihdMizMqmIjfL1GIqRRNHAP8o8LkF7A18IS26mKzAwJ9KaNM6SdUnXgwfxlaNl2al0JwWYbbBOn1GWOupFKkU3GjgqgKr+wGLU+oEZHVHC5ZMcfpE5VR94sWcy9fVB3VdULMNVokZ4SiymqIzswkSPYE9gNsj4hUASZcD27TSzj7ARRGxHCDvNOQwSeulUrTSzlDguYh4Mr2/mKwU3Lnp/eT052xgbKufLps93l3ktGihwt8FS/tExPnA+ZBVlinhuFYmW/TtuTbxYtrKbdcuH9i3J38aty7xonv37owbN65oO7169WpxfZ8+fVpc379//+brly6E805bVx90zUrPCs02UCWuEdZ6KsXhFDgtmiwC+qZapwBbAi+W0KZ1oqpOvHBahFnZVWIgrNlUCkl9gD2BawqtTykTt7GuGPgxxba1yqnqxAunRZiVXaefGq3VVIrkIOCmiHgjf6Gk64EvR8SLwGnAZen07UNkaRRWZao28cJpEWZl5/SJLsDpE2Zmbef0CTMzsxK46HYbOZUAym5zAAAgAElEQVTCzKy2eCBsI6dSmJnVFp8aNTOzuuaB0MzM6lpJA6GkQyRtll7/QNJkSTt3bNfMNlzVp0jkc6KEWUWUOiP8YUQslfRxYF+yMmQuFJ1IulrSbEmPSjo+LTtO0pOSbpc0IRcrJendkq6SNDP9jKhs72tXl0iRyOdECbOKKPVmmVy14c8Cf4qIaySd3jFd6pK+FBGvSepJVkP1n2QP9u9MVtHmVuDhtO15ZDFMd0l6P1kt1A9XotO1rupTJEaOZKuttmLhwoVM++e1sGAZ45woYdbpSp0RviDpL8ChwPWSNm3DvvXgFEkPA/cBg0hl2iLitYhYBVyRt+0+wO8lNZIlYmyeO+2cz+kTG67qUyTyLZ6/7rVrh5p1qpIqy+RFCz0SEU9Jeh+wQ0Tc1NEdrHaS9gJ+Bnw6IpZLup1s1jcmIo5J25wCbBMRJ0laBAyKiML/ShfgyjLtM+LMW9emSOQb2Lcnd4/fu8AeFZK7Nrj6rXXLuvWAb8zxrNBsA5S1skyKOnoZ+HhatBp4qv3dqyl9gP+mQXBbYDegF1kI8DtT/dLP521/E3BS7k2ql2odoKpTJPI5UcKsokq9a/THZMWiv5cWdQf+1lGd6mKmAd0kzQF+SnZ69AXg/8iKf98CPAYsSdufAgyXNEfSY8AJnd/l+lDVKRL5nChhVlGlnhptBHYCHoyIndKyORHxkQ7uX5clqXdELEszwinAXyOi+R0VJfCpUTOztit30e2VKUsvUuPv2JDO1YnT0y8Qc4HngKsr3B8zMyug1McnJqW7RvtK+grwJeCCjutW1xcR36l0H8zMrHUlDYQRcbakTwGvA0OBH0XEzR3aMzMzs05Q0kAo6ZcRcRpwc4FlZmZmXVap1wg/VWDZfuXsiJmZWSW0OCOU9DXg68BW6fGAnM2AuzuyY2ZmZp2htRnh34EDyEqBHZD3s0tEHLWhB5c0LleMulwkjZG0Xd77n0jap4ztbyvpXkkrJH2nybq+kq6U9ISkxyXtXmB/SfqtpKfTs4RO8agR7Uq6cOKEWcW1OBBGxJKIeD4ijoiIfwFvkj1C0TsVjK5GY4C1A2FE/Cgibilj+6+RPRR/doF15wHTImJbYEfg8QLb7AdsnX6OxykeNaHdSRdOnDCruFJvljkAOAfYgqzU2gfI/pHfvpX9jiIbNDYhq7LydeBosgo1LwFPAivSthOBqRFxZXq/LCJ6p9ffJStk/TZwQ0SMT49xHJ/afjqtbwAOJCtv9gOy0mY/zLUraRTZANYNmAl8LSJWSHqeLFrqALKqOYdExLq4gjwR8TLwsqTPNvmsmwMjgXFpu5XAymYNwOeAS9JzmfelWeT7IuKllr5Lq27tSrpYs9KJE2ZVoNSbZX5GVkPzyYgYAoyilWuEkj4MHAaMiIgGsiino4AzgBFkN+BsV7yFte3sRzbL+1hE7AiclVZNjohd07LHgeMi4h6y07inRkRDRDyT104PYCJwWETsQDYYfi3vUIsiYmeyGVp7ngHcCngFuEjSQ5IuKFJ4YCCQFzXAgrSs6ed2+kQX0q6kCydOmFWFUh+oXxURr0raSNJGEXGbpNb+rx0F7EKWzwfQE9gDuD0iXgGQdDmwTSvt7ANclAp/ExGvpeXDJP0M6Av0Jsv1a8lQ4LmIeDK9vxg4ETg3vZ+c/pwNjG2lrUK6keUPnhwR90s6DxhPNiPNpwL7NqtzFxHnA+dDVmKtHf2xTrRF355rky6mrdx27fKBfXvyp3Hrki66d+/OuHHj1l0bjJQ4sWalZ4VmFVLqjHCxpN7ADODS9I/86lb2EXBxmpk1RMRQ4HQK/KOfrM71R9nIuUleO4X2mQiclGZ3ZwA9SuhPS1akP9dQ+i8I+RYACyLi/vT+SrKBsdB2g/Lebwm82I7jWRVpc9KFEyfMqkapA+HnyG6U+RZZ2sIzZNfTWjIdOFjSewAkvQt4CNhLUj9J3YFD8rZ/nmwGmTte9/T6JuBLKRMx1w5kj3C8lNo5Mq+dpWldU08AgyV9KL3/InBHK5+hZBGxEJgvKfcv3yiy1ImmrgWOTneP7gYs8fXBrq/NSRdOnDCrGqWWWHsj7+3FJe7zWLph5SZJGwGryE5Fng7cS3azzINA7tfoCcA1kh4gG0TfSO1MS5l9syStBK4Hvk92yvF+4F/AI6wb/C4DJqQw3IPz+vOWpGOBK1IixEzgz6V8lnySBgCzgM2BtyV9E9guIl4HTiabMW8CPAscm/Y5IfXhz6n/nyG7wWd5bhvr+sbsNLD0iKf/KXgvlplVQIsxTJLuioiPS1rK+qcnBUREbN7RHTTHMJmZtUepMUwtzggj4uPpz0KnGs3MzLq8Vq8RpjtF53ZGZ6qJpGMlNTb5+UOl+2VmZuXV6jXCiHhb0sOS3h8R/+6MTlWDiLgIuKjS/TAzs45V6mMC7wMeTTeyrL1xJiIO7JBemZmZdZJSB8IzOrQXZmZmFVLSc4QRcUehnw09eA2mTzwv6ZF0PbHgbZ5On+j62pUyUYzTJ8wqrqSBUNJukmZKWiZppaQ1kl7v6M61UyXTJwA+mSrpFLtl1+kTXVi7UyaKcfqEWcWVemr098DhwBXAcLIEia1b26me0ifawOkTXVi7UibyNDQ00NDQwPLly5n0j785fcKsCpRaYo2IeBrYOCLWpDsq92pp+zpMn4Cs6MBNkmZLOr7INk6f6MLalTJRjNMnzKpCqTPC5alsWKOks8hmc4UihvLVW/oEZIP+i6m+6s2SnoiIGU22cfpEF9bmlIkieq15nXFvnu/0CbMqUOqM8Itp25PIHp8YRHbasSX1lj5BRLyY/nwZmAJ8tMBmTp/owtqcMlGM0yfMqkapA+HOZLVFX4+IMyLi2+lUaUvqKn1C0jskbZZ7DXwaKFSRx+kTXVibUyaKcfqEWdUodeZzIHCupBlk6Q43RkSLeYT1lj4B9AempNPA3YC/R8S0tI/TJ2pIm1IminH6hFnVaDF9Yr0Ns5nXfmQ3wHwcuDkivtyBfbPE6RNmZm1XlvSJfBGxStINZNfrepKdvvRAaGZmXVqpD9SPTs/5PU12uvECsvqjNcvpE2Zm9aHUGeE4smtvX42IFa1sWxOcPmFmVh9KGggj4nBJ7wU+lW4GeSA9ImBmZtallXpq9BDgAbLHHQ4F7pd0cMt7mZmZVb9SnyP8AbBrRBwTEUeTPSj+w47qVFdMpUht7pWuJT4qqeAzipKGSLpf0lOSLk8Ve6xKlS1pwikTZlWr1IFwoyanQl9tw77VokNTKST1Bf4IHBgR27N+sYB8vwR+ExFbA/8FjitXH6y8ypo04ZQJs6pV6s0y0yTdCPwjvT+M7MHwdqnFVArgC2SFwP8Na8usNf3cAvZO25LaPh1HMVWlsiVNvPwck2Y7ZcKsWrU4q5P0IUkjIuJU4C/AR4AdySrDnN+eA9ZwKsU2wDsl3Z7SJ44usE0/YHFeVZ6CyROpX06fqLCyJU3ck3eW3/VEzapOazPCc8nKmRERk0kJDZKGp3UHtOOYtZpK0S19rlHpM90r6b68dqHE5Alw+kQ1KEvSxNKF9Jr7N8Y5ZcKsarV2nW9wRMxpujAiZgGD23nMWk2lWABMi4g3ImIRMINs9pxvEdA31ToFJ09UtbIkTThlwqzqtTYQtjSI9GznMWs1leIa4BOSuqU+fYzs1OxaKZX+NtYVAz8m7WdVqCxJE06ZMKt6rZ0anSnpKxExIX+hpOPIThW2Wa2mUkTE45KmAXPIbt65ICLmAki6Hvhyyis8Dbgsnb59CLiwrceyzrPBSRNOmTCrei2mT6RqMlOAlawb+IaTnZ48KCIWdngPzekTZmbtUJb0iYj4D7CHpE8Cw9Lif0bErWXoo5mZWcWVWmv0NrJrW3UvnUb9RpPFd0fEiZXoj5mZbZiS8wgt41QKM7Pa0tXKpJmZmZWVB0IzM6trHghbIOnOvHT6FyVd3WT9rpLWtCWSqhrSJ8qWqGDFOW3CrMuo+4FQmYLfQ0R8IlcBh+wZx8l5+21MliTRWgm3piqaPlHWRAUrzmkTZl1GXd4sI2kwcAPZnbC7A+dKOgHYFHgGODYiluVtvxlZasSxec2cDFwF7NqG41Y8faJsiQrLlzNp0qRm64cPH86wYcNYsmQJU6ZMabZ+9913Z+jQoSxatIipU6c2Wz9y5Ei22morFi5cyLRp05qtHzVqFIMGDWL+/PlMnz692frRo0czYMAAnn32WWbMmNFs/f7770///v2ZN28e9957b7P1Bx10EH369GHu3LkUenbz0EMPpVevXjQ2NtLY2Nhs/ZFHHkn3t15l5oON7Oq0CbMuoZ5nhEOBS8jSLo4D9kkpE7OAbzfZ9iBgekS8DiBpYFq2XgUaSZvlnUpt+rMdVZA+UbZEBSvujrNYWwbXdUXNql6LlWVqVZoR3hYRQyTtT1aoe0FavQlwb0Qcl7f9DWQl065K768Afh0R9zXNS2zluO9ObX8ovR8EXJ8KhBdVzsoyI868dW2iQr6BfXty9/i9C+xhbZK7Nrj6rXXLuvWAb8zxrNCsk5WlskyNeyP9KeDmiDii0EaS+gEfJZsB5gwnqxcK0B/4jKTVZLVQ7yxyvC+QFeHuK6lbmhV2evrEqfsO5XuTH1l7ehTakahgxbWUNrH/OZXpk5m1qJ4Hwpz7gD9I+lBEPJ2SI7bMyxE8hGzGt/ZX/IgYknudNyPM3VHa0NLBJOXSJy6jAukTuQLSv7pxHi8ufpMt+vbk1H2HblhhaVunpbQJD4RmVanuB8KIeEXSOOAfkjZNi38A5AbCw4Ezy3jIiqdPbHCighXntAmzLqcurxF2NU6fMDNru1KvEdbzXaNmZmYeCM3MrL55IDQzs7rmgdDMzOqaB0IzM6trHgjLTNI9bdz+dEnf6aj+lJNTK0rk5AmzLsUDYTtJ6tbk/cYAEbFHZXrUsZxa0QZOnjDrUurugfpUZ3QacBewG/AwcBFwBvAe4Mi06blAT+BNsjSKeenB+88CPYB3SPoJ8GPgJbKKMttJWhYRvdOxTgUOJUu1mBIRP07L/xc4GpgPvALM7tAPXQZOrSgxtWLWXcyavYxxTp4w6zLqdUb4IeA84CPAtmR1QD8OfAf4PvAEMDIidgJ+BPxf3r67A8dERK5C9UeB/42I7fIPIOnTwNZpfQOwi6SRknYhq1azEzCWIjFOHZU+0V5OrSjRY9ete+3kCbMuoe4qy6QZ4c0pGBdJlwA3RsSlkrYiC989APgt2UAWQPeI2DbNCPeMiGPTvnsBP46IT+a1vywieks6m6ym6OK0qjfwC2Az4F0R8aO0/TnAixFxdrE+V0NlGadWlMDJE2ZVxZVlWrYi7/Xbee/fJjtd/FOymKZhZINij7zt32B9Td/nCPhFLuE+Ij4UEbm6ol3ut49T9x1Kz+4br7fMqRVNtJQ8YWZVq14Hwtb0AXJ3gYxrZxs3Al+SlLteOFDSe4AZwEGSekrajGygrXpjdhrIL8buwMC+PRHZTPAXY3dw8e58LSVPmFnVqrubZUp0FnCxpG8Dt7angYi4SdKHgXtTbuEy4KiIeFDS5UAj8C+K5xdWHadWtMLJE2ZdUt1dI+yKquEaoZlZV+NrhGZmZiXwQGhmZnXNA6GZmdU1D4RmZlbXPBCamVld80DYQSTdKakx/bwo6eq0fC9JS/LW/agS/XOSRDs5WcKs5vg5wg2g7AFBRTQtJwIR8Ym87a4CrslbfWdE7N8JXSwolySRK6KdS5IA/Jxga/KTJfY/p9K9MbMy8EDYRqlW6Q3AbWQFuM+VdAJZwsQzZEkVy/K23wzYGzi20ztbhJMkSkySmDuX9Z7fXLMSFjhZwqzW+NRo+wwFLgE+BRwH7BMROwOzgG832fYgYHpEvJ63bHdJD0u6QdL2hQ7QkekTTpJop8Xz1712DVGzmuHKMm2UZoS3RcQQSfsDE4EFafUmwL0RcVze9jcAF0TEVen95sDbEbFM0meA83JJGMWUu7KMkyTawckSZl2OK8t0rFzihMginXIJE9s1GQT7keUR/jO3LCJez506jYjrge6S+ndi350k0R5OljCrWR4IN8x9wAhJHwKQ1EvSNnnrDwGmRsTaaYSkAekmGyR9lOzv4NVO7LOTJNrDyRJmNcs3y2yAiHglhfX+Q9KmafEPgCfT68OBM5vsdjDwNUmrgTeBw6MC56edJNFGTpYwq1m+RtgFOH3CzKztfI3QzMysBB4IzcysrnkgNDOzuuaB0MzM6poHQjMzq2tVORBKGifp92Vuc4yk7fLe/0TSPmVsv6RUCUlDJN0v6SlJl0vapFx9sDJaupA3zhrGAb+4ygkdZjWuKgfCDjIGWDsQRsSPIuKWMh/jzrwqMz8pss0vgd+ksmr/JatValXm2at+TM83FnDIG38nWJfQ4cHQrPZUZCCUdJSkB9LM6S+SNpZ0rKQnJd0BjMjbdqKkg/Pe5yc7fFfSI6mA9Zlp2VckzUzLrkrVXvYADgR+lY75wfx2JY2S9FBq66+5h+MlPS/pDEkPpnXbbuDnFlkSxZVp0cVkA7RVk6ULGfj8ZC7Rwazu1o93sxiAN1et4Vc3zqtw58ys3Dp9IJT0YeAwYERENABrgKOAM8gGwE+RN3NroZ39yAaRj0XEjsBZadXkiNg1LXscOC4i7gGuBU5Ns7Vn8trpQVY4+7CI2IGs2s7X8g61KCVL/An4Tivdai1Voh+wOCJWp/cLgILlXToyfcJasV5d0eDkbpPXriqW3GFmXVclSqyNAnYBZqaSmz2BPYDbI+IVAEmXA9sUbSGzD3BRRCwHiIjX0vJhkn4G9AV6Aze20s5Q4LmIyJVFuxg4ETg3vc/9KzgbGNtCOw8CH8hLlbgaaJoqoQL7FSztExHnA+dDVlmmlc9g5bJ0ITReyqZazTiuAODNjTfhd6vH8gp92aJvzwp30MzKrRKnRgVcnHctbShwOkUGBGA1qZ/p1GLu5hIV2WcicFKa3Z0B9CihPy1Zkf5cQwu/OJSYKrEI6Csp186WwIutHN86U4GUiY14m5O7TXZCh1mNqsRAOB04WNJ7ACS9C3gI2EtSP0ndyVIbcp4nm0ECfA7onl7fBHxJUq+8dgA2A15K7RyZ187StK6pJ4DBuQQJ4IvAHW39UKWkSqTi2reRFd4GOAa4pq3Hsg5UIGViU61mdLcHndBhVqM6/dRoRDwm6QfATZI2AlaRnYo8HbgXeInsNGMuMG8CcI2kB8gG0TdSO9MkNQCzJK0Erge+D/wQuB/4F/AI6wa/y4AJkk5h3UBERLwl6VjgijRTmwn8uR0frWiqhKTrgS9HxIvAacBl6fTtQ8CF7TiWdZQiKRPvwXc1mdUqp090AU6fMDNrO6dPmJmZlcDBvG2UTqN+o8niuyPixEr0x8zMNowHwjaKiIuAiyrdDzMzKw+fGjUzs7rmgdDMzOpaRQfCLpoysa2keyWtkPSdJuv+KullSXNb2F+SfivpaUlzJO1crr61x9UPvcCIM2+t34SFpQvhvB1h6X8q3RMzq5BanBF2dMrEa8ApwNkF1k0ERrey/35kpde2Bo4nq2FaEVc/9ALfm/wILyx+s34TFu44Cxb/G+74ZaV7YmYV0qE3y0g6imzQ2ITsIfevA0cD3yN7cP5JUgkzSROBqRFxZXq/LCJ6p9ffJav48jZwQ0SMl/QVsoFkE+DptL6BLGViz/TQ/ufJHrCfGhFXShpFNoDlHpz/WkSskPQ8WY3RA8gq1xwSEQWfrI6Il4GXJX22wLoZkga38rV8DrgkPWx/n6S+kt4XES+1sl/Z/erGeby5ag0AozdZ93Fvu+4JFj/cl+23355dd92VVatWcemllzbbv6GhgYaGBpYvX86kSZOarR8+fDjDhg1jyZIlTJkypdn63XffnaFDh7Jo0SKmTp3abP3IkSPZaqutWLhwIdOmTWu2ftSoUQwaNIj58+czffr0ZutHjx7NgAEDePbZZ5kxY0az9ft/cjf6N16alVRrvBT2PA02e2+z7cystnXYjLDGUyY2xEBgft77ggkUnZE+USxJYeXqNR1yvKoz84J1dUXjbc8KzepUh1WWkXQSWcmzl9OinmSlx+ZExNFpm1OAbSLipGIzQkm/Bp6IiAlN2t8TWC9lIiJOKNDORGAq8BTwu4gYmZaPAk6MiLFpRjgiIl6Q9DHg5xHR4nVFSacDyyLi7CbLB6fjDyuy3z+BX0TEXen9dOC7ETG72LE6qrLMiDNv5YUCg+HAvj25e/zeZT9eVcldG1z91rpl3XrAN+Z4VmhWI6qhskxNpkyUwQJgUN77iiVQnLrvUHp233i9ZXWTsFAgZcKzQrP61JEDYU2mTJTBtcDR6e7R3YAllbg+CDBmp4H8YuwODOzbE5HNBOsmYaFAygRrVmbLzayudNjMp1ZTJiQNAGYBmwNvS/omsF1EvC7pH8BeQH9JC4AfR8SFkk5Iffhz6v9nyG7wWQ4c29Y+lNOYnQbWx8DXVJGUCTOrP06f6AKcPmFm1nbVcI3QzMys6rnodhFOmTAzqw8eCItwyoSZWX3wqVEzM6trHgjNzKyuVeVA2EVTKfpIuk7Sw5IeTdcYC223i6RHUvrEb1PxANtQ7UyRqPv0DTOrzoGwg3R0KsWJwGOp9ulewK8lbVJguz+RFQvPJVC0llZhpWhHioTTN8wMKnSzTC2mUpCVgdsszfB6k8U1rW7yud8HbB4R96b3l5AN0De042u0nKULofFSJsbnYfYy+M/5sHH2O0hLCRoP/XsxA99+J0/Tf+2yN1et4Vc3zqvPIgNmdarTZ4Q1nErxe+DDZHVDHwG+EdG0mCUDyWqN5hRMnkj96vD0iZrRtG7o4vnFt81TLGWjWCqHmdWmSswIR5HVFJ2ZLo/1BPYAbo+IVwAkXQ5s00o7+wAXRcRygIh4LS0fJmm9VIpW2hkKPBcRT6b3F5Od5jw3vZ+c/pwNjG2hnX2BRmBv4IPAzZLujIjX87YpdD2wYGmfiDgfOB+yyjKtfIb6lWaDrFnJOK7Ivs03e8Dn10+R6N69O+PGjVtv1wlF0je26NuzgzttZtWkEtcIazWV4liy2WhExNPAc8C2TbZZQJY2kVOx5ImasQEpEnWdvmFma1ViIKzVVIp/k812kfRespnms/kbpJSJpZJ2S4P60cA17TiW5WxAikRdp2+Y2Vqdfmq0VlMpgJ8CEyU9QjbLPC0iFgFIakzXQyG7/jiR7JTwDfhGmQ2zgSkSdZu+YWZrOX2iC3D6hJlZ2zl9wszMrAQuut1GTqUwM6stHgjbyKkUZma1xadGzcysrnkgNDOzuuaBsBWSBkuaW2D57ZJavRupwH5lT9aoKe1MkSjEyRJmVgoPhFZd2pEiUYiTJcysVL5ZpjTdJF0M7ESWjHF0/kpJfwJ2JXtI/sqI+HFavitwHvAOslJto5rs91ngB8ABuYfv61o7UyQAGhoaaGhoYPny5UyaNImH/r2YPbUGNoFpK7NKd06WMLNCPCMszVDg/Ij4CPA6WWxUvv9ND21+hCzq6SMpi/ByshSKHcmKhK+t8CzpIGA88JlCg2Bdpk+0M0WiECdLmFmpPCMszfyIuDu9/htZlmK+QyUdT/Z9vo8sRiqAlyJiJkAuhSIlbnwSGA58ukk6xVp1lz6xASkS+Xr16sW4ceOcLGFmJfOMsDRNB6K17yUNIcspHJVmjP8kS7wolo4BWTHuzWg9aqp+bECKRCFOljCzUnkgLM37Je2eXh8B3JW3bnOyQuBLUurEfmn5E8AW6TohkjZLRb0hKwg+FrhE0vYd3vuuYANSJApxsoSZlcqnRkvzOHCMpL8AT5Gl1R8AEBEPS3oIeJRspnd3Wr5S0mHA7yT1JLs+uE+uwYiYJ+lIstSLAyLimU79RNVmA1MkCnGyhJmVwukTXYDTJ8zM2q7U9AnPCM3MymTVqlUsWLCAt956q9JdqSs9evRgyy23pHv37q1vXIAHQjOzMlmwYAGbbbYZgwcPzt0hbh0sInj11VdZsGABQ4YMaVcbvlnGzKxM3nrrLfr16+dBsBNJol+/fhs0C/dAaGZWRh4EO9+GfuceCM3MrK55IKxnZUx6yOfUB7Ou5ctf/jKPPfZYi9uMGzeOK6+8stny559/nr///e9tPmax9iqhKgfCjogqkjRG0nZ5738iaZ+W9mlj+++UNEXSHEkPSBpWZLshku6X9JSky1NN0sooU9JDPqc+mJWuWn5pvOCCC9huu+1a37CA9g6E1aSe7hodA0wFHgOIiB+Vuf3vA40RcZCkbYE/0CRtIvkl8JuIuEzSn4HjyB7Q71xlTHrI99C/FzPg7X48x7vWLnPqg1lzuV8a31yVFYjP/dIItPv/lbPOOosePXpwyimn8K1vfYuHH36YW2+9lenTp3PRRRdx9NFH8+Mf/5gVK1bwwQ9+kIsuuojevXuz1157cfbZZzN8+HAuvPBCfvnLX7LFFluw9dZbs+mmm/L732fzkhkzZnDOOeewcOFCzjrrLA4++GDGjx/P448/TkNDA8cccwynnHIK48eP5/bbb2fFihWceOKJfPWrXyUiOPnkk7n11lsZMmQI1fQMe0VmhJKOSrOmRkl/kbSxpGMlPSnpDmBE3rYTJR2c935Z3uvvSnpE0sOSzkzLviJpZlp2laRekvYADgR+lY75wfx2JY2S9FBq66+SNk3Ln5d0hqQH07ptW/hY2wHTASLiCWBwKrmW/7kF7A3kzgdcTDZAF/qOOjZ9ooxJD/mc+mBWml/dOG/tIJiT+6WxvUaOHMmdd94JwKxZs1i2bBmrVq3irrvuYocdduBnP/sZt9xyCw8++CDDhw/nnHPOWW//F198kZ/+9Kfcd9993HzzzTzxxPoVn1566SXuuusupk6dyvjx4wE488wz+cQnPkFjYyPf+ta3uPDCC+nTpw8zZ85k5syZTJgwgeeee44pU6Ywb948HnnkESZMmMA999zT7qfqHq4AABVYSURBVM9Zbp0+I5T0YeAwYERErJL0R+Ao4AxgF2AJcBvwUCvt7Ec2iHwsIpZLyk1BJkfEhLTNz4DjIuJ3kq4FpkbElWldrp0ewESyotlPSroE+BpwbmpvUUTsLOnrZMW1v1ykSw+T1Q+9S9JHgQ8AWwL5F+D6AYsjYnV6vwAo+Ktfh6ZPlDnpIZ9TH8xKU+yXww35pXGXXXZh9uzZLF26lE033ZSdd96ZWbNmceedd3LggQfy2GOPMWJENs9YuXIlu++++3r7P/DAA+y55568613ZP6eHHHIITz755Nr1Y8aMYaONNmK77bbjP/8pfG/BTTfdxJw5c9Ze/1uyZAlPPfUUM2bM4IgjjmDjjTdmiy22YO+992735yy3SpwaHUU24M1Mg1FPYA/g9oh4BUDS5bSezLAPcFFELAeIiNfS8mFpAOwL9AZubKWdocBzEZH7274YOJF1A+Hk9OdssoGumDOB8yQ1Ao+QDeSrm2xT6B7fzj8/0FLSw/7nFN7n/7d37tFVlVcC/20gEDAIWBisoAasIAgxD8JLEqCWpIWIteBMobSCrjWotXXo0FFrebm00ymZpRarURnFdlmpLY/6wDZFCYEMFFNMwQIK1FuLOkZSVB7yzJ4/vu/CzfXmeW9yk9z9W+us+53vfK+z77l3n/195+zdQL6fP7jGdA9Y1AfDiMRFPbvG/KYxKSmJ1NRUnnrqKcaOHUtaWhobNmxg//79DBgwgEmTJvHss8/WWr++6couXbrUW1ZVWbZsGfn5+TXy161b12pfLYnH1KgAT6tqut8GA4upXSGcxo/TTy0GHy6pLczRCuB2VR2OszKTGzCeujjhP89Qx42Dqn6iqnNUNR0Xwb4P8HZYsYNAz5AoFP2B9+rpP/bEONJDKBb1wTAaRnOFCsvNzaWwsJDc3FxycnIoKioiPT2d0aNHU1ZWxr59+wA4duxYDWsPYOTIkWzcuJFDhw5x+vRpVq1aVW9/3bt35/Dhw2f38/PzefTRRzl16hQAb731FkePHiU3N5eVK1dy5swZ3n//fTZs2BDVecaSeFiErwC/FZEHVLXST2m+jrOmPoeLAH8DbqoRIICzIJ8DrgOCzuSKgYUi8svg1Ki3CrsD74tIEvANIPgY1mF/LJzget4XVHUf8E1gY2NPSkR6AsdU9SRu+rQ0POiuqqqIbACmAyuBG4HfNravqGmGSA+hWNQHw6if4G9k6e/f5L2PPuWinl35fv7gqH87OTk53H///YwZM4bzzjuP5ORkcnJy6NOnDytWrGDGjBmcOOHu7++77z4GDTo3+davXz9+8IMfMGrUKC666CKGDh1Kjx496uwvLS2NTp06cdVVVzF79mzuuOMOAoEAmZmZqCp9+vRh7dq1XH/99bz66qsMHz6cQYMGMX78+KjOM6aoaotvuDXCCmAHbspxNDAHeAunhB4CHvZl+wJbgW3AfwJHQtq5C/cUaAXwI593K84SKwGWASt8/tW+7OvAZTjLcbo/do3P3wk8CXTx+QGgt0+PwE3f1nZOY3AhmvbgplN7hRxbB1zk0wP9uewDfh3sq64tKytLDcNo/ezatSveQ4iaw4cPq6rqqVOntKCgQFevXh3nETWMSLIHyrUBOsnCMLUBLAyTYbQNdu/ezZAhQ+I9jKiYP38+69ev5/jx4+Tl5fHQQw+12rW9UCLJ3sIwGYZhGI2msLAw3kNocUwRNhIRmQPcEZZdpqrfjsd4DMMwjOgwRdhIVPUp4Kl4j8MwDMOIDa3S16hhGIZhtBSmCJsJEdnk3blViMh7IrI27Hi2iJwJdR/XrDQi0kRrcQRsGIbREpgijAJxRJShquaodxoAbOGchxpEpCPO+XZ9Xm9iRwMjTVj0CMMwGsrYsWPrLZOamsrBgwc/k19SUtIkf6O1tRcNtkbYSEQkFXgZ5w91DPCgiNwCdAH2A3NUNdQxeHeco+05Ic18B1gFZLfIoBsRaeL1dz5ivJyBzvC7k87HuEWPMIxm5PD/wZP5cFNxDV+/bYFoHGeXlJSQkpLSIGXa3JhF2DQGAz8HJuHCKH1JVTOBcuB7YWWvB15R72VGRPr5vKK6Oohp9IlGRJqw6BGG0cLEOC5oIBBg2LBz4VALCwtZvHgxEyZM4M4772TkyJEMGjTobJSKyZMns2PHDgAyMjK49957AViwYAHLly8HYOnSpWRnZ5OWlsaiRYvOtp2SkgJAdXU1t912G1deeSUFBQVMnjy5RtDdZcuWkZmZyfDhw9mzZw+BQICioiIeeOAB0tPT2bRpEx9++CHTpk0jOzub7OxsysrKAKiqqiIvL4+MjIyz4ZxijVmETeNvqrpVRApw4ZfK/AunnXHToKHMAJaH7D8I3KmqZ+p6SVVjFX2ikZEmLHqEYbQgwd+nVrvP8Xc2q1V4+vRptm3bxrp161iyZAnr168/G7opNTWVTp06nVVAmzdvZtasWRQXF7N37162bduGqjJ16lRKS0vJzc092+7q1asJBALs3LmTyspKhgwZwk033XT2eO/evdm+fTuPPPIIhYWFLF++nFtuuYWUlBTmz58PwMyZM5k3bx7jxo3jnXfeIT8/n927d7NkyRLGjRvHwoULeemll3j88cdjLhezCJvGUf8pwB/0nAPxoap6c7CQ9506EngppO4IYKWIBHA+Rx8RkYgxCWNCXZEmItBcjoANw4hA6O+zjt9lrPja11wAnaysLAKBAOB8k5aWlrJ582amTJnCkSNHOHbsGIFAgMGDB1NcXExxcTEZGRlkZmayZ88e9u7dW6PdzZs3c8MNN9ChQwcuvPBCJk6cWG+/4axfv57bb7+d9PR0pk6dyieffMLhw4cpLS1l1qxZAEyZMoVevXrFUCIOswijYyvws6DDbhHpBvTXcyGdbsDFQDwerKCqA4JpEVnhj9d4ojSm1BVpIkLIpeZyBGwYRhghszWA+4yBVdipUyeqq8/d/B4/fvbv52wYpY4dO3L6tIsSl52dTXl5OQMHDmTSpEkcPHiQJ554gqysLMD5o7777ruZO3durX3WN10Zqd9wqqur2bJlC127fnb2qbldvJlFGAXq4ifOBp4VkR04xRgaxf7rQO3Bv1qCf98Diz/+7FZHBIqvZvSj7K4v8vaPp1B21xdNCRpGc9DI2ZqG0rdvXyorK6mqquLEiRO8+OKLdZbv3LkzF198Mc899xyjR48mJyeHwsJCcnJyABdW6cknn+TIEfcM4LvvvktlZWWNNsaNG8eqVauorq7mgw8+oKSkpN5xhodvysvL4+GHHz67X1FRAbiwUs888wwAL7/8MocOHapfCI3EFGEjUdWAqg4L2X9VVbNVNc1vz4ccm6Cqv6ujrdmq+pvajhuG0Y5pprigSUlJLFy4kFGjRlFQUMAVV1xRb52cnBz69u1Lt27dyMnJ4cCBA2cVYV5eHjNnzmTMmDEMHz6c6dOn11BgANOmTaN///4MGzaMuXPnMmrUqHrDN1177bWsWbPm7MMyP/3pTykvLyctLY2hQ4dSVOSeJ1y0aBGlpaVkZmZSXFzMJZdc0kTJ1I5Fn2gDWPQJw2gbtIfoE03lyJEjpKSkUFVVxciRIykrK+PCCy9ssf4t+oRhGIYRVwoKCvjoo484efIkCxYsaFElGC2mCA3DMIyoaci6YGvF1ggNwzBiiC03tTzRytwUoWEYRoxITk6mqqrKlGELoqpUVVWRnJzc5DZsatQwDCNG9O/fnwMHDhC1W0SjUSQnJ9O/f/8m1zdFaBiGESOSkpIYMGBA/QWNVoVNjRqGYRgJjSlCwzAMI6ExRWgYhmEkNOZZpg0gIh8Cf4uiid5AbEM6tw9MLrVjsomMySUyrVUul6pqn/oKmSJMAESkvCFuhhINk0vtmGwiY3KJTFuXi02NGoZhGAmNKULDMAwjoTFFmBg8Hu8BtFJMLrVjsomMySUybVoutkZoGIZhJDRmERqGYRgJjSlCwzAMI6ExRdjGEZEvi8ibIrJPRO6KcLyLiPzKH/+jiKSGHLvb578pIvktOe7mpqlyEZFUEflURCr8VtTSY29OGiCXXBHZLiKnRWR62LEbRWSv325suVE3P1HK5UzI9fJ8y426ZWiAbL4nIrtEZIeIvCIil4YcaxvXjKra1kY3oCOwHxgIdAb+DAwNK3MbUOTTXwd+5dNDffkuwADfTsd4n1MrkEsq8Ea8zyGOckkF0oCfA9ND8i8A/uo/e/l0r3ifU7zl4o8difc5xFk2E4FuPn1ryG+pzVwzZhG2bUYC+1T1r6p6ElgJXBdW5jrgaZ/+DXCNiIjPX6mqJ1T1bWCfb689EI1c2jP1ykVVA6q6A6gOq5sP/EFV/6Gqh4A/AF9uiUG3ANHIpb3TENlsUNVjfncrEIyH1GauGVOEbZt+wN9D9g/4vIhlVPU08DHwuQbWbatEIxeAASLyuohsFJGc5h5sCxLNd57o10tdJItIuYhsFZGvxnZocaexsrkZeLmJdeOGxSNs20SyYMLfh6mtTEPqtlWikcv7wCWqWiUiWcBaEblSVT+J9SDjQDTfeaJfL3Vxiaq+JyIDgVdFZKeq7o/R2OJNg2UjIrOAEcD4xtaNN2YRtm0OABeH7PcH3qutjIh0AnoA/2hg3bZKk+Xip4qrAFT1T7j1kUHNPuKWIZrvPNGvl1pR1ff851+BEiAjloOLMw2SjYh8CbgHmKqqJxpTtzVgirBt8xpwuYgMEJHOuIc+wp9aex4IPq01HXhV3Ur288DX/dOTA4DLgW0tNO7mpslyEZE+ItIRwN/hX45b5G8PNEQutfF7IE9EeolILyDP57UHmiwXL48uPt0buBrY1WwjbXnqlY2IZACP4ZRgZcihtnPNxPtpHdui24DJwFs4y+Uen3cv7qIESAZ+jXsYZhswMKTuPb7em8BX4n0urUEuwDTgL7in47YD18b7XFpYLtm4O/mjQBXwl5C6N3l57QPmxPtcWoNcgLHATn+97ARujve5xEE264EPgAq/Pd/WrhlzsWYYhmEkNDY1ahiGYSQ0pggNwzCMhMYUoWEYhpHQmCI0DMMwEhpThIZhGEZCY4rQaPeIyPUioiJyRUjeBBF5MQZtrwiPRhChzAQRGdvIdruJyDMislNE3hCRzSKSIiI9ReS26EZdo59UEXkjyjYWi8j8WvJVRL4QkjfP542Ips96xhOMIPK6iOwWkW2hkQ9EZGowioJ/b/SPvmyOiNzg62xorvEZrQ9ThEYiMAPYjHsZOB5MwL1v1hjuAD5Q1eGqOgznw/EU0BMXOSMuBJ0NNIKd1JT7dFrmhfP9qpqhqkN8//NEZA6Aqj6vqj/25a4B9viym3Byvk1VJzakkybIw2iFmCI02jUikoLz9nEzn1WE54vIGh9LrUhEOohIR2/lveGtsXm+nXTvVHmHr9MrQl8B710EERkhIiXi4hzegvsjrvBWRx8RWSUir/nt6ghD/zzwbnBHVd9U57rqx8Blvq2l3kp8RVysvJ0icp3vP9VbNk+IyF9EpFhEuvpjWSLyZxHZAnw7ZPypIrLJt7U9aMV6i3aDiPwSp9gQkXvExahbDwyu4ytYi49W4D31fAx8GNJnnohs8f392n9fiMhCL5s3RORxERcZxMv0v7yV95Y0wCm6Otdn3wO+69uYLSIPi0g68BNgspfnImAcUORl29F/vua/97l1yGOWH1OFiDwm57wTHRGR+728t4pIX5/f119Hf/bb2LraMZqZeL/Rb5ttzbkBs4D/8en/BTJ9egJwHBdnrSMuRMx0IAsXOiZYv6f/3AGM9+l7gQd9egU+Ph0QAHr79AigxKcXA/ND2vwlMM6nLwF2Rxh3OlAJbAHuAy73+amExEvEOc4/36d74zx4iC93Gkj3x54DZkU4l6XB9oBuQLJPXw6Uh8jqKDDA72fhFEA34Hzf5/wI57AYmA+sBobhPBndiPPHOcKPtxQ4z5e/E1jo0xeEtPMLvIcfX/e/fXoysD5CvzVkFPwegU99ejbwcHg6pP0RPv2vwA99ugtQjovdGS6PIcALQJLffwT4lk9ryNh/EtLer4B/8+mOOF+3tbZjW/NuFn3CaO/MAB706ZV+f7vf36bOWkBEnsVZA68AA0VkGfASUCwiPXAKcaOv9zTOPVtT+RIwVM6FPzxfRLqr6uFghqpWeAsqz5d/TUTGAJ+GtSXAj0QkFxcrrx/Q1x97W1UrfPpPQGqEc/kF8BWfTgKCltIZajob36YubiVADrBGfQw6qT8q+0qcNZ6Pm4qc4/NH4wJEl3lZdMYpfoCJIvIfOGV7Ac7t3Qv+2OrQc6qn7yBNiTWZB6TJuTXgHrgbhJPUlMc1uJuD1/x5dMXdxODLBtei/wRM8ukvAt8CUNUzwMci8s062jGaEVOERrtFRD6H+8MZJiKKu/NW/wcLnw0Jo6p6SESuwv1pfxv4Z2BeA7s8zbnlhuQ6ynUAxqhquFILH8wR3J/+ahGpxllAq8KKfQPoA2Sp6ikRCYT0fSKk3BncH6tQeyiceTifkVf5MR4POXY0fHh1jT2MF3CWZ7mqfhJyAyA463tGaGERScZZQyNU9e8ispia8gye1xka/h+WAexuxJiD4/uOqtZwFC0iE6gpDwGeVtW7I7RxSr15R/3jrasdoxmxNUKjPTMd+LmqXqqqqap6MfA2zvIDGCnOq34H4F+AzX6Nr4OqrgIW4KZSPwYOhaxHfRPYyGcJ4O7owTnvDnIY6B6yXwzcHtzxFlgNROTq4DqkOK//Q4G/RWirB1DpleBE4NK6BKKqH+Gsj6AMvhHW1vuqWu3Psbb1qVLgehHpKiLdgWvr6fNT3LTn/WGHtgJXi3+qVNyTsoM4p/QO+jXDOp/KrQ9x67SFwLJGVv09cKuIJPl2BonIeRHKvQJMF5F/8uUuEJE6vwdf51ZfvqOInN/EdowYYIrQaM/MANaE5a0CZvr0FtzDJ2/gFOQa3NRiiYhU4Nb/gnfnNwJLRWQHbv3u3gj9LQEeEpFNuLv/IC/gFEeFV6bfBUb4BzB24R6mCecyYKOI7ARex61PrVIXK7HMP0SyFHjGt1WOU2p7GiCXOcDPxD0sE2qVPgLcKCJbcdOi4VYgAKq6HbfGVYGT56b6OlTVlb5eaN6HuDW6Z71ctwJXeGX9BG4dci0uFFBjuUz86xO49dFlqvpUI9tYjnvCdbu4V0weI4JFp6q7gB/iptF34NabP19P23fgpn934qZMr2xiO0YMsOgThmEYRkJjFqFhGIaR0JgiNAzDMBIaU4SGYRhGQmOK0DAMw0hoTBEahmEYCY0pQsMwDCOhMUVoGIZhJDT/D/JVcWUQm8llAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "from causallib.evaluation import PropensityEvaluator\n", "\n", "evaluator = PropensityEvaluator(ipw)\n", "evaluation_results = evaluator.evaluate_simple(X, a, y, plots=[\"covariate_balance_love\"])\n", "\n", "fig = evaluation_results.plots[\"covariate_balance_love\"].get_figure()\n", "fig.set_size_inches(6, 6) # set a more compact size than default\n", "fig;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that prior to applying IP-weighting (orange triangles), the maximal absolute mean-difference across treatment groups was around ~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 percentage of people owning a high-school diploma (`no_degree`). \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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.7" } }, "nbformat": 4, "nbformat_minor": 2 }