Skip to content

Instantly share code, notes, and snippets.

@allisonmorgan

allisonmorgan/gender_subtopic.ipynb

Last active May 3, 2021 23:04
Show Gist options
  • Save allisonmorgan/914ec8237a381b8ed124d0e695c95894 to your computer and use it in GitHub Desktop.
Save allisonmorgan/914ec8237a381b8ed124d0e695c95894 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
gender_subtopic.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's say you wanted to learn the fraction of women across a variety of subtopics. Your data contains for each researcher, whether or not they are a woman and their primary subtopic. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Fraction female across dataset: 0.31\n",
"Topic distribution: {'ml': 0.3756, 'nlp': 0.2112, 'hci': 0.2078, 'ai': 0.2054}\n"
]
}
],
"source": [
"N = 5000 # Number of researchers\n",
"\n",
"# Subset of topics with different probabilities of men and women specializing in them\n",
"topics = ['nlp', 'ai', 'ml', 'hci'] \n",
"women_topic_weights = np.array([0.5, 0.2, 0.1, 0.2])\n",
"men_topic_weights = np.array([0.1, 0.2, 0.5, 0.2])\n",
"\n",
"data = pd.DataFrame({'is_woman': np.random.choice([0, 1], N, p=[0.7, 0.3]),\n",
" 'department_id': np.random.choice(range(1, 100), N)})\n",
"# This next line is just assigning a primary topic based on some knowledge about the likliehood of women or men to work in that are. Assuming men and women are unevenly distributed across topics.\n",
"data['primary_topic'] = data['is_woman'].apply(\n",
" lambda x: np.random.choice(topics, 1, p=women_topic_weights)[0] if (x == 1)\n",
" else np.random.choice(topics, 1, p=men_topic_weights)[0])\n",
"\n",
"print('Fraction female across dataset: %.2f' % (sum(data.is_woman)/N))\n",
"print('Topic distribution: ', dict(data.primary_topic.value_counts(normalize=True)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's see if we can recover the same distribution of fraction female in each topic via regression and data analysis."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import statsmodels.formula.api as smf\n",
"import statsmodels.api as sm"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.489842\n",
" Iterations 6\n"
]
},
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>Logit Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>is_woman</td> <th> No. Observations: </th> <td> 5000</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>Logit</td> <th> Df Residuals: </th> <td> 4996</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>MLE</td> <th> Df Model: </th> <td> 3</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Tue, 15 Dec 2020</td> <th> Pseudo R-squ.: </th> <td>0.2042</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>09:09:24</td> <th> Log-Likelihood: </th> <td> -2449.2</td> \n",
"</tr>\n",
"<tr>\n",
" <th>converged:</th> <td>True</td> <th> LL-Null: </th> <td> -3077.7</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>nonrobust</td> <th> LLR p-value: </th> <td>3.268e-272</td>\n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> -0.9088</td> <td> 0.069</td> <td> -13.178</td> <td> 0.000</td> <td> -1.044</td> <td> -0.774</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(primary_topic)[T.hci]</th> <td> 0.1307</td> <td> 0.096</td> <td> 1.361</td> <td> 0.173</td> <td> -0.057</td> <td> 0.319</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(primary_topic)[T.ml]</th> <td> -1.4649</td> <td> 0.108</td> <td> -13.609</td> <td> 0.000</td> <td> -1.676</td> <td> -1.254</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(primary_topic)[T.nlp]</th> <td> 1.7870</td> <td> 0.097</td> <td> 18.508</td> <td> 0.000</td> <td> 1.598</td> <td> 1.976</td>\n",
"</tr>\n",
"</table>"
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" Logit Regression Results \n",
"==============================================================================\n",
"Dep. Variable: is_woman No. Observations: 5000\n",
"Model: Logit Df Residuals: 4996\n",
"Method: MLE Df Model: 3\n",
"Date: Tue, 15 Dec 2020 Pseudo R-squ.: 0.2042\n",
"Time: 09:09:24 Log-Likelihood: -2449.2\n",
"converged: True LL-Null: -3077.7\n",
"Covariance Type: nonrobust LLR p-value: 3.268e-272\n",
"===========================================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"-------------------------------------------------------------------------------------------\n",
"Intercept -0.9088 0.069 -13.178 0.000 -1.044 -0.774\n",
"C(primary_topic)[T.hci] 0.1307 0.096 1.361 0.173 -0.057 0.319\n",
"C(primary_topic)[T.ml] -1.4649 0.108 -13.609 0.000 -1.676 -1.254\n",
"C(primary_topic)[T.nlp] 1.7870 0.097 18.508 0.000 1.598 1.976\n",
"===========================================================================================\n",
"\"\"\""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mod = smf.logit('is_woman ~ C(primary_topic)', data)\n",
"res = mod.fit()\n",
"res.summary()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'nlp': 0.7064393939393939,\n",
" 'ai': 0.2872444011684516,\n",
" 'ml': 0.08519701810436642,\n",
" 'hci': 0.3147256977863331}"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dict([(t, res.predict(pd.DataFrame({'primary_topic': [t]}))[0]) for t in topics])"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'ai': 0.2872444011684518,\n",
" 'hci': 0.314725697786333,\n",
" 'ml': 0.08519701810436635,\n",
" 'nlp': 0.7064393939393939}"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dict(data.groupby(['primary_topic']).mean()['is_woman'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Those two distributions looks exactly the same :check:. Alright, now let's try and run this as a regression. The unit of analysis is departments."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"departments = pd.DataFrame(data.groupby(['department_id']).mean()['is_woman']).reset_index()\n",
"\n",
"departments['nlp_researchers'] = departments['department_id'].apply(\n",
" lambda x: data[data.department_id == x].primary_topic.value_counts()['nlp'])\n",
"departments['ai_researchers'] = departments['department_id'].apply(\n",
" lambda x: data[data.department_id == x].primary_topic.value_counts()['ai'])\n",
"departments['ml_researchers'] = departments['department_id'].apply(\n",
" lambda x: data[data.department_id == x].primary_topic.value_counts()['ml'])\n",
"departments['hci_researchers'] = departments['department_id'].apply(\n",
" lambda x: data[data.department_id == x].primary_topic.value_counts()['hci'])"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.546718\n",
" Iterations 5\n"
]
},
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>Logit Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>is_woman</td> <th> No. Observations: </th> <td> 99</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>Logit</td> <th> Df Residuals: </th> <td> 98</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>MLE</td> <th> Df Model: </th> <td> 0</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Tue, 15 Dec 2020</td> <th> Pseudo R-squ.: </th> <td>-77.09</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>09:09:24</td> <th> Log-Likelihood: </th> <td> -54.125</td>\n",
"</tr>\n",
"<tr>\n",
" <th>converged:</th> <td>True</td> <th> LL-Null: </th> <td>-0.69315</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>nonrobust</td> <th> LLR p-value: </th> <td> nan</td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> -0.8246</td> <td> 0.218</td> <td> -3.777</td> <td> 0.000</td> <td> -1.252</td> <td> -0.397</td>\n",
"</tr>\n",
"</table>"
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" Logit Regression Results \n",
"==============================================================================\n",
"Dep. Variable: is_woman No. Observations: 99\n",
"Model: Logit Df Residuals: 98\n",
"Method: MLE Df Model: 0\n",
"Date: Tue, 15 Dec 2020 Pseudo R-squ.: -77.09\n",
"Time: 09:09:24 Log-Likelihood: -54.125\n",
"converged: True LL-Null: -0.69315\n",
"Covariance Type: nonrobust LLR p-value: nan\n",
"==============================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept -0.8246 0.218 -3.777 0.000 -1.252 -0.397\n",
"==============================================================================\n",
"\"\"\""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"simple_mod = smf.logit('is_woman ~ 1', departments)\n",
"simple_res = simple_mod.fit()\n",
"simple_res.summary()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"null_model_predictions = simple_res.predict()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's make a fancier model where we know how many researchers from each topic were hired. does that help us predict the fraction of women at the department level? This model is like yours but we are learning"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.541322\n",
" Iterations 5\n"
]
},
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>Logit Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>is_woman</td> <th> No. Observations: </th> <td> 99</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>Logit</td> <th> Df Residuals: </th> <td> 94</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>MLE</td> <th> Df Model: </th> <td> 4</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Tue, 15 Dec 2020</td> <th> Pseudo R-squ.: </th> <td>-76.32</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>09:09:24</td> <th> Log-Likelihood: </th> <td> -53.591</td>\n",
"</tr>\n",
"<tr>\n",
" <th>converged:</th> <td>True</td> <th> LL-Null: </th> <td>-0.69315</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>nonrobust</td> <th> LLR p-value: </th> <td> 1.000</td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> -0.8307</td> <td> 1.589</td> <td> -0.523</td> <td> 0.601</td> <td> -3.946</td> <td> 2.285</td>\n",
"</tr>\n",
"<tr>\n",
" <th>nlp_researchers</th> <td> 0.0336</td> <td> 0.059</td> <td> 0.567</td> <td> 0.571</td> <td> -0.083</td> <td> 0.150</td>\n",
"</tr>\n",
"<tr>\n",
" <th>ai_researchers</th> <td> 0.0094</td> <td> 0.075</td> <td> 0.125</td> <td> 0.901</td> <td> -0.138</td> <td> 0.157</td>\n",
"</tr>\n",
"<tr>\n",
" <th>ml_researchers</th> <td> -0.0287</td> <td> 0.052</td> <td> -0.549</td> <td> 0.583</td> <td> -0.131</td> <td> 0.074</td>\n",
"</tr>\n",
"<tr>\n",
" <th>hci_researchers</th> <td> 0.0084</td> <td> 0.073</td> <td> 0.115</td> <td> 0.908</td> <td> -0.135</td> <td> 0.152</td>\n",
"</tr>\n",
"</table>"
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" Logit Regression Results \n",
"==============================================================================\n",
"Dep. Variable: is_woman No. Observations: 99\n",
"Model: Logit Df Residuals: 94\n",
"Method: MLE Df Model: 4\n",
"Date: Tue, 15 Dec 2020 Pseudo R-squ.: -76.32\n",
"Time: 09:09:24 Log-Likelihood: -53.591\n",
"converged: True LL-Null: -0.69315\n",
"Covariance Type: nonrobust LLR p-value: 1.000\n",
"===================================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"-----------------------------------------------------------------------------------\n",
"Intercept -0.8307 1.589 -0.523 0.601 -3.946 2.285\n",
"nlp_researchers 0.0336 0.059 0.567 0.571 -0.083 0.150\n",
"ai_researchers 0.0094 0.075 0.125 0.901 -0.138 0.157\n",
"ml_researchers -0.0287 0.052 -0.549 0.583 -0.131 0.074\n",
"hci_researchers 0.0084 0.073 0.115 0.908 -0.135 0.152\n",
"===================================================================================\n",
"\"\"\""
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fancier_model = smf.logit(\n",
" 'is_woman ~ nlp_researchers + ai_researchers + ml_researchers + hci_researchers', departments)\n",
"fancier_res = fancier_model.fit()\n",
"fancier_res.summary()"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"fancier_model_predictions = fancier_res.predict()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I always struggle with: how do I interpret these parameters on logistic regression? This [link](https://stats.idre.ucla.edu/other/mult-pkg/faq/general/faq-how-do-i-interpret-odds-ratios-in-logistic-regression/) is helpful."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Intercept -0.830695\n",
"nlp_researchers 0.033579\n",
"ai_researchers 0.009386\n",
"ml_researchers -0.028682\n",
"hci_researchers 0.008449\n",
"dtype: float64"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fancier_res.params"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.03414924055928"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Holding everything else constant, adding one NLP researcher, increases the female:male ratio\n",
"np.exp(fancier_res.params['nlp_researchers'])"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0094304471581492"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Holding everything else constant, adding one AI researcher, increases the female:male ratio\n",
"np.exp(fancier_res.params['ai_researchers'])"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9717250353098399"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.exp(fancier_res.params['ml_researchers'])"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0084845088020287"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.exp(fancier_res.params['hci_researchers'])"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2deZhU1dG435phgAFkUIGfsoNBkH1fBBRxwQQFRRaXJJC4K9FPExSTKGBMQDGamGgUTWJMXEDigoof+gm4RmVwQSUqCCibAsoMDAzMVr8/eqGnp7vndk/f7ts99T5PPzN9+i517z331Dl16lSJqmIYhmEYXiMn3QIYhmEYRiRMQRmGYRiexBSUYRiG4UlMQRmGYRiexBSUYRiG4UlMQRmGYRiexFUFJSJnishnIrJBRGZF+H26iOwSkQ/8n0vclMcwDMPIHBq4dWARyQXuBU4HtgKrRWSpqq4L23SRqs5wSw7DMAwjM3FzBDUE2KCqG1W1DHgCmODi+QzDMIwswrURFNAW2BLyfSswNMJ254nIScDnwHWquiV8AxG5DLgMoGnTpgO7d+/ugriGkXo+2lYc9bfebQtSKIlhpI81a9bsVtVW4eVuKignPAc8rqqHRORy4B/AmPCNVHUhsBBg0KBBWlhYmFopDcMlRsxfwbai0hrlbVvk8+asGq+CYWQlIvJlpHI3TXzbgPYh39v5y4Ko6reqesj/9SFgoIvyGIbnmDm2G/l5udXK8vNymTm2W5okMgzv4KaCWg10FZHOItIQOB9YGrqBiBwb8nU88F8X5TEMz3FO/7bMm9ibti3yEXwjp3kTe3NO/7bpFs0w0o5rJj5VrRCRGcByIBf4m6p+IiK3AoWquhS4RkTGAxXAd8B0t+QxDK9yTv+2ppAMIwKSaek2bA7KMAwjuxCRNao6KLzcIkkYhmEYnsQUlGEYhuFJTEEZhmEYnsQUlGEYhuFJTEEZhmEYnsQUlGEYhuFJTEEZhmEYnsQUlGEYhuFJTEEZhmEYniTd0cyNesYz729jwfLP2F5USpsW+cwc283C/BiGERFTUEbKeOb9bdz01EeUllcCsK2olJue+gjAlJRhGDUwE5+RMhYs/yyonAKUlleyYPlnaZLIMAwvYwrKSBnbIyTmi1VuGEb9xhSUkTLatMiPq9wwjPqNKSgjZVj2WMMw4sGcJIyUEXCEMC8+w8gM0u11awrKSCmWPdYwMgMveN2aic8w0s3axXB3L5jTwvd37eJ0S2QYnvC6tRGUYaSTtYvhuWug3O/JWLzF9x2gz5T0yWXUe7zgdWsjKMNIJ6/celg5BSgv9ZUbRhrxgtetKSgju/G6+ax4a3zlhpEivOB1ayY+I3vJBPNZQTufXJHKDSONeMHr1hSUkb3EMp95RUGdekt1JQqQl+8rN4w0k26vW1NQRtaixVuROMrTQkBRvnKrz6xX0M6nnLyiQI2YpHudUDLw8jWYgjKylm9oyTHsilLuIfpMMYWUgXhhnVBd8fo1mJOEkbXMK5vMAW1YreyANmRe2eQ0SWRkE15YJ1RXvH4NpqCMrKWw+enMKr+ErVUtqVJha1VLZpVfQmHz09MtmpEFeGGdUF3x+jWYic/IWmaO7cZNT5WxtGxksCw/L5d5FpzWSAJtWuSzLUJDnknR+b1+DTaCMrKWc/q3Zd7E3rRtkY8AbVvkM29ib0/Y1o3MxwvrhOqK16/BRlBGVpNuN1kje/HCOqG64vVrEFVNtwxxMWjQIC0sLEy3GIZhGEaSEJE1qjoovNxMfIZhGIYnMRNfluDlxXaGYRiJYAoqC/D6YjvDMLKPVHSKzcSXBXh9sZ1hGNlFoFO8ragU5XCn+Jn3tyX1PKagsgCvL7YzDCO7SFWn2BRUFuCFxGKGYdQfUtUpdlVBiciZIvKZiGwQkVkxtjtPRFREargZGrXj9cV2hmFkF6nqFLvmJCEiucC9wOnAVmC1iCxV1XVh2x0BXAu845Ys2Y7XF9sZ2YEbk+LmfZqZ+MKIfVTNzOdGp9hNL74hwAZV3QggIk8AE4B1Ydv9BrgdmOmiLFmPRUww3MQNT9FM9T41pZq6TrGbCqotEJrLeiswNHQDERkAtFfVF0TEFJRheIm1i4OJFIfRktMrJ7OUw4F3A5PiiTZKsSbavdrgZ6pSdYNUdIrT5iQhIjnAXcDPHWx7mYgUikjhrl01E9AZhpFk1i72paIv3gIox7CL+XkPMT7njWqb1WVSPBO9T21JR2pxcwS1DWgf8r2dvyzAEUAvYJWIABwDLBWR8apaLdieqi4EFoIvFp+LMseNDfeNrOSVW6G8uqJoImXc0GBxtfQldZkU93qqh0h4WalmY1vk5ghqNdBVRDqLSEPgfGBp4EdVLVbVlqraSVU7AW8DNZSTl0nVYjXDSDnFWyMWt5Fvg//XdVI8E71PvbqkI1vbItcUlKpWADOA5cB/gcWq+omI3Coi4906byqx4b6RtRS0i1i8U1omLbdWJubr8qpSzda2yNVYfKq6DFgWVnZLlG1HuymLG3h5uG8YdeLUW3xzUKFmvrx8jjn7d2zqMy5pp8k079Nw77Vpzd7lhrxFNHn2a1jVznff+kxJuVzZ2hZZsNg6kIk2dCM9ZNz8QJ8prN68h/bvLaC17mantGRL75kMTkPj6zWCSnXtYnjuASj1twHFW3xKHVKupLK1LbJQR3XAq8N9w1tk4vzAM+9v48erOzLs4B/pcuhRhh38Iz9e3dHTMqecCI4klJf6ylNMtrZFpqDqQCba0I3Uk4nzA5koc8qJ4kgStdxFsrUtMhNfHck0G7qRejJxfiATZU45Be3868QilKeBbGyLbARlGC7jlmvyM+9vY8T8FXSe9QIj5q9IqvnNq+7UnuLUWyAv7H7k5fvKjaRgCipbWLsY7u4Fc1r4/q5dnG6JDD9uzA+4Pa/l2pxGNtXTPlPg7HugoD0gvr9n35MWL75sxUx82UAgLE15+r2JjJq4EVjT7Th2rgQDTWY9DYkTSEH63LvpM8XeMRcRVU9FDqqVQYMGaWFhxgSbSA1394piC28P132cenkM9/A3zFVFW9iuLbmjYgpLqw6HHhJg0/zkrVNKKsmqp+GKDnymNRu9ZCwiskZVa+QDNBNfNuAhbyLDRUICuOYItMvZXSOAq6fniJJVTz3k3m24iymobCCa11CavIkMl4gRwBUyYN1LsuqpdcjqDaagIuCmd5QrRPAmqshtzJz959XtGrJpQjsbiBHANd3rXhy9M8nyesvGDlkGvGvpaBfNSSKMjExIFrC7+yeND+Qfwy37z2NJ2RAgwWswxwvvEWXdTU6Ldrx53Zg0COTD8TsTVk8Tdm6IEicwY927M+BdS1e7aE4SYYyYvyJiTKu2LfJ5c1b6GoF4qO0aHMWFM8cL7+FR54C0vDMp9OJzPY5iBrxrbj/jaE4SNoIKIxtW0Me6Bsc9Ibfs/F5xD85EPBrANS3vTIrcu1MycsiAObV0tYs2BxVGNqygj3UNjmOsuWHnD0sjHjRleNDe7kWSFsA1yfMdaX9nXJy/SUlMwgyYU0vXM65VQYlInohcIyJL/J+fiUieq1KlkWyIChzrGhz3hCJNaOfkQdn+xBsCcw+uE0lpLF3oJKT1nXG505OSkUMGhExK1zN2MoL6CzAQuM//GeAvy0qyISpwrGtw3BMKD+OSfxSIQOl3JNwQZIApw8skpbF0oZOQ1nfG5U5PSkYOGRAyKV3PuFYnCRH5UFX71laWKrIlkkS6EtiF29TB1xOqtbIlYyI3AyaDvUxSJqrntAAivfMCc4rqJF9acPl6En5fjLioSySJShE5LuRAXYDKGNsbtZDOBHaJ9oQ0yignWnlEMsCU4WWSYmbJgPmOuHD5erLBopLJOPHimwmsFJGN+EJ9dQR+4qpUWY7bgT5rIzRvTGAkd92iD2KO5L6hJcewK0q5Q5K1DiaceuIZmJQArmlYQ+SqtSCO60lUjmzMs5Qp1KqgVPUVEekKBLppn6nqIXfFym62F5UyPucNbmiwmDayOxj087mikbXvnETicaGdVzaZeXkP0UTKgmUHtCHzyifzx3hOmmz34CQtckyXyTURueq07qQunYQEOgKuu2k7vJ6MXIBvRFdQIjJGVVeIyMSwn74nIqjqUy7LlrVMa/YuN5QfbuzbiS/o51F5DYHURaKOZyRX2Px0Zu3Fr1S/ZbsezR0VU1jT/PSUyRuRWJPk4Y1nlAbWq42Xa3Il0klIsCOwYPlnnF75Kjc0rN4ZW7C8Ye3X4FQhOriedFstjMSINYI6GVgBnB3hNwVMQSXIDXmLaFJRVq2siZRxQ94iYG7K5IjHK2zm2G7c9FQZS8sOj/Ly83KZl273e6eegTEa2AXLW3qy8fJUoxpPRyCEQXtfrjbyDnTGbtoLEGMkmOTwP9mwAL8+ElVBqeps/7+3quqm0N9EpLOrUqWBZJh4wo9xSvdWrPx0V41jNin9OuL+0cqTjr9n+kXjrWyvOrpGTqFILrSuJLAjCfc9Sny6GpPkMRrY7UW3Rzx0uhsvTzWqCS4RuKnhkzShZmfspoZPAvOi75igQoxGmxb5ET0gM2kBfn3EiZPEv/GtfQplCb61UVlBMkwpkY7xr7e/Cv5e7ZhOG1U3COmZ5nA4pxDlsLRqZEyvsGRPFifFhOV0kjxGA+vVxstTciVYZ/8fu+MqD5LkNXM+C0BNd/FMWoDvOVLgnBTVzVxEuovIeUCBiEwM+UwHGidVijSTjBX6kY4RTvCY6XS3jpFTKNUutEmJjOB0kWMMd2SvRg/xlFwJ1lmJct+jlQdJsvt4Ut3FXQitlHEpflIUtizWCKobcBbQgurzUPuAS5MqRZpJhinF6bbbikoZsawlf+g9l8Ff/Cli78Ox2SuRHkyUHmi7nG9THq09aSYsJ5P+MUZa5/Rxx3xZV9wyqyZEot5/ibq1J8sdPuQdOaegHef8oI69fBdSY3jVSScmSTbBRiPWHNSzwLMiMlxV/5O0M3qQZJhSoh0jEtuKSvnx6o7Mm7i8RgV0XFkTfVHSaV4MI6UmrFoaWFfWuiTBBJI0uZJhjknE+y9RxZaMNXNu5FlyoWFOuzNMEju6yQ5b5iTUUWPgYqAnIaY9Vf1pUiVxiBuhjpIRziTSMWojUogax+FsEg0b5KGcQrHuO7g/cnB17ZOH7rOnZEklboTWciG0UudZL0Q7Ipvmu7zsJNG6keR7W5dQR/8EjgHGAq8C7fCZ+bKGuOzTUezP5/RvyyODv+TtxteysdGFvN34Wu44/lPaxhgNRDJlOTZ7JdqD8VBgymj3HXA9FJTr4abSHLk9dE7j66d+WT+jyLvRy3chtFJa05VEq6dPXRp7fi1F8+hOvPi+p6qTRWSCqv5DRB4DXk+qFB7AkSkllskAGPzRbKAUBI5hF1N2LGDK2fcwYllLx6Ysp2avA/nH0KR0R43tDuQfQ5PYV5GyZG9OiHTfR8xf4brJIx6zSkIjrbo0jnU0x4WPTFvrLl93PIyqoq2Mmr/CE/NtruCGOduFUFFp9TCMVR/DTaLh9bLvhbD+pfR48YVQ7v9bJCK9gAKgdVKlyBRi9Ypj/BaPN5bTbe8on8oBbVit7IA25I7yqXFelPdIxfofp+dIeKSVaE87Cd5R4cp3u7aMuN12PTqlgYpTjhu9fBcsEGkNSFtbfQy0b5Hq5YeP+e7lnCKfWc+FTq+TEdRCETkS+DWwFGgG3Jx0STKBRHrFxVvj8sZyuu0/SobwXU5ZjdBDzx0awpyELq5uJHM+JxXOE07PkfAEdqI97SRMwocr2TsqpjA/QhzFOyqmOL+eTMSt4MTJsECEjUbOOfUWzpmVBqtGpHoaTvHWlHntheMkWOxD/n9fA7oAiEgH1yTyMrWZDGL8Fo83lpNt27TIZ2nRyGqhh4CYc15ukWw32VSYPJyeI+HRXKKNYxLmTcKV79KqkVAON+Yt5lgOd2ZCo4ekO2pGbSTcAfKQOTuIG96FcVD9XoYueYnQfoG/3UtPstGYCkpEhgNtgddUdaeI9AFmAaOA9q5K5kVq6xWnMI1B0hrxGPMdThuFZLvJpmL9j9Nz1Gk0l0jjmIR5k0h1Y2nVSJYeih4tP91RM2IRqQM088kPmfvcJxQdKPfMujVw+M6kaTQSkC/8XgaXvOS+Gb0Ne+XWtCxPiRXNfAG+hbofADeKyHLgEnwBtNLiYp52nPSKU5SXKCmNeIye3DOVIxyPiuo0ZxRFQaYiB4+Tc6R8AjsJk/DhdSNHhMoYy0m8EDUjFpE6QOVVyp4DvulxxyP2RJ1PHO7n2JKQptEI1NKZnFVL+5biPGIQYx2UiKwDBqjqQf8c1Bagl6pudlWiWkhLyvdsTYgXYy3DiEP3OE4vnnAq8gxZnxOPeWn10gdo/94CWusudkortgyYyeDxl9esQ13PiO4BFbpt/pG+stI9Cde9aOtswPeM/tBjfdSoJukkcN+dLoCPWd8SrWtx7Bf+HgRzvuV8S07ofU3SGqJEzJ51WnPlYjsYbR1ULBPfQVU9CKCqe0RkfbqVU1pIgb04bcnyYvTkth+MNxVHAqOMNJo6ahDj5XM6mlu99AF6rfk1+VIWXGpQsObXfPHtexy3/dnqdajwr4d3DK9TgU+S6l40M2XbFvm8+YPd8NzstM2HRCORhe8xR+yJ1rU49tseppyqOaaE3tckjJITnfdNucm6jsQaQRXhc4wIcFLod1UdX+vBRc4E/gjkAg+p6vyw368ArgYqgRLgMlVdF+uYdR1BzX3uE9Zt3+t8h62roSJCAuEGjaDd4ITlCLC75BCbdu+nKuQxtJRiOufuIqeyzHeeIztB01Z1PlcNYlzb+1VdKausqvFTw9wc+ndoUaN8d8khtnxXSlllFQ1zc2h/VD4tmzWKff7Nb0T/rVMKswvv3wXfroeqkOvNyYGju0a+7/t3wZ7NvnsX8nzKNr9Dw+CqjMMogkQdw4QQXqeSVPci1bEcgc4tm9KyaK2r9TtR3v+qqFr9O5pi2stOGlHOIfLYoq35loJq+0Srm0DidS2O/UJl7ifraRShLgTva5Q65JTw+xMg5j2glrpQ2/sahR5tmjP77J4J7RsgkRHUhLDvv4/zhLnAvcDpwFZgtYgsDVNAj6nq/f7txwN3AWfGcx7XifTyxiqPky3flVarLEdTTGd2kBOofBWHfI0nJF9JHdkpcsN8ZCfaa37Eitz+qMg9rZbNGsVfwRs0it44ppI9m6vfA/B937O55j0PV2YhzyeScgKcKafAsWJ9r608CoHnErEDsdvd+l0r0ZR9mHLqIjvIwVfWiHK6yA5QgkoqVt0EEq9rcezX/qjD70xE5QSHj9W0VdT3OVZnL/S3cI6mmPZVO2FzOVW5DfmqqjXfVDWvdoyYdSESdVSkdSVWsNhX63jsIcAGVd0IICJP4FN6QQWlqqFDmaZEDnKVVBxp+lBzT+Mc0AhmhoL2cPnMOsvTedYL1b4/2/Bm2uVEyJUj7eHyOGJcObUXR9zONzh23fS4dksM+/7w5J2nNuZ8H3KjWOYvD4uvdncvyIswfyDt+brRQY5hV42fKsihATUblBqE16m7L40+V5GEuufaOeKpe89dA7mlPhsLQGk+nHYPI/Ycjr4S7Z34mlYMP/hHZ3Uz0boW536Bd+b3B+6P/B7Xcl8DpruAAiqrrGJ70UF+NqYrQLXfQqlhUsS31m1W+SUsrRwZPEZCTlQRnk+qTH21BotN+MAik4AzVfUS//cfAUNVdUbYdlcD1wMNgTGquj7CsS4DLgPo0KHDwC+//NIVmYHIk6LhJHEiP3xidWOjC8mJFJZGhVH5TzlTEulwPnDZQ8pV4pm0jhEsdPWA2w/PQfkp1YZs7zSx+hxUJCI9n1Q8x2SfI57jxbjvz4xeHpxjifZOxB2gNZV1NMH7GsvhCIjqMPJGw2siKsStVS0ZWXZP8Bhvzhrj/HrcCLYbhURMfClBVe8F7hWRC/FFq5gWYZuFwELwzUG5KlCkSVEAyQWt4kD+MdxRPpV/PNaUNstWOPeAilIpwh0MtmtL2knNihYalgZqcadNtfNBXSbzvbCQsrZJ69BnJ9FG1O0YPP5yVoPfi283O6UlWwYGvPhOdezFd3jk2pRpzS7nhvxFNCn9Ora3X11SaEDyOgnx1L0YTjrVXOUPRH4n4l6Dk2hdS2GakUSXbLTJ+TZyuRwu315UGt+7mkZ3+ABujqCGA3NUdaz/+00AqjovyvY5wB5VLYj0ewDX3cxj9JCfmfBJNWUyPucNbs97qFqPuVovKdiAbMHnyKkRtws1pU1r9i6/1vtpUHkwuGlwqO5f+Z8rQpVqdNOGCykBYpLCnpZrRGvsUzyidpz6xasu+rXUvdC6/p/G10Y0idaoN1691joQzXyeyAiqbYt83mx0TcR3sMYIKsp2Ed9VD4ygag0WKyLHi8iDIvKSiKwIfBycczXQVUQ6i0hD4Hx8sfxCj9015Os4oIZ5L+XECPIZvsjthgaLqysniBJcEWq8tCHpDs7p35Y3Z41h0/xxzPn1XBpM+BMUtKdKha1VLaspJ4BK1diBS11ICRCTqD2tLUlLie2YRNNx95nie+nCA1/GGlG7kK4k1kLKaqQ5nUdUYtS98MC7vyubTGlYwOOI7tZupYhxIXW7E2IFII4VLDpmIOkIgXFD4y0Gt4tnVJSilBqxcGLiexK4H3gQnzu4I1S1QkRmAMvxTbH9TVU/EZFbgUJVXQrMEJHT8EVM30ME817KiWHu2f5Y9QahTSSzA0QPrhhpu0j4TQqjovSmQokYUsiFlAAxiRaeB3zlz1wFL95Yp8WmjnBjzVq0Z6RVroxGXc8HFoGkOMPUZi049RYWLPusRvglyuGXDZ/kGHbHrhtxmNkcXU+kuvLUZb48SAXtE6+jDsyusTohgYXGseSP/Ft1k2JgKuK5Q0NoG7rdqjhCabkVbDcOnGTUXaOqA1MkT62kJJJElEoWPvyONjFJQXt/Q1GL+bSWobLTxYoRV4Gn0vnAiRkslLqaZ6JdmxsmiRSbL13PqBxGpDoWUC9tnSqriM/ff5SQxr7WKAYJ1tlQhVSQn8f+sgrKKw+fKaKJNNr9C5BIHa3NFOm/vqqiLWzXljUC9no6g67L1CWj7nMicpWIHCsiRwU+LsjoHaKYe8KH2HdUTIluoqjNnOZgRBOeJyZXIroyRV4FHs1k5QbVTDAOqIspKla+JDcmdVNs5nCcOyxJckXqzQeadse5oiJaC/SwsvTXvZiZYxPMgxVuLisqLa+mnCCKibS2OpFIHY1ldg25vhyBdjm7mZ/3EONzDi8ETknAXg9l1HaCEwU1DZgJvAWs8X9SHAwvvQTSZ1+36AMaNcjhyCZ5CLCm+el8PPC2yA87UgMSSGsaR6UInZ/6/ZS+jhMfppyAQnSqpBJVGrEaATfm3lL8QjtOXhdJrr4X+u5DHHMqtXmHRWzcw3HYMYipfBOcU4ukYCNR4zqd1Il462is+xDh+pqIL58bhL3HsebGkjFvlsrOax1xkg+qcyoE8SrhJpCi0nLy83K5e2o/f6MxBri85o4u2G9TkYaizjhJgAaJK41YjcDEhfHNvTk1KaXYFd5xJPdQuRKcf4sWmy2UWqPSO0wRErP+PpvY6NdpHqsaoxMH9fRrWjJ81gvO37NY9yHKdbSRb6ubUmM9R3D+jL2wvjAJ1KqgRKQJvoW0HVT1Mr/nXTdVfd516VJJlAdap1xHLjRsqUhDUSfCFXP+kVBWApVhrvi1maKivWCxGoF4OgVpThqXdBJc+xYp0G84tZqe4nDKiVp/E8yD5UTBRrQyVKsrNR07yjSXhpTyRaML2X6gJX94+nzgqsSzKEfJp5TToh1vXhcyt1jbSNLJM86iuu3ExPd3oAw40f99G3CbaxKlgxj27zrlOqqvhJoQbtwEE+6Nz0QWaz6itrkXp+aLZLlpp8lVuQYJzr+FmhQhaIQO4siEnAwzaIJzapHMhuc2eJO3Gl/DxkYX8nbja3lk8JeRFUuwrhT7Rt9++Ys4AhCOkpLgfNGtspAPXlgY+xpi3Qen1xfrOTp9xl5dgpAATrz4ClV1kIi8r6r9/WUfqmrflEgYRp28+BLw/oonL5KRJGrzUEuG+SIZi5lr8YhKaRqVNOYYShpJ8OKLtNA9Xi+1rbccFzVsULtbv3B8OTVwcn2xniM4e8ZR6zb4FKf3TH7RvPicKKi3gFOBN1V1gIgcBzyuqkPcETU2CSuoSI1JTh40OgJKv4uyU83oERDFbdVIHqmIhJGMBt1hLLkArtabZLkPZ/rcRRKea9WcFuREqH9VCDluRGIJJdZzBEfu/LW60Ice0yPPti5u5rOB/wXai8ijwCvADUmWz30iDXurymMoJ6CgnXOvKiN5pCISRjLctGOYXBxHhEgWyTCzxePq7RXTZjhJWGpwMP+YuMqTSqznWGM5R8i8WW1m8HAyxOTnxIvvZRF5DxiG745cq6pRQih4mHhdRkMaK9cdE9zutWZarzgVkTCS4WUZY2J/+zdpmLusq1OOU0cLL0/CJ+hsEUqT799KxbM/q2YmrMhtTJPvp6hBj/UcA79FGiUFnlVgpBio29HMfSkM+pooURWUiAwIK9rh/9tBRDqo6nvuieUCscLx1Ni2DqFO4sXtl93LjUk0UhViJcEGPTDnMWjv2cxv+FfyCUlo51ekbZbVIbV2ukjGJHy661QyOjd9pvgaxpD618BrnbranlVo3Y5q9nQpNmcSiTWCKgQ+BgKjpVAHH8W3AChzcLw+J8URuN1+2b3cmMQixWuPnBK6Lm4bI9EyuDFvMW3kWyQ0jUpl5KjknlhUHQ2now8PpGGISrI6Nx6tf0HiGSmmOjZnEomloK4HJgGlwBPA06pakhKp3CBZ63OSjdsvu5cbkwwkfG5padVIlh4a6fPqDFnPkvZF1YmYdZ02ZEkwo7mK15VLMohH6Xgg6GuiOPHi64IvVcYE4Evgd6r6QQpki0hSg8V6YW7GZZfqA7d3p0npjprl+cfS5MZPExZ79dIH/In5drFTWrFlgD8xX13wwvOohWgBTyGO4Koucdj0+HJk06MTpwknz8CjAUfrHRnwvjglYTdz/8498SmpHwE3qGraXHZSEs08ldIrtQ0AACAASURBVMTrVhpnQzDnttncUH4fTULyVh3QhtyRdxVzfj03IZFXL30gYmrzjwfelriSypBGL1q08QDpWoIQanqMGWU/WebrLGocjfQTt5u5iHQRkV+KyDvAXOBD4IR0KqesJJZbaRJWhP+jZAizyi9ha1XLagkQ/1GS+DK29u8tqJGoMV/KaP/egoSPWeu1esStOVLkglBcdSWPQajpMWaesmSRQQFHDecEAmN3nvUCI+avqD2SvcvEmoPaAKwFngX2Ah2AK8Wf8kFV73JduvpCNJt5EuaP2rTIZ2nRSJaWjaxW3rYO3mStdVfNmDhA67qsPoh1rR7yRAydW4o2kkpHGKzQc27XlrSLpKS8MkdkxEeKRqvhgbED6VaAtJmtYy3UvRV4GqgCmgFHhH0Mt0nCglXH+YXiYKe0ilLeMuFjxrzWdMQWizFiC6RAiabknbiSJ7unGnrOOyqmcMBJKnWX8FovPKNJME9WIqR8cbkDoiooVZ2jqnOjfVIpZL0l0orwnDwo2+/Y1OVGJIwtA2bWSNRYqg3ZMmBmwseMFdlBo4yuopXXGYeNQqLKPzzJnuPEgDEIlWVp1UhmlV/CNm2JpjgpXV2uLdWKLSMUaQo7Z14MjF1rJAkjjURzjQ+EZ3Jo6kp2JIzB4y9nNfi9+HazU1qyZWAdvfhiuMJ+89QvOYZdNXb5hpa4EnzG4dqxRF3J65TCJQrhsqxpfjqrx85IuWkm0WtLtXnJi+asiKRwmUi01CXpXFxuCsrrhK8ID48dmKZFt4PHXw5+hXSM/xMg4YjYUebi5pVNZl7eQzU8EeeVT+aPdbyOiMTRKCSi/N3qqXohV1ii1+aG0vbS+RImhWvOIuUGS/ficifBYg2vkAGLbt0wXxU2Pz2iJ2Jh89Mjnr/OZhuXg9VG65F6OgySQxK9tlSbl7xozopIMoIaO8SLgbGdZNQtAOYAo/xFrwK3qmqxi3IZkUhFb6qOHkNu9Ex9Pbuyap6I+Xm5zAvr2SXNbONyaBgv9lTjJdooOdFrS7V5Kdr5FN9at3QuuK5GiqNAeGEUHooTE9/f8MXkC9yRH+HLsjvRLaGMKLgdUysJ7txu9EydzvUkTTm63CikPQxSHXHSEYj32lKttGOlundrPirZpu/6gJNQRx+oar/aylJF1kWSiBc310QkIdlbtEgLbmYgDrz40dYlCbBp/jhXzl0fcesZpzqbb231Jpl1NlypgyU+DSVaJAknI6hSERmpqm/4DzQCXwBZIx242ZtKwhxXqnvCkV78cLJhbsdLZIuTR+B80eIrJnM+KmOcMjyGEwV1BfCIfy4KYA8wzT2RjLSRhDmuVJuvIr34oWTa3E4m4EV35LqQiuvJGKcMj+FEQe1V1b4i0hxAVfeKSGeX5TLSQZLmuFLZE471gqc7uni2kg1OHqGk4nqyTamnCidu5v8Gn2JS1b3+siXuiWSkjViBaz1KtBc8MH9gyin5eNEduS6k4nrcCDlWH4jqJCEi3YGewB1AaAyb5sBMVe3pvng1qfdOEllMIpPkNvlsZAqpdgLJJBJxkugGnAW0AM4OKd8HXJpc8Yz6TqJrmDLdZduoP3htjVEm4MTNfLiq/idF8tSKjaCyk3S4pxuG4Q0SdjP3knIyshfzcjLcwMxqmY0FizU8gXk5Gckm3RHLTTnWHQsWa3gC83Iykk06E/C5ETS5PuIkWGwj4DygU+j2qupiOlOjvuE1Zwfr/WY+6TQbW+SI5ODExPcsUAysAQ65K45Rn/GKl1O6TUNGckin2djmVJODEwXVTlXPdF0Sw/AI1vvNDtIZ8cLmVJODkzmot0SkdyIHF5EzReQzEdkgIrMi/H69iKwTkbUi8oqIdEzkPE5JSjI7I+ux3m92kM6IFzanmhycjKBGAtNFZBM+E58Aqqp9Yu0kIrnAvcDpwFZgtYgsVdV1IZu9DwxS1QMiciW+qBVTE7iOWjGzjeEU6/1mD+kyG3ttTjVTcaKgvp/gsYcAG1R1I4CIPAFMAIIKSlVXhmz/NvDDBM9VK2a2MZySbcFQjfTglTnVTMbJQt0vRWQk0FVV/y4irYBmDo7dFgjN3bAVGBpj+4uBFyP9ICKXAZcBdOjQwcGpa2JmG8OpZ571fg3DGzhxM58NDMIXm+/vQB7wL2BEsoQQkR/6z3FypN9VdSGwEHyhjhI5h5lt6jfxmnit92sY6ceJk8S5wHhgP4CqbgeOcLDfNqB9yPd2/rJqiMhpwK+A8arqmhu7TVrWb9K5aNMw6kJ9du5yMgdVpqoqIgogIk0dHns10NWf3HAbcD5wYegGItIfeAA4U1V3Ohc7fsxsU78xE6+RidR35y4nCmqxiDwAtBCRS4GfAg/WtpOqVojIDGA5kAv8TVU/EZFbgUJVXQoswDef9aSIAHylquMTvJZaMbNN/cVMvEYmUt+du5w4SdwpIqcDe/HNQ92iqi87ObiqLgOWhZXdEvL/afGJaxiJYZ55RiZS30f+jqKZq+rLIvJOYHsROUpVv3NVMiMtZGsMOjPxGplIfR/5O/HiuxyYCxwEqvAv1AW6uCuakWqy3d5tJl4j06jvI38nI6hfAL1Udbfbwhjppb7buw3Da9T3kb8TBfUFcMBtQdJBtpqzEqW+27sNw4vU55G/EwV1E76Ase8Qkm5DVa9xTaoUkO3mrESo7/ZuwzC8hZOFug8AK/DFylsT8slobOFmTWwxs2EYXsLJCCpPVa93XZIUY+asmtR3e7dhGN7CiYJ60R+s9Tmqm/gy2s3czFmRqc/2bsPIVjJ1vt2Jie8C/PNQHDbvFbopVCowc5ZhGPWBwHz7tqJSlMPz7ZkQ089JJInOqRAk1Zg5yzCM+kAmLx9xslA3D7gSOMlftAp4QFXLXZQrJZg5yzCMbCeT59udmPj+AgwE7vN/BvrLDMMwDI8TbV49E+bbnSiowao6TVVX+D8/AQa7LZhhGIZRdzJ5vt2JF1+liBynql8AiEgXoLKWfQzDNTLVI8kw0kEmz7c7UVAzgZUishFfoNiOwE9clcowomARQAwjfjJ1vj2qghKR/8HnWv4q0BVfLiiAz9xMzW4YschkjyTDMOIj1giqHfAHoDvwEfAmPoW1lZAFu4ZRG8k0yWWyR5JhGPERVUGp6i8ARKQhMAg4EZ9pb6GIFKlqj9SIaGQyyTbJWQQQw6g/OPHiyweaAwX+z3bgHTeFMrKHZAflzWSPJMMw4iPWHNRCoCewD59Cegu4S1X3pEg2IwtItkkukz2SDMOIj1hzUB2ARsB6YBu+uaeiVAhlZA9umOQy1SPJMIz4iGriU9Uz8S3IvdNf9HNgtYi8JCJzUyGckfmYSc4wjESJuQ5KVRX4WESKgGL/5yxgCDDbffGMTMdMcoZhJEqsOahr8HnunQiU45uDegv4Gz63c8NwhJnkDMNIhFgjqE7Ak8B1qrojNeIYhmEYho9Y66CyLs27YRiGkTk4WQdlGIZhGCnHFJRhGIbhSUxBGYZhGJ7EFJRhGIbhSUxBGYZhGJ7EFJRhGIbhSUxBGYZhGJ7EFJRhGIbhSUxBGYZhGJ7EFJRhGIbhSUxBGYZhGJ7EVQUlImeKyGciskFEZkX4/SQReU9EKkRkkpuyGIZhGJmFawpKRHKBe4HvAz2AC0SkR9hmXwHTgcfcksMwDMPITNwcQQ0BNqjqRlUtA54AJoRuoKqbVXUtUOWiHClh69atTJgwga5du3Lcccdx7bXXUlZW5vp5mzVrBsDmzZvp1auX6+czDMNIFW4qqLbAlpDvW/1lWYeqMnHiRM455xzWr1/P559/TklJCb/61a/qfOyKiookSJhcVJWqqozvUxiG4XFipnz3CiJyGXAZQIcOHWJuO/e5T1i3fW9Sz9+jTXNmn90z6u8rVqygcePG/OQnPwEgNzeXu+++m86dOzN37lzGjBnDX//6V3r29B1j9OjR3HnnnZxwwgn87Gc/4+OPP6a8vJw5c+YwYcIEHn74YZ566ilKSkqorKzkhRdeYMKECezZs4fy8nJuu+02JkyYEFWeUEpKSiLuO2vWLNq3b8/VV18NwJw5c2jWrBm/+MUvWLBgAYsXL+bQoUOce+65zJ07l82bNzN27FiGDh3KmjVrWLZsGfPnz2f16tWUlpYyadIk5s6dC8CyZcu4/vrradq0KSNGjGDjxo08//zz7N+/P+L1GoZhRMLNEdQ2oH3I93b+srhR1YWqOkhVB7Vq1SopwiWTTz75hIEDB1Yra968OR06dGDDhg1MnTqVxYsXA7Bjxw527NjBoEGD+O1vf8uYMWN49913WblyJTNnzmT//v0AvPfeeyxZsoRXX32Vxo0b8/TTT/Pee++xcuVKfv7zn6OqjmSLtm+oTACLFy9m6tSpvPTSS6xfv553332XDz74gDVr1vDaa68BsH79eq666io++eQTOnbsyG9/+1sKCwtZu3Ytr776KmvXruXgwYNcfvnlvPjii6xZs4Zdu3YFzxHreg3DMMJxcwS1GugqIp3xKabzgQtdPB9AzJFOupgyZQpnnHEGc+fOZfHixUya5HNYfOmll1i6dCl33nknAAcPHuSrr74C4PTTT+eoo44CfCa1X/7yl7z22mvk5OSwbds2vvnmG4455phazx1t3/79+7Nz5062b9/Orl27OPLII2nfvj1//OMfeemll+jfvz/gG4GtX7+eDh060LFjR4YNGxY89uLFi1m4cCEVFRXs2LGDdevWUVVVRZcuXejcuTMAF1xwAQsXLox5vSeccEIybrNhGFmGawpKVStEZAawHMgF/qaqn4jIrUChqi4VkcHA08CRwNkiMldVvadhaqFHjx4sWbKkWtnevXv56quv+N73vkeTJk04+uijWbt2LYsWLeL+++8HfMrj3//+N926dau27zvvvEPTpk2D3x999FF27drFmjVryMvLo1OnThw8eNCRbLH2nTx5MkuWLOHrr79m6tSpQZluuukmLr/88mrH2bx5czWZNm3axJ133snq1as58sgjmT59eq0yRbtewzCMSLi6DkpVl6nq8ap6nKr+1l92i6ou9f+/WlXbqWpTVT06E5UTwKmnnsqBAwd45JFHAKisrOTnP/8506dPp0mTJgBMnTqVO+64g+LiYvr06QPA2LFj+dOf/hQ0173//vsRj19cXEzr1q3Jy8tj5cqVfPnll45li7Xv1KlTeeKJJ1iyZAmTJ08OyvS3v/2NkpISALZt28bOnTtrHHfv3r00bdqUgoICvvnmG1588UUAunXrxsaNG9m8eTMAixYtCu7j9HoNwzDAIkkkBRHh6aef5sknn6Rr164cf/zxNG7cmN/97nfBbSZNmsQTTzzBlClTgmU333wz5eXl9OnTh549e3LzzTdHPP5FF11EYWEhvXv35pFHHqF79+6OZYu1b8+ePdm3bx9t27bl2GOPBeCMM87gwgsvZPjw4fTu3ZtJkyaxb9++Gsft27cv/fv3p3v37lx44YWMGDECgPz8fO677z7OPPNMBg4cyBFHHEFBQUFc12sYhgEgTifbvcKgQYO0sLAw3WIYMSgpKaFZs2aoKldffTVdu3bluuuuS7dYhmF4FBFZo6qDwsttBGUknQcffJB+/frRs2dPiouLa8xnGYZhOMFGUIZhGEZasRGUYRiGkVGYgjIMwzA8iSkowzAMw5OYgjIMwzA8iSmoJJGbmxv0XOvbty+///3vUxLx++GHH2b79u1Zcx7DMIwAGRHNPNk88/42Fiz/jO1FpbRpkc/Msd04p3/dMoHk5+fzwQcfALBz504uvPBC9u7dG4zw7QaVlZU8/PDD9OrVizZt2rh2HiBl5wlQUVFBgwb1snoahuGn3o2gnnl/Gzc99RHbikpRYFtRKTc99RHPvJ9QoPWItG7dmoULF/LnP/8ZVaWyspKZM2cyePBg+vTpwwMPPADAqlWrOOmkkxg3bhzdunXjiiuuCI66rrzySgYNGkTPnj2ZPXt28NidOnXixhtvZMCAATz++OMUFhZy0UUX0a9fP0pLS+nUqRM33XQT/fr1Y9CgQbz33nuMHTuW4447LhgDEGDBggVBeQLH37x5MyeccAKXXnopPXv25IwzzqC0tJQlS5bUOE8oDz74IIMHD6Zv376cd955HDhwgOLiYjp27Bi8nv3799O+fXvKy8v54osvgpEmRo0axaeffgrA9OnTueKKKxg6dCg33HAD7777LsOHD6d///6ceOKJfPbZZwAcOHCAKVOm0KNHD84991yGDh1KYOnBSy+9xPDhwxkwYACTJ08OhmyaNWsWPXr0oE+fPvziF79I2rM2DMNFVDWjPgMHDtS6cOK8V7Tjjc/X+Jw475U6Hbdp06Y1ygoKCvTrr7/WBx54QH/zm9+oqurBgwd14MCBunHjRl25cqU2atRIv/jiC62oqNDTTjtNn3zySVVV/fbbb1VVtaKiQk8++WT98MMPVVW1Y8eOevvttwfPcfLJJ+vq1auD3zt27Kj33Xefqqr+z//8j/bu3Vv37t2rO3fu1NatW6uq6vLly/XSSy/Vqqoqrays1HHjxumrr76qmzZt0tzcXH3//fdVVXXy5Mn6z3/+M+J5Qtm9e3fw/1/96ld6zz33qKrq+PHjdcWKFaqq+sQTT+jFF1+sqqpjxozRzz//XFVV3377bT3llFNUVXXatGk6btw4raioUFXV4uJiLS8vV1XVl19+WSdOnKiqqgsWLNDLLrtMVVU/+ugjzc3N1dWrV+uuXbt01KhRWlJSoqqq8+fP17lz5+ru3bv1+OOP16qqKlVV3bNnT8TrMAwjPeALIF6jva93NpTtRaVxlSeDl156ibVr1wYjnhcXF7N+/XoaNmzIkCFD6NKlC+BLTfHGG28wadKkiKksAkFmA5HHozF+/HgAevfuTUlJCUcccQRHHHEEjRo1oqioiJdeeilqSo3OnTvTr18/AAYOHBgM+hqLjz/+mF//+tcUFRVRUlLC2LFjg3IuWrSIU045hSeeeIKrrrqKkpIS3nrrrWBwWoBDhw4F/588eTK5ubnB+zRt2jTWr1+PiFBeXg7AG2+8wbXXXgtAr169gvfl7bffZt26dcG4gGVlZQwfPpyCggIaN27MxRdfzFlnncVZZ51V6zUZhpF+6p2CatMin20RlFGbFvlJPc/GjRvJzc2ldevWqCp/+tOfgg13gFWrViEi1cpEpNZUFqFpLyLRqFEjAHJycoL/B75XVFTETKkRun1ubm4Nc14kpk+fzjPPPEPfvn15+OGHWbVqFeBTlL/85S/57rvvWLNmDWPGjGH//v20aNEiOF8XTui13XzzzZxyyik8/fTTbN68mdGjR8eUQ1U5/fTTefzxx2v89u677/LKK6+wZMkS/vznP7NixYpar8swjPRS7+agZo7tRn5ebrWy/LxcZo5NXo6iXbt2ccUVVzBjxgxEhLFjx/KXv/wlOAL4/PPPg5lk3333XTZt2kRVVRWLFi1i5MiRUVNZROKII46IGG08Fk5Tajg9z759+zj22GMpLy/n0UcfDZY3a9aMwYMHc+2113LWWWeRm5tL8+bN6dy5M08++STgUyoffvhhxOMWFxfTtq3PeeXhhx8Olo8YMSKYDXjdunV89NFHAAwbNow333yTDRs2AL55r88//5ySkhKKi4v5wQ9+wN133x31fIZheIt6N4IKeOsl24uvtLSUfv36UV5eToMGDfjRj37E9ddfD8All1zC5s2bGTBgAKpKq1ateOaZZwAYPHgwM2bMYMOGDZxyyimce+655OTkBFNZtG/fPmiyikTAsSA/P5///Oc/jmQ944wz+O9//8vw4cMBnyL517/+FTStOTlPfv7hEedvfvMbhg4dSqtWrRg6dGg1RTZ16lQmT54cHFWBL4nilVdeyW233UZ5eTnnn38+ffv2rXHOG264gWnTpnHbbbcxbty4YPlVV13FtGnT6NGjB927d6dnz54UFBTQqlUrHn74YS644IKg2fC2227jiCOOYMKECRw8eBBV5a677nJ0nwzDSC8WLDaNrFq1ijvvvJPnn38+3aJkFJWVlZSXl9O4cWO++OILTjvtND777DMaNmyYbtEMw0iAaMFi690Iysh8Dhw4wCmnnEJ5eTmqyn333WfKyTCyEBtBGYZhGGnF0m0YhmEYGYUpKMMwDMOTmIIyDMMwPIkpKMMwDMOTmIJKEoF0G7169WLy5MkcOHAg4WNNnz49GBbpkksuYd26dVG3XbVqFW+99Vbw+/33388jjzyS8LlDj5vKkECdOnVi9+7dST3mnDlzuPPOOyP+duKJJ0bdL9Zv8VBYWMg111wT1z6/+93v6nTO2uqLYWQS9VNBrV0Md/eCOS18f9curvMhA+k2Pv74Yxo2bFgtcjj40kckwkMPPUSPHj2i/h6uoK644gp+/OMfJ3Su+kToPQsQeEaRfkuEQYMGcc8998S1T10VVG31xTAyifqnoNYuhueugeItgPr+PndNUpRUgFGjRrFhwwZWrVrFqFGjGD9+PD169IiadkNVmTFjBt26deO0006rFnZo9OjRwVQS//u//8uAAQPo27cvp556Kps3b+b+++/n7rvvpl+/frz++uvVRg0ffPABw4YNo0+fPpx77rns2bMneMwbb7yRIUOGcPzxx/P6669HvI69e/fGlQokUkqLXbt2cd555zF48GAGDx7Mm2++CcC3337LGWecQc+ePbnkkkuIttzh8ccfp3fv3vTq1Ysbb7wxWN6sWTN+9atf0bdvX4YNG8Y333wTcf9169YxevRounTpUk1ZNGvWDKDGMwr/bfTo0UyaNInu3btz0UUXBeVctmwZ3bt3Z+DAgVxzzTURR5uho9A5c+bw05/+NKIsofcvEJHkoosuAuCuu+6iV69e9OrViz/84Q+AL2ZiQJ4TTjiBSZMmBUfsseoLwKuvvkq/fv3o168f/fv3jztMlmGklEghzr38qWu6Db2rp+rs5jU/d/Ws02ED6TbKy8t1/Pjxet999+nKlSu1SZMmunHjRlXVqGk3/v3vf+tpp52mFRUVum3bNi0oKAim3Qikudi5c6e2a9cueKxAOo7Zs2frggULgnKEfu/du7euWrVKVVVvvvlmvfbaa4PHvP7661VV9YUXXtBTTz21xvXEmwokWkqLCy64QF9//XVVVf3yyy+1e/fuqqr6s5/9TOfOnauqqs8//7wCumvXrmoybNu2Tdu3b687d+7U8vJyPeWUU/Tpp59WVVVAly5dqqqqM2fODN7XUGbPnq3Dhw/XgwcP6q5du/Soo47SsrKyas8r/BmF/9a8eXPdsmWLVlZW6rBhw/T111/X0tLSas/i/PPP13HjxkW8h4HyWLKEEpq2pbCwUHv16qUlJSW6b98+7dGjh7733nu6adMmBfSNN95QVdWf/OQnwWdeW30566yzgvvt27cvmM7EMNIJUdJt1L8RVPHW+ModEuj5Dho0iA4dOnDxxRcDMGTIEDp37gz40m488sgj9OvXj6FDh/Ltt9+yfv16XnvtNS644AJyc3Np06YNY8aMqXH8t99+m5NOOil4rKOOOiqmPMXFxRQVFXHyyScDMG3aNF577bXg7xMnTgRip9QIpALJzc0NpgIBWLx4MQMGDKB///588sknrFu3rlpKi6eeeoomTZoA8H//93/MmDGDfv36MX78ePbu3UtJSQmvvfYaP/zhDwEYN24cRx55ZI3zr169mtGjR9OqVSsaNGjARRddFLyGhg0bBkcnsa5h3LhxNGrUiJYtW9K6deuII63QZxTpt3bt2pGTk0O/fv3YvHkzn376KV26dAnuc8EFF0TcNxFZQnnjjTc499xzadq0Kc2aNWPixInB0W5ojMYf/vCHwWcTIFp9GTFiBNdffz333HMPRUVFlrXY8DT1r3YWtPOb9yKU14HQlO+hhKaP0ChpN5YtW1ancydCIK1Gbm5u1PmxeFKBNGjQIGJKi6qqKt5++20aN26cVPnz8vKC8sW6hvD0IZG2i5W+xMn+TknmsSI9GyfMmjWLcePGsWzZMkaMGMHy5cvp3r17wnIYhpvUvxHUqbdAXljup7x8X7nLREu7cdJJJ7Fo0SIqKyvZsWMHK1eurLHvsGHDeO2119i0aRMA3333HRA9DUZBQQFHHnlksMf9z3/+Mziacko8qUCipbQ444wz+NOf/hQ8ZkCJn3TSSTz22GMAvPjii8H5sVCGDBnCq6++yu7du6msrOTxxx+P+xrcoFu3bmzcuDE4alu0aFHSjp2XlxesH6NGjeKZZ57hwIED7N+/n6effppRo0YB8NVXXwWj1z/22GOMHDmy2nGi1ZcvvviC3r17c+ONNzJ48GA+/fTTpMluGMmm/o2g+kzx/X3lVp9Zr6CdTzkFyl0kWtqNc889lxUrVtCjRw86dOgQTIMRSqtWrVi4cCETJ06kqqqK1q1b8/LLL3P22WczadIknn322WqKAOAf//gHV1xxBQcOHKBLly78/e9/j0veeFKB7Nu3L2JKi3vuuYerr76aPn36UFFRwUknncT999/P7NmzueCCC+jZsycnnngiHTp0qHH+Y489lvnz53PKKaegqowbN44JEybEdQ1ukJ+fz3333ceZZ55J06ZNGTx4cNKOfdlll9GnTx8GDBjAo48+yvTp0xkyZAjgqz/9+/dn8+bNdOvWjXvvvZef/vSn9OjRgyuvvLLacaLVlz/84Q+sXLmSnJwcevbsyfe///2kyW4YycaCxRpGApSUlNCsWTNUlauvvpquXbty3XXXpeTcmzdv5qyzzuLjjz9OyfkMw20sWKxhJJEHH3yQfv360bNnT4qLi7n88svTLZJhZB02gjIMwzDSio2gDMMwjIzCFJRhGIbhSUxBGYZhGJ7EFJRhGIbhSVxVUCJypoh8JiIbRGRWhN8bicgi/+/viEgnN+UxDMMwMgfXFJSI5AL3At8HegAXiEh4HoCLgT2q+j3gbuB2PzokBQAAB4xJREFUt+QxDMMwMgs3R1BDgA2qulFVy4AngPAwABOAf/j/XwKcKk6DihmGYRhZjZuhjtoCoVFZtwJDo22jqhUiUgwcDVRLrSoilwGX+b+WiMhnCcjTMvy4ht2TCNg9qYndk8jYfalJovekY6TCjIjFp6oLgYV1OYaIFEZaCFafsXtSE7snNbF7Ehm7LzVJ9j1x08S3DWgf8r2dvyziNiLSACgAvnVRJsMwDCNDcFNBrQa6ikhnEWkInA8sDdtmKTDN//8kYIVmWuwlwzAMwxVcM/H555RmAMuBXOBvqvqJiNyKL73vUuCvwD9FZAPwHT4l5hZ1MhFmKXZPamL3pCZ2TyJj96UmSb0nGRcs1jAMw6gfWCQJwzAMw5OYgjIMwzA8SdYrqNrCLdUHRKS9iKwUkXUi8omIXOsvP0pEXhaR9f6/R6Zb1lQjIrki8r6IPO//3tkfdmuDPwxXw3TLmGpEpIWILBGRT0XkvyIyvL7XFRG5zv/ufCwij4tI4/pYV0TkbyKyU0Q+DimLWDfExz3++7NWRAbEe76sVlAOwy3VByqAn6tqD2AYcLX/PswCXlHVrsAr/u/1jWuB/4Z8vx242x9+aw++cFz1jT8C/6uq3YG++O5Pva0rItIWuAYYpKq98Dl9nU/9rCsPA2eGlUWrG98Huvo/lwF/ifdkWa2gcBZuKetR1R2q+p7//334Gpy2VA819Q/gnPRImB5EpB0wDnjI/12AMfjCbkH9vCcFwEn4PGxR1TJVLaKe1xV8Hs/5/vWaTYAd1MO6oqqv4fO4DiVa3ZgAPKI+3gZaiMix8Zwv2xVUpHBLbdMkiyfwR4zvD7wD/D9V3eH/6Wvg/6VJrHTxB+AGoMr//WigSFUr/N/rY33pDOwC/u43fT4kIk2px3VFVbcBdwJf4VNMxcAarK4EiFY36tz+ZruCMkIQkWbAv4H/UdW9ob/5F0jXmzUHInIWsFNV16RbFo/RABgA/EVV+wP7CTPn1cO6ciS+0UBnoA3QlJpmLoPk141sV1BOwi3VC0QkD59yelRVn/IXfxMYcvv/7kyXfGlgBDBeRDbjM/2OwTf30sJvxoH6WV+2AltV9R3/9yX4FFZ9riunAZtUdZeqlgNP4as/9b2uBIhWN+rc/ma7gnISbinr8c+t/BX4r6reFfJTaKipacCzqZYtXajqTaraTlU74asXK1T1ImAlvrBbUM/uCYCqfg1sEZFu/qJTgXXU47qCz7Q3TESa+N+lwD2p13UlhGh1YynwY7833zCgOMQU6IisjyQhIj/AN9cQCLf02zSLlHJEZCTwOvARh+dbfolvHmox0AH4EpiiquEToFmPiIwGfqGqZ4lIF3wjqqOA94EfquqhdMqXakSkHz7HkYbARuAn+Dqz9bauiMhcYCo+j9j3gUvwzafUq7oiIo8Do/Gl1fgGmA08Q4S64Vfmf8ZnDj0A/ERVC+M6X7YrKMMwDCMzyXYTn2EYhpGhmIIyDMMwPIkpKMMwDMOTmIIyDMMwPIkpKMMwDMOTmIIyDD8iUikiH/ijVn8oIj8XEdffERGZLiJt0nkeEXlYRCb5/1/lzwCw1h/R/M8i0sJt+QwjHFNQhnGYUlXtp6o9gdPxRWOe7eYJ/RH3p+MLoeM28ZznIlXtA/QBDlF/F6EaacQUlGFEQFV34ksRMMO/Ej5XRBaIyGr/yOJy8C3yFZHXROQF/6jj/sCoS0T+IiKF/hHZ3MCxRWSziNwuIu8BFwCDgEf9o7d8/+/z/N8LRWSAiCwXkS9E5IqQ48wMkWeuv6yT+HI4Peg/70v+Y04KP4/D+1CGL6BuBxHpm5SbaxgOMQVlGFFQ1Y34IpC0xpfrp1hVBwODgUtFpLN/0yHAz/DlHDsOmOgv/5WqDsI3CjlZRPqEHP5bVR2gqv8CCvGNWPqpaqn/969UtR++CCAP4wupMwwIKKIz8OXZGQL0AwaKyEn+fbsC9/pHgkXAeaq6JMp5nNyHSuBDoLvTfQwjGTSofRPDMIAzgD6BeRqgAJ8iKAPe9SuzQCiYkfiCrE4RkcvwvWfH4lNga/37L6rlfIGYkR8Bzfx5vPaJyCH/fNAZ/s/7/u2a+eX5Cl9g0w/85WuATgldcXUkCccwjLgwBWUYUfDH5avEF51ZgJ+p6vKwbUZTM72A+kdXvwAGq+oeEXkYaByyzf5aTh+I6VYV8n/gewO/PPNU9YEweTqFbV8JODLnRcM/T9ab6pmHDcN1zMRnGBEQkVbA/cCf/TlulgNX+tOWICLH+xP5AQzxR8zPwRdQ9A2gOT4lVCwi/w+fw0U09gFHxCnicuCn/hxfiEhbEWldyz5xn8d/vfOALaq6trbtDSOZ2AjKMA6TLyIfAHn4olb/EwikJ3kIn6nsPX+U5l0cTm29Gl/U5u/hS8HwtKpWicj7wKf4soq+GeO8DwP3i0gpMNyJoKr6koicAPzHJw4lwA/xjZgcnaeWeahHReQQ0Aj4P3wJ+wwjpVg0c8OoA6GpOtIti2FkG2biMwzDMDyJjaAMwzAMT2IjKMMwDMOTmIIyDMMwPIkpKMMwDMOTmIIyDMMwPIkpKMMwDMOT/H/BE3asTkoj6QAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.scatter(departments.department_id, departments.is_woman, label='Department averages')\n",
"plt.plot(departments.department_id, null_model_predictions, label='Overall average')\n",
"plt.scatter(departments.department_id, fancier_model_predictions, label='Prediction based on hiring in topics')\n",
"plt.legend(frameon=False)\n",
"plt.ylabel('Women to Men Ratio')\n",
"plt.xlabel('Department ID')\n",
"plt.ylim(0, 0.5)\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Is this fancier model better than the simpler? An ANOVA test will test whether one of two nested models are better. The p-value says the fancier model results in a significant increase to the R^2."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>df_resid</th>\n",
" <th>ssr</th>\n",
" <th>df_diff</th>\n",
" <th>ss_diff</th>\n",
" <th>F</th>\n",
" <th>Pr(&gt;F)</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>98.0</td>\n",
" <td>0.487476</td>\n",
" <td>0.0</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>94.0</td>\n",
" <td>0.338255</td>\n",
" <td>4.0</td>\n",
" <td>0.149221</td>\n",
" <td>10.36697</td>\n",
" <td>5.343899e-07</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" df_resid ssr df_diff ss_diff F Pr(>F)\n",
"0 98.0 0.487476 0.0 NaN NaN NaN\n",
"1 94.0 0.338255 4.0 0.149221 10.36697 5.343899e-07"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"simple_res = smf.ols('is_woman ~ 1', departments).fit()\n",
"fancier_res = smf.ols(\n",
" 'is_woman ~ nlp_researchers + ai_researchers + ml_researchers + hci_researchers', departments).fit()\n",
"sm.stats.anova_lm(simple_res, \n",
" fancier_res)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Alright let's do this by hand based on [its definition](https://en.wikipedia.org/wiki/F-test#Regression_problems). The F-statistic is $\\frac{(RSS_{1} - RSS_{2})/(p_{2} - p_{1})}{RSS_{2}/(n-p_{2})}$ where $RSS_{1}$ and $RSS_{2}$ are the residual sum squared for both models, $p_{1}$ and $p_{2}$ describe the number of parameters in each model, and $n$ is number of observations."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"rss_1 = sum((simple_res.resid)**2)\n",
"rss_2 = sum((fancier_res.resid)**2)\n",
"\n",
"p_1 = len(simple_res.params)\n",
"p_2 = len(fancier_res.params)\n",
"\n",
"n = len(departments)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"10.366970002659706\n"
]
}
],
"source": [
"f_stat = ((rss_1 - rss_2)/(p_2 - p_1))/(rss_2/(n - p_2))\n",
"print(f_stat) # Notice this the same as in the table above!"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"5.343899306398114e-07\n"
]
}
],
"source": [
"import scipy\n",
"p_value = 1-scipy.stats.f.cdf(f_stat, (p_2 - p_1), (n - p_2))\n",
"print(p_value)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice these same values are provided on the regression result for the fancier model. The F-statistic printed out always compares the model to an intercept only model, and the probability of that F-statistic is used to figure out whether our fancy model significantly reduces the residual sum squared compared to a cimple intercept only model."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>OLS Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>is_woman</td> <th> R-squared: </th> <td> 0.306</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>OLS</td> <th> Adj. R-squared: </th> <td> 0.277</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>Least Squares</td> <th> F-statistic: </th> <td> 10.37</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Tue, 15 Dec 2020</td> <th> Prob (F-statistic):</th> <td>5.34e-07</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>09:09:25</td> <th> Log-Likelihood: </th> <td> 140.64</td>\n",
"</tr>\n",
"<tr>\n",
" <th>No. Observations:</th> <td> 99</td> <th> AIC: </th> <td> -271.3</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Df Residuals:</th> <td> 94</td> <th> BIC: </th> <td> -258.3</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Df Model:</th> <td> 4</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>nonrobust</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>t</th> <th>P>|t|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> 0.3045</td> <td> 0.044</td> <td> 6.987</td> <td> 0.000</td> <td> 0.218</td> <td> 0.391</td>\n",
"</tr>\n",
"<tr>\n",
" <th>nlp_researchers</th> <td> 0.0073</td> <td> 0.002</td> <td> 4.376</td> <td> 0.000</td> <td> 0.004</td> <td> 0.011</td>\n",
"</tr>\n",
"<tr>\n",
" <th>ai_researchers</th> <td> 0.0021</td> <td> 0.002</td> <td> 1.009</td> <td> 0.316</td> <td> -0.002</td> <td> 0.006</td>\n",
"</tr>\n",
"<tr>\n",
" <th>ml_researchers</th> <td> -0.0062</td> <td> 0.001</td> <td> -4.271</td> <td> 0.000</td> <td> -0.009</td> <td> -0.003</td>\n",
"</tr>\n",
"<tr>\n",
" <th>hci_researchers</th> <td> 0.0018</td> <td> 0.002</td> <td> 0.880</td> <td> 0.381</td> <td> -0.002</td> <td> 0.006</td>\n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <th>Omnibus:</th> <td> 0.661</td> <th> Durbin-Watson: </th> <td> 2.097</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Prob(Omnibus):</th> <td> 0.718</td> <th> Jarque-Bera (JB): </th> <td> 0.697</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Skew:</th> <td>-0.189</td> <th> Prob(JB): </th> <td> 0.706</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Kurtosis:</th> <td> 2.838</td> <th> Cond. No. </th> <td> 192.</td>\n",
"</tr>\n",
"</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: is_woman R-squared: 0.306\n",
"Model: OLS Adj. R-squared: 0.277\n",
"Method: Least Squares F-statistic: 10.37\n",
"Date: Tue, 15 Dec 2020 Prob (F-statistic): 5.34e-07\n",
"Time: 09:09:25 Log-Likelihood: 140.64\n",
"No. Observations: 99 AIC: -271.3\n",
"Df Residuals: 94 BIC: -258.3\n",
"Df Model: 4 \n",
"Covariance Type: nonrobust \n",
"===================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"-----------------------------------------------------------------------------------\n",
"Intercept 0.3045 0.044 6.987 0.000 0.218 0.391\n",
"nlp_researchers 0.0073 0.002 4.376 0.000 0.004 0.011\n",
"ai_researchers 0.0021 0.002 1.009 0.316 -0.002 0.006\n",
"ml_researchers -0.0062 0.001 -4.271 0.000 -0.009 -0.003\n",
"hci_researchers 0.0018 0.002 0.880 0.381 -0.002 0.006\n",
"==============================================================================\n",
"Omnibus: 0.661 Durbin-Watson: 2.097\n",
"Prob(Omnibus): 0.718 Jarque-Bera (JB): 0.697\n",
"Skew: -0.189 Prob(JB): 0.706\n",
"Kurtosis: 2.838 Cond. No. 192.\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fancier_res.summary()"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.3061085772316079\n"
]
}
],
"source": [
"ssr = sum((departments.is_woman - fancier_res.predict())**2)\n",
"sst = sum((departments.is_woman - departments.is_woman.mean())**2)\n",
"r2 = 1 - (ssr/sst)\n",
"print(r2)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0\n"
]
}
],
"source": [
"ssr = sum((departments.is_woman - simple_res.predict())**2)\n",
"sst = sum((departments.is_woman - departments.is_woman.mean())**2)\n",
"r2 = 1 - (ssr/sst)\n",
"print(r2) # This should be zero because the sst == ssr. The simplest model is just predicting the average"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment