{ "cells": [ { "cell_type": "markdown", "id": "bd13e811", "metadata": {}, "source": [ "# GG274 Homework 8: Hypothesis Testing\n", "\n", "Reisdential instability is one component of the [Ontario Marginalization Index](https://www.publichealthontario.ca/-/media/Documents/O/2017/on-marg-technical.pdf?la=en≻lang=en&hash=EED54DF437EDEDA2DFE1A00A4B14A50A) that includes indicators of types and density of residential accommodations, and certain family structure characteristics, such as living alone and dwelling owndership. ([see OCHPP](https://www.ontariohealthprofiles.ca/canmargCAN.php))\n", "\n", "In this homework you will explore the following question:\n", "\n", "> **Are mental health visits different in Toronto neighbourhoods with higher \"residential instability\"?**" ] }, { "cell_type": "code", "execution_count": 22, "id": "498ebc38", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd" ] }, { "cell_type": "markdown", "id": "dd16cb9a", "metadata": {}, "source": [ "## Step 1 - Read the Neighbourhood Instability data into a `pandas` `DataFrame`\n", "\n", "a) The data is stored in `1_marg_neighb_toronto_2006_OnMarg.xls` - a Microsoft Excel file format with file extension `.xls`.\n", "\n", "Use the `pandas` function `read_excel` to read the sheet `Neighbourhood_Toronto_OnMarg` into a `pandas` `DataFrame` named `marg_neighb`. \n" ] }, { "cell_type": "code", "execution_count": 23, "id": "0e70a55c", "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", " \n", " \n", "
Neighb idNeighbourhood namePOPULATIONINSTABILITYINSTABILITY_QDEPRIVATIONDEPRIVATION_QETHNICCONCENTRATIONETHNICCONCENTRATION_QDEPENDENCYDEPENDENCY_QONMARG_COMBINED_Q
01West Humber-Clairville32252-0.663910.162032.45485-0.202132.4
12Mount Olive-Silverstone-Jamestown32127-0.108111.019553.74335-0.597512.4
23Thistletown-Beaumond Heights9928-0.313110.346041.622040.284552.8
34Rexdale-Kipling107250.186620.470441.239630.273452.8
45Elms-Old Rexdale9879-0.015020.804051.99114-0.352722.6
\n", "
" ], "text/plain": [ " Neighb id Neighbourhood name POPULATION INSTABILITY \\\n", "0 1 West Humber-Clairville 32252 -0.6639 \n", "1 2 Mount Olive-Silverstone-Jamestown 32127 -0.1081 \n", "2 3 Thistletown-Beaumond Heights 9928 -0.3131 \n", "3 4 Rexdale-Kipling 10725 0.1866 \n", "4 5 Elms-Old Rexdale 9879 -0.0150 \n", "\n", " INSTABILITY_Q DEPRIVATION DEPRIVATION_Q ETHNICCONCENTRATION \\\n", "0 1 0.1620 3 2.4548 \n", "1 1 1.0195 5 3.7433 \n", "2 1 0.3460 4 1.6220 \n", "3 2 0.4704 4 1.2396 \n", "4 2 0.8040 5 1.9911 \n", "\n", " ETHNICCONCENTRATION_Q DEPENDENCY DEPENDENCY_Q ONMARG_COMBINED_Q \n", "0 5 -0.2021 3 2.4 \n", "1 5 -0.5975 1 2.4 \n", "2 4 0.2845 5 2.8 \n", "3 3 0.2734 5 2.8 \n", "4 4 -0.3527 2 2.6 " ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "marg_neighb = pd.read_excel('1_marg_neighb_toronto_2006_OnMarg.xls',\n", " sheet_name='Neighbourhood_Toronto_OnMarg', header=1)\n", "\n", "marg_neighb.head()" ] }, { "cell_type": "markdown", "id": "a9b80a4d", "metadata": {}, "source": [ "Use `marg_neighb` to create a another `DataFrame` called `instability_df` that has three columns: `'Neighb id ', 'Neighbourhood name ', 'INSTABILITY'`." ] }, { "cell_type": "code", "execution_count": 24, "id": "260275b5", "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", "
Neighb idNeighbourhood nameINSTABILITY
01West Humber-Clairville-0.6639
12Mount Olive-Silverstone-Jamestown-0.1081
23Thistletown-Beaumond Heights-0.3131
34Rexdale-Kipling0.1866
45Elms-Old Rexdale-0.0150
\n", "
" ], "text/plain": [ " Neighb id Neighbourhood name INSTABILITY\n", "0 1 West Humber-Clairville -0.6639\n", "1 2 Mount Olive-Silverstone-Jamestown -0.1081\n", "2 3 Thistletown-Beaumond Heights -0.3131\n", "3 4 Rexdale-Kipling 0.1866\n", "4 5 Elms-Old Rexdale -0.0150" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "instability_df = marg_neighb[marg_neighb.columns[[0, 1, 3]]]\n", "\n", "instability_df.head()" ] }, { "cell_type": "markdown", "id": "ba21291a", "metadata": {}, "source": [ "b) Rename the column names of `instability_df` using the following table. The DataFrame with the new column names should be called `instability_df` (i.e., don't change the name of the DataFrame).\n", "\n", "Original column name | New column name\n", "----|----\n", "Neighb id | Neighbid\n", "INSTABILITY | INSTABILITY\n", "Neighbourhood name | name\n" ] }, { "cell_type": "code", "execution_count": 25, "id": "911d50e5", "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", "
NeighbidnameINSTABILITY
01West Humber-Clairville-0.6639
12Mount Olive-Silverstone-Jamestown-0.1081
23Thistletown-Beaumond Heights-0.3131
34Rexdale-Kipling0.1866
45Elms-Old Rexdale-0.0150
\n", "
" ], "text/plain": [ " Neighbid name INSTABILITY\n", "0 1 West Humber-Clairville -0.6639\n", "1 2 Mount Olive-Silverstone-Jamestown -0.1081\n", "2 3 Thistletown-Beaumond Heights -0.3131\n", "3 4 Rexdale-Kipling 0.1866\n", "4 5 Elms-Old Rexdale -0.0150" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "colnames = {'Neighb id ': 'Neighbid',\n", " 'INSTABILITY' : 'INSTABILITY',\n", " 'Neighbourhood name ': 'name'}\n", "\n", "instability_df = instability_df.copy()\n", "\n", "instability_df.rename(columns = colnames, inplace=True)\n", "\n", "instability_df.head()" ] }, { "cell_type": "markdown", "id": "cf32e981", "metadata": {}, "source": [ "## Step 2 - Read the mental health visit data into a `pandas` `DataFrame`.\n", "\n", "a) In this step you will read in data on rates of mental health visits stored in `2_ahd_neighb_db_ast_hbp_mhv_copd_2012.xls` into a `pandas` `DataFrame` named `mentalhealth_neighb`.\n" ] }, { "cell_type": "code", "execution_count": 26, "id": "0243bcf1", "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", " \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", "
Unnamed: 0Unnamed: 1MaleFemaleBoth sexesMale.1Female.1Both sexes.1Male.2Female.2...Female.12Both sexes.12Rate Ratio**, Both sexes.4H/ L/ NS, Both sexes.4(95% CI) LL, Male.4(95% CI) UL, Male.4(95% CI) LL, Female.4(95% CI) UL, Female.4(95% CI) LL, Both sexes.4(95% CI) UL, Both sexes.4
01West Humber-Clairville938116821061391514046279616.68.2...7.27.20.85L6.18.36.28.36.48.0
12Mount Olive-Silverstone-Jamestown866113019961225613082253387.08.6...8.58.20.96NS6.69.47.210.07.39.3
23Thistletown-Beaumond Heights2754106854124445385776.49.2...8.17.90.93NS5.99.96.410.06.79.3
34Rexdale-Kipling3284537814130447086007.710.0...9.09.11.07NS7.211.67.311.07.810.6
45Elms-Old Rexdale2873966833787402878157.49.6...8.17.50.88NS4.89.26.210.56.19.1
\n", "

5 rows × 81 columns

\n", "
" ], "text/plain": [ " Unnamed: 0 Unnamed: 1 Male Female Both sexes \\\n", "0 1 West Humber-Clairville 938 1168 2106 \n", "1 2 Mount Olive-Silverstone-Jamestown 866 1130 1996 \n", "2 3 Thistletown-Beaumond Heights 275 410 685 \n", "3 4 Rexdale-Kipling 328 453 781 \n", "4 5 Elms-Old Rexdale 287 396 683 \n", "\n", " Male.1 Female.1 Both sexes.1 Male.2 Female.2 ... Female.12 \\\n", "0 13915 14046 27961 6.6 8.2 ... 7.2 \n", "1 12256 13082 25338 7.0 8.6 ... 8.5 \n", "2 4124 4453 8577 6.4 9.2 ... 8.1 \n", "3 4130 4470 8600 7.7 10.0 ... 9.0 \n", "4 3787 4028 7815 7.4 9.6 ... 8.1 \n", "\n", " Both sexes.12 Rate Ratio**, Both sexes.4 H/ L/ NS, Both sexes.4 \\\n", "0 7.2 0.85 L \n", "1 8.2 0.96 NS \n", "2 7.9 0.93 NS \n", "3 9.1 1.07 NS \n", "4 7.5 0.88 NS \n", "\n", " (95% CI) LL, Male.4 (95% CI) UL, Male.4 (95% CI) LL, Female.4 \\\n", "0 6.1 8.3 6.2 \n", "1 6.6 9.4 7.2 \n", "2 5.9 9.9 6.4 \n", "3 7.2 11.6 7.3 \n", "4 4.8 9.2 6.2 \n", "\n", " (95% CI) UL, Female.4 (95% CI) LL, Both sexes.4 (95% CI) UL, Both sexes.4 \n", "0 8.3 6.4 8.0 \n", "1 10.0 7.3 9.3 \n", "2 10.0 6.7 9.3 \n", "3 11.0 7.8 10.6 \n", "4 10.5 6.1 9.1 \n", "\n", "[5 rows x 81 columns]" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mentalhealth_neighb = pd.read_excel('2_ahd_neighb_db_ast_hbp_mhv_copd_2012.xls',\n", " sheet_name='2_MentalHealthV_2012', header=11)\n", "mentalhealth_neighb.head()" ] }, { "cell_type": "markdown", "id": "637787f9", "metadata": {}, "source": [ "b) Create a new DataFrame `mhvisitrates` by selecting the columns in `mentalhealth_neighb` that corresponds to Neighbourhood ID, Neighbourhood Name, and 'Age-Standardized rate of Mental Health Visits (2012), All Ages 20+' rename this column in `mhvisitrates` to `mhvisitrates_mf`. When you rename this column don't change the name of the DataFrame `mhvisitrates`." ] }, { "cell_type": "code", "execution_count": 27, "id": "449b31ab", "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", "
Neighbidnamemhvisitrates_mf
01West Humber-Clairville7.4
12Mount Olive-Silverstone-Jamestown7.8
23Thistletown-Beaumond Heights7.8
34Rexdale-Kipling8.9
45Elms-Old Rexdale8.5
\n", "
" ], "text/plain": [ " Neighbid name mhvisitrates_mf\n", "0 1 West Humber-Clairville 7.4\n", "1 2 Mount Olive-Silverstone-Jamestown 7.8\n", "2 3 Thistletown-Beaumond Heights 7.8\n", "3 4 Rexdale-Kipling 8.9\n", "4 5 Elms-Old Rexdale 8.5" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mhvisitrates = mentalhealth_neighb[mentalhealth_neighb.columns[[0, 1, 10]]]\n", "\n", "colnames = {'Unnamed: 0': 'Neighbid',\n", " 'Both sexes.2' : 'mhvisitrates_mf',\n", " 'Unnamed: 1' : 'name'}\n", "\n", "mhvisitrates = mhvisitrates.copy()\n", "\n", "mhvisitrates.rename(columns = colnames, inplace=True)\n", "\n", "mhvisitrates.head()" ] }, { "cell_type": "markdown", "id": "6e4619af", "metadata": {}, "source": [ "## Step 3 - Merge mental health visits and instability\n", "\n", "In this step you will merge the `mhvisitrates` with `mentalhealth_neighb`.\n", "\n", "a) Merge `mhvisitrates` with `instability_df` and name this DataFrame `mhvisitinstab`." ] }, { "cell_type": "code", "execution_count": 28, "id": "5debecca", "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", "
Neighbidnamemhvisitrates_mfINSTABILITY
01West Humber-Clairville7.4-0.6639
12Mount Olive-Silverstone-Jamestown7.8-0.1081
23Thistletown-Beaumond Heights7.8-0.3131
34Rexdale-Kipling8.90.1866
45Elms-Old Rexdale8.5-0.0150
\n", "
" ], "text/plain": [ " Neighbid name mhvisitrates_mf INSTABILITY\n", "0 1 West Humber-Clairville 7.4 -0.6639\n", "1 2 Mount Olive-Silverstone-Jamestown 7.8 -0.1081\n", "2 3 Thistletown-Beaumond Heights 7.8 -0.3131\n", "3 4 Rexdale-Kipling 8.9 0.1866\n", "4 5 Elms-Old Rexdale 8.5 -0.0150" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mhvisitinstab = mhvisitrates.merge(instability_df, on = ['Neighbid', 'name'])\n", "mhvisitinstab.head()" ] }, { "cell_type": "markdown", "id": "b748c831", "metadata": {}, "source": [ "## Step 4\n", "\n", "a) Create a new column in `mhvisitinstab` named `instab_HL` that categorizes neighbourhoods. The new columns should have two possible values:\n", "\n", "- `'High'`, if the neighbourhood's INSTABILITY value is greater than or equal to the mean\n", "- `'Low'`, if the neighbourhood's INSTABILITY value is less than the mean" ] }, { "cell_type": "code", "execution_count": 29, "id": "bcf2d441", "metadata": {}, "outputs": [], "source": [ "mean_instab = mhvisitinstab['INSTABILITY'].mean()\n", "\n", "mhvisitinstab.loc[mhvisitinstab['INSTABILITY'] >= mean_instab, 'instab_HL'] = 'High'\n", "\n", "mhvisitinstab.loc[mhvisitinstab['INSTABILITY'] < mean_instab, 'instab_HL'] = 'Low'" ] }, { "cell_type": "markdown", "id": "05edf48a", "metadata": {}, "source": [ "b) Compute the frequency distribution of `instab_HL`. Save the results in `instab_HL_frequencies`. " ] }, { "cell_type": "code", "execution_count": 30, "id": "12d2561d", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "instab_HL\n", "Low 70\n", "High 66\n", "Name: count, dtype: int64" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "instab_HL_frequencies = mhvisitinstab['instab_HL'].value_counts()\n", "\n", "instab_HL_frequencies" ] }, { "cell_type": "markdown", "id": "66a0de88", "metadata": {}, "source": [ "c) Is there evidence that Toronto has many neighbourhoods that have residential instability? Briefly explain. __(1 mark)__" ] }, { "cell_type": "markdown", "id": "10873e1a", "metadata": {}, "source": [ "> There is more than one possible way to answer the question. e.g., \n", ">\n", "> i. The data shows that about half of the neighbourhoods in Toronto, 66 out of 136, having residential instability measures that are higher than the city-wide mean. It is difficult to provide a definite answer to the question based on the result since it only provide comparison within the city.\n", ">\n", "> ii. (a bit mor. naive answer) About half of the neighbourhoods in Toronto have high residential instability.\n", ">\n", "> Any answer with a sensible reasoning based on the data is acceptable." ] }, { "cell_type": "markdown", "id": "c240f8cb", "metadata": {}, "source": [ "## Step 5 - Do neighbourhoods with high residential instability have more mental health visits compared to neighbourhoods with low residential isntability?\n", "\n", "a) Use the `DataFrame` `describe` method to compute the distribution of `mhvistrates_mf` in `mhvisitinstab` **grouped by** `instab_HL`. Store the results in `median_table`." ] }, { "cell_type": "code", "execution_count": 31, "id": "26535b3d", "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", "
countmeanstdmin25%50%75%max
instab_HL
High66.08.3621210.9373885.97.8258.459.07510.5
Low70.07.8671431.0005395.77.0257.808.7759.9
\n", "
" ], "text/plain": [ " count mean std min 25% 50% 75% max\n", "instab_HL \n", "High 66.0 8.362121 0.937388 5.9 7.825 8.45 9.075 10.5\n", "Low 70.0 7.867143 1.000539 5.7 7.025 7.80 8.775 9.9" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "median_table = mhvisitinstab.groupby('instab_HL')['mhvisitrates_mf'].describe()\n", "\n", "median_table" ] }, { "cell_type": "markdown", "id": "2cd0db75", "metadata": {}, "source": [ "Use `median_table` to compute the difference in medians between neighbourhoods with high and low instability. Store this value in `median_diff`. " ] }, { "cell_type": "code", "execution_count": 32, "id": "978a9570", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/var/folders/0j/ybsv4ncn5w50v40vdh5jjlww0000gn/T/ipykernel_11385/223760169.py:1: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n", " median_diff = median_table['50%'][0] - median_table['50%'][1]\n" ] }, { "data": { "text/plain": [ "0.6499999999999995" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "median_diff = median_table['50%'][0] - median_table['50%'][1]\n", "\n", "median_diff" ] }, { "cell_type": "markdown", "id": "5da0b797", "metadata": {}, "source": [ "## Step 6 - Set up a simulation in python to test if the medians are equal\n", "\n", "a) In this step you will write a function `random_shuffle_median` that returns a simulated value of the median difference (a simulated value of the test statistic) of mental health visit rates in neighbourhoods with high versus low residential instability assuming that there really is no difference in mental health visit rates between these types of neighbourhoods.\n", "\n", "A step-by-step explantion of a similar function was given in lecture, and you can follow this example to help guide you through this step.\n", "\n", "The function `random_shuffle_median` is started for you below. Your task is to complete the function by filling in the `...`.\n", "\n", "Try writing a meaningful docstring for `random_shuffle_median`. The `pandas` [docstring guide](https://pandas.pydata.org/docs/development/contributing_docstring.html) has some great examples and guidelines. (NB: this will not be graded)" ] }, { "cell_type": "code", "execution_count": 33, "id": "81c53baa", "metadata": {}, "outputs": [], "source": [ "# def random_shuffle_median():\n", "# \"\"\"\n", "# Put your docstring here (optional)\n", "# \"\"\"\n", "\n", "# # shuffle the column of mhvisitinstab that corresponds to high/low instability\n", "\n", "# instab_HL_shuffle = mhvisitinstab[...].sample(frac=1, replace=False).reset_index(drop = True)\n", " \n", "# # calculate the median visit rate for high and low instability neighbourhoods\n", "\n", "# visitrate_low_shuffle = mhvisitinstab.loc[instab_HL_shuffle == ..., ...].median()\n", " \n", "# visitrate_high_shuffle = mhvisitinstab.loc[instab_HL_shuffle == ..., ...].median()\n", " \n", "# shuffled_diff = ... - ...\n", " \n", "# return shuffled_diff\n" ] }, { "cell_type": "code", "execution_count": 34, "id": "20b977f6", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def random_shuffle_median():\n", " \"\"\"\n", " randomly shuffles the column of mhvisitinstab that corresponds to \n", "\n", " high/low instability and calculates the difference in median visit rates \n", " \n", " for high and low instability neighbourhoods\n", " \"\"\"\n", "\n", " instab_HL_shuffle = mhvisitinstab['instab_HL'].sample(frac=1, replace=False).reset_index(drop = True)\n", " visitrate_low_shuffle = mhvisitinstab.loc[instab_HL_shuffle == 'Low', 'mhvisitrates_mf'].median()\n", " visitrate_high_shuffle = mhvisitinstab.loc[instab_HL_shuffle == 'High', 'mhvisitrates_mf'].median()\n", " shuffled_diff = visitrate_high_shuffle - visitrate_low_shuffle \n", " return shuffled_diff\n", "\n", "random_shuffle_median()\n" ] }, { "cell_type": "code", "execution_count": 35, "id": "4b22ef64", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " randomly shuffles the column of mhvisitinstab that corresponds to \n", "\n", " high/low instability and calculates the difference in median visit rates \n", " \n", " for high and low instability neighbourhoods\n", " \n" ] } ], "source": [ "print(random_shuffle_median.__doc__)" ] }, { "cell_type": "markdown", "id": "2cd897a2", "metadata": {}, "source": [ "b) Explain the purpose of \n", "\n", "```python\n", "mhvisitinstab[...].sample(frac=1, replace=False).reset_index(drop = True)\n", "```\n", "\n", "in 1-2 sentences." ] }, { "cell_type": "markdown", "id": "1727b2fb", "metadata": {}, "source": [ "> The line randomly samples 100% (`frac=1`) of the selected column (`'instab_HL'`) without replacements (`replace=False`) and removes the old ordering (`reset_index(drop=True)`). It mimicks random asingment of the high vs low labels to each neighbourhood." ] }, { "cell_type": "markdown", "id": "3b012943", "metadata": {}, "source": [ "## Step 7 - Compute the distribution of simulated values of the median difference assuming the null hypothesis is true\n", "\n" ] }, { "cell_type": "markdown", "id": "976af74e", "metadata": {}, "source": [ "We will use your student number to generate data for this homework. Complete the assignment statement below by typing your student number as an `int`. In other words assign your student number as an integer to the variable `student_number`.\n" ] }, { "cell_type": "code", "execution_count": 36, "id": "f50aed80", "metadata": {}, "outputs": [], "source": [ "# # Replace the ... with your student number\n", "student_number = 12345\n", "\n", "# # This checks that you correctly typed in your student_number as an int.\n", "# # Make sure there's no error when you run this cell!\n", "# assert type(student_number) == int" ] }, { "cell_type": "markdown", "id": "32fb4d54", "metadata": {}, "source": [ "a) Write a function called `shuffled_diffs` that returns a list. The function should use a `for` loop that iterates the function `random_shuffle_median` an arbitrary number of times. The number of times that the `for` loop iterates should be controlled by a function parameter named `number_of_shuffles`. " ] }, { "cell_type": "code", "execution_count": 37, "id": "9c843e7c", "metadata": {}, "outputs": [], "source": [ "def shuffled_diffs(number_of_shuffles):\n", " shuffled_diffs = []\n", " for _ in range(number_of_shuffles):\n", " shuffled_diffs.append(random_shuffle_median())\n", " return shuffled_diffs" ] }, { "cell_type": "markdown", "id": "f06d8325", "metadata": {}, "source": [ "b) Use `shuffled_diffs` to compute 10000 simulated median differences between high and low instability neighbourhoods assuming that there is no difference in median mental health visit rates between high and low instability neighbourhoods. Store the values in `shuffled_diffs_10000`." ] }, { "cell_type": "code", "execution_count": 38, "id": "2edcc0f2", "metadata": {}, "outputs": [], "source": [ "np.random.seed(student_number)\n", "\n", "shuffled_diffs_10000 = shuffled_diffs(10000)" ] }, { "cell_type": "markdown", "id": "c2ec390e", "metadata": {}, "source": [ "c) Plot the distribution of the 10,000 simulated values stored in `shuffled_diffs_10000` using a `matplotlib` histogram. Name the plot `nullhypothesis_distribution_plot`. Label the horizontal axis as `'Difference in median visit rates for high and low instability neighbourhoods'` and the vertical axis as `'Frequency'`." ] }, { "cell_type": "code", "execution_count": 39, "id": "add94761", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Frequency')" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGwCAYAAABIC3rIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABO+klEQVR4nO3deXxM5/4H8M9kGyKZSSKyVSTWSIgoLqbUUiEhXEpvbZWo0KuNfSlp1drbqK3VVrlVEr2XKi2lQoglFLGFWIIgjYYrE2syYoksz++PvnJ+RhbJSDIT5/N+veZV55xnnvmeJ2fOfHrOmTMKIYQAERERkYyZGbsAIiIiImNjICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItmzMHYB1UFBQQFu3LgBW1tbKBQKY5dDREREZSCEwP379+Hm5gYzs9KPATEQlcGNGzfg7u5u7DKIiIjIANeuXUPdunVLbcNAVAa2trYA/hpQlUpl5GqIiIioLHQ6Hdzd3aXP8dIwEJVB4WkylUrFQERERFTNlOVyF15UTURERLLHQERERESyx0BEREREssdARERERLLHQERERESyx0BEREREssdARERERLLHQERERESyx0BEREREssdARERERLLHQERERESyx0BEREREssdARERERLLHQERERESyx0BEREREsmdh7AKIiIzNc3q0sUsot6vzg4xdAtFLhUeIiIiISPaMGoiWL1+OFi1aQKVSQaVSQaPRYMeOHdLyLl26QKFQ6D1Gjx6t10daWhqCgoJgbW0NJycnTJ06FXl5eXpt4uLi0KpVKyiVSjRq1AhRUVFVsXpERERUTRj1lFndunUxf/58NG7cGEIIrFmzBn379sWpU6fQrFkzAMCoUaMwd+5c6TnW1tbSv/Pz8xEUFAQXFxccPnwY6enpCA4OhqWlJT777DMAQGpqKoKCgjB69GisXbsWe/bswciRI+Hq6oqAgICqXWEiIiIySQohhDB2EU9zcHDAwoULERoaii5duqBly5b48ssvi227Y8cO9O7dGzdu3ICzszMAYMWKFZg2bRpu3boFKysrTJs2DdHR0Th37pz0vEGDBiEzMxMxMTHF9puTk4OcnBxpWqfTwd3dHVlZWVCpVBW3skRkEngNEdHLSafTQa1Wl+nz22SuIcrPz8f69evx4MEDaDQaaf7atWvh6OiI5s2bIzw8HA8fPpSWxcfHw9fXVwpDABAQEACdToekpCSpjb+/v95rBQQEID4+vsRaIiIioFarpYe7u3tFrSYRERGZIKN/y+zs2bPQaDR4/PgxbGxssHnzZvj4+AAAhgwZAg8PD7i5ueHMmTOYNm0akpOTsWnTJgCAVqvVC0MApGmtVltqG51Oh0ePHqFmzZpFagoPD8ekSZOk6cIjRERERPRyMnog8vLyQmJiIrKysvDzzz8jJCQE+/fvh4+PD9577z2pna+vL1xdXdGtWzekpKSgYcOGlVaTUqmEUqmstP6JiIjItBj9lJmVlRUaNWqE1q1bIyIiAn5+fli6dGmxbdu1awcAuHLlCgDAxcUFGRkZem0Kp11cXEpto1Kpij06RERERPJj9ED0rIKCAr0Lmp+WmJgIAHB1dQUAaDQanD17Fjdv3pTaxMbGQqVSSafdNBoN9uzZo9dPbGys3nVKREREJG9GPWUWHh6Onj17ol69erh//z7WrVuHuLg47Ny5EykpKVi3bh169eqF2rVr48yZM5g4cSI6deqEFi1aAAB69OgBHx8fDBs2DAsWLIBWq8WMGTMQFhYmnfIaPXo0vvnmG3z44YcYMWIE9u7diw0bNiA6uvp9q4SIiIgqh1ED0c2bNxEcHIz09HSo1Wq0aNECO3fuRPfu3XHt2jXs3r0bX375JR48eAB3d3cMGDAAM2bMkJ5vbm6Obdu24f3334dGo0GtWrUQEhKid9+i+vXrIzo6GhMnTsTSpUtRt25dfP/997wHEREREUlM7j5Epqg89zEgouqH9yEiejlVy/sQERERERkLAxERERHJHgMRERERyR4DEREREckeAxERERHJHgMRERERyR4DEREREckeAxERERHJHgMRERERyR4DEREREckeAxERERHJHgMRERERyR4DEREREckeAxERERHJHgMRERERyR4DEREREckeAxERERHJHgMRERERyR4DEREREckeAxERERHJHgMRERERyR4DEREREckeAxERERHJHgMRERERyR4DEREREckeAxERERHJHgMRERERyR4DEREREckeAxERERHJHgMRERERyR4DEREREckeAxERERHJHgMRERERyR4DEREREckeAxERERHJHgMRERERyR4DEREREckeAxERERHJnlED0fLly9GiRQuoVCqoVCpoNBrs2LFDWv748WOEhYWhdu3asLGxwYABA5CRkaHXR1paGoKCgmBtbQ0nJydMnToVeXl5em3i4uLQqlUrKJVKNGrUCFFRUVWxekRERFRNGDUQ1a1bF/Pnz0dCQgJOnDiBN954A3379kVSUhIAYOLEifjtt9+wceNG7N+/Hzdu3ED//v2l5+fn5yMoKAhPnjzB4cOHsWbNGkRFRWHmzJlSm9TUVAQFBaFr165ITEzEhAkTMHLkSOzcubPK15eIiIhMk0IIIYxdxNMcHBywcOFCvPXWW6hTpw7WrVuHt956CwBw8eJFeHt7Iz4+Hu3bt8eOHTvQu3dv3LhxA87OzgCAFStWYNq0abh16xasrKwwbdo0REdH49y5c9JrDBo0CJmZmYiJiSlTTTqdDmq1GllZWVCpVBW/0kRkVJ7To41dQrldnR9k7BKITF55Pr9N5hqi/Px8rF+/Hg8ePIBGo0FCQgJyc3Ph7+8vtWnatCnq1auH+Ph4AEB8fDx8fX2lMAQAAQEB0Ol00lGm+Ph4vT4K2xT2UZycnBzodDq9BxEREb28jB6Izp49CxsbGyiVSowePRqbN2+Gj48PtFotrKysYGdnp9fe2dkZWq0WAKDVavXCUOHywmWltdHpdHj06FGxNUVERECtVksPd3f3ilhVIiIiMlFGD0ReXl5ITEzE0aNH8f777yMkJATnz583ak3h4eHIysqSHteuXTNqPURERFS5LIxdgJWVFRo1agQAaN26NY4fP46lS5di4MCBePLkCTIzM/WOEmVkZMDFxQUA4OLigmPHjun1V/gttKfbPPvNtIyMDKhUKtSsWbPYmpRKJZRKZYWsHxEREZk+ox8helZBQQFycnLQunVrWFpaYs+ePdKy5ORkpKWlQaPRAAA0Gg3Onj2LmzdvSm1iY2OhUqng4+MjtXm6j8I2hX0QERERGfUIUXh4OHr27Il69erh/v37WLduHeLi4rBz506o1WqEhoZi0qRJcHBwgEqlwtixY6HRaNC+fXsAQI8ePeDj44Nhw4ZhwYIF0Gq1mDFjBsLCwqQjPKNHj8Y333yDDz/8ECNGjMDevXuxYcMGREdXv2+VEBERUeUwaiC6efMmgoODkZ6eDrVajRYtWmDnzp3o3r07AOCLL76AmZkZBgwYgJycHAQEBODbb7+Vnm9ubo5t27bh/fffh0ajQa1atRASEoK5c+dKberXr4/o6GhMnDgRS5cuRd26dfH9998jICCgyteXiIiITJPJ3YfIFPE+REQvN96HiOjlVC3vQ0RERERkLAxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7Rg1EERER+Nvf/gZbW1s4OTmhX79+SE5O1mvTpUsXKBQKvcfo0aP12qSlpSEoKAjW1tZwcnLC1KlTkZeXp9cmLi4OrVq1glKpRKNGjRAVFVXZq0dERETVhFED0f79+xEWFoYjR44gNjYWubm56NGjBx48eKDXbtSoUUhPT5ceCxYskJbl5+cjKCgIT548weHDh7FmzRpERUVh5syZUpvU1FQEBQWha9euSExMxIQJEzBy5Ejs3LmzytaViIiITJeFMV88JiZGbzoqKgpOTk5ISEhAp06dpPnW1tZwcXEpto9du3bh/Pnz2L17N5ydndGyZUvMmzcP06ZNw+zZs2FlZYUVK1agfv36WLx4MQDA29sbBw8exBdffIGAgIAifebk5CAnJ0ea1ul0FbG6REREZKJM6hqirKwsAICDg4Pe/LVr18LR0RHNmzdHeHg4Hj58KC2Lj4+Hr68vnJ2dpXkBAQHQ6XRISkqS2vj7++v1GRAQgPj4+GLriIiIgFqtlh7u7u4Vsn5ERERkmox6hOhpBQUFmDBhAjp06IDmzZtL84cMGQIPDw+4ubnhzJkzmDZtGpKTk7Fp0yYAgFar1QtDAKRprVZbahudTodHjx6hZs2aesvCw8MxadIkaVqn0zEUERERvcRMJhCFhYXh3LlzOHjwoN789957T/q3r68vXF1d0a1bN6SkpKBhw4aVUotSqYRSqayUvomIiMj0mMQpszFjxmDbtm3Yt28f6tatW2rbdu3aAQCuXLkCAHBxcUFGRoZem8LpwuuOSmqjUqmKHB0iIiIi+TFqIBJCYMyYMdi8eTP27t2L+vXrP/c5iYmJAABXV1cAgEajwdmzZ3Hz5k2pTWxsLFQqFXx8fKQ2e/bs0esnNjYWGo2mgtaEiIiIqjOjBqKwsDD897//xbp162BrawutVgutVotHjx4BAFJSUjBv3jwkJCTg6tWr2Lp1K4KDg9GpUye0aNECANCjRw/4+Phg2LBhOH36NHbu3IkZM2YgLCxMOu01evRo/PHHH/jwww9x8eJFfPvtt9iwYQMmTpxotHUnIiIi02HUQLR8+XJkZWWhS5cucHV1lR4//fQTAMDKygq7d+9Gjx490LRpU0yePBkDBgzAb7/9JvVhbm6Obdu2wdzcHBqNBu+88w6Cg4Mxd+5cqU39+vURHR2N2NhY+Pn5YfHixfj++++L/co9ERERyY9CCCGMXYSp0+l0UKvVyMrKgkqlMnY5RFTBPKdHG7uEcrs6P8jYJRCZvPJ8fpvERdVERERExsRARERERLLHQERERESyx0BEREREssdARERERLLHQERERESyx0BEREREssdARERERLLHQERERESyx0BEREREssdARERERLLHQERERESyx0BEREREssdARERERLLHQERERESyx0BEREREssdARERERLLHQERERESyx0BEREREssdARERERLLHQERERESyx0BEREREssdARERERLJnUCD6448/KroOIiIiIqMxKBA1atQIXbt2xX//+188fvy4omsiIiIiqlIGBaKTJ0+iRYsWmDRpElxcXPDPf/4Tx44dq+jaiIiIiKqEQYGoZcuWWLp0KW7cuIHVq1cjPT0dHTt2RPPmzbFkyRLcunWrouskIiIiqjQvdFG1hYUF+vfvj40bN+Lzzz/HlStXMGXKFLi7uyM4OBjp6ekVVScRERFRpXmhQHTixAl88MEHcHV1xZIlSzBlyhSkpKQgNjYWN27cQN++fSuqTiIiIqJKY2HIk5YsWYLIyEgkJyejV69e+OGHH9CrVy+Ymf2Vr+rXr4+oqCh4enpWZK1ERERElcKgQLR8+XKMGDECw4cPh6ura7FtnJycsGrVqhcqjoiIiKgqGBSILl++/Nw2VlZWCAkJMaR7IiIioipl0DVEkZGR2LhxY5H5GzduxJo1a164KCIiIqKqZFAgioiIgKOjY5H5Tk5O+Oyzz164KCIiIqKqZFAgSktLQ/369YvM9/DwQFpa2gsXRURERFSVDApETk5OOHPmTJH5p0+fRu3atV+4KCIiIqKqZFAgGjx4MMaNG4d9+/YhPz8f+fn52Lt3L8aPH49BgwZVdI1ERERElcqgb5nNmzcPV69eRbdu3WBh8VcXBQUFCA4O5jVEREREVO0YdITIysoKP/30Ey5evIi1a9di06ZNSElJwerVq2FlZVXmfiIiIvC3v/0Ntra2cHJyQr9+/ZCcnKzX5vHjxwgLC0Pt2rVhY2ODAQMGICMjQ69NWloagoKCYG1tDScnJ0ydOhV5eXl6beLi4tCqVSsolUo0atQIUVFRhqw6ERERvYRe6Kc7mjRpgn/84x/o3bs3PDw8yv38/fv3IywsDEeOHEFsbCxyc3PRo0cPPHjwQGozceJE/Pbbb9i4cSP279+PGzduoH///tLy/Px8BAUF4cmTJzh8+DDWrFmDqKgozJw5U2qTmpqKoKAgdO3aFYmJiZgwYQJGjhyJnTt3vsjqExER0UtCIYQQ5X1Sfn4+oqKisGfPHty8eRMFBQV6y/fu3WtQMbdu3YKTkxP279+PTp06ISsrC3Xq1MG6devw1ltvAQAuXrwIb29vxMfHo3379tixYwd69+6NGzduwNnZGQCwYsUKTJs2Dbdu3YKVlRWmTZuG6OhonDt3TnqtQYMGITMzEzExMc+tS6fTQa1WIysrCyqVyqB1IyLT5Tk92tgllNvV+UHGLoHI5JXn89ugI0Tjx4/H+PHjkZ+fj+bNm8PPz0/vYaisrCwAgIODAwAgISEBubm58Pf3l9o0bdoU9erVQ3x8PAAgPj4evr6+UhgCgICAAOh0OiQlJUltnu6jsE1hH8/KycmBTqfTexAREdHLy6CLqtevX48NGzagV69eFVZIQUEBJkyYgA4dOqB58+YAAK1WCysrK9jZ2em1dXZ2hlarldo8HYYKlxcuK62NTqfDo0ePULNmTb1lERERmDNnToWtGxEREZk2gy+qbtSoUYUWEhYWhnPnzmH9+vUV2q8hwsPDkZWVJT2uXbtm7JKIiIioEhkUiCZPnoylS5fCgMuPijVmzBhs27YN+/btQ926daX5Li4uePLkCTIzM/XaZ2RkwMXFRWrz7LfOCqef10alUhU5OgQASqUSKpVK70FEREQvL4NOmR08eBD79u3Djh070KxZM1haWuot37RpU5n6EUJg7Nix2Lx5M+Li4or8HEjr1q1haWmJPXv2YMCAAQCA5ORkpKWlQaPRAAA0Gg3+9a9/4ebNm3BycgIAxMbGQqVSwcfHR2qzfft2vb5jY2OlPoiIiEjeDApEdnZ2ePPNN1/4xcPCwrBu3Tps2bIFtra20jU/arUaNWvWhFqtRmhoKCZNmgQHBweoVCqMHTsWGo0G7du3BwD06NEDPj4+GDZsGBYsWACtVosZM2YgLCwMSqUSADB69Gh88803+PDDDzFixAjs3bsXGzZsQHR09ftmCREREVU8g752X2EvrlAUOz8yMhLDhw8H8NeNGSdPnowff/wROTk5CAgIwLfffiudDgOAP//8E++//z7i4uJQq1YthISEYP78+dJdtIG/bsw4ceJEnD9/HnXr1sUnn3wivcbz8Gv3RC83fu2e6OVUns9vgwNRXl4e4uLikJKSgiFDhsDW1hY3btyASqWCjY2NQYWbKgYiopcbAxHRy6k8n98GnTL7888/ERgYiLS0NOTk5KB79+6wtbXF559/jpycHKxYscKgwomIiIiMweAbM7Zp0wb37t3T+5bWm2++iT179lRYcURERERVwaAjRL///jsOHz5c5IdcPT098b///a9CCiMiIiKqKgYdISooKEB+fn6R+devX4etre0LF0VERERUlQwKRD169MCXX34pTSsUCmRnZ2PWrFkV+nMeRERERFXBoFNmixcvRkBAAHx8fPD48WMMGTIEly9fhqOjI3788ceKrpGIiIioUhkUiOrWrYvTp09j/fr1OHPmDLKzsxEaGoqhQ4cW+1MYRERERKbMoEAEABYWFnjnnXcqshYiIiIiozAoEP3www+lLg8ODjaoGCIiKpvqeDNJgDeUJNNlUCAaP3683nRubi4ePnwIKysrWFtbMxARERFRtWLQt8zu3bun98jOzkZycjI6duzIi6qJiIio2jEoEBWncePGmD9/fpGjR0RERESmrsICEfDXhdY3btyoyC6JiIiIKp1B1xBt3bpVb1oIgfT0dHzzzTfo0KFDhRRGREREVFUMCkT9+vXTm1YoFKhTpw7eeOMNLF68uCLqIiIiIqoyBgWigoKCiq6DiIiIyGgq9BoiIiIiourIoCNEkyZNKnPbJUuWGPISRERERFXGoEB06tQpnDp1Crm5ufDy8gIAXLp0Cebm5mjVqpXUTqFQVEyVRERERJXIoEDUp08f2NraYs2aNbC3twfw180a3333Xbz++uuYPHlyhRZJREREVJkMuoZo8eLFiIiIkMIQANjb2+PTTz/lt8yIiIio2jEoEOl0Oty6davI/Fu3buH+/fsvXBQRERFRVTIoEL355pt49913sWnTJly/fh3Xr1/HL7/8gtDQUPTv37+iayQiIiKqVAZdQ7RixQpMmTIFQ4YMQW5u7l8dWVggNDQUCxcurNACieTMc3q0sUsot6vzg4xdAhFRuRkUiKytrfHtt99i4cKFSElJAQA0bNgQtWrVqtDiiIiIiKrCC92YMT09Henp6WjcuDFq1aoFIURF1UVERERUZQwKRHfu3EG3bt3QpEkT9OrVC+np6QCA0NBQfuWeiIiIqh2DAtHEiRNhaWmJtLQ0WFtbS/MHDhyImJiYCiuOiIiIqCoYdA3Rrl27sHPnTtStW1dvfuPGjfHnn39WSGFEREREVcWgI0QPHjzQOzJU6O7du1AqlS9cFBEREVFVMigQvf766/jhhx+kaYVCgYKCAixYsABdu3atsOKIiIiIqoJBp8wWLFiAbt264cSJE3jy5Ak+/PBDJCUl4e7duzh06FBF10hERERUqQw6QtS8eXNcunQJHTt2RN++ffHgwQP0798fp06dQsOGDSu6RiIiIqJKVe4jRLm5uQgMDMSKFSvw8ccfV0ZNRERERFWq3IHI0tISZ86cqYxaiOglUB1/boSIyKBTZu+88w5WrVpV0bUQERERGYVBF1Xn5eVh9erV2L17N1q3bl3kN8yWLFlSIcURERERVYVyBaI//vgDnp6eOHfuHFq1agUAuHTpkl4bhUJRcdURERERVYFyBaLGjRsjPT0d+/btA/DXT3V89dVXcHZ2rpTiiIiIiKpCua4hevbX7Hfs2IEHDx4Y/OIHDhxAnz594ObmBoVCgV9//VVv+fDhw6FQKPQegYGBem3u3r2LoUOHQqVSwc7ODqGhocjOztZrc+bMGbz++uuoUaMG3N3dsWDBAoNrJiIiopePQRdVF3o2IJXXgwcP4Ofnh2XLlpXYJjAwEOnp6dLjxx9/1Fs+dOhQJCUlITY2Ftu2bcOBAwfw3nvvSct1Oh169OgBDw8PJCQkYOHChZg9eza+++67F6qdiIiIXh7lOmVWeJTm2XmG6tmzJ3r27FlqG6VSCRcXl2KXXbhwATExMTh+/DjatGkDAPj666/Rq1cvLFq0CG5ubli7di2ePHmC1atXw8rKCs2aNUNiYiKWLFmiF5yelpOTg5ycHGlap9MZuIZERERUHZQrEAkhMHz4cOkHXB8/fozRo0cX+ZbZpk2bKqzAuLg4ODk5wd7eHm+88QY+/fRT1K5dGwAQHx8POzs7KQwBgL+/P8zMzHD06FG8+eabiI+PR6dOnWBlZSW1CQgIwOeff4579+7B3t6+yGtGRERgzpw5FbYOREREZNrKFYhCQkL0pt95550KLeZZgYGB6N+/P+rXr4+UlBR89NFH6NmzJ+Lj42Fubg6tVgsnJye951hYWMDBwQFarRYAoNVqUb9+fb02hReBa7XaYgNReHg4Jk2aJE3rdDq4u7tX9OoRERGRiShXIIqMjKysOoo1aNAg6d++vr5o0aIFGjZsiLi4OHTr1q3SXlepVEpHwYiIiOjl90IXVVe1Bg0awNHREVeuXAEAuLi44ObNm3pt8vLycPfuXem6IxcXF2RkZOi1KZwu6dokIiIikpdqFYiuX7+OO3fuwNXVFQCg0WiQmZmJhIQEqc3evXtRUFCAdu3aSW0OHDiA3NxcqU1sbCy8vLyKPV1GRERE8mPUQJSdnY3ExEQkJiYCAFJTU5GYmIi0tDRkZ2dj6tSpOHLkCK5evYo9e/agb9++aNSoEQICAgAA3t7eCAwMxKhRo3Ds2DEcOnQIY8aMwaBBg+Dm5gYAGDJkCKysrBAaGoqkpCT89NNPWLp0qd41QkRERCRvRg1EJ06cwKuvvopXX30VADBp0iS8+uqrmDlzJszNzXHmzBn8/e9/R5MmTRAaGorWrVvj999/17u+Z+3atWjatCm6deuGXr16oWPHjnr3GFKr1di1axdSU1PRunVrTJ48GTNnzizxK/dEREQkPwrxondXlAGdTge1Wo2srCyoVCpjl0My4jk92tglEFWoq/ODjF0CyUh5Pr+r1TVERERERJWBgYiIiIhkj4GIiIiIZI+BiIiIiGSPgYiIiIhkj4GIiIiIZI+BiIiIiGSPgYiIiIhkj4GIiIiIZI+BiIiIiGSPgYiIiIhkj4GIiIiIZI+BiIiIiGSPgYiIiIhkj4GIiIiIZI+BiIiIiGSPgYiIiIhkj4GIiIiIZI+BiIiIiGSPgYiIiIhkj4GIiIiIZI+BiIiIiGSPgYiIiIhkj4GIiIiIZI+BiIiIiGSPgYiIiIhkj4GIiIiIZI+BiIiIiGSPgYiIiIhkj4GIiIiIZI+BiIiIiGSPgYiIiIhkj4GIiIiIZI+BiIiIiGSPgYiIiIhkj4GIiIiIZI+BiIiIiGSPgYiIiIhkz6iB6MCBA+jTpw/c3NygUCjw66+/6i0XQmDmzJlwdXVFzZo14e/vj8uXL+u1uXv3LoYOHQqVSgU7OzuEhoYiOztbr82ZM2fw+uuvo0aNGnB3d8eCBQsqe9WIiIioGjFqIHrw4AH8/PywbNmyYpcvWLAAX331FVasWIGjR4+iVq1aCAgIwOPHj6U2Q4cORVJSEmJjY7Ft2zYcOHAA7733nrRcp9OhR48e8PDwQEJCAhYuXIjZs2fju+++q/T1IyIioupBIYQQxi4CABQKBTZv3ox+/foB+OvokJubGyZPnowpU6YAALKysuDs7IyoqCgMGjQIFy5cgI+PD44fP442bdoAAGJiYtCrVy9cv34dbm5uWL58OT7++GNotVpYWVkBAKZPn45ff/0VFy9eLFNtOp0OarUaWVlZUKlUFb/yRCXwnB5t7BKIKtTV+UHGLoFkpDyf3yZ7DVFqaiq0Wi38/f2leWq1Gu3atUN8fDwAID4+HnZ2dlIYAgB/f3+YmZnh6NGjUptOnTpJYQgAAgICkJycjHv37hX72jk5OdDpdHoPIiIienmZbCDSarUAAGdnZ735zs7O0jKtVgsnJye95RYWFnBwcNBrU1wfT7/GsyIiIqBWq6WHu7v7i68QERERmSyTDUTGFB4ejqysLOlx7do1Y5dERERElchkA5GLiwsAICMjQ29+RkaGtMzFxQU3b97UW56Xl4e7d+/qtSmuj6df41lKpRIqlUrvQURERC8vkw1E9evXh4uLC/bs2SPN0+l0OHr0KDQaDQBAo9EgMzMTCQkJUpu9e/eioKAA7dq1k9ocOHAAubm5UpvY2Fh4eXnB3t6+itaGiIiITJlRA1F2djYSExORmJgI4K8LqRMTE5GWlgaFQoEJEybg008/xdatW3H27FkEBwfDzc1N+iaat7c3AgMDMWrUKBw7dgyHDh3CmDFjMGjQILi5uQEAhgwZAisrK4SGhiIpKQk//fQTli5dikmTJhlprYmIiMjUWBjzxU+cOIGuXbtK04UhJSQkBFFRUfjwww/x4MEDvPfee8jMzETHjh0RExODGjVqSM9Zu3YtxowZg27dusHMzAwDBgzAV199JS1Xq9XYtWsXwsLC0Lp1azg6OmLmzJl69yoiIiIieTOZ+xCZMt6HiIyF9yGilw3vQ0RV6aW4DxERERFRVWEgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItmzMHYBREQkH57To41dQrldnR9k7BKoCpj0EaLZs2dDoVDoPZo2bSotf/z4McLCwlC7dm3Y2NhgwIAByMjI0OsjLS0NQUFBsLa2hpOTE6ZOnYq8vLyqXhUiIiIyYSZ/hKhZs2bYvXu3NG1h8f8lT5w4EdHR0di4cSPUajXGjBmD/v3749ChQwCA/Px8BAUFwcXFBYcPH0Z6ejqCg4NhaWmJzz77rMrXhYiIiEyTyQciCwsLuLi4FJmflZWFVatWYd26dXjjjTcAAJGRkfD29saRI0fQvn177Nq1C+fPn8fu3bvh7OyMli1bYt68eZg2bRpmz54NKyurYl8zJycHOTk50rROp6uclSMiIiKTYNKnzADg8uXLcHNzQ4MGDTB06FCkpaUBABISEpCbmwt/f3+pbdOmTVGvXj3Ex8cDAOLj4+Hr6wtnZ2epTUBAAHQ6HZKSkkp8zYiICKjVaunh7u5eSWtHREREpsCkA1G7du0QFRWFmJgYLF++HKmpqXj99ddx//59aLVaWFlZwc7OTu85zs7O0Gq1AACtVqsXhgqXFy4rSXh4OLKysqTHtWvXKnbFiIiIyKSY9Cmznj17Sv9u0aIF2rVrBw8PD2zYsAE1a9astNdVKpVQKpWV1j8RERGZFpM+QvQsOzs7NGnSBFeuXIGLiwuePHmCzMxMvTYZGRnSNUcuLi5FvnVWOF3cdUlEREQkT9UqEGVnZyMlJQWurq5o3bo1LC0tsWfPHml5cnIy0tLSoNFoAAAajQZnz57FzZs3pTaxsbFQqVTw8fGp8vqJiIjINJn0KbMpU6agT58+8PDwwI0bNzBr1iyYm5tj8ODBUKvVCA0NxaRJk+Dg4ACVSoWxY8dCo9Ggffv2AIAePXrAx8cHw4YNw4IFC6DVajFjxgyEhYXxlBgRERFJTDoQXb9+HYMHD8adO3dQp04ddOzYEUeOHEGdOnUAAF988QXMzMwwYMAA5OTkICAgAN9++630fHNzc2zbtg3vv/8+NBoNatWqhZCQEMydO9dYq0REREQmSCGEEMYuwtTpdDqo1WpkZWVBpVIZuxySker4MwdELxv+dEf1VZ7P72p1DRERERFRZWAgIiIiItljICIiIiLZYyAiIiIi2TPpb5kRVSReoExERCXhESIiIiKSPQYiIiIikj0GIiIiIpI9BiIiIiKSPQYiIiIikj0GIiIiIpI9BiIiIiKSPQYiIiIikj0GIiIiIpI9BiIiIiKSPQYiIiIikj0GIiIiIpI9BiIiIiKSPQYiIiIikj0GIiIiIpI9C2MXQEREZMo8p0cbu4Ryuzo/yNglVDs8QkRERESyx0BEREREssdARERERLLHQERERESyx0BEREREssdARERERLLHQERERESyx0BEREREsscbM5JBquONyoiIiErCI0REREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHsMRERERCR7DEREREQkewxEREREJHuyujHjsmXLsHDhQmi1Wvj5+eHrr79G27ZtjV0WERFRhaqON8+9Oj/IqK8vmyNEP/30EyZNmoRZs2bh5MmT8PPzQ0BAAG7evGns0oiIiMjIFEIIYewiqkK7du3wt7/9Dd988w0AoKCgAO7u7hg7diymT59e6nN1Oh3UajWysrKgUqkqvLbqmOSJiIgqUmUcISrP57csTpk9efIECQkJCA8Pl+aZmZnB398f8fHxRdrn5OQgJydHms7KygLw18BWhoKch5XSLxERUXVRGZ+xhX2W5diPLALR7du3kZ+fD2dnZ735zs7OuHjxYpH2ERERmDNnTpH57u7ulVYjERGRnKm/rLy+79+/D7VaXWobWQSi8goPD8ekSZOk6YKCAty9exe1a9eGQqEwYmUvRqfTwd3dHdeuXauUU3/VEcdEH8dDH8ejKI6JPo6HPlMbDyEE7t+/Dzc3t+e2lUUgcnR0hLm5OTIyMvTmZ2RkwMXFpUh7pVIJpVKpN8/Ozq4yS6xSKpXKJDZUU8Ix0cfx0MfxKIpjoo/joc+UxuN5R4YKyeJbZlZWVmjdujX27NkjzSsoKMCePXug0WiMWBkRERGZAlkcIQKASZMmISQkBG3atEHbtm3x5Zdf4sGDB3j33XeNXRoREREZmWwC0cCBA3Hr1i3MnDkTWq0WLVu2RExMTJELrV9mSqUSs2bNKnI6UM44Jvo4Hvo4HkVxTPRxPPRV5/GQzX2IiIiIiEoii2uIiIiIiErDQERERESyx0BEREREssdARERERLLHQPSSuXv3LoYOHQqVSgU7OzuEhoYiOzu7xPZXr16FQqEo9rFx40apXXHL169fXxWr9ELKOx4A0KVLlyLrOnr0aL02aWlpCAoKgrW1NZycnDB16lTk5eVV5qpUiPKOx927dzF27Fh4eXmhZs2aqFevHsaNGyf9vl+h6rR9LFu2DJ6enqhRowbatWuHY8eOldp+48aNaNq0KWrUqAFfX19s375db7kQAjNnzoSrqytq1qwJf39/XL58uTJXoUKVZzxWrlyJ119/Hfb29rC3t4e/v3+R9sOHDy+yLQQGBlb2alSY8oxHVFRUkXWtUaOGXpvqvn0A5RuT4vafCoUCQUH//8OtJruNCHqpBAYGCj8/P3HkyBHx+++/i0aNGonBgweX2D4vL0+kp6frPebMmSNsbGzE/fv3pXYARGRkpF67R48eVcUqvZDyjocQQnTu3FmMGjVKb12zsrKk5Xl5eaJ58+bC399fnDp1Smzfvl04OjqK8PDwyl6dF1be8Th79qzo37+/2Lp1q7hy5YrYs2ePaNy4sRgwYIBeu+qyfaxfv15YWVmJ1atXi6SkJDFq1ChhZ2cnMjIyim1/6NAhYW5uLhYsWCDOnz8vZsyYISwtLcXZs2elNvPnzxdqtVr8+uuv4vTp0+Lvf/+7qF+/vkmu/7PKOx5DhgwRy5YtE6dOnRIXLlwQw4cPF2q1Wly/fl1qExISIgIDA/W2hbt371bVKr2Q8o5HZGSkUKlUeuuq1Wr12lTn7UOI8o/JnTt39Mbj3LlzwtzcXERGRkptTHUbYSB6iZw/f14AEMePH5fm7dixQygUCvG///2vzP20bNlSjBgxQm8eALF58+aKKrVKGDoenTt3FuPHjy9x+fbt24WZmZnejm/58uVCpVKJnJycCqm9MlTU9rFhwwZhZWUlcnNzpXnVZfto27atCAsLk6bz8/OFm5ubiIiIKLb922+/LYKCgvTmtWvXTvzzn/8UQghRUFAgXFxcxMKFC6XlmZmZQqlUih9//LES1qBilXc8npWXlydsbW3FmjVrpHkhISGib9++FV1qlSjveERGRgq1Wl1if9V9+xDixbeRL774Qtja2ors7GxpnqluIzxl9hKJj4+HnZ0d2rRpI83z9/eHmZkZjh49WqY+EhISkJiYiNDQ0CLLwsLC4OjoiLZt22L16tUQJn4LqxcZj7Vr18LR0RHNmzdHeHg4Hj58qNevr6+v3k09AwICoNPpkJSUVPErUkEqYvsAgKysLKhUKlhY6N/X1dS3jydPniAhIQH+/v7SPDMzM/j7+yM+Pr7Y58THx+u1B/76Wxe2T01NhVar1WujVqvRrl27Evs0FYaMx7MePnyI3NxcODg46M2Pi4uDk5MTvLy88P777+POnTsVWntlMHQ8srOz4eHhAXd3d/Tt21dvH1Cdtw+gYraRVatWYdCgQahVq5befFPcRmRzp2o50Gq1cHJy0ptnYWEBBwcHaLXaMvWxatUqeHt747XXXtObP3fuXLzxxhuwtrbGrl278MEHHyA7Oxvjxo2rsPormqHjMWTIEHh4eMDNzQ1nzpzBtGnTkJycjE2bNkn9PnuH88Lpso6zMVTE9nH79m3MmzcP7733nt786rB93L59G/n5+cX+7S5evFjsc0r6WxeOV+F/S2tjqgwZj2dNmzYNbm5ueh+YgYGB6N+/P+rXr4+UlBR89NFH6NmzJ+Lj42Fubl6h61CRDBkPLy8vrF69Gi1atEBWVhYWLVqE1157DUlJSahbt2613j6AF99Gjh07hnPnzmHVqlV68011G2EgqgamT5+Ozz//vNQ2Fy5ceOHXefToEdatW4dPPvmkyLKn57366qt48OABFi5caJQPvMoej6c/7H19feHq6opu3bohJSUFDRs2NLjfylJV24dOp0NQUBB8fHwwe/ZsvWWmtH1Q1Zg/fz7Wr1+PuLg4vQuJBw0aJP3b19cXLVq0QMOGDREXF4du3boZo9RKo9Fo9H4g/LXXXoO3tzf+/e9/Y968eUaszDSsWrUKvr6+aNu2rd58U91GGIiqgcmTJ2P48OGltmnQoAFcXFxw8+ZNvfl5eXm4e/cuXFxcnvs6P//8Mx4+fIjg4ODntm3Xrh3mzZuHnJycKv/Nmqoaj0Lt2rUDAFy5cgUNGzaEi4tLkW9ZZGRkAEC5+q0oVTEe9+/fR2BgIGxtbbF582ZYWlqW2t6Y20dJHB0dYW5uLv2tCmVkZJS4/i4uLqW2L/xvRkYGXF1d9dq0bNmyAquveIaMR6FFixZh/vz52L17N1q0aFFq2wYNGsDR0RFXrlwx6UD0IuNRyNLSEq+++iquXLkCoHpvH8CLjcmDBw+wfv16zJ0797mvYzLbiLEvYqKKU3jR7IkTJ6R5O3fuLPNFs507dy7y7aGSfPrpp8Le3t7gWqvCi45HoYMHDwoA4vTp00KI/7+o+ulvWfz73/8WKpVKPH78uOJWoIIZOh5ZWVmiffv2onPnzuLBgwdlei1T3T7atm0rxowZI03n5+eLV155pdSLqnv37q03T6PRFLmoetGiRdLyrKysanPRbHnHQwghPv/8c6FSqUR8fHyZXuPatWtCoVCILVu2vHC9lc2Q8XhaXl6e8PLyEhMnThRCVP/tQwjDxyQyMlIolUpx+/bt576GqWwjDEQvmcDAQPHqq6+Ko0ePioMHD4rGjRvrfa36+vXrwsvLSxw9elTveZcvXxYKhULs2LGjSJ9bt24VK1euFGfPnhWXL18W3377rbC2thYzZ86s9PV5UeUdjytXroi5c+eKEydOiNTUVLFlyxbRoEED0alTJ+k5hV+779Gjh0hMTBQxMTGiTp061eZr9+UZj6ysLNGuXTvh6+srrly5ovc12by8PCFE9do+1q9fL5RKpYiKihLnz58X7733nrCzs5O+MThs2DAxffp0qf2hQ4eEhYWFWLRokbhw4YKYNWtWsV+7t7OzE1u2bBFnzpwRffv2rTZfqy7veMyfP19YWVmJn3/+WW9bKLxFx/3798WUKVNEfHy8SE1NFbt37xatWrUSjRs3Nun/WShU3vGYM2eO2Llzp0hJSREJCQli0KBBokaNGiIpKUlqU523DyHKPyaFOnbsKAYOHFhkvilvIwxEL5k7d+6IwYMHCxsbG6FSqcS7776rdz+h1NRUAUDs27dP73nh4eHC3d1d5OfnF+lzx44domXLlsLGxkbUqlVL+Pn5iRUrVhTb1tSUdzzS0tJEp06dhIODg1AqlaJRo0Zi6tSpevchEkKIq1evip49e4qaNWsKR0dHMXnyZL2voZuq8o7Hvn37BIBiH6mpqUKI6rd9fP3116JevXrCyspKtG3bVhw5ckRa1rlzZxESEqLXfsOGDaJJkybCyspKNGvWTERHR+stLygoEJ988olwdnYWSqVSdOvWTSQnJ1fFqlSI8oyHh4dHsdvCrFmzhBBCPHz4UPTo0UPUqVNHWFpaCg8PDzFq1Kgi9+YxZeUZjwkTJkhtnZ2dRa9evcTJkyf1+qvu24cQ5X/PXLx4UQAQu3btKtKXKW8jCiFM7LuxRERERFWM9yEiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljIHoOhUKBX3/9VZq+ePEi2rdvjxo1akg/zlfcvJdJVFQU7OzsjF1GmT1b7+zZs03i7zJ8+HD069evwttWJ1qtFt27d0etWrUqfJsqy3ZqyLh6enriyy+/NLiuivS8dbx69SoUCgUSExOrrKZn95HGUJb1jouLg0KhQGZmJgDT3U8U59nay6Is69OlSxdMmDDhhWozlDG21UIlvY9k+Wv3w4cPx5o1awAAFhYWcHBwQIsWLTB48GAMHz4cZmb/nxPT09Nhb28vTc+aNQu1atVCcnIybGxsSpz3Mhk4cCB69epl7DIMNmXKFIwdO9bYZWDp0qUo643hn23bpUsXtGzZstI+mGfPno1ff/210ndOX3zxBdLT05GYmAi1Wl2pr1Wc8vwNqGye3Ue+iLi4OHTt2hX37t2r8MD82muvIT09vcTt7tn9xPDhw5GZmWn0sAc8v3aqGLIMRAAQGBiIyMhI5OfnIyMjAzExMRg/fjx+/vlnbN26FRYWfw2Ni4uL3vNSUlIQFBQEDw+PUueV15MnT2BlZWXw8ytTzZo1UbNmTWOXYTAbGxuTCKrl2ZlV1I7P1LarlJQUtG7dGo0bNza4jxdZJ36gVLxn95GmysrKqtRaTWU/UZzn1V7dPHnyxNglFEu2p8yUSiVcXFzwyiuvoFWrVvjoo4+wZcsW7NixA1FRUVK7pw8HKxQKJCQkYO7cuVAoFJg9e3ax8wDg2rVrePvtt2FnZwcHBwf07dsXV69elfotPHT/r3/9C25ubvDy8irX8xYtWgRXV1fUrl0bYWFhyM3Nldrk5ORg2rRpcHd3h1KpRKNGjbBq1Spp+blz59CzZ0/Y2NjA2dkZw4YNw+3bt0scq5IOLf/nP/+Bp6cn1Go1Bg0ahPv37z+3j23btsHLywvW1tZ466238PDhQ6xZswaenp6wt7fHuHHjkJ+fr7cuU6ZMwSuvvIJatWqhXbt2iIuLK9J3vXr1YG1tjTfffBN37tzRW/7soePjx4+je/fucHR0hFqtRufOnXHy5Em95ygUCnz//fd48803YW1tjcaNG2Pr1q0lrt9HH32Edu3aFZnv5+eHuXPnAih6uubnn3+Gr68vatasidq1a8Pf3x8PHjwo0nb48OHYv38/li5dCoVCAYVCobdNPM3T0xPz5s1DcHAwVCoV3nvvPQDAtGnT0KRJE1hbW6NBgwb45JNPpG0mKioKc+bMwenTp6X+C98DmZmZGDlyJOrUqQOVSoU33ngDp0+fll7v9OnT6Nq1K2xtbaFSqdC6dWucOHGixNp++eUX/PDDD1AoFBg+fDgAIC0tDX379oWNjQ1UKhXefvttZGRkSM8r/Pt9//33qF+/PmrUqFHi3wEAdu7cCW9vb9jY2CAwMBDp6enSsmf/Bvfv38fQoUNRq1YtuLq64osvvij2NMLDhw8xYsQI2Nraol69evjuu+9KrSEmJgYdO3aEnZ0dateujd69eyMlJUVaXni6YNOmTejatSusra3h5+eH+Ph4vX6et22Xxf79+9G2bVsolUq4urpi+vTpyMvLAwBs27YNdnZ20nsuMTERCoUC06dPl54/cuRIvPPOOyX2//Q+sizr9eeff6JPnz6wt7dHrVq10KxZM2zfvh1Xr15F165dAQD29vZ628jzxrPQxYsX8dprr6FGjRpo3rw59u/fLy173mmnp/cTs2fPxpo1a7BlyxbpPREXF4c33ngDY8aM0XverVu3YGVlhT179pTab2n7y4KCAkRERKB+/fqoWbMm/Pz88PPPP5da+8qVK+Hu7i5tG0uWLCn2qNrz9tN5eXkYM2YM1Go1HB0d8cknn+gdRb137x6Cg4Nhb28Pa2tr9OzZE5cvXy523Ap9+eWX8PT0lKZL+rwDgD/++KPU98Avv/yCZs2aQalUwtPTE4sXL9Zb/rz6gOe/jwr3Y7L8tfuQkBDRt2/fYpf5+fmJnj17StMAxObNm4UQQqSnp4tmzZqJyZMni/T0dHH//v1i5z158kR4e3uLESNGiDNnzojz58+LIUOGCC8vL5GTkyPVYGNjI4YNGybOnTsnzp07V+bnqVQqMXr0aHHhwgXx22+/CWtra/Hdd99JNb/99tvC3d1dbNq0SaSkpIjdu3eL9evXCyGEuHfvnqhTp44IDw8XFy5cECdPnhTdu3cXXbt2LXG8IiMjhVqtlqZnzZolbGxsRP/+/cXZs2fFgQMHhIuLi/joo49K7cPS0lJ0795dnDx5Uuzfv1/Url1b9OjRQ7z99tsiKSlJ/Pbbb8LKykqqVQghRo4cKV577TVx4MABceXKFbFw4UKhVCrFpUuXhBBCHDlyRJiZmYnPP/9cJCcni6VLlwo7O7si9fr5+UnTe/bsEf/5z3/EhQsXxPnz50VoaKhwdnYWOp1O7+9et25dsW7dOnH58mUxbtw4YWNjI+7cuVPs+p07d04AEFeuXCky7/Lly9LfrnC7u3HjhrCwsBBLliwRqamp4syZM2LZsmXSL88/3TYzM1NoNBoxatQokZ6eLtLT00VeXl6xdXh4eAiVSiUWLVokrly5ItUzb948cejQIZGamiq2bt0qnJ2dxeeffy6E+OvXpydPniyaNWsm9f/w4UMhhBD+/v6iT58+4vjx4+LSpUti8uTJonbt2tI4NGvWTLzzzjviwoUL4tKlS2LDhg0iMTGx2Npu3rwpAgMDxdtvvy3S09NFZmamyM/PFy1bthQdO3YUJ06cEEeOHBGtW7cWnTt31vv71apVSwQGBoqTJ0+K06dPF9t/4Tbm7+8vjh8/LhISEoS3t7cYMmSI1ObZ9/7IkSOFh4eH2L17tzh79qx48803ha2trRg/frzemDo4OIhly5aJy5cvi4iICGFmZiYuXrxYbB1CCPHzzz+LX375RVy+fFmcOnVK9OnTR/j6+or8/HwhhBCpqakCgGjatKnYtm2bSE5OFm+99Zbw8PAQubm5QoiybdvPKuz31KlTQgghrl+/LqytrcUHH3wgLly4IDZv3iwcHR2lX6fPzMwUZmZm4vjx40IIIb788kvh6Ogo2rVrJ/XZqFEjsXLlyhJf8+l9ZFnWKygoSHTv3l2cOXNGpKSkiN9++03s379f5OXliV9++UUAEMnJydI2Up7xrFu3rvj555/F+fPnxciRI4Wtra24ffu2EEKIffv2CQDi3r17Qoji92uF+4n79++Lt99+WwQGBkrviZycHLF27Vphb28vHj9+LD1vyZIlwtPTUxQUFBQ7PmXZX3766aeiadOmIiYmRqSkpIjIyEihVCpFXFxcsbUfPHhQmJmZiYULF4rk5GSxbNky4eDgUO79dOfOnYWNjY0YP368uHjxovjvf/9b5PPk73//u/D29hYHDhwQiYmJIiAgQDRq1Eg8efKkyLgV+uKLL4SHh4c0XdznXVm2lRMnTggzMzMxd+5ckZycLCIjI0XNmjVFZGRkmesry/uocD/GQPSMgQMHCm9vb2n66Te7EH8FpsKdSUnz/vOf/wgvLy+9N0hOTo6oWbOm2Llzp1SDs7OzFHTK8zwPDw+9D8R//OMfYuDAgUIIIZKTkwUAERsbW+z6zZs3T/To0UNv3rVr16SdUHGK23FYW1vrBYipU6fq7USL6+PZwPDPf/5TWFtbSyFACCECAgLEP//5TyGEEH/++acwNzcX//vf//T66tatmwgPDxdCCDF48GDRq1cvveUDBw4sNRA9Kz8/X9ja2orffvtNmgdAzJgxQ5rOzs4WAMSOHTtK7MfPz0/MnTtXmg4PD9cbk6e3u4SEBAFAXL16tdi+nt1GO3furPchXRIPDw/Rr1+/57ZbuHChaN26tTRd3Bj9/vvvQqVS6e38hRCiYcOG4t///rcQQghbW1sRFRX13Ncr1LdvXxESEiJN79q1S5ibm4u0tDRpXlJSkgAgjh07JtVmaWkpbt68WWrfxW1jy5YtE87OztL00+Oq0+mEpaWl2Lhxo7Q8MzNTWFtbFwlE77zzjjRdUFAgnJycxPLly8u83rdu3RIAxNmzZ4UQ//8B/v333xdZ7wsXLgghyrZtP+vZQPTRRx8V2acsW7ZM2NjYSGGiVatWYuHChUIIIfr16yf+9a9/CSsrK3H//n1x/fp1AUD6H5DiFBeISlsvX19fMXv27GL7evaDvyQljef8+fOlNrm5uaJu3bpS8C9PIBKi+M+JR48eCXt7e/HTTz9J81q0aFHi+hT2W9r+8vHjx8La2locPnxY73mhoaFi8ODBxdY+cOBAERQUpNd+6NCh5d5Pd+7cWXh7e+ttH9OmTZM+Ay9duiQAiEOHDknLb9++LWrWrCk2bNggvU5ZAtGzn3dl2VaGDBkiunfvrtf31KlThY+PT5nrK8v7qHA/JttTZiURQkChULxQH6dPn8aVK1dga2srnZd2cHDA48eP9Q7z+vr66l0LUdbnNWvWDObm5tK0q6srbt68CeCvQ97m5ubo3LlzibXt27dP6t/GxgZNmzYFgGIPQZfE09MTtra2xdZQEmtrazRs2FCadnZ2hqenp955e2dnZ6mfs2fPIj8/H02aNNGrd//+/VKtFy5cKHKqSqPRlFpHRkYGRo0ahcaNG0OtVkOlUiE7OxtpaWl67Vq0aCH9u1atWlCpVKWu49ChQ7Fu3ToAf21HP/74I4YOHVpsWz8/P3Tr1g2+vr74xz/+gZUrV+LevXul1l1Wbdq0KTLvp59+QocOHeDi4gIbGxvMmDGjyPo+6/Tp08jOzkbt2rX1xj81NVUa/0mTJmHkyJHw9/fH/Pnzy7UNAX/9/dzd3eHu7i7N8/HxgZ2dHS5cuCDN8/DwQJ06dZ7b37PbWGnb5R9//IHc3Fy0bdtWmqdWq/UO5xd6eltQKBRwcXEpdVu4fPkyBg8ejAYNGkClUkmnD0rbxlxdXQFA6teQbftZFy5cgEaj0dundejQAdnZ2bh+/ToAoHPnzoiLi4MQAr///jv69+8Pb29vHDx4EPv374ebm1u5r/kqbb3GjRuHTz/9FB06dMCsWbNw5syZ5/ZX1vF8enwsLCzQpk0bve3oRdWoUQPDhg3D6tWrAQAnT57EuXPnpFN7JSltf3nlyhU8fPgQ3bt313uf/fDDDyW+n5KTk/W2WwBFpp/3uoXat2+vt31oNBpcvnwZ+fn5uHDhAiwsLPS2w9q1a8PLy6vc4/rs512h570HOnTooNe+Q4cO5aqvLO+jwv2YbC+qLsmFCxdQv379F+ojOzsbrVu3xtq1a4sse3qnXqtWLYOeZ2lpqbdMoVCgoKAAAJ578XN2djb69OmDzz//vMiywo2xLEqroTzPKa2f7OxsmJubIyEhQS8AAnihix9DQkJw584dLF26FB4eHlAqldBoNEUu9CvvOg4ePBjTpk3DyZMn8ejRI1y7dg0DBw4stq25uTliY2Nx+PBh7Nq1C19//TU+/vhjHD169IW3v2e3q/j4eAwdOhRz5sxBQEAA1Go11q9fX+Rc/LOys7Ph6upa5JotANK1CrNnz8aQIUMQHR2NHTt2YNasWVi/fj3efPPNF1qHZz27TiUp7m8mKuBbZeXdFvr06QMPDw+sXLkSbm5uKCgoQPPmzUvdxgo/lJ73PqpoXbp0werVq3H69GlYWlqiadOm6NKlC+Li4nDv3r0S/+eqNKWt18iRIxEQEIDo6Gjs2rULERERWLx4canfBC3reFaFkSNHomXLlrh+/ToiIyPxxhtvPPcLNc/bzwFAdHQ0XnnlFb12SqXyhWo1ZD9dXmZmZkXeY09f01qopPewKbwHCvdjPEL0lL179+Ls2bMYMGDAC/XTqlUrXL58GU5OTmjUqJHeo7RvuRj6vKf5+vqioKBA72LCZ18jKSkJnp6eRV6jrB86VeXVV19Ffn4+bt68WaTWwm9ceHt74+jRo3rPO3LkSKn9Hjp0COPGjUOvXr2ki/VKu6i8rOrWrYvOnTtj7dq1WLt2Lbp37w4nJ6cS2ysUCnTo0AFz5szBqVOnYGVlhc2bNxfb1srKSu9i8/I4fPgwPDw88PHHH6NNmzZo3Lgx/vzzz+f236pVK2i1WlhYWBQZf0dHR6ldkyZNMHHiROzatQv9+/dHZGRkmWvz9vbGtWvXcO3aNWne+fPnkZmZCR8fH4PWt6waNGgAS0tLHD9+XJqXlZWFS5cuvVC/d+7cQXJyMmbMmIFu3brB29vboKN/hmzbxfURHx+v94F16NAh2Nraom7dugCA119/Hffv38cXX3whhZ/CQBQXF4cuXbqUu/bncXd3x+jRo7Fp0yZMnjwZK1euBADpCMLT22J5xvPp8cnLy0NCQgK8vb0NqrGk95yvry/atGmDlStXYt26dRgxYoRB/Rfy8fGBUqlEWlpakffZ00dOn+bl5aW33QIoMl1WxW1jjRs3hrm5Oby9vZGXl6fXpvDvUfj+rFOnDrRard42VlG37/D29sahQ4f05h06dAhNmjQpc31lfR81adJEvoEoJycHWq0W//vf/3Dy5El89tln6Nu3L3r37o3g4OAX6nvo0KFwdHRE37598fvvvyM1NRVxcXEYN26cdJi6Ip/3NE9PT4SEhGDEiBH49ddfpT42bNgAAAgLC8Pdu3cxePBgHD9+HCkpKdi5cyfeffddgz9wK0uTJk0wdOhQBAcHY9OmTUhNTcWxY8cQERGB6OhoAH8dfo+JicGiRYtw+fJlfPPNN4iJiSm138aNG+M///kPLly4gKNHj2Lo0KEVdluBoUOHYv369di4cWOJp8uAv3ZCn332GU6cOIG0tDRs2rQJt27dKnHn7enpiaNHj+Lq1au4fft2uf4PqnHjxkhLS8P69euRkpKCr776qkjw8vT0RGpqKhITE3H79m3k5OTA398fGo0G/fr1w65du3D16lUcPnwYH3/8MU6cOIFHjx5hzJgxiIuLw59//olDhw7h+PHj5foA8vf3h6+vL4YOHYqTJ0/i2LFjCA4ORufOnYs99VeRbG1tERISgqlTp2Lfvn1ISkpCaGgozMzMXui0ub29PWrXro3vvvsOV65cwd69ezFp0qRy92PItv2sDz74ANeuXcPYsWNx8eJFbNmyBbNmzcKkSZOk+63Z29ujRYsWWLt2rRR+OnXqhJMnT+LSpUsGHSEqzYQJE7Bz506kpqbi5MmT2Ldvn7TNeHh4QKFQYNu2bbh16xays7PLNZ7Lli3D5s2bcfHiRYSFheHevXsGBxZPT0+cOXMGycnJuH37tt5Rj5EjR2L+/PkQQrzw0VBbW1tMmTIFEydOxJo1a5CSkoKTJ0/i66+/lu6X96yxY8di+/btWLJkCS5fvox///vf2LFjh0HbbVpaGiZNmoTk5GT8+OOP+PrrrzF+/HgAf+07+vbti1GjRuHgwYM4ffo03nnnHbzyyivo27cvgL/C861bt7BgwQKkpKRg2bJl2LFjh+ED8pTJkydjz549mDdvHi5duoQ1a9bgm2++wZQpU8pc3/PeR0/vx2QbiGJiYuDq6gpPT08EBgZi3759+Oqrr7Bly5Yip2fKy9raGgcOHEC9evWk8/GhoaF4/PgxVCpVhT/vWcuXL8dbb72FDz74AE2bNsWoUaOkr3O7ubnh0KFDyM/PR48ePeDr64sJEybAzs5O74aUpiIyMhLBwcGYPHkyvLy80K9fPxw/fhz16tUD8Nf575UrV2Lp0qXw8/PDrl27MGPGjFL7XLVqFe7du4dWrVph2LBhGDduXKlHcsrjrbfewp07d/Dw4cNS74isUqlw4MAB9OrVC02aNMGMGTOwePFi9OzZs9j2U6ZMgbm5OXx8fFCnTp3nXv/ztL///e+YOHEixowZg5YtW+Lw4cP45JNP9NoMGDAAgYGB6Nq1K+rUqYMff/wRCoUC27dvR6dOnfDuu++iSZMmGDRoEP788084OzvD3Nwcd+7cQXBwMJo0aYK3334bPXv2xJw5c8pcm0KhwJYtW2Bvb49OnTrB398fDRo0wE8//VTmPl7EkiVLoNFo0Lt3b/j7+6NDhw7w9vZ+7lf7S2NmZob169cjISEBzZs3x8SJE7Fw4cJy92PItv2sV155Bdu3b8exY8fg5+eH0aNHIzQ0tEg/nTt3Rn5+vhSIHBwc4OPjAxcXl2KvqXoR+fn5CAsLg7e3NwIDA9GkSRN8++23Ur1z5szB9OnT4ezsjDFjxpRrPOfPn4/58+fDz88PBw8exNatW/WOZpbHqFGj4OXlhTZt2qBOnTp6RyoGDx4MCwsLDB48+IW2lULz5s3DJ598goiICGlcoqOjSzx93qFDB6xYsQJLliyBn58fYmJiMHHiRINqCQ4OxqNHj9C2bVuEhYVh/Pjx0u06gL/2wa1bt0bv3r2h0WgghMD27dulU13e3t749ttvsWzZMvj5+eHYsWNSYHlRrVq1woYNG7B+/Xo0b94cM2fOxNy5c/Wu2Xpefc97Hz29H1OIijjBTkT0Enjw4AFeeeUVLF68GKGhocYuh0zU1atX0bBhQxw/fhytWrUydjkA/gpwFy9exO+//27sUqotXlRNRLJ16tQpXLx4EW3btkVWVpZ0E83Cw+1ET8vNzcWdO3cwY8YMtG/f3qhhaNGiRdLvAu7YsQNr1qyRjrSRYRiIiEjWFi1ahOTkZFhZWaF169b4/fffDT7NQi+3Q4cOoWvXrmjSpInenaSN4dixY1iwYAHu37+PBg0a4KuvvsLIkSONWlN1x1NmREREJHumdxUtERERURVjICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZYyAiIiIi2WMgIiIiItljICIiIiLZ+z8Ip4ugbr8/1gAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "nullhypothesis_distribution_plot = plt.hist(shuffled_diffs_10000)\n", "\n", "plt.xlabel('Difference in median visit rates for high and low instability neighbourhoods')\n", "\n", "plt.ylabel('Frequency')" ] }, { "cell_type": "markdown", "id": "e699095a", "metadata": {}, "source": [ "## Step 8 - Compute the p-value\n", "\n", "a) Compute the number of simulated differences in medians in `shuffled_diffs_10000` that are greater than or equal to the observed median difference (`median_diff`). Store this value in `rightextreme`." ] }, { "cell_type": "code", "execution_count": 40, "id": "1a34886f", "metadata": {}, "outputs": [], "source": [ "rightextreme = shuffled_diffs_10000 >= median_diff\n", "\n", "rightextreme = rightextreme.sum()" ] }, { "cell_type": "markdown", "id": "8391adaf", "metadata": {}, "source": [ "b) Compute the number of simulated differences in medians in `shuffled_diffs_10000` that are less than the observed median difference (`median_diff`). Store this value in `leftextreme`." ] }, { "cell_type": "code", "execution_count": 41, "id": "50e101da", "metadata": {}, "outputs": [], "source": [ "leftextreme = shuffled_diffs_10000 < -median_diff\n", "\n", "leftextreme = leftextreme.sum()" ] }, { "cell_type": "markdown", "id": "117f46b6", "metadata": {}, "source": [ "c) Use `rightextreme` and `leftextreme` to compute the p-value. Store the p-value in `pvalue`. " ] }, { "cell_type": "code", "execution_count": 42, "id": "1cebcf55", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0064" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pvalue = (leftextreme + rightextreme) / 10000\n", "\n", "pvalue" ] }, { "cell_type": "markdown", "id": "b40e7f61", "metadata": {}, "source": [ "## Step 9 - Communicate what you did in the steps above\n", "\n", "\n", "a) In a few sentences introduce the question that you explored in this homework (see the [beginning](#gg274-homework-8-hypothesis-testing) of this homework). For example, why do you or others think this is an important question? (__1 mark__)" ] }, { "cell_type": "markdown", "id": "904406ca", "metadata": {}, "source": [ "> We investigated whether neighbourhoods whose residential statuses are relatively instable visit mental health services more or less frequently compared to those with relatively stable residential statuses. The analysis may tell us whether the marginalized neighbourhoods need more support for their mental health." ] }, { "cell_type": "markdown", "id": "dd1f4fd1", "metadata": {}, "source": [ "b) Briefly describe the data sources that you used to answer the question. Which statistical variables did you use and why did you use these varaibles? (__1 mark__)" ] }, { "cell_type": "markdown", "id": "331f2bb5", "metadata": {}, "source": [ "> We used the residential instabilty measure from the Ontario Marginalization Index and the age-standardized rates for mental health visits among those above age 20. We computed the mean of the instability measure to dicotomize the neighbourhoods to those whose measures are above the mean and those whose measures are below the mean. This allowed us to label the neighbourhoods into two groups - high vs. low residential instability. \n", ">\n", "> We then computed the difference in the median mental health visit rates between the two groups. The medians provided the \"typical\" (or centre of the distributions) visit rates of the two groups." ] }, { "cell_type": "markdown", "id": "1a773e29", "metadata": {}, "source": [ "c) What computational and statistical methods or analyses did you use to answer the question? Briefly describe these methods and how they were used to answer the question. (__1 mark__)" ] }, { "cell_type": "markdown", "id": "4d0390e8", "metadata": {}, "source": [ "> We simulated the distribution of the difference in median mental health visits rates under the assumption that there is no difference in the median rates betwee the two groups using the bootstrap. We then located the observed difference in the simulated distribution to assess high likely it would be to observe the value when there is actually no difference. When the location of the value in the distribution indicates that it is unlikely, we can reasonably conclude that the median mental health visit rates are different between the two groups.\n", ">\n", "> Note: Full marks if the data and method used are reasonably well explained between part b) and c)." ] }, { "cell_type": "markdown", "id": "a6bc5e37", "metadata": {}, "source": [ "d) Briefly describe the results of your statistical analysis in a few sentences. (__1 mark__)" ] }, { "cell_type": "markdown", "id": "2affe526", "metadata": {}, "source": [ "> After locating the observed difference in the simulated distribution, the resulting p-value is 0.0064 which is very small. It is reasonable to reject the hypothtesis that there is no difference between the two groups and conclude that the mental health visit rates are different between high and low residential instability neighbourhoods. Specifically, the visit rates are higher for those in neighbourhood with higher residential instability." ] }, { "cell_type": "markdown", "id": "edd25c8b", "metadata": {}, "source": [ "e) What conclusions can you draw about the question you set out to answer that is supported by the data and statistical analysis of the data? State at least one limitation of your conclusions. (See the [USC Research Guide section on study limitations](https://libguides.usc.edu/writingguide/limitations)) (__1 mark__)" ] }, { "cell_type": "markdown", "id": "8d121d9a", "metadata": {}, "source": [ "> The data shows that the typical mental health visit rates in Toronto neighbourhoods with relatively unstable residential statuses is higher. Our statistical analysis provides a reasonably strong evidence that there is a systematic difference and we didn't observe the difference by chance. \n", ">\n", "> Limitations may include...\n", ">\n", "> - the mental health visit rate data are from 2012 whereas the Ontario Marginalization Index is from 2006; a significant change in the neighbourhood composition may have altered the residential stability.\n", "> - the mental health visit rates are based on OHIP only (this requires some diggin on the data source: https://www.ontariohealthprofiles.ca/o_documents/aboutTheDataON/1_AboutTheData_AdultHealthDisease.pdf); other types of mental health visits may be relevant.\n", "> - ...\n", ">\n", "> Any reasonable limitation identified is acceptable." ] } ], "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.10.13" }, "vscode": { "interpreter": { "hash": "440dd12f919b48e435ef15d7652bb5c9f2f802a3e9de582e9da805c841a6f459" } } }, "nbformat": 4, "nbformat_minor": 5 }