{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "89Q75QiBTuAR"
},
"source": [
"# Chi Square Test\n",
"\n",
"Notebook created for Regression in Psychology [PSYCH–GA.2229](https://docs.google.com/document/d/10AW7g92O6BtX61kXVIkHtL4j_k3P_G5f/edit?usp=sharing&ouid=100340169590558171318&rtpof=true&sd=true) graduate level course at New York University by [Dr. Madalina Vlasceanu](https://www.mvlasceanu.com/)\n",
"\n",
"This content is Open Access (free access to information and unrestricted use of electronic resources for everyone).\n",
"\n",
"Sources:\n",
"- Navarro, D. (2013). Learning statistics with R: https://learningstatisticswithr.com/\n",
"- Gureckis, 2018\n",
"https://teaching.gureckislab.org/fall22/labincp/intro.html"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"id": "xIiedSn8zTAK"
},
"outputs": [],
"source": [
"# import libraries\n",
"\n",
"import numpy as np\n",
"import statsmodels.api as sm\n",
"import pylab as py\n",
"import pandas as pd\n",
"import seaborn as sns\n",
"import matplotlib.pyplot as plt\n",
"import scipy.stats as stats\n",
"from scipy.stats import wilcoxon\n",
"from scipy.stats import chisquare\n",
"from statsmodels.stats.contingency_tables import mcnemar"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Unnamed: 0
\n",
"
PRE_male
\n",
"
PRE_female
\n",
"
POST_male
\n",
"
POST_female
\n",
"
POST_male_salary
\n",
"
POST_female_salary
\n",
"
POST_male_friendly
\n",
"
POST_female_friendly
\n",
"
POST_male_intelligent
\n",
"
...
\n",
"
Hire_m
\n",
"
Hire_f
\n",
"
Choice_m
\n",
"
Choice_f
\n",
"
itemnum
\n",
"
partnum
\n",
"
Gender
\n",
"
Ide
\n",
"
Political
\n",
"
Edu
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
0
\n",
"
1
\n",
"
1
\n",
"
1
\n",
"
1
\n",
"
57
\n",
"
88
\n",
"
95.0
\n",
"
73.0
\n",
"
90.0
\n",
"
...
\n",
"
1
\n",
"
1
\n",
"
1
\n",
"
1
\n",
"
0
\n",
"
0
\n",
"
Female
\n",
"
100
\n",
"
Democrat
\n",
"
College Degree
\n",
"
\n",
"
\n",
"
1
\n",
"
1
\n",
"
1
\n",
"
0
\n",
"
0
\n",
"
0
\n",
"
63
\n",
"
44
\n",
"
66.0
\n",
"
79.0
\n",
"
74.0
\n",
"
...
\n",
"
0
\n",
"
0
\n",
"
0
\n",
"
0
\n",
"
1
\n",
"
0
\n",
"
Female
\n",
"
100
\n",
"
Democrat
\n",
"
College Degree
\n",
"
\n",
" \n",
"
\n",
"
2 rows × 21 columns
\n",
"
"
],
"text/plain": [
" Unnamed: 0 PRE_male PRE_female POST_male POST_female POST_male_salary \\\n",
"0 0 1 1 1 1 57 \n",
"1 1 1 0 0 0 63 \n",
"\n",
" POST_female_salary POST_male_friendly POST_female_friendly \\\n",
"0 88 95.0 73.0 \n",
"1 44 66.0 79.0 \n",
"\n",
" POST_male_intelligent ... Hire_m Hire_f Choice_m Choice_f itemnum \\\n",
"0 90.0 ... 1 1 1 1 0 \n",
"1 74.0 ... 0 0 0 0 1 \n",
"\n",
" partnum Gender Ide Political Edu \n",
"0 0 Female 100 Democrat College Degree \n",
"1 0 Female 100 Democrat College Degree \n",
"\n",
"[2 rows x 21 columns]"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# import data downloaded from https://github.com/mvlasceanu/RegressionData/blob/main/data6.xlsx\n",
"#df = pd.read_excel('data6.xlsx')\n",
"\n",
"# Or you can read the Excel file directly from the URL\n",
"url = 'https://github.com/mvlasceanu/RegressionData/raw/main/data6.xlsx'\n",
"df = pd.read_excel(url)\n",
"\n",
"df.head(2)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "XXRt6DJl6rQ2"
},
"source": [
"## One way chi square"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "6-hpd7_FPwDY"
},
"source": [
"Chi square goodness-of-fit tests if a proportion is different from a baseline proportion.\n",
"\n",
"Useful when both outcome (DV) and predictor (IV) are categorical."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "M_yJlrh6v9SO"
},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "JwyUm0Gwxnou"
},
"source": [
"### Example"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "1GiR7Xi1xtwY"
},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "F9l3VJEb24fL"
},
"source": [
"The variable \"PRE_female\" encodes whether participants chose a man (PRE_female==0) or a woman (PRE_female==1) at pretest, in the female condition of this experiment. Since the choice is binary, the chance level of choosing a woman is 0.5. Let's see if participants chose women at chance (Null hypothesis) or if there was a bias towards choosing men (in which case we expect the mean of PRE_female to be < 0.5) or a bias towards choosing women (in which case we expect the mean of PRE_female to be > 0.5)."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "7xtsF5ivyBaQ"
},
"source": [
"To test whether the proportion of men/women choices at pretest in Condition 1 is at chance level (comparing PRE_female to 0.5) we can run a (one way) chi square test."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "3XRfBCygzlEN",
"outputId": "92061cac-bb18-4abe-ea20-3b6e2912fb30"
},
"outputs": [
{
"data": {
"text/plain": [
"0.38846153846153847"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# first let's look at the mean of the variable, which indicates percent women chosen\n",
"\n",
"df['PRE_female'].mean()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "xmdm3Ay10hcp",
"outputId": "04c4e73d-3bdd-4890-90ea-9165407137ac"
},
"outputs": [
{
"data": {
"text/plain": [
"0.4883404432118925"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# we can also compute the standard deviation of the proportion\n",
"\n",
"df['PRE_female'].std()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 497
},
"id": "r9k_1TPx0kYz",
"outputId": "70d94988-e14a-45dd-d918-c8e7cc01f154"
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/var/folders/53/6ypqgtmd38ggqwb6vvtqdz6h0000gp/T/ipykernel_97624/3444304390.py:5: UserWarning: \n",
"\n",
"`distplot` is a deprecated function and will be removed in seaborn v0.14.0.\n",
"\n",
"Please adapt your code to use either `displot` (a figure-level function with\n",
"similar flexibility) or `histplot` (an axes-level function for histograms).\n",
"\n",
"For a guide to updating your code to use the new functions, please see\n",
"https://gist.github.com/mwaskom/de44147ed2974457ad6372750bbe5751\n",
"\n",
" sns.distplot(df[\"PRE_female\"], ax=ax[1])\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAJOCAYAAACqbjP2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABfnUlEQVR4nO3deXyU9bn///dMJjtZSEIWIEBQNokiggoqoKIoHC2i/Wm1bqfaI3WrIsdKPcdaW0sXy0GrgLYKpR4t3xa0tVArp4Kg4AIEtaLIHgwJISzZSCbJzP37Y3IPxIRwJ8zMPcvr+XjMo81kJrmGgeTt9bnuz8dhGIYhAAAAdMppdwEAAACRgNAEAABgAaEJAADAAkITAACABYQmAAAACwhNAAAAFhCaAAAALCA0AQAAWOCyuwA7eL1e7du3T2lpaXI4HHaXA8QUwzBUW1ur3r17y+mM3f9u4+cQYK/u/CyKydC0b98+FRYW2l0GENP27t2rvn372l2Gbfg5BISHrvwsisnQlJaWJsn3B5Wenm5zNUBsqampUWFhof/fYazi5xBgr+78LIrJ0GS2wtPT0/lhBdgk1pek+DkEhIeu/CyK3YECAACALiA0AQAAWEBoAgAAsIDQBAAAYAGhCQAAwAJCEwAAgAWEJgAAAAsITQBi3uzZs3XuuecqLS1Nubm5uuaaa7R169ZOn7N69Wo5HI52ty+++CJEVQMINUITgJj3zjvv6J577tH777+vlStXqqWlRZMmTVJ9ff1Jn7t161aVl5f7b4MGDQpBxQDsEJM7ggPA8d588802Hy9cuFC5ubnauHGjxo8f3+lzc3NzlZmZGcTqAIQLOk0A8DXV1dWSpKysrJM+duTIkSooKNDEiRO1atWqEz7O7XarpqamzQ1AZCE0AcBxDMPQjBkzdNFFF6m4uPiEjysoKNALL7ygpUuXatmyZRoyZIgmTpyoNWvWdPj42bNnKyMjw38rLCwM1ksAECQOwzAMu4sItZqaGmVkZKi6upqDMoEQC/d/f/fcc4+WL1+ud999V3379u3Sc6+++mo5HA799a9/bfc5t9stt9vt/9g8YT1c/xyAaNedn0V0mgCg1X333ae//vWvWrVqVZcDkySNGTNG27Zt6/BziYmJSk9Pb3MDEFkYBAcQ8wzD0H333afXXntNq1evVlFRUbe+TklJiQoKCgJcHYBwQWgCEPPuuecevfLKK/rLX/6itLQ0VVRUSJIyMjKUnJwsSZo1a5bKysq0ePFiSdLcuXM1YMAADR8+XE1NTXr55Ze1dOlSLV261LbXASC4CE0AYt78+fMlSRdffHGb+xcuXKjbb79dklReXq7S0lL/55qamjRz5kyVlZUpOTlZw4cP1/LlyzVlypRQlQ0gxBgEZ64ACCn+/fnw5wDYi0FwAACAICE0AQAAWMBMEyLK97//fR04cECS1KtXLz399NM2VwQAiBWEJkSUAwcOaP/+/XaXAQCIQYQmAEDIvPJB6ckfFAA3nd8vJN8HsYWZJgAAAAsITQAAABYQmgAAACwgNAEAAFhAaAIAALCA0AQAAGABoQkAAMACQhMAAIAFhCYAAAALCE0AAAAWEJoAAAAsIDQBAABYQGgCAACwgNAEAABgAaEJAADAAkITAACABYQmAAAACwhNAAAAFhCaAAAALCA0AQAAWEBoAgAAsIDQBAAAYAGhCQAAwAJCEwAAgAWEJgAAAAsITQAAABYQmgAAACwgNAEAAFhAaAIAALCA0AQAAGABoQkAAMACQhMAAIAFhCYAAAALCE0AAAAWEJoAAAAsIDQBAABYQGgCAACwgNAEAABgAaEJAADAAkITAACABYQmAAAACwhNAAAAFhCaAAAALCA0AQAAWEBoAgAAsIDQBAAAYEFYhKZ58+apqKhISUlJGjVqlNauXWvpee+9955cLpfOPvvs4BYIAABinu2hacmSJXrggQf06KOPqqSkROPGjdPkyZNVWlra6fOqq6t16623auLEiSGqFAAAxDLbQ9OcOXN0xx136M4779SwYcM0d+5cFRYWav78+Z0+76677tJNN92ksWPHhqhSAAAQy2wNTU1NTdq4caMmTZrU5v5JkyZp3bp1J3zewoULtWPHDv3oRz8KdokAAACSJJed37yqqkoej0d5eXlt7s/Ly1NFRUWHz9m2bZseeeQRrV27Vi6XtfLdbrfcbrf/45qamu4XDQAAYpLty3OS5HA42nxsGEa7+yTJ4/Hopptu0o9//GMNHjzY8tefPXu2MjIy/LfCwsJTrhkAAMQWW0NTTk6O4uLi2nWVKisr23WfJKm2tlYbNmzQvffeK5fLJZfLpSeeeEIff/yxXC6X3n777Q6/z6xZs1RdXe2/7d27NyivBwAARC9bl+cSEhI0atQorVy5UtOmTfPfv3LlSk2dOrXd49PT0/Xpp5+2uW/evHl6++239ec//1lFRUUdfp/ExEQlJiYGtngAABBTbA1NkjRjxgzdcsstGj16tMaOHasXXnhBpaWlmj59uiRfl6isrEyLFy+W0+lUcXFxm+fn5uYqKSmp3f0AAACBZHtouuGGG3Tw4EE98cQTKi8vV3FxsVasWKH+/ftLksrLy0+6ZxMAAECw2R6aJOnuu+/W3Xff3eHnFi1a1OlzH3/8cT3++OOBLwoAAOA4YXH1HAAAQLgjNAEAAFhAaAIAALCA0AQAAGABoQkAAMACQhMAAIAFhCYAAAALwmKfpmgy6j8X211CVEs/XOdP+uWH6/jzDrKNv7rV7hIAIGzQaQIAALCA0AQAAGABoQkAAMACQhMAAIAFhCYAAAALCE0AAAAWEJoAAAAsIDQBAABYQGgCAACwgNAEAABgAaEJAADAAkITAACABYQmAAAACwhNAAAAFhCaAAAALCA0AQAAWEBoAgAAsIDQBAAAYAGhCQAAwAJCEwAAgAWEJgAAAAsITQBi3uzZs3XuuecqLS1Nubm5uuaaa7R169aTPu+dd97RqFGjlJSUpIEDB2rBggUhqBaAXQhNAGLeO++8o3vuuUfvv/++Vq5cqZaWFk2aNEn19fUnfM6uXbs0ZcoUjRs3TiUlJfrhD3+o+++/X0uXLg1h5QBCyWV3AQBgtzfffLPNxwsXLlRubq42btyo8ePHd/icBQsWqF+/fpo7d64kadiwYdqwYYOeeuopXXfddcEuGYAN6DQBwNdUV1dLkrKysk74mPXr12vSpElt7rviiiu0YcMGNTc3B7U+APag0wQAxzEMQzNmzNBFF12k4uLiEz6uoqJCeXl5be7Ly8tTS0uLqqqqVFBQ0OZzbrdbbrfb/3FNTU1gCwcQdHSaAOA49957rz755BO9+uqrJ32sw+Fo87FhGB3eL/mGzTMyMvy3wsLCwBQMIGQITQDQ6r777tNf//pXrVq1Sn379u30sfn5+aqoqGhzX2VlpVwul7Kzs9s9ftasWaqurvbf9u7dG9DaAQQfy3MAYp5hGLrvvvv02muvafXq1SoqKjrpc8aOHas33nijzX1vvfWWRo8erfj4+HaPT0xMVGJiYsBqBhB6dJoAxLx77rlHL7/8sl555RWlpaWpoqJCFRUVamho8D9m1qxZuvXWW/0fT58+XXv27NGMGTP0+eef66WXXtKLL76omTNn2vESAIQAoQlAzJs/f76qq6t18cUXq6CgwH9bsmSJ/zHl5eUqLS31f1xUVKQVK1Zo9erVOvvss/WTn/xEzzzzDNsNAFGM5TkAMc8c4O7MokWL2t03YcIEbdq0KQgVAQhHdJoAAAAsIDQBAABYQGgCAACwgNAEAABgAaEJAADAAkITAACABYQmAAAACwhNAAAAFhCaAAAALCA0AQAAWEBoAgAAsIDQBAAAYAGhCQAAwAJCEwAAgAWEJgAAAAsITQAAABYQmgAAACwgNAEAAFhAaAIAALCA0AQAAGABoQkAAMACQhMAAIAFhCYAAAALCE0AAAAWEJoAAAAsIDQBAABYQGgCAACwgNAEAABgAaEJAADAAkITAACABYQmAAAACwhNAAAAFhCaAAAALCA0AQAAWEBoAgAAsIDQBAAAYAGhCQAAwAKX3QUAXeGNT+3w/wMAEGyEJkSUuiGT7S4BABCjWJ4DAACwgNAEAABgAaEJAADAAkITAACABYQmAAAACwhNAAAAFhCaAAAALCA0AQAAWEBoAgAAsIDQBAAAYAGhCQAAwAJCEwAAgAWEJgAAAAsITQAAABYQmgAAACwgNAEAAFhAaAIAALCA0AQAAGABoQkAAMACQhMAAIAFhCYAAAALCE0AAAAWEJoAAAAsCIvQNG/ePBUVFSkpKUmjRo3S2rVrT/jYd999VxdeeKGys7OVnJysoUOH6n/+539CWC0AAIhFLrsLWLJkiR544AHNmzdPF154oZ5//nlNnjxZW7ZsUb9+/do9PjU1Vffee6/OOusspaam6t1339Vdd92l1NRU/cd//IcNrwAAAMQC2ztNc+bM0R133KE777xTw4YN09y5c1VYWKj58+d3+PiRI0fqxhtv1PDhwzVgwADdfPPNuuKKKzrtTgEAAJwqW0NTU1OTNm7cqEmTJrW5f9KkSVq3bp2lr1FSUqJ169ZpwoQJJ3yM2+1WTU1NmxsAAEBX2Bqaqqqq5PF4lJeX1+b+vLw8VVRUdPrcvn37KjExUaNHj9Y999yjO++884SPnT17tjIyMvy3wsLCgNQPAABih+3Lc5LkcDjafGwYRrv7vm7t2rXasGGDFixYoLlz5+rVV1894WNnzZql6upq/23v3r0BqRsAAMQOWwfBc3JyFBcX166rVFlZ2a779HVFRUWSpDPPPFP79+/X448/rhtvvLHDxyYmJioxMTEwRQMAgJhka6cpISFBo0aN0sqVK9vcv3LlSl1wwQWWv45hGHK73YEuDwAAwM/2LQdmzJihW265RaNHj9bYsWP1wgsvqLS0VNOnT5fkW1orKyvT4sWLJUnPPfec+vXrp6FDh0ry7dv01FNP6b777rPtNQAAgOhne2i64YYbdPDgQT3xxBMqLy9XcXGxVqxYof79+0uSysvLVVpa6n+81+vVrFmztGvXLrlcLp122mn6+c9/rrvuusuulwAAAGKAwzAMw+4iQq2mpkYZGRmqrq5Wenp6QL/2qP9cHNCvB9hp469uDfjXDOa/v0gSq38Or3xQevIHBcBN57ffHBk4Xnf+DYbF1XMAAADhjtAEAABgAaEJAADAAkITAACABYQmAAAACwhNAAAAFhCaAAAALCA0AQAAWEBoAgAAsIDQBAAAYAGhCQAAwAJCEwAAgAUuuwsAAACd46Dj8ECnCQAAwAJCEwAAgAWEJgAAAAsITQAAABYQmgAAACzg6jkAQFB5vIZ+t3anXispU3pSvKacWaAEF//NjsjD31oAMW/NmjW6+uqr1bt3bzkcDr3++uudPn716tVyOBztbl988UVoCo4wz6/Zodl//0JfVNTqw92H9PIHe2QYht1lAV1GaAIQ8+rr6zVixAg9++yzXXre1q1bVV5e7r8NGjQoSBVGrj0H6zV35TZJ0kWn5yjO6dD2yjp9ub/O5sqArmN5DkDMmzx5siZPntzl5+Xm5iozMzPwBUWRVz4oVZPHqwtOy9Yf7jhPt774odZur9I/v9ivIflpdpcHdAmdJgDoppEjR6qgoEATJ07UqlWrOn2s2+1WTU1Nm1u0a/F4taykTJJ02wUD5HA4NG5wLzkd0leHG3Swzm1zhUDXEJoAoIsKCgr0wgsvaOnSpVq2bJmGDBmiiRMnas2aNSd8zuzZs5WRkeG/FRYWhrBie7y7vUoHat3KTk3QpUNzJUk9El0amNNDkvSvfdEfHBFdWJ4DgC4aMmSIhgwZ4v947Nix2rt3r5566imNHz++w+fMmjVLM2bM8H9cU1MT9cFpzZdVkqRJw/MVH3fsv9GL+2Ro+4E6/ausWhMG97KrPKDL6DQBQACMGTNG27ZtO+HnExMTlZ6e3uYW7d7dfkCSNG5QTpv7hxX4Zpn2HWlQQ5Mn5HUB3UVoAoAAKCkpUUFBgd1lhI39NY36cn+dHA5p7MDsNp9LS4pXTo9EGZJ2H6y3p0CgG1ieAxDz6urqtH37dv/Hu3bt0ubNm5WVlaV+/fpp1qxZKisr0+LFiyVJc+fO1YABAzR8+HA1NTXp5Zdf1tKlS7V06VK7XkLYeX/nQUlSce8M9UxNaPf5opwUVdW5tauqXsMKor/rhuhAaAIQ8zZs2KBLLrnE/7E5e3Tbbbdp0aJFKi8vV2lpqf/zTU1NmjlzpsrKypScnKzhw4dr+fLlmjJlSshrD1clpUckSaMH9Ozw80U5qfpo92HtqqLThMhxyqFp+/bt2rFjh8aPH6/k5GQZhiGHwxGI2gAgJC6++OJOd6hetGhRm48ffvhhPfzww0GuKrJ98tURSdKIvpkdfr5/dqokqby6Qc0eb5tBcSBcdftv6cGDB3XZZZdp8ODBmjJlisrLyyVJd955px566KGAFQgAiCzNHq8+a91OYERhZoePyUyOV0pCnLyGb/4JiATdDk0PPvigXC6XSktLlZKS4r//hhtu0JtvvhmQ4gAAkWdrRa3cLV6lJ7k0IDulw8c4HA71zkyWJJUfITQhMnR7ee6tt97SP/7xD/Xt27fN/YMGDdKePXtOuTAAQGT6tKxaknRW38xOxzV6ZyRre2WdyqobdG6oigNOQbc7TfX19W06TKaqqiolJiaeUlEAgMi1taJW0rH9mE6kd2aSJN9+TUAk6HZoGj9+vP/yW8nXavV6vfrVr37V5ioUAEBs+XK/LzQNzjtJaMrwLc9VVDfK28kgPhAuur0896tf/UoXX3yxNmzYoKamJj388MP67LPPdOjQIb333nuBrBEAEEHM0DQkv/PQlNUjQS6nQy1eQ0eONiurg/2cgHDS7U7TGWecoU8++UTnnXeeLr/8ctXX1+vaa69VSUmJTjvttEDWCACIEFV1blXVNcnhkAbldh6anA6HeqX5xjkquYIOEeCU9mnKz8/Xj3/840DVAgCIcF+2zjP1z0pRckLcSR/fKy1R5dWNqqx1ayin0CDMdSk0ffLJJ5Yfe9ZZZ3W5GABAZNtWWSdJGnSSeSZTblqSpGpV1tJpQvjrUmg6++yz5XA4Ot05V/INhXs8nFwNALFm5wFfaDqtVw9Lj881l+dq3UGrCQiULoWmXbt2BasOAEAU2Nl6ltzAnFRLj89NN2ea3BzDhbDXpdDUv3//YNUBAIgC5gG8Rb2shaas1AQ5JDV5vKp1tyg9KT6I1QGn5pQP7N2yZYtKS0vV1NTU5v5vfOMbp/qlAQARpLHZo7LWjSqLLHaaXE6neqYm6FB9kw7WNRGaENa6HZp27typadOm6dNPP20z52S2VplpAoDYUnroqAxDSktyKbsLey5l+0OT23LYAuzQ7X2avv/976uoqEj79+9XSkqKPvvsM61Zs0ajR4/W6tWrA1giAHSMOcvwYg6BD8xJ7dJskrmp5cH6ppM8ErBXt0PT+vXr9cQTT6hXr15yOp1yOp266KKLNHv2bN1///2BrBEAOnT66afrkksu0csvv6zGRi5Zt9vug0clSQO62C3K7uEbBic0Idx1OzR5PB716OG7pDQnJ0f79u2T5BsW37p1a2CqA4BOfPzxxxo5cqQeeugh5efn66677tKHH35od1kxq/SQLzT1y2p/mHtncsxOUx3bDiC8dTs0FRcX+ze7PP/88/XLX/5S7733np544gkNHDgwYAUCwIkUFxdrzpw5Kisr08KFC1VRUaGLLrpIw4cP15w5c3TgwAG7S4wpe1tDU2EXQ1NWj2PLcyfbBxCwU7dD03/913/J6/VKkn76059qz549GjdunFasWKFnnnkmYAUCwMm4XC5NmzZN/+///T/94he/0I4dOzRz5kz17dtXt956q8rLy+0uMSZ0t9OUleILTU0tXh1t4iIihK9uXz13xRVX+P//wIEDtWXLFh06dEg9e/ZkczIAIbVhwwa99NJL+uMf/6jU1FTNnDlTd9xxh/bt26fHHntMU6dOZdkuyDxeQ2WHfdsNdLXT5IpzKi3RpVp3i44cbVZq4invhgMERUD/ZmZlZQXyywFAp+bMmaOFCxdq69atmjJlihYvXqwpU6bI6fQ10YuKivT8889r6NChNlca/cqrG9TiNRQf51B+elKXn5+ZEq9ad4sOH21Sn57JQagQOHXdDk2NjY36zW9+o1WrVqmystK/VGfatGnTKRcHAJ2ZP3++vvOd7+jf//3flZ+f3+Fj+vXrpxdffDHElcUec2mub88UxTm7vtqQmZKgvYcbdKShOdClAQHT7dD0ne98RytXrtQ3v/lNnXfeeSzJAQi5lStXql+/fv7OkskwDO3du1f9+vVTQkKCbrvtNpsqjB3dHQI3Zab4dgI/cpRtBxC+uh2ali9frhUrVujCCy8MZD0AYNlpp52m8vJy5ebmtrn/0KFDKioq4mSCEDLnmfp2c2kts3UY/PBROk0IX92+eq5Pnz5KS0sLZC0A0CUnujy9rq5OSUldn6tB95Ud8W0u2ieze6GpZzKdJoS/bneafv3rX+sHP/iBFixYoP79+weyJgDo1IwZMyT5zrp87LHHlJJybEnI4/Hogw8+0Nlnn21TdbGpvNrXaSrI6F5YNTtNR+g0IYx1OzSNHj1ajY2NGjhwoFJSUhQf3/Zk6kOHDp1ycQDQkZKSEkm+TtOnn36qhIRjh8MmJCRoxIgRmjlzpl3lxaTyal+nqSCju8tzvt8hDc0euZs9SoyPC1htQKB0OzTdeOONKisr089+9jPl5eUxCA4gZFatWiVJ+vd//3c9/fTTSk9Pt7mi2GYYhvYd8XWaurs8lxQfp+T4ODU0e3S4oVn5hCaEoW6HpnXr1mn9+vUaMWJEIOsBAMsWLlxodwmQdKi+Se4W37YzeRmJ3f46mSnxaqj26MjRpm7t9QQEW7dD09ChQ9XQ0BDIWgDgpK699lotWrRI6enpuvbaazt97LJly0JUVWwzl+ZyeiQq0dX9DlFmSoLKqxuZa0LY6nZo+vnPf66HHnpITz75pM4888x2M020ywEEQ0ZGhn8cICMjw+ZqIEll/qW5U+sOsVcTwl23Q9OVV14pSZo4cWKb+w3DkMPhYH8UAEFx/JIcy3PhofyIeeXcqR1/Ym47wF5NCFfdDk3mICYA2KWhoUGGYfi3HNizZ49ee+01nXHGGZo0aZLN1cUOc3mudzeHwE3Hth2g04Tw1O3QNGHChEDWAQBdNnXqVF177bWaPn26jhw5ovPOO08JCQmqqqrSnDlz9L3vfc/uEmOCuTzXO1DLc5w/hzDV7R3BJWnt2rW6+eabdcEFF6isrEyS9Ic//EHvvvtuQIoDgM5s2rRJ48aNkyT9+c9/Vn5+vvbs2aPFixfrmWeesbm62HGqezSZzE5TbWOLWjzekzwaCL1uh6alS5fqiiuuUHJysjZt2iS32y1Jqq2t1c9+9rOAFQgAJ3L06FH/cU5vvfWWrr32WjmdTo0ZM0Z79uyxubrYUR6gTlNqQpxcTt+Qf01jyynXBQRat0PTT3/6Uy1YsEC//e1v21w5d8EFF2jTpk0BKQ4AOnP66afr9ddf1969e/WPf/zDP8dUWVnJFbwh0uLxqqImMDNNDodD6a3D4DUs0SEMdTs0bd26VePHj293f3p6uo4cOXIqNQGAJY899phmzpypAQMG6Pzzz9fYsWMl+bpOI0eOtLm62FBZ65bXkFxOh3J6dH9jS1Nakm/UtqaR0ITw0+1B8IKCAm3fvl0DBgxoc/+7776rgQMHnmpdAHBS3/zmN3XRRRepvLy8zekEEydO1LRp02ysLHaYB/XmZyQpznnqx2mlJ7V2mlieQxjqdmi666679P3vf18vvfSSHA6H9u3bp/Xr12vmzJl67LHHAlkjAJxQfn6+8vPz29x33nnn2VRN7Nl3pHVp7hSHwE0ZLM8hjHU7ND388MOqrq7WJZdcosbGRo0fP16JiYmaOXOm7r333kDWCAAdqq+v189//nP985//VGVlpbzetldc7dy506bKYse+AA2Bm1ieQzjrUmj65JNPVFxcLKfTNwr15JNP6tFHH9WWLVvk9Xp1xhlnqEePHkEpFAC+7s4779Q777yjW265RQUFBf7jVRA65hB4XkZgQtOxQXCW5xB+uhSaRo4cqfLycuXm5mrgwIH66KOPlJ2drdGjRwerPgA4ob///e9avny5LrzwQrtLiVmVNb7tZvLTAxSa/DNNdJoQfrp09VxmZqZ27dolSdq9e3e7VjgAhFLPnj2VlZVldxkxrbLW12nKTQtMaDp+pskwjIB8TSBQutRpuu666zRhwgR/G3z06NGKi4vr8LHMEgAItp/85Cd67LHH9Pvf/95//hxCa39rpykv/dS3G5COzTS1eA01NnuVnNDx7xjADl0KTS+88IKuvfZabd++Xffff7+++93v+nfjBYBQ+/Wvf60dO3YoLy9PAwYMaLPRriQ22g0ywzC035xpCtDyXHycU8nxcWpo9qi6sZnQhLDS5avnrrzySknSxo0b9f3vf/+koemrr75S7969/cPjABAo11xzjd0lxLSahha5W3xjGr3SAtNpknxLdA3NHtU0NAdsVgoIhG5vObBw4UJLjzvjjDO0efNmNrwEEHA/+tGP7C4hppnzTBnJ8UqKD1xHKC3JpYoaqZZhcISZoLd/GOQDEExHjhzR7373O82aNUuHDh2S5FuWKysrs7my6BfoeSaTue1ANdsOIMx0u9MEAHb75JNPdNlllykjI0O7d+/Wd7/7XWVlZem1117Tnj17tHjxYrtLjGqBnmcyse0AwhWDRgAi1owZM3T77bdr27ZtSko69ot78uTJWrNmjY2VxYb9rctzgZxnkqT0ZN9/z9dylArCDKEJQMT66KOPdNddd7W7v0+fPqqoqLChothS6V+eC1anieU5hJeghyaONQAQLElJSaqpqWl3/9atW9WrVy8bKoot5iB4XsA7TRzai/DEIDiAiDV16lQ98cQTam72/XJ1OBwqLS3VI488ouuuu87m6qLf/qB1mnzLc3XuFnm8/A5B+Ah4aDIMQ5WVlf6Pt2zZov79+wf62wCAnnrqKR04cEC5ublqaGjQhAkTdPrppystLU1PPvmk3eVFPXMQPDfAV8+lJrrkdEiGfMEJCBddvnouJSVFe/bs8be+r7zySi1cuFAFBQWSpMrKSvXu3Vsej0eSVFhYGMByAeCY9PR0vfvuu1q1apU2btwor9erc845R5dddpndpUU9wzBUWevrNAXq3DmT0+FQj0SXahpbVNvY7D+PDrBbl0NTY2NjmyW39957Tw0NDW0ew5IcgGDzer1atGiRli1bpt27d8vhcKioqEj5+fkyDIN5yiCrbmhWU+tu4IHuNElSWlK8ahpbVMcwOMJIUGaa+GEFIJgMw9A3vvEN3XnnnSorK9OZZ56p4cOHa8+ePbr99ts1bdo0u0uMeuY8U8+UeCW6An8+XI/EY3NNQLhgc0sAEWfRokVas2aN/vnPf+qSSy5p87m3335b11xzjRYvXqxbb73Vpgqjn3+eKcBLc6YercPgtYQmhJEud5ocDkebTtLXPwaAYHv11Vf1wx/+sF1gkqRLL71UjzzyiP73f//Xhspih3+eKQhLc5KU1tppqmV5DmGky6HJMAwNHjxYWVlZysrKUl1dnUaOHOn/eOjQocGoEwD8PvnkE1155ZUn/PzkyZP18ccfh7Ci2BOsI1RMZqepjqNUEEa6vDy3cOHCYNQBAJYdOnRIeXl5J/x8Xl6eDh8+HMKKYs8B/5VzQeo0te4KzvIcwkmXQ9Ntt90WjDoAwDKPxyOX68Q/vuLi4tTSwi/bYDJDU06P4IQm/yA4y3MIIwEfBC8vL9eTTz6pZ599NtBfGgAk+cYEbr/9diUmdvwL2+12h7ii2HOgrjU0Ba3TxCA4wk+3QtOWLVu0atUqxcfH6/rrr1dmZqaqqqr05JNPasGCBSoqKgp0nQDgZ6XjzZVzwVVlhqYeCUH5+uYgeFOLV+4WT1C2NQC6qsuh6W9/+5uuu+46/1lPv/zlL/Xb3/5W119/vYqLi/WnP/1JV111VcALBQATs5X2q2pdnusVpOW5BJdT8XEONXsM1TW2KLEHoQn26/LVc08++aSmT5+umpoaPfXUU9q5c6emT5+upUuXatWqVQQmAIhy7haPalpnjXoFaXnO4XD4h8HZ4BLhosuh6fPPP9c999yjHj166P7775fT6dTcuXM1fvz4YNQHAAgzB+uaJEnxcY6gngvXg72aEGa6HJpqamqUmZkpSXK5XEpOTtbgwYMDXRcAIEyZ80zZqYlB3dyYYXCEm24PgldUVEjyXcWydetW1dfXt3nMWWedderVAQDCjn8IPC04Q+Amth1AuOlWaLr00kvbfGzOMTkcDv/p4h6P59SrAwCEnapa3/JcsPZoMpmdpjo3u4IjPHR5eW7Xrl0nvO3cudP/v10xb948FRUVKSkpSaNGjdLatWtP+Nhly5bp8ssvV69evZSenq6xY8fqH//4R1dfBgCgm/x7NAU5NPVIbN0VnE4TwkSXQ1Nubq5++ctf6oILLtC5556rWbNmKTU1Vf37929zs2rJkiV64IEH9Oijj6qkpETjxo3T5MmTVVpa2uHj16xZo8svv1wrVqzQxo0bdckll+jqq69WSUlJV18KAKAbqkIUmo51mghNCA9dDk2PPfaYFi1apH/7t3/Tt771La1cuVLf+973ul3AnDlzdMcdd+jOO+/UsGHDNHfuXBUWFmr+/PkdPn7u3Ll6+OGHde6552rQoEH62c9+pkGDBumNN97odg0AAOuq6szludDMNNFpQrjo8kzTsmXL9OKLL+pb3/qWJOnmm2/WhRdeKI/Ho7i4rm0+1tTUpI0bN+qRRx5pc/+kSZO0bt06S1/D6/WqtrZWWVlZXfreAIDu8W9sGaQ9mkz+TlNji39eFrBTlztNe/fu1bhx4/wfn3feeXK5XNq3b1+Xv3lVVZU8Hk+708rz8vL8V+edzK9//WvV19fr+uuvP+Fj3G63ampq2twAAN0TupkmX2jyGIYamrm4CPbrcmjyeDxKSGjbknW5XKd0ovjX/+vB6n9RvPrqq3r88ce1ZMkS5ebmnvBxs2fPVkZGhv9WWFjY7VoBINaFaqbJFedUcrxvBYMlOoSDLi/PdXS6eGNjo6ZPn67U1FT/fcuWLTvp18rJyVFcXFy7rlJlZWW77tPXLVmyRHfccYf+9Kc/6bLLLuv0sbNmzdKMGTP8H9fU1BCcAKAbmj1eHTnq2wIg2DNNktQjyaWGZo/q3C3q/LcCEHxdDk0dnS5+8803d+ubJyQkaNSoUVq5cqWmTZvmv3/lypWaOnXqCZ/36quv6jvf+Y5effVV/du//dtJv09iYmKbkAcA6B7zCJU4p0M9U4IfmtISXTpQ66bThLDQ5dAU6NPFZ8yYoVtuuUWjR4/W2LFj9cILL6i0tFTTp0+X5OsSlZWVafHixZJ8genWW2/V008/rTFjxvi7VMnJycrIyAhobQCAtsyluazUBDmdwR/M7uEfBmeDS9ivWzuCB9INN9yggwcP6oknnlB5ebmKi4u1YsUK/15P5eXlbfZsev7559XS0qJ77rlH99xzj//+2267TYsWLQp1+QAQU0I1BG5KS+T8OYQP20OTJN199926++67O/zc14PQ6tWrg18QAKBD5nYDoZhnkqS0JN+u4Jw/h3DQ5avnAACxy9zYMth7NJl6sCs4wgihCQBgmTnT1CvUy3N0mhAGCE0AAMtCtUeTyew0MdOEcEBoAgBY5g9NaaGZaTJ3BT/qbpHHa4TkewInQmgCAFhWVWse1huaTlNqoksOSYak+ia6TbAXoQkAYFmol+ecDoe/28QVdLAboQkAYEmLx6tDR0PbaZK4gg7hg9AEALDk0NEmGYbkdPh2BA8VOk0IF4QmAIAl5jxTVmqC4kJwhIrJH5roNMFmhCYAMW/NmjW6+uqr1bt3bzkcDr3++usnfc4777yjUaNGKSkpSQMHDtSCBQuCX6jNQj3PZOrh36uJ8+dgL0ITgJhXX1+vESNG6Nlnn7X0+F27dmnKlCkaN26cSkpK9MMf/lD333+/li5dGuRK7XWg1qbQxEwTwkRYnD0HAHaaPHmyJk+ebPnxCxYsUL9+/TR37lxJ0rBhw7RhwwY99dRTuu6664JUpf2OdZpCN88ksTyH8EGnCQC6aP369Zo0aVKb+6644gpt2LBBzc3Ru4Rk2/IcnSaECTpNANBFFRUVysvLa3NfXl6eWlpaVFVVpYKCgnbPcbvdcrvd/o9ramqCXmegmYf15oTosF4TV88hXNBpAoBucDjaXj1mGEaH95tmz56tjIwM/62wsDDoNQaa3YPgR5s8HKUCWxGaAKCL8vPzVVFR0ea+yspKuVwuZWdnd/icWbNmqbq62n/bu3dvKEoNqGOD4KGdaTr+KJWjHKUCG7E8BwBdNHbsWL3xxhtt7nvrrbc0evRoxcfHd/icxMREJSaGtkMTaP7luRB3mpwOh1ISXap3t6jO3aK0pI7/jIFgo9MEIObV1dVp8+bN2rx5syTflgKbN29WaWmpJF+X6NZbb/U/fvr06dqzZ49mzJihzz//XC+99JJefPFFzZw5047yQ8LjNXSo3tdpyg3xTJMkpTHXhDBApwlAzNuwYYMuueQS/8czZsyQJN12221atGiRysvL/QFKkoqKirRixQo9+OCDeu6559S7d28988wzUb3dwOGjTfIakiPER6iY2HYA4YDQBCDmXXzxxf5B7o4sWrSo3X0TJkzQpk2bglhVeDGHwHumJMgVF/pFCnPbgVo6TUHR0ORRZW2jahqblc7y5wkRmgAAJ2WeOxfqIXATnabg2bTnsN74ZJ/cLV4tWrdbP7/uTE0b2dfussISM00AgJOya7sBE6EpOLZX1mnppq/kbvFKktwtXs34fx9r1ReVNlcWnghNAICTIjRFnxaPV8s2fSVD0jn9MvXTa4p1w+hCGYb033/5lxqbPXaXGHYITQCAkzpgd2hK4uq5QPtw9yEdaWhWepJLU8/uI6fDoR994wwVZCTpq8MN+sP6PXaXGHYITQCAk/LPNKUx0xQNvIah97ZXSZIuGZqr+Nbh/pQEl+6fOEiS9If397AD+9cQmgAAJ2X78lxrp6ne3SJvJ1c6wpqdB+p1+GizkuKdOqdfzzafu+bsPkpPcqn00FG98yWzTccjNAEATsoMTb1sCk2pCccfpcKszanauOeQJGlE30x/l8mUnBCnb47ynY34l837Ql5bOCM0AQBO6ti5c/aEpjinQ8kJcZKYazpVLR6vPq+olaR2XSbTVSMKJEn//LySgfDjEJoAAJ3yeg0drLd3pklirilQdlbVq6nFq/Qkl/r0TO7wMWf3zVRBRpLq3C1a8+WBEFcYvghNAIBOHWlo9g8EZ6fad+jwsV3Bm22rIRpsKa+RJA3NT5fT4ejwMU6nQ1cMz5ckrdrKXJOJ0AQA6JQ5z5SRHK8El32/Nug0nTrDMPRl69Lc0IK0Th87YXAvSdKaL6s6PWYolhCaAACdqvLPM9m3NCdJaYSmU3aovklHGpoV53BoYE6PTh97/sAsxcc5VHakQbsPHg1RheGN0AQA6JTdG1ua/J0mBsG7beeBeklSYVbySbuGKQkujervGxRfu425JonQBAA4iao63xB4rzSbQ1NSvCQ6TadiR1WdJGlgr867TKYLT8uRJH2w61DQaookhCYAQKfs3tjSxEzTqdtd5es0DeyVaunx5xZlSZI27D7EXJMITQCAkzBnmuzvNBGaTkV1Q7NqGlvkdEh9M1MsPce3+aVD+2vc+upwQ5ArDH+EJgBAp451muwdBDc7TRyl0j1fHfYNc+elJ1m+CjI5IU7FfTIkSR/tZomO0AQA6JQ50xQuy3NeQ2rgKJUuMztFfXta6zKZRrcOg28qPRzwmiINoQkA0KlwmWmKczqUHN96lApLdF22t7XTVHiCXcBP5Ky+mZKkf5XVBLqkiENoAgCckGEYOmh2mmyeaZKYa+our2GorJudJnN57vPyGrV4vAGvLZIQmgAAJ1TT0KKm1l+U2an2zjRJx5boatmrqUsO1LrlbvEqIc6p3PSuhd/+WSlKS3TJ3eLV9gN1QaowMhCaAAAnZG5smZbkUlLr0pid2Hage8x5pt6ZySc8b+5EnE6HzuidLkn69KvqgNcWSQhNAIATMueZetk8z2TyL8/RaeqSr7o5z2Q6s3WJ7rN9sT3XRGgCAJzQgdrwGAI3cf5c9/ivnMvq2jyTyZxr+lcZnSYAADrkv3Iuzf55Jun45blmmyuJHB6voYrqRklSn8zudZqKj+s0ebyxu0cWoQkAcELhst2Aiavnuq6qzi2PYSjR5VTPlPhufY2inFSlJMSpodmjXVWxOwxOaAIAnFBVbXhsbGnyd5qYabJsf42vy5SXniRHF4fATXFOh4abw+AxvERHaAIAnFDYdZr8R6l4OEDWogp/aDq193B4b3OuKXaHwQlNAIATCpdz50xmaPIYhhqaOUrFiv3VxzpNp2JYQZok6cv9tadcU6QiNAEATsg8d65XGOwGLkmuOKeS4n2/uliis8bsNOWfYmgalEdoIjQBADpkGIZ/c8twWZ6TpB6JvmHmWobBT8rd7NHho74rDU85NOX2kCTtr3GruiE2r14kNAEAOlTrblFTi+8IlXDpNEnsCt4V+2uP7eie0vrn1l1pSfHqneELXttitNtEaAIAdKiq9Rduj8TwOELFxK7g1pnzTKfaZTIdW6KLzW0HCE0AgA6Z80zhMgRuotNkXUVNYIbATYPzfEt0sTrXRGgCAHQo3LYbMBGarAvUELjJ7DRtqyQ0AQDgF66hKY3lOcv8G1tmBKrTxPIcAADtmDNN4XLunIlOkzX17hYdbfLIIalXgIKveQXdgVq3jhxtCsjXjCSEJgBAhw7UhdcRKiZCkzUHWkNvZkq8ElyB+XWfmujyH/obi90mQhMAoEPhujx3/KG9HKVyYsHaYyuWh8EJTQCADpmdirALTeZRKl5Djc1em6sJX+b7F+g9tsy5pljcq4nQBADokPlLN/cUD3oNtPg4pxJbl5tq3bG5M7UVwQpNp7fONe04UB/QrxsJCE0AgHaOP0IlUEPEgcRc08kF6/0b2MsXmnYeYKYJAADVNIbnESomdgXvXLPHq8P1wTlseWBOqiRpX3WjjjbF1p8/oQkA0M6B484sC6cjVEx0mjp3sL5JhqSkeKf/zypQeqYmqGeK79DkXVWxtURHaAIAtBOseZhAITR1zv/+9UiUw+EI+Nc/tkRHaAIAxLhwnmeS2BX8ZIJ95WNR6xIdnSYAQMwL/06Tb3mITlPHzD22gvX+DezlC02xNgxOaAIAtBP+oYnluc4E+/0bmNO6PEenCQAQ68I+NCURmk4kFNtFnObvNNXH1K7shCYAQDvhPtPk7zQ1cpTK19W6fdtFOCRl9QjOYcv9slPkdPhCqxmwYwGhCQDQTth3mlpDU4vXkLuFo1SOd7D1oOXMlHi5nMH5NZ/oilPfnimSYmuJjtAEAGgn3ENTgsuphNajVLiCrq1D9b73LjvIXcKBxy3RxQpCEwCgDY/X8P/iDdfQJB3rNtUy19SG2WnKTg3O0pzJPwweQ1fQEZoAAG0crHfLa0hOh5SdGv6hiWHwtqrqQxOaisxOE8tzAIBYZS7NZaUmKs4Z+N2kA4XQ1LFDdaFZnivK9oWmPQcJTQCAGBXu80wmDu1tzzAMHQxRp6l/tm8QfO+hBnm8sXEFI6EJANBGpISmNDpN7dQ3eeRu3W6gZ5BDU+/MZCXEOdXk8aq8uiGo3ytcEJoAAG2E+x5NJja4bO9g63uXkRyv+Ljg/oqPczpUmJUsSdpz8GhQv1e4IDQBANqIlE7TsQ0um22uJHyYV84Fa1PLrxvQOte0O0bmmghNAIA2Ii400WnyO9i6VUROiK567G+Gphi5go7QBABoIxJDE0ep+PiHwEPVacrxDYPvZnkOABCLIm2mqdljqImjVCSFbmNLU/8Y23aA0AQAaCNSOk2JrjjFx/n2kWKJztxuoHWPrRAF3gGt2w7sOXhU3hjYdoDQBADwa2z2qLZ136NwD00Sc03HO9rkUWOzr+MWqk5Tn8xkuZwOuVu82l/bGJLvaSdCEwDAz+wyJbicSm9d/gpn/vPn2ODSP88Uiu0GTK44p/r29G07sLsq+ueaCE0AAL/j55kcjvA9QsWUlhQviU6TdGyPpqwQdZlMsTTXRGgCAPhFyjyTKa21G1bDXk0hOz7l68y5pli4go7QBADwi7TQlJ7s6zSxPCcd8m83ENr3jk4TACAmRVpoSvPPNNFpqmpdngt1p6kox9wVnE4TACCGRMoeTSY6Tcf492gK0caWpv7+bQfqo36TUUITAMAv4jpN5kxTQ2x3mo42taih2SNJyg7RESqmvj1T5HT4tjww//5EK0ITAMCvMuJCk6/TVN/kkScGNlc8EXOeKS3JpQRXaH+1J7ic6mNuOxDlS3SEJgCA34Ea3waFeelJNldiTUpCnJytOyPE8lxTqI9P+boB5sG9UT4MTmgCAEiSvF7D32nKS4+MTpPT4fB3m2J5rqmq9fiUUF85Zzp+rimaEZoAAJKkQ0eb1OI15HBIOREyCC4dm2uK5U7TobDpNLE8BwCIAftbl+ayUxNDdgxHIJidppoY7jQdtGmPJlOs7NUUOf8qAABBtd8/zxQ5XSZJ/jPyYrnTdNCmPZpM5q7ge6qORvW2A4QmAIAkaX+NOc8UGUPgpmPLc7HZaWps9qi+ydxuwJ7QVJiVIodDqnW3+K/ki0aEJgCApOM7TZEVmtL9y3Ox2Wkyr5zrkehSYnycLTUkxcepoPXvTTTPNRGaAACSju80RdbyXKx3mg7W27s0ZxqQE/1zTYQmAIAkqTJCO02xPgh+bAjc3tDUPwauoCM0AQAkSftrI3MQ3Ow0HXW3xOSu4ObyXFaIj0/5ugExsFdTWISmefPmqaioSElJSRo1apTWrl17wseWl5frpptu0pAhQ+R0OvXAAw+ErlAAiGLm8lxuWmR1mlITXXI6JENSnTv2uk3+K+fCpdNURWgKmiVLluiBBx7Qo48+qpKSEo0bN06TJ09WaWlph493u93q1auXHn30UY0YMSLE1QJAdGrxeFVVF5lXzzkdDvVIjN1tB8zluRy7O005vk4Ty3NBNGfOHN1xxx268847NWzYMM2dO1eFhYWaP39+h48fMGCAnn76ad16663KyMgIcbUAEJ0O1LllGJLL6bB9oLg70pNb55oaYqvT5G72+LtrWTa/b/2zfJ2m6oZmHY7SbQdsDU1NTU3auHGjJk2a1Ob+SZMmad26dQH7Pm63WzU1NW1uAIBjji3NJcppnoAbQdJaO02xtu2A2WVKSYhTcoI92w2YkhPilN/apdwVpXNNtoamqqoqeTwe5eXltbk/Ly9PFRUVAfs+s2fPVkZGhv9WWFgYsK8NANHA3KMpN8KW5kz+TlOMhqZw6Q6aS3TROgxu+/KcJDkcbf+rxjCMdvedilmzZqm6utp/27t3b8C+NgBEg8oIPULFlOFfnoux0NQ6hxYuByybB/fuqorOuSaXnd88JydHcXFx7bpKlZWV7bpPpyIxMVGJieHxFwoAwlGkHqFiMkNTdayFptZOU5bNV86Zon2DS1s7TQkJCRo1apRWrlzZ5v6VK1fqggsusKkqAIg9kXqEiindH5piaxDc3KMp2+Yr50wDonzbAVs7TZI0Y8YM3XLLLRo9erTGjh2rF154QaWlpZo+fbok39JaWVmZFi9e7H/O5s2bJUl1dXU6cOCANm/erISEBJ1xxhl2vAQAiHj7a48NgkeijKRjy3OGETsbXB4KkyNUTNG+7YDtoemGG27QwYMH9cQTT6i8vFzFxcVasWKF+vfvL8m3meXX92waOXKk//9v3LhRr7zyivr376/du3eHsnQAiBr7q32dpvyMyO40NXm8crd4ba4mNJpavP6jY+ze2NL09W0HeoZJmAuUsBgEv/vuu7V792653W5t3LhR48eP939u0aJFWr16dZvHG4bR7kZgAnAqunIywerVq+VwONrdvvjiixBWHFjl1Q2SInd5LsHlVHK875L7WJlrMg/qTY6PU0qC7T0QSdG/7UBYhCYAsFNXTyYwbd26VeXl5f7boEGDQlRxYNW7W/wdi4II7TRJsTcM7p9nCpMukymatx0gNAGIeV09mcCUm5ur/Px8/y0uzt7NBburvHVpLi3RpbTW2aBIlJ7cusFlzISm8NpuwFSUE73bDhCaAMS0UzmZYOTIkSooKNDEiRO1atWqTh8bzicTmEtzBZmR22WSYq/TVFUfnp0m8+BeOk0AEGW6czJBQUGBXnjhBS1dulTLli3TkCFDNHHiRK1Zs+aE3yecTyYoP+LrNBVkJNtcyalJj7XQZHaawmS7AVM0bzsQHpNjAGCzrpxMMGTIEA0ZMsT/8dixY7V371499dRTbS5kOd6sWbM0Y8YM/8c1NTVhE5z2mZ2mCJ5nko7bdiBGjlIJ95mmaNx2gE4TgJgWqJMJxowZo23btp3w84mJiUpPT29zCxfR0mk6dpRK9G9w2djsUZ3b9zrDbabp69sORBNCE4CYFqiTCUpKSlRQUBDo8kJiX5TMNMXS8px5fEpqQpyS4sPrAoTkhDh/1zLath1geQ5AzOvqyQRz587VgAEDNHz4cDU1Nenll1/W0qVLtXTpUjtfRrdVtF491ztKOk0NzR4dbWoJm72LgsG8ci47zLpMpv7ZKSqvbtSeg/U6p19Pu8sJmOj9GwUAFnX1ZIKmpibNnDlTZWVlSk5O1vDhw7V8+XJNmTLFrpdwSswtByK905TocirB5VRTi1cV1Y0a2KuH3SUFTVXrPFNOmM0zmYpyUvX+zkNRt+0AoQkA5DuZ4O677+7wc4sWLWrz8cMPP6yHH344BFUFX01js382JtIHwR0OhzKS4nWgzh31oSn8O03Rue0AM00AEMPMIfCM5PioWM4yl+jM7lm0qgrTjS1N0brtAKEJAGJYtGw3YEr3h6YGmysJLnMQPDtMD8SN1m0HCE0AEMP8Q+CZkT0EbspM8YWmsiPRG5oamjw62uSRFH57NJmiddsBQhMAxLDyI9HVacps7TR9dTh6Q5O5NJeW5FKiK7y2GzBF67YDhCYAiGH7zCvnoiQ09WxdrormTtPB+tYh8DA7PuXr+mf7luiiaRic0AQAMcx/WG+E79FkMjtN+440yDAMm6sJjnDfbsBUlONbooumbQcITQAQw/xHqET4Hk2mjJR4OSQ1Nnv9w9LRpirMtxswReO2A4QmAIhRhmH4r56L9N3ATS6nU2lJvq0TyqJ0rulghHSaonHbAUITAMSoqromNTZ75XBEz9VzkpSZEr1zTYZhHJtpCvNOk7k8t7OqPmqWSglNABCj9h72zZrkpycpwRU9vw782w5EYaepvsmjxmavpPDdo8nUPztFDodU29jin8OKdNHzrwQA0CXmZfmFPVNsriSwMpOjt9NkHp+SkRyv+Ljw/hWeFB/n/7u140CdzdUERnj/iQMAguar1k5T357RszQnHes0ReNeTQdqfaGpV5gvzZkG9mpdojsQHXNNhCYAiFF7D/lCRd+s6Oo09YziXcEPtHaaeqVFRmg6rfXQZDpNAICIFr2dptblucPRsz+Qyd9pIjTZgtAEADHKHJSOvtDk6zTVNLaoprHZ5moCK9JCE8tzAICI5/UaUTsInuiKi8or6Fo8Xh0+6rsKLVJmmsxO097DR9XY7LG5mlNHaAKAGHSgzq0mj1dxTkfUnDt3vD6t+05FU2g6WN8kryEluo5t4BnucnokKC3JJcOQdkfBzuCEJgCIQXsP+eZ9CjKS5ArzS9e7wx+aomgY/PilOYfDYXM11jgcDn+3KRqW6KLvXwoA4KS+itJ5JlPf1iVHMxxGA/+VcxGyNGcy55p2VEb+MDihCQBikHnlXLTNM5n6Z/te155oCk0RNgRu8neaouAMOkITAMQg/x5NUR6aSg8SmuwWTdsOEJoAIAZ9dSQ692gy9c/2LQntORQdh8UahhGxy3OnHbc8F+nvBaEJAGKQf7uBKNsN3NQnM1lxTocam72qbO3QRLKaxhY1tXjldEhZPcL7oN6v65+dKpfTofomj/ZVN9pdzikhNAFAjGnxeKN2Y0tTgsup3pm+rRT2RMESnbk0l5WaIJczsn51J7icKsrxdZu+3F9rczWnJrL+5AEAp6zsSINavIYSXU7lp0ffHk2mAa1LdNGwP9D+Gl+HJjctMt+vwXlpkqQvKwhNAIAIYl7FNCA7VU5nZOz30x39sqJnGLyy1hea8tIja57JNCjPNwz+5f7IHgYnNAFAjNlthqac6JxnMkVXp8m3PJcXoZ3BIa2dpm2VdJoAABHkWGhKtbmS4OpnbjsQ4Xs1GYbhX56L1NA0yAxN++vk9UbuFXSEJgCIMbtal6uKsqM7NPk7TRG+qeK+6ka5W7yKcziUHWFXzpkGZKcoIc6phmZPRB9tQ2gCgBhjhoiiaO80tc401TS26MjRJpur6T5zeDq7R+RdOWdyxTn9x6lE8hV0kfmnDwDolqYWr/8IlWgPTckJcf7B6UjedsAMGZG6NGcyl+gieRic0AQAMWTv4aPyGlJqQlzEHcfRHf2zIn8YfGuUhKbBub4r6LbRaQIARIJdB3zhoX92qhyO6N1uwOQ/uDcKOk35EbrdgMnfaYrgK+gITQAQQ8yOS7QvzZnMKwR3RegwuMdraFvrclakd5qG5PtC0/bKOnki9Ao6QhMAxJBdMbJHk+n01iWh7ZWROUdTeuio3C1euZwO9UyNzCvnTP2yUpQcH6fGZm/ELpcSmgAghpi/rAZE+XYDpuNDUyTuD7S19cq53LREOSN8OTXO6dDQAl+3acu+Gpur6R5CEwDEkN1VsXHlnKl/Vori4xxqaPZoX3Xk7Q+0ZV+1JCk/IzoOVj6jIF2StKWc0AQACGMNTceCQ7TvBm5yxTn9ATESl+jMcNE7M7LnmUzDzNBEpwkAEM62V9bJMKSs1ATl9IjsK7G6IpLnmj5rDRcF0dJp6k2nCQAQAcxL1we3njgfK07vFZmh6VB9k8qrfWfOFWRER6dpaH6aHA7pQK1blbWNdpfTZYQmAIgRx0JTms2VhNZpEdpp+qx1nmlAdoqS4uNsriYwUhJc/uXSz8sjb78mQhMAxIhYDU2Dcn2vd1tlnQwjcq6gM5fmhvfOsLmSwDojgueaCE0AECPMM79iLTQN7JUqh0OqbmhWVV3kHNxrhiZzDihaRPJcE6EJAGJAbWOzyo74rpyLtZmmpPg4Ffb0beYZSUt05vLc8GgLTf5OU7XNlXQdoQkAYsC21rCQm5aozJTI3lm6O/xX0B2IjNBU727x794edctzrSFwZ1W9jja12FxN1xCaACAGfNm6s7R5/lesGWSGpv2RMXz8RUWNDMMXcnulRdf2ELlpScpPT5JhSP8qi6wlOkITAMQAc57JHIqONWZYjJQrtj79KjqX5kwjCn3ds817D9tcSdcQmgAgBphXzg3Jj615JpO5xLWlvCYizqDbvPeIJOnswp72FhIkIwozJUkf742suSZCEwDEADM0DYqxK+dMp/VKVYLLqTp3i0oPHbW7nJMqaQ1NI/tl2lpHsJzdGprMcBgpCE0AEOUO1TepstYt6dhsT6xxxTk1tHWJ7rMw3x/oYJ1bew76gp3ZkYk2Z/bJkMMhlR1p0IHWv5uRgNAEAFHu0zLfEkhRTqrSkuJtrsY+5nzQZ2F+qbvZfTk9t4cykqPz/UpLivcH+I8jqNtEaAKAKPfpV0ck+f7rPpad0TrXFO6dppLSI5KOLWFFqxF9MyVF1hIdoQkAotwnrVdindU3tkPTsU5TeIemzVE+z2TyD4O3hvpIQGgCgChnLs/FeqdpWH66nA6pqs6typpGu8vpkMdrHAtNUXrlnOn4YfBIuKJRIjQBQFQ7UOtWeXWjHA5peIyHpuSEOA3s5ZujCddu044DdapztyglIS7qj7sZkp+mpHinahtbImandkITAESxf7V2mU7r1UM9El02V2O/cB8G37THt9njWX0z5IqL7l/R8XFOndPP1037YOdBm6uxJrrfEQCIcf55phjvMpmKW4fBzT+XcPPBrkOSpNH9s2yuJDTOL8qWJL3f+rrDHaEJAKLYp2VHJElnxvgQuMkcrt5UeliGEV5zNIZhaP0OX8dl7GnZNlcTGmMG+sLhBzsPhd370RFCEwBEMa6ca6u4T4YS4pyqqmvybyAZLnYfPKqKmkYlxDk1qn90D4GbRhRmKsHlVFWdWzur6u0u56QITQAQpcqrG1RZ65bTIZ1RQGiSpKT4OH+A/Gh3eC0JmV2ms/tlKik+zuZqQiMpPk4jW6+i+2BneL0fHSE0AUCUMn8JDe+doeSE2PglbMWoAb4uzsbWoetwsW5HlSRp7MDYWJoznd/6ej/YFf7D4IQmAIhS5i+h84tiY6jYKnPIekMYhSbDMPR+a8iNlXkm05iiyJlrIjQBQJQyO01jYqxzcTLmvND2yjodrm+yuRqf7ZV1qqpzK9HljPqdwL9uZL+eio9zqKKmUbvDbM7s6whNABCFKmsatbOqXg6HdC6dpjayUhN0Wq9USeGzRLe+dZ+i0QN6KtEVW0upyQlxOneA7+/oqi8qba6mc4QmAIhC5n4/w/LTlZEcb3M14Sfclujebg0LF56eY3Ml9rh0aK4kadVWQhMAIMT880wD6TJ1xOy+rQ+DnaiPNrVoXeuVc5cNy7O5GntcPMQXmj7YeUj17habqzkxQhMARCFzqNjccRltjRvk6+h88tURHaxz21rL2m1VamrxqjArWYNyo/u8uRM5rVeq+mWlqMnj1Xvbq+wu54QITQAQZarq3Npe6TsAlSvnOpaXnqRhBekyDF9osdM/P98vSZo4NE8Oh8PWWuzicDgiYomO0AQAUWbNlwckSWcUpKtnaoLN1YSvi4f0kiSttvGXtNdr6O0vfO9XrC7NmS4xQ9MXB8J26wFCEwBEmVVbfb+ELxnay+ZKwtvFg31/Pmu2VcnjteeX9MdfHVFVnVtpiS6dF+NdwfOLspQcH6eKmkb9q6zG7nI6RGgCgCjS4vH6O02XtA7XomPn9O+ptESXDtU36ZOvjthSw1tbfEtz4wf3UoIrtn8lJ8XH+Zfo3vhkn83VdCy23yEAiDIb9xxWdUOzMpLjNbJfbBz62l3xcU5d1DoQvrq1OxdKXq+hv5SUSZImn5kf8u8fjr5xdm9J0l8375PXpu5fZwhNABBF/v6vCknSxGG5inPG5lBxV5jduJWtHZ9Qen/XQe2rblRakivm55lMFw/ppbQklypqGvVhmB2oLBGaACBqGIahf3zmC02TiwtsriYyXH5GnlxOh7aU1+jL/bUh/d6vbfJ1ma46q0BJ8bG1C/iJJLriNLnY13X7y+bwW6IjNAFAlCjZe0Tl1Y1KTYjz70OEzvVMTfBvrPh661JZKDQ0ebTi03JJ0rXn9A3Z940EU8/uI0la8Wm5mlq8NlfTFqEJAKKE2bmYNDyfzkUXTBvp+yX9lxDO0by1pUL1TR4VZiVrdH9mz443ZmC2ctMSVd3QbMuyaWcITQAQBZpavP4rjswQAGsmDstVWqJLZUcaQjZH8/L7eyRJ147sG7MbWp5InNOhG84tlCT9ft1ue4v5GkITAESB//t8v44cbVZuWmLMHvraXUnxcZpypm8GzOzWBdOm0sP6aPdhxcc59O3z+wX9+0Wim8f0l8vp0Ie7D+lfZdV2l+NHaAKAKGB2Lq4fXchVc90w7Rxfd+6NT/apuqE5qN/rd2t3SpKuObuPctOTgvq9IlVeepI/yC58b7e9xRyH0AQAEW57Za3W7Tgop0O6kc5Ft5xflKUheWk62uTR//tob9C+z56D9XqzdVuI744fGLTvEw1uv3CAJOmNj/fpQK29hyqbCE0AEOFeWOPrXFw2LE99MpNtriYyORwO/y/phe/tCtpVW8+v2Smv4duPaHBeWlC+R7Q4p19PnV2YqSaPV/NX77C7HEmEJgCIaBXVjXqt9VL5uyacZnM1kW3ayD7KTUvUvupG/XnjVwH/+jsP1GlJaxfre7xXlsy4fLAk6Q/v79beQ0dtrobQBAAR7Zm3t6nZY+i8oiyN4tL1U5IUH6fprWHmN29vU0OTJ6Bff/bfv5DHa+jSobk6f2B2QL92tBo3KEcXnp6tZo+hOSu/tLscQhMARKrjOxf/ecUQm6uJDjed3099MpNVXt3oX/YMhH9+vl8rt+yXy+nQI5OHBuzrRjuHw6FHrhwmSXp9c5ltByubCE0AEIEMw9CP/vqZv3Nx7oAsu0uKCknxcf5QM2/1du04UHfKX/PI0Sb98LVPJUl3XFTELFMXndk3Q9ec3VuGIc3808dqbA5sB7ArCE0AEIH++vE+rd1WpQSXU49ddYbd5USVq84q0PjBveRu8Wrmnz4+paFwr9fQf/75E+2vcWtgr1Q9cNngAFYaO/77qjOU0yNRX+6vs3WZjtAEABGm7EiD/vv1f0mS7r3kdA3ISbW5oujicDg0+9ozlZbkUknpET25fEu3v9b//N+XWrllvxLinHr6hpFKTuB4m+7I7pGon197piTpt2t36p0vD9hSB6EJACJIY7NH0/+wUTWNLRpRmKm7L+YqrGDok5msOdefLUn6/fo9mrd6e5e/xvPv7NBv3vY976fTinVm34xAlhhzLjsjT986t1CGId398kZbdgonNAFAhGj2eHXvK5v0aVm1eqbE69kbR8oVx4/xYLn8jDz9cIpvvumXb27V7L9/rhbPyZfqmlq8evyvn2n237+QJD10+WBdP7owqLXGih9PHa6xA7NV3+TRvy/6SLuq6kP6/fnXBgARoKHJo++9vEn/93mlEl1OLbh5lAqzUuwuK+r9x/jT/FcmPv/OTv1/z6/Xx3uPdPhYwzD03vYqfePZd7Wo9aDZH1w5VPdeenqIqo1+ia44PX/rKA3NT9OBWremzXtPH+4KzSHLkuQK2XcCAHTL1opa3fvKJm2rrPMHJvb5CZ17LjldhVkp+uGyT1VSekRTn3tPZxdmavygHPXNSpHXa2hXVb1Wba3Ul/t9V9tlpsTrl9edpUnD822uPvqkJ8XrD3ecrzt//5E+/qpa3/7d+/rPK4boOxcWBb3zGhadpnnz5qmoqEhJSUkaNWqU1q5d2+nj33nnHY0aNUpJSUkaOHCgFixYEKJKAUSrcP059McPS/WNZ9/Vtso69UpL1B/uOF+XDM0NyvfCiX1jRG+tnDFe157TRy6nQ5v3HtEzb2/Xw3/+RI8s+1TPr9mpL/f7Qu1tY/vr7YcuJjAFUa+0RP3xP8Zqypn5avYY+tmKL3TNvPf03vYqGYYRtO9re6dpyZIleuCBBzRv3jxdeOGFev755zV58mRt2bJF/fq1P3hy165dmjJlir773e/q5Zdf1nvvvae7775bvXr10nXXXWfDKwAQ6cL559CRhma5W7yaMLiXfn39COX0SAzo14d1BRm+4fBHrhyqlZ/v18d7j6iy1i2nw6HemUka1b+nLh2ap4zkeLtLjQnJCXF67qZz9KcNX+mny7foX2U1+vbvPtCZfTJ0zyWn6crigoB/T4cRzEhmwfnnn69zzjlH8+fP9983bNgwXXPNNZo9e3a7x//gBz/QX//6V33++ef++6ZPn66PP/5Y69evt/Q9a2pqlJGRoerqaqWnp5/6izjOqP9cHNCvB9hp469uDfjXDOa/v+4K559DXq+h5Z+W69/OLJDT6ejCqwpPr3xQGpLvc9P57cNuJOPPrXOVtY167u3tWrJhrxqbvbrqrAI9e9M5nT6nOz+LbF2ea2pq0saNGzVp0qQ290+aNEnr1q3r8Dnr169v9/grrrhCGzZsUHNzc9BqBRCdwv3nkNPp0NUjekdFYAKCJTctST+eWqz3fnCp7p84SN8L0lYcti7PVVVVyePxKC8vr839eXl5qqio6PA5FRUVHT6+paVFVVVVKiho345zu91yu93+j6urfXs71NTUnOpLaMfjbgj41wTsEox/I+bXtLnJ7ReNP4fC2dH62pB8n2j7c+XPzZp4SXee75slO9lr6c7PIttnmiTf7qvHMwyj3X0ne3xH95tmz56tH//4x+3uLyxk3wygMxm/mR60r11bW6uMjPDZ7I+fQ9Hlu3YXEKFi8c+tKz+LbA1NOTk5iouLa/dfc5WVle3+K86Un5/f4eNdLpeyszu+BHfWrFmaMWOG/2Ov16tDhw4pOzu70x+KCE81NTUqLCzU3r17w2YmBtYZhqHa2lr17t3b7lIkhffPoWj7u87rCW/R9HqsvJbu/CyyNTQlJCRo1KhRWrlypaZNm+a/f+XKlZo6dWqHzxk7dqzeeOONNve99dZbGj16tOLjO75iITExUYmJba84yczMPLXiYbv09PSI/4cdq8KpwxQJP4ei7e86rye8RdPrOdlr6erPItv3aZoxY4Z+97vf6aWXXtLnn3+uBx98UKWlpZo+3bcsMGvWLN1667EreKZPn649e/ZoxowZ+vzzz/XSSy/pxRdf1MyZM+16CQAiHD+HAFhh+0zTDTfcoIMHD+qJJ55QeXm5iouLtWLFCvXv31+SVF5ertLSY5daFhUVacWKFXrwwQf13HPPqXfv3nrmmWfYowlAt/FzCIAVtu/TBHSV2+3W7NmzNWvWrHbLHUA0iba/67ye8BZNrydYr4XQBAAAYIHtM00AAACRgNAEAABgAaEJAADAAkITAISJw4cP65ZbblFGRoYyMjJ0yy236MiRI50+5/bbb5fD4WhzGzNmTGgK/pp58+apqKhISUlJGjVqlNauXdvp49955x2NGjVKSUlJGjhwoBYsWBCiSq3pyutZvXp1u/fB4XDoiy++CGHFJ7ZmzRpdffXV6t27txwOh15//fWTPiec35+uvp5AvT+EJgAIEzfddJM2b96sN998U2+++aY2b96sW2655aTPu/LKK1VeXu6/rVixIgTVtrVkyRI98MADevTRR1VSUqJx48Zp8uTJbbZqON6uXbs0ZcoUjRs3TiUlJfrhD3+o+++/X0uXLg1x5R3r6usxbd26tc17MWjQoBBV3Ln6+nqNGDFCzz77rKXHh/v709XXYzrl98cAANhuy5YthiTj/fff99+3fv16Q5LxxRdfnPB5t912mzF16tQQVNi58847z5g+fXqb+4YOHWo88sgjHT7+4YcfNoYOHdrmvrvuussYM2ZM0Grsiq6+nlWrVhmSjMOHD4egulMjyXjttdc6fUy4vz/Hs/J6AvX+0GkCgDCwfv16ZWRk6Pzzz/ffN2bMGGVkZGjdunWdPnf16tXKzc3V4MGD9d3vfleVlZXBLreNpqYmbdy4UZMmTWpz/6RJk05Y+/r169s9/oorrtCGDRvU3NwctFqt6M7rMY0cOVIFBQWaOHGiVq1aFcwygyqc359TcarvD6EJAMJARUWFcnNz292fm5vb7nDg402ePFn/+7//q7ffflu//vWv9dFHH+nSSy+V2+0OZrltVFVVyePxtDvgOC8v74S1V1RUdPj4lpYWVVVVBa1WK7rzegoKCvTCCy9o6dKlWrZsmYYMGaKJEydqzZo1oSg54ML5/emOQL0/th+jAgDR7PHHH9ePf/zjTh/z0UcfSZIcDke7zxmG0eH9phtuuMH//4uLizV69Gj1799fy5cv17XXXtvNqrvn63WerPaOHt/R/XbpyusZMmSIhgwZ4v947Nix2rt3r5566imNHz8+qHUGS7i/P10RqPeH0AQAQXTvvffqW9/6VqePGTBggD755BPt37+/3ecOHDjQ7r/4O1NQUKD+/ftr27ZtXa61u3JychQXF9euC1NZWXnC2vPz8zt8vMvlUnZ2dtBqtaI7r6cjY8aM0csvvxzo8kIinN+fQOnO+0NoAoAgysnJUU5OzkkfN3bsWFVXV+vDDz/UeeedJ0n64IMPVF1drQsuuMDy9zt48KD27t2rgoKCbtfcVQkJCRo1apRWrlypadOm+e9fuXKlpk6d2uFzxo4dqzfeeKPNfW+99ZZGjx6t+Pj4oNZ7Mt15PR0pKSkJ6fsQSOH8/gRKt96fUxojBwAEzJVXXmmcddZZxvr1643169cbZ555pnHVVVe1ecyQIUOMZcuWGYZhGLW1tcZDDz1krFu3zti1a5exatUqY+zYsUafPn2MmpqakNb+xz/+0YiPjzdefPFFY8uWLcYDDzxgpKamGrt37zYMwzAeeeQR45ZbbvE/fufOnUZKSorx4IMPGlu2bDFefPFFIz4+3vjzn/8c0rpPpKuv53/+53+M1157zfjyyy+Nf/3rX8YjjzxiSDKWLl1q10too7a21igpKTFKSkoMScacOXOMkpISY8+ePYZhRN7709XXE6j3h9AEAGHi4MGDxre//W0jLS3NSEtLM7797W+3u0RakrFw4ULDMAzj6NGjxqRJk4xevXoZ8fHxRr9+/YzbbrvNKC0tDX3xhmE899xzRv/+/Y2EhATjnHPOMd555x3/52677TZjwoQJbR6/evVqY+TIkUZCQoIxYMAAY/78+SGuuHNdeT2/+MUvjNNOO81ISkoyevbsaVx00UXG8uXLbai6Y+Yl91+/3XbbbYZhRN7709XXE6j3x2EYrZNdAAAAOCG2HAAAALCA0AQAAGABoQkAAMACQhMAAIAFhCYAAAALCE0AAAAWEJoAAAAsIDQBAABYQGgCAKAThmHoP/7jP5SVlSWHw6HNmzfbVsvu3bttryGWEZoAAGHr9ttvl8PhkMPhUHx8vAYOHKiZM2eqvr7eHyDMW0ZGhsaMGdPuoNlFixa1eZx5S0pKslTDm2++qUWLFulvf/ubysvLVVxcHIyXigjgsrsAAAA6c+WVV2rhwoVqbm7W2rVrdeedd6q+vl4/+MEPJEn/93//p+HDh+vIkSOaN2+errvuOm3atKlNuElPT9fWrVvbfF2Hw2Hp++/YsUMFBQW64IILAveiEJHoNAEAwlpiYqLy8/NVWFiom266Sd/+9rf1+uuv+z+fnZ2t/Px8DR06VE8++aSam5u1atWqNl/D4XAoPz+/zS0vL++k3/v222/Xfffdp9LSUjkcDg0YMECSb8nul7/8pQYOHKjk5GSNGDFCf/7zn/3PW716tRwOh/7xj39o5MiRSk5O1qWXXqrKykr9/e9/17Bhw5Senq4bb7xRR48e9T/vzTff1EUXXaTMzExlZ2frqquu0o4dOzqtccuWLZoyZYp69OihvLw83XLLLaqqqrLwJ4uuIjQBACJKcnKympub293f3Nys3/72t5Kk+Pj4gHyvp59+Wk888YT69u2r8vJyffTRR5Kk//qv/9LChQs1f/58ffbZZ3rwwQd1880365133mnz/Mcff1zPPvus1q1bp7179+r666/X3Llz9corr2j58uVauXKlfvOb3/gfX19frxkzZuijjz7SP//5TzmdTk2bNk1er7fD+srLyzVhwgSdffbZ2rBhg958803t379f119/fUBeP9pieQ4AEDE+/PBDvfLKK5o4caL/vgsuuEBOp1MNDQ3yer0aMGBAu9BQXV2tHj16tLnvggsu0FtvvdXp98vIyFBaWpri4uKUn58vyRds5syZo7fffltjx46VJA0cOFDvvvuunn/+eU2YMMH//J/+9Ke68MILJUl33HGHZs2apR07dmjgwIGSpG9+85tatWqVf6nxuuuua/P9X3zxReXm5mrLli0dzlLNnz9f55xzjn72s5/573vppZdUWFioL7/8UoMHD+709aFrCE0AgLD2t7/9TT169FBLS4uam5s1depU/eY3v/Evay1ZskRDhw7Vl19+qQceeEALFixQVlZWm6+RlpamTZs2tbkvOTm5W/Vs2bJFjY2Nuvzyy9vc39TUpJEjR7a576yzzvL//7y8PKWkpPgDk3nfhx9+6P94x44d+u///m+9//77qqqq8neYSktLOwxNGzdu1KpVq9oFQvNrEZoCi9AEAAhrl1xyiebPn6/4+Hj17t3bv/S2e/duSVJhYaEGDRqkQYMGqUePHrruuuu0ZcsW5ebm+r+G0+nU6aefHpB6zCCzfPly9enTp83nEhMT23x8/DKheQXg8RwOR5ult6uvvlqFhYX67W9/q969e8vr9aq4uFhNTU0nrOXqq6/WL37xi3afKygo6NoLw0kRmgAAYS01NdVy4JkwYYKKi4v15JNP6umnnw5KPWeccYYSExNVWlraZinuVB08eFCff/65nn/+eY0bN06S9O6773b6nHPOOUdLly7VgAED5HLxKz3YGAQHAESVhx56SM8//7zKysr89xmGoYqKina3Ew1YdyYtLU0zZ87Ugw8+qN///vfasWOHSkpK9Nxzz+n3v/99t+vu2bOnsrOz9cILL2j79u16++23NWPGjE6fc8899+jQoUO68cYb9eGHH2rnzp1666239J3vfEcej6fbtaBjhCYAQFS56qqrNGDAAD355JP++2pqalRQUNDuVllZ2a3v8ZOf/ESPPfaYZs+erWHDhumKK67QG2+8oaKiom7X7XQ69cc//lEbN25UcXGxHnzwQf3qV7/q9Dm9e/fWe++9J4/HoyuuuELFxcX6/ve/r4yMDDmd/IoPNIdhGIbdRQAAAIQ7YigAAIAFhCYAQMwqLS1Vjx49TngrLS21u0SEEZbnAAAxq6Wlxb91QUe4Kg3HIzQBAABYwPIcAACABYQmAAAACwhNAAAAFhCaAAAALCA0AQAAWEBoAgAAsIDQBAAAYAGhCQAAwIL/H+3s9nilnDmqAAAAAElFTkSuQmCC",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# we can also visualize the proportion, by plotting a bargraph and a histogram of the data\n",
"\n",
"fig, ax = plt.subplots(1,2, figsize=(6,6))\n",
"sns.barplot(y=\"PRE_female\", data=df, ax=ax[0])\n",
"sns.distplot(df[\"PRE_female\"], ax=ax[1])\n",
"plt.tight_layout()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "q1V5IndN2gVP",
"outputId": "c3e39663-fc83-4f05-b998-d8b5b5bb49c0"
},
"outputs": [
{
"data": {
"text/plain": [
"0 1\n",
"1 0\n",
"2 1\n",
"3 1\n",
"4 0\n",
" ..\n",
"255 0\n",
"256 1\n",
"257 1\n",
"258 0\n",
"259 0\n",
"Name: PRE_female, Length: 260, dtype: int64"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# the variable \"PRE_female\" is a Series of 1s and 0s\n",
"\n",
"df[\"PRE_female\"]"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "9pURyyhXfnK8",
"outputId": "36856905-7b95-4806-be20-002ddc8b353e"
},
"outputs": [
{
"data": {
"text/plain": [
"numpy.int64"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"type(df[\"PRE_female\"][0])"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "ECNW3vdGVMA2",
"outputId": "36bf88a8-f382-47a7-8a40-43b7d8b6c13f"
},
"outputs": [
{
"data": {
"text/plain": [
"{'Unnamed: 0': dtype('int64'),\n",
" 'PRE_male': dtype('int64'),\n",
" 'PRE_female': dtype('int64'),\n",
" 'POST_male': dtype('int64'),\n",
" 'POST_female': dtype('int64'),\n",
" 'POST_male_salary': dtype('int64'),\n",
" 'POST_female_salary': dtype('int64'),\n",
" 'POST_male_friendly': dtype('float64'),\n",
" 'POST_female_friendly': dtype('float64'),\n",
" 'POST_male_intelligent': dtype('float64'),\n",
" 'POST_female_intelligent': dtype('float64'),\n",
" 'Hire_m': dtype('int64'),\n",
" 'Hire_f': dtype('int64'),\n",
" 'Choice_m': dtype('int64'),\n",
" 'Choice_f': dtype('int64'),\n",
" 'itemnum': dtype('int64'),\n",
" 'partnum': dtype('int64'),\n",
" 'Gender': dtype('O'),\n",
" 'Ide': dtype('int64'),\n",
" 'Political': dtype('O'),\n",
" 'Edu': dtype('O')}"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# this is how you see what type of variables are in your data (in each column)\n",
"\n",
"df.dtypes.to_dict()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Y2E8XjuY2Eiy",
"outputId": "f6ac7ea6-ccf8-41c4-a52d-b205d6a57ca2"
},
"outputs": [
{
"data": {
"text/plain": [
"array([159, 101])"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# the chisquare test expects as input the observed frequencies in each category\n",
"# to compute this we can use numpy's function bincount, which counts how many 0s and how many 1s are in the \"PRE_female\" Series\n",
"\n",
"np.bincount(df[\"PRE_female\"])"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "vVRAAYtBTuPh",
"outputId": "d4b150fd-6816-43ec-e42b-9db9933fba05"
},
"outputs": [
{
"data": {
"text/plain": [
"0 159\n",
"1 101\n",
"Name: PRE_female, dtype: int64"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# the chisquare test expects as input the observed frequencies in each category\n",
"# to compute this we can use panda's function value_counts, which is the same as bincounts but can handle strings not just numebrs\n",
"\n",
"df[\"PRE_female\"].value_counts()"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "5je5-cDi1zZ6",
"outputId": "b7a279bd-1f60-497d-d68a-0609c13ca6a6"
},
"outputs": [
{
"data": {
"text/plain": [
"Power_divergenceResult(statistic=12.938461538461539, pvalue=0.00032189944909632355)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# now we can run the chisquare test on the count of the categories computed above\n",
"# by default, chisquare compares the counts to the Null in which the categories are equally likely\n",
"# so in this case, it tests whether the 159 0s and 101 1s are significantly different from 130 0s and 130 1s, or in other words, a 50% incidence rate\n",
"# it reports the chi square statistic and the p value\n",
"\n",
"chisquare(np.bincount(df[\"PRE_female\"]))"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "B1xYG_GnUR6a",
"outputId": "c0e0b790-a706-406a-c138-7141a6eab14d"
},
"outputs": [
{
"data": {
"text/plain": [
"Power_divergenceResult(statistic=12.938461538461539, pvalue=0.00032189944909632355)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# now we can run the chisquare test on the count of the categories computed above\n",
"# by default, chisquare compares the counts to the Null in which the categories are equally likely\n",
"# so in this case, it tests whether the 159 0s and 101 1s are significantly different from 130 0s and 130 1s, or in other words, a 50% incidence rate\n",
"# it reports the chi square statistic and the p value\n",
"\n",
"chisquare(df[\"PRE_female\"].value_counts())"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "dnzyE1-h7cvz",
"outputId": "a55eb2dc-30b3-44a9-dbd8-aa9b92e566d2"
},
"outputs": [
{
"data": {
"text/plain": [
"0.2230769230769231"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# if you need the effect size of the test you can use the following formula:\n",
"\n",
"res = chisquare(np.bincount(df[\"PRE_female\"]))\n",
"np.sqrt(res.statistic/len(df))"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "sk2lcr9tyMbw"
},
"source": [
"Reporting result:\n",
"\n",
"\"At pretest in Condition 1, participants chose women (Mean=38%) significantlty less than expected by chance (50%), χ2=12.93, w=0.22, P<0.001.\""
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "7uvo14rL6xEN"
},
"source": [
"## Independent proportions chi square"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "E4nLRkNKQGRY"
},
"source": [
"Chi square test of independence tests if two proportions are different from each other"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "lCrOF7Jq8uub"
},
"source": [
"The variable \"PRE_female\" encodes whether participants chose a man (PRE_female==0) or a woman (PRE_female==1) at pretest, in the female condition of this experiment.\n",
"\n",
"The variable \"PRE_female\" encodes whether participants chose a man (PRE_male==0) or a woman (PRE_male==1) at pretest, in the male condition of this experiment.\n",
"\n",
"To test whether the proportion of men/women choices are different in the male versus female conditions (comparing the PRE_female proportion to the PRE_male proportion) we can run a (two way) chi square test (or a chi square test of independence)."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 441
},
"id": "csdn0PjU8NSt",
"outputId": "3cf23cfe-9a96-45ea-a421-18577d5e4f54"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAJOCAYAAACz06ChAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAgWklEQVR4nO3de3BU9d3H8c9CSEKQhEogQo0hcQRhIsokRRMaLCpRvLTjjagFtIIljcgl8tCmcSgwQhwvNFCFGCsgTrWphbGjZsTMFGpi6B/EMFjjMKWAoWExF2oSGEkk2ecPHrbPlwTJZbNns3m/ZnaGPTkbvmkPvnP2/HbX5fF4PAIA4P8McnoAAEBgIQwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwQpwewAnt7e06fvy4hg8fLpfL5fQ4CGAej0fNzc0aO3asBg0KvN+jOJbRVd05lgdkGI4fP67Y2Finx0A/cuzYMV155ZVOj9EBxzK6qyvH8oAMw/DhwyWd+x8oMjLS4WkQyJqamhQbG+s9ZgINxzK6qjvH8oAMw/lT7sjISP4xoUsC9WkajmV0V1eO5cB70hQA4CjCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMEKcHgAAfGXJkiWqq6uTJI0aNUobNmxweKL+iTAACBp1dXX66quvnB6j3+OpJACAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGCFODwDnLVmyRHV1dZKkUaNGacOGDQ5PBMBJhAGqq6vTV1995fQYAAIETyUBAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwAiIMmzZtUnx8vMLDw5WUlKTS0tIuPe6TTz5RSEiIbrjhhr4dEAAGEMfDUFRUpKVLlyo3N1eVlZVKS0vTrFmzVF1d/Z2Pa2xs1Lx583Trrbf6aVIAGBgcD8P69es1f/58LViwQBMnTlR+fr5iY2O1efPm73zcwoUL9cgjjyglJcVPkwLAwOBoGFpbW1VRUaH09HSzPT09XeXl5Rd93NatW/Wvf/1Lv/nNb7r097S0tKipqcncgP6IYxn+4GgY6uvr1dbWppiYGLM9JiZGJ06c6PQx//znP/WrX/1Kf/jDHxQS0rVPJs3Ly1NUVJT3Fhsb2+vZASdwLMMfHH8qSZJcLpe57/F4OmyTpLa2Nj3yyCNavXq1xo8f3+Xvn5OTo8bGRu/t2LFjvZ4ZcALHMvyha79y95Ho6GgNHjy4w9lBbW1th7MISWpubta+fftUWVmpRYsWSZLa29vl8XgUEhKijz76SLfcckuHx4WFhSksLKxvfgjAjziW4Q+OnjGEhoYqKSlJJSUlZntJSYlSU1M77B8ZGanPPvtM+/fv994yMzM1YcIE7d+/XzfeeKO/RgeAoOXoGYMkZWdna+7cuUpOTlZKSooKCwtVXV2tzMxMSedOnWtqarR9+3YNGjRIiYmJ5vGjR49WeHh4h+0AgJ5xPAwZGRlqaGjQmjVr5Ha7lZiYqOLiYsXFxUmS3G73JV/TAADwHcfDIElZWVnKysrq9Gvbtm37zseuWrVKq1at8v1QADBABcSqJABA4CAMAACDMAAADMIAADAIAwDAIAwAAIMwAACMgHgdAwD/Svqf7U6P0Cci/3PK+9uu+z+ngvbnrHhhXp9+f84YAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgBHi9AD9SdL/bHd6hD4R+Z9T3t8Q3P85FbQ/Z8UL85weAegXOGMAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAEaI0wMAgK+0DxnW6Z/RPYQBQNA4NWGW0yMEBZ5KAgAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABgBEYZNmzYpPj5e4eHhSkpKUmlp6UX3LSsr07Rp0zRy5EgNHTpU1157rX7729/6cVoACG4hTg9QVFSkpUuXatOmTZo2bZpeffVVzZo1S1VVVbrqqqs67D9s2DAtWrRIkydP1rBhw1RWVqaFCxdq2LBh+vnPf+7AT9D/tQ8Z1umfAQxMjodh/fr1mj9/vhYsWCBJys/P165du7R582bl5eV12H/KlCmaMmWK9/64ceO0c+dOlZaWEoYeOjVhltMjAAggjj6V1NraqoqKCqWnp5vt6enpKi8v79L3qKysVHl5uW6++ea+GBEABhxHzxjq6+vV1tammJgYsz0mJkYnTpz4zsdeeeWVqqur09mzZ7Vq1SrvGUdnWlpa1NLS4r3f1NTUu8EBh3Aswx8C4uKzy+Uy9z0eT4dtFyotLdW+fftUUFCg/Px8vf322xfdNy8vT1FRUd5bbGysT+YG/I1jGf7gaBiio6M1ePDgDmcHtbW1Hc4iLhQfH6/rrrtOTzzxhJYtW6ZVq1ZddN+cnBw1NjZ6b8eOHfPF+IDfcSzDHxwNQ2hoqJKSklRSUmK2l5SUKDU1tcvfx+PxmNPrC4WFhSkyMtLcgP6IYxn+4PiqpOzsbM2dO1fJyclKSUlRYWGhqqurlZmZKencb0g1NTXavn27JOmVV17RVVddpWuvvVbSudc1vPjii3rqqacc+xkAIJg4HoaMjAw1NDRozZo1crvdSkxMVHFxseLi4iRJbrdb1dXV3v3b29uVk5OjI0eOKCQkRFdffbWee+45LVy40KkfAQCCiuNhkKSsrCxlZWV1+rVt27aZ+0899RRnBwDQhwJiVRIAIHAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAccOjQIe3atUvffPONJMnj8Tg8EfBfhAHwo4aGBt12220aP3687rzzTrndbknSggUL9PTTTzs8HXAOYQD8aNmyZQoJCVF1dbUiIiK82zMyMvThhx86OBnwXyFODwAMJB999JF27dqlK6+80my/5ppr9OWXXzo0FWBxxgD40enTp82Zwnn19fUKCwtzYCKgI8IA+NH06dO1fft2732Xy6X29na98MILmjFjhoOTAf/FU0mAH73wwgv60Y9+pH379qm1tVUrVqzQ559/rpMnT+qTTz5xejxAUi/PGEpLSzVnzhylpKSopqZGkvTmm2+qrKzMJ8MBwWbSpEk6cOCApk6dqpkzZ+r06dO67777VFlZqauvvtrp8QBJvThj2LFjh+bOnauf/vSnqqysVEtLiySpublZ69atU3Fxsc+GBILJFVdcodWrVzs9BnBRPQ7Ds88+q4KCAs2bN09//OMfvdtTU1O1Zs0anwwHBIMDBw50ed/Jkyf34SRA1/Q4DAcPHtT06dM7bI+MjNTXX3/dm5mAoHLDDTfI5XJd8tXNLpdLbW1tfpoKuLgeh2HMmDE6dOiQxo0bZ7aXlZUpISGht3MBQePIkSNOjwB0S4/DsHDhQi1ZskRbtmyRy+XS8ePHtXfvXi1fvlwrV6705YxAvxYXF+f0CEC39DgMK1asUGNjo2bMmKEzZ85o+vTpCgsL0/Lly7Vo0SJfzggEnaqqKlVXV6u1tdVs//GPf+zQRMB/9ep1DGvXrlVubq6qqqrU3t6uSZMm6bLLLvPVbEDQOXz4sO6991599tln5rqDy+WSJK4xICD0+pXPERERSk5O1tSpU4kCcAlLlixRfHy8vvrqK0VEROjzzz/Xxx9/rOTkZO3Zs8fp8QBJ3TxjuO+++7q8786dO7s9DBDs9u7dq7/+9a8aNWqUBg0apEGDBumHP/yh8vLytHjxYlVWVjo9ItC9MERFRfXVHMCA0NbW5j2zjo6O1vHjxzVhwgTFxcXp4MGDDk8HnNOtMGzdurWv5gAGhMTERB04cEAJCQm68cYb9fzzzys0NFSFhYUs80bA4E30AD965plndPr0aUnn3j3g7rvvVlpamkaOHKmioiKHpwPO6VUY/vznP+tPf/pTp8vuPv30014NBgSj22+/3fvnhIQEVVVV6eTJk/re977nXZkEOK3Hq5I2btyon/3sZxo9erQqKys1depUjRw5UocPH9asWbN8OSMQ1C6//HKigIDS4zOGTZs2qbCwUA8//LDeeOMNrVixQgkJCVq5cqVOnjzpyxmBoHHmzBn97ne/0+7du1VbW6v29nbzdc60EQh6HIbq6mqlpqZKkoYOHarm5mZJ0ty5c3XTTTfp5Zdf9s2EQBB5/PHHVVJSogceeEBTp07lTAEBqcdhuOKKK9TQ0KC4uDjFxcXp73//u66//nodOXLkku8iCQxUH3zwgYqLizVt2jSnRwEuqsfXGG655Ra99957kqT58+dr2bJlmjlzpjIyMnTvvff6bEAgmHz/+9/X8OHDnR4D+E49PmMoLCz0Pj+amZmpyy+/XGVlZbrnnnuUmZnpswGBYPLSSy/pl7/8pQoKCnjXVQSsHofh/Mv5z5s9e7Zmz57tk6GAYJWcnKwzZ84oISFBERERGjJkiPk6CzcQCHr1OoYzZ87owIEDna6u4O2DgY4efvhh1dTUaN26dYqJieHiMwJSj8Pw4Ycfat68eaqvr+/wNT6iEOhceXm59u7dq+uvv97pUYCL6vHF50WLFunBBx+U2+1We3u7uREFoHPXXnutvvnmG6fHAL5Tj8NQW1ur7OxsxcTE+HIeIKg999xzevrpp7Vnzx41NDSoqanJ3IBA0OOnkh544AHt2bNHV199tS/nAYLaHXfcIUm69dZbzXaPx8NTsAgYPQ7Dyy+/rAcffFClpaW67rrrOqyuWLx4ca+HA4LN7t27nR4BuKQeh+Gtt97Srl27NHToUO3Zs8esrnC5XIQB6MTNN9/s9AjAJfX4GsMzzzyjNWvWqLGxUUePHtWRI0e8t8OHD/tyRiColJaWas6cOUpNTVVNTY0k6c0331RZWZnDkwHn9DgMra2tysjIMC9yA/DdduzYodtvv11Dhw7Vp59+qpaWFklSc3Oz1q1b5/B0wDk9/q/6o48+yidOAd307LPPqqCgQK+99pq5LpeamspbbiNg9PgaQ1tbm55//nnt2rVLkydP7nDxef369b0eDgg2Bw8e1PTp0ztsj4yM1Ndff+3/gYBO9DgMn332maZMmSJJ+sc//mG+xsv8gc6NGTNGhw4d0rhx48z2srIyJSQkODMUcIEeh6Gry+7+/e9/a+zYsVyLACQtXLhQS5Ys0ZYtW+RyuXT8+HHt3btXy5cv18qVK50eD5DUyzfR64pJkyZp//79/DYESFqxYoUaGxs1Y8YMnTlzRtOnT1dYWJiWL1+uRYsWOT0eIMkPYeDT3DDQHThwQImJid6z5rVr1yo3N1dVVVVqb2/XpEmTdNlllzk8JfBfPL8D9LEpU6Z434U4ISFBDQ0NioiIUHJysqZOnUoUEHAIA9DHRowYoSNHjkiSjh492uGzS4BA0+dPJQED3f3336+bb75ZY8aMkcvlUnJysgYPHtzpvrxrAAJBn4eBpasY6AoLC3Xffffp0KFDWrx4sZ544gkNHz7c6bGAi+LiM+AH599uu6KiQkuWLLlkGFjmDSf5/KjzeDyqra313q+qqlJcXJyv/xqgX9q6dWuXzhYmTZqko0eP9v1AQCe6HYaIiAjV1dV5799xxx1yu93e+7W1tRozZoz3fmxs7EWfTwXQOc604aRuh+HMmTPmoP3kk086fIYtBzUA9F998gQmF5wBoP/iyhYAwOh2GFwuV4eP8eQMAfAt/k3BSd1erurxeDR+/HjvgXvq1ClNmTLFu6yO6wtA7/HvCE7qdhi2bt3aF3MAA5rH41FdXZ1Gjx4t6dwy77Fjxzo8FQaqbofh0Ucf7Ys5gKAWERGhL7/8UqNGjZJ0bpn31q1bvUu7a2trNXbsWLW1tUk6t8wbcIrPLz673W7eVx64AMu80Z/06C0xqqqqtHv3bg0ZMkSzZ8/WiBEjVF9fr7Vr16qgoEDx8fG+nhMIelxwRqDo9hnD+++/rylTpuipp55SZmamkpOTtXv3bk2cOFH79+/XO++8o6qqqr6YFQDgB90Ow9q1a5WZmammpia9+OKLOnz4sDIzM7Vjxw7t3r1bd999d1/MCfRrLPNGf9Ltp5K++OILvfHGG7rsssu0ePFirVixQvn5+Zo+fXpfzAcEBZZ5oz/pdhiampo0YsSIcw8OCdHQoUM1fvx4X88FBBWWeaM/6fHF5xMnTkg695vOwYMHdfr0abPP5MmTez8dECRY5o3+pEdhuOWWW8z989cVXC6XPB6PXC6Xdz02gEtzu91au3atXn75ZadHAbofhvMfag6ge1jmjf6i22EYPXq0li9frnfffVfffvutbrvtNm3cuFHR0dF9MR8QFN5//33df//9+vbbbyVJzz//vF577TXNnj1biYmJeuedd1jRh4DR7eWqK1eu1LZt23TXXXfpoYceUklJiX7xi1/0xWxA0GCZN/qTbp8x7Ny5U6+//roeeughSdKcOXM0bdo0tbW18RGewEWwzBv9SbfPGI4dO6a0tDTv/alTpyokJETHjx/36WBAMGGZN/qTbp8xtLW1KTQ01H6TkBCdPXvWZ0MBwYhl3ugvevRBPY899pjCwsK8286cOaPMzEwNGzbMu23nzp2+mRAIEizzRn/hk89jmDNnjk+GAYIVy7zRn/AJboAfsMwb/YnPP6inJzZt2qT4+HiFh4crKSlJpaWlF913586dmjlzpkaNGqXIyEilpKRo165dfpwW6D6WeaM/cTwMRUVFWrp0qXJzc1VZWam0tDTNmjVL1dXVne7/8ccfa+bMmSouLlZFRYVmzJihe+65R5WVlX6eHOi688u8CwsLtXHjRn3wwQd69913uaaAgOR4GNavX6/58+drwYIFmjhxovLz8xUbG6vNmzd3un9+fr5WrFihH/zgB7rmmmu0bt06XXPNNXrvvff8PDnQdSzzRn/iaBhaW1tVUVGh9PR0sz09PV3l5eVd+h7t7e1qbm7W5Zdf3hcjAj7BMm/0Jz16d1Vfqa+vV1tbm2JiYsz2mJgY73rvS3nppZd0+vRpzZ49+6L7tLS0qKWlxXu/qampZwMDPeSrZd4cy/AHR8Nw3oUfcXh+TfelvP3221q1apX+8pe/aPTo0RfdLy8vT6tXr+71nEBP+WqZN8cy/MHRMERHR2vw4MEdzg5qa2s7nEVcqKioSPPnz9c777yj22677Tv3zcnJUXZ2tvd+U1OTYmNjez440E2+WubNsQx/cPQaQ2hoqJKSklRSUmK2l5SUKDU19aKPe/vtt/XYY4/prbfe0l133XXJvycsLEyRkZHmBvRHHMvwB8efSsrOztbcuXOVnJyslJQUFRYWqrq6WpmZmZLO/YZUU1Oj7du3SzoXhXnz5mnDhg266aabvGcbQ4cOVVRUlGM/BwAEC8fDkJGRoYaGBq1Zs0Zut1uJiYkqLi5WXFycpHMfefj/X9Pw6quv6uzZs3ryySf15JNPerc/+uij2rZtm7/HB4Cg43gYJCkrK0tZWVmdfu3C/9jv2bOn7wcCgAHM8Re4AQACC2EAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAARkCEYdOmTYqPj1d4eLiSkpJUWlp60X3dbrceeeQRTZgwQYMGDdLSpUv9NygADACOh6GoqEhLly5Vbm6uKisrlZaWplmzZqm6urrT/VtaWjRq1Cjl5ubq+uuv9/O0ABD8HA/D+vXrNX/+fC1YsEATJ05Ufn6+YmNjtXnz5k73HzdunDZs2KB58+YpKirKz9MCQPBzNAytra2qqKhQenq62Z6enq7y8nKHpgKAgS3Eyb+8vr5ebW1tiomJMdtjYmJ04sQJn/09LS0tamlp8d5vamry2fcG/IljGf7g+FNJkuRyucx9j8fTYVtv5OXlKSoqynuLjY312fcG/IljGf7gaBiio6M1ePDgDmcHtbW1Hc4ieiMnJ0eNjY3e27Fjx3z2vQF/4liGPzgahtDQUCUlJamkpMRsLykpUWpqqs/+nrCwMEVGRpob0B9xLMMfHL3GIEnZ2dmaO3eukpOTlZKSosLCQlVXVyszM1PSud+QampqtH37du9j9u/fL0k6deqU6urqtH//foWGhmrSpElO/AgAEFQcD0NGRoYaGhq0Zs0aud1uJSYmqri4WHFxcZLOvaDtwtc0TJkyxfvniooKvfXWW4qLi9PRo0f9OToABCXHwyBJWVlZysrK6vRr27Zt67DN4/H08UQAMHAFxKokAEDgIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDACIgwbNq0SfHx8QoPD1dSUpJKS0u/c/+//e1vSkpKUnh4uBISElRQUOCnSQEg+DkehqKiIi1dulS5ubmqrKxUWlqaZs2aperq6k73P3LkiO68806lpaWpsrJSv/71r7V48WLt2LHDz5MDQHByPAzr16/X/PnztWDBAk2cOFH5+fmKjY3V5s2bO92/oKBAV111lfLz8zVx4kQtWLBAjz/+uF588UU/Tw4AwcnRMLS2tqqiokLp6elme3p6usrLyzt9zN69ezvsf/vtt2vfvn369ttv+2xWABgoQpz8y+vr69XW1qaYmBizPSYmRidOnOj0MSdOnOh0/7Nnz6q+vl5jxozp8JiWlha1tLR47zc2NkqSmpqaujVvW8s33dofgaW7/3///8d4PB5fj9MjHMuQ+v5YdjQM57lcLnPf4/F02Hap/Tvbfl5eXp5Wr17dYXtsbGx3R0U/FvW7zB4/trm5WVFRUT6cpmc4liH1/bHsaBiio6M1ePDgDmcHtbW1Hc4Kzrviiis63T8kJEQjR47s9DE5OTnKzs723m9vb9fJkyc1cuTI7wzQQNLU1KTY2FgdO3ZMkZGRTo8TMDwej5qbmzV27FinR5HEsdwVHMud686x7GgYQkNDlZSUpJKSEt17773e7SUlJfrJT37S6WNSUlL03nvvmW0fffSRkpOTNWTIkE4fExYWprCwMLNtxIgRvRs+SEVGRvKP6QKBcKZwHsdy13Esd9TVY9nxVUnZ2dn6/e9/ry1btuiLL77QsmXLVF1drczMc6dKOTk5mjdvnnf/zMxMffnll8rOztYXX3yhLVu26PXXX9fy5cud+hEAIKg4fo0hIyNDDQ0NWrNmjdxutxITE1VcXKy4uDhJktvtNq9piI+PV3FxsZYtW6ZXXnlFY8eO1caNG3X//fc79SMAQFBxeQJluQUc1dLSory8POXk5HR4qgLoTziWe48wAAAMx68xAAACC2EAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYPwvRRndHpWZ/OkAAAAASUVORK5CYII=",
"text/plain": [
"
"
],
"text/plain": [
"PRE_female 0 1\n",
"PRE_male \n",
"0 104 63\n",
"1 55 38"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# the chisquare test expects as input the observed frequencies in each category\n",
"# to compute this we can use Panda's function crosstab, which counts how many 0s and how many 1s are in the \"PRE_female\" Series and how many are in the \"PRE_female\" Series\n",
"\n",
"pd.crosstab(df[\"PRE_male\"], df[\"PRE_female\"])"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "RVP5OPOiSxad",
"outputId": "6d1df6d4-d638-4534-b8f7-caf75532ffdd"
},
"outputs": [
{
"data": {
"text/plain": [
"Chi2ContingencyResult(statistic=0.13285927016061816, pvalue=0.7154856772979994, dof=1, expected_freq=array([[102.12692308, 64.87307692],\n",
" [ 56.87307692, 36.12692308]]))"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# main chi square test of independence\n",
"# reports test statistic, p value, df\n",
"\n",
"stats.chi2_contingency(pd.crosstab(df[\"PRE_male\"], df[\"PRE_female\"]))"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "waHOyONa90Jq",
"outputId": "a1f26845-eeca-40a1-e116-522cd528a9bf"
},
"outputs": [
{
"data": {
"text/plain": [
"SignificanceResult(statistic=1.1405483405483405, pvalue=0.6907062188741226)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# do this fisher_exact test instead if your data has small frequencies (<5 obs/cell)\n",
"# now we can run the two way chi square on this frequency table\n",
"# output is chi square statistic and p value\n",
"\n",
"stats.fisher_exact(pd.crosstab(df[\"PRE_male\"], df[\"PRE_female\"]))"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "q7V_Mfrn1dCb",
"outputId": "ef4c0360-ad73-473a-b712-b654518cfcc0"
},
"outputs": [
{
"data": {
"text/plain": [
"0.38846153846153847"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df['PRE_female'].mean()"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "VeWyV5gK1goE",
"outputId": "95eccb49-bbd1-4e69-c661-ff2773ffec8d"
},
"outputs": [
{
"data": {
"text/plain": [
"0.3576923076923077"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df['PRE_male'].mean()"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Z_3WfvuC1pQ-",
"outputId": "110346f1-2801-4403-9c6b-2c4b7474b5e0"
},
"outputs": [
{
"data": {
"text/plain": [
"0.48024533379886347"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df['PRE_male'].std()"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "y5mf5ECK13w6",
"outputId": "5d9ead1e-e1bf-4a8e-8355-f0c134d5bf5f"
},
"outputs": [
{
"data": {
"text/plain": [
"0.4883404432118925"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df['PRE_female'].std()"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "CNiOMf7I0aso"
},
"source": [
"Reporting result:\n",
"\n",
"\"The proportion of women chosen at pretest in the male condition (M=35.7%) was not significantly different (χ2=1.14, P=0.69) from the proportion of women chosen at pretest in the female condition (M=38.8%).\""
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "F6O6J-wnIuGf"
},
"source": [
"### Another example:\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "WgUXbU1S-bKd"
},
"source": [
"## Repeated measures chi square -- McNemar test"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "1hC14ZmzQLgZ"
},
"source": [
"McNemar test of marginal homogeneity is a chi square test of independence for non-independent observations (e.g., test-retest / pre-post / before-after designs)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "bZ4og1XA-jHc"
},
"source": [
"The variable \"PRE_female\" encodes whether participants chose a man (PRE_female==0) or a woman (PRE_female==1) at pretest, in the female condition of this experiment.\n",
"\n",
"The variable \"POST_female\" encodes whether the same participants chose a man (POST_female==0) or a woman (POST_female==1) at posttest, in the female condition of this experiment.\n",
"\n",
"To test whether the proportion of men/women choices are different in the pretest versus posttest measures for the same participants (comparing the PRE_female proportion to the POST_female proportion) we can run a repeated measures chi square test (or a chi square test of independence for repeated measures, or a McNemar tests)."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 441
},
"id": "SA3-89l2-eO9",
"outputId": "51676b0c-c313-4171-c9b8-ca3477117461"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAJOCAYAAACz06ChAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAofUlEQVR4nO3df1SVBZ7H8c8VBPwROIoSGiJomkrNuLBO6KijFi41M021k5OlOcqMRP4kc2Rt0xyNNo3BOmFyyhymcqjRZteJM8VulpjOmY1w7aymU/7AEENwRtRRULj7h+ttvl1QQLgPXN+vc+453ec+z73fu/s4b57n/nK53W63AAD4f52cHgAA0L4QBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiBTg/ghPr6eh09elTXXXedXC6X0+OgHXO73Tp16pT69u2rTp3a399R7Mtoqubsy9dkGI4ePaqoqCinx0AHcuTIEd1www1Oj+GFfRnN1ZR9+ZoMw3XXXSfp4v+BQkNDHZ4G7Vl1dbWioqI8+0x7w76MpmrOvnxNhuHSIXdoaCj/mNAk7fU0Dfsymqsp+3L7O2kKAHAUYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiBTg8AAK1l3rx5On78uCSpd+/eWrNmjcMTdUyEAYDfOH78uL788kunx+jwOJUEADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMBoF2HIyclRTEyMQkJCFB8fr6KiosuuX1NToyVLlig6OlrBwcEaOHCg1q9f76NpAcC/BTo9QH5+vubPn6+cnByNHj1a69atU3Jysvbs2aP+/fs3uM19992nL7/8Ui+//LIGDRqkiooKXbhwwceTA4B/cjwMWVlZmjlzplJSUiRJ2dnZeuedd7R27VplZmZ6rf+HP/xBH3zwgQ4cOKCePXtKkgYMGODLkQHArzl6Kqm2tlbFxcVKSkoyy5OSkrRjx44Gt/mP//gPJSQk6JlnnlG/fv00ePBgLVy4UGfPnvXFyADg9xw9YqisrFRdXZ0iIiLM8oiICB07dqzBbQ4cOKDt27crJCREb731liorK5WWlqYTJ040+jpDTU2NampqPNerq6tb70kAPsS+DF9oFy8+u1wuc93tdnstu6S+vl4ul0uvvfaaRo4cqTvuuENZWVnasGFDo0cNmZmZCgsL81yioqJa/TkAvsC+DF9wNAzh4eEKCAjwOjqoqKjwOoq4JDIyUv369VNYWJhn2dChQ+V2u/XFF180uE1GRoZOnjzpuRw5cqT1ngTgQ+zL8AVHwxAUFKT4+HgVFhaa5YWFhRo1alSD24wePVpHjx7V6dOnPcv279+vTp066YYbbmhwm+DgYIWGhpoL0BGxL8MXHD+VlJ6erpdeeknr16/X3r17tWDBApWWlio1NVXSxb+Qpk2b5ll/ypQp6tWrl37yk59oz5492rZtmx577DHNmDFDXbp0ceppAIDfcPztqpMnT1ZVVZWWL1+u8vJyxcXFqaCgQNHR0ZKk8vJylZaWetbv3r27CgsLNWfOHCUkJKhXr1667777tGLFCqeeAgD4FcfDIElpaWlKS0tr8LYNGzZ4Lbvpppu8Tj8BAFqH46eSAADtC2EAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAY7eIrMQD4VvxjeU6P0CZC/3La89du+V9O++3zLF417corXQWOGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIAR6PQAANBa6jt3a/C/0TyEAYDfOD0k2ekR/AKnkgAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIDRLsKQk5OjmJgYhYSEKD4+XkVFRY2u+/7778vlcnldPv30Ux9ODAD+y/Ew5Ofna/78+VqyZIlKSko0ZswYJScnq7S09LLb7du3T+Xl5Z7LjTfe6KOJAcC/OR6GrKwszZw5UykpKRo6dKiys7MVFRWltWvXXna7Pn366Prrr/dcAgICfDQxAPg3R8NQW1ur4uJiJSUlmeVJSUnasWPHZbcdMWKEIiMjNXHiRG3duvWy69bU1Ki6utpcgI6IfRm+4GgYKisrVVdXp4iICLM8IiJCx44da3CbyMhI5ebmatOmTdq8ebOGDBmiiRMnatu2bY0+TmZmpsLCwjyXqKioVn0egK+wL8MX2sUvuLlcLnPd7XZ7LbtkyJAhGjJkiOd6YmKijhw5otWrV2vs2LENbpORkaH09HTP9erqav5BoUNiX4YvOBqG8PBwBQQEeB0dVFRUeB1FXM6tt96qV199tdHbg4ODFRwc3OI5gfaCfRm+4OippKCgIMXHx6uwsNAsLyws1KhRo5p8PyUlJYqMjGzt8QDgmuT4qaT09HRNnTpVCQkJSkxMVG5urkpLS5Wamirp4qFzWVmZ8vLyJEnZ2dkaMGCAhg8frtraWr366qvatGmTNm3a5OTTAAC/4XgYJk+erKqqKi1fvlzl5eWKi4tTQUGBoqOjJUnl5eXmMw21tbVauHChysrK1KVLFw0fPlxvv/227rjjDqeeAgD4FZfb7XY7PYSvVVdXKywsTCdPnlRoaKjT46Ada+/7Skvni38srw2nQlsrXjWt2ds0Z19x/ANuAID2hTAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAACPQ6QHgvHnz5un48eOSpN69e2vNmjUOTwTASYQBOn78uL788kunxwDQTnAqCQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgtIsw5OTkKCYmRiEhIYqPj1dRUVGTtvvwww8VGBiob33rW207IABcQxwPQ35+vubPn68lS5aopKREY8aMUXJyskpLSy+73cmTJzVt2jRNnDjRR5MCwLXB8TBkZWVp5syZSklJ0dChQ5Wdna2oqCitXbv2stvNmjVLU6ZMUWJioo8mBYBrg6NhqK2tVXFxsZKSkszypKQk7dixo9HtXnnlFX3++edaunRpkx6npqZG1dXV5gJ0ROzL8AVHw1BZWam6ujpFRESY5RERETp27FiD2/z5z3/W4sWL9dprrykwMLBJj5OZmamwsDDPJSoq6qpnB5zAvgxfcPxUkiS5XC5z3e12ey2TpLq6Ok2ZMkVPPvmkBg8e3OT7z8jI0MmTJz2XI0eOXPXMgBPYl+ELTfuTu42Eh4crICDA6+igoqLC6yhCkk6dOqWPPvpIJSUlmj17tiSpvr5ebrdbgYGBevfddzVhwgSv7YKDgxUcHNw2TwLwIfZl+IKjRwxBQUGKj49XYWGhWV5YWKhRo0Z5rR8aGqpPPvlEu3bt8lxSU1M1ZMgQ7dq1S9/+9rd9NToA+C1HjxgkKT09XVOnTlVCQoISExOVm5ur0tJSpaamSrp46FxWVqa8vDx16tRJcXFxZvs+ffooJCTEazkAoGUcD8PkyZNVVVWl5cuXq7y8XHFxcSooKFB0dLQkqby8/IqfaQAAtB7HwyBJaWlpSktLa/C2DRs2XHbbZcuWadmyZa0/FABco9rFu5IAAO0HYQAAGIQBAGC0i9cYOor4x/KcHqFNhP7ltOcvhPK/nPbb51m8aprTIxjnzp1TSEiI02MAXjhiAHyovr5ev/jFL9SvXz91795dBw4ckCT967/+q15++WWHpwMuIgyAD61YsUIbNmzQM888o6CgIM/ym2++WS+99JKDkwFfIQyAD+Xl5Sk3N1cPPPCAAgICPMtvueUWffrppw5OBnyFMAA+VFZWpkGDBnktr6+v1/nz5x2YCPB21WH47LPP9M477+js2bOSLn4zKoCGDR8+vMGfrn3zzTc1YsQIByYCvLX4XUlVVVWaPHmy3nvvPblcLv35z39WbGysUlJS1KNHDz377LOtOSfgF5YuXaqpU6eqrKxM9fX12rx5s/bt26e8vDz9/ve/d3o8QNJVHDEsWLBAgYGBKi0tVdeuXT3LJ0+erD/84Q+tMhzgb77//e8rPz9fBQUFcrlceuKJJ7R3715t2bJFt99+u9PjAZKu4ojh3Xff1TvvvKMbbrjBLL/xxht1+PDhqx4M8FeTJk3SpEmTnB4DaFSLjxjOnDljjhQuqays5IdEAKADa/ERw9ixY5WXl6df/OIXki7+PGd9fb1WrVql8ePHt9qAQEf3jW98o8Gfqm3IiRMn2nga4MpaHIZVq1bpu9/9rj766CPV1tZq0aJF+t///V+dOHFCH374YWvOCHRo2dnZTo8ANEuLwzBs2DDt3r1ba9euVUBAgM6cOaN77rlHjzzyiCIjI1tzRqBDe+ihh5weAWiWq/oSveuvv15PPvlka80CXFPOnj3r9aG20NBQh6YBvtKsMOzevbvJ695yyy3NHgbwd2fOnNHPf/5zvfHGG6qqqvK6va6uzoGpAKtZYfjWt74ll8t1xU83u1wudnCgAYsWLdLWrVuVk5OjadOm6YUXXlBZWZnWrVunp59+2unxAEnNDMPBgwfbag7gmrBlyxbl5eXpu9/9rmbMmKExY8Zo0KBBio6O1muvvaYHHnjA6RGB5oUhOjq6reYArgknTpxQTEyMpIuvJ1x6e+p3vvMdPfzww06OBnhc9S+47dmzR6WlpaqtrTXLf/CDH1ztXQN+JzY2VocOHVJ0dLSGDRumN954QyNHjtSWLVvUo0cPp8cDJF1FGA4cOKC7775bn3zyiXnd4dIHeXiNAfD2k5/8RP/zP/+jcePGKSMjQ3feeaeef/55XbhwQVlZWU6PB0i6ijDMmzdPMTEx+s///E/FxsbqT3/6k6qqqvToo49q9erVrTkj4DcWLFjg+e/x48fr008/1UcffaSBAwfqm9/8poOTAV9pcRh27typ9957T71791anTp3UqVMnfec731FmZqbmzp2rkpKS1pwT8Ev9+/dX//79nR4DMFochrq6OnXv3l2SFB4erqNHj2rIkCGKjo7Wvn37Wm1AwN/86U9/0vvvv6+KigrV19eb2zidhPagxWGIi4vT7t27FRsbq29/+9ueHzfPzc1VbGxsa84I+I2nnnpKjz/+uIYMGaKIiAjz5XpN/aI9oK21OAyPP/64zpw5I0lasWKFvve972nMmDHq1auX8vPzW21AwJ+sWbNG69ev1/Tp050eBWhUi8Pw9z80Ehsbqz179ujEiRPN+oph4FrTqVMnjR492ukxgMtq8Q/1NKRnz55EAbiMBQsW6IUXXnB6DOCyWnzEcO7cOT3//PPaunVrgy+iffzxx1c9HOBvFi5cqDvvvFMDBw7UsGHD1LlzZ3P75s2bHZoM+EqLwzBjxgwVFhbqn//5nzVy5EiOFIAmmDNnjrZu3arx48erV69e/LtBu9TiMLz99tsqKCjgfCnQDHl5edq0aZPuvPNOp0cBGtXi1xj69eun6667rjVnAfxez549NXDgQKfHAC6rxWF49tln9fOf/1yHDx9uzXkAv7Zs2TItXbpUf/vb35weBWhUi08lJSQk6Ny5c4qNjVXXrl29XkS79HXCAL7y3HPP6fPPP1dERIQGDBjg9e+GN22gPWhxGO6//36VlZXpqaee8voEJ4CG/fCHP3R6BOCKWhyGHTt2aOfOnXwjJNAMS5cudXoE4Ipa/BrDTTfdpLNnz7bmLMA14a9//ateeuklZWRkeE65fvzxxyorK3N4MuCiFh8xPP3003r00Ue1cuVK3XzzzV7nSkNDQ696OMDf7N69W7fddpvCwsJ06NAh/fSnP1XPnj311ltv6fDhw8rLy3N6RKDlYfinf/onSdLEiRPNcrfbLZfLxS+4AQ1IT0/X9OnT9cwzz5i3eycnJ2vKlCkOTgZ8pcVh2Lp1a2vOAVwT/vu//1vr1q3zWt6vXz8dO3bMgYkAby0Ow7hx41pzDuCaEBISourqaq/l+/btU+/evR2YCPB2Vd+uWlRUpAcffFCjRo3yvHD261//Wtu3b2+V4QB/c9ddd2n58uU6f/68pIs/zlNaWqrFixfr3nvvdXg64KIWh2HTpk2aNGmSunTpoo8//lg1NTWSpFOnTumpp55qtQEBf7J69WodP35cffr00dmzZzVu3DgNGjRI1113nVauXOn0eICkqziVtGLFCr344ouaNm2afvOb33iWjxo1SsuXL2+V4QB/Exoaqu3bt+u9997Txx9/rPr6ev3DP/yDbrvtNqdHAzxaHIZ9+/Zp7NixXstDQ0P117/+9WpmAvxKz549tX//foWHh2vGjBlas2aNJkyYoAkTJjg9GtCgFp9KioyM1Geffea1fPv27YqNjb2qoQB/Ultb63nB+Ve/+pXOnTvn8ETA5bX4iGHWrFmaN2+e1q9fL5fLpaNHj2rnzp1auHChnnjiidacEejQEhMT9cMf/lDx8fFyu92aO3euunTp0uC669ev9/F0gLcWh2HRokU6efKkxo8fr3Pnzmns2LEKDg7WwoULNXv27NacEejQXn31Vf3yl7/U559/LpfLpZMnT3LUgHatWWHYvXu34uLi1KnTxTNQK1eu1JIlS7Rnzx7V19dr2LBh6t69e5sMCnRUERERevrppyVJMTEx+vWvf61evXo5PBXQuGa9xjBixAhVVlZKkmJjY1VVVaWuXbsqISFBI0eOJArAFRw8eLBJUbj55pt15MgRH0wEeGtWGHr06KGDBw9Kkg4dOqT6+vo2GQq+Vd+5m+qD/v/SuZvT40AX/31d+hAc4GvNOpV07733aty4cYqMjJTL5VJCQoICAgIaXPfAgQOtMiDa3ukhyU6PAKAdaVYYcnNzdc899+izzz7T3Llz9dOf/tR8QyQAoONr9ruSLn3ddnFxsebNm3fFMHzxxRfq27ev5wVrAED71uL/tX7llVeadLQwbNgwHTp0qKUPAwDwsTb/M97tdrf1QwAAWhHndwAfmDBhQrO+Q2zdunWKiIhou4GAy2jxJ58BNN3777+v2traJq/Pz3zCSRwxAACMNj9icLlcbf0QQIdw6tQphYSEXHad0NBQH00DNK7Nw8CLz8BFgwcPbvQ2t9stl8uluro6H04ENKzVw+B2uz0/XShJe/bsUd++fVv7YYAO57e//a169uzp9BjAFTU7DF27dtXhw4fVu3dvSRc/8PbKK68oMjJSklRRUaG+fft6/vKJiopqxXGBjmv06NGeP5iA9qzZLz6fO3fOnB768MMPdfbsWbMOp48AoONqk3cl8YIzYEVHRzf6hZNAe8PnGAAfuPR19V/3wQcf6MyZM0pMTNQ3vvENH08FNKzZYXC5XOaI4OvXAXhbtWqVTp8+rSeffFLSxdOtycnJevfddyVJffr00X/9139p+PDhTo4JSGrBqSS3263BgwerZ8+e6tmzp06fPq0RI0Z4rt90001tMSfQoW3cuFHDhg3zXP/tb3+rbdu2qaioSJWVlUpISPBEA3Bas48YXnnllbaYA/BrBw8e1C233OK5XlBQoHvvvVejR4+WJD3++OP60Y9+5NR4gNHsMDz00ENtMQfg186fP6/g4GDP9Z07d2revHme63379vX8njrgtFZ/V1J5eblmz57drG1ycnIUExOjkJAQxcfHq6ioqNF1t2/frtGjR6tXr17q0qWLbrrpJv3yl7+82rGBNjVo0CBt27ZNklRaWqr9+/dr3Lhxntu/+OIL9erVy6nxAKNF70ras2ePtm7dqs6dO+u+++5Tjx49VFlZqZUrV+rFF19UTExMk+8rPz9f8+fPV05OjkaPHq1169YpOTlZe/bsUf/+/b3W79atm2bPnq1bbrlF3bp10/bt2zVr1ix169ZNP/vZz1rydIA29/DDD2v27NkqKirSH//4RyUmJprXHN577z2NGDHCwQmBrzT7iOH3v/+9RowYoTlz5ig1NVUJCQnaunWrhg4dql27dunNN9/Unj17mnx/WVlZmjlzplJSUjR06FBlZ2crKipKa9eubXD9ESNG6P7779fw4cM1YMAAPfjgg5o0adJljzIAp82aNUtr1qzRiRMnNHbsWG3atMncfvToUc2YMcOh6QCr2WFYuXKlUlNTVV1drdWrV+vAgQNKTU3Vpk2btHXrVn3ve99r8n3V1taquLhYSUlJZnlSUpJ27NjRpPsoKSnRjh07zGE50B7NnDlTb731ltauXavrr7/e3JaTk6O7777bockAq9mnkvbu3atf/epX6t69u+bOnatFixYpOztbY8eObfaDV1ZWqq6uzuuXqiIiInTs2LHLbnvDDTfo+PHjunDhgpYtW6aUlJRG162pqVFNTY3nenV1dbNnBVpDWVmZNm3apP3798vlcmnw4MG655571K9fvyZtz74MX2h2GKqrq9WjR4+LGwcGqkuXLpf9OuGm+PoH5C59BfHlFBUV6fTp0/rjH/+oxYsXa9CgQbr//vsbXDczM5P3iMNxOTk5Sk9PV21trcLCwuR2u1VdXa3HHntMWVlZSktLu+J9sC/DF1r84vOlv+jdbrf27dunM2fOmHX+/j3bjQkPD1dAQIDX0UFFRcUVf+/20gvcN998s7788kstW7as0TBkZGQoPT3dc726uppvfYVPvf3225o7d67mz5+vRx991PNtxOXl5Vq1apXmzZunAQMG6I477rjs/bAvwxdaFIYJEyaY65deV3C5XM36wZGgoCDFx8ersLDQnF8tLCzUXXfd1eR53G63Obz+uuDgYPMecsDXnnnmGS1evFgrVqwwyyMjI5WVlaWuXbvq3/7t364YBvZl+EKzw9DYl4G1VHp6uqZOnaqEhAQlJiYqNzdXpaWlSk1NlXTxL6SysjLl5eVJkl544QX179/f89Ub27dv1+rVqzVnzpxWnQtoTSUlJcrNzW309qlTp2rNmjU+nAhoXLPD0KdPHy1cuFC/+93vdP78ed1222167rnnFB4e3qIBJk+erKqqKi1fvlzl5eWKi4tTQUGBoqOjJV081C4tLfWsX19fr4yMDB08eFCBgYEaOHCgnn76ac2aNatFjw/4Qn19vTp37tzo7Z07d+Z3TNBuNDsMTzzxhDZs2KAHHnhAISEh2rhxox5++GG9+eabLR4iLS2t0RfeNmzYYK7PmTOHowN0OMOHD9e///u/a8GCBQ3e/rvf/Y5vVkW70ewwbN68WS+//LJ+/OMfS5IefPBBjR49WnV1dfwQCdCItLQ0PfzwwwoODtbPfvYzBQZe/Kd34cIFrVu3To8//rhycnIcnhK4qNlhOHLkiMaMGeO5PnLkSAUGBuro0aO8OwJoxEMPPaRPPvlEs2fPVkZGhgYOHChJ+vzzz3X69GnNnTtX06dPd3ZI4P81Owx1dXUKCgqydxIYqAsXLrTaUIA/Wr16tX70ox9p48aN2r9/vyRp7Nix+vGPf6xbb73V4emArzQ7DG63W9OnTzdvmTt37pxSU1PVrVs3z7LNmze3zoSAH/jb3/6mxx57zPOmjYkTJ+r5559v8Zs2gLbUKr/H8OCDD7bKMIC/Wrp0qedNG126dNHrr79+1W/aANoKv+AG+MDX37TxwAMP8KYNtFut/kM9ALxd7k0bQHtDGAAf4E0b6Eha9F1JAJqHN22gIyEMgA/wpg10JIQB8AHetIGOhNcYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAABGuwhDTk6OYmJiFBISovj4eBUVFTW67ubNm3X77berd+/eCg0NVWJiot555x0fTgsA/s3xMOTn52v+/PlasmSJSkpKNGbMGCUnJ6u0tLTB9bdt26bbb79dBQUFKi4u1vjx4/X9739fJSUlPp4cAPyT42HIysrSzJkzlZKSoqFDhyo7O1tRUVFau3Ztg+tnZ2dr0aJF+sd//EfdeOONeuqpp3TjjTdqy5YtPp4cAPyTo2Gora1VcXGxkpKSzPKkpCTt2LGjSfdRX1+vU6dOqWfPnm0xIgBccwKdfPDKykrV1dUpIiLCLI+IiNCxY8eadB/PPvuszpw5o/vuu6/RdWpqalRTU+O5Xl1d3bKBAYexL8MXHD+VJEkul8tcd7vdXssasnHjRi1btkz5+fnq06dPo+tlZmYqLCzMc4mKirrqmQEnsC/DFxwNQ3h4uAICAryODioqKryOIr4uPz9fM2fO1BtvvKHbbrvtsutmZGTo5MmTnsuRI0euenbACezL8AVHwxAUFKT4+HgVFhaa5YWFhRo1alSj223cuFHTp0/X66+/rjvvvPOKjxMcHKzQ0FBzAToi9mX4gqOvMUhSenq6pk6dqoSEBCUmJio3N1elpaVKTU2VdPEvpLKyMuXl5Um6GIVp06ZpzZo1uvXWWz1HG126dFFYWJhjzwMA/IXjYZg8ebKqqqq0fPlylZeXKy4uTgUFBYqOjpYklZeXm880rFu3ThcuXNAjjzyiRx55xLP8oYce0oYNG3w9PgD4HcfDIElpaWlKS0tr8Lav/4/9+++/3/YDAcA1rF28KwkA0H4QBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGAQBgCAQRgAAAZhAAAYhAEAYBAGAIBBGAAABmEAABiEAQBgEAYAgEEYAAAGYQAAGIQBAGC0izDk5OQoJiZGISEhio+PV1FRUaPrlpeXa8qUKRoyZIg6deqk+fPn+25QALgGOB6G/Px8zZ8/X0uWLFFJSYnGjBmj5ORklZaWNrh+TU2NevfurSVLluib3/ymj6cFAP/neBiysrI0c+ZMpaSkaOjQocrOzlZUVJTWrl3b4PoDBgzQmjVrNG3aNIWFhfl4WgDwf4FOPnhtba2Ki4u1ePFiszwpKUk7duxotcepqalRTU2N53p1dXWr3TfgS+zL8AVHjxgqKytVV1eniIgIszwiIkLHjh1rtcfJzMxUWFiY5xIVFdVq9w34EvsyfMHxU0mS5HK5zHW32+217GpkZGTo5MmTnsuRI0da7b4BX2Jfhi84eiopPDxcAQEBXkcHFRUVXkcRVyM4OFjBwcGtdn+AU9iX4QuOHjEEBQUpPj5ehYWFZnlhYaFGjRrl0FQAcG1z9IhBktLT0zV16lQlJCQoMTFRubm5Ki0tVWpqqqSLh85lZWXKy8vzbLNr1y5J0unTp3X8+HHt2rVLQUFBGjZsmBNPAQD8iuNhmDx5sqqqqrR8+XKVl5crLi5OBQUFio6OlnTxA21f/0zDiBEjPP9dXFys119/XdHR0Tp06JAvRwcAv+R4GCQpLS1NaWlpDd62YcMGr2Vut7uNJwKAa1e7eFcSAKD9IAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAAzCAAAwCAMAwCAMAACDMAAADMIAADDaRRhycnIUExOjkJAQxcfHq6io6LLrf/DBB4qPj1dISIhiY2P14osv+mhSAPB/jochPz9f8+fP15IlS1RSUqIxY8YoOTlZpaWlDa5/8OBB3XHHHRozZoxKSkr0L//yL5o7d642bdrk48kBwD85HoasrCzNnDlTKSkpGjp0qLKzsxUVFaW1a9c2uP6LL76o/v37Kzs7W0OHDlVKSopmzJih1atX+3hyAPBPjoahtrZWxcXFSkpKMsuTkpK0Y8eOBrfZuXOn1/qTJk3SRx99pPPnz7fZrABwrQh08sErKytVV1eniIgIszwiIkLHjh1rcJtjx441uP6FCxdUWVmpyMhIr21qampUU1PjuX7y5ElJUnV1dbPmras526z10b409//ff7+N2+1u7XFahH0ZUtvvy46G4RKXy2Wuu91ur2VXWr+h5ZdkZmbqySef9FoeFRXV3FHRgYU9n9ribU+dOqWwsLBWnKZl2Jchtf2+7GgYwsPDFRAQ4HV0UFFR4XVUcMn111/f4PqBgYHq1atXg9tkZGQoPT3dc72+vl4nTpxQr169Lhuga0l1dbWioqJ05MgRhYaGOj1Ou+F2u3Xq1Cn17dvX6VEksS83Bftyw5qzLzsahqCgIMXHx6uwsFB33323Z3lhYaHuuuuuBrdJTEzUli1bzLJ3331XCQkJ6ty5c4PbBAcHKzg42Czr0aPH1Q3vp0JDQ/nH9DXt4UjhEvblpmNf9tbUfdnxdyWlp6frpZde0vr167V3714tWLBApaWlSk29eKiUkZGhadOmedZPTU3V4cOHlZ6err1792r9+vV6+eWXtXDhQqeeAgD4FcdfY5g8ebKqqqq0fPlylZeXKy4uTgUFBYqOjpYklZeXm880xMTEqKCgQAsWLNALL7ygvn376rnnntO9997r1FMAAL/icreXt1vAUTU1NcrMzFRGRobXqQqgI2FfvnqEAQBgOP4aAwCgfSEMAACDMAAADMIAADAIAwDAIAwAAIMwAAAMwgAAMAgDAMAgDAAAgzAAAIz/AxM+T5DYrmM7AAAAAElFTkSuQmCC",
"text/plain": [
"
"
],
"text/plain": [
"POST_female 0 1\n",
"PRE_female \n",
"0 85 74\n",
"1 24 77"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# the chisquare test expects as input the observed frequencies in each category\n",
"# to compute this we can use Panda's function crosstab, which counts how many 0s and how many 1s are in the \"PRE_female\" Series and how many are in the \"POST_female\" Series\n",
"\n",
"pd.crosstab(df[\"PRE_female\"], df[\"POST_female\"])"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "2AUUJITO_uR_",
"outputId": "9d32e466-91d4-45e5-e829-40d171adf0a7"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"pvalue 4.2159207369916054e-07\n",
"statistic 24.0\n"
]
}
],
"source": [
"# now we can run the McNemar test on this frequency table\n",
"# output is chi square statistic and p value\n",
"\n",
"print(mcnemar(pd.crosstab(df[\"PRE_female\"], df[\"POST_female\"])))"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "RcsYvXka2gH8",
"outputId": "0f2d6e60-30be-4a58-c643-ad21b24068b4"
},
"outputs": [
{
"data": {
"text/plain": [
"0.38846153846153847"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df[\"PRE_female\"].mean()"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Unb_zQ3f2kdX",
"outputId": "6a1be4c0-f08e-403b-a8de-464005792f01"
},
"outputs": [
{
"data": {
"text/plain": [
"0.4883404432118925"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df[\"PRE_female\"].std()"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "bf1AXL8B2kfk",
"outputId": "8da17925-9da2-43b6-f057-967e99a98f1b"
},
"outputs": [
{
"data": {
"text/plain": [
"0.5807692307692308"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df[\"POST_female\"].mean()"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "siMfz-nY2khy",
"outputId": "fb3de747-a106-4e53-d0de-a6969d8a7635"
},
"outputs": [
{
"data": {
"text/plain": [
"0.4943848646716376"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df[\"POST_female\"].std()"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "27NK885t2K7z"
},
"source": [
"Reporting result:\n",
"\n",
"\"The proportion of women chosen at posttest in the female condition (M=58%) was significantly higher (χ2=24, P<0.001) than the proportion of women chosen at pretest in the female condition (M=38.8%).\""
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "ODHos2hurTWC"
},
"source": [
"## Power"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "YxaQheK5rWF5"
},
"source": [
"the probability of detecting a significant effect, given that the effect is real\n",
"\n",
"The power is affected by at least three factors:\n",
"- Signicance level (α, typically 0.05): the higher the significance level, the higher the power\n",
"- Sample size (n): the greater the sample size, the greater the power\n",
"- Effect size (ES): the greater the effect size, the greater power\n",
"- Other: the tests methods, distribution of predictors, missing data\n",
"\n",
"Compute power: WebPower: https://webpower.psychstat.org/wiki/models/index"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "0HhFCl_UrZgQ"
},
"source": [
"Report power analysis:\n",
"\n",
"\"For a power analysis we used the software webpower (Zhang & Yuan, 2018), and we calculated that in order to detect an effect of at least 0.2, at a significance level of 0.05, in a two sided comparison, with a power of 0.95, we need a sample size of 325 observations (participants).”\n"
]
}
],
"metadata": {
"colab": {
"provenance": []
},
"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.11.7"
}
},
"nbformat": 4,
"nbformat_minor": 4
}