{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Doubly Robust Methods\n",
"### AIPW\n",
"\n",
"#### Background\n",
"\n",
"Augmented Inverse Propensity Weighting (AIPW) is a modification of standard Inverse Propensity Weighting to achieve double robustness. We first consider basic IPW, which considers a sample weight, or propensity score $\\hat \\pi (X_i)$, in the model.\n",
"\n",
"$$\\widehat{ATE}_{IPW} = \\frac{1}{N} \\sum_{i=1}^N \\left[\\frac{D_iY_i}{\\hat \\pi (X_i)} - \\frac{(1-D_i)Y_i}{1-\\hat \\pi (X_i)}\\right]$$\n",
"\n",
"The augmenteed IPW, AIPW, as presetned by Glynn and Quinn, 2009 includes the outcome model in such a way that ensures doubly-robust estimation. This equations below reqrites the AIPW formulation such the basic IPW is seen clearly first along with the ourcome model augmentation. \n",
"\n",
"$$\\widehat{ATE}_{AIPW} = \\frac{1}{N} \\sum_{i=1}^N \\left(\\left[\\frac{D_iY_i}{\\hat \\pi (X_i)} - \\frac{(1-D_i)Y_i}{1-\\hat \\pi (X_i)}\\right]-\\frac{(X_i - \\hat \\pi (X_i))Y_i }{\\pi (X_i)(1-\\pi (X_i))} \\right) [(1-\\hat \\pi (X_i))\\hat{\\mathbb{E}}(Y_i|D_i=1,X_i) + + \\hat{\\pi} (X_i) \\hat{\\mathbb{E}}(Y_i|D_i=0,X_i)]$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"# Hide Cell\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"pd.set_option(\"display.precision\", 3)\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('seaborn')\n",
"import seaborn as sns\n",
"\n",
"from sklearn.linear_model import LogisticRegression, LinearRegression"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Dataset\n",
"\n",
"We will use a simulated dataset based on The National Study of Learning Mindsets. This was a randomized study conducted in U.S. public high schools, the purpose of which was to evaluate the impact of a nudge-like, optional intervention designed to instill students with a growth mindset. The study includes measured outcomes via an achievement score, a binary treatment of a growth mindset educational intervention, and 11 other potential confounding factors that could be parents of both the treatment and outcome. The first 5 rows of the dataset are shows in the table below. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
schoolid
\n",
"
intervention
\n",
"
achievement_score
\n",
"
success_expect
\n",
"
ethnicity
\n",
"
gender
\n",
"
frst_in_family
\n",
"
school_urbanicity
\n",
"
school_mindset
\n",
"
school_achievement
\n",
"
school_ethnic_minority
\n",
"
school_poverty
\n",
"
school_size
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
76
\n",
"
1
\n",
"
0.277
\n",
"
6
\n",
"
4
\n",
"
2
\n",
"
1
\n",
"
4
\n",
"
0.335
\n",
"
0.649
\n",
"
-1.311
\n",
"
0.224
\n",
"
-0.427
\n",
"
\n",
"
\n",
"
1
\n",
"
76
\n",
"
1
\n",
"
-0.450
\n",
"
4
\n",
"
12
\n",
"
2
\n",
"
1
\n",
"
4
\n",
"
0.335
\n",
"
0.649
\n",
"
-1.311
\n",
"
0.224
\n",
"
-0.427
\n",
"
\n",
"
\n",
"
2
\n",
"
76
\n",
"
1
\n",
"
0.770
\n",
"
6
\n",
"
4
\n",
"
2
\n",
"
0
\n",
"
4
\n",
"
0.335
\n",
"
0.649
\n",
"
-1.311
\n",
"
0.224
\n",
"
-0.427
\n",
"
\n",
"
\n",
"
3
\n",
"
76
\n",
"
1
\n",
"
-0.122
\n",
"
6
\n",
"
4
\n",
"
2
\n",
"
0
\n",
"
4
\n",
"
0.335
\n",
"
0.649
\n",
"
-1.311
\n",
"
0.224
\n",
"
-0.427
\n",
"
\n",
"
\n",
"
4
\n",
"
76
\n",
"
1
\n",
"
1.526
\n",
"
6
\n",
"
4
\n",
"
1
\n",
"
0
\n",
"
4
\n",
"
0.335
\n",
"
0.649
\n",
"
-1.311
\n",
"
0.224
\n",
"
-0.427
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" schoolid intervention achievement_score success_expect ethnicity \\\n",
"0 76 1 0.277 6 4 \n",
"1 76 1 -0.450 4 12 \n",
"2 76 1 0.770 6 4 \n",
"3 76 1 -0.122 6 4 \n",
"4 76 1 1.526 6 4 \n",
"\n",
" gender frst_in_family school_urbanicity school_mindset \\\n",
"0 2 1 4 0.335 \n",
"1 2 1 4 0.335 \n",
"2 2 0 4 0.335 \n",
"3 2 0 4 0.335 \n",
"4 1 0 4 0.335 \n",
"\n",
" school_achievement school_ethnic_minority school_poverty school_size \n",
"0 0.649 -1.311 0.224 -0.427 \n",
"1 0.649 -1.311 0.224 -0.427 \n",
"2 0.649 -1.311 0.224 -0.427 \n",
"3 0.649 -1.311 0.224 -0.427 \n",
"4 0.649 -1.311 0.224 -0.427 "
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_mindset = pd.read_csv('learning_mindset.csv')\n",
"# print(df_mindset.info())\n",
"df_mindset.head()\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"## Hide cell for book\n",
"# Convert categorical values to binary indicators (one-hot)\n",
"categ = [\"ethnicity\", \"gender\", \"school_urbanicity\"]\n",
"cont = [\"school_mindset\", \"school_achievement\", \"school_ethnic_minority\", \"school_poverty\", \"school_size\"]\n",
"\n",
"df_categ = pd.concat([\n",
" df_mindset.drop(columns=categ), # dataset without the categorical features\n",
" pd.get_dummies(df_mindset[categ], columns=categ, drop_first=False) # categorical features converted to dummies\n",
"], axis=1)\n",
"\n",
"# df_categ.head()\n",
"\n",
"T = 'intervention'\n",
"Y = 'achievement_score'\n",
"X = df_categ.columns.drop([T, Y])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Understanding the data and potential bias\n",
"\n",
"We begin by visualizing the achievement scores of treated and untreated cohorts with no control or consideration for the other variables. It is clear form the plot below there is an impact of the treatment as the average of the treated group's achievement scores is clearly higher. But we can intuit a positive bias in this measurement. We should note again the intervention was an option to take a growth mindset course. So although the option was offered in a random fashion, *it is highly likely students who opt-in to the treatment are likely to have the features to provide higher achievement scores regardless*. Thus, we might hypothesize controlling for this bias would decrease the ATE from the *naive ATE* (meaning no adjustment or simple difference of means of the treated and untreated groups)."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEfCAYAAABI9xEpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABgqklEQVR4nO2dZ3hcxdWA321aaVe9N8tyHXcbV2xMh1ASeggJJNSQUD5ISCUVkhCS0AIpJJCEHieEZlMChGLANu69aWy5qffetn8/7l1rLWRLlu5au9K8z7PP7t5y5uzde+fMnDlzxhQIBFAoFArFyMM81AooFAqFYmhQBkChUChGKMoAKBQKxQhFGQCFQqEYoSgDoFAoFCMUZQAUCoVihGIdagWMRAjxd+Am/etOKeU0g+QGY2W3Siln9eP4j4DT9a8pUsomI/QYaQghTgecUsr/9vP4xcDtwGIgE3ADB4B3gIeklDXh0jWSEUIUol0HgGVSykvDXN4ZwPLjOGWMlPJgeLTpGyFEBnCtlPLhodJhMAghxgDnSSn/erznDpsegBAiDrgyZNNUIcT8odJHMXCEENlCiH8BHwET+3nOncAnwJeBfCAGiAemA98HtgshJoVFYUXUIoT4FrAX+NpQ63K8CCFsQohfALuA8wciYzj1AC4DEntsuwFYNwS6/BRI1z+3D0H50c75aBV5v9BbuA8DJqAJeBzYASQBNwOz0XoET9DdM1OEjx1oz2Mor+nvtcA3euwbyp7Zo0NY9mDJA34+GAHDyQBcp797gUYgA/iKEOI7UsrOE6mIlHLliSxPwYV038s/kFL+LbhDCLEE2IdmkE8TQqRKKRuGQMcRg5SyDlgauk0IEfzYIaVciiIiGBYGQAiRB5yjf30f2A3chdYCvAxYcpTzZgE/AU7Tjy0D3gPul1KWHuWcXOAhtErHpJf3fSnl/pBjPuIoYwBCiCuB7wAzAQ+wEc0//V99fzZQivbfrJVSntyj/J8Dv9C/XiqlXKZvPwm4V/8tsYAEngb+JKX0hZx/EBgNvAt8E63l/Dk0w/myrlsM8ABwOWABVgJ3Sin39dClv2UGr8c+YIZ+ztVAGrAduDfk9z9DtzEH+L0Q4vfAmVLKj+gdZ8jnc4QQLwSNvpSyRQjxVWCUvt8XeqIQYizwM+BctEZDJbAKuE9KubvHsQ7g22i9k7H6NdsEPC6lfLnHsaG/+ZvAk2iuqQ1SylOO8/rFAz9C+z8K0f6TSrR79R4pZflRrkuvCCFOBh4E5gINwH+An0spW4UQVrTnIAtoATKllK6Qc7+Pdm+A5jd//njKPoo+Z9A9ZvB1tOtxJdAJfENK+YoQIgb4Ltq9MQZoBj4EftHL/5QC/ADt2c/XN1cCb+jHN/cyTjFTH+t7Vkp5fch96AJS0Z65r6J5GVYBt6GNq/wAuAXt3tkJ/FRK+W4PfbL08y9Ca4hUoPWI7gttjAgh7gXu0b+OQqtjvgWMQ6sT/go8IqUMCCGuR7tXglyi6/8LKeW9n73KvTNcxgC+RvdvWQL8M2Tfjb2dIIS4AFgDfBHNPWBHu9C3AJ8KIfJ7OS0NWA18Bc1gJKI9lB/oD84xEULcg/awnQzE6eefCbwlhLgDQEpZBbytn7JAr6BC+Yr+XgsEK83P6XpdDCSjVSYz0bq3/z6KOvnAWuAKIAFIQXOXvAisQBtMT9F1vFD/jXEhv2UgZcagGZ7vo3VfY4F5wOtCiJlHOac/hPa4vgSUCyGeEUJcI4TIllK+K6X8u/5qDvkNs9EM8PW6PjFoxvFqYLUQYkbIsZlo1+vXaOMKTrR74EzgJSHE40fRLR2tNTxWl79Hl9ev6yeEMAFvAT8GJunH2YACtP9ouRAitf+XiuloFd9iXVYumlF7XwgRI6X0Av/Sj01EaxyEcrn+3km3W8dI7gWuRXs+UoGN+rP1JnA/INCuYwZwFbBOCLEweLIQwoZ2ve7Wj3Xqr/FojcI3jlMfM1oQwfeAbMCB1lh4H+063Y/2X8ShGdQ3e9w3+cB6tEZArq57oa7Lp0KItKOU+wCay3IKWt00Hq3heftx6t/njxsOXKu/dwJLpZQb0VpTAGcJIUaHHiyEcALPoV1YH/AYmnUP3vj5wH29lJOvH38bmqEIViaFwAXHUlAIMZdu674DzQ96I1oFBPBQSGUfatkP+8KFEHPQKgGAJVJKj14pP6//lnq0G/8aNEMD8EW9BdyTqWhRMreitYCDkU4XAhPQWiw3A9X69tFolR2DKHMUmj/+52hGe4++3aKXBfAH4I8h5zyF1pLb0Ys8AKSUqzmyl5eC1np7AagQQqwUQlwUeo5esT6PVvkCPIt2DwTLTtJ1CfIkEIwqW6Pr+wOgTt92qxCit8ZGElor8la0VvzTx3n9FqC1iEFr8V+nHxuMjJqAdi/2l7FoLdBb0Xp7Qf3no93XoD0bQa4IftB7vwv0r8uklG3HUW5/yUe7Dl8DfqNHB92JVukCvK7v+y7avRkPPCuECNZllwJBg/AuWoPpG2gDvQCn6mNGPccpDurfQ/9z0IztfLR79pt0j1eMRuul/FPXZ7W+3Yp2HwV5HO2+9wO/R2tc/A6t9yj0z73xFeAf+vFvhWwP/tcf6voEWafrf7TGV69EvQtIj/SZrH99U0rZqn/+J/BLNDfN9XS7TQAuoXuQ9o9Syrt0WUvQWhaddBuQnlwopSzSj09A60qDVjEfq3Vxk64LwAVSyjJdxjK01nwM2sN9D1prp1bX5StorQzQboYgz+rvF6P1YAC+LaV8Qf+8RAgxEZiFVlm9wGe5Rkq5QtfjNLofsgeD3UjdJfUrfXvQjTKYMr8bDFcTQpSiRfqA1vtCSrkptAUFbO+nz/h6oBitZZUQst0EnILWywjtHi9Ea12B1mi4Xv/8T92FkAFs1l0PhWj3DGgRF6dJKT36b3gXzQ1kQavgn+pFt3tDQ/SEEFfR/+sXHyJnLfCSlLJTCPEyWsNnN5rrob+4gbODYZdCiM10u0KuAh6VUm4WQmxH6y1cLISw6b/3Mrrv4X8SHqrR7ktvyLZg42AvcJmU0q/rvh+tFzIBzd22HC0S7DI0V+PDUsp2/Vgb8GddTq6U8lNgacjYRPMx7rNfSyl/pctJR+sFAnwqpfyqvv0gWs8Z9OdEf3a+oG97Wkr5Hf3zv4QQsWjunauFEHdKKTt6lPmilPLrupzXgSq0eyH4nJQIIf4XcnzlQMZWot4AcKS/OLQVGDQAANcLIX4ppQy2cueEHPdO8IO+/1yOTnmw8tc5FPI5rufBPZgb8rk05MYL5WRdD48Q4p9oXfNpQohpaBVPsDewTUq5uRe5zwshevPJzhNCWEL9ymgtktUh36tDPn8U8rky5HPsIMsErRUbJPRaxjII9ArqHiHEA2i9sfPQ/svQ3t89QoilUsotHOUe0GUdERKo+4uDPBWs/PVjtwkhVqO5VMbrLqeqHuqt6PG939cPzd9cilap/BT4nl7eCuBtKeWaXs49FptDY+6llB8JIRrQ3C2TQ457Dq1xkwKchdaaDrp/6vXv4WB1aOWvj38Ee70TAN8xnp3lUspqtIr9fWCRPt5xCkdGf9mOU6dQF+PxPCdz6DaYNwkhbuKzxKG5/lb32H74OZFStgshyuh2ARpGVLuA9NZZaLjga0KIgD4YEjpgWYjuvtBJCvlcR//pGT3iCfnc17VM6Yf83JDPoW6gr6DdwMH9zxyn3Di6XR1BOnq0skIr6sajbA/ezAMtE7QwzSCukM8mBoEQwimEGCWlbJdSviylvFlKWYj2EG4KOTTY7T+eeyDUT9vbgGvott788T3l9/v66YPZZwDL0NwGsWj38s/Rxik2672G/tJbBFTQrRE6mP5Puv/7L+q+6qAr6j+hRtBgel6r5H6elwuHY+MfQ/tN/0NrBE5G6ykFOd57rSXks9HPCRz53Adp6vHd1csxgybaewAX0/sD1xs3ovnNoNt3D0c+3MFZde1HmTXqP24Nuwm6pnxovsPeVuI5PGdAb1luQvOZX0r3w+nlyO53a8jnn6D1FHqjZxezZ8s8FO8x9g2mzL7KHRBCiCI0f6pHCJEROtCru5TuRfMdQ7fr5Vj3QJ5+brBiD70X8npRIXRbbS/7e4YhH9f10yPMLtUHe88BTkVrlU9Bcxe9gOan7g+9VTbB3x963SqFEO+hzcm4FM39FKwvenPtGcWxrtV2jh73HuyN/xRtzAC0cbwnpJRleuv77wPU6Wj37PE8J/+ie4ynJxt72Wb4c9Ib0W4AQt0/j9P7w3c32mDb5UKIJL1yCG0RXojWUgjyEjBHCFEOjJdSdhmk62a0h9UC1IX43oWu4w6O7FKC1guYjfagZ+vb3u5hnDaHfLYE/YB65MSTaINbmwyeCxHuMkMNbX9aa5vQDIANbQzg3h77Tw35fCDknCAXoukd5A9o90sz2sBvaPf8eiHEY8Hek+6eCw467pZS9nYP9mw49Pv6CSHOQQsfnIoWbvwf9IpECLET7d6Yq0fwuHspuyczhBBTpZQ7dRmL0MY7ALb2OPY5NAOQTvcY2gHdfx4ujrhWesjmAbTQzzzgo2BYtRDiCrTe0E66K9FgD8+LFhIZrKRPOkp5AbR7bFA90KOwJeRzVqiPXgjxU7Sxqp0cnxcilON9Tj5D1BoAPSwvOP25ErgjODjU47jpaL7LODR30RNoYXn1aC2fO/UBolW6vKBveJWBlT9olfkN+ueXdV91LVo0Q3DQ80sc2RpYghb6Zae7p/NMD7lL0bqiKcDPhBBJuowr6X4YHkcbWDaKcJcZ2nO4SK+IV0gp9x7l+D+g/bcm4OdCiAloPlQTcDbdg+duuseJPkVzC0xGi6FeghZtMZduX7cMGaxfjlbZTAM+EUL8A+33/wDNqEP34GBfLKX/1y+R7hbtK0KIR9Dum7loRg+0san+VP6gXZO3hRC/QXNb/jRkX8+B3aVo7o9EunsOvc6pCTNPo7lyUoEPhRB/RgvH/AXaNXTTHYARdNdYgRf0IIvTOXL2sT3kcwda73qcEOIGoHEgg6m9IaU8JIT4EK23dpY+cP8a2j30Q7T/YidaRNhACH1O5gghvgZUSCk/6K+AaB4DuIZuA/ZSb5W/Tmh39UbQBlXQokbcaH/CbXSHc4Hm0/22kcrqLf5H9a+ZaBX7s3RX/i8Br/Q4p4Fu1wVoRuvNHse0oE2e8aO1gL+L9pAGK5IiusNPDeEElBnaQj4TLRxu3jH0WYNWEQdbc1ejVRpPod0nJrQu9TeCFbo+4H8d3RXGV9DulW/r31s4cg7JtXSHrS5Ecyc8SHfr+VEpZb8iY47z+r1G9z08Bi1M9d9ocekWXUZQ5/7wFtr99zjwJ7p7lu/To3Gh9+Be6nF+uKJ/jsVDdPfCTkK79n+g28f+3eD/ypFRWFehXddvcqS7JnSOT/Bec+rn3oax3EK3Z+IKtP/ybrR7sg34ekhwynEhtRnXwQmreWg9tquOR0Y0G4BQ98/RfGugxUsHB2vmCyGmAkgp3wQWoT1g9UAX2gP+R2CelLKyF1mDQg83/Spab6MVzee/FS0c7JqjGLHQweB/9dbSk1K+ivZblqF1J11oIZGPAKfqN4qhhLNMqc04vhMt7C8ot7GPcx5Ci8T5N9pD4UZrIe1Hu4ZzpZTP9jhnPVpL+nm0MDu3fvxzwJygm0Q/tgyt8vkx2n/WgWYklqOFJt51nL+xX9dPrxyuRzNAK9GiULz6+1K0kNRX6D9r0Fxiy9H87RVosegXH+X+C322Nsses25PBLohOhttvspONL0b0FymX5BS/ink2KfQwkaDxx1Cd+vS7ZMPnRNyO1roaBtaPbAHA9F7rSeheR6C92U52nVdOIAorp5cB2xA+61VaO7DfmMKBAZkfBQKxQhAHJl65Pu6oVUME6J2DEChUIQHoaXNHo8Wd363vtlDeKN/FEOAMgAKhaIns/msr/+xXia4KaIcZQAUCkVP9qKNMSQBJWhjIr8ZUo0UYUGNASgUCsUIJWp6ALW1rYZbqpQUB42NvU1WjVyiUWeITr2jUWdQep9IokHnjIyEo04Si+Yw0EFjtVr6PijCiEadITr1jkadQel9IolGnUMZ0QZAoVAoRjLKACgUCsUIRRkAhUKhGKEoA6BQKBQjFGUAFAqFYoSiDIBCoVCMUJQBUCgUihGKMgAKhUIxQlEGQKFQKELYtGkD99zzoyO2/eUvf+S//32j1+OrqqpYufKTQZf7yisv9vvYY+lzPERNKgiF4nhYumJ/r9svPXXsCdZEMdzZtGk9hw4dZPHi0wYl59lnn+KKK45rQa9BowyAQqGISP7zYTHri2oMlTlvUiZfOmv8gM//7nfvxGazUllZwVlnnctdd93BCy88Q1dXF9Onz+Df//4nyckptLa28uCDj/Lww7+lrKwUv9/PzTffyuzZc1m+/H1effUlgok477vvAZYte4WWlmYeeui3fPvb3+PBB+//zHkfffQBzz77D5KTU/B4PIweXTjo66FcQAqFQtEPTCYT1dWV3HffA/z1r0+zZMlzWCwWvvrV6zn33PNZvPh0AM4993wee+xx3nrrdZKSkvnzn//Gb3/7MI888gAApaUlPPjgY/zpT09SUDCadetWc911N5GYmMT3vnc3b7yxtNfzHn/8Dzz66OM88sifiI2NNeQ3qR6AQqGISL501vhBtdYHit0ei9vtOWJbZ2cHMTF2xo4dj9VqxWq1Yrf3XgkXFIwGYN++YrZt28yuXTsA8Pm8NDc3kZKSyn333YPD4eDQoYNMmzbjiPN7O6+hoR6n00lSUjLAZ84ZKMoAKBQKRQiFhYXs3Supq6sjPT0dl8vF1q2bmTBBYOolsbLJZCIQ8B/+bjZrjpXRowvJzMzk2mtvxOXq4tlnn8JisfKPfzzBK6+8CcBdd91+2BUUfO/tvISERNra2mlsbCQlJYWiol1kZmYN+rcqA6BQKBQhOJ3x3HHHXfzgB9/Cbo/F6/VwxRVXkZeXz4YNaz9z/Lhx43nuuaeYOHHSEdsvueRyfve7+/i///sG7e1tXHbZlTidTqZPn8mNN36VuLg4EhISqKurBaCwcAy//OXPuPvun33mPJvNxo9//HO++93/IyEhCavVmKo7alYEC8eCMBkZCdTWthotNqxEo85w4vU2IgpIXesTSzTqHQ06qwVhFAqFQvEZlAFQKBSKEYoyAAqFQjFCUYPAimGNz++noq6DpjYXXp82jDRzfDqF2QmYegvpUChGEMoAKIYtFXXtbN/fgMvjO7yttKaN11cdZExOAlefM5FxeUlDqKFCMbQoA6AYdgQCAXYfamRvWTMWs4nxeYnkpDmxWc1MLUxl5fZKNu+t49fPb+SiRYVcsngMZrPqDShGHsoAKKKa3sI995Y1sbesGWeslXmTM0l0xBzed9LEDE6amMGe0ib+/uYu3vj0IGW1bdxyyVRsVsuJVF0Rgfzxj79Hyt00NNTT1dVFbm4eyckp3Hff7wYkb9myV/n85y/uV9z+0qUvU19fz003fXNAZQ0EZQAUw4rK+nZ2H2oiLsbComnZxNl7v8Unjkrmnhvm8fhrO9i8t45HX9rGt744gxibMgIjmTvuuAuA//73DQ4dOsitt94xKHnPP/8055//ecMmbhlN2LQSQpiBx4GZgAv4upSyuJfjngQapJR3h0sXxcigy+VlS3E9ZrOJBVOyjlr5B3HG2vj2lTP4y9KdbCmu46/LdnL75dOwmFVwXCRw76c/5Y19Sw2VedG4S7l30X3Hdc6vf30vzc3NtLQ088ADj7JkyXNs3boJvz/AzTffxNy5i9m8eSNPP/03ALq6uvjpT3/Btm2baWio5957f8xvfvMwf/3rnw6fd9VV13DWWeewdesWHnvsIRITEzGbLUydOs3Q39sX4TRLlwKxUsqFQoiTgYeBS0IPEEJ8E5gOfBxGPRQjgEAgwNZ99Xi8fqaPTSXRGdPrcTe+/P1eTjZD3By2FMPNf32Wi6drWR2dTjvt7S61hoCCOXPmctVV17B69SoqK8v5y1+ewuVycfvtN/H738/kwIH9/PznvyI9PYPnnnuK5cvf57rrbuKZZ/7Bvffe/5nzvvnNG5g3bwF//OMj3HvvrykoGM1DD/3mhP+ucBqAxcA7AFLKNUKIuaE7hRALgZOBJ4BJnz39SFJSHFjD4KPNyEgwXGa4iUadITx6O512AMpq2qhu7CQr1cHUcelHDfGMsfV+ywdGbcN7aB6BltFUNnQyflTyYfnReL2jUWc4Uu8/X/IYf+axIdEjISEWhyOGjIwEYmNtTJ8+mYyMBKqrSyku3sN3vnMbAF6vF7e7hXHjCvjLXx7F4XBQXV3N7NmzychIwGIx93oe+HG7W6ivr2XOHK3Vv2jRAkpKSk7ofxdOA5AINId89wkhrFJKrxAiB7gXuAz4Un+ENTZ2GK5gNOTx6Ek06gzh07u93YXP72fj7mpMJpgyOpmODvdRj3d7vEfZ44XsjVC6iI1F1TjsFnIzE2hvd0Xd9Vb3yOBpbe2io8NNbW0rXV0eWlq6qK1tJT09hxkzZvPDH/4Ev9/PSy89T1xcMj/5yXX85z/LcDic3HffPYfvG78fampaPnPeM8/8nbi4ZFJSUlm/fhuFhWNYt24TCQnGX4NjGZRwGoAWILRks5Qy+PRdCaQD/wWyAYcQokhK+UwY9VEMUw5VtdHh8jI2N5EER++un35h64Ssbfgr57JpTy2ZaU7jlFQMC0455TQ2b97Ibbd9nc7ODs4//zwcDifnnXch3/jG9SQkJJCSknY4w+fMmbP43vfu5I9/fOKI80477UwcDic/+9mv+PWv78HhcOJwOEhIOLE9t7BlAxVCXAFcJKW8Xh8DuEdKeUEvx10PTOprEFhlA9WIRp0hfHq//FExH2wsx+vzc87cfOx9RPG8Xv2XPmWOaf8SBypbmViQzKRRyVE3BqDukRNHNOh8rGyg4ewBvAacK4T4FDABNwghrgbipZRPhrFcxQjiUFUrLo+PCflJfVb+/WXK6BTqmrrYU9JESvwgehQKRYQTNgMgpfQDt/TYXNTLcc+ESwfF8Mbl9rG3vBmrxcS43ETD5FosZmZPTGfFtiq2FtfT6fL2GVIaygPr7h+0Dj+Y/+NBy1Ao+kIFPCuilg83leH2+Bmbm2j4BK6keDtTxqTS5faxdMUBQ2UrFJFCZE5PUyj6wOXx8fbaEsNb/6FMGZPKgYpm3t9YyqJp2YzO7v8AXau7leKmPZS0HKKus5YuXxcWk5X0uHTGJI1lWtp07NbeFxVXKE4UygAoopKV2ypp6/QwMT8pbDl8LBYzM8alsXpnNc++U8RPr517zKRx/oCfj0uXs7T4FfY1FRNAi1uItcTisDnx+D0cbDnAwZYDrCpfwemjzmRG+iyVlloxZCgDoIg6/P4A/1tfgtVipjAnPK3/IBnJcZw8JYs1u6pZvrmcs+fkf+aYmo4aXpL/5rldT3GgWUtOl+nIYnr6DMYmjScxJvFwJd/uaWdH3TbWVa3hvUPvcrD5AJ8fezFWs3oUFSceddcpoo5Ne2qpberi9Fm5xMaEP3nbVWdPYOu+epau2M+CKVnEx9mo7ajl04oVvLznRd4/9D98AR+xlliuEldjNVnJdub02rJ32pwsyFnIlLRpvLX/dfY27eGVvf/h8glXYjPbwv5bFIpQlAFQRBWBQIC315ZgAj43bxRrd1WHrSyv302tu4xtjaUki818JHdy5UuP02ItZm/TnsPHzciYxZfF1Vwx8UukxKb2KwooISaBL0686rARePvAW1w09hLlDlKcUJQBUEQVe0qbOFDZwkkT0skxeKZuIBCgyVvDkvL7ke0bqHId0Pz4Uj/ADrSBw+rkzFFnszD3FM4dfT5T0weWwdFqtvL5sRfz8p4X2dNYxOrKdBblLjbs9ygUfaEMgCKqeHddKQDnLygwVG6Lt54drato8FQCYDfHMTFhNimWHNJickmz5WBxpXHwgI04uxVz1wbcPjdv7X+dt/a/PuByrWYrl4y/nOd2Pc3qilUUJo4hNz7PqJ+lUBwTZQAUUUNFXTtbiusYl5fIhPxkw+SWdhaxrXUFAfxkxhTwldy7EfHzSIqPp73ddcSxlqZqqho6oS0bbZmLwRNnjePCMV/gRbmEtw68wXVTbjRErkLRF2oimCJqeHddCQDnzx9tmMzi9i1sbf0Yq8nGvKTzmZ98AVMTFmE19T4gO7UwFfBD/ST8PuMen1EJBczLXkCzq4k1lZ8aJlehOBbKACiigqY2F6t3VpGVEsdJE9INkVnWuYei9rXEmuM5JeVSsux9GxZnnA1SDoA3jvpDuYboEWRRzmISYxLZUL0O2fCZrCkKheEoF5AiKvhgYxleX4Dz5hccczJWf2nx1LOt9WNsphgWJF9IvDW5/yen7IOWPBpKcknKqSUmzhhXkM1i4+yCc3mt+BWue/srXD7hyqMeq3IFKYxA9QAUEU+X28tHm8tJcNhYNC170PJ8AR+bWz7Ej59ZiWeTYE05PgFmH6RLAn4ztcXGDkaPTRpPfvwo9jfvo7y1zFDZCkVPlAFQRDwrtlbS3uUlN83Jf9ccYumK/YdfA6G4fTOtvgZGx00hyz7ACjy+grikVlpr02hvMG42sslk4rT8MwBYUf4R4VqvQ6EAZQAUEY7X5+d/60uwmE0U5gx+taROXyv7OrZgNzuY7Dx54IJMkDXhIBCgZm8hAf+gVTtMbnwe45LGU9ZWxoGWgRk5haI/KAOgiGjW7a6mvsVFQVa8IQu+7G5bhx8fk5zzsQ4y9UJsYjtJOTW42h00VWQNWrdQFuedBsDKso9VL0ARNpQBUEQs/kCAt9eUYDaZGJebNGh55V3FVLiKSbKmkx870QANIWNcKWarl9r9o/C6jYupyHBkMiV1KjWdNexplH2foFAMABUFpIgYevr0qxo6KK9rJz/DiSN28Lfq2zV/B2Cic84xc+4E1w2OsVlxe7zHlGmN8ZI+poyavYXUHRhFtjBu8ZiFuYvZ3bCLtVWrmZgiVJ4gheGoHoAiYikuawZgfN7gW//VrkOsb/4fidY0MmOMm0gGkJJXTYyjg6byTLpaHcbJjU1hYoqgpqOaQy0HDZOrUARRBkARkdS3dNHQ6iIrJY5E5+AXZn+/7gUC+BnvOMmwlnRRSSNFJY3IsgbcyTsAEwe352Oky35+tjZQvbZqtXFCFQodZQAUEcnh1n/+4Fv/Hb4WVje+QZoth2z7mEHL6xVHPTiroSuVluo0w8RmObMpTBxDaWsJlW0VhslVKEAZAEUE0tLuprqxk5QEO2mJg183d1XDMtyBLk5P+xJmUxhv+fTdYPJRs3c0Po9xC9UEewHrqtYYJlOhAGUAFBFIcbnW+p9gQOvfH/DzUcN/sJliOTXl8kHLOya2Tkjdi88TQ02xceMMoxIKyHbmsLdpD/Wd9YbJVShUFJAiomjr9FBe206Cw0ZWStxxnx+M4AlS5y6n1l1KfuxEPqhfYpSaRyf5IHbXKJorM0nMrsOZ0jJokSaTiQXZJ7Ns32usr17L+YUXGqCoQqF6AIoIY29ZEwFg4qhkQwZrSzp3A1AQO3nQsvqFKUD2pP1AgGo5Br/PmAHn8ckTSY1NZVf9DlrdgzcqCgUoA6CIINo6PZTVaK3/3LTBh1O6/V1UuQ4Qb0kmxWbsTN1jEZfYTsqoKtwdcdQfMmZ1L5PJxLzsBfgDfjZUrzdEpkKhXECKiGFPqdb6Fwa1/iu69uHHz6i4SSd8ElXGmFJaa1KpP5RLYmY99vhOQAsd7Y1JBX1nJJ2SOo1V5SvZVruFhq56UmONizZSjExUD0AREVQ3dFCm+/5zDGj9A5S7igHIs483RF5/KSppZE9FPd7U7RAwc2BbAUWHGo9a+fcXi9nCvOz5ePwe/rbtrwZpqxjJKAOgiAje+PQgYFzrv8PXSqOninRbHrEW56DlDQhnLcRXQFcKtBizbsD09JnEWeP4+/YnaHO3GiJTMXJRBkAx5FQ3dLB6Z5Whrf+KLq31nxt7Ylv/nyF9N5g9UCfAO/gZzTGWGOZkzqPZ1cQzO58yQEHFSEYZAMWQ88anBwkEjGv9g5b504yZnHDN/O0vVjek7YGAFRomGCJyVuZs4m0J/GXLH+nydhkiUzEyUQZAMaRU1rezemcVeRlOw1r/Ld56Wn0NZNoLsJnthsgcFImlYGuDllHgjh+0uFhrLDdOu5nazhr+VfSCAQoqRirKACiGlFc+3k8gAJefOtbQ1j9Art2YFvegMQUgvQgwaa4gA/jGzNuItcTy582P4fF5DJGpGHkoA6AYMorLm9m0p5bxeUnMmpBuiMxAIEBFVzFWk23g6/2GA0ctxNVDRyZ0DD58M9ORyTVTrqWk9RCvFb9sgIKKkYgyAIohIRAI8PJyraX+xTPGGdb6b/RU0+lvI9s+Bospgqa5mIC0IiAA9QIMSBl926w7sZqt/GHTI/iNXJRYMWJQBkAxJGzbV8+esmZmjU9n4qhkw+SWu/YCJz72v1/EtkB8JbiSoGPwPZ5RCQV8ceJV7GmULC1+xQAFFSMNZQAUJxy/P8DLH+/DBFx++ljD5HoDHiq79hNjiiMtxpgUDIaToi972TjOEHHfnftDbGYbv117nxoLUBw3ygAoTjird1ZRXtvOounZ5GcMPiomyO62tbgDXeTGjgtv3v/BYG8FRw10pUJn3+kf+mJ0YiFfm3I9B1sOsKToeQMUVIwkIvQpUQxXPF4fS1fsx2oxc+li41r/AOua3gYgb6gnf/XF4V6AMb//rrk/wGF18PCG39Hp7TREpmJkEDYDIIQwCyH+KoRYLYT4SAgxvsf+K4QQ64UQ64QQXw+XHorI4sNN5dS3uDh7Th5pSYNf7SuIy9/JlpYPcZgTSLZmGiY3LMQ1QmyDFhHkShi0uCxHFjfPuJWq9kr+sf1JAxRUjBTC2QO4FIiVUi4E7gYeDu4QQliA3wLnAAuB7wshjIkDVEQMS1fsP/xa8m4RLy0v5tVP9hNnt/L5hYWGlrWt5RNc/k5yY8ef8MyfA8LgXsDts+4kyZ7MHzY9TENngyEyFcOfcBqAxcA7AFLKNcDc4A4ppQ+YLKVsBtLQguTawqiLIgIoLm/G4/VTmB3P+xtKjzAQS1fsH5TsbvdPhEz+6gtHLcS0Qlu2ITmCkmNT+Nbs79LkauIXH/3CAAUVI4FwBkonAs0h331CCKuU0gsgpfQKIS4H/gy8BRwzhCElxYHVatxC20EyMgbfBT/RRIvOTmd3GoZOl5f9FS3E2S1MG5+B1WJc26PN28SOtlUUOCaRFpdhmFyAGFv4HhFfSjm+6klY2gtwOo+v/dPbPfDjs77PEvksf17/Z74595tMyZhilKonjGi5t0OJRp2DhNMAtAChV8YcrPyDSClfFUIsBZ4BrgWePpqwxsYOwxXMyEigtja6UupGk87t7a7Dn3cdasTnDzAxPxlXlwfXMc47Xj6pfwtfwMPcxPNwe4xLjhZjs+L2ePs+cKA4SsE0Hl9jHm1t2zgez9X33/pRr9unp86kuKGY29+4gxe/8Fp0uMN0ouneDhINOh/LQIXTBbQKuBBACHEysD24QwiRKIT4WAhhl1L6gXZATWUcprR1ethX3kx8nJVRWcaFfQZZ1/RfTJhYkBxli6VbvJBQAV4H7fXJhogcmzSesSlj+aj0Q9479I4hMhXDl3AagNeALiHEp8DvgbuEEFcLIb4hpWwB/gl8IoRYiTYxXqU1HKYUHWokEIBJo1MwG9wirXOXs7djExOdc0/our+GkVQCQGO5MbqbTCbOG3ceFpOFn6/6MW6f2xC5iuFJ2FxAesv+lh6bi0L2PwmomLVhTmOri4r6DtKSYslJNSbdcyjBwd8FyZ83XPYJwd4KsY201yfj7rQTEzd451imM5Mbpn2dv29/gr9vf4LbZt1hgKKK4YiaCKYIK8F1cGdOyDDcHx0IBFjT9BZWUwyzk842VPYJJbEEMNFUYdz8he/P+xEp9hQe3vA7ajpqDJOrGF4oA6AIGw0tXdQ2dZGeFEtWGFr/JV1FVLkOMDPxdByW6I3EIL4Ks9VLc2UGRiX1TIlN5Qfzf0Kru4Xfrv2VMUIVww5lABRhQ5Y2AdpSj+FgXdN/gSh2/wQx+0nMqsPnjqGtfvD5gYJcN/VGJqVO5p+7n2Nb7RbD5CqGD8oAKMJCcVnz4da/kSkfgvgCXtY2vY3TksS0+FMMl3+iSc7V3DTNlcbNY7Cardy3+HcECPCTlT8kEDBgEQLFsEIZAEVYWLZSm9krCpLDIn9760pavHXMSzoPq9kWljJOJLEJHdgT2mirT8HrMu73nJZ/BheOuYi1latZVvyqYXIVw4MIWjJJMVwoLmtm58FGrfWfaHzrH2BFg7YAyqmpV4RF/lCQnFNL9Z54mqsySBtdMShZD6y7//Dn3Pg8LCYL3/34W+ys24HNYuMH8388WHUVwwDVA1AYzlurDwLh8/3XuyvZ0bqKMXHTGRVnzCLrkUBiVh0ms5+migyM9NYk25OZmzWfVncLG2vWGydYEfUoA6AwlLKaNrbuq2d8flJYfP8AqxqXEsA/rFr/ABabj4SMejydcXQ2GRvVND/nZOKscayrWqvWDFAcRrmAFIMmNJPnxj21AKQm2I92+KDwBbysalxKrNnJvOTzwlLGUJKUW0tLdQZNlZk4UozLMWO32Dk5ZxHLSz9gbeVqw+QqohtlABSG0d7loaK2nQSHjayUuLCU8VTpT2j0VDM6bgrv1j4TljKGEkdyC7a4LlprUvFNPIjF6jNM9syMk9hYvZ7NNRspay0lP2GUYbIV0YlyASkMY195CwFgQl5S2LJQHurcDUBB7OSwyB9qTCZIyqkh4LfQUp1mqGyr2copeafhC/h4YP39fZ+gGPYoA6AwBJfHR0lNGw67ldwMZ1jKqHYdosZdQrI1kyTb8FpArqik8fCrzl8MBKg+mGp4OZNTp5Ael8GLRUvY0yANl6+ILvplAIQQ3xdCZIdbGUX0crCqFb8/wNjcRMMzfgZ5r+55AMY6ZoRFfsRgdWkrhrmS6WozNoWG2WRmce6pBAjw6KaHDJWtiD762wNwAB8JId4SQlwphIj+mTcKw/D5AxysbMFqMVEQhnz/AC2eej5tfB2HOZFs+5iwlBFRJJYB0Fxh7ApnAOOSJzA5dSqv7n2JA82DW4pTEd30ywBIKX8hpZwE/AY4E9gqhPiTEGJWOJVTRAflde24PH5GZyUYutRjKMvr/4034GasYwZm0wjwXDprwOKiuSodv9/YHpXJZOKuOd/DH/Dzx02/N1S2IrrodxSQEMIJjAHGoq3e1QA8JoT4VErZ+/p0imFPIBBgf0ULAGNyEsNSRpevg+UNLxJvSWFU3MSwlBFxmAKQUI6/aSxttakkZtUbKr6oYTcp9lSWFD2Pw+og0Z70mWPUbOHhT3/HAF4AioEzgPuklNOklD8HPgd8M3zqKSIdWdJES7ub3DQHjtjwRBWvalxKh6+FM9OuwmIaQd5H3Q3UFAY3kNlkZkHOQvwBP+ur1xouXxEd9Lcv/SEwXkp5k5RyJYAQIkZK6QKmhE07RcTz3oZSAMbmhqf17/Z38b+6Z7GZYjkj7aqwlBGxxLQTl9RCR2MS7k7jJ9ZNTp1CYkwS22q30u5pN1y+IvLprwG4WUp5+A4RQpiBjQBSyqpwKKaIfKobO9iyt47k+BhSwjTz98P6f9HoqeastC+TYDUuV360kJRTC5gMTRMdxGK2MC97Pr6Ajy01mwyXr4h8jmkAhBAfCiH8wAIhhD/4AroAFUQ8wvlgQxkBYFxuYlgmfrV7m3mn5ikclkTOz7zRcPnRQGJWPWarl6aKTMMHgwGmpk0n1hLLltrNePwew+UrIptjOm2llGcBCCEek1J+68SopIgGOrq8rNheSUqCnZy0/k38er36L30ec3HWrYc/v137Dzr8rXwx+zs4LeFxMUU6Zouf5JwaGkpzaa1JIym7zlD5MZYYZmTMYl3VGnbX72JGxkxD5Ssim2MaACHEF6SUbwKbhBDX9twvpXwubJopIpqV2ytxuX18YeFoPF6DFrINod5dyYf1/ybVlsOZI83334Pk/GoaSnNoLM023AAAnJQ5hw3V69hQvY7p6TPClsZDEXn0NQYwT38/Ay3+P/R1Rti0UkQ0gUCA5ZvLsVrMnD4rLyxlvF79Z7wBN5dk3Y7NHJ7xhWghJs5FfHojXa3xdDYbP9EuISaBSSmTaeiq52DLAcPlKyKXvlxA9+jvNwS3CSESgVFSyp1h1k0RoRQdaqS6oYOFU7OJjzM+LLO0U7Km6S3yYyeyIPlCw+VHIyn5VbTVpdJYlk1cUrHh8udkzWNXw042VK1jTNJYw+UrIpP+zgO4SQjxjBAiA9gFvCyEULNERijLN5cDcObs8LT+X636AwECXJ79rZEx67cfOFJaiHF20FKTisfANYODZDmzGZVQwKHWg9R11houXxGZ9Pfpug34EfAVYBkwHbg8XEopIpemNheb99YxKjOecWGI/S9qW8fOtlVMcs5navwiw+VHG8EMobK0EbdjPwTM7NuVHJay5mTOBWBT9YawyFdEHv1uXkkpK4ELgbeklF4gPCt+KCKaT7ZW4PMHOOOkPMMHCwOBAK9UPQrAFTnfVoORPUkoB4sLmgvweS2Gix+bPJ4kezK76nfS4ekwXL4i8uivAdgphHgTLQ/Q+0KIFwG1uvQIw+f38/GWCuwxFk6ekmW4/ErXPg517mJe0vmMjlMTzD+D2Q/JB8Fvo6nM+OtvNpmZnTkHb8DLtrothstXRB79NQA3Ag8AJ0sp3cALwE1h00oRkWwrrqex1cWiqdnE2Y3N++MP+ChqW4/FZOXS7P8zVPawIrEEzB4aSrPx+4zvIU1Ln0GMOYYtNZtw+9yGy1dEFv19iuPR/P6nCyGCd91JwC/DopUiIjk8+HuS8YO/hzp30+FvoTBuGqsb3zBc/rDB4oWkEnyN42iuzCQlv9pQ8XaLnenpM9hYs4E39i3liolfMlS+IrLobw/gJbTYfwtgCnkpRgg1jR3sONBAaoKdDbKGpSv2H34NFq/fzd72jVhNNiY4Zxug7TAn6SAms5/6klwCYUgPcVKWNhj8xNY/EwgEDJeviBz62wPIllKeG1ZNFBHNx1sqACjMTjBc9sHOnbgDXUx0zsVuVrEFfWJ1k5RbQ1NZNs1V6STnGhu2mWxPZnzyBLbUbmZ91Trm5ywwVL4icuhvD2CzEGKYL8SqOBpen59VO6qwWc3kpBu7Rq3X72ZfxzZsJjtj4qYZKns4kza6HJPZT93B/LD0AuZkaUkAntz2uOGyFZFDf3sA09CMQDVaJlATEJBSqimDI4Ad+xtoaXczJicBi9nYiVkHO3fiCXQhnPNGfMqH48Fm95CcV01jaQ5NFcaPBeTHj2Ja+gze3L+M0tYSRiUUGCpfERn092m+DC0EdCHdeYDODJNOighjxTbN/VOQaWweGq31vxWbyU5h3FRDZY8EtF6Aj/qDeYZHBJlMJr4x41b8AT9Pbf+bobIVkUN/F4U/BJwCfAOoBU7XtymGOS3tbrbtq6cgM56keGNb6Frr38VYxwzV+h8A1hgvqaOq8LpjaCo3fl7AZRO+SHpcBi/sfpY2T5vh8hVDT39zAf0WbRbw5WhuoxuEEA+HUzFFZLB6ZxU+f4DFM3IMlesNeA77/lXrf+CkFlRgtnipP5SH32use85usXPDtK/T7GriP/JfhspWRAb9vWPOA74GdEkpW4BzgQvCppUiIggEAqzcVonVYuLkqdmGyi7v2osn0EVh3FTV+h8EFpuP1IJKfB4bDWXG/kcA1029iRhzDH/b9hf8AePXfVAMLf01AD3/eXsv2xTDjINVrZTXtTNrQoahaZ8DgQAHOrZjwqxSPhhAyqgqzFYPDSW5+DzG5gjKdGRy+cQr2ddUzIcl7xkqWzH09DcK6D/Ai0CKEOLbwLXAkmOdoC8c/zgwE3ABX5dSFofs/wrwbcAHbANuk1IqozLEhE7s2ravHgCLCUMmfAWpdZfS5msiL3YCsZb+LSepOJKiksYjNyTth3rB3h0pzFncamhZN8+4lX8X/ZMntj7OOaPPM1S2Ymjpbw/gLeANoA44FfiZlPL+Ps65FIiVUi4E7gYOjxkIIeKA+4AzpZSLgCTgC8enuiKc+P0BKuraibGZyUgxdnLWgc7tAIyJm26o3BFN0iEtU2hTIR6XsXmapqfP4JTcU/m4bDm76tU6UMOJYxoAIUSmEOIT4GPgdsALnAXcLoRI6kP2YuAdACnlGmBuyD4XsEhKGcw5a0WbX6CIEGqbOnF7/eSlOzEbmJa51dtIrbuMVFs2ybYMw+SOeMw+SNkHAStVxcZHBN06S0vQ94dNjxguWzF09NVU+A2wEjhbSukBEELY0JLAPQZcf4xzE4HmkO8+IYRVSunVXT3Vurw70JLNHdPBmJLiwGo1Pgd6RobxqQ3CTTh1djq1Adkq3f0zoSDl8LbBEmOzUtKmtSAnJMwixmZsSzUcRIOOQQJpFXiaxlBzII1GcxEmm+uI/dPHpx+XvND77Or0K3lg469ZWvwKvzv/fsanjjdE52OVGS1Eo85B+rq7F0kpJ4dukFJ69OUgt/RxbgsQemXM+kIywOExggeAicAVUspjZp1qbDR+gYqMjARqa431l4abcOvc3u7C6/VTVtOGM9aK3WKivd3V94l94HTaaXd1UNIhiTPHk24Zhdvj7fvEISTGZo14HT9DajHUTMdTUwiZR7prjvd/7Hmf3THzO9z8v+v5xfv38ciZfxyspp9BPY/h4VgGqq8xgF7dMnpl3deA7Sq0uQMIIU4GtvfY/wQQC1wa4gpSRACVDR34/AHyM+INXZWrwrUPH14K4iZjUmv9hoeEcohph5Z88Bg7dvOFsZcwPnkCL8ollLeWGSpbMTT09RQeq1XeV57Y14AuIcSnwO+Bu4QQVwshviGEmI22oMx04EMhxEdCiMv6rbUirJTVarM+8zKMjdAp7SwCID92oqFyFSGYAljS9wNmaDDWTWMxW7hz9nfw+D08vuUPhspWDA19uYCmCiF6i/8zAcecGqr7+W/psbko5LNqAkYgXW4ftU1dJMfHGBr7X9axlyZvDRkxo4izGJtTSHEk5sQqfHWjoTUPUvZrPYIB8MC6zwb6+fw+EmMSeXrH37lzznfJchg/4Kw4cfRVCU9ES/rW83UGIMKqmWJIqKjTKov8DGMr6Y9qXwVgVOwkQ+UqPovJBKTtBUzQMMFQ2RazhQU5C/EGvPxho8oGE+0cswegEr6NPMpq2zABeenGuX+8fg8ra5dhM8WSZR9tmFzFMXDUgL0Z2rLB7RxwL6A3pqXNYF3lGp7d+RS3zrqD/IRRhslWnFiUG0ZxmKqGDpra3GQkx2GPMS7kdlvrJ7R6G8mPnYDFZHwor6IXTEBKsfah0dhlOyxmC4tyF+P2u3lkwwOGylacWJQBUBxmzc4qAPINHvxd2fgaAKPilNfwhOKsgZhWaM01PCJoctpUJiRP5F9FL7C/qbjvExQRiTIACkBL0LZmZzUWs4nsNOOWfWz0VLOz9VPGOqeTaE0zTK6iH5jQZgdjNrwXYDaZ+eH8n+AL+Hhw/W8Nla04cUTPNEdFWNlf2UJNUyd56U6sFuPaBasb3yCAn9MzL6e2vcowuYp+El+pDQS35OPpqsMW6zZM9O76XWTGZfLK3v8Qb4snw5F5xP4fzP+xYWUpwoPqASgAWLNTW1PWSPdPIBBgVeMybCY7i9I+b5hcxXEQ0gtoKDF2UR+TycTi/NMB+KTsI0NlK04MygAo8Pn9rN9dTXycjYxk43zFe9o3UusuZXbSOTis0ZsvJepJqABrB00VWXjdxs3tABiTOJZRCQUcaNlPSYsKGow2lAFQsPtgIy0dHuZNzsRsNi71w6eNSwFYnHKpYTIVA8AUgJQDBPzh6QWcnn8mAB+XLScQ6CtBgCKSUAZAwWrd/bNwinFLCnb4WtnQ/D7pMflMcM4xTK5igCSUYY1x01SeZfiqYdnOHETKZKo7qpCNRX2foIgYlAEY4bg8PjbtrSU9KZZxeYmGyd3Q/C6eQBenpFyCWSV+G3rMflIKKvH7LDSWG5++YXHeaZhNZlaWf4zP7zNcviI8qCdzhLNlbx0ut48FU7IMzfy5smEpJswsSrnYMJmKwZGcW4PZ6qWxNBu/z7j/GiAlNoWZGbNocjWxtW6LobIV4UMZgBFOcPLXyVONc/+Ude7hYOcOpiYsIsWmkoVFCharj+S8anyeGJqrjF+N7eScU7CZY1hdsQq3b/BrSCjCjzIAI5jWDjc7DjRQkBlvaO6f4MzfU1MuN0ymwhhS86swmf00lORg9Hit0+ZkfvYCOr0drK9aZ6xwRVhQBmAEs6GoBp8/YGjr3+3vYk3jmyRa05meeKphchXGYLV7SMquxdMZR2ttquHy52TNw2F1sqF6HdUd1YbLVxiLMgAjlKUr9vPWai1uu7G1i6Ur9rN0RW9LPxwfm5s/oMPfyqKUi7GajI05VxhDakElEKDhUK7hvYAYSwyLchfj8Xt4SKWIiHiUARihtHV6aGh1kZ4US5zduIwgKxq1vP8q9j/yKCpppKikkf11leCsoqs1HimNrwKmp88gxZ7KC7ueYV/TXsPlK4xDGYARSmmNtuxjQaZxC79Uuw6xp30jwjmPTHuBYXIVYSDlgPZucJI40NJFn5p/Or6Aj1+v+aXh8hXGoQzACMTn91Na04bVYiLHwMyfKxq01v+pqWrwN+KJbYa4euhMp6vVuHsgyITkiczJmseb+5exQQ0IRyzKAIxAduxvoMvtIz8jHotBmT9d/k5WNr5GgiWFkxLPMkSmIswka2M+9YdyDRdtMpm4Z+GvAPjl6p+rFBERijIAI5CV2yoBY90/qxtfp8PXwulpV2Iz2w2TqwgjjjqIaaG1Jg13p/H/2cm5iziv8ALWVH7Ke4feMVy+YvAoAzDCaGl3s6W4jkSHjaT4GENk+gN+3q/7J1aTjTNSrzJEpuIEYAJS9gMmw5PEBfnJyfdiNpm5b829KkVEBKIMwAhj9c4qfP4ABVkJhqV+2N66ghp3CfOTLyDRplb9iiriq7DFdtFcmYnXbfz6UJNSJ/NlcQ1FDbt5US4xXL5icCgDMILwBwJ8vKUCq8Vk6MIv79e9AMDZadcYJlNxgjAFSC2oJOA301hm3ITAUH4w/8fEWmL53bpf0+ntDEsZioGhDMAIYueBBqoaOlgwOYsYmzEpgUs7i5Dt65nknK8WfY9SknJqsdg8NJZl4/caXyXkxufxjRm3Udlewd+2/dVw+YqBo9YEHkG8v6EMgLPn5rNlb50hMt+tfRaAc9K/yuvVfznqcTE2datFKnvK6yHhIDRMYM+uBEg+CMCkghTDyrhj9rd5ftfT/GHTI3x1yrWkxipXYSSgegAjhKqGDrbvr2d8fhKF2cbk/a/sOsD65nfIj53ItITFhshUDBFJh8DkhaZCCBibKhogyZ7Mt+d8nxZ3M49ufNhw+YqBoQzACOGDjVrr/5w5+YbJfKvmCQIEuCjzFrXoS7Rj8UBiGXjjoDU8EUE3Tr+ZUQkFPLX9SUpbS8JShuL4UE/tCKCjy8vK7ZWkJNiZPdGYPPClnZL1ze8yKlYwK/FMQ2QqhpjkA4AfmsZCGOZt2S127p7/U9x+N79de5/xBSiOG+WYHQGs3F6Jy+3jCwtHYzVg5m8gEODlqkcIEODy7G8ZupKYYgixdUFCJbTmQUcG4B2UuAfW3f+ZbYFAgIy4TF7a829unXUH09KnD6oMxeBQPYBhjsfr5911JcTYzJw205gp/ztaV7K7bS1T4xcxNWGRITIVEYKeHiIcSeJASxFxWv4ZAPxq9c/DUoai/ygDMMxZtb2SxlYXZ56UR4Jj8DN/Xf5OllT8BjMWrsi5ywANFRGFvQ0cNdCVSkezcalCQilMHENBwmiWl37AJ2UfhaUMRf9QBmAY4/X5+e+aQ1gtZs6fb0x65ter/0K9p4LPZVxLfuwEQ2QqIowUrRfQEIYkcRDsBWjjRr9afQ/+gD8s5Sj6RhmAYcyandXUNXdx+sxckuIHn+zrUOdu3q97gYyYUXwh85sGaKiISGIbwd5IW10qXW1xYSki25nN5RO+yNbazSwrfjUsZSj6Rg0CD1P8/gBvrT6IxWzigpMH3/r3Bbz88eAdBPAzzjGDd2qfNkBLRURiAlL3QeVcavcVMGqmDEsxd8//GW/sW8b9a3/J58deHJYyFMdGGYBhytrd1VQ3djI6K55PtlYMWt7bNU/R4q0jP3Yi6THGzSVQRCiOWuKSW2ivT6GjKQFHcqvhRRQmjeH6qTfxt+1/5bmdT/Gj7O8bXobi2CgDMAzx+vy89sl+LGYTE/KT+n3e0VI5NLir+LTpdWLN8UyJX2iUmopIxgSZ40o4tHEatcUFFMzZidHRvg+sux+bOYYYcwy/WnMPLYEGvK4jj/nB/B8bW6jiCNQYwDDk4y0V1DV3ceZJeThibYOS5fa72NzyAQAnJZ5FjDnWCBUVUUBcUhvxGfV0tiTQVmdcXqBQHDYH87NPptPbyarSVWEpQ3F0VA9gmLB0hRa54fX5eX9jGRazCYt5cE22QCDAttaP6fS3MdE5h7SY8KQIUEQuGWNLaatLpXZfAfFpTZjMxk8Rnp01l821m1hTtoapyTOJjwlP+Knis6gewDBjX0ULbo+f8XlJ2GMGl/K5pGs3Va4DpNpymOCYbZCGimjC7uwiObcGd0ccDWFaLyDGEsOi3MV4/B5WV64MSxmK3glbD0AIYQYeB2YCLuDrUsriHsc4gPeAm6SUReHSZaTg8vjYV95MjM3MuNzBZfxs8dazs/VTbCY7JyWehUklexuxZIwtpbUmlfoD+SRm1WGzewwvY3r6DDbXbmBb7VbmZM1T6aJPEOF8qi8FYqWUC4G7gSNywAoh5gKfAOPCqMOIQpY04fUFmJifjNU68L/WG/Cwqfl9/PiYlXgGcRbVJR/JWGxeMsaV4PdZqNlbGJYyzCYzZ485mwABVpR9HJYyFJ8lnGMAi4F3AKSUa/QKPxQ7cBnwfH+EpaQ4sFqNWcUqlIyMBMNlhpvedPYG4FBVK4nOGKaOS8c8AP9/cNGWbY0f0+ZrYrxzJqPixw9a357yo4lo1BmM0dvp7J486BjfQmt1O601afg6GknMMD4sVDgEoxJHsbdpDw3eGkYljYqK5zMadDwa4by7E4HmkO8+IYRVSukFkFKuAhCif8sINjZ2GK5gRkYCtbXG38jh5Gg6r99VRQCYXJBMZ6d7QLLdHi9lnXso6SgiyZrOBMc83J7BZYQMEmOzGibrRBGNOoNxere3HxmTmTF+H+3rp3Nwaz5j5m/DbDE2hYPTaWdxzun8q+UF3i3+H18W10T88xkNdcixDFQ4XUAtQGjJ5mDlrzCWbfvqqW3qIiM5lsyUgU/db/M2sb1tBVaTjdmJ52AxGd/jUkQvsQkdpBZU4umMpabYmNxSPclLyGd88gTK28rY11zc9wmKQRHOHsAq4CLgP0KIk4HtYSxrxOL1+Xnxw70ATC1MHXBufl/Ay+aWD/EFvMxOPBuntf8TyBTDk6KSxs9utDZjdybTVJ5NfHoj8WnNnz1mkJyadzr7mor5pOwjvH4vVnN0uuGigXD2AF4DuoQQnwK/B+4SQlwthPhGGMsccXy8pYLK+g5GZ8WT6Bx4uud3ap+m2VtLXuwEcmON8/srhhlmPzlTisHkp2r3OHwe4yvntLh0pqXPoKGrnheLlhguX9FN2EyrlNIP3NJj82dCPaWUZ4RLh+FOe5eHpSv2E2e3MKlg4DM1yzr38GbNE9jNDqbGn2KghorhSGxCBxljyqjdX0Dl7rHkTd9jeJqIU3IXs7thJ79b/2sum/BFHDaHsQUoADURLKpZtvIA7V1evrCwcMCTvrx+D0+X/QxfwMuMhNOJMQ8+bbRi+JM6ugJHcjNtdak0lBi/bkB8TAJzMudR1V7J37b1nqNKMXiUAYhSymrb+HBjOZkpcZwzd9SA5bxT+zSlXZJTUi4lyx6egT3F8MNkgtxpe7HaXdTuG0V7w+AmHvbG/OwFpMam8ofNv6ehq95w+QplAKKSQCDAkvf24A8EuPqcCdgGOOmrxlXCf2v/TpI1gytzvmuwlorhjjXGS960vWAKULFzAp7OwS85GordGstdc75Pq7uF32940FDZCg1lAKKQ9UU1FJU0MXNcGjPGpQ9IRiAQYEnFb/AG3FyV8z0cluidzKIYOuKS2siacAifx0bZtkn4vMaGDl8/7esUJBbyjx1Psq9pr6GyFcoARB1dLi8vfliM1WLiy+cMfE3ejc3/Y1fbaqbEL2RO0ucM1FAx3CkqaTziVe0rgqSDuNodlG+fQMBv3Iiw3WLn3oX34fV7uWfVTwyTq9BQAbZRxNIV+9lX0UJjq4sJ+Ums3lE1IDmdvjZerHwIqymGq3N/NOC5AwoFoC0hmb4bPHF0NGZRJceQPWm/YZFBnx97Eafknsr/Dr3D8pIPOLPgbGMEK1QPIJpo6/Sw+2AjcTGW41rpqyfLqv9Ms7eWCzNuIlMN/CqMwARkbyU2oY3mykzqD+YZJ9pk4leLf4sJEz9f9SO8fpVQwCiUAYgSAoEAOw404A8EmDImFatlYH/doc7dLK9/kcyYAs7LuMFgLRUjGrOPrrS1YO2g7sAoinbEHXYTDZZp6dP56pTrkI1FPLPj7wYoqwDlAooa1hfVUNPYSXaag9y0458U83r1XwgE/KxqXEYAP2Mc03m79h9h0FQxorG6IXcDlC2EmulgcYOzzhDRd8//GcuKX+O3637NReMvI8uRZYjckYzqAUQB7V0elry3B7PZxNzJWQP22Zd0FdHkrSHXPo6MmHyDtVQodGLaIWcjmAJQdRJ0GZNXKsORwU9OvocWdzP3rFKLxRuBMgBRwEvLi2np8CBGJZPgGFistcvfQVHbWqymGKbELzRYQ4WiB3GNkL0FAmaomIurPdYQsddOuYHZmXN4de9LfFy63BCZIxllACIcWdLIJ1sryc+IH9Qyj7va1uAJuBHOecRanAZqqFAcBWcNZO4AfwylWybj6Rr8RDGL2cKDpz+K2WTmh598hy5vlwGKjlzUGEAE0+ny8o+3dmMywXUXCLbvG9h0eNm2gfKuvSRZ0ymMm2KwlgrFMUgsB58db72gdOskRs/eicXm6/fpD6y7v9ftszJms6lmA49ufJC7F/zMKG1HHKoHEMEseX8Pdc1dfH7haMblDsyP6vV7WFKhPUTTE05Vi7srTjzJ+yHpAO52B3vXj6foYPOgo4NOyTuVhJhEHtv0CJurNxqo7MhC1QYRykZZw6rtVYzOTuDiU8YMWM57dc9T6drP6LgpJNsyDdRQoegnJiC9CBLKoStFGxgODG6WmN1i54LCC/EFfPzfB9+k09tpjK4jDGUAIpB/v7+Hv72xC7PZxNicBN789CBLV+w/bjl17nLerHmSBGsqk5zzw6CpQtFPTEDmdnDUQkcm1EyDwOBEFiQWcvP0W9jbtIf71/zCEDVHGsoARBhen5/1sha318/UwpQBR/34A36eKfs5nkAXX8r5HjaV518x1JgCkL0Z7E3Qmg/1YtAif3LyvYxLHs8T2x5nRdnHg9dxhKEMQITx7w/20tjqIi/dSWH2wDN0Lq//F3vaN3JS4tnMT7rAQA0VikFg9mkTxWxt0DSW+pKcQYlz2Bz8+ewnsZqtfPO9G6loKzdI0ZGBigKKIFZtr+TDTeUkOGzMHJ824AlfVa6DvFr1B+ItKVyT9xOV7E0RWVg8kLseyhZSWzwaq81DUs7AZgsHo4ROyzuDD0vf54JXzubL4hosZi0t9Q/mqwljx0L1ACKEPaVNPPeuxGG3Mn9S5oBz/fgCXp4u/RmegIuv5v2URGuqwZoqFAZg64Lc9ZitXiqLxtFWlzwocSdlzmFy6hQq2ytYXvqBMTqOAJQBiABKa9p47OVt+P0Bbrl0Ks4424BlvVr1GAc6tzM/6QJmJ6m0uYoIxt5G/swiTCY/5Tsm0NYw8AmKJpOJc0efT3pcBltqN7G1dotxeg5jlAEYYmqbOnnkP1vodHm56fOTmTYmbcCy1je9y3t1z5NtH8M1eWrxDEXk40hqI2/aXgIBM8XrxuJqixuwrBhLDJeMu5w4axzvH3qXPY3SQE2HJ8oADCF1zZ08/O8tNLe5+crZEzh5avaAZVV0FfNc+b3YzQ5uLXiYOEu8gZoqFOEjPr2JnEn78HmslG6dNKiUESmxKVwx4UtYzTbe2v86q8pXGKjp8EMNAg8R1Y0dPPivzTS0uLj4lELOnTdqwLJavQ08uP8mXP5OZieey/rmd6HZQGUVijDRPRu4EUumH2/NRPatn8j4BUVYYwa28Eu2M4dLx1/OK3v/w9f++2VevOhV5mUvME7pYYTqAQwB5XXt/Pafm2hocXHF6WO59NSxA5bV4WvlsQO30+5rZpxjFrmxA5elUAwllrRDWtoIT/ygF5gfnVjI58dcTKe3gytfv4QPS94zUNPhgykQGOR0vBNEbW2r4YpmZCRQW9tqtNjPEDqLt66pk/VFtXh8fqaNSWXscWb4dDrttLe7AHD5O3nswG0Ud2xmVOwkZiScFrEhnzE2K25PdC3lF406Q5Tr7fZqC8m05hOboA0SD7QnADAz8yRufvc6vAEvfzr7CS6fcKWBGp+4OmQwZGQkHLVSUD2AE0hJdSurd1Xj9fs5aUL6cVf+oXT5Onj80F0Ud2xmbtJ5zEg4NWIrf4Wi35iAzB0k5dTQ1RpPyaapgxoTOK/wAv5z0VIcVie3vHcTv1n7S3z+/mcjHe4oA3ACCAQC7DrYwJbieqwWM4umZjMqc+CDtE2eGh7cfyO729YwI+F0bhx1n8ryqRg+mAJkT9pPakEF7o44Dm2cSlfLwENET85dxLJL36YgsZDfb3yIL715GbUdtQYqHL2oWiPMdLq8bJC1FJe34Iy1cuqMHNKSBr46UknHHn6z71pKu4o4NeVybhn9EFbTwOcNKBSRiMkEmeNLyBh/CK8rhkObptJcmT5geVPTp/H+Fz/m/MILWVH2EWf95xT+u/9NAzWOTpQBCCMVde3c99wGKus7SEu0c+qMHOIHOMkrEAjwcf3L3LPjKho9VVyWdSdfzfuZqvwVw5q0gkryZ0hMZj+Vu8dTtacQv29g1VZybArPXvAv7ll4H41dDVz/ztVc//Y1VLZVGKx19KAMQJhYt7uaXz2rVf5jcxNZODWbGNvAohpavQ08fugu/llxHzaTnVsKHuaCzBuVz18xIohPb6Jw7g5inB00lWVzYN10OpoH5kI1mUzcftKdLL/qU07OWcR/D7zBwiVz+O3aX9HiGnmx0yoKyOARfK/Pz0vL9/HehlLsMRZuvHAy5bVtA5LlD/j4pOEVllX/mXZfM2m2XOalnYvVb8wC2yeSaIxMiUadYRjr7TdDw0RoKgQgOa+a9DFlx4wSOlYyOH/Azz93P8dv195HbWcNyfZkbpt1J9dOvYHU2P7NyI/2KCBlAAz88yrq2nnyjZ2UVLeRk+bg9sumk5vuPO7FXAKBAEXt63ip8mHKuvYQa3Yy1jGDMXHTscfYhufDHYFEo84w/PUuSCygqmgs7o44zBYfqQUVpIyqwmIdWHSP2+dmc81G1lWtweVzYTVbmZI6lVmZs8mIy+SHC46eViXaDYCaCWwAfn+ADzeV8dJH+/B4/Zw6I4evnDOB2Jjju7yBQADZvp43qv/K3o5NACxKuZjLsu7k44aXwqG6QhF1OJJbGTN/G00VmdQdyKfuwCgaSnJJzK4lJb8Ku7PruOTFWGJYkLOQWRknsb1uG5trNrKtbivb6raSFpuO3Wrn4nGXMjZ5fJh+0dChegCDtN6Hqlp57l3JgcoW4uNsXHf+JOaIjCOO6asH4A142NT8Ae/XvcDBzh0AzEg4jS9kfpNCx1QAXq/+CzD8W3eRRDTqDMNf70kFKYc/+7xmmsqyaSzPwuvSVr2LTWgjIauehIwGYuJcx62HP+Bnf/M+dtZtZ3/zPnwBrWcxJmks5xR8jsX5pzMvewHpcelR3wNQBmCAf15DSxdvrj7Ex1vKCQQgL93J1DEpx9Xqr3OX82nj67xf9wJd/nYAsmJGM8E5h2RbRq/nDPeHO5KIRp1hhOodMEF7FjTnQ2cah+NbbO0kZ7ThSGkhNqENW5yL44md6PJ2MT5lAv87+A4fly2n3dM9njcueTynFZ7K1KRZTE+fwaS0KcRZB57NNFwoA3AUBmIAaps6eXddCZ9srcDrC5CV6uBrn5vIntKmfp3v8neypflDVjUuo6h9HQAWk41RsYIxcdNwWpOOef6IfLiHiGjUGZTe+GzQlqUtPt+ZCv7uUGmz1UtsQhtxie3EJrZhj+/AFntsoxAcSHb73KyrWsOaik9ZV7WGDVXrafN01x8mTKTGppHpyCTDkUVmXCYZjkx+ccqvB/+bBoEyAEehvwago8vDluI6Vm6rpKikSTs3OZYvLCpk4dRsrBbzMd08Hb4WtresYFPLh+xsXYU7oPkoJzjncErKJVR1HcRq7l88/4h/uE8g0agzKL2PIGBidHI+nc3xdLXG09nixNPZo5Vu9kBMK9hbD79PHGfBbPEfU7Q/4KedZg7WlVDTWU1NRw21HTW4/e4jjst25jAtbTpT06czJW0qU9OmMzZ5HFbziRmCVQbgKBzNADS3uzlU1cr+imZ2HWxkf0ULfv06TRyVzGkzc5g/OeuIZRtDDYA/4KOkswjZvp5dbWsoaltHAO1mclqSyLGPZVSs6LO13xvq4T5xRKPOoPTuE58VXEnQlQTuBHAlgseJlogoSABbXBf2+A5i4zuwx3cQ4+jEFuvCbOmuikKTM4IWyNHsajpsEGo6qunydlHRfuRi9XaLHZE6mSlpU/XXNKakTSM9buCznY/GkBgAIYQZeByYCbiAr0spi0P2XwT8HPACT0kp/3YseUYZAJ/fT0u7h5Z2NyabhUPlTdQ1dVHb3EldUxc1TZ20tB9pwVMS7GQmx5GX4TxiJq8/4KPRU0Odu5x6TzllnXs51LWb0s6iwz59gCRrOtn2MWTbxxBvSR7UBC71cJ84olFnUHoPCL8Z3PGHDYLDlEJXmwO/t5eeuaULbJ1g7cRsd+E3dzBmrOeY0Ucdng5qO2uo66yltqOG2s5a6jprDw8wB8l0ZDE5dQpjk8dRmDiWMUnaa0LKRMwDzPc1VAbgcuBiKeX1QoiTgR9JKS/R99mA3cA8oB1YBVwkpaw6mryBGoC6pk5eeG8P9c1dNLe7ae/0cDRBJhPE2a0kOmwkOmNIcsSw0/869d5SOv1tdPra6PC10ulro93XTKOnBj9H3rAmTGTZCxnvmIWIn4dwzuOThpcHonqvqIf7xBGNOoPS2wgmFaQQCIDXbcPV5sTV6sDdGUtzkxk8DvDGEppIwe5sZ8yC7cdVhj/gp7GrkdrOGiakTGBX/U521e+ktLXkM8feMvP/+OUp9w/otwzVPIDFwDsAUso1Qoi5IfsmA8VSykYAIcRK4FTA8GD3hlYXOw80YDaZsMdYSE20Y7dZsMdYSHDaMRMgzm7FGWsl1m7FHNI6b/U28NLuBwj0MBkWrFjNMSRaU3FYEnFYEnBYEoi3pJBoTTvsz692HaLadcjon6RQKMJM90plOmbAqb8AAoA3FivxeLts5BQcf+vcbDKTFpdGWpw263h6+kymp8/E7XPR5GqisauRJlcTLe5mzi+8cDA/56iE0wAkcuTChD4hhFVK6e1lXytwTIf4saxYH+dxyuyBL7f4HY49EKRQKBTRSjiTwbUACaFl6ZV/b/sSgKYw6qJQKBSKHoTTAKwCLgTQxwBCHWS7gQlCiFQhRAxwGrA6jLooFAqFogcnIgpoBlp81Q3AbCBeSvlkSBSQGS0K6M9hUUShUCgUvRI18wAUCoVCYSxqQRiFQqEYoSgDoFAoFCMUZQAUCoVihDJiF4QRQjiBJUAq2mzkr0kpa4dWq74RQiQBL6DNpYgBviOljJoIKiHEZcCVUsqrh1qXo9FXGpNIRgixAPidlPKModalP+hZAZ4CCgE7cJ+U8vUhVaofCCEswN8AAfiAG6SU+4ZWq+NnJPcAbgY2SilPBf4N/HSI9ekv3wE+kFKeDlwPRE30lBDiMeA3RP59dykQK6VcCNwNPDy06vQPIcQPgL8D0bRo9FeBev05vAD40xDr018uApBSnoIWzfjI0KozMCL9QQwbUspHgWCi7gKgeui0OS5+Dzyhf7YCx7f+3dDyKXDrUCvRD45IYwLMPfbhEcM+4PKhVuI4eQn4Wcj3yEgG1AdSyqXAN/Svo4me+uMIRoQLSAhxE3BXj803SCnXCyE+BKYD5554zY5NH3pno7mCvn3CFeuDY+j9ohDijCFQ6Xg5VhqTiEVK+YoQonCo9TgepJRtAEKIBOBloqcnjpTSK4R4FrgM+OJQ6zMQRoQBkFL+A/jHUfadJYSYBLwFjDuhivXB0fQWQkxHc1t9T0r58QlXrA+Odb2jhGOlMVEYjBBiFPAa8LiUcslQ63M8SCmvE0L8EFgrhJgipWzv86QIYsS6gIQQPxJCfE3/2o42kBPxCCGmoHWbr5ZSvj3U+gxTjpXGRGEgQogs4H/AD6WUTw21Pv1FCPE1IcSP9K8dgJ8oqUNCGRE9gKPwFPCs7q6woKWqiAZ+gzbI95gQAqA5uM6CwjBeA84VQnxKdxoTRXj4MZAC/EwIERwLuEBK2TmEOvWHV4GnhRCfADbg21LKaBqPA1QqCIVCoRixjFgXkEKhUIx0lAFQKBSKEYoyAAqFQjFCUQZAoVAoRijKACgUCsUIRRkAhUKhGKGM5HkACoMRQnwOeFf/eoqU8tPjPD8AuKSUvSYzE0JcDzwNPCGlvGUwukYqQggT8HXgZSll41GOmYU2H2QhWiOuGPijlPLpE6WnYnigegAKIwlN8fy1ox41cEqAZcCWMMiOFN4EnkRLjfwZhBA5wHLgfLSKfx0wFXhKCDEsjaIifCgDoDAEIUQcWiZKt77pKiFEjJFlSCk/lFJeKqX8q5FyI4ypfey/AkhGy5szV0p5DhCcCX5zOBVTDD+UC0hhFBehJVB7A8hBS6H8BbQp88DhxT/uR+sdJKDl2PmelHJlqCAhxKVoaa8zgBeBW6SUnt5cQEKIr6OlE8jV5X1XSvmJEOJa4FngeSnltfqxAigCNkkp5wghEoE/oGVzDACvA9+SUjbqWTUPoPU4NgHfAjzAd4F6tLz1OWh5mb4hpXQfSx993zPAdWiZI28DgnmGbpdSbhRCfISWWhigUghxg5TymR7X2aO/nyGEWKS72f4HnEq38UUIkQs8BpyHlqPmI+AOKWWZvr8QbZ2Ds9DqgffQFhc6qO8/CKQDD6Jldn1LSnmNEGIR8CjaYjllwENSyr+giEpUD0BhFNfo76/RXelf2+OY+4HvoSXOWgvMA/6rV1ZBYoBn0CoXK3Aj8JXeChRCfB5tVaYk4BNgMvCuEGIC8Apakr+LQ3oiX9LfX9Df/4FWIZegGYavoaUkDuU84E5gP5Cln/MKUKrvv05/9aVPKE8BTrSU0wvQjBDASiCYA+cdXa+evAk0AVOAVUKIcuCvgEdKuU7Xw4T2P3xR13M/2iI3y/T9qXpZl+tl7EczgiuEEGkhZTmAHwLbgPV64rZ3gVnACjQ31eMhSRUVUYYyAIpBo1co56O1NN+g2wBcGKxQhBAOtIq0BZgupTwL+BHwATAqRJwJuFxfIeo+fduCoxT9Q/19oZTyc2gVWixwq56W9zW0yvhz+nFf0nX8lxBiHFoF+SkwQ0p5MvAv4CwhxIyQMmxoLfX5aAuu2IEHpZRnog3EAkzrS58eev9PL+90/fscACnlT4EafdsNUsoPe/5gKWU52rXeqG/KRXP9rBFCfF/fdrau70q0az1H/207hRDpuj55wBIp5Uwp5Uy05VHz0XomQUzA/0kpT9MXULodiNe3nYO2jkYHn137QRElKAOgMIIr0VruK6WUdVJKCexCqzy/rB8zXj9mq5SyHkBK+YCU8jIp5doQWQEguMbBAf3dcZRyg/5yqUcQBSOQgit4PRfUT1/zYRracppVaC1ogEWAXz//Kz3OByiTUhZLKQNAhb7tE/09+D04YNuXPkGCq43tRTNIvQ74Hg0p5Vop5Vy0Hsb3gQ36rvt0Yxz8bcullH79nKullNdKKevQDBrA8yFig59P5khC15sO/r4n9N/XgPbfzNDde4ooQ40BKIwg6P45Xa8YQvka2rrFJv374XtOd814g5WUjltKGcyrHnw30TvBSmdZj+2H9PcPgHK08Yn9+rag+yd4bgmwucf59SGfQ1P8BvXs6PG9v/oE6Qj57EVLR94vhBB3A6cBP5FSbgaKhBCPADuBScAEer/WsSHpinvqTcg5Pf+/lpDPwd+3kiOvEWg9HQ+KqEL1ABSDQghRgLaGrh+tEgp9ASwQQkxEC1n0ADN1XzJoa8G2CiFuHGDxu/T3+6SUl6INHO8FlgLohmUJWr75u9Aq3tf0c3br743AF/Xz16KFVYb2SAzTJ4Rj5WAPVs5HezYz0BZP/57u6wctKihF/1wRosfZQoigcXlPCFEvhJhKt8H7aojc4OfQFj8cuchJUO4y/ffdjGZgX5NSth7jNykiFGUAFIPlarTW47tSymmhL7TIE4Cv6T75x9FcBlv1tZh/hFbhfTDAsoODp8uFEB+gDZB+jyNdRkE3UBKwNLgGrZRyN1rky0xgtxBiHdog9S0cuR6w0fr0RZ3+/roQ4upe9j+KNgh8NZqr6V1gD9oA9UtSylLgfbRKfj6wXQixAc1Il6JV4n8FaoFrhBBbhBBbdXklaP/R0XgScAEPCCE+RpuPcRvdYyCKKEMZAMVgCbp/elvO7wn9/at6a/WHaKGJZrSB3fXA56WUPV0k/UJfP/b/gGq0Cq4CuF5K+VbIMTvonjj2Qg8RX0Hzfaej+c3/B5w70NWo+qNPP/gNUAkIejEcegW/GK0n4wTOQHPTPABcrx8TQJsbsBQoAMbqx18spQzoYyAL0QbrC/X9rwKnHm32sS53P1rvYz3a/2cGHkIz5IooRK0IplAoFCMU1QNQKBSKEYoyAAqFQjFCUQZAoVAoRijKACgUCsUIRRkAhUKhGKEoA6BQKBQjFGUAFAqFYoSiDIBCoVCMUP4f1qQURZWr7AIAAAAASUVORK5CYII=",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"sns.histplot(df_mindset['achievement_score'][df_categ['intervention'] == 0], kde=True, stat=\"density\", linewidth=0)\n",
"sns.histplot(df_mindset['achievement_score'][df_categ['intervention'] == 1], kde=True, stat=\"density\", linewidth=0,color='g')\n",
"plt.title(\"Achievement Scores by Treatment\", fontsize=20, fontweight = 'bold')\n",
"plt.xlabel('Achievement Score', fontsize=14, fontweight = 'bold', labelpad=5)\n",
"plt.legend(['Untreated', 'Treated']);\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Method Implementations\n",
"\n",
"The following code block implements the naive ATE, the standard IPW, and finally the AIPW methods as python functions. Note that propensity score, or the exposure model, is constructed as a *Logistic Regression problem*, and the outcome model is generated as a *Linear Regression problem*.\n",
"\n",
"We do this to allow us to readily run many iterations of each method. We will use a bootstrap subsample method, where we will sample 1% of the original data (~100 data points), 100 times. This will allow us to generate a distribution of ATEs with an empirical standard deviation. Thus we can report our results comparing each of the three methods using various exposure and outcome models with 95% confidence intervals as well."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"#### Define Estimation methods ####\n",
"\n",
"### Linear Regression T on Y\n",
"def naive_ATE(df, T, Y):\n",
" return df[Y][df_categ[T] == 1].mean() - df_categ[Y][df_categ[T] == 0].mean()\n",
"\n",
"### IPW\n",
"def IPW(df, X, T, Y,true_ps = True):\n",
"\n",
" if true_ps:\n",
" p_scores = LogisticRegression(C=1e6, max_iter=1000).fit(df[X], df[T]).predict_proba(df[X])[:, 1]\n",
" else:\n",
" p_scores = np.random.uniform(0.1, 0.9, df.shape[0])\n",
"\n",
" df_ps = df.assign(propensity_score=p_scores)\n",
"\n",
" weight = ((df_ps[\"intervention\"]-df_ps[\"propensity_score\"]) / (df_ps[\"propensity_score\"]*(1-df_ps[\"propensity_score\"])))\n",
"\n",
" weight_t = 1/df_ps.query(\"intervention==1\")[\"propensity_score\"]\n",
" weight_nt = 1/(1-df_ps.query(\"intervention==0\")[\"propensity_score\"])\n",
"\n",
" y1 = sum(df_ps.query(\"intervention==1\")[\"achievement_score\"]*weight) / len(df_ps)\n",
" y0 = sum(df_ps.query(\"intervention==0\")[\"achievement_score\"]*weight_nt) / len(df_ps)\n",
"\n",
" return np.mean(weight * df_ps[\"achievement_score\"]), p_scores, df_ps\n",
"\n",
"### AIPW\n",
"def AIPW(df, X, T, Y,true_ps = True,true_mus = True):\n",
" if true_ps:\n",
" p_scores = LogisticRegression(C=1e6,max_iter=500).fit(df[X], df[T]).predict_proba(df[X])[:, 1]\n",
" else:\n",
" p_scores = np.random.uniform(0.1, 0.9, df.shape[0])\n",
"\n",
" if true_mus:\n",
" mu0 = LinearRegression().fit(df.query(f\"{T}==0\")[X], df.query(f\"{T}==0\")[Y]).predict(df[X])\n",
" mu1 = LinearRegression().fit(df.query(f\"{T}==1\")[X], df.query(f\"{T}==1\")[Y]).predict(df[X])\n",
" else:\n",
" mu0 = np.random.uniform(0, 1, df.shape[0])\n",
" mu1 = np.random.uniform(0, 1, df.shape[0])\n",
"\n",
" return (\n",
" np.mean(df[T]*(df[Y] - mu1)/p_scores + mu1) -\n",
" np.mean((1-df[T])*(df[Y] - mu0)/(1-p_scores) + mu0)\n",
" ), p_scores, mu0, mu1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Experiments and Results\n",
"\n",
"The following block shows our bootstrap sampling method results displayed (100 iterations for ~100 samples of 1% of dataset). In this initial experiment, we correctly specify both the exposure and outcome models. The results are displayed in the plot and table below."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Hide cell\n",
"bootstrap_sample = 100\n",
"AIPW_ates = []\n",
"IPW_ates = []\n",
"naives_ates = []\n",
"for iSample in range(bootstrap_sample):\n",
" df_bootstrap = df_categ.sample(frac=1)\n",
" ate, ps, mu0, mu1 = AIPW(df_bootstrap, X, T, Y)\n",
" AIPW_ates.append(ate)\n",
" ate, ps,_ = IPW(df_bootstrap, X, T, Y,true_ps=False)\n",
" IPW_ates.append(ate)\n",
" naives_ates.append(naive_ATE(df_bootstrap, T, Y))\n",
"\n",
"sns.histplot(AIPW_ates, kde=True, stat=\"density\", linewidth=0,color='g')\n",
"sns.histplot(IPW_ates, kde=True, stat=\"density\", linewidth=0,color='b')\n",
"sns.histplot(naives_ates, kde=True, stat=\"density\", linewidth=0,color='r')\n",
"# plt.vlines(np.percentile(AIPW_ates, 2.5), 0, 20, linestyles=\"dotted\")\n",
"# plt.vlines(np.percentile(AIPW_ates, 97.5), 0, 20, linestyles=\"dotted\", label=\"95% CI\")\n",
"plt.title(\"ATE Bootstrap Distribution\", fontsize=20, fontweight = 'bold')\n",
"plt.xlabel('ATE', fontsize=14, fontweight = 'bold', labelpad=5)\n",
"plt.legend(['AIPW', 'IPW', 'Naive'])\n",
"\n",
"Results = {\"AIPW\":{\"Mean ATE\":np.mean(AIPW_ates), \"Std Dev\":np.std(AIPW_ates), \"[.025\":np.percentile(AIPW_ates, 2.5), \".975]\":np.percentile(AIPW_ates, 97.5)},\n",
"\"IPW\":{\"Mean ATE\":np.mean(IPW_ates), \"Std Dev\":np.std(IPW_ates), \"[.025\":np.percentile(IPW_ates, 2.5), \".975]\":np.percentile(IPW_ates, 97.5)},\n",
"\"Naive\":{\"Mean ATE\":np.mean(naives_ates), \"Std Dev\":np.std(naives_ates), \"[.025\":np.percentile(naives_ates, 2.5), \".975]\":np.percentile(naives_ates, 97.5)}}\n",
"\n",
"df_results = pd.DataFrame(Results)\n",
"df_results.T"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From the results it is clear both IPW and AIPW account for a positive bias we hypothesized. They estimate the ATE at ~$0.39$, up from the naive ATE estimate of ~$0.47$. We also note the IPW and AIPW methods agree with very close estimates and with very similar 95% confidence intervals. This is unsurprising considering the exposure model is correctly specified using logistic regression for both methods.\n",
"\n",
"Now that we have propensity scores, we can also perform a quick positivity check visualizing the distribution of our propensity scores to ensure we meet the positivity/overlap assumption which the plot below demonstrates."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAE1CAYAAAAcUKCZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABZt0lEQVR4nO2dd3gVVfrHP7fkpgcChAChF490EEVBBcQCYkHsvbvqquv6s6y6umJX7KtrXXsDsYB1FUFEepXOofcWCDX9lt8fZy65CTfkBnIzN8n7eZ48mTtzZuY77Z0z7znnfR2BQABBEAShduG0W4AgCIJQ9YhxFwRBqIWIcRcEQaiFiHEXBEGohYhxFwRBqIWIcRcEQaiFuO0WIFQ/SqnhwCNlZgeAYmAfoIGXtNZfVqOmYJ/c+VrrHmWWdQc6a60/C5k3Eehv/UzXWu+O9j6rEqVUCnA7MAw4CkgCNgMTMed+QZnyHwDXWD97aq3/jIauilBKtQbWWD/Haq3Ps0OHUDFScxeCOAAP0BDoC4xWSt1lpyClVKpS6jVgDtC7tuzTenEsBZ629lEfc+5bA9cCc5RSf4nGvoW6g9TchVeBCYALSAMuA063lj2plPpMa72tGnQMs/7vCZnXC7itnPIPAY2s6dxq2ucRo5RqBvxKifb/AaOt6fOAczDP5etKqUVa66nR0iLUbsS4C3O11mOCP5RSHwLTgeOAROA04NNoiwjVEGH5ydW9zyriWUoM+/Na63tDlr2nlHoJ+DvmZXs/cG71yhNqC2LchVJorf1Kqd8wxh2gRehypVQX4B5gIJAJbAfGAU9prVeWKdsDU8M+EWPQCjH+/I+BV7TWgZCypfzfYdoF7lRK3Qlcp7X+oKzPHWgPzLJ+j9JaX1pGy3vAddbPHlrr+ZXZJ3AHcAymbaK51npzyLYvAr6wfv5La/04YVBKpQOXWD9zrHNTlkeBzsBkYHy47QBOpdTdwF+BLIyL5xmt9agy+/MAd2N89W0wXygTgEe11kvD6OsP/AM4AUgA1gNjgBFa65xytATXfTzkeGYDA7TWh/tFJVQB4nMXwtEmZPqAy0IpdT7GgF6DMfoeoDnG+M1TSp0aUrYbMAW4AGiCqUgkYwzkS8DzVSlYaz0bWGj9PEcplRyiJR443/o5X2s9/zB28ZH134E5plDOD5k+VANsHyDOmv5Na11YtoDWerfW+gyt9WNa6ynlbOffmPPXFogHegAjlVKnBQsopdzA98BTgMJcqwzMy2WmUqpP6AaVUjdiDP+ZmJdlorXeP4AJSqm08g5KKXUJJYZ9NXCWGHb7EeMuAKCUciqlMpRSd1DaeE20lrcAPsDU6ALA68DVwLtWuRRMI2zQ5fBXTA8QP/AExqjcAWyylt9pbbM8RgIPh/z+HuMjn3CIdd63/icBQ0Pmnw3Us6Y/OMx9fgZ4rfkHzo/14hhi/ZyutV51iO23C5lefYhyFXEC8BymfeT3kPm3hkz/jZK2k2+BqzC1+G2Ya/WhUsoJB67tfzD2oAB4HLgS+MVavzvwf+GEKKWOpeS87wTO1FpvP4JjE6oIccsI7yul3i9vWcjn+y1AqjX9kNb6KWv6Y6VUDnAvpsb3F0xtMcVang98p7WeCaCUmozp+rcY2FKeKK31MqtskFUR+Mg/wfi04zCGL1iLvtz67+UQNeuK9qmU+gnT4HmyUirTamg+DdMQDRW3TYTWfgsqKHsontda329pGo9xjQEcHVLmJuv/CmCY1tpvlV8NfAN0wLi1fsMYfo9V/n6t9StW2S+BaZgXUbhr1RQYi6nl5wPnaK2XH8FxCVWIGHchHAWYmtz9IfNOCZl+q0z5tzDGHaAfxriPBK7AuGJmKKU2A39gvgS+DfVZVxVa62yl1A+YXieDlFINAB9wllXkxyOsVX6EMe5Oax9vUeKS8QKjwq92gLyQ6eRyS1XMr8EJ65jzMF8riXCgD33Q0HcAfEqpcNs5AWPce4XM+1/ItgsxbrTyCO0qWoDx/QsxgrhlhFcxroehwGDMg95Ia32P1tobUq6h9b9Qa72zzDY2hUw3ANBaf49xxSyz5jezfr8BrFNKvRF0C1Qxwa+QOIz75HyMXxoO7ZKJhO+AXdb0hUopFyW9WcZprbMrWH9dyHTb8goppS4OcW+Fo2zjZrH1P3g+61egI0gz63+9kHk7Ily3LOnAg4e5rhAFpOYuzI3A3QHm0/8oIF4p1bCMgc8KmT5g4LTWXwBfKKU6YdwXJwJnYIzPLZgG0NePSP3B/IjxK2diatdBw7cT+OFINqy1LlRKjcJoH4B5KQaN8CcRbGIKpg3CCfRXSsWXbVRVSnXAfAH4lFJjtdZlG2+xtnEo9oVMLwT+VU654MsmtJ9/Q8y5Cuo5CsjWWu8iPK9hGop7AXcopV7VWm+oQJ9QDUjNXYiUaSHTN5VZFjqa8nellEMpdY9S6kOl1BfAUq31v7XWlwCDQsr2rWCfoUbMEYlI62sjaGhPxbxUAD7TWhdFsImK9hnsNeMGXrGmczG+54q0bcN0LQRT0x0eulwp5QBesH66OMxatNZ6DyUhArKAiVrrMdZL3IU5J01Dtj83ZPVg43CwK+VEIEcpNSfMrmYCd1JSY0/ANMYKMYDU3IVI+S+mt0sCZuRqFjAD42MPGvts4B2tdUAp1Q/jnwZwK6W+wrhKrg3Z5hoOTaiP+hSl1NWYF8Ws8laweB/TMySeyrtkDrlPrfU0pdQKjC876NYYU4muf/dgXjr1gPuVUp2BrzANmlcBJ1vl9gJPRrjNcLwPPIZxk01QSv0H45d/FPNiKcK4mcC8DP9laXhWKdUY0+B9OeYlACF+/hC2WA21v1hjI04BrlJKvVg2No5Q/UjNXYgIqxfETRij4MQEvfqYEsO+F7gw5PP9Tkp88cMwBuR9SgYerQZermC3yyjpVdIV+JCS7n2H0rqYkgFNAIu01nPLK38Y+/yozO+IR/Bqrddg2jaC7qtzMC+etykx7LnAJVrr9ZFuNwzPU/K11RPzcv43xrAD3K213mhpWofpOhnAGPgHMNcrWItfQMU18ges/07gmSPQLVQRYtyFiNFaf4Lxrb4PbMAY+k3Ae5hRn5NCyq4BjsX0x14C7KdkhOqLQO8wDbNl97cfM2BqEcbgbuAQ3SfLENq984MI14l0n1+ETGdjRuhGjNZ6OmaA0OMYw5mLOTcrML2Uumit/1f+FiLaRz7mC+FhTC08H9MQOxE4W2v9Wpnyb2FeYr9gfPB5mOv2BNDfOi+H2t8MSlxTZyqlBh6JfuHIcQQCgYpLCYJwAMtV86H18z9a69vt1CMI4RCfuyBEgDWKsyfGBx3qCy9vAJgg2IoYd0GIjBYc3CPmS611uF4kgmA74nMXhMhYi8mUVGhNv4CJrSMIMYn43AVBEGohUnMXBEGohYhxFwRBqIVIg6oNhMn4A2YASTEmLogGXtJaf1nN0moEZc7fsDJheT2YNHWvV9Q3+zD268DEOb8aE+M8HXO9lmDyoL4RYYiDGoFS6lrC9wYqBHYDf2K6gn5XtkDZLFdHoKENMEhr/ebhbqOuIjX32MGBGR3YEBNzZbRS6i57JdUslFKnYwYfPUsVV1ysCJajMKNTT8NkNXJjDPyJmNG2E0MzQNVi4jGB2QYB3yqlXqzqHSil4pRSj2JenIOrevt1Aam528+rmEw/Lkwyh8soGe7+pFLqMyvglFDCSEytEUx8myD/xMR8iQaXAxdZ08sx4QI2YlIS3olJJdgHuI+Dv8pqA6Mw592FScTSG7gRUyG5Syk1X2v9YUj5Ydb/PRweWZQfzVKIADHu9jO3jFvhQ2A6JkF1IqaWGHHskrqA1noZJXHiq4tzQ6bP1lqvCP5QSk3ChPMFE2a4Nhr3ZWVCQ3+olPoeE2IZ4FGl1Cdaax9AhGGkhSgixj3G0Fr7rQh7x1mzWgAopQZgsuaAqTH1w9Qk84G/aK2/ssqdCNyFce00xMR++R54Wmt9IEZKme1dh4k78i+gM6Y/9zvAc8GHNWS9gZgacm+MK2kRxu/6cZlya4FWwHhMko6nMYYvBVPbvt+KRxIsH2fpvhxT+/ZgYsj/AQy3DHqw7HDK+NxDfLxBdiml1mEiWX5rzXtba31zGZ2LgU6Y2DeZWus8whPqbjlbKfWy1joAoLWeqpS6DvM8HbS+Uqo/JtH0CZiomusxoX9HaK1zypRthYkceRYm6uQeTJ7UZ8oGPws55ncx7ot/Ylwmr2ut77PKXITJf9od06YzB5Om70eOEK31T0qpKRi3VCvMPfdHGW2lfO5KqQsw16Q7Jm3jHsxX2AtBTWF8/UOt7T2qtR5emXulLiM+99ikTch0uM/a4ZhGvURMSNc5AEqp24FJmAxETTE3fRvMwzRPKdWtnP1djhl92QtjfNpijHGpxjTLgI0DBmKMdDJwPPCRUur5cradjqnV3oTxUydiEl1MUEplhpR7D+Mr744JTeumJHvTJKVUuZmLKuAnSiIwnmdlTwoez1EYww7wzSEMO0BobtUXgZVKqReVUmcppVK01h9orf+rtS6Vo1UpdSPG7XYm5lwkYoKG/QNzDtJCyp4IzMNE3GyDMdSNMS/xGUqpK8rRNggzqKoB5postrb3CCbI2QnWftMwYXl/sBKhVwWTQqaPP1RBpdQNwJeYyKD1MS6eBpj76Tul1HkR7jNa90qtQox7jKCUciqlMqyHLjT7zsQwxZtjHtqrMDXytVYW+pcw17QQE3b1WuBra51M4GurN0lZTsckbLgRU3vPt+ZfpZQaZOlrhsma5MRk8Pkb5gUz3ip7t1LqZA7mGMwXxD3ADZgsSWAeyqusbTfB9EIBmG2Vu5SSpBsZlM7nGo5hWEbN4krMF40X+Nya15iSsLpQ+jxX5Pp6DdOLKUhbTO3xe2CHUmq0UqpL6ApWPJr/YM5ZASYK5JWYyItgjNP/WWVTMT7tYEjezzDRKV/E5Gd1A+8ppTqG0dYcU/u9CpNAZIx1PwS/bhZhEqpcj1URAJ6vIiMYmgu3WbmlDME8u9swPZouxsScD4aRftbqkTQBCP3Cmom5viOr6F6pE4hbxn7eV0qVF3zqfa11uKTD24AryuQ4vYuS63l9SA0ymA3pIqAdcCHGcISyAxigtd4HoJTaiKkdgUly/TPGcCRY867UWk+2yo7EuH4yMC+HP8LovSoYwlYp5aMkBG87639KSNkFmJgte4FRSqnxmJfJojDbPYDlmvl7yKwftNa7g+cA8zICc/wTrelgcuvthE9GEbr9fVYCkpcwNURXyOJ4a7tnK6WGhYTrvQrz9QTGDfUKgFLqS0ys9dWUhBO+FGOkwSQ8OZDdSim1CvOS8GCuc2jmqyBXaa0XYRk5q5YczCR1ZjB2u1JqLOZLxoN5eRxp+0B+yHRqBWWD13kT5ktpPaZX2AxL6xLL1bVeKfVLyHpbgj58pVT7kPmHda/UFaTmHpsUYD6zwz3EANPKGHYwn9tg+l1/XmbZWyHT/cJs76egYbcYHTLd1fp/bMi8P5RSAcsPWoQx7GA+/8sSoKR2D6UbQhMAtNYrMbUwMLXLHUqpqUqpZzH+6UkRJJ8uF8tXHXzgh1lpAFuEHNPIsm0L5Wxnu9b6Ckw7yM2YDEqhuUUTgA9Cvo56hSw7EJ9da12otT5Ga32hFUcdSq4flL5eYBJtBPWFu367LcMeSuj12hByvXZS8tyHu16VJTFk+lBuLTBfJmC+5tYppRYrpd7CGP3frRwAhyTa90ptQoy7/byK+eQciunP2wtopLW+J4wBDxIut2ZD6//mYENfCJtCphuEWbdUo541+CeYjShY20qnYsJ9ludqrYtDfocmhA7NUToE04e8AJOOL9itcDyglVJHaoiCDb7NMA1/54csi6g3klIqXSnVWGu9RWv9ttb6QsyL7TxKkkpnYhoYwaTSC1JRPtSGIdOh1wtrYFRw/XDXL9y2D/d6VZYmIdNbKyh7P6a9KGh8O2EqMKOAjVb7RCRE+16pFYhbxn5KdYWMkPww87ZjPuubKaUcZQx8Vsh0uFpNqYdcKZVCiQtmt/U/tGZ/E+ENSrgodBXWiAGs2tY1SqnbMLXY/tbfsRj3zTdKqeaR1LDL4RPgKYw75QJKatUrtNYzD7Wi1YNFY9wvszA9hYK6fcBYy2Xzf9bsxtb/0MbwhpS8AIKNudkhaQm3h5TNIsRQWl8Cwa+jcNcv3P0QvF4+jEsu3LWJNO/roQhtRJ1+qIJWZeVRpdSTmK+GAZg2kFMwL8K3lFIztNYLK9hOtO+VWoHU3Gsm/jDzgvkyUzE+4VBC3Tu/h1n3TKVUaM3xvJDp+db/eSHzCrXWY6yX0u8YY9mRkhdBpVBKHauUel4p9SMmpdt31pfLcZjeLmBqiG3K3wpQ+ryEfhWgtd5MiXvockpq15HU2tdTYpiPVWVSyFk9cPqGzAq6F0K7Lg4JKe/B+P1zlFLBBs5pIWVvojQ3UPKshrt+4e6H4PVyATtCrtdSzFdiO8p8sVUWq/dV0J20jtI9Z8qWbaKUesxqb3hAaz1Za/2E1noQ8JBVzEnJyyLstazCe6XWIzX32sPrmEY9B8bv2xVT2zyXkl4hyyjpPRNKCqYL2cuYz/5/hiwLJoP+BPMQxgGvKxPzYyXG9xz0A99HeONTET7gbmv6WKXUcxhj0ZGShNoFVPzZH+rzvUMptb1MTJKPgDMwrpMgZRuXD0JrHVBKvYLpfufADLl/A2NA62EanYOugNDk3J9geh95MD1BGlvLL8d0VYWShtyRwKNAI+BmZcIY/AJ0w/QsAXMOIh3q/z5m/ALAl0qpEZha/93WNsH0Vok02cjRIV0VE4EuwK2UNCw/WkFNeRdwC5Yby+oCOhvzlXNLSLngizH0WvZSSl2F6ZmTQ9XcK7UeqbnXErTWEzEZ6P0Y98GDmF4iQcO+GTi/HD/+Ykzf67cxXSiDvR7e1FbSa631akq6sqVgDNGnlBj2PzDdBQ9H+zxKUtdlACMwftjhmC6TYGp7FQUCC/26eBR42epaF+QbSruXZoaONK2AF631wfQlvwdz/K9T8hWwC7gsZHBTsMtoAGPgH8AY/GAtfgGmeyTWYKbLgL3WsisxL6N7MJWwIuBqq0GxQrTWf2Di3YAxoM9j7oegYR+NaRCOlEswx/8N5oX4ICV+/de01odMN6i1LsR0zS3GvBDuwbzQ/k1JLXus1nq8VX4HJjk5GDfVR8AlVXiv1HrEuNcitNbPYoztF5gudkWYmtArQM9yulWCqT0OwfSVLgJWYR6+28ps/xWr3K8YQ1aA+Rp4BBiitQ7n+41U+0MYd8E4TIOiF+MK+QU4R2v9cgSbeQFz7DswI05nUdJ2gDVI6fuQ8hGHdbBeihdQ0k8929K4B2OkRwBHl/UXW71hTrfW2YOpkS4BnsC4FfaHlP0VY3xfxXSTLMT44r8AemutQ3sxRaL5LkvvFMxLLRfjZrsT05U2nDsnEooxNePvgLO01hENiLJGoB6HOe/B4wtep79jvjxDuQZTu8+39rfW2k5V3Cu1HsnEVEdRpcMPvKK1/rt9aqoHpZQb85VyFMYVlKUlKJtQSxGfu1DrsUbZxmNqsUdZs38Swy7UZsS4C3WBuykJowzGHTDcHimCUD2Iz12oC8zD+JyDERYHa60j7SUiCDUS8bkLgiDUQmLGLZOdvc/Wt0x6ehK7dlUUGqP6EV2VQ3RFTixqAtFVWTIyUh3h5otbxsLtdlVcyAZEV+UQXZETi5pAdFUVYtwFQRBqIWLcBUEQaiFi3AVBEGohUWtQtZLYfgi0xowGvEkS1wqCIFQP0ay5DwHcWuu+mDyJT1ZQXhAEQagiomnclwNupZQTk3W9uILygiAIQhURtUFMVo7KsZjwsI2As7XWU8sr7/X6AjWtq5EgCEIMELafezSN+4uYjD0PWIZ+AtBVa10Qrrzdg5gyMlLJzt5XccFqRnRVDtEVObGoCURXZbFjENMuSnJI5mAy+EjVXBAE25g7dzaPPPJAqXlvvPEqP/74XdjyW7duZfLkcrMHRsxXX42KuOyh9FSGaIYfeAl4Tyn1ByYLzYNa66pIyCsIh2TECE/Y+ffdV1TNSoSazty5s1i3bi0nndSv4sKH4MMP3+OCC8qmNo4uUTPuVoaZi6O1fUEQajbDh8fz3XcHmyCnE/z+5MPa5jnneBk+vPCw1r377r8RF+dmy5bNDBx4OldeeS2ffPIBBQUFdO3aja+/HkVSUir79u3juede5oUXnmHjxg34/X5uuulWjjnmWH777Ve+/no0QXf3E0+MYOzYr9i7dw/PP/8Mf//7PTz33FMHrTdx4ng+/PBd6tdPp7i4mFatWh/WMYQig5gEQajzOBwOtm3bwhNPjODNN9/ns88+wuVyceWV13L66YM56SSTe/v00wfzyiuv88MP31KvXn3+8593eOaZF3jxxREAbNiwnueee4XXXnubli1bMXPmNK655gbS0upxzz338913Y8Ku9/rr/+bll1/nxRdfIyEhoVydlSFmokIKglC3GD68MGwt2zRcRseDGx+fQFFR6V7Z+fl5eDzxtG3bHrfbjdvtJj4+vIFt2bIVAKtWrWTBgnksWbIIAJ/Py549u0lPb8ATTzxCUlIS69atpUuXbqXWD7deTs5OkpOTqVevPsBB6xwuYtwFQagztG7dmhUrNDt27KBRo0YUFhYyf/48OnRQOML0OXE4HAQCJXnEnU7j7GjVqjWNGzfm6quvp7CwgA8/fA+Xy827777FV1+ZHOx33XXbAfdM8H+49VJT09i/P5ddu3aRnp7OsmVLaNw484iPVYy7IAh1huTkFO644y7uu+9O4uMT8HqLueCCS8jKas7s2TMOKt+uXXs++ug9jjrq6FLzhw49n2effYLbb/8Lubn7GTbsIpKTk+natTvXX38liYmJpKamsmNHNgCtW7fhscce5v77Hz5ovbi4OB588F/cffftpKbWw+2uGrMcM5mYpJ97eERX5cjISOXee8M3qNnZWyYWz1csagLRVVkkWYcgCEIdQoy7IAhCLUSMuyAIQi1EjLsgCEItRIy7IAhCLUSMuyAIQi1E+rkLglBnePXVl9B6KTk5OykoKKBZsyzq10/niSeePaztjR37NWeddW5EfdPHjPmSnTt3csMNNx/WviqLGHdBEOoMd9xxFwA//vgd69at5dZb7zii7X388fsMHnxWlQ08qkpiT5EgRAkJBRxbDJ/6EN+tGnPQfKfTgd9/eGMaz2l3HsP7PlGpdZ58cjh79uxh7949jBjxMp999hHz58/F7w9wySVXMHDgacybN4dPPnmP4mIfBQUFPPTQoyxYMI+cnJ0MH/4gTz/9Am+++dpB682f/yevvPI8aWlpOJ0uOnfucljHdTiIcRdqFcMnDmeKx1vu8hOLHq5GNUJNoVevY7nkkiuYNm0KW7Zs4o033qOwsJCbb76O4447njVrVvPcc8/hdCbx0Ufv8dtvv3LNNTfwwQfvMnz4U+Wu9+qrLzJ8+JO0bNmK559/ulqPSYy7IAi2MLzvE2Fr2XYM8w9Ge1y9eiVaL+P22/8CgNfrZevWLWRkZPDkk0/idMaRnb2drl27l1q/vPWys7cf2HbXrt3ZuHFDtR2TGHdBEOo8DkdJtMeePY/lH//4J36/nw8++C9ZWVncdddfGT9+PPn5AZ544pFS6wUCgXLXa9iwIWvXrqF16zYsXbqE1NTUajsmMe6CIAgWJ57Yj3nz5vDXv95Ifn4e/fqdQlJSMoMGDeHiiy8mKSmZ9PSGB6I9du/eg3vu+RuvvvpW2PUefvhxnnzyEZKSkklKSqpW4x61qJBKqWuBa62fCUAPoInWene48hIVMjyiq3L8Z/EL/DKucj736mhQjcXzFYuaQHRVlvKiQkYzh+oHwAcASqn/AO+VZ9gFIVp4C+PZNud4ti/oicMRYGdDN717+2jSJDZCXQtCtIi6W0YpdSzQWWt9W7T3JQihFOclsviTG8nfkYkzrhCn28fynS6WL3fRo4ePU07x4nLZrVIQokN1+NwfBB6tqFB6ehJut71PWkZG9fnDKoPoqhyeODe+YjeLRl9D/o5MmvaaTdvTxuGKL6T12uH89BP8+aeLnTtdXH45ZGTEV4uuWDxfsagJRFdVEFXjrpSqDxyttf6torK7duVFU0qFxLA/TXRVkqJiLxsnn8S+Tc1p1Hk+rc4Yi98RwO+Fxo0Luewy+PlnN1q7eP99P1dfnUdGRnTdNLF4vmJRE4iuylLeCyfagcP6Ab9GeR+CUIqifalsmtaPuKT9tBn8LQ5HacPt8cDZZ3vp3t3H9u1Ohg5NZPPmsG1SglBjibZbRgGro7wPQSjFhj8G4i/20Pr0H3DHh8+n6nDAaad58XgCzJrl5pxzkvjyyzzatCl5EUi4AqEmE9Wau9b6Oa31y9HchyCEkrcnkexFPUhI30HjbnMPWdbhgH79fDzwQCEbNjg599wkli2TKNhC7UDuZKFWMe+nngS8cWT2monDWbEf3eGAu+4q4sknC9i2zcl55yXy55/yWAg1H7mLhVqDzwezxx6HM66Ixl0PXWsvy003FfPKK/ns3u1g6NAkRo6UwdtCzUbuYKHGUtYnvm6dg91bU2ncYzbuxIJKb++yy7w0aJDPbbcl8re/JdKtm4+BA73EYKhuQagQuW2FWsOyZWacRKPOfx72NgYN8jFuXC433JDIggUuNm+BM87aQ8NGPjykMtXzBCNm+sKue1/vBw97v4JQ1YhxF2oFPh+sWOHEk7KPtBbrjmhbDZrt4pqXP+ax0T+yI3Umn6VsByDB35AkGtFq37E0T21RFbIFIWqIcRdqBevWOSkocJDVe0lEDall8fl9/L7xN0Yt+5Qf13xPoa8QmkJiUUsKVw/G73VR3HQxOamakVrTvn4HBrUeQqI7MQpHIwhHjhh3oVawfLnpG5DReVGl1stxLOfJ6R/yhf6cLbmbAWiQ0IDeTU5g/5rOJJBOUf0UVv0wjN2rFM52E6h/xV9ZuVuzbcn7XNDhYholNqry4xGEI0WMu1DjCQRg1SonSUkB0lpspLj8iL8AFLIH7f6SRXEfs9k1HeZCmqceV3e6Hq+/mKbJzXA4HExZY3z4npT9HH3xx2ybdxzrfj2TnMcX0/jav7K91duM0p9x8VGXkpHUuBqOVBAiR4y7UOPZutVBfr6DLl18B4UaCKWAXfwSfxtL3J/hdeRDwEEr76k8cOblnNnmbBLdiYyY+VTYdR0OaHLMLOq1WsOW8Rey7f23SD69Pbkn/oMvV4zi8qOvitbhCcJhIcZdqPGsWmVcMu3a+dkdZrkfL5tcU9jknErA4aeevw1di66hk/cK0gItOL9D5OEEEhvu4LLn3uT3d4cw7/t7cea5yT39//hqxWj+0fufpMXXq5qDEoQjRIy7UONZvdqJyxWgVauDjXsxuSxzj2a/cxOeQCqnF7yK8l6Ek8MPL+2O83HqLd/RvMsafn7lVopS1pPT52VuG3crH531KQ5H6SBkI0Z4SE6G3NzS/fIlRo0QTWSEqlCj2bcPtm930rx5AE+ZOF+F7GVh3Afsd26ika8z3YtvpqP30iMy7KGokxZx1SuvkbHyLlgzgJ/Xf8/j41+tkm0LwpEiNXehRrNunamftGnjLzW/mDyWxn1OoWM3Wb4TaeHrj4OqCes7ZWroy2E3bS94j8Cke9nRSPPaskeYM7+IE08uqZVP8bjw4OY4HqiS/QtCJEjNXajRrF1rbuHWrUuMewA/y93fkO/YQVNf7yo17OFwur10GDiVE9w3ATDN+w5jXj6NvD1JUdunIFSEGHehxhIImJp7SkqAhg1LeslsdE5mr3Mt6f4OtPKdFlXDHspJJxdxXL3BkLqFlVlP8P7tt7Niaqdq2bcglEWMu1Bj2b7ddIFs1cpPsA1zLxvZ6JqMJ5BGO+851WbYg/Q7qgvt6h0FbSZScNxTjH3qKpaPuYjiPBnJKlQv4nMXaixlXTI+itB8C44AHYqHEsfBBnWK5/GD5pUXCOxwcDgcDGkzhI+XZrO77/PUD7Rm57jbmLWuHY1Pd9K+vb/ijQhCFSDGXaixBBtTW7Y0BnNW3IvkObLJ9B1DWqBlxNsp3UB65MS7Eziv/QV8tvQj9p/0d5okBNj24y2MGeOmVy8v/fr5cFXtLgXhIMS4CzWS3FzYtMlB48Z+kpNhv2ML0z0jiAsk09J3it3yaJTYiLPansM3K78i55h/0q3FfjaMvI85c9xkZzs599xiydEqRJWo+tyVUg8opaYppeYopW6I5r6EusX06S58PscBl8wUz2N4HXm04RTcJNisztCufgdOataPIsdeVmf8h4uv3EX79j7Wr3cycmQcubl2KxRqM1GruSulBgB9gROBJOCeaO1LqHtMnGhu3dat/WQ7F7HI/SENfZ1o4uxJMbHj1z6+aR+Wbchhh2sRP6ZdyHnnjuGP31KYN8/FF1/EcfHFxSQnm7LB9gBJBiJUBdF0ywwCFgLfAGnAvYcqnJ6ehNttryMyIyPV1v2Xh+g6mMmTwe2GDh08jHI/RMDhZ5DzeTYxA0+cPZ3AkpPDP06dnMNYEvCy3j2Rn1Kv4OJzviY+3sX06U7Gjo3nmmvA4wGP9TiWt51onG+5typHrOoKRzSNeyOgFXA20Ab4Vil1tNY6bNi+XbvyoiilYjIyUsnO3merhnCIroPZvNnB4sUptG7tR3t/ZGXc/2jlHUjTglPYlDyDoopi/kaJX8aVv6xj3AVs8M5khfsHvvBdwpATP2TfvmQWL3YxerSPoUO9FMUb3bm54WvuVX2+5d6qHLGsKxzRrOLsBH7WWhdprTVQAGREcX9CHWHCBMsl06aY3z0PQsBB/6Knq71Pe2Vw4mZowShaePuxwj2GsYkXcsoZe2nZ0s/KlS5mzpTuM0LVEk3jPhkYrJRyKKWaAckYgy8IR8SECcYQert8RLZrIZ29V9LY391mVRUTRxLnF4ylrfdM1rrH8U3K2Zx+zg5SUgJMnuxi74ZWdksUahFRM+5a6++BecBM4DvgNq111Y0WEeokxcXw++9umrfdx5/pw3EHEjmp6BG7ZUVMHIkMLRiFKr6ITa6pfN9gMKcN3WyySX0/DF9xnN0ShVpCVPu5a63vi+b2hZpHeZmOglTUI2T2bBf79jlod8NzbHRu5oSi+0kNNK9KiVHHhYezCj/AQwoL495nUrtT6Xrizyyc0ooNv58K/X+wW6JQC5DYMkKNYvx4F6RuYmn6CyT7m9C7qGb2sHXi4ozC1+lVdCc5Ts26gf3xtJ7Flll92LFO8rEKR44Yd6FGMWGCG+dpD1EYyOOkokfwkGK3pMPGgYMBRc/Qt/Bh9rrW4b/8TEjO5rf/nk2g/FSwghARYtyFGsO2bQ4W7ViAv9uHdGrYhc7eq+2WdMQ4cNC3+J/0LXwYr2cnrmtPZ92iLNbMVnZLE2o4YtyFGsP48U4YdDc4Ajza98kqS5cXC/QpfpAMX3d8jRbB4DuZ8tlpUnsXjggx7kKN4eN530Cb3+jbaDD9W9gfHKwqceCgrW8wyf5M6PVftnmmsHrW0XbLEmowYtyFGsH2/TuZm3kXDm8CL57xtN1yooITNx28w3DjgbP+yuSve0rtXThsxLgLNYLbvn2QQFI2vfY9Qtv67eyWEzUSachJzU+CpByy2/6bjYtb2y1JqKGIcRdingnrx/H77s9hcy/uPemvdsuJOsdkHkt9RxYc8y5TfkuzW45QQ5FkHULMMGWqixETSyewKGIfXzX5O/jcpP72X056JHbjx1QVToeTQR36M2r5Z2zMfJddmwaTniWRO4TKITV3Iab5w/MvNu3fAJPv56xenYmrA6Pzp0x1sX5RG5L2doe24/lmjKfKUwEKtR+puQsxyybnVObFvUl9r2L3pIcY8n5xqfAFUzwHGzxPLbql23iOZTHzyWk2Cr/3dLvlCDWM2vMkCLUKLwX8nHArAK7v3yEt2cMpp+xnwXybhR0BwUxLkZLmbEJ89vEUtpjKZn1WlFQJtRUx7kJMMt3zDDlOTfudf2XlnydzxRVFxMfbrar6ae7qwSpmsCVuGiNG3B22jCTUFsIhPnch5sh2LmRm3POk+lsQP9m4YS64wJ7sSnaTkdIQZ87ReFv/xObc9XbLEWoQYtyFmMKPl//F34Lf4eXUvP+wamkaTZr46dOnbqYCcOCg4e5Twenjj7wP7ZYj1CDEuAsxxZy4V9nmmkOn4ivwLhtMQYGDCy7w4qrDnUVaNGoABWlsyvgAX6BufsEIlUeMuxAzFLKXqZ4nSAw04pTCESxcaCz65ZcX26zMXuIT/MRvOBN/yibm5/5itxyhhiANqkLMsN71G8WOXAYWPE/RnoasW+cgK8tPhw5+u6XZTiZdWM8opjheJN8zp9SyETONy6qiLFZC3SKqxl0pNQ/YY/1co7W+Lpr7E2oum/ZtZIdrEZm+XnTxXsOUhS7AQdeu4oYAaJpVzPptXSlsNIPi4uOJcybYLUmIcaJm3JVSCQBa6wHR2odQO/AH/IzfMA6ADH8XJrueYe7Ce3AluMnpOoIRMwtsVmg/TreP5B39yM1cyKbsHbRuULPyxgrVTzR97t2BJKXUL0qpCUqpE6K4L6EGs3DHArbnbSPD15XUQHNylnWmODeVxt3m4Iqr2/72UJolG4O+w7nEZiVCTSCabpk84Hngv0AH4CellNJah/3OTk9Pwu22t0tERkaqrfsvj9qkKzm59EikYl8xUzf/gcfloV3gdDxON9vmmnpAi+Pn4olzk5wc/jb1lBNnxhMXm01JR6qractCVm7uTXHTWfi9p5DgTgQ4cH4O53rUpnurOohVXeGI5lOwHFiptQ4Ay5VSO4GmwIZwhXftyouilIrJyEglO3ufrRrCUdt05eYWlvo9c+sMcotzOaFpXxzrEtm1NYO9G1pSv53GlbqdomLIzQ3fx72oOExsmTg3RcWx56evKl2pu/qwt9lM1uzYTLuGrYCS81PZ61Hb7q1oE8u6whFNt8z1wAsASqlmQBqwJYr7E2oYhb5CZm2dTrwrnmMzewOwdc7xADTpNcNOaTFLs7Qm4HeSE7fAbilCjBNN4/4uUF8pNRkYBVxfnktGqJvMz/6TfG8+x2b2JsGdQHF+IjsWdyO+fg71266wW15MUr/xfhwbT8bbaAF53v12yxFimKi5ZbTWRcDl0dq+ULPx+X3M3TaLOKeHno17AZC9oCd+r4cmx8zA4ZTkoeFwOKDe3t7s5nc27dlKh4bt7ZYkxCgyQlWwhaU5i9lfvJ9uGd1JcCcQ8DvYNrc3DlcxGd3m2i0vpslKbwheD7sS5tktRYhhIjLuSql7lVJNoi1GqBsEAgFmbZ2J0+GkV+ZxAKyb346CXY1o1GkhcUn5NiuMbVIb7cO5oT++9OXsL4q9Bj4hNoi05p4ETFRK/aCUukgpVQeSnQnRYvWeVews2MHRDTqS5jEJoOf/aBpSM4+Zaae0GoHDAfXyegKwaU+2zWqEWCUi4661flRrfTTwNHAKMF8p9ZpSqkc0xQm1k1lbTU+Y4zKNQd+3I42VMzqSnLmZlGYb7ZRWY8hqWA98bvYkSK8ZITwRN6gqpZKBNkBbwA/kAK8opaZqrR+Ikj6hlrEldwsb92+gdVpbMpIaA7Dg5+MI+F1kHjMTh+PgdSQ59MGkNsjDuekkfC0nsj3nNBo3SLRbkhBjROpz/wRYCQwAntBad9Fa/ws4A7g5evKE2sa87Sai4bGWr93ndbLw5+PwJBXQqHMNTpBqA/XzewAwfVGOvUKEmCTSmvsE4GatdW5whlLKo7UuVEp1io40obaxM38nOmcp6fENaJXWGoDVM49mf049ep4zFZdH4shUhqz0dHL8LtYVLASy7JYjxBiRNqjeVMawO4E5AFrrrdEQJtQ+Pl36Eb6Aj3q5vZg6zc2UqS4mfWn6uPsy51SwtlCWlPpFuLYcT2GDP9m2Q8YHCqU5ZM1dKTUB44pBKRWaMcELfBs9WUJtw+f38eHid3E748jwdwOgOD+R3auOIiljK8mNt9mssGZSL78HOUxlxsIcGGK3GiGWOKRx11oPBFBKvaK1vrN6JAm1kV/X/8KGfevp1qgH7s0m0UTOss4EfG4adfnTXnE1mKz0BuQEHKwtlF4zQmkqqrmfrbX+HpirlLq67HKt9UdRUybUKt5b+DYAPRsfw/LNZl72ou6An0adF9onrIaTkhrAk30cRRmzWLhmG13bZNotSYgRKvK5H2f9H4Dp3x76NyBqqoRaxerdK/ltw3iOb9rnQPfHwj312behDWkt1xKftqeCLQiHonVCN3AEePnnH+yWIsQQFbllHrH+H8h9qpRKA1porRdHWZtQS3h/0X8BuL7LTazYtRyAHYuN371RF+n+eKQc3zWd5evgt21jMJG2BSHyfu43KKU+UEplAEuAL5VSkmpdqJDc4lw+X/YpGYmNOavtuQAEAsYl43B5aXi01BGOlMwMF/G7erC/4SSmzNtptxwhRoi0K+RfgQeAy4CxQFfg/GiJEmoPX68Yzd6iPVzd+To8Lg8AedubkL8jk/T2y3AnSPLrqqB9Wkdw+nn+p+/sliLECBGH/NVab8F0tvrBSroh452FQxIIBHhv4Tu4HC6u7nTAs2c1pIpLpirp07UhBBxMzxtFgbwvBSI37ouVUt9j4sr8qpQaBcyKniyhNjB9y1QW71zI2W2H0jSlGQB+n4OdS7rhSsgnvd1ymxXWHuonpZBV3B9fsyl89F3YNMVCHSNS4349MAI4wcqw9AlwQ9RUCbWCdxa8CcCN3W45MG/jojYU7atHw6MX4XSHT3wtHB7X9LwUgHdmfGmzEiEWiNS4p2D87Lcppf4F9AT+WdFKSqnGSqkNSqmjj0CjUAPZuG8DP675jm4ZPejd5PgD85f8ZuKQZ4hLpsq5oc/ZOHwJrEv7jHXr7FYj2E2kxn00pm+7C3CE/JWLldDjLUDS6tRB3l/0X/wBPzd2vRmHFcc3Px+WT+2CJ203qS3E+lQ1qZ40eiScDY00L42SgWF1nUijQjbRWp9eyW0/D7yJ6WUj1CHyivP4eMn7NEpsxHntLzgwf9w4N0V5CTTrMwOHQxJgR4PbTr6IGyd8yTerRvJsYVfi4+1WJNhFpMZ9nlKqm9Y6ogAWSqlrgWyt9c9KqYiMe3p6Em63vUkZMjJSbd1/edQ0Xe/MGcnuwt08dPJDtGiacWD+d1YvvWY9FuOJizhPTKWJ5raPhGjpSk42283ISOXqBsO4/deG5LcfybgJL3Dd1YfOiFnT7i27iVVd4Yj0buuCMfDbgAKMSyagtW5bTvnrgYBS6jSgB/CRUurcQ4UH3rUrL3LVUSAjI5Xs7NhLNlzTdPkDfl6c+hJup5uL2lx5oMyuXfDjjylktNlCXPpmiqIUut0T56aoOPbC30ZTV26uaZgOnutz2lzA6HVv83+v/Y/FC844KLvVffcVATXv3rKbWNYVjkiN+7DK7Exr3S84rZSaCNwicd/rBuPW/cyynKVceNQlB7o/Anz7bRzFxQ46DvgTSckRXf5y3BWMXvc2uzu8xebNg8jKEhdYXSTSBNnrgBOBvwDZQH9rniCU4tW5LwFwe8+/l5r/5ZduHI4AHftLL5lo071xT45KOg46/Mi0ZfKY1lUijS3zDGZ06vmY2v51SqkXIllXaz1Aa73s8CUKNYXpW6Yxc+t0zmg1mE4NOx+Yv369gxkz3Jx4oo/URnttVFh3uLPPTeAIsDbjLfbKKa+TRNoVchBwFVCgtd4LnA6cGTVVQo3k1bkvAnDHMf9Xav7XX5tGvQsvFIdMdXFu+2F4ijOg53vMnl9otxzBBiI17v4yv+PDzBPqMEt2Lmbcup85vmkfjm96woH5gQCMHu0mPj7A2WfHXkNnbSXeFU9333WQuIsF/tEUFdmtSKhuIjXuXwCjgHSl1N+BP4DPoiVKqHkEfe13lPG1z57tZMUKF0OGeElLs0FYHaan7yYIOPEe8x8WLoo4RqBQS4j0iv8AfAfsAE4GHtZaPxU1VUKNYuWuFXyz8ks6NujMaa0GlVo2cqRxyVx2mbhkqpu0QAvaFJ4DTecxc/NM/PKtXaeoKIdqY+BLoDOwAvACA4FEpdRkrbXkRxN4Yfaz+AN+7j3uAZyOkvrCE094GDUqjtTUANOnu5g508UUj70D1eoavf23sYax5PZ8lpUrv+Koo8TC1xUqqrk/DUzGhB84QWt9AtAYmA+8Em1xQuyzPEfz9YrRdG7YlSFtzy61bMUKJ0VFDjp39uEUr4AtNPefTOOCPqC+Y9rqRXbLEaqRih65vlrrB7XWB76prekHMZEhhTrOC7OfIUCA+3o/WKrWDrBokamld+kioX3twoGDkwL3A5B99DNs3nzIeH9CLaIi4x42p4vWOoD0lqnzLN25hDErv6ZbRg8Gtx5SatmaNQ42bHDSooWf+vXt0ScY2vjOID2/J3QezRS9wm45QjVRkXE/1LhlGdNcx3k+WGs/7oEDYX2DjBplGlKl1m4/DhyczP3gCLCuxQjWrpXae12gotgynZVSq8PMdwBNo6BHqCEs2LaA71aNoWfjYzi91eBSy4qK4JNP4vB4AtKAVw1MmWrcXyMmesot08F3Dil5ndjf7VNeeO8BXn0sq7rkCTZRkXE/qlpUCDWGETNND9jv144BoG299jw36+kDy+/r/SDffedm+3YnvXp5iTt0xFmhmnDg5GT+wU/Oa/hy+wge3/0KGRkVryfUXA5p3CU4mBCObblbWbZjGc2Ss2id1qbUshEzn+KTl/4Kjiz8vV9himeXTSqFsnT0X8jv+58lr8tHvPLxPbz2WCe7JQlRJDazGggxzdTNkwE4Mevkg3ztm5a2ZOvyFrQ7fgkJ6WLYY4lpnqdo6u3KKucS3lr6OOnjOlFQVJJH4b7eD9qoTqhqxLgLlWJL7hZW7VlJy7SWtExtddDyGV8MAODYYX+wXqIRxhwZ7nas29sOr/qCT997iSZtQ4KK9bZPl1D1iHEXDiLoVw/H1E1/ADCg9YCDau3bVzdh9ayOZHVaS/POa1k/TUajxhoOHLTwn8Qa5yo2uCeQGehyUKYmoXYg4waFiNm8fxNr9q6mRWpL2qS3OWj5jNGnAHD8RRPFYMQwmUlZuHI64m33LTuyxQTUVuTKChEzZbOptfdtdtJBy7avboL+oxuZ7TfS5lhd3dKESuDAQbP808ERYH1gmt1yhCghxl2IiE37NrJu71papbamRWrLg5ZP/vgMAE6++heptdcAmjVKwbmtJ0UtfmV3XuwlfRaOnKj53JVSLuAdQAE+4Dqt9apo7U+ILjO2mhpen2YnHrRsw6LWrJ7VkeZdVtOqpwxvrwk4HQ6a7B7K5sx5rCmaTc+kU+yWJFQx0WxQPQdAa32iUmoA8CIwNIr7E6JEdn42q/esollyFs1TW5Ra5vc5mfDmuQD0v+4nqbXbyBTP45Uq36a1i80bT6Cg+VT25Uuf99pG1NwyWusxwF+sn62AbdHalxBdZm+dAUDvkPR5Qeb/rzfZa5vS5fRZNFUbq1uacAS4PV4abL4EgNVFc21WI1Q1Ue0KqbX2KqU+BIYBFx6qbHp6Em63vV3nMjJSbd1/eVS3ruTk+APTewv3sjRnCY2SGtGtWedS3R8d3nSmfHIG8ckFDL51Yqn1ADw2hR7wxMVmD99Y1NWufYCctf3Ja/07qwqWcEKL4+2WBMizWBVE/W7TWl+jlPoHMEMp1UlrnRuu3K5deeFmVxsZGalkZ8dew5IdunJzSwa2TNowGX/AT6/Gx5GXV5JlOTk5nl/e6kfBvkROuel78OSQW+bKFhVX/8vaE+emqDj2EnHHoi5PnBt3yk5S5l7L/ta/c9fXj/Pt5Z/bLUuexUpS3gsnam4ZpdRVSqkHrJ95mPjvEv+1BlHgLWBB9nxS4lLo2KBzqWVbVjRh/v9607DlNnqcJd3pajItWvpgQx+m7/6BZTlL7ZYjVBHR7Ar5NdBTKTUJ+Bn4u9Y6bPIPITZZsGM+xf4ijml8LG5nyUdeIAA//XsIBJwM/Mt3uNwS1rcmU6/VGjzzbwFgxNSX7RUjVBlRc8tY7peLo7V9Ibr4A37mb5+L2+mmW0aPUsuWTuzBhkUt6dB3Ea16rDoQT1yomTgckNUolTXbO/ND4As27Ptn2LEMQs1CBjEJYVmzZzV7ivbQsUEnEtwJB+YX5Xn4/f0zcXuKGXDDDzYqFKqSjC4LSJx7HwGHj9fmvGq3HKEKiL3meyEm+HO76RrXo/ExpeZP/+IUcnPS6H/NROpl7rZBmRANXJ5iOjZrztzdLflw8fskxsWTFJdUqoyEBK5ZSM1dOIjdBbtYs3c1zZKzyExqUjJ/SwNmjzmJtMa7OPGyyTYqFKLBMUNmw7S78TuKmLt9tt1yhCNEjLtwEH9mzwMOrrX//v5g/F43/a77ibj42OrWJxw59Zvsoq2rL+Q2Yu6WeRT5CiteSYhZxC0jlCKvOI95WxfgJomduhNTMI2le9e3YsXUrqQ2X8cO5xImTnTb0o9diC7HnjWP1T/fSdHAh5mf/SfHNYmNQU1C5ZGau1CKb1d9g89RQKa/B07r3R8IOFg7/kwAWp0q8WNqMy26rabBlguhMIVZm+fg88vQlJqKGHehFJ8v+wSAxr6eB+btWNyN3C3NadhpPqlZEj+mNuNwQJ+hf8Kcv5Dn38vSnMV2SxIOEzHuwgHW7FnNtM1TSPO3IoH6APiK41j/2xk4XMW0OmWcvQKFakGdvIB6a64Cv4sZm+YQCATsliQcBmLchQN8oU1ckcb+bgfmbZnZl6J99Wjaeyrx9XbbpEyoTpyuAH3OXgGLLmVX8TbW7F1ttyThMBDjLgBmROoX+nOS41Jo4D8agKL9KWya1g930n6y+k6yWaFQnXTs/ycpy24CYPqGOTarEQ4HMe4CAFM2/cGGfesZ2m4YLjwAbJh0Kv6ieFr0G487XrrF1SVcbj99B+6EVaexuWA1W3O32i1JqCTSFVIAShpSL+14JU8vmETu9ky2z+9FYqNtZPaQmltdoGyMIH/yfNy/34q33a9MWzfPJlXC4SI1d4F9RXv5YfW3tKnXluObnEAgAOvGnwkBJ61O/R8Op0R9rIs4XT5atfHC1m6s2r+I9XvX2S1JqARi3AXGrvyGfG8+l6orcDgc7F7dgT1r2lOvzQrS20nC67pMRueFpK2+Hpw+Hv/1dbvlCJVA3DICI5d9igMHF6vL8Hph3a9ngsNPq1P/Z7c0wWYcjgCDToljdHYW3xZ/SNYUD0lxiWHLSmCx2EJq7nWcVbtXMHPrdPo1H0BWanM+/jiO/J2Nadx9DsmNJae5AK26bqDh5ksJxOUxYc5mu+UIESLGvY4zapnp237p0Vewaxc8+6wHl6eAFv3G26xMiCWG9K8PBWnoot8pkHxqNQJxy9RhfH4fo/RnpHrSOLPN2TzxSDw5OU5aDpyIJ2W/3fKEGCKzRS6ZC85iW+bn/G98IeedFX9Q75oRE00X2vvuKwq3CaGakZp7HWbSxolsyd3Mee0vYP2qFN5/P462bf00PU4SXgsHM+SkJuBzs9L1I3u2pdotR6iAqNXclVJxwHtAayAeeEJr/W209idUnlH6UwAuUZfz0B3x+HwOHnssn/8skUiAwsE0rBdPlq8vmxpN4odvU2nVOc9uScIhiGbN/Upgp9b6ZOBM4LUo7kuoBCNmPsVjU//FtyvH0CChAe+M3sSkSW5a99LMq/e43fKEGObU7h0A2Jz+JXvWtLVZjXAooulzHw18GfL7kKl70tOTcLvtTf6QkRGbn5pVrSs5OZ6lmxfiDXjp3rgXk/59Dk6Xj7P+No6UlHg8cZFtxxMXm002oityItGUnFxSpk1yc5qv7sTG1pNYNfoBerfbiNNlBrklx8UDkJERf8S66sqzGE2idrdprfcDKKVSMUb+oUOV37XL3k+8jIxUsrP32aohHNHQlZtbyNzN83DgIHfy1eza3IBeQyeT2HATublElGHJE+emqDj2Uu2JrsiJVFNubmk3XZ9W3Ri9fAmFnd5l/fTraNZ7qilXZOIPZWcfWYNqXXoWq4LyXjhRrUoopVoA3wCva60/i+a+hPIZMcJT6vc4zy62eDbTIvEo5n12EYlp++lzmXR9FCKjZWorGsU3ZUfHr9nw1qM06jQfT0qu3bKEMkTN566UygR+Af6htX4vWvsRKk+2awEA3pk3UlwQT79r/0dCinReFiLD4XBwfLNjwenH3+t11k8cZLckIQzRbFB9EEgHHlZKTbT+wo9bFqoNP16ynQtx+pLYMvYOmqr1dDltrt2yhBrGUelH4wmkQc93yV7Zkn2bmtstSShDNH3udwJ3Rmv7wuGx1vUrxY79uBZfh8Pv4bS/jsHhlDRqQuVwOV008x3P2rhxcOwbrPn5evyXg1NGzsQMcinqGIvcHwPgm34rPYZMJ7PdFpsVCTWVxv7uuALxOPq+RO6OBixaJOYkloi9vllC1Mgnh5Xu7yC7I67dR3PilS/aLUmoAZQNMxDEhYtMfy82J0zF0etd/vjjNo46SkIPxAryqq1DLI0bhd9RBPOuo/XAX6QRVThimvqOxRFw4RzwBPmFXqZPt3esilCC1NzrCAECzPa/Cz43ydkDyOg/ttwamSBEiodUGvt7sC1xDgnHfca82Vezdq2D1q2lHcdupOZeR9gYmMXexEWgh9Ku/3QcDrsVCbWFLF8fHAEnrgFP4fP7eeqpIx+hKhw5YtzrCBP3vg9A+z03SBIOoUqJpx6N/F3JTVxOvb6jGTMmjtmzxbTYjVyBOkBO7l62ZYzCsbs1g48aYLccoRaS5euLI+DEcfJTQIBHHkkgIJ4ZWxHjXgf438YvwZNHu73XkRAvl1yoehJpgPJeyO6EhfS6/BtmzXLx/ffSpGcn8qTXchYudLK56bvgdzEw/Sq75Qi1mBOK74eAg729huNy+3j88XiKpGekbYhxr8UEAnDPC4uh2Rya7R9CmqOZ3ZKEWkwjfyc6ei9hxb6F9Lt1JGvXOvnggwjjRwtVjhj3Wsyvv7qYF/cmACfE3WCzGqEucGLRw7idbla3eoTUesW88EI8u3fbrapuIsa9luL1wsPP7oKun5JW3IE2vjPsliTUAeoH2nFFx2tYt38VA+58l127HLz0knSNtAMx7rWUjz+OY3WDt8FdRG/fHTjkUgvVxP/1upcEVwJzUp+kees83n03jrVrZWBFdSNPfC3k8Ykv8Miz+6D368Q7E9nr3MAUz+MH/gQhmnz8emu65v+VzbkbSR34GkVFDq66SqJ9Vzdi3GshM74YQEHrMZC0gx6Ne+DCU+E6glCV9C66h4RAA1a1eIKMthvR2iUDm6oZOdu1jA0bHMwe2xfnSS/gdDjp2fgYuyUJdZBEGtCv8HGKHPtIPO9uAO6/PwFvbKWRrdWIca9FjBjh4corE/G3Go+/gaaBtzPzZ9e3W5ZQR+nqvY4mvuNYn/YFrQaMY8ECF2+8IV+R1YUY91rE5s0Oli514jztIQCa+XvbrEioyzhwclrhK2Zg00l30CizgBEjPCxdKmanOojqWVZKHa+UmhjNfQgGvx8mTHBDu1/wN5lLuv8okgNN7JYl1HGa+I+hR/Ff2OXW9P/n0xQWOrj55gTy8uxWVvuJWvAHpdR9wFVAbrT2IZQwerSbrVsdJNz+LwqAFr5+dksS6hjl9cSKI5kUfxZjdz/F0FvOZOybJ3D//Qm88kqBhJ6OItGsua8Czo/i9gWL/fvhiSficXYdTUGjmTTwH01yINNuWYIAgJtEzix8G6/fyyJ1NV177WHkyDjeeENCE0QTRyCKcTmVUq2BkVrrEyoq6/X6Am63ZAY6HP75T3jq2UIS7ulIUcJGjuUWEmlgtyxBOMAAhrOvz928OP1Fzml7MbPvH8nWLQ4++QQuv9xudTWesN8/MROTc9cue51wGRmpZGfvs1VDOCrStXixkxEjkkgb9BJ7E9fQq+gOXKRRRHT7nHni3BQVx16/NtEVOdWpKbeokLu7/5PJa6fy3eovuPGpLnzxt4e4+mrIyytg6NASHTX1WbSLjIzUsPOl2boG4/XC3/+egDdlLQXHP0aiP4MTih6wW5YghCXOFce7gz6iaXIz/rv2X9z82gckJsJNNyXyxhtxktyjihHjXoN588045s930uSGWykK5HNK0QhxxwgxTWZyEz4/+ytSPWm8vPZm/vHeaDIz/TzySAIPPRSPz2e3wtpDVN0yWuu1QIX+dqHyrFrlYMSIeJIHvM7W1J8Z0GIgHZdearcsQaiQTg078+mQL7j0+wt4bNmVPPluEe/dcxnvvONh5Uono0ebciNGhB/wdN99kgEkEmLG5y5ETnEx3H57IgWpS4kbcA/p8em8csrrfLxU+pUJsUtpYz2Ac5xj+DpxGPfNvop7nllF038/yG8T4jjmGHjnHXEqHClyBmsgTz4Zz5zFe0m5cSjFFPDiKa/RNEWyLAk1ixb+k7ksfwLNkrN4/s/hZN50PXfft48NG+Ccc5KYP98pfvgjQGruNYyxY928/hYk3nAx++NX8ree/8dZbc+xW5YgHJJDhZo+t90wpm+Zwkj9CT1bL+Gdr0dy7/XtGTcujs2bfZx2mpc46RJfacS41yDmzXNy+9/icFxwOfnNfqWd9yzcfzzOiD9kfIBQc5k/ux79i8ZTFH8H87Z/yq3betLn6hdYNvJGFi92sX27g6FDi6lf326lNQtxy9QQtHZy6RUuCs+6kkCn0TT3ncjZBR/jRAy7UPOZ6XmBeoHWdPCehx8vk9JvoeimTjTq8x3Z2U4+/tjDqlViriqDnK0agNZOLriigF1nngtdRpHl68uw/G+II8luaYJQpTTyd+ZYbiHN34rd7uXknHEBDa+/Hq/fxzffxPHHHy6JCR8hYtxjnGnTYMj1K9h+3vHQbhxntBrMhfnfE0+a3dIEISokUJ9O3ito5z0bBy52tnyftHt6kdzlN2bMcHPJJYlkZ0vPsIoQ4x6jBALw4ccOTv7Hc+y75HhosIo7j7mbD8/8XGrsQq3HgYPG/u70LL6FDF83dsUvIPfCgSRfP4w/lqzmtNOSJG1fBUQ1cFhlyM7eZ6uQWIobsWOHgxsfm8PU+n+DJguo58rgjcGvc1qrQQAMfflZmxXGZqwUEF2VIRY1QXhdbb1nMjH+Xja5puEIOAksuBL3lH/y+F0tuf764moJHRxLNiKUjIzUsEcvxt0iFi5cQQG89vEWnp/3OP4unwDQqPB4Wjv6xlxtvSYZhlggFnXFoiYIr+vEoocJEGCFawwrsh5nac4S8LtgyYWc4LyVroOmUi9zb7nbvK/3g0esKxZsRDjKM+7SFTIG2LcPXvsghzcXP09+p7ehSxHJxa1p7e9PI0/rmHwABaG6ceDgKN8w3r7kTL5fNZZnp49gRZdRTGcU02f1oK3zJE7r7yEtWfK0ghj3aiV0+LXPB2vXOlmwYQ3rsl7C2+UD6FZAmq81xzY+hh4t2jNtmozcEISyPP9cAnAJ53IxG5yT+L3gDbY1/pbVzj95e5Gb+vm96Nkmi6ObZZEcl2y3XNsQ415N+P3Gl75hg4O5G1awpzAHf+9/w5CvwOknvjiTPpmn0jOrEy6n9F0XhIpw4KClvz9XefozLu8h1mzKZW+9aexuMoPfsuG37Q4S8trRs00LlrRbTMcGnXDUobx+YtyjxJYtDubOdTFvnpN581z8+aeLfQX50Gk0nPwWtJgGQEJhC049+hhUA4XTIa3/gnA4JHni6dwmHl/ReWyadxPbcrfjbf4bBS0nM237SgaM6kNDV0vO6jCIIe3P5KSsfnhctdt9I8a9Cti1C/7808X8+SXGfOvWoKEOQOZC6p31Ds62n+D37IaAg/r+djTzHU+aozU52sE0Ow9AEGoJLk8xLTtvpkUA8rZfRPbEe9gT2ERe49/Z2eEnPlr2Dh8te4ckZ33ObDuYczsMZUCLgSS6E+2WXuVIbxmLQ7WEh/rK/X7Yvt3Bxo1O6tcPMG+ei3XrSte4k1P8NDh6PnT+kt3NvmJfvDbz/U3p6r2GQvaRQP2IdNWkHg2xgOiKnFjUBFWvq28fH2cmPsTXYwN8OWMG2xt8Cx2/hnobAYh3JHNGm0Gc2/5cTm11BilxKWG3U9N6y4hxtyjvwvl8cN998WzY4GTDBmPUCwtLzmV6eoDu3X106rmL+KMmkZ3yGz+s+IVdzhUAuAOJtPENplPxpbT1DcFF3CEj5JWlrjyAVYXoipxY1ATR1RUIwP4tWexY0pm9uUXkNR0Hnb6CBqsAiHMkMKDFKQxqM5gzWg+mSXLTA+uKcT9MYsW4FxfDwoVOpk93MX26i2nT3OzZU3Lu6tcP0LxFMQ06LMeVNY99qX+yyTWZrc45BBwmR5g7kEhb35kc5T2ftt7BeChdExDjHj1EV+TEoiaoPl2BgIP9m5qTvaQLO3eCt/XP0PEraLzkQJmM+Ga0b9CGu4/7B4M6n0L+ntiwl6FUu3FXSjmB14HuQCFwo9Z6ZXnl7TDuxcWwerWTpUudrF2byKRJXubMcZGfD3hyIWkHjTusJa7xGhKbrsHZYA374leQ7VyI15F3YDvOgJsm/mNp6RtAS98A1rl+w1lFzRl1/QGsLKIrcmJRE9ijK+B3sHdDa3Yu6cKenAQKmkwE9R20+h1cRosj4KK5uyvt0zrTKUPRpYmibaNmZKU1pVFiI9s6RNhh3M8HztVaX6uUOgF4QGs9tLzyh2vc9+6FefNc+Hzg9UJREUzZ9zk7izdR6CukyFdMbmEh+UXF5BcVUVBcRG5BMXmFheQXF4OrEFxF5s+zH3fqTvyJO/E7CsPuzxFwkRhoSHIgk6RAJsmBTFICTXERfzjyK0QewMohuiInFjWB/br6nuBn+5omrJvXgbVLM9havJKihnOh5WRoOgfcYXK4+l0485rgLErH5U3B6UsiMz2JHp0TSY1PJt4Vj9sZR5wzjjinG5fTTZwzDrczjoYJDbngqIsPu/eOHSNUTwL+B6C1nq6UOjYaO/nb3xL48ceQwT711sFdN5Yu5ADirb9ycARcOHHjIpH4QAPiAkm4A4l4SCMhUI/4QH3zR1qV1coFQYg9HM4Ame22kNluC70xfvrcXfWZ8sut7J+TRK53PwWu7RQnbMbn2YUvfgf+xGz8SdvxJ27C69kPrmLWAmtXRbbP9ukdOK7J8VV6HNG0UmnAnpDfPqWUW2sd9pVc3tunIn74oeycVkDs+cUEQajhXGW3gMoRTSfRXiA1dF/lGXZBEAShaommcZ8CDAGwfO4Lo7gvQRAEIYRoumW+AU5XSk3FeL2vi+K+BEEQhBBipp+7IAiCUHVIpCpBEIRaiBh3QRCEWogYd0EQhFpInRiNE0koBKVUEjAOuEFrvcyaN4+SvvprtNZV2ihckS6l1GXA3wEfsAD4q7Uo4rAO1aVLa+2P5vmKQNMFwP2YQQ5va63/W9kQGNWly5pv670VUu5tIEdrfX8snK9wuqzfdj+L/wfcAGRbs24GVkRyLHZRJ4w7cB6QoLXuY3XLfAE4EArBGj37JtA8ZF4CgNZ6gB26lFKJwBNAV611nlLqc+BszDUr91js0qWU+gWier4OpckFPAMcC+wHliilxgD9ylvHZl37wb57K4hS6magK/B7pOvYocvuZ9HiGOBqrfWcEJ3nV7COrdQVt0ypUAiYhy2UeGAYsCxkXncgSSn1i1JqgnXxqlNXIdBXax2MUOYGCipYx05d0T5f5WrSWvuAjlrrPUBDTNfb/RUch5267L63UEr1AU4A3op0HRt12X6+gF7AA0qpyUqpByJcx1bqinEPGwoh+ENrPUVrvaHMOnnA88Ag4Bbg09B1oq1La+3XWm8DUErdAaRg3EaHPBYbdUX7fFV0Db1WTWo+MAkormgdG3XZem8ppZoCw4HbIl3HZl22ni+Lkda+BwInKaXOjmAdW4kZIVHmcEIhLAdWaq0DwHKl1E6gKVD2JRA1XZYfcARwFHCB1jqglKqOsA6Hoyva56vC49Zaf225PT4Aro5kHZt0fYa999ZFQCPgR6AJpla8LJJjsUnX59h4vpRSDuBl6wsMpdQPQM8KjsV26krN/XBCIVyP8aGhlGqGeUtvqWZdbwEJwHkhbpDqCOtwOLqifb7K1aSUSlNK/a6Uitda+4FcwB/Bcdily9Z7S2v9b611L8uH/Qzwmdb6g0OtY7Muu5/FNGCRUirFMvQDgTkVrGM7dWKEakhLeDdKQiEcA6Rord8OKTcRuEVrvUwp5cHUtFpiejr8Q2s9tbp0AbOtvz8oCXP5CjC27DrB3j026/qBKJ6viq6hUuovmN4MxZgePHdYOmw7V4fQ5cLGe6vMPX8tcHSZ3jK2na9ydNn6LFrX8Srgb5g2p/Fa60eq43wdCXXCuAuCINQ16opbRhAEoU4hxl0QBKEWIsZdEAShFiLGXRAEoRYixl0QBKEWIsZdEAShFlJXRqgKVYhSKrT/bAAoAjYCj2mtP7JHVdVhjXfoD/SxYoaglBoGLNFa68Pc5inAo5iRjV5MHKOntdbfVoloQSiD1NyFw8WPGVD1A6CBdsAHSqkBdoqqIiZjji0HQCn1H+BrIP1wNqaU6oYJMHUSJsbMAqA3MFYpdWZVCBaEskjNXThcirXW5wV/KKVGARcD1wATbdJUJWitHyozq/MRbvJqwAPcp7V+DkApdQvwBnAj8NMRbl8QDkKMu1BV/IEx7s2glGvjQeBuYKnW+mSlVEPgOeBcIAkTn+MerfV8a70PMC+Ii4C7MKFWZwLXBxMhWPFFXgdOx0QM/By4V2tdaC2/FxPBLwvYAXwCPGglFLkWeB94GDOc/TJMnJLHtdYfl9Hex9pOf+sYpymlHgUuBzoAXbTWi611vgfOAk7TWo8vc26Krf/nKKV+1lovAD4EFgH7goWUUkdjQjmcjIlD8wNwZ0jAqh7As0Bfa5vfAndrrXdaywPAEuucXgq8bg3fPxd4GmgPrAQe1lp/ba3T1NrnACAZ4y56UGv9M0KNRtwywhFjxdgIuhc2lln8KMbgTLdihIzHxO3IwbgoTgMmKaXal1nvPYzxX48xdl8rpZxW4KZvMEkRFmEM8x0YY48VinWEte4EzD3+D8yLIpR7gVOAtRhD/aZSql6Yw5sL7LSmJ2GM3+fW7/OsfSYBpwLbCf/VMgpjjE8G5iulVmGMbU7ISy0R+Bk4A1iMeSldY50HlFLtMC/QM6zjDi7/1TqvQTpiXrILgVlKqa4Yl1IrTPKLTGC0Uir4wnoV8yJdD0zFxEkZo5TKCnMcQg1CjLtwuMQppcZYNdYlmOh4PuCdMuWe0Vr301rfizEi3TFGspPWug/GyKUBD5RZbwImcFMXjLHriqlBn4LxV3+utT5ea90NY5SuVko1ANpa67+LeQH0xaQE/K3M9ndZWnoC6zAvg45lD1Jr/W+MMQUTsGokJmQvWMYdY3ATgC+tBB1lt/GnVXa5NastcCfG0F9mzbsC8yXxmdb6OEvbz8BOpVQcJlVfCvCUdd46Yc5jD+CSkN05gPO11idqrb/CvMRc1rwzrPPhtPYf1OLFfBmdDlwI/AUTIEuowYhxFw4XJ8Z4DsGkJ5wJnBPsXRLCtJDpYAadz0PiXn9UZlmQ77TWAa11EfCLNe8oSvzflymlApYroi/GxdgDU0vdDPwT83XwGib7UdlwrH9orfMsHcG8l/EVHjVg9ZiZC/RSSjWnJLXaqEOs86PWWmHcTP/CfAG4gZetL59OVtHxVvkirfVgrfVftNbFlJyfj63lXkq+IMqeu9BzHjxfP1vnKtjbJ5g16DWM8Z+vlFqD+QJbrrXeUcFpEGIc8bkLh0uh1johgnJ7Q6b9YZY7rP9lw5N6yha01o+zppdRYqiC5GutNyqlOmNcE4MxeVTPsv7ODymbFzIdfNE4iJzPMF8Ww6xtb8b0sjkIpdSLmC+Qq7XWc4G5SqnnMS6lxkBGyL7dIeslaq3zrZ+RnjtvyDpQcr5+AULn5wNord9TSs3E1P4HYsIT36yUOk9rPbacYxdqAFJzF6JNqJtinvX/0pB0ZFdZ/0NrmwDDLB97HMafDcaYL7GmV1q9dYZZ8ycDC5VSVwL/BVZprc8HWlsaTi+z/crEug4a1tDnZaQ1/wGMcR5tJeQIRytr/3eGzGuKceXkY1xEweM6A8A67sVKqc2Wuyl47q6ylrsxjaZQ+tyVdQsFt/uOdb4exLQzjFJKeZRSrwCPAE9qrU/ExCwH0xYi1GDEuAvVyRcYQ9wfWKKUmoYxjnswvvdQBlLSJ7wbxg3yB6YGugw4WykVXH4fcLHWej/GWF6A6UP+P2AWxu0w6Qh0B10Ubyml7gbQWm+yttnUWlauSwZ4CtOger9Sar5Sapx1bPHAa5br6VNM7f8CpdRsa3kbYLbWOgeTQzQfeNA6b0swXyVzMee1PF7DvMg+VUqNx5zDO4Esa79tMX72RUqpnyi5DkdyvoQYQIy7UG1YKfn6YbLqpGMaDX8FTtZary5T/CFgN8bA/Q5caPng/Ri/8FiMYWoBfIXpWonViHgtsMbaV32MX/+6I5D+MrAaM1ArLWR+0Oe9Hijb1nAArfUcTFfDXzC5QU/G9Cp60PoLnptBmIbkzphcoh9iesSgtV5krfcLxsXTCHMeT7d88uXtezLG5bIcM4hqP6Zh+D9WkSsxaRMTLY1bgNu11qPLPx1CTUAyMQkxRUg/98usnikxi1LqZUwt+AWt9T02yxGEUkiDqiBUEiuf5iWYLwgfxscvCDGFuGUEofLEYxp5NwA3xlJSZEEIIm4ZQRCEWojU3AVBEGohYtwFQRBqIWLcBUEQaiFi3AVBEGohYtwFQRBqIf8Pp1NmkNxHjj0AAAAASUVORK5CYII=",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Hide Cell\n",
"sns.histplot(ps[df_categ['intervention']==0], kde=True, stat=\"density\", linewidth=0,color='b')\n",
"sns.histplot(ps[df_categ['intervention']==1], kde=True, stat=\"density\", linewidth=0,color='g')\n",
"plt.title(\"Positivity Check\\nPropensity Score Dists\", fontsize=20, fontweight = 'bold')\n",
"plt.xlabel('Propesnity Scores', fontsize=14, fontweight = 'bold', labelpad=5)\n",
"plt.legend(['Untreated', 'Treated']);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the second experiment, we specify a bad exposure model. Instead of using logistic regression, we simply sample a uniform random distribution: \n",
"\n",
"$$ \\hat \\pi (X_i) \\sim U(0.1,0.9) $$\n",
"\n",
"As we can see from the results below, the AIPE method is effectively stable, estimating a slightly lower ATE of about ~$0.38$. The standard deviation also increases slightly. On the other hand, the IPW method does far worse than the naive method, which again makes sense as we are feeding it random noise for the propensity scores. This is the first example of a doubly robust method showing how, since the outcome model is correctly specified, the estimation is still robust even to random noise for the exposure model.\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Hide cell\n",
"bootstrap_sample = 100\n",
"AIPW_ates = []\n",
"IPW_ates = []\n",
"naives_ates = []\n",
"for iSample in range(bootstrap_sample):\n",
" df_bootstrap = df_categ.sample(frac=1,replace=True)\n",
" ate, ps, mu0, mu1 = AIPW(df_bootstrap, X, T, Y,true_mus=False)\n",
" AIPW_ates.append(ate)\n",
" ate, ps,_ = IPW(df_bootstrap, X, T, Y)\n",
" IPW_ates.append(ate)\n",
" naives_ates.append(naive_ATE(df_bootstrap, T, Y))\n",
"\n",
"sns.histplot(AIPW_ates, kde=True, stat=\"density\", linewidth=0,color='g')\n",
"sns.histplot(IPW_ates, kde=True, stat=\"density\", linewidth=0,color='b')\n",
"sns.histplot(naives_ates, kde=True, stat=\"density\", linewidth=0,color='r')\n",
"# plt.vlines(np.percentile(AIPW_ates, 2.5), 0, 20, linestyles=\"dotted\")\n",
"# plt.vlines(np.percentile(AIPW_ates, 97.5), 0, 20, linestyles=\"dotted\", label=\"95% CI\")\n",
"plt.title(\"ATE Bootstrap Distribution\\nRandom Outcome Model\", fontsize=20, fontweight = 'bold')\n",
"plt.xlabel('ATE', fontsize=14, fontweight = 'bold', labelpad=5)\n",
"plt.legend(['AIPW', 'IPW', 'Naive'])\n",
"\n",
"Results = {\"AIPW\":{\"Mean ATE\":np.mean(AIPW_ates), \"Std Dev\":np.std(AIPW_ates), \"[.025\":np.percentile(AIPW_ates, 2.5), \".975]\":np.percentile(AIPW_ates, 97.5)},\n",
"\"IPW\":{\"Mean ATE\":np.mean(IPW_ates), \"Std Dev\":np.std(IPW_ates), \"[.025\":np.percentile(IPW_ates, 2.5), \".975]\":np.percentile(IPW_ates, 97.5)},\n",
"\"Naive\":{\"Mean ATE\":np.mean(naives_ates), \"Std Dev\":np.std(naives_ates), \"[.025\":np.percentile(naives_ates, 2.5), \".975]\":np.percentile(naives_ates, 97.5)}}\n",
"\n",
"df_results = pd.DataFrame(Results)\n",
"df_results.T"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the final experiment, we show the impact of a bad outcome and exposure model: $$\\mu_d(X_i) \\sim U(0,1), \\hat \\pi (X_i) \\sim U(0.1,0.9) $$\n",
" \n",
"In this experiment, we see that AIPW performs very poorly, vastly over-estimating the ATE. In this instance, either naive or IPW would perform better, although the naive without any consideration for random models does best. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Mean ATE
\n",
"
Std Dev
\n",
"
[.025
\n",
"
.975]
\n",
"
\n",
" \n",
" \n",
"
\n",
"
AIPW
\n",
"
1.054
\n",
"
0.037
\n",
"
0.986
\n",
"
1.132
\n",
"
\n",
"
\n",
"
IPW
\n",
"
0.575
\n",
"
0.032
\n",
"
0.515
\n",
"
0.630
\n",
"
\n",
"
\n",
"
Naive
\n",
"
0.475
\n",
"
0.017
\n",
"
0.439
\n",
"
0.508
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Mean ATE Std Dev [.025 .975]\n",
"AIPW 1.054 0.037 0.986 1.132\n",
"IPW 0.575 0.032 0.515 0.630\n",
"Naive 0.475 0.017 0.439 0.508"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAE1CAYAAADJbraRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABVrklEQVR4nO2deZwUxfXAvzN7n7DAAnJf+lS8AFG8UdQoSgSNmnhFY2J+xiRqNMYYTdAYjcQYo8ZbNEZNvFHx5hBBvEEUkFJu5FzOhQWW3Z35/VHdTO8wszs7O9fuvu/ns5/t6a7uetPTXa/eq1evfMFgEEVRFEVpKv50C6AoiqK0TFSBKIqiKHGhCkRRFEWJC1UgiqIoSlyoAlEURVHiQhWIoiiKEhfZ6RagtSEijwKXOh/nGWMOCDv+BPDjGC83zRgzXEQuBh6P8ZzjjTHvNSDfWOBPEQ4FgGpgDTAFGGuM+S7GOhOGiOQCVwH3G2O2NeM65cBFxpi/J0q2ZNDA87ADWAvMAP5mjPky7LzhwFTn4z+NMVc1Q4bjgCJjzBtNOMeN/59jjDnE2Xcxoef0amPM3fHK1EjdxcBVxphbPftSUrdSH7VAEoiIFABne3YNFJHD0iVPE/EDBUBfrAL8RES6plIAETkJmAvcQTM6NyJyJfAtcGGCREsHBUAf4ALgUxE5u+HiTUdEuorIf4H3gH0Sff1kICI/BL4Brk23LIpaIIlmDFAatu8S4BPP53uACZ7PJwC/crafBf7nObY+Qh3hZcKZG4ugDvdirQ2ALKAdcAUwGNgLuB5rDaSKPwB7J+A6dyfgGungJuzvlw2UAaOcv1zg3yLylTFmgVN2LvZ5A1gUZ32nAD+M81y37i1xnh8vf8U+m+H1TiEk01cplagNowoksbiuiFpgE1AO/EhEfmOM2QFgjJkFzHJPEJH2nvMXGGMmNFJHLGViZVb4tURkKrDY+Tg0QfUosTEjzP34iIjcAVyHtUhuxFokGGPWU78jklIS+AwmBGPMcmB5uuVoa6gCSRAi0h040fk4CfgauBrbqx8DPJMm0ZrKZs/2xvCDItIb6z44DeiG7QlOA/7qKMfw8p2A3wBnAr2B7cBHwN+NMVM85cJz6mwSkWXGmD7O8eOB32KVWhlQhe2FP2iM+Y9TZjihcQGAg53r/tsYc3HY+M9RwDjnemuBY4wxyxyX4++BQ4HOTj1fe+tx6vJeqz/wA+D/sL3j+cCfE9TI/gVroRYAZ4nIJcaYmmhjICKSg33uzsNac7nAOmA6dlxrgVPuCeqPvfxDRP6BM4YmIkuxv9dkYDzwN6AD8Iox5oeRxkDCEZEfY63YvsAS4H7gPmNM0Dnex9mPc93RnnPvBq50Ph4PLPWUBWjnyBBpnLDeGIiIZAM/d76vYN2184DHgMeMMQFPWfe+1AF5WOX9U6A71i36N2PMk5G+b1tEx0ASx4WE7uczwNOeYz9JvThNQ0T8zpjH3Z7d48PKHAXMBn6JbRTysI3s2cDHInJ+WHnBWlu/x764+dhGaCQwWUR+F6NsJwPvAqcCnbDutlLgSOBJEflVA6dH41msEskF6hzlMRQ7HjAa6OEcK/PUc2nkS/EAdtymL/Y7DgZeFpFm/+7GmErgC+djPnBwI6eMd2Q5GCjEdhK7AecC74tIvyaKMBB40rlGPlaZxsJPgSeAfbHPyb5Y9+1dTay/WTjjkpOA+7CdhVKgGDgceBiY4CiYSPwXuA3oh/0OB2BdiacnW+6WgiqQxHGR838HMMEY8zlgnH0nOD33RPAnEQlG+ZvQxGs97p6L7XGtxn6PWuAPxpiX3YIiUoIdeylzdj2D7and5ZTPBsaLyH5OeT/2BezplH8LOx50M9YKAbhdREY422OwvUKXC4DLnO3fYJVGFbZHew62Z+j6wW8RkULqjwuA7bWOwTZc4XQH/gz8DOtXBxiL7enXOnKe6/x3e6jRBuVPBl5w7sejnv33iEiHKOc0hVWe7W7RCjkdgAucj59hgyF+CDzl7CvH3j+w9+Rez+njsfcqfAytK7AC2wn6Cw2Pv3kZiLVeLsZaerXO/qtEZFCM1/CyzpGvwvm83fl8YyPn3Q4c52wvAH6Btehca2YU8McI52VhreY7sdbcR55j/9dE2Vst6sJKAI7bYz/n40RjzFZn+2ngFsCHfZFuTr10ceEHhohIZ2PMOmffD7G9coBHjDFu446ILAL+he2xX41t+EcAbkPxDjDS47r4BHgde19+B0w2xkwQkas8MrxujNnsbBc7/zcBrxpjvnau8xHWIpkH7DTGbMf2KN1rbGnAjfSSMSa84fgD8CJQ63VTiMhIbO81WuP9tjHGjZJ60glFvggowloz46OcFys7PNslDZQr9mx/CbzgWDDPishkYBmOgjDGzBKRgzzlv2rgXv3KGDOxiTIvAL5njKkDEJEdhFx+52At2Zjx/LZ3O7tqGnMRikg7Qp2QCuAI95kSkRew0VwlwNUi8hdjTHXYJe40xlzvlJ+J7ZCAdVkqqAWSKLy+ZO9Yh9eNdbGI+BJQ17PYnlekv9ubeK17PeeejX3ZXsM+F2cC74pIllP2eM95D4Vd51GsBQNwbITyD7vKA8CZb+DOMTkmBjndXm8PYL6ILHJ81f2wPvBvvH7sGJkevsMY84UxZjwwTUTOF5F/isgcQsEEOVGu9VzY5+c92wc2Ua5IFHi2t0crZIxZiLU8wFoM60VkpjMQvxx43xhTEe38BtjjXsXAK67ycHjRs71feOEkcTihe/esp0OCMWYN8IrzsZhQZ8fLu57yy4Cdzsf8hEvaQlELpJk4vU1vKOTLnh6wlz7YRnVKpINNIKlRWNjIn6nAcOAg7HjFa0BHT5mV3hOMMbtEZD3QBTvGQUPlPft6APkiUuj0MCNijLnPcbP9DusS6+f8/RjYKSK3G2NuaeiLRmCPEGnH/fYIdmwErNvlC6z7pDPWYopEeLCBt5Eupvl45+OsaaTsSKzb5RxsQ3eE83cdsEhELjDGfNTA+eHUGmPiCdUNvyfrPNtFEcqH39tEtE2xPIMukVyNm8M+V2PvaSI6gq0CtUCaz/eJ/PBFIuMH0x0+9my7vUVvA9DdW9hRouXOx4rGyoft296Q8nAxxvwLGxV0KLYxfBXbG88HbhaR0xq7Rhhet5AbqTMRqzzWYH/XdsaYoVhXR0OEu7bKPdubmyhXPUQkj9DA+S4acf0YYyqMMT92ZPg+8HdCVkl/bAcnK9r5EdjReJGIhN8Tb2PuKiRv5F24ddeQqy5WYn0Gob7Sd6mLsE/xoBZI8/G6r+4n8oN4PTaK40wRaRdnjy6VeN0uG5z/HxIaoP0Z8LmnzKWEOiPTPOXxlN/twhCRUwmNp0zzlPO6oXxO2SLsuMpAYKMx5gqn7r+JyDlYlx7YXvbrznbQOb+hnmK4y+sgrFUDdkzjNaf+QmwEWUNcICIPeNx0oz3H5jRybmP8lJAV87wxZme0giJyKNYa3h/4l/Md3O/xBjaKrSs2WmwhEe53BJrqGnQZIyK/d+c/EfmeVHn2hSucaAPtrjyxWAGfY5VuLnCOiNzovnsi0gU4wym3hVCkm9IEVIE0AxHpjJ3NCzaC6VeRfPEiciB2TKEA+4KHjyE0hX1FZHQDxxd4Zis3xmAR2ez5XISNKBrpfN4FvOls/w8bBNAJ+LnTsL+DbXivcsrsJBSm+Q42bn5v4Hsi8jp2rKA31oIA29B7x228lsivRGSdMeZBsWk8DoLdPfJJ2HDMX3jKe+cIbHe+S38RuQTYFIPbr9KzfY6IzMW6YX5ByKLIi3LuEcBrIvIccBghS3MLIT97LBztTCz1Yb/fsdjgC7D39rZGzq8DrnG2DxWRv2EHzvcjFIm0k5AbzHu/R4nIFmC6MebbJsgcjV7YMbSHsGlS3JDtOkJjWhuxFlp77Jyd32I7FJcRfezIlblERK7ABko8FamgMWaziDyNjf7rDHwoIvdgI6x+Q8jK+ZsxpiaeL9nWUQXSPM4ndA+fb2Ag9ymsAgHbuDRHgZzr/EXjZmw4aiz8ilAalUj80RizCsAYs1FEfoS1JEqx1sgFnrK7sMkLFzrla0XkB9jw3b2wSmmkp3wQuMYY4x2gnY2doOh+j2qnAfoJduyoFGvthM/H+IRQqKp7naOxSmQ8djB0QgPfE2w6kPexjXYBduKcSw3WxdJJRPIjWAFzHLm9brQgcLknIi8W/hxl/y7gJ8aY+Q2dbIyZLSJ/wUaTlWPDZ8P5vSdJpdcddrzzdz5W8TeXl7DP/FFh+29xFZQxJuAEQ1zlHHPlDWJdlN+PcN3ZWGvUh53b8S31f/twrsRaZIdjFekDYcdfIBTGrTQRHQNpHl73VXgkjpc3sCGoAIeJyMDkidQsarE98ZnA+caYO7wHjTGTsJbAvdh0J9VYP/NzwGHGmOfDyn/plL8dG9a5E9vrfB074/kfYfX/3bnWemAb8CmQ78ypGQQ8iB2P2OH8fYUNDT0+LATzCqwy2IZ1wTU2hoHjfjoLG1G2ktBM95uxLjiwnYVTIpw+FtujXebck8+AUcaY/zZWbwPsxFpV44EhsV7LGHMj1jXzLvZ71GLvwTuOTHd7yi4Cfo1thKuxbq1NJIZ/YDsYc51rzwd+GiHY4TpsA74K+5vOwFrBzxOZP2Ct4i1Y66UxpboV25n4JXZsrxL7236E7YicExYtpjQBXzAYnkFCUZTGkPqpTMZkWm4oRUkFaoEoiqIocaEKRFEURYkLVSCKoihKXOgYiKIoihIXaoEoiqIocaEKRFEURYmLpEwkDAtxDGcXNt3HHGCcMWZalHIpxbMC2xZjTPv0StMwzmzsn2Inah2AXaNjCzbm/kVsuvXw1NTx1tUXm5b7wURcT2kcqb9a4CBjzBcNlB1L9HctnL7GmKXNkU3ZE9lzNc3Lw98XEXkWm+DSZZoxZngSZBlLM8PLm9IWpsMCycUmMRsJTBWRUWmQocUiIj2xE9XuA07ApmjIwaYYGY6d5DdXRPZtZj05InIzdqJWpMlziqJE5vgI+46LsK/Fk4pUJvcSSmHuw6aJ+Al2wSEfdi3mid71IpTIiEgxNjXI/s6uZdi0KIuwyeguxmZuHYBdMnaQZ0GoptKdyCu1KZmL912LRLzPgtI0hns/OMsEdEmPKMklFQpkjzUnRORF7DKZ5dgU0wNITP6d1s71hJTHx8BJ3lxLInIfNhXIGKxCuQObSE5pG0Ra30VJHeuxnoDOInKAMcZdHti1SAJYV3NZpJNbImlJpmiMqXaWQXWznLZ3j4lIPja52nlYP1wusBab2+dPbnI/p6xrtTyGXTr2b9g8OtnAVOC3xhh3XXL3nP7YrKYnYV0/72LzGEVFRI7CphQ/EruuwUrs2hG3G2NWe8pdDDzufDwCaw1cjV1MagE2kd2bInImtne/L3Zlvvu8OYqiyOCj/noiPwtP1OckMPwZNm13PnC+iFztZCX1yna1tz4R+cKRFWOML6wswBnOvb7ZGDPWOScXm0fpQmzG3e3YNeDvi5S3yVmv45fAEGxSxCXY8Zq/edPbh/lwe2KXhr0M24ObDVxpjPlURC7DZp7tjc3hdLsxxrsCpLuO+41Y33M37Njbm873+I4YEJFzHbkFaIfNFfUJcJt3YSYReQ/rpliEzf81FvsMd8Tm7BrrrMTovXaxU+5cbMMzh9Ca5UnBySA9z6kP4Pue1PV3YZ9XgGeMMeeLSB9CmY5vxmZCvg27LssmbCLDsZ607W49hdj3+IfYNPm1wCzgfmPMC2Fli4HfY8f0+mCz5a4m9M6v9JR9ggjjQ04WYzeP1+7xhbDn6Shswsah2DblGGPMMie1+83Y9dE7YfNyvQzcaowJXxirIWY61/BhlYarQIY7/+dg1w6KqEBifUecsjnYDMcXY5dG+AbbBkZFRE7A5hI7zJFxLjbt/39i/4r1SUsUlpMK3LvGwmLP9lPY5HsHYm9iPraR+Cl2qdFIKbX7YhPvnYNVRsXYH3KKt7yjPD52ypU55cYAH2Abh0iy/hKbmO8sbFbZXKe+XwGzpf660l7+gk3+J9g04AcDr4jIP7APxcHO/v5YN94volzHZaBTP9iU7V9FKmSM2YBVnmAVZHg21GbjPLxvYBX2QVi3ZEesgn1GRP4UVv4OrMI9BdtpyMMqzz8An4pIDyLzb+x97E1odb13ReRhrOtuH+daA4GnxK5d7tZZgk3Mdx22UXLH3n4KfCYiA2L4nr/Eph4/2pE7F6vIRgHvicj+EU7LBd4GfuvUl49tsF4VEXdhKFcBT8YqwR5OucOxSQ+PaEy2eHFcmr/07PqniOQ7a4n82tm3mshZmo/FPlvHYH/zbtj7+7p3kSpHSX2M/e0OxGZFbodtVJ8Xkfs9ZX3Y5Jo3YJ+JfOxz2wub7HCqiMS6YFtjPIt9H3KBOkd59MC2HT93vk8u9nm5GpgpIh2jXCsSmwgpDe84iDv+8X60E+N4R57DZm/u75Q9EJuA8uwo178Eq5BPwLZ7Rdjn7UkRuTO2r7cnqVAgg0VktPP3AxH5OfYlcbXwa06jh4gMxTbUYHt5F2I17KfOvgFEfrlOwPaAr8C+HO6aAd2wL7vLXYRWRvsCm2X199gb2j78os5L9Q/sfarGZg29GJuqGmxj8pLTGESS6Qlsb8ld4S8H2yv7BNuzftJT/uII1/DSy7P9dSNlvVZX30bKRmIK9oVy+QSraN11HK7GjmGB7c1ehm1w1jr7/igi/QCctUvc9T8qsS/ETwgpub2pv468l+OxGXp/Qiijbjvs7/YWNtvr657yF3u2/4Kzhgj2dzgfuAnYiv3dGowqcyzhm52Pa7HP1nkeWfOoH1Xj0hMYjLUwL/TInUUoqy/O9Q5ztjcA12Lv+SKsYoyHx0UkGOXvbreQMeZZQtlu+zqyPurICNa6jdTzPh67tvrlWKt9g2f/Tz3lHsZGB4LNevsz7DPgLiN8uYi41vThWMUEtoH7Mfa3cq21vYH/i/H7N0Z3bKP7M0Ip3O/H/mYB7Lt+Htb1W4vt/N2x52UaxFUSx4mI38m83TnsWD2a+o6IyBhCC3Rtx1pYl2Lv9R6dGhHphv2efuy46a+x7c9kp8g1InJME77jblLhwmpozYnZ1F8UaCF2TYVDgH+7pquIVBBqKMJXLnM51V1ISUQ6EVoTo7+zr5DQehQbgeGe1ckWYE3WcK4mdI9+Yoxxf8R/i1086Gzn+j9gz0ZwijHmEuf6C7FWDthU0qcYYzaJyH8JrV3dM8r3cvGavVVRS+15PKJl1RDGmOUi8o5n1+ow37rbWOwETjTGrAEQkZXYnuRsj7y/9Zw30hjzgVP231iT/3DgGBE50hgzM0yUx40x1zrlqwHXRbUcm5q8VkSmYcfTwLmHYpenvdjZN9n9HZxj27ANxQgR6WeM8Vq/XgJYpXkI8IUx5n3n/FewjQxEfxavccM4RWQF8J6zv7+njLenOMZdF0VEXsNa5PlRrp0ofoHtGXfGdqJcHjfGvB75FHYBJxhjlsFu16c7aH8+8JCI7ENopb/5wLHuYk0i8jbWjZXl1Dme+mvGf4xdV2eHiLyAbeS+xnZSEsFLxpjdgSEi0hU43fn4uDHGdWX/1+lAXAmcJyK/NjEsu+wwHds56ID1MhwVdiwSTX1HvM/O/7kuKKc9WYxdddLLhYSepwuMMTOc8v/DuuPLse90NPmikq4Fpb7F+qZf9S7OY4zZBLwhIpOBw0XkfKxbxGsOhq+dDPBd2Cp83m33xgmh7zspzKf4KrbHEX4/3Hq3AuF+/YcI/ZDHsqcCmeHZXuvZnuN8T3fMYj0hF0ZDeJVCcdRSFu960pVRS8WB4xra2/n4has8AIwxL+NRxI770LUY57ovhlM2ICKPYF8OsPcwXIFEu4czjDG1zvZqz37vb+3egxGyZ5y+y+HUd5/uxhizC3hfRKYDB4nI5U75EZ5ikZ5FsD1pl0jPIli3G8A676JaxpjVIjITa8E2lYaisBZ5Pxhj1jvf6UXP7hWExkAi8bGrPJxrTBWRDVir3l1BcLin/HjjWenPGPOliHyIdQkOcBrwD5x6e2LbhGudMtOBN73jTAkgvIEcQmhp3EtFJHyhMrCuuoOpv0RzQ3itjOOx7RfA18aYCpH6qyPH+Y541xN60VN+h9ili71jpWDHq1ymh8vgMKyhLxWNVCiQS7CN72FYn3ZfbAN0HHY1sHqIyI1Yc97tOa/F9mTcLxhpLeTNYZ+9k+jc8kWefRs82+4PtZ49Nbfr7loVIcx4pWc7ko/W23B7F6wJX7DHPdbYGs9LPdsHRCsU4fiSCMfD62rKc+C1aNZHLWXp4KlrZYTjzb6Hxpg6zwvh1hVrlEs0CwIAEfk+8E+sTxzsAlUfE1rPPdpvttmzHelZhNDzWO9ZdIg33LapUVhvYAML3GCW2eGDtWFEcmtVYN8Tt1PjHTNo9Dc3xqwRkeFY9/JpWCXrro74R8fKOdcYE2lRMO/9jOUZDn9eE/KceHE6AAux7vYTCDX+0cY/4nlH3GdnRwTLKNKzE8v3jPk7ekmJBeLMip7uDHJ+DhRiTegNeOYaOAM97rKej2Bnqi8UkRHY6I9oxLKimPfFKPceEBE/9R98l3XYxqKbiPjClEh3z3ZFE2SqjbK/Mb7ENpxl2N7bYcaYT8ILiUh3Qj7lGkKuM6/s4T3nEmLHex/r3TNncF2Ab5wefIVTr4/698slWffQG502GTvpMhJRx5LEaqUXsPfqA+xg92cx1h/r89iBsGfRoXOEfcnglrD6vy8iZxljXoxSPlIj456/2fnvbcBi+s0dN+JoZ7D8ROwg/QlYf/4h2MAad7wo2nMcyzO8I+yz9zn5L9FXFf08hmt7mY5VIKcQGleKpkDieUfcd7BARIpNaIliiPzseL/nz4jc8YtrHl5Ko7AcN9O1nl03iMjhns9jPNu3GWd9bexyps2uHuuzBzhRREo9x84ksjvCNVtL2HMd8ss820lPx2LseuvetdQfCY9OcczhJwh1DP7rusuo7wLr5jmnM5EfXO/67rt7ek7osPu7DBYbAulyCjZktUpEfuO4mdxG9wAR2R0A4Sht78BrIu/hAkKNRQ/gdWPMBKd33hVr/XakYQvqNELPxKPGmI+NXfo0Ec8i2CAOsOus7x7AdKJtEh45F44TsOL6/NcQakDubyDyaIjYSXHuNY4i1ImY4/z3unoudsaj3PIHEHLXuC6dE0XknyIyCehtjHnOGPMrY8xAQsvVHuoJVIn4HGMDFxojEPb5C892F/cZcZ4TV9ZSGre0w3GVRVaEffWI8x3xyj3GU76I0JiOl9me7WrPd5yGDVrajz29ODGRjjGQB7FCj8De4EdEZLBzI70ui0dFZDzW/3iVZ3+kMN5GMcbscgaNLsZGXE0WkXuxjecNUU67HztA7gOeEJEDsYro+4SixRYQispKNrc59e6NjTD6SkQewEb69MD2LtwUJmuwceIuXh/4pSIyBdszuZX6D7qL1zQeIiIXYl15k7Hzbm7H/hbvi50/kIWNHsHZfsvZvg/rugQ7vjXOke0CQm7JSZGsqXgxdp7Rf7G+YAHeFjt/oBs2uCIP25traJ1x77P4R7Hhph2o3wGK61l0eJLQOMfzIvIXrLvr2mZcd7CIbG7g+CwnQCIXO8/H/d0vwAYG/ATbg70P+FGE8/3YMOq/Yt8J73vzJIAxZr6ITMW6oA7APh+PYS3n6zx1/sX5X0oofPhF51mqwPrtXd/kSseihfrP8S0e1/M9DXzviDhhvFOwv8MJzsD9y47cv3O+4zygqfMkwpXFkkbmHTX1HXmS0DjHv5woq3XY6LhIFshT2PGlHGwHoS+2E/hzQt6K64ijE5fyeSCOG+hSQmbVgVj3ANiH2u0JjcBG3VxHfV9ntDkDsXA9NowN7AP6b2yjvJ1QqLBX1vew0SIB7Et9g3OOqzxWAWd6BnSTitP7PwnrzgLbIP4ZG9/+d0LKYwme6CiHOdgIGLD+6lewA677EaF3ZIxZTyi6qTv2oXWtsLsIPWz7YDsF/yLUI/yTMWa+c50nsYoYrOK+DRt94zaeX2OjRBLN9YSyGxyP/d1cpVcLXNZIZM2rhHqefbEy34ltrNyB4eY8i08SClXtgm0AH8JGar0V7aRG+BW2AYz2597zmwgNxD7ldAq8YbY/FBsqGs5CrCK415HXHTN8k/qN7EWEwpePwIYI/42Qu+tuE5r0+TK2gQN7n+/Fhotfi1U2Aep3IJ/HjkWBbTumO/tWEnkMoTH+j5Br6CxHluuxbc424KcRxj8bxHHJeWWJOv/DKd+kd8TYBLRuGHoJNiR5PHa+0R4RdI48bqRXMTY8/WlCymM60d28DZKWiYROJIe3J/dHEenrPMhnYn2O27EN9BtYk959IONOvmiMWYuNivgP1mSrxEYxHI4NDY10zh3YG/0cNuJnF7aB/id2JmxjczISinPvhmBN23exD38tdoBzOjaKZqAxZl7YeQFsGPOzWB/qFqyPfyjR08j8GGte78D2iJY619oFfA+rXOdhe84VOBMujTF/9l7EGHMF1rR+3SlXjbXc/gwcHqboEoIxpgL7u96Jbfiqsb20N7Ah3BMaOX8d9gV+Czv2tAkbZ38+oVn6w5yQ8XjkC2Kf9T9jsxHsxI61nEiERiBRiMghhGa7b8JxYxk7F8v7Tj4Q7iLF/tbHOXJWYxvJ27BhyLsbWae3PQjb4ZqDfZcrsfMaxhhjrvaUDWK9Ahdho+7WYp/ntcAEbBjwi2HXPhH7rO/AthF3YTsJuyM6Y8UY860j60PYDtMu53s9BxzRjCgwb8RXgwrEkaOp78gvsL+d+2x/gY0KjTi/yRjzT+z7Pwn7u+90rv8nbOhw+PhQTOiKhIqiREXqpzJ5xRgzOn3SKJmGLiilKIqixIUqEEVRFCUuVIEoiqIocaFjIIqiKEpcpCsXVkxUVGzdQ7uVlRWyaVOsec3Sh8qZWFTOxKJyJpZMk7O8vKSx1EgJocW5sLKzI815yzxUzsSiciYWlTOxtBQ5E02LUyCKoihKZqAKRFEURYkLVSCKoihKXKgCURRFUeJCFYiiKIoSF6pAFEVRlLhQBaIoiqLEhSoQRVEUJS5UgSiKoiSBp556gjPO+B7V1dUAPPbYQ0yY8AIAw4cP45e/vIxf/ernXH75pdxxx63U1tZy7rmj2bRpIwDr16/n2GMPY+rUSbuvec45Z1BZuWXPytJERqcyySQKx91W7/P266KtgqsoigLvvvsWI0aczOTJ7zByZP118EpL23HffQ/v/vzHP/6ejz76gEMPPYw5c2YzfPgIPvpoBsOHj+DDDz/g+ONPZNWqlZSVdaC0tF2qv0pUVIEoitIqGTvzRl5bNCGh1xzVfzRjj7y10XKzZn1Gt249GD36LG655Y97KBAvtbW17NixnYKCQoYOPZw5c77YrTh++tPL+cMffkswGGT27M85/PAjEvl1mk3SFIiIZAGPAALUAZdg1xl+Arvu+VzgCmepVUVRlFbDxImvMGrUaHr16kNOTg7z5s2td7yycgu//OVl+Hw+fD4fw4YdyZAhQ6msrOTpp/9NbW0tq1evom/ffvTrNwBjFjB79ueMGXN2mr5RZJJpgYwCMMYcJSLDsesW+4AbjTHviciDwBnAy0mUQVGUNsrYI2+NyVpINJWVlXz44Qds2rSRF154lqqqbbz00rN069Zjd5lwF1ZofylZWdl89NFMDjzwYACGDTuSr76aw+LFi9hvv/1T9j1iIWkKxBgzQUQmOh97A2uB04Bpzr43gZNpQIGUlRVGzHJZXl6SWGFjoSiv/scYZEiLnHGgciYWlTOxtDQ533prAmef/QN+97vfAbBjxw5GjBjB975XRnFxd8rLS/D7fVG/19FHH8nzzz/Nz3/+c8rLSxg58iSuvvpq9tlnAF26ZM74ByR5DMQYUysi/wbGAD8ATjfGuGt8bAUavBuR8uuXl5dQUbE10aI2SmFVdb3P2xuRIV1yNhWVM7GonImlJcr53/8+y0033VJP7mOOOZ7nnnuOq666loqKrQQCwajfa+DAQYwf/zj9++9PRcVW/P5CtmzZyimnjIr5XqRK6aZkRUIR6Qp8DJQaY8qcfWcAJxljfhntvEgLSqVNgTQxCqslPviZjMqZWFTOxJJpcrb4BaVE5EIR+b3zcTsQAD5zxkMATgWmJ6t+RVEUJbkk04X1EvC4iLwP5ABXAV8Dj4hIrrP9QhLrVxRFUZJIMgfRq4BzIhw6Lll1KoqiKKlDU5koiqIocaEKRFEURYkLVSCKoihKXGguLEVRlASyevUq/vSnG+jduw/ffLOAkpJSfD4fdXV1XHvt75k2bQq5ubmcf/6PAbjmml+TleVn3Li7AZu1t6SklHPO+VEav0VsqAWiKIqSJC6//Nfcd9/D3HvvQ1x44SU8+uiDDB06jDlzvgBg586dVFVtY+3aNVRX7wRg9uzPGTYss5ImRkMtEEVRWiVjx+bx2muJbeJGjapl7NjqxgtGYOvWSgoKCthvv/1ZsmQxwWCQzz77hEGDhlBVtY1Zsz5j8OChbNq0kV69+iRU7mShCkRRFCVJPPDAPTz11BNkZWXRqVMnLr/8SrKysth7731YvHgRH300k5NO+h5VVVV89NFM8vLyGTRoSLrFjhlVIIqitErGjq2O21pIFJdf/muGDTtyj/123Y/ZzJv3JVdddS21tbU8+eR4SkvbZdyaHw2hYyCKoigpZujQw5k06W169OhFdnY2+fn5lJSUMGvWZwwZMjTd4sWMKhBFUZQU06NHTyoqKjjiiKN27xs69HD8fj+FhUVplKxpqAtLURQlgey1VzcefviJRss9//wr9T6fc855nHPOeUmSKjmoBaIoiqLEhSoQRVEUJS5UgSiKoihxoQpEURRFiQtVIIqiKEpcqAJRFEVR4kIViKIoSoKYNeszTjllOGvXrtm974EH7uWNN16LWP6NN15jxoxpqRIv4agCURRFSSDZ2TncdtstBIPBRsuOHDmKo49uuat860RCRVFaJUVjbyTvtQkJvWb1qNFUjb21wTJDhhxKIBDkpZee46yzzt29/8EH72PBgvls376dPn36csMNf+Kxxx6iY8eOrFixnAED9uHUU09nw4b1/Pa3VzF+/FM8+OB9zJkzi0AgyLnnns8JJ5yY0O/TXNQCURRFSTDXXns9zz77DCtWLAegqqqKkpIS7r77fh58cDzz5n1FRcW63eVHjRrDm29OBODtt9/gtNNG8eGHH7B69UoeeGA899zzIE8+OZ6tW7em5ftEQy0QRVFaJVVjb23UWkgW7dq159e/vobbbhvLgQceTF5eHmvXruFPf7qBwsJCduzYQW1t7e7yffr0pa6ujjVrVjN58rvcfff9vPrqSxizgF/+8jIAamtrWbNmNSUlJWn5TpFQC0RRFCUJHH30sfTs2Zs33phIdXU169at5eabb+Oyy66gunrnHmMkp59+Bvfffw99+vSlpKSE3r37MGjQodx338Pcc8+DnHDCiXTv3j1N3yYyqkAURVGSxJVXXkNeXh47d+5k1aqVXHbZxdx00+/o1q0769dX1Ct7/PEn8sknHzJq1GgAjjrqWAoLC/jFL37KpZdegM/ny7hMvb5YIgXSRUXF1j2EKy8voaIi9X7AwnG31fu8/bobGiyfLjmbisqZWFTOxKJyxkd5eYkvFfUkbQxERHKA8UAfIA+4FfgOeA341in2gDHm2WTJoCiKoiSPZA6iXwBsMMZcKCIdgdnALcBdxpi/J7FeRVEUJQUkU4E8D7zg+VwLDAFERM7AWiFXGWMyx+5TFEVRYibpYyAiUgK8CjyCdWV9aYz5XET+AJQZY66Ndm5tbV0wOzsrqfLFzNixDX9WFEXJHFr2GAiAiPQEXgbuN8Y8IyLtjTGbncMvA/c2dP6mTdv32Je2QfSq6nqftzciQ6YNqkVD5UwsKmdiUTnjo7w8NXNFkhbGKyJdgHeA3xljxju73xaRw5ztEcDnyapfURRFSS7JtEBuAMqAm0TkJmffb4C7RWQXsAa4LIn1K4qiKEkkaQrEGHMlcGWEQ0cmq05FURQldehMdEVRFCUuVIEoiqIocaEKRFEURYkLVSCKoihKXKgCURRFUeJCF5RqAjkz3id77ldUf390ukVRFEVJO6pAYmX7dnI/mglAzocfpFkYRVGU9KMurBjJWrUytL1iOdTUpFEaRVGU9KMKJEb8jgIJtC/DV1tL9tfz0iyRoihKelEFEiP+dWsBqB14AABZixamUxxFUZS0owokRvybNxMsLCTQdS9AFYiiKIoqkFiorcVXuYVA+zICZR0AyFr4bSMnKYqitG5UgcSAf9VKfIEAwXbtCZaWEvT7yVq2JN1iKYqipBVVIDGQtWwpAIF27cDvJ1hUjH/NmvQKpSiKkmZUgcRA1lJrbQTbl9n/xcX4166BQCCdYimKoqQVVSAx4F/5HQCB0lLAKhBfbS2+DRvSKZaiKEpaUQUSA/6KdQAEi4rt/2L7379mddpkUhRFSTeqQGLAnQMSLCqy/4vtgvVZa1WBKIrSdlEFEgP+dWsJZmVBXh7gtUB0IF1RlLaLKpAY8K9bZ60Pnw+AgGOB+FevSqdYiqIoaUUVSGMEg9YCccY/AIKFhQD4N6xPl1SKoihpRxVII/g2b8JXU7N7/AOAAqtAfBs1CktRlLaLKpBG8K+rH4EFECwosMc2bEyLTIqiKJmAKpBGCI/AAiAri0BJKX61QBRFacOoAmkE/1obaVVPgQDBDh3UhaUoSptGFUgjhFxY9RVIoGNHa4EEg+kQS1EUJe0kbU10EckBxgN9gDzgVmA+8AQQBOYCVxhjMjqhVMiFVVxvf6BDR3zV1VBVBcXFkU5VFEVp1STTArkA2GCMOQY4FbgPuAu40dnnA85IYv0JYbcCKQx3YXW0x9WNpShKGyWZCuR54CbP51pgCDDN+fwmcGIS608IUV1YqkAURWnjJM2FZYzZBiAiJcALwI3AncYYd9BgK9CuoWuUlRWSnZ21x/7y8pLECtsQGyugfXuK2tVXIHSwmXnL6nZAFHlSKmczUDkTi8qZWFTOzCVpCgRARHoCLwP3G2OeEZFxnsMlwOaGzt+0afse+8rLS6io2JpIMRuk4+rVBMo7U11VXW9/oF0nSoDKxSuojiBPquWMF5UzsaiciUXljI9UKbOkubBEpAvwDvA7Y8x4Z/dsERnubJ8KTE9W/Qmhpgb/hg0EOnfZ45C6sBRFaesk0wK5ASgDbhIRdyzkSuAeEckFvsa6tjIW//oKAAKdO+9xLNihA6DpTBRFabskcwzkSqzCCOe4ZNWZaNwIrIgWSKkdvvFv3pxKkRRFUTIGnUjYALsVSPmeCiTYvj0Avi1bUimSoihKxpDUQfSWjhvCG+jcmazly+odC7RrD4Bvy+YUS6UoSiyM++S2qMeuO+yGFErSelELpAHcPFiRXFgUFhLMyVEXlqIobRZVIA2w24XVpeueB30+gu3aqwWiKEqbRRVIA/jXRh9ELxx3GwTqyFq10m4riqK0MWJSICLyWxGJ0A1v3fjXrSWYlUWwY8eIx4N5+bBzp2bkVRSlTRLrIHoh8J6ILMJm051gjKlJmlQZgn/dWgLlncEfWc8G8/PxBQJQ0+pvhaIoyh7EZIEYY242xuwL3A4cD8wRkftE5JBkCpdWgkGrQCINoLvk5wPgq96ZIqEURVEyh5jHQESkCOgL9AMCwEbgnyJye5JkSyu+bVvx7dhBoEt0BRLMswqEnapAFEVpe8TkwhKRp4ARwBvArcaYGc7+PGA18PukSZgmGhpAdwm6FogqEEVR2iCxjoFMAX5ujKlyd4hIrjGmWkT2T45o6SUUwtuQBZIHqAtLUZS2SawurJ+FKQ8/8DmAMWZNMgRLNw2lMdmNWiCKorRhGrRARGQKMNzZ9q5dXgu8mjyx0s/uWeiRJhE6uC4sqqujllEURWmtNKhAjDEnAIjIP53sum2G3WMg5XumcndxB9HVAlEUpS3SmAVyujFmIjBLRC4KP26MeTJpkqUZ/3fLAQj06BG9kLqwFEVpwzQ2iD4UmIjjxgojCLRaBZK1fBnBnBwCXfeKWiao80AURWnDNObC+pPz/xJ3n4iUAj2NMfOSLFtayVq+jLoePSErK2qZ3WMgaoEoitIGiTUX1qUi8oSIlAPzgRdEpPUm1N+2Df/69QR692m4XE4uQZ9PXViKorRJYg3j/QV2suCPgFeAA4EzkyVUuslaYcc/6nr1abigzwf5+apAFEVpk8ScysQYsxoYCbxujKkFCpImVZpxVx+s69W70bLBvHzQMRBFUdogsSqQeSIyEZsHa5KIPAt8mjyx0kvW8qUABHrHoEDUAlEUpY0SqwL5CTAOGGaM2QU8BVyaNKnSjH/ZUiBGCyQ/H19dHezYkWSpFEVRMotYc2EVY8c9jhMRn7NvEHBLUqRKMyEXVp/GCzuTCf1bNhMoaLVePUVRlD2I1QJ5HrsOSBbg8/y1SrKWLSNQVEywQ4dGy+6eC7J5c5KlUhRFySxitUC6GmNOSqokmUIwiH/5MgK9etsoq8aKqwJRFKWNEqsFMltEDkqqJBmCb+NG/FXbqGtsDoiLx4WlKIrSlojVAjkAq0TWAjux7qugMaZfQyeJyOHAHcaY4SIyGHgN+NY5/IAx5tk45U44heNuA8C/ejUAdTFEYAEEC1wLZFNyBFMURclQYlUgY5p6YRG5DrgQcNcRGQzcZYz5e1OvlUp8jiURiCECC0IZef2qQBRFaWPE5MIyxiwDjgIuAyqA45x9DbGI+rPVhwCnicj7IvKYiJTEI3CycV1RMUVgoWMgiqK0XWJdE/2vQA+sErgDuEREDjbGXBPtHGPMiyLSx7PrE+BRY8znIvIH4E/AtQ3VW1ZWSHb2nskMy8uToHuK7PK0VG0FoN0h+4O3Hvd4OGWl9nB1FUVhciVFziSgciYWlTOxxCtnUbR3thnXbIiWcj8TSawurO9hXVCzjDGVInIS8CUQVYFE4GVjzGZ3G7i3sRM2bdq+x77y8hIqKrY2odrYKKyyqwrmb9hIFlBR1BE89bjHw/EFsygEdq5ex1ZP+WTJmWhUzsSiciaW5shZFeWdBRL+3TPtfqZKmcUahRUI+5wXYV9jvC0ihznbI3DWVM80fFu2ECwohOLimMqHXFg6BqIoStsiVgvkOeBZoExErgIuAp5pYl2XA/eJyC5gDXY8JbMIBPBVbiHQuUvs5+TkEPT78esYiKIobYxYFcjrwCpsMsVjgJuMMa83dpIxZikwzNmeBRwZn5ipwbdtK75AgGD79k04yUnprhaIoihtjMbWRO8MvAAMxM7fqAVOAApEZIYxZkvyRUwdvi326wTatW/SecH8fJ1IqChKm6OxMZDbgRnYVCbDjDHDgM7AHOCfyRYu1bhKINiuXZPOC+YX4Nu0CYLBJEilKIqSmTTmwjrSGLOfd4cxpsZZzvaLpEmVJuK2QPJsSndf1TaCxW0vlE9RlLZJYxZIxJWSjDFBmh6FlfH44rRAKNDJhIqitD0aUyAN+WRanb/Gv2ULQZ+PYElpk85z05n4NulAuqIobYfGXFgDRWRxhP0+YK8kyJNWfJVbrAsqa8/Z7w3hzgXxb9lMXTIEUxRFyUAaUyD7pESKTKCuDt+2bQS6dW/yqbsnE6oFoihKG6JBBRJDwsRWg2/bVnzBIMFSO/7hpnePiXy7lG0mhPKOG5fb4PHrrtuVIkkURWntxJrKpNXjq6wEIFDatPEPUAtEUZS2Sawz0Vs9/kobwhtshgLJBAtEUZT6BIIBFmz8ml111ezXcSB5WdGz9CpNQxWIg2uBuC6spqAWiKJkJoFggFcXvczCzXYh1NnrZnGunJdmqVoP6sJy2O3CamIILwCuAlELRFEyii/WzWLh5m/pWdKLAzsdzIad65m8/N10i9VqUAvEwbfVtUDicGG5y9qqBaIoGUNNXQ0zV80gLyuP0/udQUF2Aet3VGA2fc0X62ZxSOfB6RaxxaMWiIO/cgvBggLIbTiKKSLZ2QQLC9UCUZQMYv7Gueys28mgzkMoyinC7/NzdPdjAbhvdqtL5ZcW1AIBCAbxVVYS6NAx7ksE2pfhz6CU7qtW+ZgzJ4usLBg0qI7y8laXOEBRGmTe+rn48HFIecjS6FXSm04F5by15HU27dxIWX6HNErY8lELBDt24autJRjjKoSRCLZrnzG5sL75xs9//5vDvHlZfPllFk8/ncPSpb50i6UoKWPVtpWsqlpJz5JeFOeG3mufz8f+HQayK7CL1xa9kkYJWweqQAD/2rUAzVIggbIyGwpcl95kJpWV8Oab2WRnw1ln7WLUqBqCQXj99Ry2bUuraIqSMl5f/CoA+5TJHsf27bg/Pny88M2zqRar1aEKBPCvcxRIUfMsEEh/JNb772dTU+PjhBNq6ds3iEiA446rZccOH++9px5LpW3w+uLXANi7bM9sTKW5pRy+1xF8vPpD1u9Yn2rRWhWqQAD/2jVA8xRIoKzMXiuN4yDG+FmwIIsuXQIccEAo2/6gQQG6dAmwYEEWX36pP7nSuqmqqeLTNR/TpbArRTmR3+mT+pxCkCBTl09KsXStC21N8LqwiuK+RrDMDsb5NmxIiEzx8MgjOQAMG1aHzzPk4fPBMcfUAnDvvXFEmSlKC+Lj1TOpCdTQu7RP1DIn9joZgMnL30mRVK0TVSAkyALpVG6vlSYFsm0bPP98Du3aBenff8+1vnr3DlJeHmDixGxWrdIBdaX1Mm3Fe4CNuIrGvh32o3txD6Ysn0RdQBdhiBdVIHjHQOK3QAIdbQiwf0N6fKqTJmWzY4eP/fevwx/hV/X5YPDgOurqfDz+eE7qBVSUFDF95TTysvLoXtIjahmfz8eIXiezuXozn6/9LIXStS5UgZCYQfTc96YAkPfS8wmRqam89podIN9nn+grDe+7b4CysiDPPJNDbW2qJFOU1LFt11bmb5jLoM5DyPE33FEa0fskAKau0HGQeFEFgnVhBQsKmrwSoZdgYSEAvu3bEyVWzFRVWQtkwIA6OnWKPmEwJwfGjKmhosLPtGnxf1dFyVS+qJhNIBhgcJdDGy17ZLej8Pv8fLByegoka52oAsEOogcL43dfAQQLHAWyI/UKZMoU674aNaq23uB5JM45pwaA555TN5bS+pjluKOGxKBA2uW158BOB/P52k/ZXpP697Y1oApkxw6bB6sZkwghZIGQBgvEdV+dfnrjfqlBgwIMGFDHm29m4yQgVpRWgzueMbhz4woE4Ojux1ITqOHTNR8nU6xWS1IViIgcLiLvOdsDRGSGiEwXkQdEJCOUlzv+EWjGADoAOTkEc3JS7sLasQPeeSebPn3qz/2Ihs8H55xTy86dPiZO1ImFSushGAwya+1ndCnsSrfi7jGdc3T3YwCYsfL9ZIrWaklaIy4i1wGPAvnOrruAG40xxwA+4Ixk1d0Uds8BacYAukuwoDDlCmTq1Gy2b/cxalRNo+4rlzPOsG6s115TN5bSeli1bSVrt69hcJdD8cX4Mhy+1xFk+bJUgcRJMrugi4Azgf84n4cA05ztN4GTgZcbukBZWSHZ2XsO9paXlyROyp12KdvcsnbkFjVzqcuSYlizhvJOVhklVM4ovOusjXPRRXmUl+fRmCFVXp5HeTkMGmTTnmzalBo5E4HKmVham5zTKuYBcGy/oygvL6GogffZvWY5JQztPpRPV35KfimU5MV/T1rK/UwkSVMgxpgXRaSPZ5fPGOOGCG0FGl07dtOmPXvz5eUlVFRsTYiMAPnfLqUE2JmbT11VdbOulZeXT3ZdHeuXrKJTv+4JlTMSO3fCK68U06tXkF69qqiogKqqhmeaV1TsAmDkyFxmz87jtdfg1FOTK2ciSPTvnixUzsTSFDmnfmujqaT4QCoqtlLVwPvsvebhnY/io+8+4vWv3mFE75OTLmcqSJUyS+U4hNdBXwJsTmHdUfGva/4sdJfdkVgVFc2+VixMm5bFtm0+Tj+98eircE4/3bqxnk/PtBVFSTiz1n6G3+fnkPJBTTrPXWRqhobzNplUKpDZIjLc2T4VyIhfK6FjIE4kVqpmo7tjGKNG1TT53P79g+y/fx3vvINGYyktntpALV9WfIGU7UdxbtN630O7Hk6OP0fHQeIglWE41wCPiEgu8DXwQgrrjsruPFjNSKTosluBrE++Atm1C956K5vu3QMMHtx49JXLuHEhF1dZWZBdu+DXv85n//3tNa67blfCZVWUZPP1xvlsr90e0/yPcApzChnSZSifrPmILdWbaZfXPvECtlKSqkCMMUuBYc72N8BxyawvHvzr1tmGP6f5WWpdK8ZVSslk+vQsKit9/OhHsUdfhbPPPgE++AC+/da/W4EoSkti3Ce3ATCnYjYA67av272vKRzV/Rg+Wj2TD1fN5JS+IxMqY2smI+ZipBP/2jUEOnch7lbYgzsZ0b9mVbOv1RgTJlj3VSyTB6PRsWOQjh1hyRI/NU33gilKxrB6m33n9iraK67zj+lu+7bTv3svUSK1Cdq2Aqmrw7++gkCXrgm5XLDY+l79q1cn5HrR2LkT3njDuq+GDm1eKmoRqK31sXx5234UlJbN6qrV5Phz6VjQKa7zD+16GIXZRUz7bmqCJWvdtOmpyP71FfgCAeoSpkCsBZK1JnkKZNy4XL75xs/WrT723beOO+9snutNBGbOhEWL/BHXEVGUTKe6rpoNO9fTs6QXfl98HaHcrFyO7HYUk5a/w6ptK2Oeyd7WadPdzt1pTDp3TswFc3MJ5uXhT6ICAViwwP5s++7b/Aa/Rw8oKAiyaJGfYPREvoqSsaypsu/bXkXdmnWdY3sOB+B9dWPFTNu2QJzB7kCXrviqmzeJ0CVYXIJ/dfLGQKqrYfFiP2VlATp3bn6L7/dD//4B5s7NYs0aXalQaXmsrmr6+EekgfaKHXb+1nsrpvDDfc9PjHCtnDauQBwLpEtXspYvS8g1g8XFZC1barMcxoE3zDYS8+f7qa21Kw8mYNwfgAEDrAJZuLBNG6RKCyWkQJpngXTK70RRTjHvf/cegWAgbndYW6JN36HdFkjnLgm7ZsBNC78q8VZIMAhz5mTh9wc56KDErePcq1eA7OygKhClxREMBlm9bTUlOSVNnkAYjs/no3dpb9bvqGD+hnkJkrB106ZbDHesIlFRWBCKxGLlyoRd02XlSh/r1/sZMCDQaNLEppCbC717B9iwwc+SJerGUloOlbsq2V5bRdfi5lkfLr1L+wIwbYVGY8VC21YgHhdWogiWltqNZYlxiXn5+GObmXjQoMRZHy4DBtgB+bffbtNeTaWFkSj3lUvvkj4ATNF10mOibSuQdWsIZmcT7NAhYdcMtGtvNxYvTtg1AVat8rFkSRY9ewbo2TPx4VL9+gWAIG+9pQpEaTm4CqRbghRIcW4xB5UfwoerZrClenNCrtmaabMKpHDcbWR98w3BgkIK7/xrwq4bbN/ebiRQgQSDdu0OgCOPjH/meUMUFcFeewX5+OMsTa6otBhWbVuJ3+enS1HivAin9BlJbaCWKcvVCmmMNqtACAbxba8imMjBBCBYUkrQ70+oApk718933/np378uKdaHS9++AerqfEybplaIkvnsrN3J2u1rKC/oTI4/catrntL3NADeWvJ6wq7ZWmm7CmTnTnx1dbtnjyeMrCwCPXomTIFs2gTvvZdNTk6QESOSY3249O1rx0EmT95zFUhFyTS+rJhDIBige4JnjQ/seAA9S3oxafm77KrT7NQN0WYViG/bNiAx64CEU9ertw3jjXMuiMuuXTZpYnW1jxEjanHH55NF165BOnUKMHlyts5KVzKez9Z+ApDwtCM+n49T+oxk665KZq6akdBrtzbargKpchVIYl1YAHW9+wCQtWJ53NcIBm3CxA0b/AweXMsBByQ/T5XPB8cfX8fatX7mzm2zj4bSQvh0zccA7FWU+LxVrhvrjcWvJfzarYk220r4XQsk0S4soK6PjSXPWrQw7mt8+GEWCxfaqKvjjkt82G40XDfZlCk6DqJkLsFgkM/WfEJRTjGluYk3zY/odhSdCjoxcfEr1NTpWgfRaLMKJGSBJEGB7LMvAFnfmrjOX7TIz8yZWZSWBhk1qoasFA5JDB9ei98fZNIkHQdRMpfvtq1g7fY1dCvqhi9ROX08ZPuzOWPAmazfsZ7pK99L+PVbC6pAkqFARADINguafG5lpXVdZWfDGWfU4KySmzI6dIDBgwN89lkWmzentm5FiZVPVn8EJH78w8uZe58NwIvfPJ+0Olo6bVeBbKsCErMWejh1vfpAXl6TLZBgECZNyqa62scJJ9TSpUt6RrJHjKjVcF4lo3EHt3uU9EpaHYd2OYxepX14Y8lEttdsT1o9LZk220L4tlYS9PsJFiZegZCdDSJkf/MNBAI2Z3oMTJyYzeLFWfTqFeDAA9OzuNO4cbm707rffXcuxtSX/brrNKxRST8zVr5PSW4pXQoTlwg1HJ/Px1l7/4B/fH4n7yx9k9F7n5W0uloqbdcC2VppEx/G2Lg3mf33x7e9Cv/K72IqXlsLt92Wh98f5KSTahKWqj0eunQJUlgYZOlSXWRKyTxWbVvJki2LGbbXEUlPuf6DfX4IwDML/pPUeloqbVOB1NTg27YtlPgwGey/PwDZ38Q2DvLCC9ksWuTngAMClJUlT6xY8PmgT58AVVU+1q3T7LxKZvHByukAHNX92KTXtXfZPhy+1xG8t2IKyyqXJr2+lkabVCD+1avwYdOOJI2BAwHI+vrrRosGAnDPPbnk5gYZNiy5s81jxSZXtKsfKkom4Y5/HNXt6JTUd8F+Pwbgma+fTEl9LYk22TpkrbJrdQSSaYEcfDAA2fO+arTotGl2zsfo0cmfbR4rvXsH8PmCLFnSJh8RJUMJBoNM/24apbntOKDTQSmpc1T/0ZTmtuOZr5+iNpAZHbxMoU0Oovu/WwFAsKR5K5g1SN++BIqKY1Igjz1ml7G99NJdvPtuZvwkBQU2O+/q1T527LCfFSXdfLPJsHzrMr7ffwxZ/uTMVYq0Xnr/9v2ZvW4W7y57m1OdWepKGiwQEZktIu85f4+nun4Av2OBJMuFVTjuNrjlFoLt25NlFsDOnVHLrl7t4913sxg8uI5Bg9ITeRWNvn0DBIM+li1TK0TJDCYteweAE3ufnNJ6D+p0CABPzhuf0noznZS2DCKSD2CMGe78XZLK+l2yHAskqS4sIFDeGV8wSLaJPg4yYUI2waCPc8/NvHQJ7jiIurGUTGHSsrcBOKHXSSmtt7ywM92KujN5+bss2vxtSuvOZFLdMhwMFIrIOyIyRUSGpbh+ALKWLQUgWNouqfUEOncGIHtudDfWhAk5ZGUFGTUq83yrnTvbcN4lSzScV0k/ldVb+HjNhwzqPJjOhZ1TXv/gLocC8OhXD6W87kwl1Q737cCdwKPA3sCbIiLGmIitZ1lZIdnZe/o5y8ubOXaxfCkUFVHUIbkWSF6vHgCULFpASQSZFy6E2bPhe9+D/fazKVWSkBy4UYqK8qIe23tvmDMHKivz6NYNysujl002zf7dU4TKmVhcOafOe5PaQC3f32/U7n0NPbuJZlDhQcxZP4v/LXiaO0feQfv89hHlbEukWoF8Ayw0xgSBb0RkA7AXsCJS4U2b9kwfUF5eQkXF1vglqKmh09KlBLruxc6q6viv0whFRXlUFbWj0Oej5rNZbIkg8/jxuUAep522g4oKq0OrqnKTJlNUORu4Dz17+pkzJ4f582tp166Oior0zERv9u+eIlTOxOKV88lZTwNwfNdTdu9r6NlNBn1K+rOicgVjnj6LQ7setnt/UVEeVwy8JqWyNESqlFmqXVg/Af4OICLdgFJgdSoFyFqxDF9dHQF37fJkkpNDsKyDdWEF9hwgf/nlbPLygowcmXnuK5c+fTScV0k/ldVbmLTsbaRsX/brsH/a5Dio/GCy/dnMWvc5gWBmBb2kg1S3Co8B7UVkBvAs8JNo7qtkkbV4EQDB9qmZ7h3o3AX/tq34ly+rt3/+fD8LFmSlZKXB5pCfD926hcJ5FSUdvLnkdarrqhmz9w+Skr49VgqyCxjY8QAqd21h4eZv0iZHppBSF5YxZhdwXirrDCdriV2rPFDWISX1Bco7w4L5ZM+byy5noSmw0VcAZ56ZudaHS9++AVauzNZwXiVtvLzwBYCMSGg4uPNQ5lR8wcerP2Tv9pJWhZZuMmPWWgpJvQXiRmJ9ya7TRtm6g/DyyzkUFQU58cSWoUBmzNC0Jkp6WF65jKnLJzOky6H0a9c/3eLQsaAjUrYvZtMCllQublSmSBMTXa477IZEi5dS2lyLkPX1fII+H4EOKbJAXAXimZE+a5afZcv8nHJKbcoXjIqHzp2DFBXZcZAIQzmKklT+M/8JggT58cBL0y3KbobtdSQAH62aSbANx7i3LQUSDJI9fy51fftBbmqinYJFxQTKO5M9b+7ufRMm5AAwZkzmTR6MhJudd8cOH3PmtK1HRkkv1bXVPP31vynLK+OMAWemW5zdlBd2pn/7AayqWtmms/S2qdbAv3oV/s2bqdv/gJTWW3vAgWStWI5v8ybq6uz4R/v2QYYPr0upHM2hf39rerzxRpvzeipp5H9z/8f6Hev50X4XUpCdWQnZjup2DADvfTelzUZktanWIHu+tQJq9x8YMaw2WdQOPJDcqZPJnjeX6cHhrF3r54ILdqXKCEoIffsGyM4OMnFiDjfcsCutC14pbYPaQC23Tr+VHH8OtXU1DY4lpIPOhV04oNNBzF3/JbNWz4LU9kszgjZlgWTP+hyA2oMOTmm9bn3ZX87hxRetzj7rrMwfPPeSk2NzYy1a5GfBgjb12Chp4qVvn2fhxoX8aN8LKc1LbtqheDm627Hk+HOZvGQya6pSOqUtI2hTLUHOJx8DUHPoYY2UTCw1Bw8CwDdrFq++mkO3bgGOOKLluK9c9tnHWm0TJ7Ypw1VJA9V11dz56V/J8edw5ZDfpFucqBTnFnNcj+HsrN3JlVN+QV2g5b3XzaHtKJDaWnI+/IBAh44UPJraZGiBPn0JtGtP7UdfsHWrjzPPrEnaUuzJpF+/AHl5QVUgStJ5aM6/WFq5hMsPvZyeJb3SLU6DHFw+iAFlA5i6YjJjP7yxTUVltZmWIPvLL/DV7KK2e/fUV+7zUXvwINq9P5V2bGbr1kLGjWtBAyAOubkwfHgdb7+dzcKFPgYMaDsvipI6Vm1byV2fjaNTQSduPv5majI8ZZfP5+PM/c7k5fkTeGjOv9hes51bj/5rxg36J4MW2A+Oj9x37ToCdX37paX+bftaN9bxpZ9RXt5yG1439PjZZ3PSLInSGgkGg1w19Qq2127npmG37JHxNlMpyCnghe+/ysCOB/Kf+Y9z1DOHctdn45ix8n2qara1WqukzVggue+8RdDvp65338YLJ4HJWw7lbOB7HT9jEcekRYZEcOqptZSWBnnuuRyuv34XWclZVVRpo4yf+zDvrZjCiF4n8cN9z0+3OE1ir+JuvHbm29z12Tge+fIB/vrJrbuPZfuzKc1tR6eCcvYp24e920vSluRNJW3CAsla8DU5X82hrncfyEv9ehaBANwz066ddZj/s5TXn0gKCmD06BpWr/YzbVrLfwGUzOHrDfO5eeZNdMjvwN3H/6tF5pgqzinmj0fcwryLF/LwSY/z60G/Ye8yoWN+R6pqtvHNpgVMXPwqj819mMWbF6Vb3GbTJiyQ/KefBKD2gIPSUv/772cxY0UfKvM60Xv952mRIZH86Ec1PPlkLs88k8MJJ7StqBMlOWzauZGL3vwhO+t28uBJ4+lS1DXdIjWL0rx2jN77LEbvfRa5WXa8MxgMsmHnBuZUzGZOxWxeWvg8PUp6csOwP+L3tcy+fMuUugn4Nm4g/6l/U9e5C3X9B6RFhnvvzQV8VB80mI5bllJUtS4tciSKwYMD7LdfHa+/ns3KlS2vl6hkFrWBWi575xKWVS7l6iHXMrLf6ekWKSn4fD46FXRiRK+TuHC/iynLK+Oe2XdxxaTLWmz4b6tWIIXjbqP04vPxV22jduABkJ16g+uTVT2ZPj2b43t/S7u8nQD0W/lByuVIJD4fXH75LurqfDz8cMuLJlMyh2AwyE0fXM+076Zycu9T+N1hN6ZbpJRQXtiZH+17IYd2OYwXv32O696/ukUOtLdqBUJVFTmff0agsIjagw5JefXBINw280QArhs2lboePQHot2J6ymVJNGPG1NKlS4D//CeHysp0S6O0VO7+/E4e++ph9uswkPtPfKTFunLioTCnkP+d/iIHdjqY/8x/gts+viXdIjWZVj0GkvvhB/hqdrHr2OE2F0eKeWvxvkxaug/Dey3k6B5LCNTtRZ0vm37fzUi5LIkmLw9+9rMabr01j0cfzeU3v0nPWulKy8Kbz+rLii94Z9lblOaWclyP4Tw45197lC8qykv5uueppDSvHf87/SVGvXwy/5z1d/q3H9Cios9arbrPWvQt2V9+QaCsLOW5rwC21+Rw/dTTyPLVMe6E12zywexsNrXrRbut36VcnmRwySW76NgxwH335bJ+vY6FKLHz7aZveHfZ2xRkF/CDfX5IcW5JukVKG+WF5Txz2vO0z2vPNe/9mg9XtRwXd6u1QIr+cgu+QIDqY44jlZMVPvjA1vWvZSezaHMnztxrBpvMBj4wdn/OgDP46JCfpUyeRBM+g/7ggwNMmZLNeecVcOKJtVx3nVoiSsMs3rKIiYtfIdufzZl7n0OH/NQs7pbJ9Gs/gPGnPMU5r43mkrfO582zptC3XXomPTeFVmmB+DZuIG/iK9Tt1Y26vSXl9b+3/gBeXTWUvoVruKTXpHrHanIK2dg+8x+MWDn44DrKygJ88YWf1avVClEaZsmWxbyy8CV8+Bg94Cz2Ktor3SJlDEd3P5Y7jr2LjTs3cuEb51JZvSXdIjVKq7RAgmUdqLz3QbK/mkOqF66YW9mbOxeeRWFWNb/f+1ly/S0rbXtTycqCk0+u5dlnc3nrrWxuv72a/Px0S6VkItNWTOWVhS8BMHrAWfQu7ZNegRJMPOuVRDpnSJehfL72U66cegWPn/JUIkRLGq3SAsHno/rc8wiWpnYNgRkr+nDj1xdRF/Qzdv/n6F1YkdL600XPnkEOOaSODRv8PPCAhvUqkfnd+78hSJDRA86iT7v0pBRqCRzX43j27zCQumDmzw1plRZIOqjYXsSZL15CTdDPH/Z5lqEdFrGrDQ0HHHdcLYGAdWkpSiTGHfcPJnz7El1b+CzzZOP3+RnZbxTXHXZDukVpFFUgCaI4p5oz9pnLAcEvObT9Qtrarc3Jsa4sTW2iROPYHsP5aNXMdIuhJJDW6cJKAwU5tTwy8nlHeSiKorR+VIEoiqIocZFSP4uI+IH7gYOBauCnxhjtsiuKorRAUm2BjAbyjTFHANcDf09x/YqiKEqCSLUCORp4C8AY8xFwaIrrVxRFURKEL5UphEXkUeBFY8ybzuflQD9jTOuebacoitIKSbUFUgl4s6b5VXkoiqK0TFKtQD4ARgKIyDDgqxTXryiKoiSIVM92exk4SURmAj7gkhTXryiKoiSIlI6BKIqiKK0HnUioKIqixIUqEEVRFCUuVIEoiqIocZGRKWNjTXkiIg8DG40x16dYRLf+BuUUkd8AlwLuwiA/N8aYDJRzKHAXNrBhDXCBMWZnJskpIl2B/3mKHwJcb4x5MJPkdI6fD1wD1AHjjTEPpFpGjyyNyXoh8FtgC/CEMeaxtAhqZTkcuMMYMzxs/yjgj0At9n4+kgbx6hFNVudYIfAucKkxZkGqZUslmWqBjKaRlCci8nPgwBTLFc5oGpZzMHCRMWa485dy5eEwmihyiogPeAS4xBjjZgronQ4haUBOY8wa9z4CvwdmYeVOB6Np+He/EzgROAq4RkTKUitePUYT/bfvBNwKDAeOA84XkT6pFxFE5DrgUSA/bH8O8A/gZKyMlzmdibQRTVbn2KHA+0D/VMuVDjJVgTSY8kREjgCGAQ+lXrR6NJaaZQjwexGZISK/T7VwHhqScx9gA3CViEwDOqRR0TWa6sZRePcClxtj0rX4SGNyfgm0wzYwPiCdoY4NydoP+MIYs9EYEwA+xb5X6WARcGaE/fsBC40xm4wxu4AZwDEplWxPoskKkAeMAVq15eGSqQqkFGtSu9SJSDaAiOwFjAWuSINc4USV0+F/wP8BJwBHi8jpqRTOQ0NydgKOxLo5TgRGiMiIFMvn0tj9BBgFzEujkoPG5ZwLfA7MAyYaYzanULZwGpL1W2CgiHRx3C4jgKJUCwhgjHkRqIlwKFz+rVjlnDYakBVjzAfGmBUpFiltZKoCaSjlydnYRu8NrEl+nohcnFrxdhNVTqenfLcxZr3Tc3odGJQGGaHh+7kB28Obb4ypwfZWh6RaQIdYUt1cADycOpEi0tDvfhBwGtAX6AN0FpGzUy5hiKiyGmM2AVcDLwLjsW7B9SmXsGHC5S8BNqdHFCWcTFUgUVOeGGPuMcYMcXzhfwWeMcY8kQ4haTg1SykwV0SKHWVyArZXmg4aknMxUCwiA5zPx2B7zukgllQ3Q4B0r4vakJxbgB3ADsfFtg5I5xhIVFkdS2QYcCxwEbCvUz6T+BrYW0Q6iEguVtYP0yyT4pCRUVhESHkiIucBxcaYdPc+vTQop4jcAEzFRr9MNsa8kaFyXgo84yi6mcaY1zNUznJgqzEm3ekTGpPzIWCGiOzC+sufSJ+ojcq6C9ux2Qn83RiTERZImIy/Ad7GdnjHG2NWple6+mRo25QSNJWJoiiKEheZ6sJSFEVRMhxVIIqiKEpcqAJRFEVR4kIViKIoihIXqkAURVGUuFAFoiiKosRFps4DUZS0ICInY+ccABxljJnp7H8C+HEDp7rLMz/eQJmyNKc1UZSEogpEUepznmf7QkKz3mcB7Z3twUBP4DPAndS2HOjlbC8F5kS4dsT8SYrSUtGJhIriICIFwFpsRtVcYBPQ1cll5i33BNYa+ZEx5n+e/RdjLZCHjDH/lyKxFSVt6BiIooQYhU3W9zbWuigD0pVBWVEyHnVhKUqI853/LwNdsWtnXAS81MTrfE9EJoTt+9gYc3vzxFOUzEIViKIAItIBOAW7DO1rQEfgNmCkiHQ0xmxowuX6OH+K0qpRBaIolrOx4x7TnIy060VkPrA/8EPgX024lo6BKG0CVSCKYnHdV8eJSHhkyYU0TYEoSptAB9GVNo+I9MKuHR7ALqbl/QM4XET2SZN4ipKxqAWiKHbuhw94yxgz0ntARKYCw7FWyE0xXi/SIDrAjcaYuc2QU1EyCrVAFCXkvhof4dhDzv8LnBUbY6EPcEaEv07NkFFRMg6dSKgoiqLEhVogiqIoSlyoAlEURVHiQhWIoiiKEheqQBRFUZS4UAWiKIqixIUqEEVRFCUuVIEoiqIocaEKRFEURYmL/wcVudAsPsCoOgAAAABJRU5ErkJggg==",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Hide cell\n",
"bootstrap_sample = 100\n",
"AIPW_ates = []\n",
"IPW_ates = []\n",
"naives_ates = []\n",
"for iSample in range(bootstrap_sample):\n",
" df_bootstrap = df_categ.sample(frac=1,replace=True)\n",
" ate, ps, mu0, mu1 = AIPW(df_bootstrap, X, T, Y,true_mus=False,true_ps=False)\n",
" AIPW_ates.append(ate)\n",
" ate, ps,_ = IPW(df_bootstrap, X, T, Y,true_ps=False)\n",
" IPW_ates.append(ate)\n",
" naives_ates.append(naive_ATE(df_bootstrap, T, Y))\n",
"\n",
"sns.histplot(AIPW_ates, kde=True, stat=\"density\", linewidth=0,color='g')\n",
"sns.histplot(IPW_ates, kde=True, stat=\"density\", linewidth=0,color='b')\n",
"sns.histplot(naives_ates, kde=True, stat=\"density\", linewidth=0,color='r')\n",
"# plt.vlines(np.percentile(AIPW_ates, 2.5), 0, 20, linestyles=\"dotted\")\n",
"# plt.vlines(np.percentile(AIPW_ates, 97.5), 0, 20, linestyles=\"dotted\", label=\"95% CI\")\n",
"plt.title(\"ATE Bootstrap Distribution\\nRandom Outcome and Exposure Model\", fontsize=20, fontweight = 'bold')\n",
"plt.xlabel('ATE', fontsize=14, fontweight = 'bold', labelpad=5)\n",
"plt.legend(['AIPW', 'IPW', 'Naive'])\n",
"\n",
"Results = {\"AIPW\":{\"Mean ATE\":np.mean(AIPW_ates), \"Std Dev\":np.std(AIPW_ates), \"[.025\":np.percentile(AIPW_ates, 2.5), \".975]\":np.percentile(AIPW_ates, 97.5)},\n",
"\"IPW\":{\"Mean ATE\":np.mean(IPW_ates), \"Std Dev\":np.std(IPW_ates), \"[.025\":np.percentile(IPW_ates, 2.5), \".975]\":np.percentile(IPW_ates, 97.5)},\n",
"\"Naive\":{\"Mean ATE\":np.mean(naives_ates), \"Std Dev\":np.std(naives_ates), \"[.025\":np.percentile(naives_ates, 2.5), \".975]\":np.percentile(naives_ates, 97.5)}}\n",
"\n",
"df_results = pd.DataFrame(Results)\n",
"df_results.T"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Concluding Thoughts for AIPE\n",
"\n",
"We clearly demonstrate AIPW's doubly robust properties using the simulated National Mindset dataset. But it is important to note, our 'incorrect' models were uniform random which would be about as poor as one could imagine. In reality, misspecified models contain more subtle biases or noise, and thus there is a whole host of literature investigating the sensitivity of doubly robust methods to various types and degrees of misspecification. For instance in the example where both models were incorrect, one could imagine scenarios where model misspecifications cancel out, and actually produce a relatively accurate ATE estimate. It is an area of active research on when doubly robust methods should be used when there might be uncertainty on both models."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### TMLE\n",
"\n",
"#### Background and Setup\n",
"Targeted Maximum Likelihood Estimation (TMLE) is a semi-parametric method with minimal assumptions on the underlying data distribution demonstrated by Van der Laan & Rubin in 2006. We will briefly walk through the steps of a TMLE estimation algorithm on the same data without diving too deep into the formulation."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"# Hide Cell\n",
"import numpy as np\n",
"import pandas as pd\n",
"pd.set_option(\"display.precision\", 3)\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('seaborn')\n",
"import seaborn as sns\n",
"from scipy.special import logit, expit\n",
"import statsmodels.api as sm\n",
"\n",
"# Super Learner Import\n",
"from math import sqrt\n",
"from numpy import hstack\n",
"from numpy import vstack\n",
"from numpy import asarray\n",
"from sklearn.datasets import make_regression\n",
"from sklearn.model_selection import KFold\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn.metrics import mean_squared_error\n",
"from sklearn.linear_model import LinearRegression\n",
"from sklearn.linear_model import ElasticNet\n",
"from sklearn.neighbors import KNeighborsRegressor\n",
"from sklearn.tree import DecisionTreeRegressor\n",
"from sklearn.svm import SVR\n",
"from sklearn.ensemble import AdaBoostRegressor\n",
"from sklearn.ensemble import BaggingRegressor\n",
"from sklearn.ensemble import RandomForestRegressor\n",
"from sklearn.ensemble import ExtraTreesRegressor"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The TMLE algorithm begins by first estimating a model by training and predicting a super learning ensemble of algorithms. In the hidden code block below, we do so using 9 models from pre-built libraries. We report the root mean squared errors for all algorithms, demonstrating the super learner ensemble performs best. "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Train (5195, 31) (5195,) Test (5196, 31) (5196,)\n",
"Meta (5195, 9) (5195,)\n",
"LinearRegression: RMSE 0.828\n",
"ElasticNet: RMSE 0.990\n",
"SVR: RMSE 0.849\n",
"DecisionTreeRegressor: RMSE 1.023\n",
"KNeighborsRegressor: RMSE 0.887\n",
"AdaBoostRegressor: RMSE 0.827\n",
"BaggingRegressor: RMSE 0.900\n",
"RandomForestRegressor: RMSE 0.898\n",
"ExtraTreesRegressor: RMSE 0.954\n",
"Super Learner: RMSE 0.808\n"
]
}
],
"source": [
"# Hide Cell\n",
"# create a list of base-models\n",
"def get_models():\n",
"\tmodels = list()\n",
"\tmodels.append(LinearRegression())\n",
"\tmodels.append(ElasticNet())\n",
"\tmodels.append(SVR(gamma='scale'))\n",
"\tmodels.append(DecisionTreeRegressor())\n",
"\tmodels.append(KNeighborsRegressor())\n",
"\tmodels.append(AdaBoostRegressor())\n",
"\tmodels.append(BaggingRegressor(n_estimators=10))\n",
"\tmodels.append(RandomForestRegressor(n_estimators=10))\n",
"\tmodels.append(ExtraTreesRegressor(n_estimators=10))\n",
"\treturn models\n",
" \n",
"# collect out of fold predictions form k-fold cross validation\n",
"def get_out_of_fold_predictions(X, y, models):\n",
"\tmeta_X, meta_y = list(), list()\n",
"\t# define split of data\n",
"\tkfold = KFold(n_splits=10, shuffle=True)\n",
"\t# enumerate splits\n",
"\tfor train_ix, test_ix in kfold.split(X):\n",
"\t\tfold_yhats = list()\n",
"\t\t# get data\n",
"\t\ttrain_X, test_X = X[train_ix], X[test_ix]\n",
"\t\ttrain_y, test_y = y[train_ix], y[test_ix]\n",
"\t\tmeta_y.extend(test_y)\n",
"\t\t# fit and make predictions with each sub-model\n",
"\t\tfor model in models:\n",
"\t\t\tmodel.fit(train_X, train_y)\n",
"\t\t\tyhat = model.predict(test_X)\n",
"\t\t\t# store columns\n",
"\t\t\tfold_yhats.append(yhat.reshape(len(yhat),1))\n",
"\t\t# store fold yhats as columns\n",
"\t\tmeta_X.append(hstack(fold_yhats))\n",
"\treturn vstack(meta_X), asarray(meta_y)\n",
" \n",
"# fit all base models on the training dataset\n",
"def fit_base_models(X, y, models):\n",
"\tfor model in models:\n",
"\t\tmodel.fit(X, y)\n",
" \n",
"# fit a meta model\n",
"def fit_meta_model(X, y):\n",
"\tmodel = LinearRegression()\n",
"\tmodel.fit(X, y)\n",
"\treturn model\n",
" \n",
"# evaluate a list of models on a dataset\n",
"def evaluate_models(X, y, models):\n",
"\tfor model in models:\n",
"\t\tyhat = model.predict(X)\n",
"\t\tmse = mean_squared_error(y, yhat)\n",
"\t\tprint('%s: RMSE %.3f' % (model.__class__.__name__, sqrt(mse)))\n",
" \n",
"# make predictions with stacked model\n",
"def super_learner_predictions(X, models, meta_model):\n",
"\tmeta_X = list()\n",
"\tfor model in models:\n",
"\t\tyhat = model.predict(X)\n",
"\t\tmeta_X.append(yhat.reshape(len(yhat),1))\n",
"\tmeta_X = hstack(meta_X)\n",
"\t# predict\n",
"\treturn meta_model.predict(meta_X)\n",
" \n",
"# create the inputs and outputs\n",
"X = df_categ[df_categ.columns.drop([Y])].to_numpy()\n",
"y = df_categ[Y].to_numpy()\n",
"# split\n",
"X, X_val, y, y_val = train_test_split(X, y, test_size=0.50)\n",
"print('Train', X.shape, y.shape, 'Test', X_val.shape, y_val.shape)\n",
"# get models\n",
"models = get_models()\n",
"# get out of fold predictions\n",
"meta_X, meta_y = get_out_of_fold_predictions(X, y, models)\n",
"print('Meta ', meta_X.shape, meta_y.shape)\n",
"# fit base models\n",
"fit_base_models(X, y, models)\n",
"# fit the meta model\n",
"meta_model = fit_meta_model(meta_X, meta_y)\n",
"# evaluate base models\n",
"evaluate_models(X_val, y_val, models)\n",
"# evaluate meta model\n",
"yhat = super_learner_predictions(X_val, models, meta_model)\n",
"print('Super Learner: RMSE %.3f' % (sqrt(mean_squared_error(y_val, yhat))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the second step we use the super learner algorithm to estimate the expected value of the outcome using the treatment and confounders as predictors. Within this, there are three steps:\n",
"1. predict with the interventions\n",
"2. predict with every sample receiving no treatment\n",
"3. predict with every sample receiving treatment.\n",
" \n",
"We can take the difference of the last two as an ATE estimate, which is effectively the g-estimation approach. We see below this provides a decent 1st estimate. "
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.4044480650676532"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_predict = df_categ.copy()\n",
"X = df_predict[df_predict.columns.drop([Y])].to_numpy()\n",
"Q_a = super_learner_predictions(X, models, meta_model)\n",
"df_predict['intervention'] = 0\n",
"X = df_predict[df_predict.columns.drop([Y])].to_numpy()\n",
"Q_0 = super_learner_predictions(X, models, meta_model)\n",
"df_predict['intervention'] = 1\n",
"X = df_predict[df_predict.columns.drop([Y])].to_numpy()\n",
"Q_1 = super_learner_predictions(X, models, meta_model)\n",
"\n",
"df_tmle = pd.DataFrame([df_categ[Y].to_numpy(),df_categ[T].to_numpy(), Q_a,Q_0,Q_1]).T\n",
"df_tmle.columns = ['Y','D','Q_a','Q_0','Q_1']\n",
"\n",
"df_tmle['Q_1'].mean() - df_tmle['Q_0'].mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the third step we obtain propensity scores (ps) and form a \"clever covariate\" from these values which will be used to refine our model. the inverse ps with indicator is added with the negative inverse of not being treated (1-ps) also multiplied with indicator if not being treated:\n",
"\n",
"$$H(D,X) = \\frac{I(D=1)}{\\hat \\pi (X_i)} - \\frac{I(D=0)}{1 - \\hat \\pi (X_i)}$$"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"T = 'intervention'\n",
"Y = 'achievement_score'\n",
"X = df_categ.columns.drop([T, Y])\n",
"ate, ps,_ = IPW(df_categ, X, T, Y)\n",
"\n",
"df_tmle['H_1'] = 1/ps\n",
"df_tmle['H_0'] = -1/(1-ps)\n",
"df_tmle['H_a'] = df_tmle['D'] * df_tmle['H_1'] + (1-df_tmle['D']) * df_tmle['H_0']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the fourth and fifth steps, we estimate the fluctuation parameter using the logit function:\n",
" \n",
"$$logit(\\mathbb{E}[Y|D,X]) = logit(\\hat{\\mathbb{E}}[Y|D,X]) = \\epsilon H(D,X)$$\n",
" \n",
"We then update out initial estimates with the fluctuation parameter adjustment.\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"eps_fit = np.polyfit(df_tmle['H_a'], df_tmle['Y'] - df_tmle['Q_a'], 1)[0]\n",
"df_tmle['Q_0_hat'] = df_tmle['Q_0'] + eps_fit * df_tmle['H_0']\n",
"df_tmle['Q_1_hat'] = df_tmle['Q_1'] + eps_fit * df_tmle['H_1']\n",
"df_tmle['Q_a_hat'] = df_tmle['Q_a'] + eps_fit * df_tmle['H_a']"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"TMLE TAE estimate: 0.3864\n"
]
}
],
"source": [
"TMLE_ate = df_tmle['Q_1_hat'].mean() - df_tmle['Q_0_hat'].mean()\n",
"print('TMLE TAE estimate: {:.4f}'.format(TMLE_ate))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see how the fluctuation adjusted outcomes estimates vastly improves the ATE to be more in line with AIPW. One major benefit of TMLE is a whole set of nice statistical and convergence properties. In this case, we can use something called the influence function to calculate a closed form standard error, unlike the empirical error estimates we gained by bootstrapping in the AIPW case.\n",
" \n",
"$$\\hat{IF} = (Y-\\hat{\\mathbb{E}}*[Y|D,X])H(D,X) + \\hat{\\mathbb{E}}*[Y|D=1,X] - \\hat{\\mathbb{E}}*[Y|D=0,X] - ATE_{TMLE}$$\n",
"$$SE = \\sqrt{var(IF)/N}$$\n",
" \n",
"Using the above method, we see we get a SE very similar to our AIPW empirical methods."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SE calculated from influcence function: 0.0166\n"
]
}
],
"source": [
"IF = (df_tmle['Y'] - df_tmle['Q_a_hat']) * df_tmle['H_a'] + df_tmle['Q_1_hat'] - df_tmle['Q_0_hat'] - TMLE_ate\n",
"print('SE calculated from influcence function: {:.4f}'.format(np.sqrt(IF.var()/df_tmle.shape[0])))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### References \n",
" \n",
" 1. Glynn, A. N., & Quinn, K. M. (2010). An introduction to the augmented inverse propensity weighted estimator. Political analysis, 18(1), 36-56.\n",
"\n",
" 2. https://matheusfacure.github.io/python-causality-handbook/12-Doubly-Robust-Estimation.html\n",
"\n",
" 3. Gruber S, van der Laan MJ. Targeted minimum loss based estimator that outperforms a given estimator. Int J Biostat. 2012 May 18;8(1):Article 11. doi: 10.1515/1557-4679.1332. PMID: 22628356; PMCID: PMC6052865.\n"
]
}
],
"metadata": {
"celltoolbar": "Tags",
"interpreter": {
"hash": "e1b97f395c729d34f085d8e45c9fd7029ba997e22632cb4ee6e2174cac7541d0"
},
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.10"
}
},
"nbformat": 4,
"nbformat_minor": 2
}