{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Around Simpson's Paradox" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Data are available in the MOOC repository. The CSV file contains data for the 1314 women that were polled in Whickham, England, in 1972-1974 and were categorized as \"currently smoking\" or \"never smoked\". Each line is related to a person and contains whether she smokes or not, whether alive or dead twenty year after the survey and her age at the time of the survey." ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of rows: 1314\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
SmokerStatusAge
0YesAlive21.0
1YesAlive19.3
2NoDead57.5
3NoAlive47.1
4YesAlive81.4
\n", "
" ], "text/plain": [ " Smoker Status Age\n", "0 Yes Alive 21.0\n", "1 Yes Alive 19.3\n", "2 No Dead 57.5\n", "3 No Alive 47.1\n", "4 Yes Alive 81.4" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data_url = \"https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/-/raw/master/module3/Practical_session/Subject6_smoking.csv\"\n", "data = pd.read_csv(data_url)\n", "print(\"Number of rows:\", len(data))\n", "data.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The CSV file does not contain any missing data." ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of rows with missing values: 0\n" ] } ], "source": [ "print(\"Number of rows with missing values:\", data.isnull().any(axis=1).sum())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "Let's visualize the number of women alive and dead after twenty years, according to their smoking habits. A heatmap is effective in this case." ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAHwCAYAAABuVI8jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XWYFdXjx/H32V526ZLu7hKQVMAAUSxsvyoGKip2d3ciovITCxFEDMRACaWREESQDsmltoPdnd8fZ3a5d3fvBizuDnxez8PD3Zk5M2fmTnzmTFzjOA4iIiIiXhFU0hUQERERKQqFFxEREfEUhRcRERHxFIUXERER8RSFFxEREfEUhRcRERHxlFIXXowx44wxT5d0PQrDGHO1MWZOSdfjeGWMuckYs9sYk2CMqVzEstnrkTGmlzHmn2NTS79p1jfGOMaYkGMwbs9sF8XNGHOeMWabux50KOn6yLFhjJlljLnO/Xy5Mebnkq5TSclvezfGvGuMeeQYTPOo9l9u2cYB+hX791lgeDHGbHYPIFE+3a4zxswqzopIySitB0VjTCjwKnC64zjRjuPsO9JxOY7zu+M4zYqxbu8ZY24orvFJgV4GRrjrwTJ3n9S/uEZujAk3xow1xmwxxsQbY5YZY87KMUw/Y8waY0ySMWamMaZecU1fcnMc5zPHcU4/FuMu7vXnv+Y4znDHcZ4q6XoURc7vM7+gU1iFbXkJAW4/mgmVBGNMcEnXoTgcizN5D6gORACrSroieTgTmFbSlTiB1KOY1gNj5dzvhQDbgD5AeeARYKIxpr5bpgrwldu9EvAH8EVx1OdIBZgPTzhB92dS3BzHyfcfsBm4H9gPVHC7XQfMcj/XBxwgxKfMLOA69/PVwFzgNeAgsBE4xe2+DdgD/M+n7DjgXWA6EA/MBur59G/u9tsP/AMMzVF2NPbAkgj0z2N+rgFWu+PeCNzo068v8C9wl1uvncA1Pv0rA98CccAi4ClgToDllrVcbgB2uOO6y6d/kLtcNwD7gIlApRxlhwFbgd/c7j2Bee5y3AZc7XYPx56dbgV2u8svsqB5cut2CEgDEoDv3O5Z9YoH/gbO86l3MPAKsBfYBIzw/f6xO/+x7nS2A08DwQGWUTjwurt8drifw4Gm7vfnuPWaEaD8JGAXEAv8BrTKsS487bsMfObtyxzjeQN4szD1B9oCK3yWxcvustgI3FLYZQE0Ama43/1e4DPc7cvt3wFY6n4HXwATsuYnj+WwBejkfr7CrUNLn2316/yWd4715F6f9WQIMBBYi93eHizi+vs/7Dq5F3gon33MIGAZdrvaBjzuU98Ed1yJ7rQ+ATKBZLffve6w3Ti8bfwJ9M2xP3oGux9KBhoXYr+3ArjAZzuZ59Mvyh1P8zzKXQQsydHtrhzfQaBttSIwFYgBDrifa+c3H9j96EZ3PdkEXF7QvBW0H3T7nwssd7+TDcCZbvdKwIfu+nMga77cftcD69115Vugpk8/B7t9rAM2ud0GAGuw2+/b2H2973FjTo7yw93yB4BRgCnMPinHfAVaf87BBuSD7nJuEWC5GeyxbI9b7xVAa599zjvAD+645wInYbezA+68dvAZVwt3WgfdaZ8TYP9VFpgJvOlO37dfXwo+Zn3nfo+Lsfuggo5ZeW63wMnAfLe+O93vLCzHd3Qbdn3aC7wEBOX8PrH76qxtOgG4GKiCXd8PYtef37PKBlyHC7GSbwb6Y888shZYUcNLOnZjCXYX3lbsyhcOnI7dgKJ9vrR4oLfb/w2fmY7C7tyuwZ4tdXQXUiufsrFAD+zONSLAjrKRuxL0AZKAjj4rQjrwJBCK3XEnARXd/hOwO+kooDX2gFTQivC5O3wb7E6pv9t/JLAAqO3O5xjg8xxlP3bLRgJ13eVyqVu3ykB7d/jXsTuLStgV/TvguULO0zhyHBSxO+Ca7jK8GLuS1XD7DccGmtrYne0v+B+wv3bnJQqohg15NwZYRk+6y6AaUBV78Hkq0HqVR/lr3fnNOigvD7Dx9+VweKnnzn85nx3fTqBbYeqPPWA/57Ms1gB13GU/s7DLAnvgGeDWvSp2g37d7ReGDSR3uN/ZhdiQGSi8fIwbjIH3sAebm3z63VGI5Z21njzqTvN67Po63l3GrYAUoGER1t/3setuOyCVwAeEvtjtIwgbDncDQ3LsFBv7/L0ZnxMToBY2QA10xzHA/buqz/5oqzsPIUBoAfu86u68Nnf/fgMYnWOYv3DDTY7u4didbwufbss4HITy21YrAxcAZdx+k/APBznnozz2oNTM7V8DnwBfwDzmtx88GbsfHeAuz1o+y+J7bJiu6K4nfdzup2H3xR3dZfAW7kmXz3c43Z3vSOzBKg67bodi1/V08g8vU4EK2H1hDIcDVb77pEDHNJ+/s06WBrh1uRcbwsLyKHsGsMSth8EGkKx94zh3GXTCthrPwIapqzh87JvpDhvqTuNB7PZ+Gnb/3sxnXE+768QifLZ9cu/bCjpmTcCuUy2xx8+Cjll5brfufHXDrnv1seF3ZI7vaKb7HdfFnvTk9336btPPYYN8qPuvF244DbgOF2Il34wNL62xK3RVih5e1vn0a+MOX92n2z4OH4jHARN8+kUDGdgDxMXA7znqNwZ4zKfsx4XZeH3Kfw3c7rMiJOeYlz3uFxaMPYA09+n3bCFWBN/hXwTGup9XA/18+tVwxx/iU7ahT/8HgCl5TMdgN7xGPt26c/jsJuA85dwQ8llGy4Fz3c8z8D+Y98/6/rE7/VTcM0m3/6W4G2we490ADMyxY9gcaL0qoI4V3OHLB9jA//UZdg5wlft5ALDB/Vxg/bFnBL18lsVwn36nH8WyGAIscz/3xp7ZGp/+8wJ9T9gWum991qvrcLchbAjKOijlt7yz1pOslqGy7rx09Rl+CW6ooHDrr2+rwSLgkkJ+l68Dr/n8XVB4uQ/4JMc4fsJt0cXuj54s5LRDsQe/MT7dxgLP5xhuLm7LZx7jGA08435uhT3rDqeAbTWP8bQHDvj87Tcf2FB8EBt4Igszf/nMt+9+cIzv8s/xHWfiHhhz9BsLvOjzd7S7PtT3+Q5P8+l/FbDA52+DbUHI72DX0+fvicD9PtthnvukAPOac/15BJjo83cQ9sS0bx5lT8MelLuRo2UAu8953+fvW4HVPn+3AQ66n3thW42DfPp/zuFWx3HA/2FD8j15TMd331bQMauZT7/CtLwUarvFnsBM8fnbwQ2U7t83A7/m8336btNPAt9QiFbRrH+FvmbqOM5f2OR7f2HL+Njt8znZHV/ObtE+f2/zmW4C9kymJvasuasx5mDWP+BybNNcrrJ5McacZYxZYIzZ75YfiD0LyLLPcZx0n7+T3LpV5fC18Sxb8ptWHvXZ4s4H7rxM8ZmP1diQVj1A2TrYg09OVbGpeonPuH50uxc0T3kyxlxljFnuM77WHF5GNXPUy/dzPezOf6dP2THYM/281MR/Gfoun3wZY4KNMc8bYzYYY+KwOyTw/y4DGY8NEgCXuX8XWH9jTAXsZct5PvUPtD4UNK5qxpgJxpjtbv0/xX8Zb3fcrTqPcec0G+hljDkJu8P6Aujh3rNRHhs+s8ab3/Le5zhOhvs52f0/0HZamPV3l8/ngOucMaarexNsjDEmFnsmXZjvMUs94KIc+4We2INtlnz3C249grCXFdKwlx6yJADlcgxeDnumnJePgMuMMQa4EntgTKWAbdUYU8YYM8a9cTgO2xpXIce9e777xkTsCd1w7Hr2vTGmeUHz6U4rv/1goH1NHWC/4zgH8ujnt265++192FabXHUnx7bjrusFfUeB1qf89kmFkbPume44auUc0HGcGdjLJaOA3e7N+77rRs7tJdD2UxPY5k4ry5Yc0xyEbQF5t4D6F+WYVZhlk+dyNsY0NcZMNcbsctfPZ8m9nQY63hXkJWxL1M/GmI3GmAJzRlFv+HoM25zsu4AT3f/L+HTzDRNHok7WB2NMNLYZagd2wcx2HKeCz79ox3Fu8inrEIAxJhyYjL3mXN1xnArY+2NMIeoUg22eq+PTrW5R5sUdfof7eRtwVo55iXAcZ7vP8L7zsg3bzJvTXuxG0cpnPOUdxwkYTnLwW17uUxTvY3feld1l9BeHl9FObPNsXvO3DdvaUMWnLuUcx2kVYNo7sAeeLL7LpyCXYa/L98ceoOtnzUIhyk4C+hpjagPncTi8FFT/M7BnElkH+J0EXh8KGtdz2GXf1nGccth7VXyXcS334JfXuP04jrMeu5O5DdtUH4/dAd2APdvJ2kEezfLOqTDrb2GNx15KqeM4Tnnszjq/7zHnNr4N2/LiW5cox3Gez6eMH3dZj8WGrwscxznk03sVtgk9a9go7LaY503EjuMswAagXtj19BO3V0Hb6l1AM2xrVzlsCxz4Lwu/+XAc5yfHcQZgg9oa7Labr0LsBwPta7YBldwQn5PfuuUuo8rYFoy86u637bjL33dbKor89kl5ybku5Kx7Vl3yXJcdx3nTcZxO2Fa1psA9Ra2wO806OW66rptjmu9jw+0036d9iyDrmFWUZZOf0dh1rIm7fj5I7u000PEuX47jxDuOc5fjOA2BwcCdxph++ZUpUnhxd5JfYHeSWd1isAv8Cvds+FryXvGLYqAxpqcxJgx7U+xCx3G2YVt+mhpjrjTGhLr/uhhjWhRyvGHY5tsYIN19HLJQj+O5B6yvgMfdM6SW2BubCvKIO3wr7L06WU8pvAs8k/XIpTGmqjHm3HzG8xnQ3xgz1BgTYoypbIxp7x6Y3gdeM8ZkndXXMsacUZj5wp4ZNPT5Owq7cce447oG2/KSZSJwuzuNCtgmewAcx9kJ/Ay8YowpZ4wJMsY0Msb0CTDtz4GH3Xmvgr3f4tNC1rssNhzswwbnZwtZLmudnYW98XCT4zirC1n/Qfg/ZTQRuM0YU9sYUxGfVslCjKss9oz+oDGmFv47wPnYnc5t7nd9PvY+hPzMxgbO2e7fs3L8DUe3vHMq6vqbn7LYM/oUY8zJ2AN+fnKus58Cg40xZ7j7oAhjTFY4LazR2PsXBjuOk5yj3xSgtTHmAmNMBHa5rXAcZ00+4/sYe4ae7jjOHMg+o89vWy2LDTcHjTGVsCeLARljqhtjznEPbKnY9SnD7Zf1zo76eRQtaD84FrjG2MfDg9w6NnfX6R+Ad4wxFd39b1bAGu+Wae+Go2ex++3NAar/PdDKGHO+sU8f3caRn/QG3CcFkHP9mQgMcuc3FBsiUzncwprNPd50dYdLxN4blZFzuEJY6Ja/112OfbEH7Qk5hhuBfTBlqjEmsigTyOOY1Rx7ue5IlcXep5TgjuumPIa5x1036mCfUA70VJ7fd2CMOdsY09gNjnHYZZrvcj2SR+2exB7gfF2P3fnuw6bRXF96EY3Hbrj7sTcJXQ42nWE3skuwiW4X8AJ2QyyQW/427Mp6ALuT/LYI9RqBbULbhb3u+GEhyszGNof9CrzsOE7Wi3recKf9szEmHnvzY9d86r4V27R7F3a5LOfw2eB97jQWuM15v2DP4ApjLNDS2Gbsrx3H+Rt75/587ArWBnt9P8v72IPyCuyNiNOwB9qsFe0q7M7xb+wy/hL/5ntfT2MfO10BrMQ+XVPYd858jG2W3O5Oa0Ehy2UZj221GZ+je571dzeqAdgzoSzvY++t+NOt+1eFGZfb7wnszY2x2B15dlnHcdKA87HXiQ9gLw3kHHdOs7E7l98C/A1Ht7xzKtL6W4CbgSfd8TyK3T7z8xw2hB00xtztnticiz0TjMG2ENxDIfdvbgC7EXuPyS5jX4aXYIzJ2u/EYO8reQb7fXTF7oPy8wk29H+So3t+2+rr2MsEe7HL80fyF4TdH+zA7hP6YJcl2DPgrO3DT0H7QcdxFmFPtF7Drp+zOdwycSX2Poo12HsrRrplfsXeOzIZ2xLSiHyWkeM4e7EPBjyPPW40wX8/UxQF7ZNyyrn+/INt+XwLu+wHY0NsWh5ly7nTO4BdvvuwLVhF4o77HOAsd5rvYO/DW5NjOAfbgroN+MYNz0UxAtsyvQu7Ln6ODWZH4m7suhKPXQZ5BZNvsPfGLcfu18YGGNfjwEfudzAU+/3/gg3g84F3HMeZlV9lsh41k2LmnvFswj7ZkJ7/0N7lnrW96zhOvQIH9jC3ReBtx3EKagERwT1L3oO9WXpdCUz/YSDGcZwx//W0S9qJsk86EsaYF4CTHMcpzFWDUk0vC5IicXfKp2LPdKpjW8imlGil/jv5NuOL+LgJWFwSwQXAcZxS99bsY+UE3yfly728E4Ztae2CfTLxuhKtVDFRy8sxcry2vBhjymCbkZtjr89/j33EMq5EKyZSShhjNmNvZBziOM6yEq7OcU/7pMCMMV2wl4pqYlsCx2Af+/f8gV/hRURERDzFk7+NISIiIicuhRcRERHxlBP2ht3IDiN0vUzkGJox6YS5Z1SkxHRvXKEwL+Y87qjlRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8JaSkKyDHpysGd+X9J68scLiBw99i5sJ/8uxXrVJZ7ry6P2f1ak2dkyqSnHqI1Rt28unUhYybMj/gOBvWqcLZfdrSu3MT2jStRfXKZUnPyGTHnljmLlvPexN/Z9nqbUc8byKlQUJcLMsW/sbfy/9gy4Y17N2zi8yMDMqWr0D9Ji3o2W8QnU7pm2fZNSuX8tfShWxev5qYXduJj40lNSWJMtHlqFW3AZ1O6UufM84lLDwi3zrs2fkv0778hL+WLeTg/n1ElomiXqOm9DlzCF16nHYM5lrEMo7jlHQdSkRkhxEn5oz/R7LCS0ZGJjEHEgIPd+9Y5i7bkKt7hxZ1+HbULVSpGA1AfGIKEWGhhIYGAzB93mouHDmGtEPpfuW6t2vIjHF3+nWLS0gmPCyE8LBQADIyMnlh7E88Nfr7o5pHyd+MSU+XdBWOa8POOYWMjIzsv0PDwgkKCiI1JTm7W9vO3bnlgecJj/APIa89fid/Lp6b/Xd4RCSAX9mq1Wty11NvcFKtunlO/8/Fcxn13IOkpaYAEFkmipSUZJzMTAB6DTiba29/GGPMUc6p5Kd74won5AJWy4scU//uPkDzQY8VqUy56AgmvzGcKhWjWbNxF8Me+Zilf28lNCSYa8/vwYt3n8+AU1rw4t3nM/K5iX5lQ0KCSU/PYNpvfzHhh8XMXryO/bGJBAUZOjSvwwt3nU+Pjo158Iaz2LpzPx99HbgFR6Q0y8jIoGHTVvTsP4jWHbtRrUYtAGJ27+C7CR/y28/fsuKP+Yx7+zluvPsJv7It23ehdcduNG3Vjmo1ahNZJgqwrTnzZ/3EpHFvE7N7B28+fS9PjxpPUJD/HQYxu3bwzvMPkZaaQpOWbRk28hFOqlWXlOQkfpj8Kd98Ppbfp0+lRu36DLyw4BZYkaLSPS9S6oy8qh81qpYnKTmNIbeOZunfWwE4lJ7BmIm/8dS70wAYdn4PGtet5ld247YY2l/wNBff9T5TflnO/thEADIzHZb8vZWzbnyLFWv/BeCea07/D+dKpHjd9+woHn3t/zht0AXZwQVsi8m1tz9E37POA2D+zB/ZF7Pbr+wZQy5lwDlDqdeoWXZwAYguV54B5wzlsuvvAGDH1k1sWLMy17S/+nQMqSnJlK9YmZGPvZrdOhMRWYbzrriBvmcOAeC7Lz4kMT6ueGdcBIUXKYUuP7srAJN+WsKWHfty9R/9+SziE1MICQnmkoGd/fpt33OQDVtjAo77UHoGE75fDECjulWpUDayGGsu8t9p0a5zvv17n35O9ufN61YXadyNmrfO/rx/7x6/fqkpySyZNxOA0waeT1R02VzlBw39HwDJSYksXTC7SNMWKQyFFylVmtSrRt0alQD4ee7feQ6TmJyWfZ9M/+4tijyNlLTD98kEB2sTkONTaFhY9ufMzIx8hsztn1XLsz9Xq1Hbr9/aVX+SlpoKQJvOp+RZvmr1mtSsUx+Av5YuLNK0RQpD97zIMVWlYjRzP7uXpvWrExxk2LU3jgV/buTDKfP5fcm6XMO3alwz+/OqDTsCjvfv9Ts4s2crmjc4qch16t25CQA7Y2LZdzCxyOVFvGDNiqXZn2vXb1zg8GmpKezfu4fFc2bwzedjAWjWugMNmvifIPy75fAN9rXqNgw4vlr1GrFj22a2b91U1KqLFEjhRY6pqMhwOrasy/7YRKIiw2hQuwoNalfh0kEn89HX87nl6c/JyMjMHr5G1fLZn3fsiQ043qx+5ctGEhUZRmJyWqHq07VtAwb3bQvAh1PmHcksiZR6iQnxfD/pIwCatmpPjdr18hzu4P59jLxyYJ792p/ci+vufCSPMnsBiIoul+spJl8VK1d1hw98GVfkSCm8yDGxMyaWp9+dxjczlrN28x7SDqUTFGQ4uU19Hh4+iH7dmvO/Id1JSknjzhcmZZeLLhOe/TkpJXAg8e1XNiqiUOGlSsVoPnruaoKDg1i3ZQ+vjpt+hHMnUnplZmby3iuPc3D/XkJCw7hi+F0Bhw0KDqJcBXuZNjkpkUNp9nJQl579OO+KG4guWz5XmZQk21oZFh6eq5+vrHfEpCQlHdF8iORH4UWOiV8XrOHXBWv8umVmOiz4cxODbx7FF69cx+BT23HDRb0Y9fmsfG+yLQ5RkWF8+fqN1KtZmbiEZC6/d2yhW2tEvOSzMa/y56I5AFx18z3Ubdg04LDlylfkzc9+AMBxHA7s28PMaVP4ccp4li6YzZXD785+akmkNNHdivKfcxyH+1+bAtgbZgf1bpPdLyEpNftzmYiwXGXz6hefmJLv9MpEhDHlrZvo2rYB8YkpnHfraFau3X6k1RcptSZ88Aa/TrUtmZdeP9LviaOCGGOoVKU6F1w1nBvvfoKM9HQ+eudFtm5c6zdchPtoddZNu4FkvbwuokyZosyCSKEovEiJ2LhtLzEH4gFoUKtydvedMYfvc6lZLXeTdc5+sfHJ+bagZAWXXp2akJCUynm3jWbe8o1HW32RUueL/3uLH6eMB+Dia2/ljCGXHvG4Ovc4lcrVauBkZvLbz9/59atQqQoAiQlxpKYEPnE4sC/GHb7qEddDJBCFFylVVq0//IRRq0Y1Aw7X0n0qac2mXQGHyQouvTs3ITHZBpe5S3P/FIGI100Y+yY/TP4UgKHXjuCsC6446nFWdEPKnp3+vwNWu16j7M/btwY+EdjuPpVUq26Do66LSE4KL1IiGtSuQtWK9uVWm31eRLduyx627twPwIAeeb/DpUxEGD062B3oL/PzfvlWmYgwvn7bBpeEpFSG3DqaOUvWF+csiJQKEz54gx+/+gywwWXgBUf/On7HcYjZbU8kIiKj/Po1bdUu+2bdlUvy/nmNvXt2smPbZgBad+x61PURyUnhRUrEcyPt68MzMjKZ9ttffv0+m2pfanXRGZ2yX1jna/jFvSkbFUF6egYTpv2Rq39WcMm6VDTk1ncUXOS4NOGDNw5fKhp2W6GCS0ZGeoHD/D79O2IP2JOK5m06+vULj4ik0ymnAjBz2lckJeb+4dVpkz4B7M8FdOzWp8DpiRSVwosUu7o1KvH7J3cz7IIe1Pe5n8UY+6j0N2/fzLn92gPwweQ5rNvi//rx1z/+lZ0xsURFhjPlrZvo0KIOAKEhwVx/UU8evXkQAGO/msv6rf5lIyNC+erN4fTq1IT4xBSGjHhHl4rkuDTx/97ODi6XXjeSs86/vFDl1q76k2fvvZG5M6axf6//bx7t2r6ViR+O4qO3nwfs23V79h+UaxznX3Ej4RGRHNy/l9efuItd2+3vj6WmJPPN+A+Y+cNXAJxzybVElS13xPMoEohxHKek61AiIjuMODFn/D9Qt0Yl/pn2ZPbfKamHiE9KoWyZCCLCQ7O75/WSuiwdWtTh21G3UKViNABxCclEhIcSFmqf7p8+bzUXjhxD2iH/s8jLzj6ZsU9dBUByShqxCfk/iXTp3e+z4E+9AfRYmDHp6ZKuwnFr355d3HXNuQCYoCDKlquQ7/BnnX959n0wq1cs4YUHbs7uFxoWTkREJKmpyX5PENVp0ITbHnmRqtXzvvfsz8VzGfXcg9lPFUVGRZOanJz9UwQ9+w9i2MhHMMYc+YxKgbo3rnBCLmC950WK3Z798dzx/ES6tm1A22a1qVIxmoply5CSdojN2/exYMVGPv56AfP/DHyz37LV2+h04TPcdc0AzurVmtrVK5CYnMbilZv5dOpCPvp6AXkF7yCfHWVkRBiR+TxuDRAaok1AvCfTORz4ncxM4g7uz3f4lJTk7M/1Gzfn+rseY82KpWxev4bYA/tIjI8lJDSMajVqU69RMzr3OJUuPU4jKDg44DjbdenB06M+4/tJH7Nq+SIO7ttLmeho6jVqRt+zzqMvksy/AAAgAElEQVRLj9OOfkZFAlDLi4gcE2p5ETn2TtSWF93zIiIiIp6i8CIiIiKeovAiIiIinqLwIiIiIp7i6fBijKltjJlijIkxxuw2xkw2xtQu6XqJiIjIsePp8AJ8CHwL1ABqAd+53UREROQ45fXwUtVxnA8dx0l3/40DAv6EqTHmBmPMH8aYP9L3rvrvaikiIiLFxuvhZa8x5gpjTLD77wpgX6CBHcd5z3Gczo7jdA6p0uo/rKaIiIgUF6+Hl2uBocAuYCdwodtNREREjlOefje64zhbgXNKuh5e1L55bQb2bkOHlnVpUrcqVSqWpVxUBHGJKazdvIsf5/zN+5N+50BcUsBxnNe/PZef3ZWOLetSuUIUh9Iz2L77IHOWrmfMF7+xYu32I6pbwzpVOLtPW3p3bkKbprWoXrks6RmZ7NgTy9xl63lv4u8sW72tUPN42xWn0btzE6pUjGZ/bBKLV27mnQmzmb14bcByHVvW5ZGbBtG9XUNCQ4L5e8MOXvvoV776ZVnAMr07N+Gn92/n+9kruXDkmCOabzn+pKak8M9f9jX8W9b/w+b1a9gXswuAcy+7jvMuvz7PcjG7d3DPtecVejo9+w/iujseLVLd3n/1Seb++n2Bw439di7BwYEPFZvXr+GnKZ+zZuUS4mMPElW2HI2at6b/4KG0bNc5YLlN61Yz5dMxrPt7BRkZ6dSq14izzr+Ck3v1C1hm9Z9/8MKDt9Du5J7c8dgrBdZdjl+eDC/GmPy2UsdxnKf+s8p41P/O7c7wSw7/VH1yShrJqWlUrhBF9/aN6N6+ESMu78tFI99j4Qr/Hy4MCw3hs5eGcXafNtnd4hNTCAsNpmn96jStX53/ndudB16bwlufzSxSvbq3a8iMcXf6dYtLSCY8LIQm9arRpF41rhzcjRfG/sRTowPveK8+rztvPnAJoaH2t1kOxidRvXJZzjmtHeec1o6n353GM2Om5SrXpXU9fnr/diIjwkhPz+BQegadW9fns5eGMfK5iYyZ+FuuMmGhIbz10CXEJ6Zwx/MTizS/cnzbuHYVrz52R5HLBQUFU65CpXyHOXQojeTEBAAaNGl5RPUD+8OMkWWiAvY3BH77/OyfvuHjUS+QkWF/jDEyKpq4g/tZOn82S+fPDhjQNqz5i+cfuJlDaakEBQUTHBLCprV/887zDxIXezf9z74oV5lDh9L4aNQLRESW4aqb7jmCOZXjiSfDC5CYR7coYBhQGVB4KcDiVVvY8uoU5i3fwD+bdhObYH+4LSoyjCH9O/DsyCFUq1SWL169nrZDniTO59eZ7x12enZwefeL33hp7E/siInFGEP75rV56e4L6NGxMc/feR5zl21g6d9bC12vkJBg0tMzmPbbX0z4YTGzF69jf2wiQUGGDs3r8MJd59OjY2MevOEstu7cz0dfz881jq5tG/DWg5cQEhLMtzP+5M4XJrF9z0EqlY/i8RGDuf7Cnjw8fCBrNu5k8nT/1pRn7ziPyIgwxk9dxK3PTCAl7RC3XNqXF+++gCdvHcxnUxeSkJTqV+a+686gaf3q3PPSl2zbdaDQ8yonhqjoctRr1Ix6jZtRr1EzPn//dWIPBLw1D4DKVavz5mc/5DvMJ6Nf5tepkwgLD6db3zOOuH4n9+rP9XcWrdUGYP3qlXz09gtkZmbQsXsfrhh+F5WqVCchLpYvPx7NrB+m8M34D6hVtwEn9+rvV/aL/3uLQ2mpdD/1TK4ecT+hoWFM/3Yin3/wOl+Oe4cepw3MFai+m/Ahu7Zv5dLrR1K52klHPL9yfPDkPS+O47yS9Q94D4gErgEmAA1LtHIeMX7qIl7/5FcWrdycHVwAEpPT+Oy7hVz70EcAVK9cjoG9WvuVvfzskwH47Y913PH8RHbExALgOA7LVm/j/NvfJT4xhaCgIIb0a1+kem3cFkP7C57m4rveZ8ovy9kfa3NqZqbDkr+3ctaNb7Fi7b8A3HPN6XmO45nbzyUkJJiVa7dz+X1j2b7nIAD7YxO57ZkJ/Dz3bzvcyCEEBfn+CnUo3ds1JD09g5HPTyQpJY3MTIe3PpvJstXbKBcdSde2Dfym1axBde66uj9L/t7KOxNmF2le5fjXrFV7Rn0xnXuffZuLr72Vbn1OJyQ0/186L4y0tFQWzPoJgE6nnEpUdNmjHmdRTfzwLTIzM6hdvzE33/8slapUByC6XHmuHnE/rTt2A+CL/3ubTLdlBuyltHWrVxAUFMxVN99DeEQkQcHBnHHepdRr1IyU5CQ2rFnpN60dWzcx7ctPqN+4OQMGD/3vZlJKLU+GFwBjTCVjzNPACmwLUkfHce5zHGdPCVftuLBo5ebsz7WqV/Drd1KV8gABW1TiElJYt8V+DdFlwos03e17DrJha0zA/ofSM5jw/WIAGtWtSoWykX7969eqTI+OjQF4/ZNfSU/PzDWOlz78GYB6NSvT0x0WoGLZMgQHB7H3YCLxiSl+ZTZstfNTtWK0X/e3H76U4KAgRjw1nsxM/VC5+AsKDj4m410ybxaJCXEA9Dnj3GMyjfzs2bmdtav+BOCs8y8nJCR3I/7ZQ/8HwL49O/nnr8MtnEkJcTiZmUSXK09kGf/tqXrNOgDExR7M7uY4DuPefp7MzEyuvvWBY7ZMxVs8GV6MMS8Bi4F4oI3jOI87jqP2+mLUo0Oj7M8b/93r12/Tdvt3hxZ18ixbLjqCJvWqAbCkCJeMCislLT37c3Cw/yrcr1vz7M/T3RaWnOYt20Cc29rUv3uL7O4H4pPIyMikSoUoykZF+JVpUMe+PijmQEJ2t2vOO4WeHRsz6vNZLF/z7xHOjUjR/fbzt4A92Ddv0/E/n/6q5QuzP7fp1C3PYZq2bEdEZBkA/lp2ePgy0eUwQUEkxMWSnJTgV2bPTrsdlSt/+IRp9k/fsHbVcgacezH1GzdHBDwaXoC7gJrAw8AOY0yc+y/eGBNXwnXzrLDQEOrWqMTwi3sz9umrAFi/dQ/fz/7Lb7j3J80BoE+Xprx2/1BqVi2f3a9989p89cZwykZFsHDFJiZMW1zs9ezduQkAO2Ni2XfQ//anVo1rArB7X5xf0PCVmemwdvNuAFo0qpHdPTnlEAv+3EhISDCv3X8RkRGhBAUZbr60D51a1iU+MSW7RapapbI8ffu5bN25nyffKfiJDZHismfndtasWAJA79OP/mHL1X8u5r7rL+S6Ib0YfuGpPHzzZXz23qvs2h74xOPfzRsBKFehYsAbi4OCg6lRpz4A27ccvuk/PCKCJi3akJmZwSejXyY1JYXMjAymf/sFm9evISKyDI2a20vVsQf2MenDUVSuehLnX3HDUc+rHD88ecOu4zheDV2l0oEFrxERHpqr+7xlG7j6wXGkHUr36/7uF79Rq1oFRl7Vj+EX92b4xb2znzYKDwtlZ0wsL/3fzzz73g9kZOS+bHM0urZtwOC+bQH4cMq8XP1ruEFqx57YfMeT1b+GT/ACeOiNb/jxvdu4/OyuXHxmZw6lZxAZYe9ReGLU1OzLSS/dcwGVykcx7OGPSUpJO7qZEimC36d/h+M4BAcH06PfoKMe3/69ewgKCiayTBTJSYn8u2UD/27ZwMzvv+KyG+7gtEEX5CpzcL+9tFuxcrV8x12xclU2+QyfZeg1t/L8Azczb8YPLJj1M8EhIRxKszfCn3/ljdmXk8a/9xqJCXFcf9djhEdE5hy9nMA8GV6keO3eF0d4WCjRZcKz71GZtegfHnrjmzyfnnEch0fe+pbVm3bx2n0XUTYqwu8yS0R4KOWiI4iKDCMl9VCx1bNKxWg+eu5qgoODWLdlD6+Om55rmKz6JxcQKLICR9kc9+QsXLGJAcNe59Gbz6Zr2/qEhgSz5O+tvPHxr0z6yZ7tDjilBUPP7Mzkn5fy4xz7MxM3XdKH6y/qScPaVdh3MJEpvyzj8VFTcz2ZJHI0MjMymPPLVADadulBhUqVj3hc9Rs3o0HTFrQ/uSeVKlcjKDiY1JQUVi6Zz8QP32bPzn/5+J0XKVu+Al16+r97JSXZvv8pLDz/e9rCwiP8hs/SuEUbHnzhXb76dAzrV68kIyOD+o2bc+b5l9Otj70Rf8Uf81n423S69OxH+5N7AjD924nMnDaZPTu3E12+Al16nMr5Vw7P91FvOT4pvAjNBz2W/blqxWguO/tk7h12Br9/cjfPf5D7fSqVK0Tx2YvD6NOlKb/MX82zY35g1YYdRIaH0rVtA56+fQg3Du3NgO4tGDDs9eynkY5GVGQYX75+I/VqViYuIZnL7x1LYvKxafH4Y9UWzrllVJ79IiNCeeOBizkYn8TdL30JwLMjh3DH//qzffcBJv64hLbNanPLZafSsWVdBlz3RrG3PsmJa8WS+RzYZ1sxjvZG3QHnXJyrW3hEBJ17nErzNh15fOT/2Lt7JxPGvkXnHqdhTOD3vRyJhs1acfdTb+bZLzUlhU/eeZHIqGguv9G+92nC2Df58avPqFi5Kl37nM7WjWuZ/u1ENq1bwwMvjM73RXpy/NHlF/ETcyCBNz6Zwbm3vIPjwIM3nMVZOR6Vfv/JK+nTpSm//bGOwTePYv6fG4lLSGH3vni+nbmCfte8SsyBeBrWqcpTtx39kxBlIsKY8tZNdG3bgPjEFM67dTQrA7y9N6ulI+tST37jBIgvYsvIQzcMpEHtKjz8xjfs2htHk3rVuP3K09i1N47ul77ADY99Ss8rXmT24rV0b9+Iq87J+2ZGkSORdaNuxcpVadup+zGbTnS58gy++BrAPi20ZcM/fv2zbsRNS81/+0lLTfEbvrC+Gf8BMbt3MPSaW6hQqQo7/93CT1PGU65CJZ5482Ouv/NRHnt9HM3bdGT96hX8Pn1qkcYv3qfwInn6Y9UW5i3fAMC1F/TI7t6sQfXsMPPGJ7/mWTbmQALjpy4C4Nx+7Y6qHlnBpVenJiQkpXLebaOZt3xjwOF3uq08NauVDziMb/+dRWgVat2kJrddcRrzl29g7OS5AAzu25agoCAmTFucfYNwenomb4+3bxYefGrbQo9fJD9xB/fz52K73vXsf/Yxf2S4cfPDb9CO2eV/slChkn367sC+/N9MkdVKlDV8YWzbtI6fvh5P4xZt6Xum/YmEZQt+w3EcTjn1zOwbhENCQjj93EsAWDpf71g60Si8SEBZN7U2qlMlu1uLhoefzsn5CLWv9e67WqIiw6lW6cheoJUVXHp3bkJisg0uc5duyLfMqvU7APtyvSo53smSJSjI0LS+faHW6g07C1UXYwyjHr4UB4dbnvo8u3v92nbZbPzX/4bErPlvUKsKIsVhzq/fk5GejjGGXgMGl2hdate37wKNO3iAuNi831KRmZHBzm2bAahVr0Gew+Qqk5nJh289B8Zwza0PZF+q2uOGp2o1avsNX72WfV1DzK4dRZ4H8TaFFwmoQS17M2BC4uGm4czMw/dv1K0R+LdXfAPLkdy0WiYijK/ftsElISmVIbeOZs6S9QWW+3XBmuzPA05pkecw3ds3pFy0fXLhl/mrC1WfG4f24uS2DXh13C+s3rgrV//I8LAcf9untxz04jopHr/9/B0Azdt2olqNWsd8euvXHH5FQtXqNf36tWrfNfvzyiW5f6IDYN3qFdk36rbu0DXPYXKa8f1kNv6zioEXXEmterlflp6W5r8vOeRetirm23HEAxReTkC+r8QPpO/JTencuh5gfwYgy/I1h3/N+fqLeuVZtkxEWPZPCKxY+2+RHyXOCi5Zl4qG3PpOoYILwObt+5i71A57+5X9CAnJvYrf7f6swJYd+5iztODx1qxansdvGcz6rXt4/oMf/fptcV/Yl7WsspzcpkF2fUSO1tpVy9n17xYA+hTDu10cJ/9QnRAfy9SJ4wCoVKUadRs18+tfrUYtmrayl4R//Go86enpOUfB95PsT4xUrlaDZq07FFinA3v3MPnj0VSvWZvBl1zj1y8rPG1cu8qv+/p/bMCqkiNcyfFP4eUEVLt6RRZMuJ9hF/Sgfq3KOfpV4O5rBjDptRsJCgpi38FE3vpsRnb/rTsPMHW2/d2Rs/u0YexTV9HAvXQSEhJEt3YN+PmD22novpH2jU9mkNN7T1xB8rK3SV72dq5+kRGhfPXmcHp1akJ8YgpDRrxT4KWinB564xvS0zNo16w2nzx/bfZL9CqWK8PrDwzlzJ6t7HCvf12oV/q/et9FlC8bya3PTCA1zX8nnbUszj2tHeeeZnfmLRqexD3XDrD9Z/n/RouceBLj44iPPZj9z3Fs62Vaaopf95yPE/ua/ZO9UTeqbDk69Ti1UNN9/9UnuXpQV64elLvVY97MH3jr6ftYPHcGcQf3Z3dPS01hyfzZPHXnsOz7XC4edhtBQbkPFUOvuZWgoGC2bVrH6Bce4sBee/9LQnwsH496kRV/2BaZi68dUaj7cz599xWSkxL53y33Exbm/wh2h672RGnJvFn8MXcmjuOwfctGvp9oA1LHbr0LsUTkeKJny05Q7ZrV5u2HLwUgNe0QcYkpRIaH+f0W0aZ/93Lp3R+we1+8X9nhj3/KN6NuoVPLulx29slcdvbJJCanEhYSQmjo4Z3Uax/9kn3jbmGd178Dfbo0BSAkOIhPXxyW7/CX3v0+C/7c5Ndt4YpN3PrsBN584BKG9GvPkH7tORCXRPnoiOyd8NPvTsv1i9J5ObtvG87t155Pv1vIrEVrc/X/Z9NuRo2fyS2XncqEV64nKTmNMpH2EtKiFZv4+Nu8m9TlxPHobVexb0/ue6t+mPwpP0z+NPvvHv0G5fnrzslJiSyeY2+OP+XUMwkthh92zMzIZMn8WSyZPwuA8IhIQsPCSEpIIDPT/ohiSGgYl153O117D8hzHI1btOF/I+7j41EvsGTeLJbMm0WZqLIkJyVkt+yce9l1uX5ROi9L589myfxZ9Og3kJbtu+TqX7NuAwacM5Tp307k7WfvJyw8PPtJp4bNWtGzhO8Bkv+ewssJaGdMLJffM5ZenZvQpU19alQpR+UKUWRkOGzduZ8Va7czddYKvvjhjzxfMrfvYCJ9rnqZKwZ35fz+HWjbrDaVypchPSODbbv2s+DPTYydPCffp4ICCTK+v/IcVuAjz6F5/CAcwLgp81m+ehu3X9mPXp0aU6ViNHv2J7BoxSbemTCb2YtzB5GcoiLDePXei9h7IIH7X/0q4HB3vzSZLTv2M+yCHjSoXYWdMbF8/etyHh/1XZ4/DClSFAtn/5z9yHHvYvoRxhbtOnHBVTexfs1Kdm7bbH9nKDGByDJRVKtZmxZtO3PqWedR9aT8L8f0OeNc6jVqxo9TxvPPyqXExx6kbPmKNG7Rhv6Dh9KyXecC65KSnMSnY14hulx5Lhl2e8DhLrvhTqpUq8nMH6cQs2s75StWpnOP07jgquF5/jCkHN9MQdc+j1eRHUacmDMu8h+ZMenpkq6CyHGve+MKJ+TtyrrnRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDxF4UVEREQ8ReFFREREPEXhRURERDylxMOLMSbYGPNLSddDREREvKHEw4vjOBlAkjGmfEnXRUREREq/kJKugCsFWGmMmQ4kZnV0HOe2kquSiIiIlEalJbx87/4TERERyVepCC+O43xkjIkE6jqO809J10dERERKrxK/5wXAGDMYWA786P7d3hjzbcnWSkREREqjUhFegMeBk4GDAI7jLAcalGSFREREpHQqLeEl3XGc2BzdnBKpiYiIiJRqpeKeF+AvY8xlQLAxpglwGzCvhOskIiIipVBpaXm5FWgFpALjgThgZInWSEREREql0tLyUs1xnIeAh7I6GGO6AItLrkoiIiJSGpWWlpevjDG1sv4wxvQG/q8E6yMiIiKlVGkJLzcCXxtjTjLGDATeBAaWcJ1ERESkFCoVl40cx1lsjLkN+Bn7UwEDHMeJKeFqiYiISClUouHFGPMd/o9ElwFigbHGGBzHOadkaiYiIiKlVUm3vLxcwtMXERERjynR8OI4zuysz8aY6kAX989FjuPsOZbTDm/Z7ViOXuSE16F+hZKugogcp0rFDbvGmKHAIuAiYCiw0BhzYcnWSkREREqjkr5slOUhoEtWa4sxpirwC/BlidZKRERESp1S0fICBOW4TLSP0lM3ERERKUVKS8vLj8aYn4DP3b8vBn4owfqIiIhIKVUqwovjOPcYYy4AegAGeM9xnCklXC0REREphUpFeAFwHGeyMWY6bp2MMZUcx9lfwtUSERGRUqZUhBdjzI3Ak0AykIltfXGAhiVZLxERESl9SkV4Ae4GWjmOs7ekKyIiIiKlW2l5omcDkFTSlRAREZHSr7S0vDwAzDPGLARSszo6jnNbyVVJRERESqPSEl7GADOAldh7XkRERETyVFrCS7rjOHeWdCVERESk9Cst97zMNMbcYIypYYyplPWvpCslIiIipU9paXm5zP3//hzd9ai0iIiI+CnRlhdjTBdjzEmO4zRwHKcB8ATwFzAV6FySdRMREZHSqaQvG40B0gCMMb2B54CPgFjgvRKsl4iIiJRSJX3ZKNjnJwAuxv6m0WRgsjFmeQnWS0REREqpkm55CTbGZAWoftjHpbOUdLASERGRUqikA8LnwGxjzF7s7xr9DmCMaYy9dCQiIiLip0TDi+M4zxhjfgVqAD87juO4vYKAW0uuZiIiIlJalXTLC47jLMij29qSqIuIiIiUfiV9z4uIiIhIkSi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKcovIiIiIinKLyIiIiIpyi8iIiIiKeElHQF5MQxcnArHr+kQ/bfFS7/tNBlX732ZK7t1xSArTEJtB35dZ7DndmhFj1aVKd9g0rUrhxF5XLhRIQGsy8+lb+2HuDrBVuYMGcTGZnO0c2MSAlLTk5myR+L+HvVKlav/pvVq1axc+cOAIbfPIKbbrk1YNk/Fi9i/ry5rPrrL7b/u40DBw+QnJRE2XLlaNSoMaf1H8AFFw4lIiIi4DgyMzP5YdpUvv3ma/5Z/TcJCQlUqFiRDh06csllV9Cpc5din2eRLAov8p9oXKMc953f5ojK9mxRnatPbVKoYR+9uAMt61TI/jsuOY2MTIcaFctQo2IZBrSrxfWnN+OiF2cSE5dyRPURKQ3+WrmCW4bfcERlx304lt9nz8r+OzKyDKGhYRzYv58/9i/ij8WL+OyTj3hnzAfUr98gV/mkpCTuvP1W5s+bA0BwcDBRUdHs27uXn3/6kek//8SNN92Sb4ASORoKL3LMGQNvXd+NyLAQFq6NoWvTqoUuGxkWzJvXdSM9M5O/Nh2kY6PK+Q7/7aKtvPvTGhaujWHLngRSDmUAcFKFSK46tTH3nd+G9g0qM3r4KVz44oyjmi+RklauXHlatGxJixYtad6yFS+/8Bx798YUWK5bt+6c0qMnHTp2om7dukRFRQNw8OABpn0/lTdefZnt//7LHbeNYPLX3xEU5H+HwVNPPMr8eXMICgpixG0jueSyy4mKiiYuLo6PPxzL+++9y7vvvE29evUZePbgYzLvcmJTeJFj7sbTm9G9WTW+mLOJTbvjixReHhnanoYnleWlr1dSq1KZAsPL81+tyLP7roPJvDhlJeGhwdx1bmv6t6tJzUpl2LE/qUjzIlJadOzUmd/nL/Lr9uZrrxSq7BVXXZ1n9woVKnLZ5VcSFhrGU088ysYN6/lz+TI6dOyUPcy6dWuZNvU7AC67/EqGXX9jdr9y5cox4vY72L17N99+M4VXX3mRAaefQWhYWBHnTiR/umFXjql6VaN4eGh79sWn8OCnfxSpbOfGVbjxjGas2xHHy1+vLJb6/LF+b/bnGhUji2WcIiUhODj4mI27bbt22Z93797t18/3ctP/rh2WZ/mrryEm8PUAACAASURBVL0OgJg9e5g3b27xV1BOeAovcky9cV03oiNCeejTJeyLTy10ubCQIN6+oRsGwx3/t5DUQ5nFUp/uzaplf968J6FYxilyvFm6ZEn25zp16vj127nD3hRctmxZqlWrnmf5uvXqERJiG/bnzvn9GNVSTmS6bCTHzFWnNqZv6xrMXLmTCXM2Fansvee1oXmtCnw0cx1zVu8uuEA+osJDqFs1ikt6NWTEwBYAfP7bhiKFKZHjXUpKCrt372L6Tz8yZvQoADp17kKr1nnfaJ+RkRFwXI7jkJlpTzjWr1tb/JWVE57CixwTNSpG8tSlHUlKTWfk2IVFKtu2XkVuP7sVuw8m8+j4ZUc0/c6Nq/DLE2fm6p6ekclnszZwz7jFRzRekePJ3pgY+vXtmWe/Pn1P5alnns/VvWatWoB94mjHju3UrFkr1zAbN2zIDi8xe/YUY41FLF02kmPitWFdKR8VxvNfrWBLTOEvzwQHGd6+oTuhIUHc9/FiYpPSjmj6h9Iz2X0wmd0Hk0lLP3yG+OGMdbzw1Yrsp5BETmRBwcFUrlyFypWrEB4ent399DPO5I677qF8hQq5yvTs3Sf78/vvjs5zvB+8d7h7QqIuz0rxU8uLFLuhPRpwZofarNi8n1HTVhep7B3ntKJt/Ur8uPRfvl649Yjr8Ofm/TS7ZTJgH9WuXy2am89qwbX9mnBJz4bc8M5cflj67xGPX+R4UKlSJWb8Zm+odRyHPbt3M+mLz/n4ow+Z8euvPPDQI1w49GK/Mk2aNOXMswbx4w/f89XkSZQpU4bLrryKatWqs3PHDsZ9OJbpP/9ESEgo6emHCDI6R5bip7VKilWVcuE8d2Un0jMyue2DBUV6k22zWuW5Z0gb4pMPcdeHiwouUEiOA5t2J3DPuMU8On4ZZSNDee/mHlSvoKeNRLIYY6h+0kmMuP0Onn3hZdLTD/HMU4/zz5o1uYZ97Imn6HZKDwA+/eQjBp7ej87tWzN44OlMnvQF7dp34KyBgwD7+LRIcVN4kWL1xCUdqVw2gnEz1rFuRxxR4SF+/0JDDq9y2d2CbbeXr+5CeGgwr3zzFwcT03KVDXZflGXM4bIhwaZI9fvgl39IScugbGQoF3avX2zzLXI86T/gdGrWrEVmZiZTvvoyV/8yUf/f3p3H6VQ9cBz/nNlX+8yUdeyM7GtS2qyFRIVUSimlffkpftosrRIivxKyRJIWqZAipcJYspTs+zIYZl/P749n5mHMDDMM0x3f9+vl9Xrm3nPuc66XZ3yfc88SyPgJH/LmyFFcf2MbKoWHU7ZsOZo0bcbzg4bw0ZRpREW5liWoFB5+kVsvlwI9NpICVSnEtVLn/W1qcn+bmmcsu/ejHgCM/3YTz09b5a77Uo+GWfZAOl2FMkHuus9PXcn477J/M8xNUko6x+KSuNwngMphQXmuJ3KpCQkNZd++vezetTPH8x4eHrRt14G27TpkO5eSnMyf69YA0KBhowvaTrk0qedFLilBfl6UCXZtNhebmFrIrRH5d7LWsnePa0xYQGBgvusvWPAdsbGxeHl50eEmbQ8gBU89L1Kgbh628IznB95aj4Hd6gHZd5XObafoTOMevJJe11TNdVdpTw9z1jE2j90U4X50db7rx4g4UWpqqnsBudx8MXeOe4+kpk2b5ev6hw8fYtTItwC45dZuhIXlvJCdyPlQeJEi4/arKtO5aUVmLN3K8s2HiDrhWoTOGIgoX4IH2takz/Wu3amX/32IRWv3FWZzRc7biePHSUs/Oe0/c22VxIQEjh076j7u6+Pr7kFZHbmKcWNHc2v322jWrAVhl13mLrdz5w7mzvmMqVMmAVChQkU633JrtvdduuQndu/ayTWtr6NsuXJ4enoSHx/Pkh8X8+47b3Po4EEqV6nCU888d0HuW0ThRYoMY6BD4/J0aFwegNjEFPfgXF/vk/vALFm/n3tGa8lycb47undl37692Y5PnjSRyZMmun/u3KUrrw4/ueBc5KqVRK5y7TXm6+tLQEAACQkJJCYmusvUrFmLd8a8h5+fX7br79q5gzdfH8Ebrw3Hy8uLgIBAYmJOYK2r57Nuvfq8O2ace7dqkYKm8CJFxver9/LYh7/RqnYYV1QsSWhxP0oE+pCQnMaOQ7Gs3naEOct3sFA9LnIJi6hTh6EjXmflij/YtGEDUVFRHD8ejbe3DxUqVKR2RAQ3tGlHm7btct38sUXLq+jV+y5WR0ZyYP8+YmJiKVOmDLUi6tCh40106HgzHh4aUikXjslMypeaEndOuzRvXOQiOTCld2E3QaTI8/Mif+tFFBGKxiIiIuIoCi8iIiLiKAovIiIi4igKLyIiIuIojg0vxpgBxphiGa8nGGP+MMbcUNjtEhERkQvLseEF6GetPWGMaQuUA/oDbxRym0REROQCc3J4yZzq3AGYZK1dxVnuxxjTzxiz0hizMnnL4gveQBERESl4Tg4va40x84FOwLfGmCBOBpocWWv/Z61tYq1t4lPt+ovSSBERESlYTl5h916gMbDFWhtvjCkD9C3kNomIiMgF5tjwYq1NM8ZUAdoAwwB/nN2TdFHVDy9F+0blaBBemqqXB1Mm2I9gf29iElLYvP84C9fsY+KizUTHJWere1WtUK6rezkNq5QmPDSI0sG+BPp6Ex2XzF97o5m3cjdTFm8hMSUth3c+s1a1w5g3uE2ey4+Ys5bXP/8zy7GKZQJpUTOUBpVLUT+8FHXDS1LM3weAeo/PZVdU3Bmv2aByKV7oXp/mNULw9vRg055oRs/byJd/7Mq1ztURYXw9qA3fRe6hx9s/5bn9UrQlJCSwauUfbNywgU2bNrJpwwb273dtT/HQwwPo/8ijZ73GkagoJn30IUuX/MiB/fvx9fWjarVqdO7Sla7dumPMuS2w+t8XBvLVl3PPWm7V2g257kKdmJjI7Fkz+WHRArZt3UpcXCx+fn5UqFiJq1pdTc9evSkTEpJj3Q3r/2Tc2NGsWR1JamoqVatVp899fWnbrkOubfnj99944L57uKb1dYwZ937eblSKJMeGF2PMWMAbuAZXeIkD3geaFma7nOLO1lXp17am++eE5FQSUlIpFexLi+BQWtQIpX/7WvR8+ydWbInKUvfRmyNo37C8++fYxBSSUtMIKe5HSPHLuDriMvq3r0W31xez9UBMvtqVnJrGweiEM5YJ8PUi2N8bgMitR7KdH9itHr2uqZqv983UuGpp5g1ug7+PF6lp6aSkpdO4ahmmPH4Nz0z+gw8Xbs5Wx8fLg5H3NicmIYVnJv9xTu8rRdP6P9fxyEP9zrn+xg3r6d+vL9HR0QAEBAQQHx/H6shVrI5cxYLvv2P0e+Px8fE55/fw9fUlKCg41/O5haN9+/by4P33sWvnDvex4OBg4uLi2LRxA5s2bmDWJ9N5d+x4GjfJ+mt53bq13N/nLpKSkvD09MTLy4sN6//k2aee4Oigo/TodWe290tOTmboKy8SEBDAC4OHnNvNSpHh2PACtLTWNjLGrAaw1h41xpz7J/gSE7k1isHTY/lt82H+2Xec4/EpAAT6etG5WUVe6dmIkOJ+TH+qNU2e/ooTCSnuukvWH2Dxuv389vchth2MITYxFYCSQT7c1rIyL/VoSHhoMNOebE3LgfPIz/ZZf/wTRc1H5pyxzMynr6V9o/LsPRrHD+v2Zzufnm7ZdiCGtTuOsmb7EYwxvNSjYZ7e/9VejfD38WLWsm08+dHvJCan81C7mgy/qwlD7mjAzJ+3ue830zNdrqB62WI8P3Ule47E5/1m5ZJQrFhxakdEULt2BLUi6vDW6yOIijp81noxMTE8+vBDREdHU7lKFYaNeIM6V9QlJTmZOZ/N5s3XR7D812W8+dpwBg156Zzb1659xyw7TufV4Of/w66dO/D29uaZ556nc5dbCAgMJCU5mV+W/cywoS9z6OBBnn3qCeYv+CHL7tTvvPUGSUlJ3NSpM4OHvIyvry8zpk3lrTdGMHrU23Tq0iXbjtQfTBjPzh07ePY/z3N52bLnfL9SNDg5vKQYYzzIGKRrjCkNpBduk5xj5rLtOR6PS0rlk5+3cSA6gbkDbyC0uD/tGpZj9q873GXGf/dXjnWPxSbzvwV/k5SSxrv3t6B2+RI0qx7C75vP/os6ry4r4c+N9V2/uGYs2UZ6DsnosQ9/z3K8Ve2wPF3b38eT5jVCSE1L55nJK4hPcj32GvfdX9zRqgr1K5eiafUQfvzzZGCqUbYYj3eqw+ptR5jw/d/nc2tSBDVq3ISfl2ftjRv9ztt5qjtl0kSiog7j5+fH2PH/o3z5CgB4+/jQo9edxMXFMnrUSOZ89il33n0P4eGVC7z9udm3by+rVq4AoO8DD2bpKfH28eHa62/APyCAfn37cORIFKtWruCqVlcDrkdpa1ZH4unpyQuDXT0pAHfd04d5X3/JX5s2snbNGlpe1cp9zW1btzJp4gdE1KlDzzvvumj3Kf9eTh4j8h4wBwgxxrwMLANeL9wmFR0rT3lUVLZUwEWreza9rqmCl6cH6emWaUu25Fgmp0CTFyUCffD08OBITBIxp/Q0AWw9eAKAMsG+WY6P6tscTw/DExN/P+f3laLL09PznOvO++pLANp16OgOLqfq2as3AQEBpKWlMX/e1+f8Puci6vDJLyQRda7IscwVdeu6XyfEn+yRPHHiBOnp6ZQoUZKgoKy9KxUrVQLg2LGj7mPWWl59eQjp6ekMeenV8/o7laLDceHFGDPfGBNurf0YGAy8BRwDbrPWzizc1hUdV9YMdb/efig2f3VrnVL3YP7GvJxN79bVAFiy4QA7D5954G1+Rcclk5aeTulgX/eYmkyVQ11jAqJiktzH7r6uGi1rhfH+93+xdsdRRArKju3b3AN7W119TY5lAgIDadS4CQDLf/3lorUNoNwpYWrjhvU5lln/p2sgvYeHB7VqR7iPFytWDA8PD6KjjxEbm/V3y57drkHxJUuWch/7/LPZRK5aSa/ed1M7ok6B3YM4m+PCCzAZWGCMGQRstta+a60dZa3N+RMkeebj5UHFMoE80KYGE/q3BGDrgRN8F7nnrHX9vD2pEhbMU53rMLRXYwB+2XSQNdsL7j/1qyPCqHKZK0R8/FPOvS7nIyE5jd83R+Hl6cGb9zTF38cTD2N4sF1NGlYpTUxCirtXKaSYHy/3aMjuqFiGf7a2wNsil7Yt//zjfl2tWo1cy1WrVh2AbVvP/fPw++/L6dSxHU0b1qVls0Z0u6UTb4wYxs5TBuKernTp0txwY1sAJn4wgZkzphMf5/oykZKSwk+Lf2DwC/8B4O577qV8hZNhx9/fn/oNGpKWlsaIYa+QkJBAWloa06d9zMYNGwgICKBe/QaAa6bVqJFvcfnlZXl4wGPnfI9S9DhuzIu19lNjzDfAEGClMWYqp4x1sdaOLLTGOdSBST3x88neFbv870M88N4yklNzHkoUWtyPzeO653ju21V76D/h1wJtZ2avy5GYRL5ZubtAr53pxU8i+XpQG3pcXYXuLcNJSUvH38f1MRk2e437cdKIu5pQMsiXB8f/6h4bI1JQDh0+5H4dGpb7mK3Mc7GxscTHxREQGJjv9zp44ACenp4EBgYRFxfLln82s+WfzXw66xOeG/gCt/folWO9l14ZSkJCPL/+sowRw15hxLBX3LON0tPTqV6jBv0ffpRbu9+Wre6TTz/L/ffezbyvvuTbb+bh5eVFUpKrV3PAY0+4Hye98dpwTpw4zrDX3nCPjREBB4aXDCm4pkb7AsFooO55OXQ8AV9vTwL9vAjycz0uWbrhAEM+iTzj7Jm0dOue1lwswNv9n/zc33Yy/LO1Oa4Rc66KB3jTuZnr29uny7bnGqjO14otUXR8dQGDbqtP0+pl8Pb0YPW2I4z5ZiOf/7YTgBvqXU73luHM/W0nC9bsBaBf25rcd2N1KocGcyQmia/+2MXQ2WuyzUwSyYvMXgwgyyyd0/n5+btfx8XnL7zUioigzhV1uebaawkLuwxPT08SEhL4ZdnPjHr7TXbv3sWwV1+mZMlStGnXPlv9YsWLM/LdsYwbO5qpUyZhrSUm5uRj4vj4eI4dO0ZaWlq2cSr1GzTkoynTeW/Mu6xbu5rU1FQi6tTh7j596dDxJgB++Xkp3337DW3bteea1tcCMGP6VGbP/ITdu3dRomRJ2rRtx4DHnsg2M0mKPseFF2NMe2Ak8BXQyFqruannqd4TX7hflynmS49WVXi6yxUsfqUDb33xJ8PnrMux3pGYpCzTmsuWCuDeG6ozoGNtbmpSnmcnr2DKjwXzeOe2qyq7w9GFeGR0qshtR+j2es57X/n7ePL2vc04HpfMwKkrAXilZyMeuzmCvUfjmLN8B3UrleSh9rVoUKUUN726kLR0DeSVf587e9+d7Zi/vz83tmlLk6ZN6Xl7N/bt3cvIt97gxrbtsq33smnTRp4Y8DCHDh3k9h69uO2OHpQvX4EjUVEs+Wkx498bw+hRbxO5agVjxk3AwyPrKIW69erx/gcTc2xbQkICw4a+THBwMM89PwiAkW+9zpRJHxEaFkaHjjfz91+bmDFtKhvXr2filGm5LqQnRZMTx7wMwjU4d6CCS8GLOpHE2Pmb6Pb6YiyW526tR7uG5fJUd9/ReIbNXku/937Bx8uTkfc144qKJQqkXXdlPDJaseUwm/YcL5Brnov/3FqP8NBgXpy5moPRCVS7vBgDOtbmYHQCrQfN5+EJy7nuv9/y88YDtKgRyp2tz22xPLm0ndqDkpiYmGu5xMSTCzoGBuT/kVFuSpQoyQP9+gOuadF/bdqY5XxcXCyPPPgABw7s54EH+/PC4CFUr14Df39/yleowJ133cPbo8ZgjGHZz0v5Yu6Z12463fvjxrJ3zx4ef+oZQkJC2bF9Gx9PnkTp0mWYOXsurw5/jemzPqNJ02asWbOaL+d+XmD3Ls7guPBirb3aWruhsNtR1EVuO8Jvf7umQ/a5rnq+6n69cje7Dsfi6eHBXddWO++21A8vRf3KrtkHH/+49byvd67qVCjBIx1q89vmQ0xe7BpQ2bFxeTw8DJ/+sp2oE65n9qlp1r0Wzk2Ns09xFTmb0JCTM/YOHTyYa7nMc0FBQec03uVM6jVo4H69Z3fWMWbffP0VR464Bq/f3ee+HOs3b3Gle5bRDwsX5Pl9N//9F9M+nkyDBg3pftsdAPy4+AestdzUqTOlS5cGwNvbm95335NxflGery9Fg+PCi1w8+465OrYqX5b/58n7j7m+EVYOy33Z8bzqfa2r9yI2MYXPl+847+udC2Ng1P3NsVie+PB39/HwUNffzelTwjO3Rcg8L5If1aqf/MKwZUv2LSlOnnOF6CpVz/9LQn5s3er6ElGyVKlsa7WcKnPdlr17zz5jESA9PZ1XXhoCGP778qvuR1V79rjCU4UKWb8MVKwYnuW8XDoUXiRXmf/xxibkf9BppZCMuokpZyl5Zn7entzWMhxwDQSOSyqcAbD331iDptVCeHfeRv7am/2xVeZ4nEx+3q4BihaNd5H8C69chcsvd60k/cuyn3MsEx8fT+Qq17irK1teVeBtWLd2jft1ufLls5zz8HCFiuhjx4iPz/3p/ZEoV+9MQB4fac36ZAZ/rlvLvff1dU8DP1ViYlKWn5OSXI/UDOe2OaU4l8LLJcgjD7vQXlPnMhpXKQPAsk0nu609Pc5et3frqlxW0jULYtnG3Lu886Jzs4qUCHStavtxAQ3+za/LS/oz+PYGbD1wgre+yLqD9c6MBfwaVS2d5XjT6mUyzhfsQnpy6bi5cxcAvvt2fo49F7M+mU58fDyenp50vLlTvq5tz7Ia9PHoaCZ+MAGAsMsuy7LIHEDt2nXc1/ns01k5XuOffzazZnUkAPVPeQSVm4MHDzJ29DtUrFiJBx56OMu5cuVc4Wn9+qyTB9atda2xdHq4kqJP4eUSVL50AD8P70if66u7e0gylSsVwBOd6jDjqdZ4eBiOxiQx7ttN7vNX1gxl/n/bcEerytmW/q8SFsyLdzTgnfuaA7DtQAwzlm7L9v7jHryS6Om9iZ7e+6xtvTtjzMymPdHZdrfOjZenoVSQr/tPsYCTq+WWCPTJcs7L8+xh7I17mlI8wIcnP/qDpJSsU7S/zVjAr1PTCnRq4urSrlWuOE92dv1yn79K3dmXuhPHj3Ps2FH3n/R017+hxISELMdPnR4NcM+9fSlTJoTEhAQG9H/QvZJtSnIyn86cwXtj3gWgW/fbc9zX6L8vDKR+nZrUr1Mz27l5X3/Jk48PYNGC7zly5OTO7ImJiSz+YRG9e93hHufy1DP/yTZT6Ma27QgJCQFgzLsj+WDCeKKjjwGuHqFv539D/wfuIzU1FW9vb3r0Ovtn/bXhrxIbG8vgF10bNZ7q2utuAOCHhQtZtHAB1lq2bPmHiR9OyHJeLh2aW3aJqlupFKP6ukJGUkoaMQkp+Pl4utd5AdhxKIa7Ry3l0PGssx1a1gqjZS3X4lgJyanEJaYS4OtFgO/Jf05/7jzKnSOXkJhy7gu4VQ4LomXGVgNT89Hr0qJGKPMGt8nx3NLhN2X5+eahC7P0LJ2uY+PydGpakU+WbmXphgPZzm/ed4L3v/uLh9rXYuqTrYlPSnX/PazYcpjpSwtvgLH8O9zRvSv79u3NdnzypIlMnnRyqnDnLl2z7O4cHBzMmHHv079fX7Zt3ULP27sRGBhIUlIyqamux7FXtmzFswNfyHeb0tPSWbxoIYsXLQTA3z8AX18fYmJiSEtzfWZ9fHx45rmBtO/QMVv9gIAA3hk9jkcfeZBjR48ydvQoxo4eRWBgIPHx8e6eHT9/f4YOe41KlcLP2J4fFy9i8aKFdOp8C81bXJntfJWqVenV+y5mTJvK0088ip+fn3sWVt169enS9dZ8/x2Isym8XIL2H0vgnneX0ioijCZVyxBWwp/Swb6kpVt2R8Wyfmc081ftZvavO7KFjzXbj/Dg+F9oVTuMBpVLEVrcn1JBviSlprHtQAxrdxzlqxW7+PL3Xee9UWHv1tXw8DAkpaTlugv2hRTo68XrdzflSEwig6ZH5lpu4NSV7DocS58bqhMeGsSBYwl89cdOhs5eS2qaxrzIuYuocwWff/kNH038gKVLfuLggf34B/hTrVo9Onfpyi23dsvWK5IXTZs359HHn2TtmjVs37aV6OPRxMbGEhgYRIWKFWnWvAXdb78jxw0hM9WtV48vvp7PpzM/4eelS9ixfTtxcbGu6dLlK9CsxZX07NU7y9YAOYmPi+O1YUMpUaIETz/3n1zLPTdwEGXLluez2TPZs3sPZcqEcGPGInXe3t651pOiyZzt2WdRVeLOaZfmjYtcJAemnP1RgYicHz+vS3O0ssa8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoyi8iIiIiKMovIiIiIijKLyIiIiIoxhrbWG3QSRPjDH9rLX/K+x2iBRV+oyJU6jnRZykX2E3QKSI02dMHEHhRURERBxF4UVEREQcReFFnETP4kUuLH3GxBE0YFdEREQcRT0vIiIi4igKL1LojDFdjTHWGFMr4+dwY8z6jNdNjDGjC7eFIs5kjEkzxqwxxmwwxqw1xjxljCmQ3/vGmJeMMc8UxLVE8kvhRf4NegLLgB6nn7DWrrTWPnbxmyRSJCRYaxtYa+sAbYCOwIuF3CaR86bwIoXKGBMEXAX0JYfwYoy51hgzzxjjYYzZYYwpccq5LcaYMGNMiDFmjjFmRcafqy7iLYg4grX2EK51XAYYF09jzJsZn5l1xpgHwfWZNMb8YIyJNMb8aYzpknkNY8wgY8zfxphFQM1CuhURvAq7AXLJuwX4zlq72Rhz1BjTCDh6eiFrbbox5kugKzDJGNMc2GGtPWiMmQG8Y61dZoypCHwP1L6YNyHiBNbabRmPjUKBLsBxa21TY4wv8IsxZgGwG+hqrT1hjCkD/GaM+QpohOsLRkNc/3dEAqsK5UbkkqfwIoWtJzAq4/XMjJ/fy6XsLGAIMAnXL9FZGcdvBCKMMZnlihljgq21MRekxSLOlvlBaQvUM8Z0z/i5OFAd2AMMN8ZcA6QD5YAw4GpgrrU2HiAj0IgUCoUXKTTGmNLA9cAVxhgLeAIWGJdLleVANWNMCK4em6EZxz2AK621CRe4ySKOZoypAqQBh3CFmEettd+fVqYPEAI0ttamGGN2AH4Zp7W2hvwr4du87gAAAvdJREFUaMyLFKbuwMfW2krW2nBrbQVgO1A+p8LWtSjRXGAksMlaeyTj1AJgQGY5Y0yDC9tsEefJCP3vA2MzPkvfA/2NMd4Z52sYYwJx9cAcyggu1wGVMi6xFOhqjPE3xgQDnS7+XYi4qOdFClNP4LXTjs0BXjhDnVnACqDPKcceA94zxqzD9W96KfBQwTVTxLH8jTFrAG8gFZiKK/wDfAiEA5HG9cz1MK4ezenA18aYlcAa4C8Aa22kMWZWxrGdwM8X8T5EstAKuyIiIuIoemwkIiIijqLwIiIiIo6i8CIiIiKOovAiIiIijqLwIiIiIo6i8CIiZ5Sxn82GjP1v1hhjmhtjnjDGBOShbp7KiYjkh6ZKi0iujDFX4loX5FprbVLGXjc+wK9AE2tt1Fnq78hLORGR/FDPi4icyeVAlLU2CSAjhHQHygI/GmN+BDDGjDfGrMzooXk549hjOZSLzbywMaa7MWZyxuvbjDHrjTFrjTFLL+L9iYgDqedFRHJljAkClgEBwCJglrV2yek9KsaYUtbao8YYT+AH4DFr7bocysVaa4MyXncHbrbW9jHG/Am0t9buNcaUsNZGX+x7FRHnUM+LiOTKWhsLNAb64Vo+flbGxn2nu90YEwmsBuoAEfl8q1+AycaYB3Bt0CkikivtbSQiZ2StTQN+An7K6CG559TzxpjKwDNAU2vtsYxHQX6nXyfzcqe8dpex1j5kjGkO3ASsMcY0OGXjTRGRLNTzIiK5MsbUNMZUP+VQA1yb8sUAwRnHigFxwHFjTBjQ4ZTyp5YDOGiMqW2M8QC6nvI+Va21v1trhwBRQIWCvxsRKSrU8yIiZxIEjDHGlMC1K/EWXI+QegLfGmP2W2uvM8asBjYA23A9Asr0v1PLAQOBecBuYH3G9QHezAhJBteYmbUX/tZExKk0YFdEREQcRY+NRERExFEUXkRERMRRFF5ERETEURReRERExFEUXkRERMRRFF5ERETEURReRERExFEUXkRERMRR/g8i0hcDYy/L1wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "count = np.array(data.groupby(['Smoker', 'Status']).count())\n", "count = np.reshape(count, (2, 2))\n", "annots = np.array([f\"{v}\\n{v/len(data):.2%}\" for v in count.flatten()]).reshape(2,2)\n", "\n", "plt.figure(figsize=(10,8))\n", "sns.heatmap(count, annot=annots, fmt=\"\", cmap='Blues', cbar=False, square=True,\n", " xticklabels=['Alive', 'Dead'], yticklabels=['No', 'Yes'], annot_kws={\"fontsize\": 25})\n", "plt.title(\"Number and percentage of alive/dead women after 20 years, according to smoking habits\")\n", "plt.xlabel(\"Status\")\n", "plt.ylabel(\"Smoker\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is possible to see that the fraction of smokers and non smokers is quite balanced (in total, 582 smokers and 732 non smokers). As expected, there are less dead than alive people (369 versus 945).\n", "\n", "We can then compute the mortality rate for the two groups. For a population proportion $p$, confidence intervals are computed as $\\hat{p} \\pm z \\cdot \\sqrt{\\frac{\\hat{p}(1-\\hat{p})}{n}}$, where $\\hat{p}$ is the sample proportion, $n$ is the sample size and $z$ is the value derived from the standard normal distribution. For 95% confidence intervals, $z=1.96$." ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mortality rate for smokers:\t23.88% ± 3.46%\n", "Mortality rate for non smokers:\t31.42% ± 3.36%\n" ] } ], "source": [ "z = 1.96\n", "\n", "num_smokers = sum(data['Smoker'] == \"Yes\")\n", "num_dead_smokers = sum(np.logical_and(data['Smoker'] == \"Yes\", data['Status'] == \"Dead\"))\n", "rate_smokers = num_dead_smokers / num_smokers\n", "ci_smokers = z * (rate_smokers * (1 - rate_smokers) / num_smokers) ** 0.5\n", "print(f\"Mortality rate for smokers:\\t{rate_smokers:.2%} \" + u\"\\u00B1\" + f\" {ci_smokers:.2%}\")\n", "\n", "num_non_smokers = len(data) - num_smokers\n", "num_dead_non_smokers = sum(np.logical_and(data['Smoker'] == \"No\", data['Status'] == \"Dead\"))\n", "rate_non_smokers = num_dead_non_smokers / num_non_smokers\n", "ci_non_smokers = z * (rate_non_smokers * (1 - rate_non_smokers) / num_non_smokers) ** 0.5\n", "print(f\"Mortality rate for non smokers:\\t{rate_non_smokers:.2%} \" + u\"\\u00B1\" + f\" {ci_non_smokers:.2%}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Surprisingly, the mortality rate is sensibly higher for women categorized as non smokers. However, we are not taking into account an important information: the age of those people at the time of the poll. This result can be expected, for example, if the average age of polled non smokers was higher than the one of smokers.\n", "\n", "---\n", "\n", "Let's now include the age in the analysis. The following age classes are considered: 18-34 years, 35-54 years, 55-64 years, over 65 years." ] }, { "cell_type": "code", "execution_count": 54, "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", "
SmokerStatusAgeBinned age
0YesAlive21.018-34 years
1YesAlive19.318-34 years
2NoDead57.555-64 years
3NoAlive47.135-54 years
4YesAlive81.4Over 65 years
\n", "
" ], "text/plain": [ " Smoker Status Age Binned age\n", "0 Yes Alive 21.0 18-34 years\n", "1 Yes Alive 19.3 18-34 years\n", "2 No Dead 57.5 55-64 years\n", "3 No Alive 47.1 35-54 years\n", "4 Yes Alive 81.4 Over 65 years" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def bin_age(age):\n", " if age < 18:\n", " return None\n", " if age < 35:\n", " return \"18-34 years\"\n", " elif age < 55:\n", " return \"35-54 years\"\n", " elif age < 65:\n", " return \"55-64 years\"\n", " else:\n", " return \"Over 65 years\"\n", "\n", "data['Binned age'] = data['Age'].apply(bin_age)\n", "data.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, let's check that no missing are present, to ensure that no women under 18 was polled." ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of rows with missing values: 0\n" ] } ], "source": [ "print(\"Number of rows with missing values:\", data.isnull().any(axis=1).sum())" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[221],\n", " [ 6],\n", " [182],\n", " [ 7],\n", " [172],\n", " [ 19],\n", " [190],\n", " [ 39],\n", " [ 81],\n", " [ 40],\n", " [ 64],\n", " [ 51],\n", " [ 28],\n", " [165],\n", " [ 7],\n", " [ 42]])" ] }, "execution_count": 56, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.array(data.groupby(['Binned age', 'Smoker', 'Status']).count())" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Age
Binned ageSmokerStatus
18-34 yearsNoAlive221
Dead6
YesAlive182
Dead7
35-54 yearsNoAlive172
Dead19
YesAlive190
Dead39
55-64 yearsNoAlive81
Dead40
YesAlive64
Dead51
Over 65 yearsNoAlive28
Dead165
YesAlive7
Dead42
\n", "
" ], "text/plain": [ " Age\n", "Binned age Smoker Status \n", "18-34 years No Alive 221\n", " Dead 6\n", " Yes Alive 182\n", " Dead 7\n", "35-54 years No Alive 172\n", " Dead 19\n", " Yes Alive 190\n", " Dead 39\n", "55-64 years No Alive 81\n", " Dead 40\n", " Yes Alive 64\n", " Dead 51\n", "Over 65 years No Alive 28\n", " Dead 165\n", " Yes Alive 7\n", " Dead 42" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.groupby(['Binned age', 'Smoker', 'Status']).count()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Barchart with one bar for each age, divided in four parts(smoker or not, dead or alive)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }