From e28b55afa40d21ee9d4be1da4e3a6cd96fcf8d74 Mon Sep 17 00:00:00 2001 From: "Alexander E. Zarebski" <36244258+aezarebski@users.noreply.github.com> Date: Tue, 22 Nov 2022 22:25:30 +0000 Subject: [PATCH] Ex4 (#68) * rename old materials * update requirements file * include the renamed notebook and md pairing * update answers notebook * include question notesbooks * update checker script * update docs --- README.org | 1 + additional-resources.org | 4 + example-4/example-4-answers.ipynb | 526 ++++++++++++++++++ example-4/example-4-answers.md | 241 ++++++++ example-4/example-4-questions.ipynb | 282 ++++++++++ example-4/example-4-questions.md | 127 +++++ {example-4 => old-example-4}/README.org | 0 .../cat-grades-by-breed.csv | 0 .../cat-weights-by-breed.csv | 0 {example-4 => old-example-4}/cat-weights.csv | 0 .../example-4-answers-grades.ipynb | 0 .../example-4-answers-weight.ipynb | 0 .../example-4-data-grades.R | 0 .../example-4-data-weight.R | 0 .../example-4-questions-grades.ipynb | 0 .../example-4-questions-weight.ipynb | 0 requirements.txt | 1 - resources/check-questions-text.py | 2 +- 18 files changed, 1182 insertions(+), 2 deletions(-) create mode 100644 example-4/example-4-answers.ipynb create mode 100644 example-4/example-4-answers.md create mode 100644 example-4/example-4-questions.ipynb create mode 100644 example-4/example-4-questions.md rename {example-4 => old-example-4}/README.org (100%) rename {example-4 => old-example-4}/cat-grades-by-breed.csv (100%) rename {example-4 => old-example-4}/cat-weights-by-breed.csv (100%) rename {example-4 => old-example-4}/cat-weights.csv (100%) rename {example-4 => old-example-4}/example-4-answers-grades.ipynb (100%) rename {example-4 => old-example-4}/example-4-answers-weight.ipynb (100%) rename {example-4 => old-example-4}/example-4-data-grades.R (100%) rename {example-4 => old-example-4}/example-4-data-weight.R (100%) rename {example-4 => old-example-4}/example-4-questions-grades.ipynb (100%) rename {example-4 => old-example-4}/example-4-questions-weight.ipynb (100%) diff --git a/README.org b/README.org index dc6f91f..24485b4 100644 --- a/README.org +++ b/README.org @@ -20,6 +20,7 @@ to pull down changes from the main branch each week for the tutorial. and hypothesis tests. - [[https://github.com/aezarebski/aas-extended-examples/tree/main/example-2][Example 2]] reviews simple linear regression. - [[https://github.com/aezarebski/aas-extended-examples/tree/main/example-3][Example 3]] reviews multiple linear regression. +- [[https://github.com/aezarebski/aas-extended-examples/tree/main/example-4][Example 4]] reviews hypothesis tests for contingency tables. - The [[https://github.com/aezarebski/aas-extended-examples/blob/main/setup.org][setup]] page gives some instructions for setting up the software and notebooks. - The [[https://github.com/aezarebski/aas-extended-examples/blob/main/additional-resources.org][additional resources]] page gives some pointers to additional resources diff --git a/additional-resources.org b/additional-resources.org index f153638..ab0112a 100644 --- a/additional-resources.org +++ b/additional-resources.org @@ -29,6 +29,10 @@ Some of these may be helpful for inspiration - [[https://en.wikipedia.org/wiki/Politics_and_the_English_Language][Politics and the English Language]] by George Orwell is great, in particular the six rules given within. - [[http://paulgraham.com/simply.html][Write Simply]] by Paul Graham +** LaTeX + +- If you want to use LaTeX, the [[https://www.overleaf.com/learn][Overleaf documentation]] can be a very good resource. + ** Tools Some of these may be helpful for mechanical corrections diff --git a/example-4/example-4-answers.ipynb b/example-4/example-4-answers.ipynb new file mode 100644 index 0000000..aff971c --- /dev/null +++ b/example-4/example-4-answers.ipynb @@ -0,0 +1,526 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a61f4b1e-5b27-4cdc-83af-1119d0f692a9", + "metadata": {}, + "source": [ + "# Example 4\n", + "\n", + "This notebook is available on github\n", + "[here](https://github.com/aezarebski/aas-extended-examples). If you find\n", + "errors or would like to suggest an improvement, feel free to create an\n", + "issue.\n", + "\n", + "As usual we will start by importing some useful libraries." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "abb9087e-7cbf-43fd-9fd3-6465f613cf03", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import statsmodels.api as sm\n", + "from statsmodels.graphics.mosaicplot import mosaic\n", + "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "import scipy\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "id": "65f9801b-b722-4312-b720-2edb27a6d175", + "metadata": {}, + "source": [ + "Today we will look at a dataset from a double-blind clinical trial of a new\n", + "treatment for rheumatoid arthritis. We will test whether treatment is correlated\n", + "with a change in symptoms using a $\\chi^{2}$-test.\n", + "\n", + "First, we need to load the data which comes bundled with `statsmodels`." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "de17ffee-8a0c-4730-949a-a7c621c5a1b5", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
IDTreatmentSexAgeImproved
057TreatedMale27Some
146TreatedMale29None
277TreatedMale30None
317TreatedMale32Marked
436TreatedMale46Marked
\n", + "
" + ], + "text/plain": [ + " ID Treatment Sex Age Improved\n", + "0 57 Treated Male 27 Some\n", + "1 46 Treated Male 29 None\n", + "2 77 Treated Male 30 None\n", + "3 17 Treated Male 32 Marked\n", + "4 36 Treated Male 46 Marked" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ra = sm.datasets.get_rdataset(\"Arthritis\", \"vcd\").data\n", + "ra.head()" + ] + }, + { + "cell_type": "markdown", + "id": "c3729f85-5d76-4768-9c94-24533dbd34f1", + "metadata": {}, + "source": [ + "### Question\n", + "\n", + "Use `pandas` to generate a cross tabulation of the treatment status and\n", + "improvement.\n", + "\n", + "[hint](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.crosstab.html)\n", + "\n", + "### Answer" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "aa6959e1-34d2-4f4c-9783-4c756dda3407", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Improved Marked None Some\n", + "Treatment \n", + "Placebo 7 29 7\n", + "Treated 21 13 7\n" + ] + } + ], + "source": [ + "outcome_tbl = pd.crosstab(ra.Treatment, ra.Improved)\n", + "print(outcome_tbl)" + ] + }, + { + "cell_type": "markdown", + "id": "e07a21a7-d039-45c8-b89f-f0051b889cab", + "metadata": {}, + "source": [ + "### Question\n", + "\n", + "Generate a mosaic plot to display this data.\n", + "\n", + "[hint](https://www.statsmodels.org/dev/generated/statsmodels.graphics.mosaicplot.mosaic.html)\n", + "\n", + "### Answer" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "07c7b300-f8a2-44e9-8bc9-0dd04f720d68", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAGzCAYAAAAhXWNYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMN0lEQVR4nO3deXgV1cHH8e/cPclNbnYSIAEEggmgbIKigigqVpQXsFhX0LpU1KIVq7alYn0rUkVstWrf962IFsWligsqWguIuKBUQCASQBJAluz73e+8fwQvXAPKngR+n+fJ89w7Z+bMmZGY35xzZsYwTdNERERE5DhnaekGiIiIiLQGCkUiIiIiKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQi++mss86iV69eh1THokWLMAyDRYsW/ei6xcXFGIbBM888c0j7FBHZXwpFIseYJ554AsMwGDRo0AFvu23bNqZOncqKFSsOf8P24fnnn+fRRx89InV/F6z256e4uPiItOH7HnjgAebNm3dU9nWoPv74Y6ZOnUp1dXVLN0XkqDD07jORY8vpp5/Otm3bKC4uZv369XTr1m2/t/3iiy845ZRTmDVrFhMmTIgpO+ussygvL2f16tUH3bZIJEIgEMDhcGCxNF2TjRw5ktWrVzcLJaZp4vf7sdvtWK3Wg9pfQ0MDr732WsyyGTNmsHXrVmbOnBmzfPTo0SQkJBzUfg6E2+3mkksuaRM9YA8//DB33nknmzZtonPnzi3dHJEjztbSDRCRw2fTpk18/PHHvPrqq9x4443MmTOHe++990e3C4VCRCKRI9Yun88XDUIul2u/tjEMY7/X3ZeEhASuvPLKmGVz586lqqqq2fI9maaJz+cjLi7ukPYvIm2Lhs9EjiFz5swhJSWFCy+8kEsuuYQ5c+Y0W+e7IaWHH36YRx99lK5du+J0OnniiSc45ZRTALjmmmuiw0rf79FYu3Ytw4YNIz4+ng4dOvCnP/0ppvy7eUNz587ld7/7HR06dCA+Pp7a2tpmc4rOOuss5s+fT0lJSXR/3/VI7G1O0Y4dO7jmmmvo2LEjTqeT7OxsRo0adchDX507d2bkyJEsWLCAAQMGEBcXx9/+9jcAqqurue2228jJycHpdNKtWzemT5/eLEQ+/PDDDB48mLS0NOLi4ujfvz+vvPJKzDqGYdDQ0MDs2bOjx/tdj9zUqVMxDIOioiKuvPJKPB4PGRkZTJkyBdM02bJlC6NGjSIpKYmsrCxmzJjR7Dj8fj/33nsv3bp1w+l0kpOTw69//Wv8fn+zdtxyyy3MmzePXr164XQ66dmzJ++++250nalTp3LnnXcC0KVLl6M+zCjSEtRTJHIMmTNnDmPGjMHhcHDZZZfx5JNP8vnnn0fDzp5mzZqFz+fjhhtuwOl0Mnr0aOrq6vj973/PDTfcwJlnngnA4MGDo9tUVVUxYsQIxowZw7hx43jllVe466676N27NxdccEFM/ffffz8Oh4PJkyfj9/txOBzN2vDb3/6WmpqamOEst9u9z+MbO3Ysa9as4dZbb6Vz586Ulpby/vvvs3nz5kMe3lm3bh2XXXYZN954I9dffz09evSgsbGRoUOH8u2333LjjTeSm5vLxx9/zD333MP27dtj5kL9+c9/5uKLL+aKK64gEAgwd+5cfvrTn/LWW29x4YUXAvDcc89x3XXXMXDgQG644QYAunbtGtOOSy+9lPz8fB588EHmz5/Pf//3f5Oamsrf/vY3zj77bKZPn86cOXOYPHkyp5xyCkOGDAGahiYvvvhiPvroI2644Qby8/P56quvmDlzJkVFRc3mMX300Ue8+uqrTJw4kcTERP7yl78wduxYNm/eTFpaGmPGjKGoqIgXXniBmTNnkp6eDkBGRsYhnWeRVs0UkWPCF198YQLm+++/b5qmaUYiEbNjx47mpEmTYtbbtGmTCZhJSUlmaWlpTNnnn39uAuasWbOa1T906FATMJ999tnoMr/fb2ZlZZljx46NLlu4cKEJmCeccILZ2NgYU8d3ZQsXLowuu/DCC81OnTo129937fyuLVVVVSZgPvTQQ/txNvZtb/vr1KmTCZjvvvtuzPL777/fTEhIMIuKimKW33333abVajU3b94cXfb9Yw0EAmavXr3Ms88+O2Z5QkKCOX78+Gbtuvfee03AvOGGG6LLQqGQ2bFjR9MwDPPBBx+MLq+qqjLj4uJi6nnuuedMi8ViLlmyJKbep556ygTMpUuXRpcBpsPhMDds2BBdtnLlShMwH3vsseiyhx56yATMTZs2NWuvyLFIw2cix4g5c+bQrl07hg0bBjQNkVx66aXMnTuXcDjcbP2xY8ce8FW/2+2OmYvjcDgYOHAg33zzTbN1x48ff1jn5MTFxeFwOFi0aBFVVVWHrd7vdOnShfPPPz9m2csvv8yZZ55JSkoK5eXl0Z/hw4cTDof58MMPY9r3naqqKmpqajjzzDP5z3/+c0DtuO6666KfrVYrAwYMwDRNfv7zn0eXJycn06NHj5jz/vLLL5Ofn8+JJ54Y09azzz4bgIULF8bsZ/jw4TG9VCeddBJJSUl7/W8pcrzQ8JnIMSAcDjN37lyGDRvGpk2bossHDRrEjBkz+OCDDzjvvPNitunSpcsB76djx44YhhGzLCUlhVWrVjVb92Dq/yFOp5Pp06dzxx130K5dO0499VRGjhzJ1VdfTVZW1iHXv7f2rl+/nlWrVu0zPJaWlkY/v/XWW/z3f/83K1asiJnD8/3z9WNyc3Njvns8HlwuV3T4as/lFRUVMW0tLCzcr7bubT/Q9N/ySAROkbZCoUjkGPDvf/+b7du3M3fuXObOndusfM6cOc1C0cH04uzr1nhzL0/2OBJ3bt12221cdNFFzJs3jwULFjBlyhSmTZvGv//9b/r27XtIde+tvZFIhHPPPZdf//rXe90mLy8PgCVLlnDxxRczZMgQnnjiCbKzs7Hb7cyaNYvnn3/+gNqxt3O8P+c9EonQu3dvHnnkkb2um5OTc8B1ihxvFIpEjgFz5swhMzOTv/71r83KXn31VV577TWeeuqpHw0qB9qrcTgc6D67du3KHXfcwR133MH69evp06cPM2bM4B//+Mdhb1vXrl2pr69n+PDhP7jeP//5T1wuFwsWLMDpdEaXz5o1q9m6R+ocd+3alZUrV3LOOecctn20xL8HkZakOUUibZzX6+XVV19l5MiRXHLJJc1+brnlFurq6njjjTd+tK7vHl54NJ9gnJCQQE1NzY+u19jYiM/ni1nWtWtXEhMTm91yfriMGzeOTz75hAULFjQrq66uJhQKAU29LoZhxMzdKi4u3uuTqxMSEo7I+R03bhzffvst//u//9uszOv10tDQcMB1tsS/B5GWpJ4ikTbujTfeoK6ujosvvniv5aeeeioZGRnMmTOHSy+99Afr6tq1K8nJyTz11FMkJiaSkJDAoEGDDvv8oD3179+fF198kV/96leccsopuN1uLrroombrFRUVcc455zBu3DgKCgqw2Wy89tpr7Ny5k5/97GdHpG133nknb7zxBiNHjmTChAn079+fhoYGvvrqK1555RWKi4tJT0/nwgsv5JFHHmHEiBFcfvnllJaW8te//pVu3bo1m2/Vv39//vWvf/HII4/Qvn17unTpclCvZPm+q666ipdeeolf/OIXLFy4kNNPP51wOMzXX3/NSy+9FH0G04Ho378/0PTohJ/97GfY7XYuuuiio/Lkb5GWoFAk0sbNmTMHl8vFueeeu9dyi8XChRdeyJw5c2Im5u6N3W5n9uzZ3HPPPfziF78gFAoxa9asIxqKJk6cyIoVK5g1axYzZ86kU6dOew1FOTk5XHbZZXzwwQc899xz2Gw2TjzxRF566SXGjh17RNoWHx/P4sWLeeCBB3j55Zd59tlnSUpKIi8vj/vuuw+PxwPA2Wefzd///ncefPBBbrvtNrp06cL06dMpLi5uFooeeeQRbrjhBn73u9/h9XoZP378YQlFFouFefPmMXPmTJ599llee+014uPjOeGEE5g0aVJ0/tOBOOWUU7j//vt56qmnePfdd4lEImzatEmhSI5ZeveZiIiICJpTJCIiIgIoFImIiIgACkUiIiIigEKRiIiICKBQJCIiIgIoFImIiIgAek7RPgUCAd577z06d+68z3cEiYiISOsSDocpLi7mvPPOw+FwHNC2CkX78N577+31AXIiIiLS+r355puMHDnygLZRKNqHzp07A/CzxERS1VMk0ibd8vHHAAx+fHALt0REjpZwZZi6uXXRv+MHQqFoH74bMku1Wsm021u4NSJyMPLz8wGwZ+p3WOR4czBTXzTRWkRERASFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERacUqZlYQLAkesfrLppYRrgkfsfqlbbG1dANEREQqZlYQaYiAARanBWdPJwnnJbR0s+Q4o1AkIiKtQvJVydg72QlVhKiZVYM13drSTZLjjEKRiIi0KrY0G/ZcO6HSUMzy4NYg9W/XEy4PY8QZxA+OJ25QXLS88dNGfMt8ROojWNOteK70YIm3ENoZon5+PaHSENYUK+6Rbuwd7NHtAl8HaFzaiBk2iR8cT/zp8QCYQZP69+rxr/VjWA1c/VzED4nHsBhH50TIUadQJCIirUqoPERwc5D4YfEEigK7CyzgHunGlmUjtD1EzbM12HJt2LPt+Fb58H7mxfMzD9YMK6EdIbCC6Tep+UcN7gvcOE50EFgXoPbFWlJvTcWwN4Ubf5GflIkpROoi1MyuwZZtw3GCg8YPGwmXhkm9JbWpnmdrsHqsuPq6WujMyJGmidYiItIq1MypoXxaOTX/qMHZx9ksfNjb27G3t2NYDOwd7Di6OwhtbupN8q3wEX9GPLZ2tqby9nYsTgv+Ij/WDCvOAieGxcCZ78SSYCG4dffk7fgz4rG4LNgybLj6uvCv9jfVudpH/FnxWOIsWJOtxA2Ow7fad/ROiBx16ikSEZFWwXOFB3sn+z7LQ6Uh6t+tb+oFCoMZMqPzjiK1EazJzecgRWoiBEuClE8rjy4zIyaRukj0u9WzezuLp2m4DSBSF4kps3qsMdvJsUehSERE2oT6t+ux59rxXObBsBvUvlILZlOZJclCuLr5rfWWJAuOrg48l3v2WW+4Jow1dVe4qolgcTcNolgSLTFl4ZowlkQNsBzL9F9XRETaBNNvYrgMsEGgJBAz38jVx4V3qZdQaQjTNAluCxLxR3DkOQjtCOEv9GOGTcygSWB9gIhvd49P49JGIr4IofIQvi99OHs6AXD2dNK4uJGIN0K4Joz3Ey+uXppPdCxTT5GIiLQJCecmUPdmHY2LGrF3tePo4YiWuU5yEWmIUPN8DWajiTXDiueKprvPPJd7qF9QT93rdWAFe46dxJzE6LbO7k6qnqhquvvs1HgcXZvqTRiaQP2Ceiofr8SwNN195jzZedSPW44ewzRNs6Ub0RoVFhZSUFDAxORkMu37HuMWkdbr3tJSADLvy2zhlojI0RIsDVL9RDVr164lPz//gLbV8JmIiIgICkUiIiIigEKRiIiICKBQJCIiIgIoFImIyDEsXBWm7L6yQ67H96WP6tnVh94gadUUikREpNWomFlB2R/LMAO7b4w2AyblD5RTMbOiBVsmxwOFIhERaVWsiVb8X/uj3/3r/NGnTB8IM6wnzsiB0cMbRUSkVXH2duJf5cd1UtPTo/2r/DhPcuL7sullrA0fNuBb7sP0mlgzrSSOTMSW1fTnrGJmBXGnxOFb4cMMmiRPSI6p2/uZF+9/vCRflYzhNKh/r55AYQAMiBsYR/yZ8UBT71Tdm3UEigJYU6w4ujuQY59CkYiItCr2LnZ8X/qINDS9iiO0M0T8mfHRUGRLt5FyQwqGy6BxcSO1r9WSelNqdHt/oZ/kCU2hJ1If+zoP/1d+kscnY4m3UDe/DtNrknprKhF/hJrnarBmWnH2cNKwqIFIfYTU21OJ1O0qS2v+wlk5tmj4TEREWhXDMHAWOPGv9uNf7ceZ7wRjd7mzwIklwYJhNYg/M57wzjCmf/dQWdypcVjcFgz77o0aFjfgX+PHM77p1R+maeL70kfCeQkYTgNrkpW4AXH41zYN2/nX+okfEo/FZcGWYcN1st55djxQT5GIiLQ6zt5O6t+pB8B9gRt2d/jgXe7F+6mXSO3uhRFvBKuzqSfHkvS9630TvJ96cf/EjSWuqcxsMCEEVX+tilnPltv0ZzFSF8Hq2d0zZPFYYOvhPEJpjRSKRESk1bF3sGN6zejn4JYg0HSLff279SRPSMaWbYMwlP+xHPaYU23s2a3UtADPlR5qnq/BkmjB0dmBEW+AFVInpWKJbz5oYkm0EK4JY01tCkaRmkizdeTYo+EzERFplZIuTSLp0qSYZWbAxDAMLAkWiEDDwob9qsvewU7ST5OofbmW4LfBprfe93FR/149EW8EM2ISKgsR3NoUvpwFThqXNBLxRQiVh/Ct9B3245PWRz1FIiLSKtkym/+JsrWz4ervourJKgy7QfzQeNjP+c+Ozg4SRyVS+0Itnqs9uM930/BBA1VPVmH6TaxpVhLOTgAgYWgCdW/WUTmzEmuKFdfJrmhgkmOXYZqmHuSwF4WFhRQUFDAxOZlMu72lmyMiB+He0lIAMu/LbOGWiMjREiwNUv1ENWvXriU/P/+AttXwmYiIiAgKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiABga+kGiIgcaaX3lrZ0E0TkKCksLKTgiYKD2lY9RSIiIiKop0hEjgO97uvV0k0QkaPEW+o96G3VUyQiIiKCQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkctjMrKigJBhs6WaIiMhBsrV0A0T2xx/LyqKfg4B9j7KbU1NJtloPeR+zqqvp53Jxsst1yHWJiEjbo1AkbcJvMzKin+8vK2Niaiop3wtCpmliAhbDOMqtExGRY4FCkbRpr9XW4jAMysJhtgaD3JSSQgiYX19PaShEitXKSLebDvamvqUPGxpY7vPhNU0yrVZGJiaSZbOxuKGBzcEgW4NB3qqr48z4eIYkJFAcCLCgoYHKcJgsm41RiYmk7gpj6/1+3q6vx2uanBoX14JnQUREDgeFImnzVvv9XOXx0M5mI2SaPF5VxQVuNyc6HKwLBHixtpZbU1OxGwbpNhs3pKTgMgwWNzbyWm0tN6WmMjQhgW+CwZjhs5pwmJdqa7k0KYkcu53PfT5eqa3lhpQUGiIRXq6rY2xiIl0dDv7d0EBtJNLCZ0JERA6FJlpLm1fgdNLebsdqGBQFAmRYrRQ4nVgMg3ynkwSLha27JkAX7PpuNQzOjI9nZziM3zT3Wu8qv58TnU46ORxYDINBcXFUh8NUhcOsDwTIttno4XRiMwzOSkhAg3YiIm2beoqkzUuy7M72NZEIJcEg08rLo8sipkndrl6c5V4vn3q9Mb063kgE514mateEw6zy+Vjj90eXhXfVVReJ4Nljvw7DIN6iawwRkbZMoUiOKUkWC10dDi73eJqVVYXDvFtfz4TkZLJtNsLAH8vL+a6f6Ps9PUkWC/3i4viJ292srspwmA17BKugadKo4TMRkTZNl7ZyTMlzONgRClHo9xM2TYKmyfpAAF8kQsA0MQyDBIuFCLCwoSFm2wSLhepwOPq9t8vFWr+fkkCAiGnij0SivUbdHQ62h0IU+f2ETJNFDQ3sfRBORETailYXip555hmSk5MPa52LFi3CMAyqq6sPa73S+rgsFi73eFjm9fJQRQWPVlSw3OsFoJ3NRn+Xiyerqni0spIUq5U9B80GxsWxwudjWnk5SxobSbFauSQxkfcaGpheUcHjVVWs2xWKEiwWxiYm8k59PQ9XVGA3jJhhPBERaXsOaPhswoQJzJ49mxtvvJGnnnoqpuzmm2/miSeeYPz48TzzzDOHs40iMabs8cyi0UlJzcqzbDbG7yNYn+92c/4ew2Gn7HErfSe7nUlpaTHrd3Y4uN7h2GtdPZxOejid0e9nJSTsV/tFRKR1OuBL25ycHObOnYt319U3gM/n4/nnnyc3N/eQGhPUKxJERESkhRxwKOrXrx85OTm8+uqr0WWvvvoqubm59O3bN7rs3Xff5YwzziA5OZm0tDRGjhzJxo0bo+XFxcUYhsGLL77I0KFDcblczJkzp9n+ysrKGDBgAKNHj8bv9xOJRJg2bRpdunQhLi6Ok08+mVdeeSVmm7fffpu8vDzi4uIYNmwYxcXFB3qYIiIicpw5qEkQ1157LbNmzYp+f/rpp7nmmmti1mloaOBXv/oVX3zxBR988AEWi4XRo0cT+d4dOnfffTeTJk2isLCQ888/P6Zsy5YtnHnmmfTq1YtXXnkFp9PJtGnTePbZZ3nqqadYs2YNt99+O1deeSWLFy+ObjNmzBguuugiVqxYwXXXXcfdd9/9o8fk9/upra2N/tTX1x/MqREREZE26qBuyb/yyiu55557KCkpAWDp0qXMnTuXRYsWRdcZO3ZszDZPP/00GRkZrF27ll69ekWX33bbbYwZM6bZPtatW8e5557L6NGjefTRRzEMA7/fzwMPPMC//vUvTjvtNABOOOEEPvroI/72t78xdOhQnnzySbp27cqMGTMA6NGjB1999RXTp0//wWOaNm0a991338GcDhERETkGHFRPUUZGBhdeeCHPPPMMs2bN4sILLyQ9PT1mnfXr13PZZZdxwgknkJSUROfOnQHYvHlzzHoDBgxoVr/X6+XMM89kzJgx/PnPf8bY9YLPDRs20NjYyLnnnovb7Y7+PPvss9GhucLCQgYNGhRT33cB6ofcc8891NTURH+WLVu23+dD2qaqcJj7ysoOuZ4vfT5m685GEZE276Af3njttddyyy23APDXv/61WflFF11Ep06d+N///V/at29PJBKhV69eBAKBmPUS9nLHjtPpZPjw4bz11lvceeeddOjQASA6pDV//vzosj23ORROpzOmDvdeHtgnLWtmRQWNkQh3pqfj2BWUA6bJwxUVxBkGt3/vzjEREZEDcdAPVhkxYgSBQIBgMNhsLlBFRQXr1q3jd7/7Heeccw75+flUVVXtf6MsFp577jn69+/PsGHD2LZtGwAFBQU4nU42b95Mt27dYn5ycnIAyM/Pb9bL8+mnnx7sYUork2i18vUer91Y5/fjPojnA4X38b4zERE5fh10T5HVaqWwsDD6eU8pKSmkpaXxP//zP2RnZ7N58+b9muz8/frnzJnDZZddxtlnn82iRYvIyspi8uTJ3H777UQiEc444wxqampYunQpSUlJjB8/nl/84hfMmDGDO++8k+uuu47ly5fruUnHkN5OJ6v8fk7a9Sb7VX4/JzmdfOnzAfBhQwPLfT68pkmm1crIxESybE3/zGdWVHDKrgc0Bk2TCd97ltFnXi//8Xq5KjkZp2HwXn09hYEABk0PdjwzPh5o6p16s66OokCAFKuV7vt4jpGIiLQth/QI3qSkJJL28vA8i8XC3LlzWb58Ob169eL222/noYceOuD6bTYbL7zwAj179uTss8+mtLSU+++/nylTpjBt2jTy8/MZMWIE8+fPp0uXLgDk5ubyz3/+k3nz5nHyySfz1FNP8cADDxzKYUor0sVupzQUoiESoSESYWcoxAl7hJJ0m40bUlK4Ky2NExwOXqutjdm+0O9nQnIyt6Smxixf2tjIlz4f45OTcVssvFdfj9c0uTU1leuTk1nl80WfZr2ooYH6SITbU1MZm5jIyl2BTERE2jbDNDWOsDeFhYUUFBQwMTmZTLu9pZsjNPX0jElKotDvJ2VX72RlOEwvp5NXamubzSkKmiZ/LC/nnvR0nIbBzIoKhick0HtXL1NVOMxfKis5Kz6edYEAV3k8xFksmLu2+2VqKkm79vOZ18u2YJDRSUk8WlHBqMREuuwKYx80NLA1GNznU7Sl5dxbWgpAr/t6/ciaInKs8JZ6+eaJb1i7di35+fkHtO1BD5+JtJTeTifv7Jp0f4HbzZ5Pvlru9fKp10vtHs/D8kYiOHeFm++/n8wEPvV6+YnbTdyusgbTJAT8dY95cCaQu2sYri4SwbPHkLHHYmHr4Ts8ERFpIQpF0uZ0sNvx7urg7GC3s2XX62GqwmHera9nQnIy2TYbYeCP5eWxb6/fddda9CtwpcfD8zU1JFosdHY4iDcMrMCk1FTi9zKJO9FioSYcJnVXMKr53gNJRUSkbdJrvaVNujQpiUu/N58tYJoYhkGCxUIEWNjQsF91dbDb+WlSEi/X1vJtMIjFMOjjcjXNK4pEiJgmZaEQW3eFrwKnkyWNjfgiEcpDIc0pEhE5RqinSNqkTFvzf7rtbDb6u1w8WVWF3TAYGh+PdS/b7k1nh4NRiYm8UFvL1R4P57vdfNDQwJNVVfhNkzSrlbN3PVNraEICb9bVMbOykhSrlZNdrmhgEhGRtksTrfdBE61F2j5NtBY5/hzKRGsNn4mIiIigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIq1Y0cwiGkoajlj9a6auIVgTPGL1S9tia+kGiIiIFM0sItQQwjAMLE4LST2TyDovq6WbJccZhSIREWkVOl3ViYROCfgr/BTPKsaZ7mzpJslxRqFIRERaFWeak/jcePyl/pjljVsb2f72dgLlAaxxVtIGp5E2KC1aXvFpBZXLKgnVh3CmO8m9MhdbvA3fTh/b52/HV+rDkeKg/cj2xHWIi25X+3UtFUsrMMMmaYPTSD89HYBIMMKO93ZQt7YOw2qQ3C+ZjCEZGBbj6JwIOeoUikREpFXxl/tp3NxI5rBM6orqossNi0H7ke1xZbnwbfdR/Gwx8bnxxGXHUb2qmorPKsj9WS7ODCe+HT4Mq0HYH6bkHyVkX5BN4omJ1K2rY8uLW+h2azcs9qZptXVFdXSd2JVQXYji2cW4sl24T3BT9mEZ/lI/3W7p1lTPsyXYPXZS+qa01KmRI0wTrUVEpFXYPGczhdMKKflHCcl9kknumxxTHtc+jrj2cRgWg7gOcbi7u2nc3AhA9Ypq0s9Ix9XO1VTePg6r00p9UT3ODCdJBUkYFoOk/CSsCVa8W73RejPOyMDqsuLMcJLcN5na1bUA1KyuIfOsTKxxVhzJDtIGp1GzuuaonQ85+tRTJCIirULuFbkkdErYZ7mv1MeOd3fg2+HDDJuYITM67yhYG8SR7Gi2TbAmSGNJI4XTCqPLzIhJsG73HWd2jz3ms2+nD4BQXahZWagudPAHKK2eQpGIiLQJ29/eTkJuArmX5WKxW9jyyhYwm8rsSXaC1c1vrbcl2XB3dZN7ee4+6w3WBHGkOqKf7e6mIGRLtDUrsyXqz+axTMNnIiLSJkT8ESwuC4bNoKGkgfqi+mhZcp9kypeW4yv1YZom3m1ewv4wiXmJeHd4qS2sxQybRIIR6tbXEfaFo9uWLy0n7AvjL/dT/WU1ST2TAPD09FC2uIywN0ywJkjFJxV4enmO+nHL0aPIKyIibUK7c9ux7c1tlC0qI6FrAok9EqNlySclE24Is/n5zYQbwzgznORe0XT3WafLO7FjwQ62vb4NrBCfE098Tnx0W3d3Nxuf2EgkHCHt1DTcXd0AZAzNYMeCHax/fD2GxSClXwrJJycf7cOWo8gwTdNs6Ua0RoWFhRQUFDAxOZlMu/3HNxCRVufe0lIAet3Xq4VbIiJHi7fUyzdPfMPatWvJz88/oG01fCYiIiKCQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiLQBRTOLKJpZhBnefcP0tje3UbqwtAVbJccaPafoR9zy8ccHfEufiIgcfhF/hKovq0gdkNrSTZFjlELRj3h85Egy4+JauhkichDuXb26pZsgh1HqaamULyknpW8KhtWIKatcVkn5x+VEAhHcXd1kX5iN1WWlYVMD297YhudkD5WfVmJxWmh/UXvc3Zoe0BhqDLHj7R3Uf1OPxWEh86xMkvskt8DRSWug4TMREWkTEjonYPfYqVpRFbO8fmM9ZR+WkXt5Lnm35REJRdjxzo5oeaA6gGEz6HFnD9LPTGfbm9uiZd++9i22RBt5v8qj0xWd2Pmvnfh2+I7aMUnrolAkIiJtRsZZGZQvKY+ZW1SzuoaU/im4Ml1YHBbandOOmjU1fPfCBovDQvrgdAyrQfJJyQRrgk3vM6sL0lDcQLvh7bDYLDgznHh6e6gtrG2pw5MWpuEzERFpM9wnuLEn2qleUR1dFqoLxbzLzO6xY4ZMwt6ml77a4m0YlqbhNoujqS8gEogQrAtihky+/tPXu3dggqe3Xvp6vFIoEhGRNiXjrAy2vbmNhC4JANgSbQRrgtHyYE0Qw2ZgjbP+YD32JDsWh4UT7z4RwzB+cF05Pmj4TERE2hR3Vzc2t426r+sA8PT0ULW8Cn+Zn0ggQum/S0nqmfSjQceeZCe+Yzyl/y4lEohghk2827z4SjWn6HilniIREWlzMs/KpOQfJQC4u7lJPzOdkjklRPy77j4bkb1f9XQY24EdC3aw/s/rMcMmzkwnWSOyjmTTpRVTKBIRkVYv7/a8mO/ubm56Tu0Z/Z42KI20QWnNtkvokkD3Sd1jlu25nS3eRsfRHQ9za6Wt0vCZiIiICApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIofNzKIiShoaWroZIiJykGwt3QCR/fHHwsLo52Akgt0wwDAAuLlrV5IdjkPex6xNm+iXksLJycmHXJeIiLQ9CkXSJvw2Pz/6+f61a5nYrRsp3wtCpmliApZdYUlERORAKBRJm/bat9/isFgo8/vZ2tjITV27EjJN5m/fTqnPR4rDwcj27ekQFwfAh2VlLK+qwhsOk+l0MrJ9e7JcLhaXlrK5sZGtXi9vbd/OmenpDMnIoLihgQU7dlAZCJDlcjGqQwdSd4Wx9XV1vL1jB95wmFNTU1vyNIiIyGGgUCRt3uqaGq7q1Il2LhehSITHN2zgguxsTkxMZF1dHS9u2cKt3bpht1hIdzq54YQTcFmtLC4r47Vvv+Wmrl0ZmpnJNw0NMcNnNcEgL23ZwqU5OeTEx/N5ZSWvbNnCDV270hAK8fLWrYzt2JGuCQn8u7SU2mCwZU+EiIgcEk20ljavICmJ9nFxWA2Dovp6MpxOCpKSsBgG+UlJJFitbPV6o+sm2GxYDYMz09PZ6fPhD4f3Wu+q6mpOTEqiU0ICFsNgUFoa1cEgVYEA6+vryXa56JGYiM1i4azMTAwN24mItGnqKZI2L8m2+59xTTBISWMj0/aYmB0xTep29eIsr6zk08rKmF4dbziM02ptVm9NMMiq6mrW1NREl4VNk7pQiLpgEI/dHl3usFiI30sdIiLSdigUSdu3Rw9Nks1GV7eby3Nzm61WFQjw7o4dTOjShWyXi7Bp8sfCQszvqvne+kl2O/1SUvhJdnazuioDATbU10e/ByMRGvfR4yQiIm2Dhs/kmJKXmMgOr5fC2lrCpkkwEmF9XR2+cJhAJIJhGCRYrURMk4WlpTHbJthsVO/Rg9Tb42FtbS0lDQ1ETBN/OBztNerudrPd56Ooro5QJMKisjJM00RERNquoxqKJkyYgGEYPPjggzHL582bp/kYcli4rFYu79SJZZWVPLRuHY8WFbG8qgqAdi4X/VNSeHLjRh5dv54UhwPrHv/uBqamsqK6mmmFhSwpKyPF4eCSjh15b+dOpn/9NY9v2MC6ujqgKUCN7diRd7Zv5+F167AbBkl7DKeJiEjbc9SHz1wuF9OnT+fGG28kJSXlaO9ejgFTCgqin0d36NCsPMvlYnznznvd9vysLM7Pyop+P2WPW+k7JSQwqXv3mPU7JyRw/Qkn7LWuHomJ9EhMjH4/KzNzv9ovIiKt01EfPhs+fDhZWVlMmzZtn+v885//pGfPnjidTjp37syMGTNiyjt37swDDzzAtddeS2JiIrm5ufzP//xPzDpbtmxh3LhxJCcnk5qayqhRoyguLj4ShyQiIiLHgKMeiqxWKw888ACPPfYYW7dubVa+fPlyxo0bx89+9jO++uorpk6dypQpU3jmmWdi1psxYwYDBgzgyy+/ZOLEidx0002sW7cOgGAwyPnnn09iYiJLlixh6dKluN1uRowYQSAQ2Gu7/H4/tbW10Z/6PSbRioiIyLGvRSZajx49mj59+nDvvfc2K3vkkUc455xzmDJlCnl5eUyYMIFbbrmFhx56KGa9n/zkJ0ycOJFu3bpx1113kZ6ezsKFCwF48cUXiUQi/N///R+9e/cmPz+fWbNmsXnzZhYtWrTXNk2bNg2PxxP9GThw4GE/bhEREWm9Wuzus+nTpzN79mwK93ieDEBhYSGnn356zLLTTz+d9evXE97jlueTTjop+tkwDLKysijddTfRypUr2bBhA4mJibjdbtxuN6mpqfh8PjZu3LjX9txzzz3U1NREf5YtW3a4DlVERETagBZ7TtGQIUM4//zzueeee5gwYcIBb2//3p0+hmEQiUQAqK+vp3///syZM6fZdhkZGXutz+l04nQ6o9/dbvcBt0lERETarhZ9eOODDz5Inz596NGjR3RZfn4+S5cujVlv6dKl5OXlYd3PJwb369ePF198kczMTJKSkg5rm6V1mllUBMAvu3eP3mb/5rZtuG02humuMBER2Q8t+vDG3r17c8UVV/CXv/wluuyOO+7ggw8+4P7776eoqIjZs2fz+OOPM3ny5P2u94orriA9PZ1Ro0axZMkSNm3axKJFi/jlL3+518ndcmzwRyJ8ueuZRCIiIgeqxV/z8Yc//IEXX3wx+r1fv3689NJL/P73v+f+++8nOzubP/zhDwc0xBYfH8+HH37IXXfdxZgxY6irq6NDhw6cc8456jk6hp2WmsqS8nL6pqTEPJQRYFllJR+XlxOIROjqdnNhdjYuq5VNDQ28sW0bJ3s8fFpZidNi4aL27em2a/i0MRTi7R07+Ka+HseuF7/2SU5ugaMTEZEj7aiGou/fVg9Nzxzy+/0xy8aOHcvYsWP3Wc/enje0YsWKmO9ZWVnMnj37YJopbVTnhAQ2NjSwoqqK/ns8lHFjfT0flpVxdadOJDscvLZ1K+/s2BF98GN1IIDNMLizRw++rKrizW3buD0vD4DXvv2WdKeTX+XlURUMMru4mCyXiyyXq0WOUUREjhy9+0yOKWdlZLCkvJzwHu8hW11TQ/+UFDJdLhwWC+e0a8eamprou8ocFguD09OxGgYnJSdTEwziDYepCwYpbmhgeLt22CwWMpxOens8FNbWttThiYjIEdTiw2cih9MJbjeJdjsrqqujy+pCIXLi46PfPXY7IdPEu+sRD/E2G5Zdw20OS9N1QiASoS4YJGSa/Onrr6PbmjS9KFZERI49CkVyzDkrI4M3t22jS0ICAIk2GzXBYLS8JhjEZhjE/cjdjEl2Ow6LhbtPPFEvLBYROQ5o+EyOOV3dbtw2G1/veqN9T4+H5VVVlPn9BCIR/l1aSs+kpB8NOkl2Ox3j4/l3aSmBSISwabLN66XU5zsahyEiIkeZeorkmHRWZib/KCkBoJvbzZnp6cwpKcG/6+6zEdnZ+1XP2A4dWLBjB39ev56waZLpdDIiK+tINl1ERFqIQpEcE767W+w73dxupvbsGf0+KC2NQWlpzbbrkpDApO7dY5btuV28zcbojh0Pc2tFRKQ10vCZiIiICApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAGBr6Qa0dre89Rb5+fkt3QwRERE5wtRTJCIiIoJ6in7U41dfTabb3dLNEJGDcO/ChS3dBBFpQxSKROSYt/re1S3dBBE5SgoLCyl4ouCgttXwmYiIiAjqKRKR48Cw+4a1dBNE5CipL60/6G3VUyQiIiKCQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkctjM/OQTSqqrW7oZIiJykGwt3QCR/fHHDz+Mfg5GItgtu/P8zQMHkuxyHfI+Zn35Jf2yszk5K+uQ6xKRw+OTmZ+QPyaf5E7JR6T+RVMXcertp+LyHPr/Q6TtUyiSNuG3Q4ZEP9+/eDETTzmFlLi4mHVM08QELIZxlFsnIofqk5mfEGwIggE2p42Mnhl0Pa9rSzdLjjMKRdKmvVZYiMNqpayxka21tdw0YAChSIT569dT2tBAisvFyLw8OiQlAfBhSQnLt23DGwqRmZDAyLw8stxuFhcXs7mmhq21tbxVVMSZnToxpFMniqurWbBhA5VeL1luN6NOPJHUXWFsfUUFb69fjzcU4tSOHVvyNIgcE0666iSSOyXTWNHIilkriE+Pb+kmyXFGoUjavNWlpVx18sm0S0ggFInw+LJlXNC9Oyemp7OuvJwX16zh1oEDsVutpMfHc0P//rhsNhaXlPBaYSE3nXIKQzt35puqqpjhsxqfj5fWrOHSnj3J8Xj4/NtveWXtWm7o35+GQICX165lbH4+XVNT+femTdT6/S18JkSODfFp8XhyPTSUNsQsr91ay/q319NY3ogtzkbO4Bw6Dtp9QbL10618u+xbAvUB4tPjOenKk7DH26nfWc/6+etpKG3AleIib2QeSR2SotuVf13OlqVbiIQj5AzOIff0XADCwTAb39tI2doyLFYL2f2y6TSkE4ZFvdHHKk20ljavICOD9omJWC0WiioqyEhIoCAjA4thkJ+RQYLdztba2ui6CQ4HVouFM3Nz2dnQgD8U2mu9q3bu5MT0dDolJ2MxDAZ17Ei1z0eV18v6ykqy3W56pKdjs1g4q3NnDA3biRwWjeWN1GyuwZ3ljlluWAzyRuZxxt1n0HNcTzb9exN12+sA2LlqJ1s/20rPS3tyxt1nkDcyD8NqEPKHWPWPVXQ8tSOn//p0Og3pxJoX1xAOhqP1VhRVcMrEU+gzoQ9bP9lK1TdVAJR8WEJDaQMDbxlI32v7snPVTnas3HH0ToQcdeopkjYvyemMfq7x+ymprmbakiXRZRHTpC4QAGD5tm18unVrTK+ONxTCaWv+q1Dj97Nq507WlJZGl4V31VXn9+PZY3K3w2ol3m4/rMclcrz5as5XYIA9zk5Wnyyy+2ZT8mFJtDyxfWL0c1KHJNK6p1GzuYbE7ER2rNhB7hm5uNu5Y9bd+dVOEjISyCjIACAjP4OSD0uo3VpLSpcUAHLPyMXmsmFz2cjqm0Xp6lJSTkihdHUpPS7ugT3Ojj3OTs7gHEpXl5LdN/tonRI5yhSK5JiS5HTSNTWVy3v3blZW5fXy7oYNTOjTh+zERMKRCH9csgTTNAGa9fQkOZ30y87mJ927N6ur0utlQ2Vl9HswHKYxGDzMRyNyfOl9Re8fvMusobSBDe9uoH5HPZFwhEgoEp135K/140pufgeZv8ZPdUk1S6btvlAyIyaBukD0+553nrk8Lhp2Ng3bBeoCMWVOjzNmOzn2KBTJMSUvLY1/ffMNhWVl5KWlETFNiquryfF4CITDGIZBgsNBxDRZWFwcs22C3U61zxf93jszk79/+SU9MzLI8XgIhsNsqKykZ2Ym3VNTeXv9eooqKjghJYVFxcXRcCUiR8b6t9fjyfXQ67JeWO1W1r6yNvp750xy4qv2NdvGmeQktWsqvS9vfqH0HV+Nj7jUuOhnh9sBgCPREVPmr/HjSHQc7sOSVuSIzykqKyvjpptuIjc3F6fTSVZWFueffz5Lly490ruW45DLZuPy3r1Z9u23PPTxxzz66acs374dgHZuN/2zs3ny88959NNPSXG5sO7ROzSwQwdW7NjBtCVLWFJSQkpcHJcUFPDexo1M/+gjHl+2jHUVFQAkOByMzc/nnfXrefjjj7FbrTHDeCJy+IX8IWwuGxabheqSaiqKKqJlWX2y2LJ0Cw2lDZimSd22OkL+EGl5adTvqKessIxIOEI4GKZifQUh3+65hFuWbiHkC9FY3siOL3eQ0bNpqC2zZyYli0sIeoP4anxs+WQLmb0yj/pxy9FzxHuKxo4dSyAQYPbs2Zxwwgns3LmTDz74gIqKih/fWGQvpgwdGv08Oj+/WXmW2834Pn32uu353bpxfrdu0e+ndOgQ/dwpOZlJp54as37n5GSu799/r3X1SE+nR3p69PtZnTvvT/NF5CB1Pbcr695cR/GiYlK6ppDWIy1a1u6kdgQaAnz1/FcEG4PEZ8Rz0hVNd5/1vrw3GxZsYN3r6zCsBp4cD54cT3Tb1O6pfP7E50TCETqe2pHUrqkAdBraiY0LNrLs8WUYFoPsftlknayHux7LDPMI9vlXV1eTkpLCokWLGLrHH7I9bd68mVtvvZUPPvgAi8XCiBEjeOyxx2jXrh0AU6dOZd68efzyl79k6tSpVFZWcvXVV/PYY48xY8YMHnnkESKRCJMmTeK3v/1tzL4nT57M66+/jt/vZ8CAAcycOZOTTz55v9peWFhIQUEBEwcMINPt/vENRKTVuXfhQgCG3TeshVsiIkdLfWk9XzzxBWvXriV/LxfOP+SIDp+53W7cbjfz5s3Dv5dnuEQiEUaNGkVlZSWLFy/m/fff55tvvuHSSy+NWW/jxo288847vPvuu7zwwgv8/e9/58ILL2Tr1q0sXryY6dOn87vf/Y7PPvssus1Pf/pTSktLeeedd1i+fDn9+vXjnHPOoXKPybF78vv91NbWRn/q6+sP78kQERGRVu2IhiKbzcYzzzzD7NmzSU5O5vTTT+c3v/kNq1atAuCDDz7gq6++4vnnn6d///4MGjSIZ599lsWLF/P5559H64lEIjz99NMUFBRw0UUXMWzYMNatW8ejjz5Kjx49uOaaa+jRowcLd10VfvTRRyxbtoyXX36ZAQMG0L17dx5++GGSk5N55ZVX9trWadOm4fF4oj8DBw48kqdGREREWpkjPtF67NixbNu2jTfeeIMRI0awaNEi+vXrxzPPPENhYSE5OTnk5ORE1y8oKCA5OZnCwsLoss6dO5OYuPv5FO3ataOgoADLHi8FbdeuHaW7niezcuVK6uvrSUtLi/ZWud1uNm3axMaNG/faznvuuYeamproz7Jlyw73qRAREZFW7Kjcku9yuTj33HM599xzmTJlCtdddx333nsvd9xxx35tb//eQ/EMw9jrskgkAkB9fT3Z2dksWrSoWV3Jycl73YfT6cS5x91Dbs0jEhEROa60yHOKCgoKmDdvHvn5+WzZsoUtW7ZEe4vWrl1LdXU1BQUFB11/v3792LFjBzabjc66I+iYV1xdzfsbN1Le2IjFMGjndjOqRw9Sdr24VUREZH8c0VBUUVHBT3/6U6699lpOOukkEhMT+eKLL/jTn/7EqFGjGD58OL179+aKK67g0UcfJRQKMXHiRIYOHcqAAQMOer/Dhw/ntNNO47/+67/405/+RF5eHtu2bWP+/PmMHj36kOqW1sUXCjF39Wou7tGDE9PTCYbDfFNVpfeQiRxDqour2fj+RhrLGzEsBu52bnqM6kFcii585PA6oqHI7XYzaNAgZs6cycaNGwkGg+Tk5HD99dfzm9/8BsMweP3117n11lsZMmRIzC35h8IwDN5++21++9vfcs0111BWVkZWVhZDhgyJ3uovx4aKxkashkFBRtPD1pw2G/m7PgfDYd7buJG1ZWVYLRb6ZWczpFMnLIbBwk2bqPB6CUcibKisJMvtZlzPniwqLmZ1aSkpcXGM69mT1F29TTvr65m/fj2lDQ2kuFyMzMujQ1LSPtslIodHyBdi9dzV9Li4B+knphMOhqn6Rhc+cmQc0ecUtWV6TlHb4AuFePTTT8lPT6dXZiYdkpJw7Xq56wfffMPmmhp+1qsX/nCYZ1eu5MzcXPpmZ7Nw0yY+3rKFK086iQ5JSTz/1VdUer0MP+EE8tPTeauoiIhpMjo/H38oxOPLlnFB9+6cmJ7OuvJy3tmwgVsHDsRutbbwGZAfoucUtX2139by1fNfcfqdpzcrCwfDbHxvI2Vry7BYLWT3y6bTkE4YFoNNCzfhrfASCUeo3FCJO8tNz3E9KV5UTOnqUuJS4ug5rmf0FR71O+tZP389DaUNuFJc5I3MI6mDLnzaolb7nCKRI81ls3FNnz6ETZNXCwt5aOlS/rl2Lf5QiNWlpZzVuTNxdjvJLheDc3JYvccb77ukpNApORmbxUJ+ejoOq5VemZlYLRYKMjLYsetZVUUVFWQkJFCQkYHFMMjPyCDBbmdrbW1LHbbIcSM+LR4zbPL1619TubEy5vUcJR+W0FDawMBbBtL32r7sXLWTHSt3RMvL15XT8dSOnHHXGVjtVr78+5ekdEnh9F+fjjvbTfHiYqDp9SGr/rGKjqd25PRfn06nIZ1Y8+IawsHw0T5caWF6Iay0ee3cbsbsuhrYVlfHS2vW8GFJCXWBAB7X7jdce5xO6gK733CdsMcdjDaLpdn3QLjpf4g1fj8l1dVMW7L7LdsR04ypS0SODJvLRp9r+rBl6RYKXy0k5AuRUZBB3sg8SleX0uPiHtjj7Njj7OQMzqF0dSnZfbMBSOmSQnKnZADS89PZ9vm26LvLMgoy+OZf3wBQUVRBQkYCGQVNQ+8Z+RmUfFhC7dZaUrqkHP2DlhajUCTHlPaJieSnp1Pa0ECiw0GNzxedF1Tj95PoOPA3XCc5nXRNTeXy3vt+y7aIHDnudm7yxzRd+NRtq2PNS2so+bCEQF0Al2f3hY/T4yRQt/tixZ6w+0LHYrM0+x4ONF34+Gv8VJdUs2Ta7gsfM2LG1CXHB4UiadPKGhrYUFlJz8xMkpxOyhsbKaqooHe7drRzu1lcUkKW200gHOaTLVs4Izf3gPeRl5bGv775hsKyMvLS0oiYJsXV1eR4PNH5SyJydCS2TyQ9P52G0gYciQ58Nb7ovCB/jR9H4oFf+DiTnKR2TaX35brwOd7p/+jSpjltNrbU1rJ0yxb8oRBxdjsFGRmckZuLaZos2LiRx5ctw2IY9MvO5uSsA3/Dtctm4/LevVmwYQOvr1uH1TDI8XjI8Xh+fGMROSQNZQ1Ubqgks2cmziQnjeWNVBRV0K53O9zt3JQsLsGd5SYcCLPlky3knnHgFz5peWl8869vKCssIy0vDTNiUl1cjSfHg82lP5PHE/3XljYtyelkXM+e+ywfmZfHyLy8ZsuHdekS871vdjZ9s7Oj37ukpDDp1FOj37Pcbsb36XPoDRaRA2Jz2qjdUsuWpVsI+UPY4+xkFGSQe0bThc/GBRtZ9vgyDItBdr9ssk4+8Asfm8tG78t7s2HBBta9vg7DauDJ8eDJ0YXP8UahSEREWi1nkpOe4/Z94ZM3Mo+8kc0vfLoMi73wye6bHZ2ADU2TsE+dtPvCx53lps/4PofeYGnTdEu+iIiICApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICgK2lG9BahcNhACq93hZuiYgcrMLCQgDqS+tbuCUicrR4K5v+bn/3d/xAKBTtQ3FxMQBz16xp2YaIyEF7oqCgpZsgIi2kuLiYXr16HdA2hmma5hFqT5sWCAR477336Ny5M1artaWbI0dAfX09AwcOZNmyZbjd7pZujogcZvodPz6Fw2GKi4s577zzcDgcB7StQpEct2pra/F4PNTU1JCUlNTSzRGRw0y/43KgNNFaREREBIUiEREREUChSI5jTqeTe++9F6fT2dJNEZEjQL/jcqA0p0hEREQE9RSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJHDETJkzgv/7rv1q6GSJtzjPPPENycvJR3edZZ53FbbfddlT3Ka2PQpG0SoZh/ODP1KlTj8h+FWREjo4JEyZEf58dDgfdunXjD3/4A6FQqKWbJscxvRBWWqXt27dHP7/44ov8/ve/Z926ddFle77HyDRNwuEwNpv+OYu0JSNGjGDWrFn4/X7efvttbr75Zux2O9nZ2S3dNDlOqadIWqWsrKzoj8fjwTCM6Pevv/6axMRE3nnnHfr374/T6eSjjz4iEokwbdo0unTpQlxcHCeffDKvvPJKtM5wOMzPf/7zaHmPHj3485//HC2fOnUqs2fP5vXXX49ewS5atAiALVu2MG7cOJKTk0lNTWXUqFEUFxfH1P2rX/2K5ORk0tLS+PWvf40eASbyw5xOJ1lZWXTq1ImbbrqJ4cOH88YbbzRbb+PGjYwaNYp27drhdrs55ZRT+Ne//hWzjt/v56677iInJwen00m3bt34+9//Hi1fvXo1F1xwAW63m3bt2nHVVVdRXl4eU0coFOKWW27B4/GQnp7OlClTYn6Pq6qquPrqq0lJSSE+Pp4LLriA9evXH+azIi1JoUjarLvvvpsHH3yQwsJCTjrpJKZNm8azzz7LU089xZo1a7j99tu58sorWbx4MQCRSISOHTvy8ssvs3btWn7/+9/zm9/8hpdeegmAyZMnM27cOEaMGMH27dvZvn07gwcPJhgMcv7555OYmMiSJUtYunQpbrebESNGEAgEAJgxYwbPPPMMTz/9NB999BGVlZW89tprLXZuRNqiuLi46O/Unurr6/nJT37CBx98wJdffsmIESO46KKL2Lx5c3Sdq6++mhdeeIG//OUvFBYW8re//S3ao1xdXc3ZZ59N3759+eKLL3j33XfZuXMn48aNi9nP7NmzsdlsLFu2jD//+c888sgj/N///V+0fMKECXzxxRe88cYbfPLJJ5imyU9+8hOCweAROiNy1JkirdysWbNMj8cT/b5w4UITMOfNmxdd5vP5zPj4ePPjjz+O2fbnP/+5edlll+2z7ptvvtkcO3Zs9Pv48ePNUaNGxazz3HPPmT169DAjkUh0md/vN+Pi4swFCxaYpmma2dnZ5p/+9KdoeTAYNDt27NisLhFpsufvWiQSMd9//33T6XSakydPbvY7vzc9e/Y0H3vsMdM0TXPdunUmYL7//vt7Xff+++83zzvvvJhlW7ZsMQFz3bp1pmma5tChQ838/PyY3/O77rrLzM/PN03TNIuKikzAXLp0abS8vLzcjIuLM1966aUDOnZpvTQJQ9qsAQMGRD9v2LCBxsZGzj333Jh1AoEAffv2jX7/61//ytNPP83mzZvxer0EAgH69Onzg/tZuXIlGzZsIDExMWa5z+dj48aN1NTUsH37dgYNGhQts9lsDBgwQENoIj/grbfewu12EwwGiUQiXH755UydOpWXX345Zr36+nqmTp3K/Pnz2b59O6FQCK/XG+0pWrFiBVarlaFDh+51PytXrmThwoUxcxG/s3HjRvLy8gA49dRTMQwjWnbaaacxY8YMwuEwhYWF2Gy2mN/ztLQ0evToQWFh4SGfC2kdFIqkzUpISIh+rq+vB2D+/Pl06NAhZr3vXgY5d+5cJk+ezIwZMzjttNNITEzkoYce4rPPPvvB/dTX19O/f3/mzJnTrCwjI+NQD0PkuDVs2DCefPJJHA4H7du33+fNEpMnT+b999/n4Ycfplu3bsTFxXHJJZdEh9ri4uJ+cD/19fVcdNFFTJ8+vVmZJnXLnhSK5JhQUFCA0+lk8+bN+7xaXLp0KYMHD2bixInRZRs3boxZx+FwEA6HY5b169ePF198kczMTJKSkvZad3Z2Np999hlDhgwBmiZsLl++nH79+h3KYYkc0xISEujWrduPrrd06VImTJjA6NGjgaaQs+eNDr179yYSibB48WKGDx/ebPt+/frxz3/+k86dO//gXarfv0D69NNP6d69O1arlfz8fEKhEJ999hmDBw8GoKKignXr1lFQULA/hyttgCZayzEhMTGRyZMnc/vttzN79mw2btzIf/7zHx577DFmz54NQPfu3fniiy9YsGABRUVFTJkyhc8//zymns6dO7Nq1SrWrVtHeXk5wWCQK664gvT0dEaNGsWSJUvYtGkTixYt4pe//CVbt24FYNKkSTz44IPMmzePr7/+mokTJ1JdXX20T4PIMal79+68+uqrrFixgpUrV3L55ZcTiUSi5Z07d2b8+PFce+21zJs3L/o7+t1NFDfffDOVlZVcdtllfP7552zcuJEFCxZwzTXXxFwEbd68mV/96lesW7eOF154gccee4xJkyZF2zBq1Ciuv/56PvroI1auXMmVV15Jhw4dGDVq1NE9IXLEKBTJMeP+++9nypQpTJs2jfz8fEaMGMH8+fPp0qULADfeeCNjxozh0ksvZdCgQVRUVMT0GgFcf/319OjRgwEDBpCRkcHSpUuJj4/nww8/JDc3lzFjxpCfn8/Pf/5zfD5ftOfojjvu4KqrrmL8+PHRobnvrmpF5NA88sgjpKSkMHjwYC666CLOP//8Zr2wTz75JJdccgkTJ07kxBNP5Prrr6ehoQGA9u3bs3TpUsLhMOeddx69e/fmtttuIzk5GYtl95/Bq6++Gq/Xy8CBA7n55puZNGkSN9xwQ7R81qxZ9O/fn5EjR3LaaadhmiZvv/02drv96JwIOeIMUzNBRURERNRTJCIiIgIKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIAP8PqIw11zW7JC4AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "mosaic(ra, ['Treatment','Improved'], title = \"Arthritis Treatment\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "ca90f04b-49ab-4ef8-bf96-44f926f7640f", + "metadata": {}, + "source": [ + "### (Bonus) question\n", + "\n", + "What fundamental error does the default plot from `pandas` make?\n", + "\n", + "### Answer\n", + "\n", + "1. It uses red and green as the only colours which can be difficult for people\n", + " with colour blindness.\n", + "2. It has not respected the implicit ordering of the response values.\n", + "\n", + "The figure below makes the pattern in the data far clearer (at the expense of a\n", + "few extra lines of code)." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "6f429058-2786-401f-9116-f4657113ec77", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAGzCAYAAAAhXWNYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOx0lEQVR4nO3deVxUVePH8c8My7AvAoKKiIoauAtqae6amlvumZWWpaU9pWWZ/dIsn8elMkuztHpUNM0tM0vTrNRcKpfcJVQUN3AHFFln5v7+4HGM0MoVoe/79ZrXa+ace8899xbw9Zwz95oMwzAQERER+YczF3YHRERERO4ECkUiIiIiKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQif1PTpk2pVq3aDbWxZs0aTCYTa9as+cttExMTMZlMzJw584aOKSLydykUiRQzH3zwASaTifr161/zvklJSYwaNYrt27ff/I5dxdy5c3n33XdvSduXgtXfeSUmJt6SPvzRmDFjWLJkyW051o3auHEjo0aNIjU1tbC7InJbmPTsM5HipWHDhiQlJZGYmMj+/fuJiIj42/tu2bKFunXrMmPGDPr27ZuvrmnTppw5c4bdu3dfd9/sdjs5OTm4urpiNuf9m6x9+/bs3r27QCgxDIPs7GxcXFxwcnK6ruNdvHiRL774Il/ZhAkTOHbsGBMnTsxX3rlzZzw9Pa/rONfCy8uLbt26FYkRsLfffpsXX3yRQ4cOER4eXtjdEbnlnAu7AyJy8xw6dIiNGzeyePFiBgwYwJw5c3jttdf+cj+r1Yrdbr9l/crKynIEITc3t7+1j8lk+tvbXo2npycPP/xwvrJ58+aRkpJSoPz3DMMgKysLd3f3Gzq+iBQtmj4TKUbmzJmDv78/7dq1o1u3bsyZM6fANpemlN5++23effddKlasiMVi4YMPPqBu3boAPPbYY45ppT+OaOzdu5dmzZrh4eFBmTJlePPNN/PVX1o3NG/ePF599VXKlCmDh4cH58+fL7CmqGnTpixbtozDhw87jndpROJKa4pOnDjBY489RmhoKBaLhVKlStGpU6cbnvoKDw+nffv2rFy5kpiYGNzd3Zk2bRoAqampDB48mLJly2KxWIiIiGD8+PEFQuTbb79NgwYNCAgIwN3dnejoaBYtWpRvG5PJxMWLF4mNjXWc76URuVGjRmEymdi3bx8PP/wwvr6+BAUFMWLECAzD4OjRo3Tq1AkfHx9CQkKYMGFCgfPIzs7mtddeIyIiAovFQtmyZXnppZfIzs4u0I9nnnmGJUuWUK1aNSwWC1WrVmXFihWObUaNGsWLL74IQPny5W/7NKNIYdBIkUgxMmfOHLp06YKrqyu9evXiww8/ZPPmzY6w83szZswgKyuL/v37Y7FY6Ny5MxcuXGDkyJH079+fRo0aAdCgQQPHPikpKbRp04YuXbrQo0cPFi1axLBhw6hevTpt27bN1/7o0aNxdXVl6NChZGdn4+rqWqAP//d//0daWlq+6SwvL6+rnl/Xrl3Zs2cP//rXvwgPD+fUqVOsWrWKI0eO3PD0Tnx8PL169WLAgAE8+eSTVKlShYyMDJo0acLx48cZMGAAYWFhbNy4keHDh5OcnJxvLdR7771Hx44d6d27Nzk5OcybN4/u3bvz9ddf065dOwBmz57NE088Qb169ejfvz8AFStWzNePnj17EhkZybhx41i2bBn//ve/KVGiBNOmTaN58+aMHz+eOXPmMHToUOrWrUvjxo2BvKnJjh07sn79evr3709kZCS7du1i4sSJ7Nu3r8A6pvXr17N48WIGDhyIt7c3kyZNomvXrhw5coSAgAC6dOnCvn37+Oyzz5g4cSKBgYEABAUF3dB1FrmjGSJSLGzZssUAjFWrVhmGYRh2u90IDQ01nnvuuXzbHTp0yAAMHx8f49SpU/nqNm/ebADGjBkzCrTfpEkTAzBmzZrlKMvOzjZCQkKMrl27OspWr15tAEaFChWMjIyMfG1cqlu9erWjrF27dka5cuUKHO9SPy/1JSUlxQCMt956629cjau70vHKlStnAMaKFSvylY8ePdrw9PQ09u3bl6/85ZdfNpycnIwjR444yv54rjk5OUa1atWM5s2b5yv39PQ0+vTpU6Bfr732mgEY/fv3d5RZrVYjNDTUMJlMxrhx4xzlKSkphru7e752Zs+ebZjNZmPdunX52p06daoBGBs2bHCUAYarq6tx4MABR9mOHTsMwJg8ebKj7K233jIA49ChQwX6K1IcafpMpJiYM2cOwcHBNGvWDMibIunZsyfz5s3DZrMV2L5r167X/K9+Ly+vfGtxXF1dqVevHgcPHiywbZ8+fW7qmhx3d3dcXV1Zs2YNKSkpN63dS8qXL0/r1q3zlS1cuJBGjRrh7+/PmTNnHK+WLVtis9n48ccf8/XvkpSUFNLS0mjUqBG//vrrNfXjiSeecLx3cnIiJiYGwzDo16+fo9zPz48qVarku+4LFy4kMjKSu+66K19fmzdvDsDq1avzHadly5b5Rqlq1KiBj4/PFf9bivxTaPpMpBiw2WzMmzePZs2acejQIUd5/fr1mTBhAt9//z333Xdfvn3Kly9/zccJDQ3FZDLlK/P392fnzp0Ftr2e9v+MxWJh/PjxvPDCCwQHB3P33XfTvn17Hn30UUJCQm64/Sv1d//+/ezcufOq4fHUqVOO919//TX//ve/2b59e741PH+8Xn8lLCws32dfX1/c3Nwc01e/Lz979my+vsbFxf2tvl7pOJD33/JWBE6RokKhSKQY+OGHH0hOTmbevHnMmzevQP2cOXMKhKLrGcW52lfjjSvc2eNWfHNr8ODBdOjQgSVLlrBy5UpGjBjB2LFj+eGHH6hdu/YNtX2l/trtdlq1asVLL710xX0qV64MwLp16+jYsSONGzfmgw8+oFSpUri4uDBjxgzmzp17Tf240jX+O9fdbrdTvXp13nnnnStuW7Zs2WtuU+SfRqFIpBiYM2cOJUuWZMqUKQXqFi9ezBdffMHUqVP/Mqhc66jGzXCtx6xYsSIvvPACL7zwAvv376dWrVpMmDCBTz/99Kb3rWLFiqSnp9OyZcs/3e7zzz/Hzc2NlStXYrFYHOUzZswosO2tusYVK1Zkx44dtGjR4qYdozD+fxApTFpTJFLEZWZmsnjxYtq3b0+3bt0KvJ555hkuXLjA0qVL/7KtSzcvvJ13MPb09CQtLe0vt8vIyCArKytfWcWKFfH29i7wlfObpUePHvz000+sXLmyQF1qaipWqxXIG3UxmUz51m4lJiZe8c7Vnp6et+T69ujRg+PHj/Pxxx8XqMvMzOTixYvX3GZh/P8gUpg0UiRSxC1dupQLFy7QsWPHK9bffffdBAUFMWfOHHr27PmnbVWsWBE/Pz+mTp2Kt7c3np6e1K9f/6avD/q96Oho5s+fz/PPP0/dunXx8vKiQ4cOBbbbt28fLVq0oEePHkRFReHs7MwXX3zByZMnefDBB29J31588UWWLl1K+/bt6du3L9HR0Vy8eJFdu3axaNEiEhMTCQwMpF27drzzzju0adOGhx56iFOnTjFlyhQiIiIKrLeKjo7mu+++45133qF06dKUL1/+uh7J8kePPPIICxYs4KmnnmL16tU0bNgQm83Gb7/9xoIFCxz3YLoW0dHRQN6tEx588EFcXFzo0KHDbbnzt0hhUCgSKeLmzJmDm5sbrVq1umK92WymXbt2zJkzJ9/C3CtxcXEhNjaW4cOH89RTT2G1WpkxY8YtDUUDBw5k+/btzJgxg4kTJ1KuXLkrhqKyZcvSq1cvvv/+e2bPno2zszN33XUXCxYsoGvXrrekbx4eHqxdu5YxY8awcOFCZs2ahY+PD5UrV+b111/H19cXgObNm/Pf//6XcePGMXjwYMqXL8/48eNJTEwsEIreeecd+vfvz6uvvkpmZiZ9+vS5KaHIbDazZMkSJk6cyKxZs/jiiy/w8PCgQoUKPPfcc471T9eibt26jB49mqlTp7JixQrsdjuHDh1SKJJiS88+ExEREUFrikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQPcpuqqcnBy+/fZbwsPDr/qMIBEREbmz2Gw2EhMTue+++3B1db2mfRWKruLbb7+94g3kRERE5M731Vdf0b59+2vaR6HoKsLDwwGYOm0OZcPCC7UvInL7+NiOsWvKQ4XdDRG5TskXbIxef/nv+LVQKLqKS1NmZcPCqVjx2m+PLyJFk7/NhVRfTZmLFG2261r6ooXWIiIiIigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUSTHVvFkdtm75+Za1f1flIE6cSLpl7YuIyO3nXNgdELkRzZvV4dzZM5jNJjw9vWjTthPDXn6jsLslIiJFkEKRFHn/nb6A6Ji7SUxM4JHeHalQoVJhd0lERIogTZ9JsREeXpHomLvZv/+3fOU7dmylW5dWRNcuT/OmtZk96+N89bNmfUTrVvWoUzuc7l3vIyXlHADx8Xt5+KEO1IuJoEvnFuzatT3ffqu+XUazJrW4t0EU//3kfUd5dnYWb4x6iXsbRNGsSS2mvP82drv91py0iIjcNApFUmwcPHiArVt+JjKyWr5yZ2dnXh/9Npu3JvDe5Bm89+5Y9u7dCcBXSxcxO/ZjJk+JZcvWg4x64y1cXFy5eDGdJ/v15JE+/fnpl3gGDnyBZ5/pS3Z2lqPdNWu+5atl65g1+0tmzviQnzb+CMAHUyaw/0A8y77ZyJzPvmbp0oUsWTL/9l0IERG5LgpFUuQ9+eSD1I2uSP8nHuSBzg/StVvvfPVVq9akatWamM1mqlevReMmLfl16yYAvlg8jyf7P0vlypGYzWaqVq2Jl5cXa1avIqLSXbRu3QEnJydatrqfEgGBbN++1dFu/wGD8fLypkLFSnTt1pvly5cAsHzZEgY98yK+vn6ULh3KY48PZNnXi2/b9RARkeujNUVS5H388TyiY+6+av3+/b8x9j+vEhe3m9zcHLKzsx3rjk6cSCI0NKzAPsnJx9i8aSN1oys6yqxWK6dOnXB8LlWqtON9SKky/PbbHgBOnTpB6dJlHHVlSody6tTJ6z9BERG5LRSKpNgb/cbLxMTczQdTZ+Pm5s7zQ/pjGAaQF2aOHz9SYJ/g4FLc26gZH0799KrtJicnERZWHoATyccJCgoGoGTJEJKSjjvqkpKOU7Jk8M0+LRERuck0fSbF3sWL6Xh7+2KxuLFl80+sXbPKUde5y4N88vFkDhyIxzAM9uzZQXp6Ok2b3Ufc3l18t2o5VquVrKxM1v34PRcunHfs+8lHk0hPv8DBgwf4/PO5tG3bCYC293fiww8mkJaWSnLycWbO+JD723W+7ectIiLX5o4LRTNnzsTPz++mtrlmzRpMJhOpqak3tV0pGoYOHcncOdOJrlOe2NhpNG/RxlHXoUNXHur9OE/17010nfK8MWoYVmsu3t4+TPtoLp/O/oSG90TSolk0C+bPztdu4yYt6dCuEY/07sgjjzxJg4ZNABg46AXKl4/g/jYNeLBHW+5v15nOnR+8recsIiLXzmRcmkf4G/r27UtsbCwDBgxg6tSp+eoGDRrEBx98QJ8+fZg5c+Z1d2jmzJkMHjz4pgaYNWvW0KxZM1JSUv524IqLiyMqKopl32ygYsXKN60vInJn87cdYtPr9xZ2N0TkOh1NszHwGxt79+4lMjLymva95pGismXLMm/ePDIzMx1lWVlZzJ07l7CwggtWr0Vubu4N7S8iIiJyva45FNWpU4eyZcuyePHlrxgvXryYsLAwateu7ShbsWIF9957L35+fgQEBNC+fXsSEhIc9YmJiZhMJubPn0+TJk1wc3Njzpw5BY53+vRpYmJi6Ny5M9nZ2djtdsaOHUv58uVxd3enZs2aLFq0KN8+y5cvp3Llyri7u9OsWTMSExOv9TRFRETkH+a61hQ9/vjjzJgxw/F5+vTpPPbYY/m2uXjxIs8//zxbtmzh+++/x2w207lz5wJ39n355Zd57rnniIuLo3Xr1vnqjh49SqNGjahWrRqLFi3CYrEwduxYZs2axdSpU9mzZw9Dhgzh4YcfZu3atY59unTpQocOHdi+fTtPPPEEL7/88l+eU3Z2NufPn3e80tPTr+fSiIiISBF1XV/Jf/jhhxk+fDiHDx8GYMOGDcybN481a9Y4tunatWu+faZPn05QUBB79+6lWrXLdxwePHgwXbp0KXCM+Ph4WrVqRefOnXn33XcxmUxkZ2czZswYvvvuO+655x4AKlSowPr165k2bRpNmjThww8/pGLFikyYMAGAKlWqsGvXLsaPH/+n5zR27Fhef/3167kcIiIiUgxc10hRUFAQ7dq1Y+bMmcyYMYN27doRGBiYb5v9+/fTq1cvKlSogI+PD+Hh4QAcOZL/njAxMTEF2s/MzKRRo0Z06dKF9957D5PJBMCBAwfIyMigVatWeHl5OV6zZs1yTM3FxcVRv379fO1dClB/Zvjw4aSlpTlemzZt+tvXQ4qXY8eOUDUy5IbbWbz4Mx7r0/WvNxQRkTvCdd+88fHHH+eZZ54BYMqUKQXqO3ToQLly5fj4448pXbo0drudatWqkZOTk287T0/PAvtaLBZatmzJ119/zYsvvkiZMnl3B740pbVs2TJH2e/3uREWiyVfG15eXjfUntw+zZvVIeXcWTb8tBcPj7z/nzIzM2jYIAo/vxL8sPrXQu6hiIgUBdd9n6I2bdqQk5NDbm5ugbVAZ8+eJT4+nldffZUWLVoQGRlJSkrK3++U2czs2bOJjo6mWbNmJCUlARAVFYXFYuHIkSNERETke5UtWxaAyMjIAqM8P//88/WephQRJYND+P67bxyfv/9+heMO09fCarXezG6JiEgRct2hyMnJibi4OPbu3YuTk1O+On9/fwICAvjoo484cOAAP/zwA88///w1tz9nzhxq1qxJ8+bNOXHiBN7e3gwdOpQhQ4YQGxtLQkICv/76K5MnTyY2NhaAp556iv379/Piiy8SHx/P3Llzb+i+SVI0tGvXha+WXv4W4ldfLqRDh26Oz1M/nEjzprWJrl2enj3aEv+/55RB3kjTJx9P5v42DbivVb0CbX86+xM6dWjK2bOnycrK5I3Xh9Ho3mo0aVSDj6a959guI+MiQ194irrRFenyQHMOJx68RWcrIiK3wg3d0drHxwcfH5+CjZrNzJs3j61bt1KtWjWGDBnCW2+9dc3tOzs789lnn1G1alWaN2/OqVOnGD16NCNGjGDs2LFERkbSpk0bli1bRvnyec+ZCgsL4/PPP2fJkiXUrFmTqVOnMmbMmBs5TSkC7r7nXvbti+PcuTOcO3eG+Pi9NGjQ2FFfoUIlFi1exc+b9tGwYROGvTQo3/6rvl3GrE+XsPybDfnK//vfKSz+fC4zZy0mICCIN8e9RlpaKitW/syCRd+y9MuFrF79LQBT3n+bs2dOs3rtdt6eMI0vlyy49ScuIiI3zTXd0fqfRHe0LjqaN6vDW299wLfffk3ZsuEAHDlyiPvbdeb5If0LrCnKzs6iVo0wtvx6EE9PL5o3q8PzL7xK+/Z534I8duwIrVvV45l/vcgP36/kk+kL8PX1wzAMatcsx8pVmwgOzluI/ensT9i9ezvjxr9Pi+bRjBk7ifr1GwLw7sQx7Ni+lRmxn9++iyE3THe0FinabuSO1te90FrkTtO+Qzf+8+9XMAyDV18dg81uc9QtmD+bWbHTOHEiCZPJhGEYpKam4OmZt6A+JLhUvrbsdjuxMz9ixMhx+Pr6AXDu3BmysjJpf3/DfNvVqZM35Xb61ElKlSrtqAspVYYd27feqtMVEZGbTKFIio3q1WuRlpq3oL96jdps374FyBv5GTvmVWZ9+iVVq9YgJyeb2jXL8ftB0ku3fbjEbDbz8X/n8/SAhygZHEzdug3w9w/A1dXCqu+34OfnX+D4QSWDSU5OIiwsbyr3RPLxW3WqIiJyC9zQmiKRO83kKTOZPGVmvrKMjIuYzCZKlAjAarUyedKbf6ut6tVrMfG9Txjy3BPs2rUds9nMA517Mn7cSM6fT8Nut5NwYB87d+RNz7Vu3YFpU98lPf0CBw8e0JoiEZEiRqFIipWIiCpERFTJV1a5ciQ9e/ahU4cmtGheh9DQMFxcXP9We3XrNuA/Y95j4NMPc+BAPMNfGY23tw8d2zemft1KDHtpEGnnUwEY9MxQ/Pz8adq4JkOf70/HTt1v9umJiMgtpIXWV6GF1iL/TFpoLVK03chCa40UiYiIiKBQJCIiIgIoFImIiIgACkUiIiIigEKRiIiICKBQJCIiIgIoFImIiIgACkUiIiIigEKRiIiICKBQJCIiIgIoFImIiIgACkUiIiIigEKRiIiICKBQJCIiIgIoFImIiIgACkUiIiIigEKRiIiICKBQJCLiYDKBk5OpsLshIoXEubA7cKfzdnfB19O1sLshIreBm7szplSnwu6GiBQShaK/8MM38cT9ml3Y3RCR28DJycT9zRWKRP6pNH0mIvI/NptBTpatsLshIoVEoUhEREQEhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIiqlXRnfnwMGdt6z9Ac83IiX11C1rX0REbj/nwu6AyI14ZXR3zqenYDaZsVjcianVnG4dBxV2t0REpAhSKJIib/CAd4ioUIOTp4/y9vvPEFIyrLC7JCIiRZBCkRQbwUFlqVShJkknEvOVHzq8h88Wv8vJU0fw8PCmVdMHad6om6P+hx8XsXr956RdOEupkuH8q/9beHn6cjwpgbmLJ5KUnEBgQGl6d3+R8LJ3OfbbtutHvl39GTablZZNetK6+UMA5OZms3DpFLbtXIuTkzP31m/P/a36YDZrtlqkOOq3NIfn73GmapB+xos6hSIpNk6cOsL+gzvo1PYJdu3d6Ch3cnLm4e5DCS1diSPH9zHxw8FElK9BWGhlftn6Ld+vW8jAx8dSKjicY0n7cXZyISs7g0kfDaVn5+eoVb0RO/dsYOqM/2P08Lm4uFgA2LV3I6+9NIu0tDNM+OBZwkIrE1k5hmWrYkk+cYhRwz4lOzuDiVMHU8I/mAb17i+sSyNSZHVflON4n2UFixOYTHmfp7R1oaSn6YaPMfz7XO6raKZZuNMNtyVFm0KRFHmTPh6K2eSEp4c3Deq2pUG9dixfNctRHxZaxfE+vOxdVI+8m4RDuwgLrcxPm7+hTfPelClVId+2m7d9R6mQ8tSp2RSAWtUbs2zVLA4e3kuViNoAtGnxCO5unri7edKwfnu2bP+ByMoxbN72PY/2GIanhzeeHt60avIgm7d9p1Akch0WdnN1vO+yIIcpbV0I9sofhAzDwADMphsPSPLPplAkRd6zT75NRIUaV61POnGIBUsmcfT4fqw2K9bcHIJLlgMgJfU0gQGlC+xzLuUk+xO2M/iVto4ym91KWtoZx+cSfiXzvT+WdACAtLQzlPAPvlznH0Jq2tnrP0ERKWDiz1bcXeBomsFvZw0mt3Ehx2Zn6lYbiakGIV4mBtV1olKJvCmtBXtsrEiwkZ4DYb55deX9zMzbbWPvGYP4szY+2GKje6QTPao6seuUnf9us3Ei3aC8n4ln6zlTyjsvdG1JsvPRr1Yu5EDHyhpdKk4UiqTY++zziVSqUIOBj4/D1dXCJ7NHgWEA4O9XkrPnkgvs4+cbRNRd9RjUb9xV2z2XeoqgwDKO974+AQD4+gZyLuXk7+pO4ucbcJPPSkR+PGznjabOlPczkW2Dp5db6V/HmbvLmNiUZDBmvZVp7VxwdTIR6mNi4n0ueLrCvD02Jv5sY1IbMw9Wc2LHSXu+6bPTFw3GbbDySkNn7go0sfyAnTd/sjLxPhfSsgze3Ghl6D3O1A4x8ekuG2cyC/lCyE2jVWFS7GVlZ+Du7oWLiyv7D+5g196fHHUN6rZl5Q9zSTpxCMMwOHIsnqysDGpENeDo8f1s3/UjNpuVnJxsdsf9QmZmumPflT/MITPrIidOHWHjpmXE1GoOQEyt5ixbFcvFjAucSznJd2vnE1O7xW0/b5HirkFZMxElzDiZTWxJshPmY6Jh2bzP94Sa8bOYiD9jOLb1dTPhbDbRI8qJxFSDzFzjiu2uOWznnjJmqpbMa6tDZSdOXTQ4mW6wJdlORX8T9cqYcXEy0auaE2bN2hUbGimSYq9L+6f4dOFbfLVyOlGV61GjakNHXb3oVpxPT2HKJ8O4cDGV0sHl+deTb+Hp6cMzT4xnwZeTiZ03DicnZyqWr07F8KqOfatF3s3rbz6K1ZpLyyY9iKwcA0C7Vn1ZuHQyr43r7fj22T0xbQv0S0RuTKD75TRyOgN2nzZ48PPLC7OtdjibmRd8VibY+DLeztlMAxNgABdywN2lYLunMwxWH7az/ujltnL/19a5TAj0uHxcN2cT3q4F25Ci6ZaHotOnTzNy5EiWLVvGyZMn8ff3p2bNmowcOZKGDRv+dQMif2LMiIV/WR5ZOYb//N/8q7bRskkPWjbpUaA8tHQEzz/93hX3mfbOOoB8X+2/xNXVQu9uQ+ndbeif9l1Ebszv11UHuEOdEBMjGhdMOSfTDT7+1cbYFs5U9DeRa4Pui3IvzaLzx/XZAe4m7qtgZkB0wT+Ryek2fk2+PMKUbTW4kFNgMymibvn0WdeuXdm2bRuxsbHs27ePpUuX0rRpU86e1cJTERG5OeqWNnMwxeCnY3ZsdoNsq8HWZDsXcwyyrAZmE/haTNjsMHe3Ld++vhYTpy5e/tyknJkNR+3sOWXHbhhk5BpsOGoHIKaUmYQUgy1JdnJtBp/tsWG/8iycFEG3NBSlpqaybt06xo8fT7NmzShXrhz16tVj+PDhdOzYEYAjR47QqVMnvLy88PHxoUePHpw8edLRxqhRo6hVqxbTp08nLCwMLy8vBg4ciM1m48033yQkJISSJUvyn//8p8Cxn3jiCYKCgvDx8aF58+bs2LHjVp6uiIgUEk9XEyObOPP1PhsPL8nlia9yWZmQF37K+ZlpXdHMv1bklQd7mnD+3V+/9pXNfH/IxoOf57Bwr40QLxMvNnBmxg4bDy3O5enlufx8LC8U+bqZGHqPM9O2WnlkSS4WJxOB7oVxxnIr3NLpMy8vL7y8vFiyZAl33303FoslX73dbncEorVr12K1Whk0aBA9e/ZkzZo1ju0SEhL45ptvWLFiBQkJCXTr1o2DBw9SuXJl1q5dy8aNG3n88cdp2bIl9evXB6B79+64u7vzzTff4Ovry7Rp02jRogX79u2jRIkSBfqanZ1Ndna243N6enqBbUREpPAs7nF58c6Quwv++SrvZ+Y/za/8b/1+tZ3pV/vy5/srXf4qfdUgMx+1z78wqHpJM2+3unJb9cqYqVfm8va9qulr+cXFLR0pcnZ2ZubMmcTGxuLn50fDhg155ZVX2Lkz7+nl33//Pbt27WLu3LlER0dTv359Zs2axdq1a9m8ebOjHbvdzvTp04mKiqJDhw40a9aM+Ph43n33XapUqcJjjz1GlSpVWL16NQDr169n06ZNLFy4kJiYGCpVqsTbb7+Nn58fixYtumJfx44di6+vr+NVr169W3lpRERE5A5zW9YUJSUlsXTpUtq0acOaNWuoU6cOM2fOJC4ujrJly1K2bFnH9lFRUfj5+REXF+coCw8Px9vb2/E5ODiYqKiofM+SCg4O5tSpUwDs2LGD9PR0AgICHKNVXl5eHDp0iISEhCv2c/jw4aSlpTlemzZtutmXQkRERO5gt+Ur+W5ubrRq1YpWrVoxYsQInnjiCV577TVeeOGFv7W/i0v+bxOYTKYrltnteXO+6enplCpVKt8U3CV+fn5XPIbFYsk3vefl5fW3+iZ3rn0J21n81YcknzqMk9mJ0NIRPNpz2BXvYC0iIlIo9ymKiopiyZIlREZGcvToUY4ePeoYLdq7dy+pqalERUVdd/t16tThxIkTODs7Ex4efpN6LUVJZmY6H04fziM9X6ZWtUbk5GYRF78Fs1lz/yIicmW3NBSdPXuW7t278/jjj1OjRg28vb3ZsmULb775Jp06daJly5ZUr16d3r178+6772K1Whk4cCBNmjQhJibmuo/bsmVL7rnnHh544AHefPNNKleuTFJSEsuWLaNz58431LYUDSdPH8XZ2ZU6NZoA4GbxoHaNxgDk5mazcOkUtu1c67i54v2t+mA2m/lqxXROnj6K1ZrL3vhNhJauyIC+/+brlTPYtO07ggJKM6Dvvwn632jT8aQE5i6eSFJyAoEBpend/UXCy95VaOctIn/tZLpB/2W5fNnzxu66+N1BG2sO2/l3syvcAVKKpFu6psjLy4v69eszceJEGjduTLVq1RgxYgRPPvkk77//PiaTiS+//BJ/f38aN25My5YtqVChAvPnX/1Ge3+HyWRi+fLlNG7cmMcee4zKlSvz4IMPcvjwYYKDg/+6ASnygoPKYrXmMmv+OPbGbyYz6/JNSJatiiX5xCFGDfuUl/71Ab/8+i0/b1nhqN+xZz0tmnTnnX8vw9XVnfGTnqZKpTq8M/prwspU5uuVM4C8x4dM+mgoLRp1Y8Lor2nXqg9TZ/wfubnZBfojItev39Icui3MIct6+YZAWVaDHoty6LdUd06Um+eWjhRZLBbGjh3L2LFjr7pNWFgYX3755VXrR40axahRo/KVzZw5s8B2f1w/5O3tzaRJk5g0adK1dFmKCXd3L4Y+8z4rf5jDjLn/JiPjAnVqNqV3t6Fs3vY9j/YYhqeHN54e3rRq8iCbt31Hg3r3A3BXpWgqVagJQO3qjVi7cYnjuWZ1ajbli2XTANi1dyOlQspTp2ZTAGpVb8yyVbM4eHgvVSJqF+yUiFy3AA/4+Zidpv97aOsvx+34u4PV9hc7/oFNd1qUP6Fnn0mxVaZUBR7vPQKAw0fj+Sh2BMu/m0Va2hlK+F8eMSzhH0Jq2uU7rHt7+Tneu7hY8Pbyz/c5OzvvkdjnUk6yP2E7g1+5/Fwzm91KWtqZW3VKIv9YjcPMrDl8ORStSbTTtJwT3x3MS0UL9thYkWAjPQfCfE0MqutEeb+8yZB+S3O4v5IT3x2ykWOFMc3zT3d9vc/GyoN2Rjd1xt0ZZuyw8dNROyYTtKvkRPeovGNmWQ3e32xjS5KdYC8T0aX0JNjiRqFI/hHKla1C7RpNSDpxEF/fQM6lnCQosAwA51JP4ucbcM1t+vkGEXVXPQb1G3ezuysif1Aj2Mx3B62kZeWN9CSmGnSPMvHdwbz6UB8TE+9zwdMV5u2xMfFnG5PaXF4h8tMxO2Oau+DhDKlZl9td/JuNtYft/KeZMz4WEx9usZKeA1PbuZBhhRGrrZTzNVGvjJnPdttIzTKY3tGFs5kwcnUuZXwUjIqTW36fIpHCcOLkYb5bO5/U/43anDx1hJ17NhBeNpKYWs1ZtiqWixkXOJdyku/WziemdotrPkaNqAYcPb6f7bt+xGazkpOTze64X8jM1N3QRW42swkalDWz7oiddUfsNChrxvy7PNKgrBlfNxPOZhM9opxITDXIzL08Vdahshl/NxMW58s7zd9jY/2Ry4HIMAy+O2Tn8dpOuLuYCHA3cX+E2fHcs/VH7fSMcsLDxURZHxPNy+tPaHGjkSIpliwWDxISd7Pyh7lkZWfg6eFDnZpNadPiYex2OwuXTua1cb0d3z67J6btXzf6B+7uXjzzxHgWfDmZ2HnjcHJypmL56lQMr3oLzkhEmoab+WirDQPoX8cp34NYVybY+DLeztlMAxNgABdywP1/M2WB7vlHdAwDvoy38VS0E16ueXVp2ZBjg0HLcx3b2Q2IDMqrT8mEQI/L7QR6mIg/qzVKxYlCkRRL/n5BDOgz+qr1vbsNpXe3oQXKO7R5PN/nBvXudyzABqgSUZt//988x+fQ0hE8//R7N6HHIvJXKpUwcyHHCkDlADO/nckbwTmZbvDxrzbGtnCmor+JXBt0X5SL8fu88odZLpMJXm/qzBs/WinhbqJaSTM+FnAxw8ftXfC2FJwW83eHMxkGpbzz6s5kKBAVNxr7ExGRIuOVe1145d78C6WzrAZmE/haTNjsMHf33/tKWqUSZoY1cGb8Riv7z9kxm/KmxP673UZ6joHdMDiaZrDvbF74aljWzIK9NjJyDY6dN/gh0X7Tz08Kl0KRiIgUGWG+JsJ884/ilPMz07qimX+tyOWJr3IJ9jTh/Df/ulUraea5es78+0crR9IMnqjthKcL/OubXB5anMs7v+QtvAboVdUJH4uJx5bm8tZPVpqF609ocWMyDEPjf1cQFxdHVFQUo4bNplRweGF3R0Ruk9aNTOya0LywuyEi1+lomo2B39jYu3cvkZGR17SvYq6IiIgICkUiIiIigEKRiIiICKBQJCIiIgIoFImIiIgACkUiIiIigEKRiIiICKBQJCIiIgIoFImIiIgACkUiIiIigEKRiIiICKBQJCIiIgIoFImIiIgACkUiIiIigEKRiIiDk5MJVzenwu6GiBQS58LuwJ2ues0QKoSHFnY3ROQ28A/0wsP5aGF3Q0QKiULRX5j6ynJ8LcGF3Q0RuQ1cXJ0Y80n9wu6GiBQSTZ+JiPxPbo6Ni+ezCrsbIlJIFIpEREREUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEJB9njxKYnS2F3Q0RKQTOhd0BEZE7icXdQkyT2uRmZxR2V0TkOnglZ8I3e65rX4UiEZHfMVkv4OZmxs3Nq7C7IiLXwfPi9e+r6TMRERERFIpEREREAIUiEREREUChSERERARQKBIREREBFIpEREREAIUiEREREUChSERERARQKBIREREBFIpEREREAIUiEREREUChSERERARQKBIREREBFIpEREREAIUiEREREUChSERERARQKBIREbkh1Yf8xE/xqYXdDbkJnAu7AyIiIldT+okfHe8vZtvxcDVjMuV9/mVcPcoGut3wMdr9ZxuPNi1Fz4YhN9yWFG0KRSIicsdK+qSx433Jx9fy87i6lAtyz7eNYRgYBpjNptvdPSlmFIpERKTIeXpaHJ5uTvx2PIPNB86zcUwMWbl2hsbuZ8/Ri4SXdOPdxypTp4IPAG8vPczMH5JIvWjlrlBP3n2sMtXCvHhzSSIb49PYnHCeITP28XzHcgztWI71cam8MvcAiacyqR7mxeQn7qJCcF4YW7XjLC/O2k/KRStP3xdamJdBbjKtKRIRkSLp859PMbpXRY5/fC8lfV3p+tZOnmodyqEPG/JSp3I8MmkPWTk2ACqX8mDNG9Ec+rAhzar589S0OABeeiCcBlV8mdyvCkmfNGZox3IcO5vFo5P3MO7hCA59eC8d6wbx+JS9AJw5n0Pf9/cy7uEI9r/fgIvZNo6fyy60ayA3l0KRiIgUSR3rBlG7vDfOTmZWbj/LXWU86VQ3CCezifYxQQT6uLA54bxj20AfV1yczbzQIYzdRy+SnmW9YrsLNp6kfXQgDar44WQ2MeC+UI6cyeLw6Uy+3XGOWuFetKkdiKuzmeFdwjGbNG1XXGj6TEREiqQyJSyO98fOZrPht1TCBqxzlOVaDU6k5AAwc3USH648RtK5bEwmMAw4l27Fy63gn8FjZ7OZv+EkSzadcpTlWA2SU3I4kZpNmYDLi7s9LE6U8HK5FacnhUChSEREiqTfj8+U9rfQonoJ5j1fvcB2h09nMvzTAyz7v1rUCvcm22qn1BPrMAwjr50/jPSU8rfQp2kp3ny0UoG2Dp7M5Pud5xyfM3NsnEvPvTknJIVO02ciIlLkta4dwM7D6Xy95TRWm53MHBvf7TxLWoaVi1k2zGYTgT6uWO0GYz5PzLdvoI8LR05nOT53v6ckSzadZmN8Kna7wYVMq2PU6L6aJdiemM6328+SY7Uz7otE7P8LV1L0aaRIRESKPF8PZxa8UJ3hcw4w6JN4XJxM1K/sS90IX6LKetG3WSkavrIZD4sTL3Yqh6vz5dGhAa3K8PRHvzH5m6MMbh/G8x3KMX1QFCM+S2B/cgYeFicaR/nzQL2SBPq48t+Bkbw4ez8p6VYGtg7NN40nRZvJMBRxryQuLo6oqChahQ/B1xJc2N0Rkdvk7dmN8Y9/pbC7ISLXKT4pnXrDtrB3714iIyOvaV9Nn4mIiIigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigr+SLiEgRteG3VEbOS2BfUgZOTiaqlfXi/SeqEF7SvbC7JkWUQpGIiBQ5aRlWek3czftPVKF9dCAZOTZW707ByaznkMn1UygSEZEi58CJDCwuJjrWDQLAy82ZDjF577NybPzf3ASWbjmNq7OZR5uU4sVO5TCbTYxdfIgDJzLJzrXzw65zVAvzYvazVRn7RSKf/3SK8JLuzHq2KuX/N9q052g6Q2P3s+foRcJLuvHuY5WpU8Gn0M5bbi2tKRIRkSInIsSDHKvBM5/8xg+7znE+8/IT79/88jC/Hb/IpvH1WPFqbeZvPMln60846pdvPcPANqEkTr0XD4sTLV/fRuMofw592JAa5bwY/0UiAOlZVrq+tZOnWody6MOGvNSpHI9M2kNWju12n67cJgpFIiJS5Ph6OPPN/9Uix2owYFocFQdu4MkP93Ih08rnP59iWOdw/D1dKBvoxjNty7Lo58tPvG9S1Z8GVfxwdTbTPjoQb3cnutQvibOTmQfqBbH7SDoAK7ad5a4ynnSqG4ST2UT7mCACfVzYnHC+sE5bbjFNn4mISJEUVdaLj57Ke4zDtkMX6DN5D28vPcyJlBxCA9wc24UFWjiRkuP4HOjj4njv7mrO99nN1Ux6Vt5I0LGz2Wz4LZWwAesc9blWI19bUrwoFImISJFXu7w3HWICiTt2kRB/V46dzaJCcN66oKNnsgnxd73mNkv7W2hRvQTznq9+s7srdyhNn4mISJGzL+kiU745SnJKNgAHkjNYse0s0RV86FK/JG8uOUzKxVyOnc1iyoqjdL275DUfo3XtAHYeTufrLaex2uxk5tj4budZ0jKsf72zFEkaKRIRkSLHy82ZX/af571lR7mQaaWEtwud6gYxpEMYNrvBK3MSqPvSJlycTDzatBQP3Rtyzcfw9XBmwQvVGT7nAIM+icfFyUT9yr7UjfC9BWckdwKTYRhGYXfiThQXF0dUVBStwofgawku7O7INVqeMI56pR4k0CP8lrS/KP5l7q8wHA8X/XIsbt6e3Rj/+FcKuxsicp3ik9KpN2wLe/fuJTIy8pr21UiRFGnLE8aRbbsImHAxWwj1rk6Nku0Ku1siIlIEKRRJkdcotB+BHuFcyDnD2iPT8Ha99rUDIiIiCkVSbHi7BhLoEc75nJP5ys9mHmH7qaVcyDmNq9mdyiUaEeHf0FG/P2UDCSkbybJewNtSkntDH8Pi5Ela9gm2nVxCWvYJPF0CqBPSmRJuoY79ktJ3E3/uR+yGjcolGlGlRBMAbPZcdp5exrELuzGbnCjvW5fIgOaYTPpeg4jInUyhSIqNCzmnOZORSNXAViSnxznKzSYn6gR3xs9SipTsJH48+jEB7uH4u5XhyPltHEjZQIMyj+LjWpLU7GTMJmes9mzWHZtOrZIdKONVlaT0OH46Pps25YfiZM67p0ly+m/cFz6ETOt5fjz6MX6WMgR7RhB39gfSsk/SuvzzWO05/Hj0Ezxc/Aj3jSmsSyMiIn+D/ukqRd76Y9P5cv8o1h+bQTnf6ALhw9+tDP5uZTCZzJRwC6WU512czTwMQGLaVqqUaIqvJQSTyYy/WxlczBaS03/DxzWYUO/qmExmynhXxc3Ji7NZRxztVgloiouTGz6WkoT7xnDswk4Ajl7YQVRgS1ydPPBw8aNyiUYcPb/j9l0QERG5LhopkiLv3tDH//RbZmnZJ9lx6itSs5OxG1bshhVv17wHR2Za0/B08S+wT0ZuKmcyD/Ll/lGOMrthI8t6wfHZw9nv8nsXP9Kyk//X5vkCdZlWPRZAROROp1Akxd72k18S6BFOwzJ9cDK78EvSZxjk3YnC3dmPjNyUAvu4u/gQ7FGZhqF9rtpuhjUVL9eAvPe5qbg5+/yvTZ8Cde7Oeqq2iMidTtNnUuxZ7dm4mN0xm5w5nXGI5Iu/OerCfaOJP7eW89knMQyDlKzj5NqzKeUZSWp2Escv7MFu2LDZczlxMZ5cW5Zj3/hza8m1ZXEh5zSJaVsI9c57FECodw3izn5Pji2DjNxU9qesp6xPzdt+3iIicm00UiTFXrWgtvx6cjF7z6wi2LMSpT0v38wrzKcW2dZ0NhyPJduajo8lmHtDH8PVyYOGZfqy4/TXbDmxCLPJiQD3cgSElHPsW8qzCt8mTsRu2Kjkfy/BnpUAiApowY5TX7Py0DuYTWbCfetSzqfObT9vERG5Nrf1jtZ9+/YlNjaWsWPH8vLLLzvKlyxZQufOnbmTbq6tO1qL/DPpjtYiRduN3NH6tk+fubm5MX78eFJSCq7jEBERESkstz0UtWzZkpCQEMaOHXvVbT7//HOqVq2KxWIhPDycCRMm5KsPDw9nzJgxPP7443h7exMWFsZHH32Ub5ujR4/So0cP/Pz8KFGiBJ06dSIxMfGqx8zOzub8+fOOV3p6+g2dp4iIiBQttz0UOTk5MWbMGCZPnsyxY8cK1G/dupUePXrw4IMPsmvXLkaNGsWIESOYOXNmvu0mTJhATEwM27ZtY+DAgTz99NPEx8cDkJubS+vWrfH29mbdunVs2LABLy8v2rRpQ05OzhX7NXbsWHx9fR2vevXq3fRzFxERkTtXoXz7rHPnztSqVYvXXnutQN0777xDixYtGDFiBJUrV6Zv374888wzvPXWW/m2u//++xk4cCAREREMGzaMwMBAVq9eDcD8+fOx2+188sknVK9encjISGbMmMGRI0dYs2bNFfs0fPhw0tLSHK9Nmzbd9PMWERGRO1ehfSV//PjxxMbGEhcXl688Li6Ohg0b5itr2LAh+/fvx2azOcpq1KjheG8ymQgJCeHUqVMA7NixgwMHDuDt7Y2XlxdeXl6UKFGCrKwsEhISrtgfi8WCj4+P4+Xl5XWzTlVERESKgEL7Sn7jxo1p3bo1w4cPp2/fvte8v4uLS77PJpMJu90OQHp6OtHR0cyZM6fAfkFBQdfVXyl6lieMA6BNhRcxm5wA+PXEF1icvaga2KowuyYiInegQr1P0bhx46hVqxZVqlRxlEVGRrJhw4Z8223YsIHKlSvj5OT0t9qtU6cO8+fPp2TJkvj46E7C/2S59mwS07ZQwa9+YXdFRETucIUaiqpXr07v3r2ZNGmSo+yFF16gbt26jB49mp49e/LTTz/x/vvv88EHH/ztdnv37s1bb71Fp06deOONNwgNDeXw4cMsXryYl156idDQ0FtxOnIHquR/L7+dXU24b4xjtOiShJSf2Jeyjlx7NiEelagd/AAuTm6cykhg64nFlPOpzf6UDbiYLdQJ6UKIZ2UAcmwZbDv5JacyDuBsthAZ0IJw3+jCOD0REbmJCv0xH2+88YZj2gvyRnkWLFjAvHnzqFatGiNHjuSNN964pik2Dw8PfvzxR8LCwujSpQuRkZH069ePrKwsjRz9w5T0qIiHix+JaVvzlZ+8uJ+4sz/QsEwf7q8wDJuRy/ZTSx31GbkpmE3OdIh4lbsCmvHricWOuk3J83Fz9uH+CsNpWKYvu0+vIDUr+badk4iI3Bq3daToj1+rh7x7DmVnZ+cr69q1K127dr1qO1e639D27dvzfQ4JCSE2NvZ6uinFTGRAS7ae/DzfaM7RCzso71cPn//drbxaUBtWJb5HTEh3AJzMLlQp0RiTyUyYT21+PfkFObZM7IaV0xkHaVDmUcwmJ3wsJSnrU4vj6bvxcytVKOcnN5fh7A1mF7DnFnZXROQ207PPpNgL9ozA3dmHw78bLcq0nifA7fJzzDyc/bAbVnLsGQBYnLwwmfIGUp3NrgBY7TlkWc9jM6x8dWC0Y1/DsBPmU/t2nIrcBiaLD57RT2BYLxZ2V0TkOrj5JgFbrmtfhSL5R4gKaMmvJ78gyKMCAO7OPmRYUx31GdZUzCZnXM0ef9qOu7MPzmZXOka8hslkupVdlkLiZLJidvMGvAu7KyJyHZw8rn+Ut9DXFIncDsGelbA4eZF0YS8Aod41OJS6ifPZp7Dac9h9+lvKetf4y6Dj7uJLgFs59pz5Fqs9B7thIyXrOOezT96O0xARkVtII0XyjxEV2JL1x6YDEOJZmbsCmrL++Ays9myCPSpRs2T7v9VOvdI92XlqGd8cfBO7YcPXEvy39xURkTuXyTAMo7A7cSeKi4sjKiqKVuFD8P3fYlwRKf4mftGdMh7xhd0NEblOcQeOUrXlAPbu3UtkZOQ17avpMxEREREUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhERuSHlG/Zl/eY9hd0NuQmcC7sDIiIiV+Md1cXx/mJGFh7uFkwmEwB7Vk0lrEzJGz5Gs57D6Pdgax7u3PyG25KiTaFIRETuWBf2Lna8d6/cid3fTiW8bHC+bQzDwDAMzGZNfsiN0f9BIiJS5Dz2wjs8M/IDWvR6Gc+7OpNwOJldvx2iSY8XKVGjBzHtn2XLzn2O7cdMmU94wz74VutKg87PszPuEACjJ81l3eY9PDnsPbyjujBmynwA1v68i5j2z+JfvTvNeg4j4XCyo61vVm+mUpN+BNTswevvzrm9Jy63lEKRiIgUSfO/+pE3X+nH+T2fExLkT9s+I3n2sU6c3vYZr/7rQbo+9R+ysnIAuKtiKJuXTuLMtvm0bFSbPs9PAGDEsw/RqG5VPh7/HBf2LuaVQT05mnSa7gP/w8SR/TmzfR5d2jak17/GAXD6bBo9nxnHxJH9Sd48h4uZWRw7cabQroHcXApFIiJSJHVp05Do6pVwdnZi2Q+bqVo5jK5t78XJyYkHWjegZIAvP2/7zbFtUIAvLi7OvDKwJzt/O0T6xcwrtjtnyWoeaN2ARvWq4eTkxL/6diTx2EkSj55k+erNRFeLoH2L+ri6ujBqcG9N2xUjWlMkIiJFUmipAMf7I0mnWPvLLvyrd3eU5VqtJJ06C8DHn63gvelLOHbiDCZMGIbB2ZTzeHm6F2j3SNJpPl38Awu/Xucoy8m1cvzkGZJPnaNs6UBHuYe7GwF+3rfi9KQQKBSJiEiRdOlbaABlQgJp3TiaLz95rcB2iUdPMuSNaayZP5461SLIzsnFK7ILhlGwnby2AniiVxsmjXqqQFsHEpNZuXar43NmVjZnUy/cpDOSwqYxPxERKfLaN6/Htj0JLFm5EavVRmZWNivWbCHt/EXSMzIxm80ElfDFarXx2sRP8+1bMsCXxGMnHZ8f6tSURcvWsW7Tbux2OxfSM1i0fD0A9zery9bdB1i+ejM5Obm8/u5c7Hb7bT1XuXUUikREpMjz9fHk6+mv837sVwRH96L8vY/x8WcrAKhWJZz+D7WlZttBlL/3McqHhuDqcnmi5Jm+HYld9B3+1bsz7oMFlC8bwmeTX+alMf8loGZPIlsMYOmqnwAICvBl7qSXePa1DwmJ6Y27myuhIYFX7JMUPSbDuDSAKL8XFxdHVFQUrcKH4GsJ/usdRKRYmPhFd8p4xBd2N0TkOsUdOErVlgPYu3cvkZGR17SvRopEREREUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSEZEirHzDvoQ37ENurtVR9tQrkxn1h3sRifwdCkUiIlKkXUjPZMbCVYXdDSkGFIpERKRIG9zvAcZOmZ9vtOiSD2Z/TcVGj1OyzoM8Mvgt0s5fBGDNTzup1KQfb7w3l4CaPSjfsG++x3ecS71A72fHExzdi4qNHid20Xe37Xyk8CgUiYhIkdbsnpqElQli5h+Cy6p1v/LvSZ/x1fRRHFo/k8ysbAa/Ps1Rn3jsJBZXF05u/Yzhg3oy4JVJjrpHh7xNqeAAjmycxbIZr/PKmzPZGXfotp2TFA6FIhERKfJGPte7wGjRvK9+5MlebYmqFIanhxv/ebEv87/+kUsPcvD0cOPFAV1xdnbi4c7NOHL8NKlp6Zw4dY41P+9k7Et9sVhcuCuiLL06NWXxig2FdXpymzj/9SYiIiJ3thYNa1EmJIDYzy+PFiWfPEuD6MuPeShXpiRZ2Tmc+99T7YNK+GI2540NeLi7AZCekUXSybNkZedSsk4vx742u52HOjW9DWcihUmhSEREioWRzz3EU6+8T9N7agDkTX8dP+2oP5J0CjeLKyX8vP+0nTIhAXh5unFu5wJMJtMt7bPcWTR9JiIixUKrRnUICfLny2/znmjfs31jPpm3grgDR7iYkcWrb8+iR7tGfxl0yoQEck/tSF59exYZmVlYrTZ+3X2AvfuP3I7TkEKkUCQiIsXGyOceckyP3de4Di8P7EG7vq8R3rAvLs5OTHyt/99q59P3XuR48hkqNupHcHQvhrwxjcys7FvZdbkDmIxLK84kn7i4OKKiomgVPgRfS3Bhd0dEbpOJX3SnjEd8YXdDRK5T3IGjVG05gL179xIZGfnXO/yORopEREREUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhEREREAnAu7A3cqm80GQHrO2ULuiYjcTvsTDnDe/Whhd0NErlPC4WTg8t/xa6FQdBWJiYkA/JQ0q3A7IiK31aoOEwu7CyJyEyQmJlKtWrVr2sdkGIZxi/pTpOXk5PDtt98SHh6Ok5NTYXdHbpP09HTq1avHpk2b8PLyKuzuiMgtoJ/z4s1ms5GYmMh9992Hq6vrNe2rUCTyO+fPn8fX15e0tDR8fHwKuzsicgvo51yuRgutRURERFAoEhEREQEUikTysVgsvPbaa1gslsLuiojcIvo5l6vRmiIRERERNFIkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkxcjMmTPx8/O7rcds2rQpgwcPvq3HFJE7W9++fXnggQcKuxtyHRSKpEjp27cvJpMJk8mEq6srERERvPHGG1it1sLumojcBJd+vq/2GjVq1C05roKMgB4IK0VQmzZtmDFjBtnZ2SxfvpxBgwbh4uJCqVKlCrtrInKDkpOTHe/nz5/PyJEjiY+Pd5T9/lllhmFgs9lwdtafMrk5NFIkRY7FYiEkJIRy5crx9NNP07JlS5YuXVpgu4SEBDp16kRwcDBeXl7UrVuX7777Lt822dnZDBs2jLJly2KxWIiIiOC///2vo3737t20bdsWLy8vgoODeeSRRzhz5ky+NqxWK8888wy+vr4EBgYyYsQIfn/7r5SUFB599FH8/f3x8PCgbdu27N+//yZfFZHiISQkxPHy9fXFZDI5Pv/22294e3vzzTffEB0djcViYf369djtdsaOHUv58uVxd3enZs2aLFq0yNGmzWajX79+jvoqVarw3nvvOepHjRpFbGwsX375pWNEas2aNQAcPXqUHj164OfnR4kSJejUqROJiYn52n7++efx8/MjICCAl156Cd3+r+hSKJIiz93dnZycnALl6enp3H///Xz//fds27aNNm3a0KFDB44cOeLY5tFHH+Wzzz5j0qRJxMXFMW3aNMe/RFNTU2nevDm1a9dmy5YtrFixgpMnT9KjR498x4mNjcXZ2ZlNmzbx3nvv8c477/DJJ5846vv27cuWLVtYunQpP/30E4ZhcP/995Obm3uLrohI8fbyyy8zbtw44uLiqFGjBmPHjmXWrFlMnTqVPXv2MGTIEB5++GHWrl0LgN1uJzQ0lIULF7J3715GjhzJK6+8woIFCwAYOnQoPXr0oE2bNiQnJ5OcnEyDBg3Izc2ldevWeHt7s27dOjZs2ICXlxdt2rRx/M6ZMGECM2fOZPr06axfv55z587xxRdfFNq1kRtkiBQhffr0MTp16mQYhmHY7XZj1apVhsViMYYOHWrMmDHD8PX1/dP9q1atakyePNkwDMOIj483AGPVqlVX3Hb06NHGfffdl6/s6NGjBmDEx8cbhmEYTZo0MSIjIw273e7YZtiwYUZkZKRhGIaxb98+AzA2bNjgqD9z5ozh7u5uLFiw4JrOXeSf5o8/06tXrzYAY8mSJY6yrKwsw8PDw9i4cWO+ffv162f06tXrqm0PGjTI6Nq1q+Pz73+3XDJ79myjSpUq+X6+s7OzDXd3d2PlypWGYRhGqVKljDfffNNRn5uba4SGhhZoS4oGTcRKkfP111/j5eVFbm4udrudhx56iFGjRrFw4cJ826WnpzNq1CiWLVtGcnIyVquVzMxMx0jR9u3bcXJyokmTJlc8zo4dO1i9enW+NQyXJCQkULlyZQDuvvtuTCaTo+6ee+5hwoQJ2Gw24uLicHZ2pn79+o76gIAAqlSpQlxc3A1fC5F/opiYGMf7AwcOkJGRQatWrfJtk5OTQ+3atR2fp0yZwvTp0zly5AiZmZnk5ORQq1atPz3Ojh07OHDgAN7e3vnKs7KySEhIIC0tjeTk5Hw/387OzsTExGgKrYhSKJIip1mzZnz44Ye4urpSunTpqy6yHDp0KKtWreLtt98mIiICd3d3unXr5hj2dnd3/9PjpKen06FDB8aPH1+gTou6RQqPp6en4316ejoAy5Yto0yZMvm2u/TA13nz5jF06FAmTJjAPffcg7e3N2+99Ra//PLLnx4nPT2d6Oho5syZU6AuKCjoRk9D7kAKRVLkeHp6EhER8Zfbbdiwgb59+9K5c2cg7xfc7xdIVq9eHbvdztq1a2nZsmWB/evUqcPnn39OeHj4n3675Y+/WH/++WcqVaqEk5MTkZGRWK1WfvnlFxo0aADA2bNniY+PJyoq6u+croj8iaioKCwWC0eOHLnqqO+GDRto0KABAwcOdJQlJCTk28bV1RWbzZavrE6dOsyfP5+SJUvi4+NzxbZLlSrFL7/8QuPGjYG8L15s3bqVOnXq3MhpSSHRQmsptipVqsTixYvZvn07O3bs4KGHHsJutzvqw8PD6dOnD48//jhLlizh0KFDrFmzxrH4ctCgQZw7d45evXqxefNmEhISWLlyJY899li+X55Hjhzh+eefJz4+ns8++4zJkyfz3HPPOfrQqVMnnnzySdavX8+OHTt4+OGHKVOmDJ06dbq9F0SkGPL29mbo0KEMGTKE2NhYEhIS+PXXX5k8eTKxsbFA3s/hli1bWLlyJfv27WPEiBFs3rw5Xzvh4eHs3LmT+Ph4zpw5Q25uLr179yYwMJBOnTqxbt06x++IZ599lmPHjgHw3HPPMW7cOJYsWcJvv/3GwIEDSU1Nvd2XQW4ShSIptt555x38/f1p0KABHTp0oHXr1gX+9fbhhx/SrVs3Bg4cyF133cWTTz7JxYsXAShdujQbNmzAZrNx3333Ub16dQYPHoyfnx9m8+UfnUcffZTMzEzq1avHoEGDeO655+jfv7+jfsaMGURHR9O+fXvuueceDMNg+fLluLi43J4LIVLMjR49mhEjRjB27FgiIyNp06YNy5Yto3z58gAMGDCALl260LNnT+rXr8/Zs2fzjRoBPPnkk1SpUoWYmBiCgoLYsGEDHh4e/Pjjj4SFhdGlSxciIyPp168fWVlZjpGjF154gUceeYQ+ffo4puYujU5L0WMytBpMRERERCNFIiIiIqBQJCIiIgIoFImIiIgACkUiIiIigEKRiIiICKBQJCIiIgIoFImIiIgACkUiIiIigEKRiIiICKBQJCIiIgIoFImIiIgA8P/+1ILQtSNL+gAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "ra['Improved'] = pd.Series(\n", + " pd.Categorical(ra.Improved,\n", + " categories=['None','Some','Marked'],\n", + " ordered=True))\n", + "\n", + "props = lambda key: {'color': 'r' if 'a' in key else 'gray'}\n", + "props = {}\n", + "props[('Treated','Marked')] = {'color': '#b35806'}\n", + "props[('Treated','Some')] = {'color': '#f1a340'}\n", + "props[('Treated','None')] = {'color': '#fee0b6'}\n", + "props[('Placebo','Marked')] = {'color': '#d8daeb'}\n", + "props[('Placebo','Some')] = {'color': '#998ec3'}\n", + "props[('Placebo','None')] = {'color': '#542788'}\n", + "\n", + "mosaic(ra, ['Treatment','Improved'], title = \"Arthritis Treatment\", properties=props)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "cfa426fe-0491-4fad-9616-cfb652a03f08", + "metadata": {}, + "source": [ + "### Question\n", + "\n", + "For this trial, what was the null hypothesis?\n", + "\n", + "### Answer\n", + "\n", + "The null hypothesis is that any change in symptoms is independent of whether the\n", + "patient recieved the treatment or a placebo.\n", + "\n", + "### Question\n", + "\n", + "Is it valid to use a $\\chi^{2}$-test for this data?\n", + "\n", + "### Answer\n", + "\n", + "We have more than 5 counts in each cell of the table, so we meet the rule of\n", + "thumb that says we can use the $\\chi^{2}$-test. Be warned, some statisticians\n", + "would prefer that you had at least 10 in each cell.\n", + "\n", + "### Question\n", + "\n", + "How many degrees of freedom are there in this data?\n", + "\n", + "### Answer\n", + "\n", + "Two, because there are two rows of three columns and $(2 - 1)(3 - 1) = 2$.\n", + "\n", + "### Question\n", + "\n", + "Perform a $\\chi^{2}$-test on the contingency table; are treatment and changes in\n", + "symptoms independent?\n", + "\n", + "### Answer\n", + "\n", + "The following code carries out the test and suggests we should reject the null\n", + "hypothesis." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "38bf2706-777a-48de-99a5-b676f4856aed", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Observed values\n", + "[[ 7 29 7]\n", + " [21 13 7]]\n", + "Expected values under null\n", + "[[14.33333333 21.5 7.16666667]\n", + " [13.66666667 20.5 6.83333333]]\n", + "Chi-squared statistic\n", + "13.055019852524108\n", + "p-value\n", + "0.0014626434089526352\n" + ] + } + ], + "source": [ + "t = outcome_tbl.to_numpy()\n", + "col_sums = np.sum(t, axis=0)\n", + "row_sums = np.sum(t, axis=1)\n", + "total_sum = t.sum()\n", + "\n", + "my_et = np.zeros((2,3))\n", + "for ix in range(2):\n", + " for jx in range(3):\n", + " my_et[ix,jx] = col_sums[jx] * row_sums[ix] / total_sum\n", + "\n", + "print(\"Observed values\")\n", + "print(t)\n", + "print(\"Expected values under null\")\n", + "print(my_et)\n", + "print(\"Chi-squared statistic\")\n", + "my_chi = (np.power(t - my_et, 2) / my_et).sum()\n", + "print(my_chi)\n", + "print(\"p-value\")\n", + "print(1 - scipy.stats.chi2.cdf(my_chi, df=2))" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "e326eb76-8b96-4091-8699-fd3fa7c15657", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "13.055019852524108\n", + "0.0014626434089526504\n", + "2\n", + "[[14.33333333 21.5 7.16666667]\n", + " [13.66666667 20.5 6.83333333]]\n" + ] + } + ], + "source": [ + "chi2, p, dof, expected = scipy.stats.chi2_contingency(outcome_tbl.to_numpy())\n", + "print(chi2)\n", + "print(p)\n", + "print(dof)\n", + "print(expected)" + ] + }, + { + "cell_type": "markdown", + "id": "542cd9fa-e404-4408-99e8-1f399ae6c079", + "metadata": {}, + "source": [ + "### (Bonus) Question\n", + "\n", + "What can we conclude from this hypothesis test? Why do we need to randomise the\n", + "treatment?\n", + "\n", + "### Answer\n", + "\n", + "- We can conclude that treatment and changes in symptoms are not independent.\n", + "- Randomisation helps avoid confounding factors which lends more credibility to\n", + " the conclusion that the treatment improves the condition.\n", + "\n", + "Note that a proper treatment of *causality* goes well beyond the scope of this\n", + "course, but recall that randomised controlled trials provide very very high\n", + "quality evidence.\n", + "\n", + "Recall from earlier notebooks the function `estimate_and_ci` which computes the\n", + "probability of success in repeated Bernoulli trials and the $95\\%$ confidence\n", + "interval on this estimate." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "bbf97501-455e-42fd-b4ca-85804deed44f", + "metadata": {}, + "outputs": [], + "source": [ + "def estimate_and_ci(num_trials, num_success):\n", + " p_hat = num_success / num_trials\n", + " z = 1.96\n", + " delta = z * np.sqrt(p_hat * (1 - p_hat) / num_trials)\n", + " return (p_hat,(p_hat - delta, p_hat + delta))" + ] + }, + { + "cell_type": "markdown", + "id": "74bf03e8-b3e0-4fb0-9d7e-3ad40b583773", + "metadata": {}, + "source": [ + "The functions `rand_small_table` and `rand_big_table` defined below return\n", + "random datasets of the same shape as our arthritis dataset under the null\n", + "hypothesis, i.e. when the outcome is independent of treatment. The\n", + "`rand_small_table` returns data from a smaller cohort and the `rand_big_table`\n", + "returns data from a larger cohort." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "52e68788-f037-4a9c-a7e2-286857eef4f7", + "metadata": {}, + "outputs": [], + "source": [ + "_, _, _, expected = scipy.stats.chi2_contingency(outcome_tbl.to_numpy())\n", + "\n", + "def rand_small_table():\n", + " x = np.array(0)\n", + " while x.min() < 1:\n", + " x = scipy.stats.poisson.rvs(mu = np.array(0.5) * expected)\n", + " return x\n", + "\n", + "def rand_big_table():\n", + " x = np.array(0)\n", + " while x.min() < 1:\n", + " x = scipy.stats.poisson.rvs(mu = np.array(1.5) * expected)\n", + " return x" + ] + }, + { + "cell_type": "markdown", + "id": "0eabddc5-8a23-414b-b8a0-90d74fe55312", + "metadata": {}, + "source": [ + "### Question\n", + "\n", + "Using the functions `estimate_and_ci`, and `rand_small_table` and\n", + "`rand_big_table`, demonstrate how the $\\chi^{2}$-test will fail if the cell\n", + "values are too small.\n", + "\n", + "### Answer\n", + "\n", + "The false positive test suggests that the test will not be powerful enough on\n", + "smaller tables." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "748f2fb8-1627-44e7-8f9d-09b2d0af4b2f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(0.0783, (0.07303459542887727, 0.08356540457112271))\n" + ] + } + ], + "source": [ + "num_trials = 10000\n", + "false_pos_count = 0\n", + "for _ in range(num_trials):\n", + " x = rand_small_table()\n", + " _, p, _, _ = scipy.stats.chi2_contingency(x)\n", + " if p < 0.1:\n", + " false_pos_count+=1\n", + "print(estimate_and_ci(num_trials, false_pos_count))" + ] + } + ], + "metadata": { + "jupytext": { + "formats": "ipynb,md" + }, + "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.10.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/example-4/example-4-answers.md b/example-4/example-4-answers.md new file mode 100644 index 0000000..7db17d3 --- /dev/null +++ b/example-4/example-4-answers.md @@ -0,0 +1,241 @@ +--- +jupyter: + jupytext: + formats: ipynb,md + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.14.1 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# Example 4 + +This notebook is available on github +[here](https://github.com/aezarebski/aas-extended-examples). If you find +errors or would like to suggest an improvement, feel free to create an +issue. + +As usual we will start by importing some useful libraries. + +```python +%matplotlib inline +import statsmodels.api as sm +from statsmodels.graphics.mosaicplot import mosaic +import matplotlib.pyplot as plt +import pandas as pd +import scipy +import numpy as np +``` + +Today we will look at a dataset from a double-blind clinical trial of a new +treatment for rheumatoid arthritis. We will test whether treatment is correlated +with a change in symptoms using a $\chi^{2}$-test. + +First, we need to load the data which comes bundled with `statsmodels`. + +```python +ra = sm.datasets.get_rdataset("Arthritis", "vcd").data +ra.head() +``` + +### Question + +Use `pandas` to generate a cross tabulation of the treatment status and +improvement. + +[hint](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.crosstab.html) + +### Answer + +```python +outcome_tbl = pd.crosstab(ra.Treatment, ra.Improved) +print(outcome_tbl) +``` + +### Question + +Generate a mosaic plot to display this data. + +[hint](https://www.statsmodels.org/dev/generated/statsmodels.graphics.mosaicplot.mosaic.html) + +### Answer + +```python +mosaic(ra, ['Treatment','Improved'], title = "Arthritis Treatment") +plt.show() +``` + +### (Bonus) question + +What fundamental error does the default plot from `pandas` make? + +### Answer + +1. It uses red and green as the only colours which can be difficult for people + with colour blindness. +2. It has not respected the implicit ordering of the response values. + +The figure below makes the pattern in the data far clearer (at the expense of a +few extra lines of code). + +```python +ra['Improved'] = pd.Series( + pd.Categorical(ra.Improved, + categories=['None','Some','Marked'], + ordered=True)) + +props = lambda key: {'color': 'r' if 'a' in key else 'gray'} +props = {} +props[('Treated','Marked')] = {'color': '#b35806'} +props[('Treated','Some')] = {'color': '#f1a340'} +props[('Treated','None')] = {'color': '#fee0b6'} +props[('Placebo','Marked')] = {'color': '#d8daeb'} +props[('Placebo','Some')] = {'color': '#998ec3'} +props[('Placebo','None')] = {'color': '#542788'} + +mosaic(ra, ['Treatment','Improved'], title = "Arthritis Treatment", properties=props) +plt.show() +``` + +### Question + +For this trial, what was the null hypothesis? + +### Answer + +The null hypothesis is that any change in symptoms is independent of whether the +patient recieved the treatment or a placebo. + +### Question + +Is it valid to use a $\chi^{2}$-test for this data? + +### Answer + +We have more than 5 counts in each cell of the table, so we meet the rule of +thumb that says we can use the $\chi^{2}$-test. Be warned, some statisticians +would prefer that you had at least 10 in each cell. + +### Question + +How many degrees of freedom are there in this data? + +### Answer + +Two, because there are two rows of three columns and $(2 - 1)(3 - 1) = 2$. + +### Question + +Perform a $\chi^{2}$-test on the contingency table; are treatment and changes in +symptoms independent? + +### Answer + +The following code carries out the test and suggests we should reject the null +hypothesis. + +```python +t = outcome_tbl.to_numpy() +col_sums = np.sum(t, axis=0) +row_sums = np.sum(t, axis=1) +total_sum = t.sum() + +my_et = np.zeros((2,3)) +for ix in range(2): + for jx in range(3): + my_et[ix,jx] = col_sums[jx] * row_sums[ix] / total_sum + +print("Observed values") +print(t) +print("Expected values under null") +print(my_et) +print("Chi-squared statistic") +my_chi = (np.power(t - my_et, 2) / my_et).sum() +print(my_chi) +print("p-value") +print(1 - scipy.stats.chi2.cdf(my_chi, df=2)) +``` + +```python +chi2, p, dof, expected = scipy.stats.chi2_contingency(outcome_tbl.to_numpy()) +print(chi2) +print(p) +print(dof) +print(expected) +``` + +### (Bonus) Question + +What can we conclude from this hypothesis test? Why do we need to randomise the +treatment? + +### Answer + +- We can conclude that treatment and changes in symptoms are not independent. +- Randomisation helps avoid confounding factors which lends more credibility to + the conclusion that the treatment improves the condition. + +Note that a proper treatment of *causality* goes well beyond the scope of this +course, but recall that randomised controlled trials provide very very high +quality evidence. + +Recall from earlier notebooks the function `estimate_and_ci` which computes the +probability of success in repeated Bernoulli trials and the $95\%$ confidence +interval on this estimate. + +```python +def estimate_and_ci(num_trials, num_success): + p_hat = num_success / num_trials + z = 1.96 + delta = z * np.sqrt(p_hat * (1 - p_hat) / num_trials) + return (p_hat,(p_hat - delta, p_hat + delta)) +``` + +The functions `rand_small_table` and `rand_big_table` defined below return +random datasets of the same shape as our arthritis dataset under the null +hypothesis, i.e. when the outcome is independent of treatment. The +`rand_small_table` returns data from a smaller cohort and the `rand_big_table` +returns data from a larger cohort. + +```python +_, _, _, expected = scipy.stats.chi2_contingency(outcome_tbl.to_numpy()) + +def rand_small_table(): + x = np.array(0) + while x.min() < 1: + x = scipy.stats.poisson.rvs(mu = np.array(0.5) * expected) + return x + +def rand_big_table(): + x = np.array(0) + while x.min() < 1: + x = scipy.stats.poisson.rvs(mu = np.array(1.5) * expected) + return x +``` + +### Question + +Using the functions `estimate_and_ci`, and `rand_small_table` and +`rand_big_table`, demonstrate how the $\chi^{2}$-test will fail if the cell +values are too small. + +### Answer + +The false positive test suggests that the test will not be powerful enough on +smaller tables. + +```python +num_trials = 10000 +false_pos_count = 0 +for _ in range(num_trials): + x = rand_small_table() + _, p, _, _ = scipy.stats.chi2_contingency(x) + if p < 0.1: + false_pos_count+=1 +print(estimate_and_ci(num_trials, false_pos_count)) +``` diff --git a/example-4/example-4-questions.ipynb b/example-4/example-4-questions.ipynb new file mode 100644 index 0000000..74be72c --- /dev/null +++ b/example-4/example-4-questions.ipynb @@ -0,0 +1,282 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a61f4b1e-5b27-4cdc-83af-1119d0f692a9", + "metadata": {}, + "source": [ + "# Example 4\n", + "\n", + "This notebook is available on github\n", + "[here](https://github.com/aezarebski/aas-extended-examples). If you find\n", + "errors or would like to suggest an improvement, feel free to create an\n", + "issue.\n", + "\n", + "As usual we will start by importing some useful libraries." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "abb9087e-7cbf-43fd-9fd3-6465f613cf03", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import statsmodels.api as sm\n", + "from statsmodels.graphics.mosaicplot import mosaic\n", + "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "import scipy\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "id": "65f9801b-b722-4312-b720-2edb27a6d175", + "metadata": {}, + "source": [ + "Today we will look at a dataset from a double-blind clinical trial of a new\n", + "treatment for rheumatoid arthritis. We will test whether treatment is correlated\n", + "with a change in symptoms using a $\\chi^{2}$-test.\n", + "\n", + "First, we need to load the data which comes bundled with `statsmodels`." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "de17ffee-8a0c-4730-949a-a7c621c5a1b5", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
IDTreatmentSexAgeImproved
057TreatedMale27Some
146TreatedMale29None
277TreatedMale30None
317TreatedMale32Marked
436TreatedMale46Marked
\n", + "
" + ], + "text/plain": [ + " ID Treatment Sex Age Improved\n", + "0 57 Treated Male 27 Some\n", + "1 46 Treated Male 29 None\n", + "2 77 Treated Male 30 None\n", + "3 17 Treated Male 32 Marked\n", + "4 36 Treated Male 46 Marked" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ra = sm.datasets.get_rdataset(\"Arthritis\", \"vcd\").data\n", + "ra.head()" + ] + }, + { + "cell_type": "markdown", + "id": "c3729f85-5d76-4768-9c94-24533dbd34f1", + "metadata": {}, + "source": [ + "### Question\n", + "\n", + "Use `pandas` to generate a cross tabulation of the treatment status and\n", + "improvement.\n", + "\n", + "[hint](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.crosstab.html)\n", + "\n", + "### Question\n", + "\n", + "Generate a mosaic plot to display this data.\n", + "\n", + "[hint](https://www.statsmodels.org/dev/generated/statsmodels.graphics.mosaicplot.mosaic.html)\n", + "\n", + "### (Bonus) question\n", + "\n", + "What fundamental error does the default plot from `pandas` make?\n", + "\n", + "### Question\n", + "\n", + "For this trial, what was the null hypothesis?\n", + "\n", + "### Question\n", + "\n", + "Is it valid to use a $\\chi^{2}$-test for this data?\n", + "\n", + "### Question\n", + "\n", + "How many degrees of freedom are there in this data?\n", + "\n", + "### Question\n", + "\n", + "Perform a $\\chi^{2}$-test on the contingency table; are treatment and changes in\n", + "symptoms independent?\n", + "\n", + "### (Bonus) Question\n", + "\n", + "What can we conclude from this hypothesis test? Why do we need to randomise the\n", + "treatment?\n", + "\n", + "Note that a proper treatment of *causality* goes well beyond the scope of this\n", + "course, but recall that randomised controlled trials provide very very high\n", + "quality evidence.\n", + "\n", + "Recall from earlier notebooks the function `estimate_and_ci` which computes the\n", + "probability of success in repeated Bernoulli trials and the $95\\%$ confidence\n", + "interval on this estimate." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "bbf97501-455e-42fd-b4ca-85804deed44f", + "metadata": {}, + "outputs": [], + "source": [ + "def estimate_and_ci(num_trials, num_success):\n", + " p_hat = num_success / num_trials\n", + " z = 1.96\n", + " delta = z * np.sqrt(p_hat * (1 - p_hat) / num_trials)\n", + " return (p_hat,(p_hat - delta, p_hat + delta))" + ] + }, + { + "cell_type": "markdown", + "id": "74bf03e8-b3e0-4fb0-9d7e-3ad40b583773", + "metadata": {}, + "source": [ + "The functions `rand_small_table` and `rand_big_table` defined below return\n", + "random datasets of the same shape as our arthritis dataset under the null\n", + "hypothesis, i.e. when the outcome is independent of treatment. The\n", + "`rand_small_table` returns data from a smaller cohort and the `rand_big_table`\n", + "returns data from a larger cohort." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "52e68788-f037-4a9c-a7e2-286857eef4f7", + "metadata": {}, + "outputs": [], + "source": [ + "_, _, _, expected = scipy.stats.chi2_contingency(outcome_tbl.to_numpy())\n", + "\n", + "def rand_small_table():\n", + " x = np.array(0)\n", + " while x.min() < 1:\n", + " x = scipy.stats.poisson.rvs(mu = np.array(0.5) * expected)\n", + " return x\n", + "\n", + "def rand_big_table():\n", + " x = np.array(0)\n", + " while x.min() < 1:\n", + " x = scipy.stats.poisson.rvs(mu = np.array(1.5) * expected)\n", + " return x" + ] + }, + { + "cell_type": "markdown", + "id": "0eabddc5-8a23-414b-b8a0-90d74fe55312", + "metadata": {}, + "source": [ + "### Question\n", + "\n", + "Using the functions `estimate_and_ci`, and `rand_small_table` and\n", + "`rand_big_table`, demonstrate how the $\\chi^{2}$-test will fail if the cell\n", + "values are too small." + ] + } + ], + "metadata": { + "jupytext": { + "formats": "ipynb,md" + }, + "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.10.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/example-4/example-4-questions.md b/example-4/example-4-questions.md new file mode 100644 index 0000000..8cf5dd5 --- /dev/null +++ b/example-4/example-4-questions.md @@ -0,0 +1,127 @@ +--- +jupyter: + jupytext: + formats: ipynb,md + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.14.1 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# Example 4 + +This notebook is available on github +[here](https://github.com/aezarebski/aas-extended-examples). If you find +errors or would like to suggest an improvement, feel free to create an +issue. + +As usual we will start by importing some useful libraries. + +```python +%matplotlib inline +import statsmodels.api as sm +from statsmodels.graphics.mosaicplot import mosaic +import matplotlib.pyplot as plt +import pandas as pd +import scipy +import numpy as np +``` + +Today we will look at a dataset from a double-blind clinical trial of a new +treatment for rheumatoid arthritis. We will test whether treatment is correlated +with a change in symptoms using a $\chi^{2}$-test. + +First, we need to load the data which comes bundled with `statsmodels`. + +```python +ra = sm.datasets.get_rdataset("Arthritis", "vcd").data +ra.head() +``` + +### Question + +Use `pandas` to generate a cross tabulation of the treatment status and +improvement. + +[hint](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.crosstab.html) + +### Question + +Generate a mosaic plot to display this data. + +[hint](https://www.statsmodels.org/dev/generated/statsmodels.graphics.mosaicplot.mosaic.html) + +### (Bonus) question + +What fundamental error does the default plot from `pandas` make? + +### Question + +For this trial, what was the null hypothesis? + +### Question + +Is it valid to use a $\chi^{2}$-test for this data? + +### Question + +How many degrees of freedom are there in this data? + +### Question + +Perform a $\chi^{2}$-test on the contingency table; are treatment and changes in +symptoms independent? + +### (Bonus) Question + +What can we conclude from this hypothesis test? Why do we need to randomise the +treatment? + +Note that a proper treatment of *causality* goes well beyond the scope of this +course, but recall that randomised controlled trials provide very very high +quality evidence. + +Recall from earlier notebooks the function `estimate_and_ci` which computes the +probability of success in repeated Bernoulli trials and the $95\%$ confidence +interval on this estimate. + +```python +def estimate_and_ci(num_trials, num_success): + p_hat = num_success / num_trials + z = 1.96 + delta = z * np.sqrt(p_hat * (1 - p_hat) / num_trials) + return (p_hat,(p_hat - delta, p_hat + delta)) +``` + +The functions `rand_small_table` and `rand_big_table` defined below return +random datasets of the same shape as our arthritis dataset under the null +hypothesis, i.e. when the outcome is independent of treatment. The +`rand_small_table` returns data from a smaller cohort and the `rand_big_table` +returns data from a larger cohort. + +```python +_, _, _, expected = scipy.stats.chi2_contingency(outcome_tbl.to_numpy()) + +def rand_small_table(): + x = np.array(0) + while x.min() < 1: + x = scipy.stats.poisson.rvs(mu = np.array(0.5) * expected) + return x + +def rand_big_table(): + x = np.array(0) + while x.min() < 1: + x = scipy.stats.poisson.rvs(mu = np.array(1.5) * expected) + return x +``` + +### Question + +Using the functions `estimate_and_ci`, and `rand_small_table` and +`rand_big_table`, demonstrate how the $\chi^{2}$-test will fail if the cell +values are too small. diff --git a/example-4/README.org b/old-example-4/README.org similarity index 100% rename from example-4/README.org rename to old-example-4/README.org diff --git a/example-4/cat-grades-by-breed.csv b/old-example-4/cat-grades-by-breed.csv similarity index 100% rename from example-4/cat-grades-by-breed.csv rename to old-example-4/cat-grades-by-breed.csv diff --git a/example-4/cat-weights-by-breed.csv b/old-example-4/cat-weights-by-breed.csv similarity index 100% rename from example-4/cat-weights-by-breed.csv rename to old-example-4/cat-weights-by-breed.csv diff --git a/example-4/cat-weights.csv b/old-example-4/cat-weights.csv similarity index 100% rename from example-4/cat-weights.csv rename to old-example-4/cat-weights.csv diff --git a/example-4/example-4-answers-grades.ipynb b/old-example-4/example-4-answers-grades.ipynb similarity index 100% rename from example-4/example-4-answers-grades.ipynb rename to old-example-4/example-4-answers-grades.ipynb diff --git a/example-4/example-4-answers-weight.ipynb b/old-example-4/example-4-answers-weight.ipynb similarity index 100% rename from example-4/example-4-answers-weight.ipynb rename to old-example-4/example-4-answers-weight.ipynb diff --git a/example-4/example-4-data-grades.R b/old-example-4/example-4-data-grades.R similarity index 100% rename from example-4/example-4-data-grades.R rename to old-example-4/example-4-data-grades.R diff --git a/example-4/example-4-data-weight.R b/old-example-4/example-4-data-weight.R similarity index 100% rename from example-4/example-4-data-weight.R rename to old-example-4/example-4-data-weight.R diff --git a/example-4/example-4-questions-grades.ipynb b/old-example-4/example-4-questions-grades.ipynb similarity index 100% rename from example-4/example-4-questions-grades.ipynb rename to old-example-4/example-4-questions-grades.ipynb diff --git a/example-4/example-4-questions-weight.ipynb b/old-example-4/example-4-questions-weight.ipynb similarity index 100% rename from example-4/example-4-questions-weight.ipynb rename to old-example-4/example-4-questions-weight.ipynb diff --git a/requirements.txt b/requirements.txt index 82fc2d9..e1a0e67 100644 --- a/requirements.txt +++ b/requirements.txt @@ -60,7 +60,6 @@ patsy==0.5.2 pexpect==4.8.0 pickleshare==0.7.5 Pillow==9.3.0 -pkg_resources==0.0.0 pkgutil_resolve_name==1.3.10 prometheus-client==0.14.1 prompt-toolkit==3.0.31 diff --git a/resources/check-questions-text.py b/resources/check-questions-text.py index 3fd0437..851caab 100644 --- a/resources/check-questions-text.py +++ b/resources/check-questions-text.py @@ -16,7 +16,7 @@ def check_file(filepath): def main(): offending_files = [ - "./example-{n}/example-{n}-questions.md".format(n=n) for n in range(4) + "./example-{n}/example-{n}-questions.md".format(n=n) for n in range(5) ] checks = [(check_file(fp), fp) for fp in offending_files] print("\n")