{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Linear regression with the Associated Press\n",
"\n",
"In [this piece from the Associated Press](https://apnews.com/66ac44186b6249709501f07a7eab36da), Nicky Forster combines from the US Census Bureau and the CDC to see how life expectancy is related to actors like unemployment, income, and others. We'll be looking at how they can write sentences like this one: \n",
"\n",
"> \"An increase of 10 percentage points in the unemployment rate in a neighborhood translated to a loss of roughly a year and a half of life expectancy, the AP found.\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<p class=\"reading-options\">\n <a class=\"btn\" href=\"/ap-regression-unemployment/simple-regression-with-census-data-statsmodels-with-formulas\">\n <i class=\"fa fa-sm fa-book\"></i>\n Read online\n </a>\n <a class=\"btn\" href=\"/ap-regression-unemployment/notebooks/Simple regression with census data (statsmodels with formulas).ipynb\">\n <i class=\"fa fa-sm fa-download\"></i>\n Download notebook\n </a>\n <a class=\"btn\" href=\"https://colab.research.google.com/github/littlecolumns/ds4j-notebooks/blob/master/ap-regression-unemployment/notebooks/Simple regression with census data (statsmodels with formulas).ipynb\" target=\"_new\">\n <i class=\"fa fa-sm fa-laptop\"></i>\n Interactive version\n </a>\n</p>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Prep work: Downloading necessary files\n",
"Before we get started, we need to download all of the data we'll be using.\n",
"* **R12221544_SL140.csv:** Employment status by census tract - ACS 2015 5-year, tract level, table B23025\n",
"* **R12221544.txt:** Data dictionary for R12221544_SL140.csv - Description of columns and tables\n",
"* **R12221550_SL140.csv:** Demographic data by census tract - ACS 2015 5-year, tract level, tables B23025, B06009, B03002, B19013, C17002\n",
"* **R12221550.txt:** Data dictionary for R12221550_SL140.csv - Description of columns and tables\n",
"* **US_A.CSV:** USALEEP national data - estimates of life expectancy at birth for most of the census tracts in the USA\n",
"* **Record_Layout_CensusTract_Life_Expectancy.pdf:** USALEEP data dictionary - None\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"# Make data directory if it doesn't exist\n",
"!mkdir -p data\n",
"!wget -nc https://nyc3.digitaloceanspaces.com/ml-files-distro/v1/ap-regression-unemployment/data/R12221544_SL140.csv.zip -P data\n",
"!unzip -n -d data data/R12221544_SL140.csv.zip\n",
"!wget -nc https://nyc3.digitaloceanspaces.com/ml-files-distro/v1/ap-regression-unemployment/data/R12221544.txt -P data\n",
"!wget -nc https://nyc3.digitaloceanspaces.com/ml-files-distro/v1/ap-regression-unemployment/data/R12221550_SL140.csv.zip -P data\n",
"!unzip -n -d data data/R12221550_SL140.csv.zip\n",
"!wget -nc https://nyc3.digitaloceanspaces.com/ml-files-distro/v1/ap-regression-unemployment/data/R12221550.txt -P data\n",
"!wget -nc https://nyc3.digitaloceanspaces.com/ml-files-distro/v1/ap-regression-unemployment/data/US_A.CSV -P data\n",
"!wget -nc https://nyc3.digitaloceanspaces.com/ml-files-distro/v1/ap-regression-unemployment/data/Record_Layout_CensusTract_Life_Expectancy.pdf -P data"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import numpy as np\n",
"\n",
"# Keep decimal numbers to 4 decimal places\n",
"pd.set_option(\"display.float_format\", '{:,.4f}'.format)\n",
"pd.set_option(\"display.max_columns\", 100)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reading in our data\n",
"\n",
"We'll start by reading in our datasets (as we must!).\n",
"\n",
"### Life expectancy data\n",
"\n",
"The first dataset is [USALEEP](https://www.cdc.gov/nchs/nvss/usaleep/usaleep.html), the U.S. Small-area Life Expectancy Estimates Project. It's from the National Center for Health Statistics at the CDC.\n",
"\n",
"We're going to rename a few columns so they make a little more sense. The columns themselves aren't too complicated, mostly pointing out the census tracts and the life expectancy."
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>tract_id</th>\n",
" <th>STATE2KX</th>\n",
" <th>CNTY2KX</th>\n",
" <th>TRACT2KX</th>\n",
" <th>life_expectancy</th>\n",
" <th>life_expectancy_std_err</th>\n",
" <th>flag</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1001020100</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>20100</td>\n",
" <td>73.1000</td>\n",
" <td>2.2348</td>\n",
" <td>3</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1001020200</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>20200</td>\n",
" <td>76.9000</td>\n",
" <td>3.3453</td>\n",
" <td>3</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1001020400</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>20400</td>\n",
" <td>75.4000</td>\n",
" <td>1.0216</td>\n",
" <td>3</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>1001020500</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>20500</td>\n",
" <td>79.4000</td>\n",
" <td>1.1768</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1001020600</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>20600</td>\n",
" <td>73.1000</td>\n",
" <td>1.5519</td>\n",
" <td>3</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" tract_id STATE2KX CNTY2KX TRACT2KX life_expectancy \\\n",
"0 1001020100 1 1 20100 73.1000 \n",
"1 1001020200 1 1 20200 76.9000 \n",
"2 1001020400 1 1 20400 75.4000 \n",
"3 1001020500 1 1 20500 79.4000 \n",
"4 1001020600 1 1 20600 73.1000 \n",
"\n",
" life_expectancy_std_err flag \n",
"0 2.2348 3 \n",
"1 3.3453 3 \n",
"2 1.0216 3 \n",
"3 1.1768 1 \n",
"4 1.5519 3 "
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"life_expec = pd.read_csv(\"data/US_A.CSV\")\n",
"life_expec.columns = ['tract_id', 'STATE2KX','CNTY2KX', 'TRACT2KX', 'life_expectancy', \n",
" 'life_expectancy_std_err', 'flag']\n",
"life_expec.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`flag` sounds a little mysterious and scary, but it's just whether the numbers were observed or predicted."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Unemployment Data\n",
"\n",
"For unemployment, we'll be using data from the Census. You can get unemployment measures from a lot of places in the US government, so the Census Bureau actually [publishes data on which ones you should use](https://www.census.gov/topics/employment/labor-force/guidance.html). In our case we're going to collect unemployment data from the American Community Survery, for the 5-year period nearest to our life expectancy dataset.\n",
"\n",
"We'll be using Table B23025005, Employment Status for the Population 16 Years and Over. Take a look at [the codebook](https://www.socialexplorer.com/data/ACS2015_5yr/metadata/?ds=ACS15_5yr&var=B23025005) for more details. We're only reading in a few columns to keep things looking clean."
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Geo_FIPS</th>\n",
" <th>total_pop</th>\n",
" <th>unemployed</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1001020100</td>\n",
" <td>1554</td>\n",
" <td>54</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1001020200</td>\n",
" <td>1731</td>\n",
" <td>116</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1001020300</td>\n",
" <td>2462</td>\n",
" <td>91</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>1001020400</td>\n",
" <td>3424</td>\n",
" <td>216</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1001020500</td>\n",
" <td>8198</td>\n",
" <td>221</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Geo_FIPS total_pop unemployed\n",
"0 1001020100 1554 54\n",
"1 1001020200 1731 116\n",
"2 1001020300 2462 91\n",
"3 1001020400 3424 216\n",
"4 1001020500 8198 221"
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"columns = {\n",
" 'Geo_FIPS': 'Geo_FIPS',\n",
" 'ACS15_5yr_B23025001': 'total_pop',\n",
" 'ACS15_5yr_B23025005': 'unemployed'\n",
"}\n",
"employment = pd.read_csv(\"data/R12221544_SL140.csv\", usecols=columns.keys(), encoding='latin-1')\n",
"\n",
"employment = employment.rename(columns=columns)\n",
"employment.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the important thing is _percent_ unemployed, not just how many people are in an area, we'll need to do a little calculation.\n",
"\n",
"> **Note that I'm multiplying by 100 here** - if you have a litle extra time, try running this notebook on your own with a 0-1.0 percentage instead of the 0-100 version."
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Geo_FIPS</th>\n",
" <th>total_pop</th>\n",
" <th>unemployed</th>\n",
" <th>pct_unemployment</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1001020100</td>\n",
" <td>1554</td>\n",
" <td>54</td>\n",
" <td>3.4749</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1001020200</td>\n",
" <td>1731</td>\n",
" <td>116</td>\n",
" <td>6.7013</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1001020300</td>\n",
" <td>2462</td>\n",
" <td>91</td>\n",
" <td>3.6962</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>1001020400</td>\n",
" <td>3424</td>\n",
" <td>216</td>\n",
" <td>6.3084</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1001020500</td>\n",
" <td>8198</td>\n",
" <td>221</td>\n",
" <td>2.6958</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Geo_FIPS total_pop unemployed pct_unemployment\n",
"0 1001020100 1554 54 3.4749\n",
"1 1001020200 1731 116 6.7013\n",
"2 1001020300 2462 91 3.6962\n",
"3 1001020400 3424 216 6.3084\n",
"4 1001020500 8198 221 2.6958"
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"employment['pct_unemployment'] = employment['unemployed'] / employment['total_pop'] * 100\n",
"employment.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Demographic data\n",
"\n",
"The other dataset we're using is also from the Census Bureau's American Community Survey. It's a collection of many more columns with very, very long names and very, very mysterious codes. The tables include:\n",
"\n",
"* [B03002](https://censusreporter.org/tables/B03002/): Hispanic or Latino Origin by Race\n",
"* [B06009](https://censusreporter.org/tables/B06009/): Place of Birth by Educational Attainment in the United States\n",
"* [C17002](https://censusreporter.org/tables/C17002/): Ratio of Income to Poverty Level\n",
"* [B19013](https://censusreporter.org/tables/B19013/): Median Household Income\n",
"\n",
"Again, we're only picking a few columns to read in. Generally they're the raw count of certain populations, as well as the total population counted of each table."
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Geo_FIPS</th>\n",
" <th>race_table_total</th>\n",
" <th>white</th>\n",
" <th>black</th>\n",
" <th>hispanic</th>\n",
" <th>edu_table_total</th>\n",
" <th>less_than_hs</th>\n",
" <th>poverty_table_total</th>\n",
" <th>poverty_100_124</th>\n",
" <th>poverty_125_149</th>\n",
" <th>income</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1001020100</td>\n",
" <td>1948</td>\n",
" <td>1703</td>\n",
" <td>150</td>\n",
" <td>17</td>\n",
" <td>1,243.0000</td>\n",
" <td>184.0000</td>\n",
" <td>1948</td>\n",
" <td>81</td>\n",
" <td>101</td>\n",
" <td>61,838.0000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1001020200</td>\n",
" <td>2156</td>\n",
" <td>872</td>\n",
" <td>1149</td>\n",
" <td>17</td>\n",
" <td>1,397.0000</td>\n",
" <td>356.0000</td>\n",
" <td>1983</td>\n",
" <td>232</td>\n",
" <td>58</td>\n",
" <td>32,303.0000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1001020300</td>\n",
" <td>2968</td>\n",
" <td>2212</td>\n",
" <td>551</td>\n",
" <td>0</td>\n",
" <td>2,074.0000</td>\n",
" <td>221.0000</td>\n",
" <td>2968</td>\n",
" <td>148</td>\n",
" <td>207</td>\n",
" <td>44,922.0000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>1001020400</td>\n",
" <td>4423</td>\n",
" <td>3662</td>\n",
" <td>162</td>\n",
" <td>464</td>\n",
" <td>2,899.0000</td>\n",
" <td>339.0000</td>\n",
" <td>4423</td>\n",
" <td>141</td>\n",
" <td>182</td>\n",
" <td>54,329.0000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1001020500</td>\n",
" <td>10763</td>\n",
" <td>7368</td>\n",
" <td>2674</td>\n",
" <td>80</td>\n",
" <td>6,974.0000</td>\n",
" <td>310.0000</td>\n",
" <td>10563</td>\n",
" <td>256</td>\n",
" <td>1064</td>\n",
" <td>51,965.0000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Geo_FIPS race_table_total white black hispanic edu_table_total \\\n",
"0 1001020100 1948 1703 150 17 1,243.0000 \n",
"1 1001020200 2156 872 1149 17 1,397.0000 \n",
"2 1001020300 2968 2212 551 0 2,074.0000 \n",
"3 1001020400 4423 3662 162 464 2,899.0000 \n",
"4 1001020500 10763 7368 2674 80 6,974.0000 \n",
"\n",
" less_than_hs poverty_table_total poverty_100_124 poverty_125_149 \\\n",
"0 184.0000 1948 81 101 \n",
"1 356.0000 1983 232 58 \n",
"2 221.0000 2968 148 207 \n",
"3 339.0000 4423 141 182 \n",
"4 310.0000 10563 256 1064 \n",
"\n",
" income \n",
"0 61,838.0000 \n",
"1 32,303.0000 \n",
"2 44,922.0000 \n",
"3 54,329.0000 \n",
"4 51,965.0000 "
]
},
"execution_count": 67,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"columns = {\n",
" 'Geo_FIPS': 'Geo_FIPS',\n",
" 'ACS15_5yr_B03002001': 'race_table_total',\n",
" 'ACS15_5yr_B03002004': 'black',\n",
" 'ACS15_5yr_B03002003': 'white',\n",
" 'ACS15_5yr_B03002012': 'hispanic',\n",
" 'ACS15_5yr_B06009001': 'edu_table_total',\n",
" 'ACS15_5yr_B06009002': 'less_than_hs',\n",
" 'ACS15_5yr_C17002004': 'poverty_100_124',\n",
" 'ACS15_5yr_C17002005': 'poverty_125_149',\n",
" 'ACS15_5yr_C17002001': 'poverty_table_total',\n",
" 'ACS15_5yr_B19013001': 'income'\n",
"} \n",
"census = pd.read_csv(\"data/R12221550_SL140.csv\", usecols=columns.keys(), encoding='latin-1')\n",
"census = census.rename(columns=columns)\n",
"census.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We'll again convert these **raw population numbers into percentages** - what percent of people are certain races? What percent of people have not finished high school?"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Geo_FIPS</th>\n",
" <th>pct_black</th>\n",
" <th>pct_white</th>\n",
" <th>pct_hispanic</th>\n",
" <th>pct_less_than_hs</th>\n",
" <th>pct_under_150_poverty</th>\n",
" <th>income</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1001020100</td>\n",
" <td>7.7002</td>\n",
" <td>87.4230</td>\n",
" <td>0.8727</td>\n",
" <td>14.8029</td>\n",
" <td>9.3429</td>\n",
" <td>61,838.0000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1001020200</td>\n",
" <td>53.2931</td>\n",
" <td>40.4453</td>\n",
" <td>0.7885</td>\n",
" <td>25.4832</td>\n",
" <td>14.6243</td>\n",
" <td>32,303.0000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1001020300</td>\n",
" <td>18.5647</td>\n",
" <td>74.5283</td>\n",
" <td>0.0000</td>\n",
" <td>10.6557</td>\n",
" <td>11.9609</td>\n",
" <td>44,922.0000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>1001020400</td>\n",
" <td>3.6627</td>\n",
" <td>82.7945</td>\n",
" <td>10.4906</td>\n",
" <td>11.6937</td>\n",
" <td>7.3027</td>\n",
" <td>54,329.0000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1001020500</td>\n",
" <td>24.8444</td>\n",
" <td>68.4567</td>\n",
" <td>0.7433</td>\n",
" <td>4.4451</td>\n",
" <td>12.4964</td>\n",
" <td>51,965.0000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Geo_FIPS pct_black pct_white pct_hispanic pct_less_than_hs \\\n",
"0 1001020100 7.7002 87.4230 0.8727 14.8029 \n",
"1 1001020200 53.2931 40.4453 0.7885 25.4832 \n",
"2 1001020300 18.5647 74.5283 0.0000 10.6557 \n",
"3 1001020400 3.6627 82.7945 10.4906 11.6937 \n",
"4 1001020500 24.8444 68.4567 0.7433 4.4451 \n",
"\n",
" pct_under_150_poverty income \n",
"0 9.3429 61,838.0000 \n",
"1 14.6243 32,303.0000 \n",
"2 11.9609 44,922.0000 \n",
"3 7.3027 54,329.0000 \n",
"4 12.4964 51,965.0000 "
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"census_features = pd.DataFrame({\n",
" 'Geo_FIPS': census.Geo_FIPS,\n",
" 'pct_black': census.black / census.race_table_total * 100,\n",
" 'pct_white': census.white / census.race_table_total * 100,\n",
" 'pct_hispanic': census.hispanic / census.race_table_total * 100,\n",
" 'pct_less_than_hs': census.less_than_hs / census.edu_table_total * 100,\n",
" 'pct_under_150_poverty': (census.poverty_100_124 + census.poverty_125_149) / census.poverty_table_total * 100,\n",
" 'income': census.income,\n",
"})\n",
"census_features.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Merge the data\n",
"\n",
"To perform our regression, we need all of our data in **a single dataframe**. Fortunately these dataframes all keep track of census tracks with [FIPS codes](https://transition.fcc.gov/oet/info/maps/census/fips/fips.txt), unique geographic codes that we can rely on to merge our datasets."
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>tract_id</th>\n",
" <th>STATE2KX</th>\n",
" <th>CNTY2KX</th>\n",
" <th>TRACT2KX</th>\n",
" <th>life_expectancy</th>\n",
" <th>life_expectancy_std_err</th>\n",
" <th>flag</th>\n",
" <th>Geo_FIPS</th>\n",
" <th>total_pop</th>\n",
" <th>unemployed</th>\n",
" <th>pct_unemployment</th>\n",
" <th>pct_black</th>\n",
" <th>pct_white</th>\n",
" <th>pct_hispanic</th>\n",
" <th>pct_less_than_hs</th>\n",
" <th>pct_under_150_poverty</th>\n",
" <th>income</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1001020100</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>20100</td>\n",
" <td>73.1000</td>\n",
" <td>2.2348</td>\n",
" <td>3</td>\n",
" <td>1001020100</td>\n",
" <td>1554</td>\n",
" <td>54</td>\n",
" <td>3.4749</td>\n",
" <td>7.7002</td>\n",
" <td>87.4230</td>\n",
" <td>0.8727</td>\n",
" <td>14.8029</td>\n",
" <td>9.3429</td>\n",
" <td>61,838.0000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1001020200</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>20200</td>\n",
" <td>76.9000</td>\n",
" <td>3.3453</td>\n",
" <td>3</td>\n",
" <td>1001020200</td>\n",
" <td>1731</td>\n",
" <td>116</td>\n",
" <td>6.7013</td>\n",
" <td>53.2931</td>\n",
" <td>40.4453</td>\n",
" <td>0.7885</td>\n",
" <td>25.4832</td>\n",
" <td>14.6243</td>\n",
" <td>32,303.0000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1001020400</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>20400</td>\n",
" <td>75.4000</td>\n",
" <td>1.0216</td>\n",
" <td>3</td>\n",
" <td>1001020400</td>\n",
" <td>3424</td>\n",
" <td>216</td>\n",
" <td>6.3084</td>\n",
" <td>3.6627</td>\n",
" <td>82.7945</td>\n",
" <td>10.4906</td>\n",
" <td>11.6937</td>\n",
" <td>7.3027</td>\n",
" <td>54,329.0000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>1001020500</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>20500</td>\n",
" <td>79.4000</td>\n",
" <td>1.1768</td>\n",
" <td>1</td>\n",
" <td>1001020500</td>\n",
" <td>8198</td>\n",
" <td>221</td>\n",
" <td>2.6958</td>\n",
" <td>24.8444</td>\n",
" <td>68.4567</td>\n",
" <td>0.7433</td>\n",
" <td>4.4451</td>\n",
" <td>12.4964</td>\n",
" <td>51,965.0000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1001020600</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>20600</td>\n",
" <td>73.1000</td>\n",
" <td>1.5519</td>\n",
" <td>3</td>\n",
" <td>1001020600</td>\n",
" <td>2855</td>\n",
" <td>190</td>\n",
" <td>6.6550</td>\n",
" <td>11.9190</td>\n",
" <td>72.9161</td>\n",
" <td>13.0615</td>\n",
" <td>17.4873</td>\n",
" <td>10.8543</td>\n",
" <td>63,092.0000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" tract_id STATE2KX CNTY2KX TRACT2KX life_expectancy \\\n",
"0 1001020100 1 1 20100 73.1000 \n",
"1 1001020200 1 1 20200 76.9000 \n",
"2 1001020400 1 1 20400 75.4000 \n",
"3 1001020500 1 1 20500 79.4000 \n",
"4 1001020600 1 1 20600 73.1000 \n",
"\n",
" life_expectancy_std_err flag Geo_FIPS total_pop unemployed \\\n",
"0 2.2348 3 1001020100 1554 54 \n",
"1 3.3453 3 1001020200 1731 116 \n",
"2 1.0216 3 1001020400 3424 216 \n",
"3 1.1768 1 1001020500 8198 221 \n",
"4 1.5519 3 1001020600 2855 190 \n",
"\n",
" pct_unemployment pct_black pct_white pct_hispanic pct_less_than_hs \\\n",
"0 3.4749 7.7002 87.4230 0.8727 14.8029 \n",
"1 6.7013 53.2931 40.4453 0.7885 25.4832 \n",
"2 6.3084 3.6627 82.7945 10.4906 11.6937 \n",
"3 2.6958 24.8444 68.4567 0.7433 4.4451 \n",
"4 6.6550 11.9190 72.9161 13.0615 17.4873 \n",
"\n",
" pct_under_150_poverty income \n",
"0 9.3429 61,838.0000 \n",
"1 14.6243 32,303.0000 \n",
"2 7.3027 54,329.0000 \n",
"3 12.4964 51,965.0000 \n",
"4 10.8543 63,092.0000 "
]
},
"execution_count": 72,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"merged = life_expec.merge(employment, left_on='tract_id', right_on='Geo_FIPS')\n",
"merged = merged.merge(census_features, left_on='Geo_FIPS', right_on='Geo_FIPS')\n",
"merged.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Running the regression\n",
"\n",
"Using the `statsmodels` package, we'll run a linear regression to find the relationship between **life expectancy** and our **calculated columns**. Note that we're using the **formula method** of writing a regression instead of the dataframes method. I find it both more readable and more useable than the dataframes method.\n",
"\n",
"> Python note: I'm also using a [multiline string](https://stackoverflow.com/posts/10660443/revisions) which is started and ended using three quotation marks. Multi-line strings and formulas are a match made in heaven!"
]
},
{
"cell_type": "code",
"execution_count": 74,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>OLS Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>life_expectancy</td> <th> R-squared: </th> <td> 0.490</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>OLS</td> <th> Adj. R-squared: </th> <td> 0.490</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>Least Squares</td> <th> F-statistic: </th> <td> 8997.</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Thu, 12 Dec 2019</td> <th> Prob (F-statistic):</th> <td> 0.00</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>21:50:03</td> <th> Log-Likelihood: </th> <td>-1.6208e+05</td>\n",
"</tr>\n",
"<tr>\n",
" <th>No. Observations:</th> <td> 65656</td> <th> AIC: </th> <td>3.242e+05</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Df Residuals:</th> <td> 65648</td> <th> BIC: </th> <td>3.243e+05</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Df Model:</th> <td> 7</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>nonrobust</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>t</th> <th>P>|t|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> 81.2365</td> <td> 0.122</td> <td> 665.628</td> <td> 0.000</td> <td> 80.997</td> <td> 81.476</td>\n",
"</tr>\n",
"<tr>\n",
" <th>pct_black</th> <td> -0.0666</td> <td> 0.001</td> <td> -56.960</td> <td> 0.000</td> <td> -0.069</td> <td> -0.064</td>\n",
"</tr>\n",
"<tr>\n",
" <th>pct_white</th> <td> -0.0386</td> <td> 0.001</td> <td> -36.707</td> <td> 0.000</td> <td> -0.041</td> <td> -0.037</td>\n",
"</tr>\n",
"<tr>\n",
" <th>pct_hispanic</th> <td> 0.0131</td> <td> 0.001</td> <td> 10.298</td> <td> 0.000</td> <td> 0.011</td> <td> 0.016</td>\n",
"</tr>\n",
"<tr>\n",
" <th>pct_less_than_hs</th> <td> -0.0862</td> <td> 0.002</td> <td> -48.979</td> <td> 0.000</td> <td> -0.090</td> <td> -0.083</td>\n",
"</tr>\n",
"<tr>\n",
" <th>pct_under_150_poverty</th> <td> -0.0596</td> <td> 0.003</td> <td> -21.738</td> <td> 0.000</td> <td> -0.065</td> <td> -0.054</td>\n",
"</tr>\n",
"<tr>\n",
" <th>income</th> <td> 4.825e-05</td> <td> 5.8e-07</td> <td> 83.217</td> <td> 0.000</td> <td> 4.71e-05</td> <td> 4.94e-05</td>\n",
"</tr>\n",
"<tr>\n",
" <th>pct_unemployment</th> <td> -0.1490</td> <td> 0.004</td> <td> -33.408</td> <td> 0.000</td> <td> -0.158</td> <td> -0.140</td>\n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <th>Omnibus:</th> <td>2114.193</td> <th> Durbin-Watson: </th> <td> 1.520</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Prob(Omnibus):</th> <td> 0.000</td> <th> Jarque-Bera (JB): </th> <td>4788.035</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Skew:</th> <td> 0.183</td> <th> Prob(JB): </th> <td> 0.00</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Kurtosis:</th> <td> 4.271</td> <th> Cond. No. </th> <td>7.06e+05</td>\n",
"</tr>\n",
"</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 7.06e+05. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: life_expectancy R-squared: 0.490\n",
"Model: OLS Adj. R-squared: 0.490\n",
"Method: Least Squares F-statistic: 8997.\n",
"Date: Thu, 12 Dec 2019 Prob (F-statistic): 0.00\n",
"Time: 21:50:03 Log-Likelihood: -1.6208e+05\n",
"No. Observations: 65656 AIC: 3.242e+05\n",
"Df Residuals: 65648 BIC: 3.243e+05\n",
"Df Model: 7 \n",
"Covariance Type: nonrobust \n",
"=========================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"-----------------------------------------------------------------------------------------\n",
"Intercept 81.2365 0.122 665.628 0.000 80.997 81.476\n",
"pct_black -0.0666 0.001 -56.960 0.000 -0.069 -0.064\n",
"pct_white -0.0386 0.001 -36.707 0.000 -0.041 -0.037\n",
"pct_hispanic 0.0131 0.001 10.298 0.000 0.011 0.016\n",
"pct_less_than_hs -0.0862 0.002 -48.979 0.000 -0.090 -0.083\n",
"pct_under_150_poverty -0.0596 0.003 -21.738 0.000 -0.065 -0.054\n",
"income 4.825e-05 5.8e-07 83.217 0.000 4.71e-05 4.94e-05\n",
"pct_unemployment -0.1490 0.004 -33.408 0.000 -0.158 -0.140\n",
"==============================================================================\n",
"Omnibus: 2114.193 Durbin-Watson: 1.520\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 4788.035\n",
"Skew: 0.183 Prob(JB): 0.00\n",
"Kurtosis: 4.271 Cond. No. 7.06e+05\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"[2] The condition number is large, 7.06e+05. This might indicate that there are\n",
"strong multicollinearity or other numerical problems.\n",
"\"\"\""
]
},
"execution_count": 74,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import statsmodels.formula.api as smf\n",
"\n",
"model = smf.ols(\"\"\"\n",
" life_expectancy ~ pct_black + pct_white + pct_hispanic + pct_less_than_hs\n",
" + pct_under_150_poverty + income + pct_unemployment\n",
"\"\"\", data=merged)\n",
"\n",
"results = model.fit()\n",
"results.summary()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It's a lot to take in, so I'll zero us in on the important parts:\n",
" \n",
"<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAT0AAACTCAIAAABzmFngAAAYZmlDQ1BEaXNwbGF5AABYhZV5Z1QTbbT1mcwkISGE3lvovXeQ3nvvCAQIvRk6WGgiRUUQUQQsqCiIWCkqdhQLIqhgQ0WwgKCiAoqg3B+gb/vWd9c9PzJn7WeffdqsmZUEQDiXnpwcj+MASEhMZXrYWtD8/ANoxFGggDhwAwCNHp6SbO7m5gQA8Pv6T5sbAgQA4IEKPTk5/r/n/1/jjGCkhAMgbgAQFpESngCAnAHACsOTmakABH0AkMpITU4FIAQBAA/Tzz8AgJAMADxRy34hAPCELfu7AICH6eVhCUBoAWCh0OnMKABqJwDQ0sOjUgGojwEIXIkRMYkAPFMABJPwaHoEgLAyACgnJCRFAAj7AYB82N90ov6hGfZHk06P+uMv9wIAACxWMSnJ8fSs/+M4/ndLiE/7nUMWACjRTDsPAOABQB7HJTl6AAAFAJlKDHNxBQAuAOR7TATAso8jR6fZeS/zcSLhKZYBAMAHgFOPoFs5AoAIAM4mMd7FaQUPi4yxsQcADgBcZkyqvddKbDEjxdpzRbOOmeTh+tuPZFqar8S20pkAK/zutDhv8xX9x9EM+9/6M9nRXr7LNaPk9BgfFwCgAqB8KXGejsscVDo72tLlN4eZ5uENANIAqCEj0dZiWR8NjmTaeKzwmQkpv/tFi6Nj7F1W/N2p0V52Kzot4XRrTwAQAEA7GYnm3r91GCl+Tr97iWBYWS/3jvYzEr1X+kVHklMtPFZivyTHu63wMTIj3tYDACQBMJGUdM+VWMwklem1siPMJTnVzWu5Tiwslu7gtlwPlglOYAlWQIM0oEEYJEEsxPRNdUwBbeXEBujAhChggMoK8jvCF+jAhESggydkwwdIBAak/ImzADowgQHpkAg//6DLnyoQCXRgQjowIAXi4C0wIQEcIR4YkAZMYEDin2w+8BqYEPOf7OGQBPGQBEyI+X9g5mAJTitI2m9dGvtvJsGaYEWwI9gQFDAhzAQzwpwwE8wMM8E0MX3M4He1f/Hxb/ED+FH8IH4E/yQkJp/5r35o4AwjkLYyKwaE/b1nTBbTxHQwC8wYM8EMgIbxYUKggmlj+pg5ZooZYTqYAViuVJ4G/9X+Rw9/m/oKj6ROwpH4SWYk+X9HUhWpOn9UGJD4jwkt1xr2Z66Wf07+nd/yb5OOgCRw/DcTLUZPoz3oFfQW2oV2AA29hHaivegFtONvd9FrYELUn2wewIBEiIN4iPlPPvpKTiYwIEX9qPqE+uLyWSojMxUAwDIpOYsZExWdSjNPTo5n0OwTw1WVaZrqGvoAfv4BtOXH1FcPQAAA4bv3F8YYB1jFAUDq/wuL3Q7QfB2Av/QvTDYQQFAZ4OT98DRm+jKGAQDggQzswAOCIAZSIA8qoAm6YARmYA0O4Ape4A/BEA7RkABMyIC1kAdFUAbbYAfshr1wAI7AMTgFHdAFV+AG3IF+GIRnMAJvYBKmYQ4WEAQhImwINyKIiCMyiBKiiegjJog14oR4IP5IKBKFJCJpyFqkAClDKpHdyH6kCTmJnEWuILeQAeQJ8gqZQL4gP3AojoLjwYniZHFqOH2cOc4R54VbjYvCrcFl4wpxW3G7cA24Flw77gruDm4QN4KbxM2igLKifKgEqoLqo5aoKxqARqJMdD1ailajDWgreg7tQR+gI+gUOo8RMG6MhqlgRpgd5o2FY2uw9dhmbDd2BGvHurEH2CtsGvuFZ8OL4JXwhnh7vB8+Cp+BL8JX4xvxbfjr+EH8G/wcgUDgI8gR9Ah2BH9CLCGHsJlQTzhOuEwYIIwRZolEoiBRiWhMdCXSianEImINsYV4iXif+Ib4nYWVRZxFk8WGJYAlkSWfpZqlmeUiy32WdywLJA6SDMmQ5EqKIGWRykkHSedI90hvSAtkTrIc2ZjsRY4l55F3kVvJ18nD5K+srKySrAas7qwxrLmsu1hPsN5kfcU6T+GiKFIsKUGUNMpWymHKZcoTylc2NjZZNjO2ALZUtq1sTWzX2F6wfadyU1Wp9tQI6gZqLbWdep/6kZ3ELsNuzh7Mns1ezX6a/R77FAeJQ5bDkoPOsZ6jluMsxyOOWU5uTg1OV84Ezs2czZy3OMe5iFyyXNZcEVyFXAe4rnGNcaPcUtyW3OHcBdwHua9zv+Eh8Mjx2PPE8pTxHOPp45nm5eLV5vXhzeSt5b3AO8KH8sny2fPF85XzneIb4vvBL8pvzs/gL+Fv5b/P/01AWMBMgCFQKnBcYFDghyBN0FowTrBCsEPwuRAmpCjkLpQhtEfoutCUMI+wkXC4cKnwKeGnIjgRRREPkRyRAyK9IrOiYqK2osmiNaLXRKfE+MTMxGLFqsQuik2Ic4ubiMeIV4lfEn9P46WZ0+Jpu2jdtGkJEQk7iTSJ/RJ9EguScpLekvmSxyWfS5Gl9KUipaqkrkpNS4tLO0uvlT4q/VSGJKMvEy2zU6ZH5pusnKyv7CbZDtlxOQE5e7lsuaNyw/Js8qbya+Qb5B8qEBT0FeIU6hX6FXGKOorRirWK95RwSrpKMUr1SgPKeGUD5UTlBuVHKhQVc5V0laMqr1T5VJ1U81U7VD+qSasFqFWo9aj9UtdRj1c/qP5Mg0vDQSNf45zGF01FzXDNWs2HWmxaNlobtDq1PmsraTO092g/1uHWcdbZpHNV56euni5Tt1V3Qk9aL1SvTu+RPo++m/5m/ZsGeAMLgw0GXQbzhrqGqYanDD8ZqRjFGTUbja+SW8VYdXDVmLGkMd14v/GICc0k1GSfyYiphCndtMF01EzKLMKs0eyduYJ5rHmL+UcLdQumRZvFN0tDy3WWl61QK1urUqs+ay5rb+vd1i9sJG2ibI7aTNvq2ObYXrbD2znaVdg9she1D7dvsp920HNY59DtSHH0dNztOOqk6MR0OueMc3Zw3u487CLjkujS4Qqu9q7bXZ+7ybmtcTvvTnB3c691f+uh4bHWo8eT2zPEs9lzzsvCq9zrmbe8d5r3VR92nyCfJp9vvla+lb4jfmp+6/zu+Av5x/h3BhADfAIaA2YDrQN3BL4J0gkqChpaLbc6c/WtYKHg+OALIewh9JDTofhQ39Dm0EW6K72BPhtmH1YXNh1uGb4zfDLCLKIqYoJhzKhkvIs0jqyMHI8yjtoeNRFtGl0dPRVjGbM75nOsXeze2G9xrnGH45bifeOPJ7AkhCacTeRKjEvsThJLykwaSFZKLkoeWWO4ZseaaaYjszEFSVmd0pnKk5qc2psmn7Yx7VW6SXpt+vcMn4zTmZyZiZm9WYpZJVnvsm2yD+VgOeE5V9dKrM1b+2qd+br965H1YeuvbpDaULjhTa5t7pE8cl5c3t189fzK/JkC34JzhaKFuYVjG203Hi2iFjGLHm0y2rS3GCuOKe4r0SqpKflVGlF6u0y9rLpscXP45ttbNLbs2rK0NXJrX7lu+Z5thG2J24YqTCuOVHJWZleObXfe3l5FqyqtmtkRsuNWtXb13p3knWk7R3Y57eqska7ZVrO4O3r3YK1F7fE6kbqSum/1EfX395jtad0rurds7499Mfse77fd394g21B9gHAg/cDbgz4Hew7pH2pqFGosa/x5OPHwyBGPI91Nek1NzSLN5UdxR9OOTrQEtfQfszrW2arSuv843/GyE3Ai7cT7k6Enh045nrp6Wv906xmZM3Vt3G2l7Uh7Vvt0R3THSKd/58BZh7NXzxmdazuvev5wl0RX7QXeC+UXyRcLLy5dyr40ezn58tSVqCtjV0OuPrvmd+1ht3t333XH6zdv2Ny41mPec+mm8c2uW4a3zt7Wv91xR/dOe69Ob9tdnbttfbp97ff07nX2G/SfG1g1cPG+6f0rD6we3Hho//DOoMvgwJD30ONHQY9GHkc8Hn8S/+Tz0/SnC89yh/HDpc85nle/EHnR8FLh5fER3ZELr6xe9Y56jj4bCx+bfJ3yevFN4Vu2t9XvxN81jWuOd03YTPS/D3z/ZjJ5cmGq6APnh7qP8h/PfDL71DvtN/3mM/Pz0pfNXwW/Hp7Rnrk66zb7Yi5hbuFb6XfB70fm9ed7fvj+eLeQsUhc3PVT4ee5X46/hpcSlpaS6Uw6AACgAICLjAT4chiAzR+Aux+AHLj8PW/FUAQABwA+iDXOHNXHBPBkAgtRncWfVEC+RCGw0akdHGTOeK7bPDq8dfwgECfYJ6wrsk10UtyMVi4xIEWWNpDxl42TS5APUrBQFFX8rHRDuUYlTtVYjU3tpfpxjVxNdy0JrQ/aZ3U26rrriei90W81yDQ0NyIbPVhVZxxhomzyxbTDbK25hQXF4qXlRatm63qbCtv1dnR7UwcBh8+OvU6tzvUu+1273MY88J6CXkLeHD6oz6Lvgj8EkAKpQWyrsdWzwaMh/aGX6afDGsNrIkoZWZFRUV7RFjHasYpxEvGCCeyJaOJM0mhy/5rzzIMpW1M3pBWlt2ViWYzsy2thnex6ww32uYF5aflbC3YU5mzU3jhWVL7JrVimhLUUynCbObfIbzUpd9nmWxFQGbDdr8pnh1e1+06XXY41trstak3qDOq19qjsVdynvt+xoeDAyCH7xpbDk02czTJHNVqMjlm1Oh/3PRFyMvpU8umMM+vb8ts3dhR3lp0tP7fjfF1X44UzF69fenR55MrQ1ePXIrsFum9er76R0RN5c/Ut39vudxx7be/a9XndW9O/b+DJA9aHaoOWQ/aPrB/rP5F5Sn06/2x8+PHzKy8OvCwYiXrlPeoy5vza9Y3rW4d3BuP84yMTpe+1349MHpnK/mD3keVj0yfbT2PTBz5nfgn+6jrjPBs7d/X7ph8dP62Wllb2r4Fi6AQ2gh8jTLOgJF1yNGsdZYSqyJ7BcYNLkDuL5yGfJn++wHMhHeEikX4xIXE/WoVEl+Sw1Kz0nMx72btyB+SZCiaKLIoPlfYqx6roqPxSvaG2Vd1XQ1zjnWarVrq2sQ6ic123VM9Vn1t/yKDGMNBI1Gh4VZ1xkImgySPTnWaB5rLmCxaDlietNlszbFbZctq+teuy3+GQ7shwCnOOdklyTXALc3f1MPJU9BL2pvrgfOZ83/kN+V8LaA2sDSpdnR0cE+IXakVXCxMIR8LfRwwyuiPbohqjq2MKY5Pi/OPNEuQS2RJnkl4lv1gzkyKRGpJWk34l43HmWNZU9vxa1nVi6+U30HIJuS/z2vLLC5iFwRu9i/w2xRQXlNSXHitr29y+5czWk+XHtjVVHKrct722aseO8uqSnfm7smqSdkfVxtTl1l/aq7DvSIPcgcqDDw7NH6YeEWqSalY8qtGid8yk1eq48wn/k/Gnik4fOHOxbaD9Rcd459dz6Hn+LqULRhfNLuldlriCuzJ6tedaW/fh67U3tvVsvJl9i3k79U5Jb1cf3711/c/vCz0wfeg1GDmU++jQ43tPZp5xDas8d3qR/HLnyPlX90dfjI2+nnyLf2cwnjcxMMk5pf5B56PsJ/ZP36fffn705fbXszP7ZzfM+XyT+zb3vWs++4fRAmXR6ufEyv5VkUlcPRqMKeCJ+M+ECeJ7llHSZ1YyRYbNnBrAnsfRwjnAtcQjw2vNF8u/UWCv4Bmh68I3RW6InhfbL55Js6D9kDgo6Sg5KVUsLSd9VSZYZl62Sk5d7rZ8lAJR4bCineI7pSJleeXrKuGqoFqvtkrtsXqaBpfGcU0nzXGtAm0x7U4dD50p3Y164nod+m764wYbDPkMjxqZG91fFb7qo3GOCdGk1lTbdMgs21zMvNPC1eKJZbTlklWDtZsNyeaa7Vo7bbv39g0OQY4CjkNOO5w9XdhdbrkWuBm5zbgf94jzlPN87bXfe7WPoM9D33I/O78l/7aA+EDpwOdB1atdVs8FV4XIhJwJNQ99Ss8Mkwx7HF4TEc2wjdSLMoi2j6HHJsTR400TOBKGEw8lJSTrJC+uucYsTXFL5U19lrY3PSJDNuNt5p4s66zh7PgcnpwHa8+vu7i+e8O13LN5TfnVBQWFSRsDi6w3KRbjix+W1JQGlEmXLWwe2XJ369nyfdvWVwRWGm4X2j5fNbTjVPXOnVt2Vdbs33269kbd4/r3exb2se2nNWgdsDsYdCipcf3hkiObm3Kb6Uf1WqgtX459aJ0/QTkpdkrztNuZnLYz7d87Dc4mn6s5f6Kr88L5i7cuzV6xvXq22/P6bE/1La3bD3u39IX22983f2gxFP+EOjw52vd+dmZ+aQlg+fc+AACCLsD2PAC/IgBvU4CKbgDZQQB+MoAbG4CXAeBkIwFH6QXE8NWf9wcCKBCAFThBAMRBDtTBEKzAFQIgElIgD8phD7TCRbgHr2AGISIiiAZii4QgGUgF0oLcRN7iCDh5nBMuBVePu4NbQvXRdPQs+guzxbZjo3gtfDH+JcGQUENYINKJt1n0WA6ThEkVZFZyCSuZdRtFiHKYTZuti2pMPceuz36ew47jGWcqFwfXMW4r7gEeL54BXlfe+3whfN/5awSMBV4IrhMSFjonHCxCEukSTRfTFvsqforGlNCRWJTskaqWjpZZJUuVHZE7LV+sEKZoriSrTFVeUPmo+lptUL1NI0dTQ/OFVrG2jvYnnU7dSr0s/QgDJ0N1I/5VVGNVk1ozJfMtFrcsP1mz2PDaCtqJ2Es7aDu6OK1x3uXS7frFXcrD13OrV48P5mvlV+TfG8gXFLa6Ofh1KIHOGUYImw1/EzHMeB/FHu0YsyP2XfyqhKrEj8kOa5pTKKlr0p5m2GR2ZqvkNK6jra/N5curKCAX5m2c3RRbPFlatjlha1sF53ahqg/VTbtCdvPV9tdv2Wu7b7ah/CDPoeLGuSNxTV+ObjtmfZzzxOdTb8+Mt092vjs31vX5Ev8Vy2vB10N7PG+Z3lG7q3BPdyDxwfdH2FPS8N6X3K8uvqGOr500/3D808IX3RmbOfK3Ld9vz4//eLPwZPHMz22/wpbUl5ZW9k8ECnCBIEiAImiBMdiBF4RCAuRACdRAE5yFO/AcphE8IoSoI7ZICJKFVCEnkD7kA44dp4ULwBXgTuHeoOJoCHoQncJ0sUJsEK+Az8MPEwwJtUQgRhMHWaxZOklqpGayArmFVZv1EsWNMsaWSSVRd7BLsJ/gsOR4xpnBxcfVwe3D/YFnHS+ZdxefCt9t/iQBfoHLgjFCPEKXhZNEpEWGRWvE/MQFxJ/Q6iUiJNWlQOqh9FGZQtkgOW15Nvn3Cr2Kp5X2KperFKiuVUtVD9cw06Ro9mmVajvq8Ot81n2i16PfbtBguNkoe1W6cYlJp+k3cy2LCMsyq0brdpvztuftLtjfcnjlhHNWdPFx3ejW4T7lKe0V4l3v88JP0j82oD2IuNo3eHfI9dAB+tWwpvDiiBiGR6RdlH90fszlOLb4sISuJKHk7DXPUyxSm9LZM5Iz72RL5KSv7V+vs+FgnnB+VSF5Y07RVDG9ZLQse4t6OW7b88qTVenV2ju/1JysTas33PNjX2OD5oH6g+8a5Q5HHznRzH+07phx64cTNacMTve10dsXOhvOuXfBhaZLTpc/X93bHXbD8KbEbezO3bvp9wj9pfcpDxoGQx45P4l/dvj5uxHxUbfXeW8vTvBPbvsoO333a9Xc5nn7Bc3FPT9f//q8sn8MSMABAiABSqALluAGwZAA66Ac9sMZuAkvYAahILKIGbIayUFqkQvIKxwJp4Wj43bg+lFelIFewESwXOw93h9/l2BJuEA0Jl5hcWJ5Tkohs5NPsPpQUEoH2xqqBvU7+3WOGs40Ln9uex4HXnc+B349AQVBHaEQ4SyRVNEwMS9xF5qzhLOkk5SztIdMiGyK3Bb5ZoWbihPKbCp6qpFqu9WHNIW0IrSP6yzouenfNSxZ5W+CN91mtmjhaFlg1WjdYdNle9Guz37B0dGp3UXVtcVd1aPdy9J7yDfBnxzQEuQTzBnKGhYSEch4HWUUXRbzNs4jvjfROen+mkDmeGpOuljGi6wbOZfX1W/wzv2Rv7/Qp0h803TJhbLNWyLLbSsEK+9URe6Y21lQw7m7oU63/u7eyP1IQ91B/UODh9OahJtvtmxotT2hdsrmzIb2hs7yc/5d/BceXaq94n+N2H3ohnbP+VvWtx/1Zvap9aMD0w/GBwceVTyRe1r/7Ndz6xelL++8Yh/1Htv3euKtxru48X0TN9+/n8J/EPmo/slq2vcz/UvEV7cZyZnZ2S1zInPN3wy+7f42/933e/s83zxzvn1+4YfZj8IftxaoC54LOxf6F1kWzRYzF08uTvyU+On/s/Ln7Z8/f2n8ivi189edX7+WNJYYS7uWepeWAFIitTQBAAChWADgXywtfZUFIFYC/KxYWlpoWFr6eQAAHQa4HL/8HxIAAIEDoG4aAOBO549c+Jf9DwXo2fFLJl3tAAAACXBIWXMAAAsTAAALEwEAmpwYAAAFUWlUWHRYTUw6Y29tLmFkb2JlLnhtcAAAAAAAPD94cGFja2V0IGJlZ2luPSLvu78iIGlkPSJXNU0wTXBDZWhpSHpyZVN6TlRjemtjOWQiPz4gPHg6eG1wbWV0YSB4bWxuczp4PSJhZG9iZTpuczptZXRhLyIgeDp4bXB0az0iQWRvYmUgWE1QIENvcmUgNS42LWMxNDUgNzkuMTYzNDk5LCAyMDE4LzA4LzEzLTE2OjQwOjIyICAgICAgICAiPiA8cmRmOlJERiB4bWxuczpyZGY9Imh0dHA6Ly93d3cudzMub3JnLzE5OTkvMDIvMjItcmRmLXN5bnRheC1ucyMiPiA8cmRmOkRlc2NyaXB0aW9uIHJkZjphYm91dD0iIiB4bWxuczpleGlmPSJodHRwOi8vbnMuYWRvYmUuY29tL2V4aWYvMS4wLyIgeG1sbnM6eG1wPSJodHRwOi8vbnMuYWRvYmUuY29tL3hhcC8xLjAvIiB4bWxuczpkYz0iaHR0cDovL3B1cmwub3JnL2RjL2VsZW1lbnRzLzEuMS8iIHhtbG5zOnBob3Rvc2hvcD0iaHR0cDovL25zLmFkb2JlLmNvbS9waG90b3Nob3AvMS4wLyIgeG1sbnM6eG1wTU09Imh0dHA6Ly9ucy5hZG9iZS5jb20veGFwLzEuMC9tbS8iIHhtbG5zOnN0RXZ0PSJodHRwOi8vbnMuYWRvYmUuY29tL3hhcC8xLjAvc1R5cGUvUmVzb3VyY2VFdmVudCMiIGV4aWY6UGl4ZWxYRGltZW5zaW9uPSIzMTciIGV4aWY6UGl4ZWxZRGltZW5zaW9uPSIxNDciIHhtcDpDcmVhdGVEYXRlPSIyMDE5LTEyLTEyVDIxOjQyOjA2LTA1OjAwIiB4bXA6TW9kaWZ5RGF0ZT0iMjAxOS0xMi0xMlQyMTo0Mjo0Ny0wNTowMCIgeG1wOk1ldGFkYXRhRGF0ZT0iMjAxOS0xMi0xMlQyMTo0Mjo0Ny0wNTowMCIgZGM6Zm9ybWF0PSJpbWFnZS9wbmciIHBob3Rvc2hvcDpDb2xvck1vZGU9IjMiIHBob3Rvc2hvcDpJQ0NQcm9maWxlPSJEaXNwbGF5IiB4bXBNTTpJbnN0YW5jZUlEPSJ4bXAuaWlkOjdlYjg4NDJiLTk3OTctNGMzNC04MDNiLTFiZGFlMTRmM2Y0YSIgeG1wTU06RG9jdW1lbnRJRD0ieG1wLmRpZDo3ZWI4ODQyYi05Nzk3LTRjMzQtODAzYi0xYmRhZTE0ZjNmNGEiIHhtcE1NOk9yaWdpbmFsRG9jdW1lbnRJRD0ieG1wLmRpZDo3ZWI4ODQyYi05Nzk3LTRjMzQtODAzYi0xYmRhZTE0ZjNmNGEiPiA8eG1wTU06SGlzdG9yeT4gPHJkZjpTZXE+IDxyZGY6bGkgc3RFdnQ6YWN0aW9uPSJzYXZlZCIgc3RFdnQ6aW5zdGFuY2VJRD0ieG1wLmlpZDo3ZWI4ODQyYi05Nzk3LTRjMzQtODAzYi0xYmRhZTE0ZjNmNGEiIHN0RXZ0OndoZW49IjIwMTktMTItMTJUMjE6NDI6NDctMDU6MDAiIHN0RXZ0OnNvZnR3YXJlQWdlbnQ9IkFkb2JlIFBob3Rvc2hvcCBDQyAyMDE5IChNYWNpbnRvc2gpIiBzdEV2dDpjaGFuZ2VkPSIvIi8+IDwvcmRmOlNlcT4gPC94bXBNTTpIaXN0b3J5PiA8L3JkZjpEZXNjcmlwdGlvbj4gPC9yZGY6UkRGPiA8L3g6eG1wbWV0YT4gPD94cGFja2V0IGVuZD0iciI/Pn4qA/QAAErNSURBVHic7Z15XI3Z/8A/z22xdG+LlotqVLRMktRQXFmTLCENBoVokMlSwyBjbBWmRItJWWKiGsmSnZAlidHcFEVSUdo3995u3e35/XGm53fndrsq15Lv8355ed3OeZ7zfM55zueczznP83w+GI7jANDY2AhfBX/++aeXl1dkZOScOXPYbLampqZ4bk1NjaqqqpKSEpEiFAqbmpoKCgokWsDMzExdXV32tfh8fkVFRc+ePXv16iXj0iKRqKKigkajUanUDtWFxWIBAI1G69BZJP8LKH5uAT4W3bp169atm0SihBoDgIKCgoqKyrNnzyorK8XTtbW136u3SkpKenp67700hULp06dP+yUnIDWWpC2wr2y+bWhoqKys1NHRUVNT+9+5NMn/Gl+b3pKQ/C9A+dwCkJCQdBhSb0lIuh6k3pKQdD1IvSUh6XqQektC0vX4dz+ZhISkC0HOtyQkXQ9Sbz8OOA6kIUPy0fhq33P8bAiFwOFAQwPweNC9O6iqQs+eoKDwucX6CKCBCcM+txz/i5B6K1eam+HVK0hPh9xc4HCARgNzcxg+HAwNQexjhq4NjgOfDywWsFggEkHPnqCmBt27kwr8KSH1Vn4IBPD8OURHw4UL8OYN4DgoKEC/fjBtGixZAt9+C5SuvyrBcaivh6dP4dEjKC4GgQB0dMDKCmxsgE4HRbI7fSLIhpYftbUQFwcxMUC87C0QwKtXcPAgUKng4wMaGp9Vvg8GKe3t2xAbC/fuQXU1AICKCgwaBLNng4sL6OuTs+6ngdRbOSEUQmEhJCSAxBcaOA4cDiQkgItLl9dbPh+ePYPYWLh4EZqb/01ksSAjA7hc0NCAmTOhg98Yk3SOrm+5fSGIRPD8Obx+LSULx6GoCAoKQCT65GK1i3fv3nl5eenr6w8ePDg6Orr1AcHBwaampjampsmbNuGpqf+vtAihEJ4/hytX3jx+PH/+fH19/UWLFj179gwABAKBv7+/qalp//79N27cyOPx5Ch2fn7+1KlT6XT6yJEj09LS2pN75cqVYcOG0en02bNnv3nz5uNJuHTp0s2bN6PfT58+nTx5MpIkNTVV4sgHDx6MHz+eTqc7OTnl5uYCwF9//TXyv7x586apqWnTpk39+/c3NTXdtWsX4F2KpKSk3bt3NzQ0SM01MTEBkH+NHj58uHv37pycHBzH3d3dAeDFixeSBzU14cHBOIbhAFL+USj4/v24QNC68DYL/IR4eXkxGIwnT55cvnwZAK5fvy6em5CQQKPRbt28mXP+/OlevURSK4hhQkvLGTo6gQEBOTk5vr6+DAYDx/GgoCATE5OHDx8+fPjQ3Nzc399fXjILhUITE5OVK1e+ePFi7969AFBVVSU799WrVwAQExOTl5fn7u7+8SQ8ffo0AKxZswb9aWRktG7duvz8fH9/fwAQ771o7IiOjs7Pz9+8ebOenh6HwykpKUltYdWqVebm5iwWy9fX187OLisr6969e0ZGRl1Mb6dPnw4Ab968kZprZGT0MfQ2LCwMAOLj4/EWNcvLy5M8qKkJ37GjTb3FMDwkBBcKWxfeZoGfEBqNdvXqVfR76dKlK1asEM+dNWuWn58fLhLhz569tbKSXkGARk3Nlf374yIRn8/n8/n5+flCoXDSpElRUVGonG3btrm4uMhL5idPngDAu3fv0J96enonT56UnfvHH38MHToUpRQXFwNAaWmp3CUsKyszMTH5+eefkd6KRCIajYaGwszMTACoqakhDj558qSRkRH6LRKJAODOnTtEbmZmpo6OTmlpKY7jRkZGRAX37Nnz+e3ks2fP2tjYbNmyxdTUlE6n+/j48Pl8ABAIBFu2bOnfv7+qqurMmTNLSkpWr1597tw5AJg4ceKVK1dkF3v48OGBAweqqqrOmTOnpqYGAG7cuDF48GAMw+h0+vbt23Ecb25uXrZsmaqqKoZho0aNevr0aetyjh49iobJdevWBQYGosSkpCR9fX1VVdU9e/YAAJfL9fHx2env39a7FjiO3755U8abGBIFFhcXOzo6Yhimqqo6f/78d+/etaMhOwmXy2WxWPr6+uhPAwODsrIy8QPevn1rYGAAGAY8nooMU5/L1afTHSdOVFJSGjhwYHl5OYVCSUxMXLRoEQAIBIILFy4MGjRIXmJXVlbSaDTClc+3335bUVEhO7eiosLY2Bil6OrqAkBFRYV8JcRxfOnSpTt37tTW1kYpGIZt377dw8Nj+/btbm5u3t7eyBsZgkqlVlVVIeO8vr4eAEpLS4nctWvX+vr69u3bFwDU1dWJClZWVn5+va2trc3MzNy+ffvw4cO1tLT27du3e/duAAgMDNy+fbumpua4cePOnDnzww8/WFhYoBmVwWCgyrRFYmKip6enioqKk5PTyZMn0SDq4uJSWFj4xx9/GBsbb9my5fr167GxsdHR0TNmzAgMDLx7966Hh0fror755htzc3MAGDhw4LfffosSN23aNGDAABaLtXbt2rdv3wYHBx+IjDTo31+GSOcuXHj7X30QR6JAPz+/69ev79ixY9asWXFxcahB5EhZWdmjFpqamgBAseURjpKSksQaj8fjKSoqokEHa2vowTAhjj+4f3+Cg8OLFy8cHR1dXV25XK6KioqysnJdXd3s2bPLysq8vb3lVQWhUNijRw/iT0VFRXGxpeby+XzCJaCCggIACASCD5RQvCX5fP7Ro0cVFRVnzpwpfkxxcXFJScmDBw/evHlTXV3N5XKJrOHDhwOAl5fXrVu3Fi5cCABE7p07d27evLl8+XL058yZM/38/JKSko4ePbp79+7Pr7eINWvWHD169Pbt2wAQExMDAMHBwTQaLTU19ezZs97e3jo6Om5ubmg4/O233ywtLWWUFhoaCgAHDx4MDAy0t7e/e/duVlYWi8VisVg4jgcGBl6+fNnU1BRNLLW1tUOHDr106RIxnYozbtw4dBsWLVrk4uKCEmNjY2/duuXq6goAhYWFw4cPj4qOnjZtmgyRRACFhYVt5UoU+PbtWwCoq6tbvHjxhQsXvv/+e5mN12HOnTs3rIXq6moAQNoLAFwuV1lZWfxgZWXlf3MpFFFbeovjFBzv27fvunXrjI2Nd+/eXVlZibamMjIyLC0tGxoaMjIyiFnow1FUVBRXAC6XK+6LT2qusrIy4Y8J2XSoph8ioXhLMpnMxYsXe3h4ZGZmvn79+u3bty9fvnz9+vW+ffuYTOalS5devHhx8eJF8a0pdXX1a9euvXr1ytPTc+jQoZaWloTrwvj4+Hnz5hG+ytauXbtixYpffvklLi4uICDg869vDx8+DACxsbHoTzS5oc40btw4iYPbub5t7Qnx6tWrISEhxJ+WlpYvX7588+aNtbU1kRgQECC12Nbr26ysLBzHN2zYAAA3btzIzMycMGrUr20s/HAAEYAPwM3/7vcgpBZ448YNHR0dJJWOjs6lS5c62bjtg0ajnTlzBv12c3NbtWqVeO68efN8fHxwkQh/8qR8wIC26shVUVk7YAAuEuE4jrQiPT09KysLAPbu3SuUtrb/ENCgUFlZieO4SCTS0dFJSkqSnRsdHW1ubo4OyMvLA4Dy8nI5SpiVlaXTArp3P/zww/Xr1wGAx+OhY8aMGbNv3z7ilIaGhpSUFJFIhOM4h8MheoJAIKDRaImJicSRmZmZr169Qr+/IL3dtWsX3iIujUZDN55Ysv/+++/Lly/ncDhIb4uLi6UWRegt2lj++++/8/LyMjIyLl26VFRUlJ6e/uzZs8uXL48bNw4A1qxZk5ubm56enpubi7YcAaCioqJ1sUhv4+Li8P9uI/n5+SE1Gzp0aDeAp3PmtNWncQBfgJspKa0Ll1rggwcPmExmRkbGihUrAIDobR8Jd3d3e3v72tpaJpMJAGfPnsVx/NChQ8nJyTiOHz58mEaj5b94UZuamtmzZ1sVFGhozAC4kZKC43hERASNRmOxWC4uLj4+Po0tEN33wxEKhTo6OuvWreNyuQcPHgSAsrKypqamgICArKwsqblIV5OSkhobG5csWYJa9SNJuGvXLrQvhXaMw8LCmpub//77bwBITU0l5GxubtbR0YmLixMIBIGBgUZGRkiHCwoKAIBQVBzHo6Ojhw4dWltbW1JSoqen96XoLY1Gi4uLW7lyJQC4ubnhOI6Mxs2bNyO1cXd3x3F81qxZAODt7Z2bm9u6KEJvf/31VwBYtGhRfHw80uGioiIAsLa2fvr0Kdr72bp16/r16wEgMjIyNzcXncvhcFoX+8cffwDAhAkTUlJSpKqZnp5eN4C3S5e2tZ8sAvhZ5nwrUSCyAu7cuXPnzh0AQE8sPh7FxcWE3bF06VI081hbWy9ZsgTHcR6PN3v2bAAYDPCERmtzbNLUvOHtjaYaHR2dy5cvo61Ucatn4cKFchQ7JSWFKP/w4cM4jtfV1RG2W+tcHMfRrUcSZmRkfDwJCb3FcTwmJoa4yrZt2yTkTExMRH1PT0/v0aNH6JSrV68CgEDswSExaQHA8uXLvxS9RXMgAAwdOrSkpATH8bdv306cOBElOjs75+fn4zh+8uRJlIKsVgkIvWWxWEgf0GR14cIFHMcjIyOJ5rOzsysrKysrK7O3tyfu2f79+6VKmJeXh7rj8uXLparZ/v37uwO8305u93yblpZGeFTX0dFJS0uTW3O3TU1NTWNjY1u5bBbr3d27OIPRpt5qaeFJSSKhUPw56scGxYLg8/ntz+XxeMh+/pQIBILy8nIZk7n4wyEZcDgcLpeL4/jn93dx5MiRJUuWHD58eP78+SwWS0tLSzwXPQJRVVUlUlBkkIqKivLycvEj+/bta2BgIJ7C5XLZbLb4TgOfzy8vL+/Zs6d44ILa2loOh0On05WVlYuKitCekESxOI6jDRvFNl6dZ1VXQ3Aw7fffpT/soVAgJKTOzS33+XPxZCqV2tYGG+pzOI737t2b8iV8kIDj8PQpLFsG9+9LP6BXLzh0CGbMIF9R/gR8Qe8nS40MIq6xCBQZpKSk5NGjR+LpNjY2Enrbo0cP8YcBAKCkpEQ8qCTo1asX8Ujt9evXUovFMKxnz54yhKfRaKCiIutbeZGorrY2PT1dPE1HR6ctve10dJKPCI7L0kkMIzX2k/H59XbmzJkMBoNOp3forFGjRo0aNUruwnxosRjWpuriuFH//j///HPnC//sYNjX8CniV8Hn11t1dfX3RtDqGgiFMpS2y3uuwTAQid7zaYTsCZlEfpDDp5zAcVBUlNVrKZSu3affO/Qgxe7SY1PXgdRbOYFWd22/Bvg1WJiyvWSJRF9DHbsIZEPLDxlTDY5/DXMRjsuyGsh9qU8Iqbfy472ru6+gW8s2lbv6wNR1IPVWfsiYUTFM1q5Vl6Cdy9evYGzqCvy7n0zGrf5QmpqURKI2Xa2KRDyhUMDlduEVII5TmpuVcbytCuAAPC5XyOGQqvsJ6LLd6EtD9uru61j7oVq0VRG0+iX5JJANLT9k25Bd/fkttLxVImMt8BWMTV0EUm/lx1e/9hOJ3mNTkHwqSL2VHzK7NdbVnwOhqsl4XwpVkNTeT8Lnf8/xKwFFFWlbM/Gv430p2S+EiVvRXbqyXzyk3sqJ974R9RX04/e+LyUSYWw2cDggFEK3bjiVCt26fQ0V//Ig9VZ+vPed+y5tJ8P7XizBcaysTCE5GSsqgqYm0NISmZuLBg3CtbW/zjCinxVyfSs/3tetv9iZh8VirV692sTExNbW9siRI/+m4jgIBCAQgEgUGho6ZMiQGS4uyN2PdESimkOHGn76SSEwUDEkRNHfXykgQCEx8UFS0kwXFwsLi19//bWhoUGOYr98+dLV1dXAwMDBwUHiw+a2cq9fvz5q1CgDAwN3d/eSkhIA4PP5W7ZssbCwsLCw2LFjh1Ao/BCRamtrvby8UEvGxcWJZ4WGhiJH3OLk5+c7/BcUMiIhIcHW1tbAwMDT0xP5SMRxPDg42MLCwsTEZN26dXLT2/Dw8P+/5XLi119/VVFRud+Wg4X2ce3aNVtbW+QqHgDc3d1txWCz2cgll4GBAZ1OX7lyJZvN7uSVZL4Rhb1vvq2oqNi7dy9y//eJ2bx589OnT8+cOePv779y5crUlBSsvFzh7l3FEycUY2PTt249HhAQvnfv7zt3NtTXt1UIzuHoPH+uzedTADAArLGRkpkpiog4s3y5x+zZCQkJmZmZERER8pJZJBJ9//33hoaGKSkpM2bMcHBwQN7tZeQWFRXNmDFj6dKl169f79atG3J3HhkZmZKScuLEicjIyIiIiPj4+A+RasmSJWw2++LFi1u2bPnxxx9RXKIzZ874+Pj4+fm17lp9+/bd0oK7u3t6erqBgcHff/+9evXqXbt2nTt3rqioaNOmTQBw6tSpiIiII0eOJCQkXL58WW528oYNGwwNDRcvXiyvAgk67Unn6dOnf//995YtW6qqqpCTS4FAcPr0aTMzs5EjR6JjFBQU/vjjj4CAgMmTJ/fq1evIkSM8Hi8qKqoTUoKCgoxPgvD3Pd4sLy//9ddfPT09J0yY0OGrfxjx8fFxcXEDBw4cOHDg0oULnx0+7EijKdy+jdXVAY4PFQpP9+1rUF2N9+3boKuL5+dLrQYmEEhOAkKhQnn5im++MRg6lG9gcPLkSTnOt7m5ufn5+ffu3aNSqd7e3mFhYbdv3yYcjkvNrampsbGxcXNzA4AtW7aYmZmVlZVVVlb6+fkNHjwYAOzs7KqqqjotkkAguHbt2uPHj42NjY2NjRkMRnZ2NoPB6NGjx5AhQxwcHFqfoqKigjycsVgsLy+v5OTkb7/9NjExcfXq1WPHjgWAyZMnI6fiLBZrwIABw4YNA4CBAwfKmm+3b9/OYDDQXGRlZYXckQNATU2Nh4cHnU43MDBAIcyQj4jCwkIGg1FZWSlRTnFxMYPB2L59OwBwuVwGg+Hl5QUAzs7O7u7uGzZsoNPpJiYmyB90fX29u7s7nU53dHRE3igRx44ds7GxodPpCxYsqK2tJcQ7fvy4gYEB8rUpQXBw8IoVK8TvBPI87uHhsXDhwlWrVoWGhvbo0ePIkSNUKvXEiRNRUVE2NjbHjx8Xd5lN4OzsPHfuXC8vLzqdPmrUqAcPHqD0xMREW1vbb/r1izl6VEZjvuNw3Nzc6HS6ra1tQkIC0SyBgYFz587V1taeP38+ACQkJMyaNWvUqFFOTk4oYExgYCCDwfjnn39kFP4hIC9c/7qhEwptlZTGPnigmJSElZYChwONjbTmZtOiIuV16xQTEnpibY49UtOVBIJedXU/TpmioaY2e/bsts/uMFVVVVQqldoSttPU1FS840nNrays7N8SUwLFu6isrPT3958yZcqDBw82b96clZVFuLbvBIqKiqWlpch/aFlZWVpamqmpKQA4OTktWLBAtqf+yMjI/v37jx8/HgBQQKaCgoKIiIjdu3cjZZk6dWp5efmKFStWrlyZkZEhS29fv37NZDJDQ0OnTp2KQjDcu3cPAObOnXvy5Mlx48bp6uqGhYUFBwePHj0anTJ27NjWfph4PB6TyXz9+jUACIVCJpOZn58PANnZ2adPn05ISBgwYEBpaemaNWsAYNWqVadPnx48eLCmpmZycjIq4fTp0ytWrFBRUXF0dExKSvrhhx8I8ZYtWybhq54gKCgoOzsbhW9AoDiF69evt7e3t7S09Pb2xnE8Pz/f1NQUua63sLAAAAmPc4js7Ozk5OT09PQZM2Y8fvx4/Pjx7969S01NXbRoUV1dnaura0VFhQxL+GB0dPK5c66urnV1dUuWLLlz5w5qloCAgOTkZE1NTTs7OwDo06cP8tpz9+7dJ0+eiESi6Ojoly9fIsHkRXl5+eMWmpuboSXuBlZbOzAvz7ymBpqaiLpgAIDjWFWVYmxsD2ktIwMcgFVXt3zx4oyMDB6PJ0c3PfKKM4L+TE9Pv337Ngo40iExxFuSz+erq6tTKJQXL15MmjRp8uTJY8aMaU8hLBZr27Zta9euFU/Mz8+/ceMGACAjpbKyks1mZ2Rk5OTkcLnc969vDx48GBERERsbCwCJiYm5ublpaWmOjo7x8fEXLlyYNm1aY2Pjjh07AMDQ0NDf35/awcjFqampyJ7Jz89/9+5dUlKSoaHhpUuX4uPjHR0d0THIg/H+/fu3bt3KYDDS0tKetzhG9PT0LCsrW7BgQeuStbS0jIyMxF03Ijf2np6eJ06cMDQ0jImJQZM8EVkD3WwZ5tyVK1eioqJ8fX0B4OrVq8iJ7MmTJ/fu3Tt18mQZBn1NZeWy5ct//vnnoKAgaImlgmAymXl5ech39OjRo9esWYP8FV+5ciU7O7uqqsrNzY3ocHLhwoULo1pA2x7Nzc1oQ7j369c92tqb4fMVWawOXQgHUNbS+s7R0XzgwM2bNycnJ4vkFAS4c3FGiETxOCMA4OPjc+fOnZkzZyKrsP2ItyQa7mNjY4cMGeLk5HT8+PF22hdXrlzR1tYeMWKEeKKTk9OZM2cOHz68ePFitHazt7d/9OjRrVu3vL2936+3yPQfMmQIABQVFaEdRXQNDQ2N+Pj41rtkMpBYrFKpVOQtEVkUKBiZtbU1cneKLgoA2dnZAGBnZ2dpaYnW+mgzEAAWLFjQfjelGzZsYLPZoaGhM2bMQEE9UAw+Vkt3RD8kfMES6Orq9u7dGwCGDh0KAG/evEEhVVHUot50uoy7hAFEHThgaWk5b948ACCWAE5OTkSQOIJJkyYBwPnz51NSUgAAuYCXI56enpwWBgwYQKVSX716BTiO1dUptL3zBNDhL2z5ioplffvienqAYWhMlJff3969e7PZbGKvNS8vD90aGbl9+vRBAze0+MGn0+k2NjZEbMdOrEfEW1JfX//48eNr1669ePHirl27pNqAUjl79qyrq6tCy9MyHx+fLVu2oN/ICistLc3JybG2tkZdfezYse/v8Sj4FVo8qKmpaWhoAAByMszj8dasWfP777+jI9uK1Y0uhk5BHZ1Awk8qKpzo00RMQeSR9N69e0wm8/bt22fOnEGjCUhz1CoDHx+fSZMmoRohFaVSqWZmZjk5OSjk18OHDwGAiO8iQWlpKYpwhSIaqqmp6ejoYACVFRWA42wuV0aXFAEsXLSIyWQymczk5GSi0SR84qHpiEqlurq6MpnMqKgobW1tW1vb9texEzg7O4eHh9c3NBSUlNTLMWYnhnF6995bXFzC4QgEgsjISAcHBwU5Pcs1NjbW1tbeu3dvU1PTsWPHqqqq7Ozsmpubg4KCcnJypOaOHDkyLy/v3LlzXC537969ZmZmOjo6AwYM2L9//9u3b9ls9pEjR8T94HcUtBCIioqytbXlcrlcLpewwyUg5ER/pqWloZkAYWJicuzYsaysLLRKolKp+vr61tbWx48fLyoq4vF4MTEx79fbjRs3njt3bt26dQAwfvx4S0tLbW3t6Ojo2NjYX3755eDBg2ico1KppaWl+/bta73ZjTyP3759OyQkZOPGjTKupa6uPnz4cCaTGRgYeOzYsePHj6N0tE944MABJpPp6enp4uIiETOunaioqNy9e9fd3T0kJGTr1q0AMGbMGGRjL1q06KeffsrPz3dzc5NRuK+v7+nTpwMCAhQAxgwe7P7dd8MBDixalHrkSOq1azLmWxHAxYsX09PTjx8/Pm3aNLRJCACEKYVMjFu3bqFaI1O5tLTU3d1dXn29LbZs2cLhcHT19KYvXSrU15fXc2ZcXV3lp5/6urqaDByopqZWUFAgHlrtA6FQKEePHj106JCmpuZPP/0UGRlJp9Obmpq2bt2ak5MjNdfY2Hjnzp3z5s3T0tK6dOkSemqwc+fOhoYGY2NjOp3O4XB27drVaZFevHjBZrPnz5+v1QLagESI28yEnADA4XCqqqr6iwVh9fDwGDNmzIgRI2g02tGjR5OSkhQVFf38/Pr27Ttw4EANDY2cnBwg4he0Zu7cuQBALDLnzp1bV1fH4XBSU1ORaUelUhcuXIj8/a9atYqQvnVRaEEILSbf8OHDORyOtra2trY2OsDGxgYAampqHj9+bGhoiA5mMBgAkJKSUlFRgYQBADMzs6SkJEI8JpMpVXiJS1+5coXD4VRWVhI2p7a2dmJiIofDqaioIBKHDx/+8uVLqeWg0QeJpARwzNeXt2ULf/z4N9rafwOcBbjQdpwRnEJ5unhxfy0tbQA6wI+zZlWVlmYxmahVUfk1NTUoGqqhoSGHw0E2HgCkpaXJrqC8KCkpqXn9ujk0VKShISNAWfv/8RcsaMzN5bDZ9fX1ZWVlH0NmNptdVFTU0NDQ/tz6+vri4mKJI9++fVteXv4xJOw01dXVr1+/lkisqKgoKSnhcDj/xhmR6u/C09MzPj6eyWT27t0bwzCJDaeamhpVVVXx/RI+n8/n81+9eiUx5ZqZmamrq9fV1YlEIvFdIhmUl5drampKbMZwuVwOhyN18VlcXCwRJb1Pnz79+vWTWjiHw6mvr+/Tp4/4wpjFYvH5fBS4QGppaM+8qKioprJSvaioZ3i4wrVrQNQUw6BbN2iJIisJhvFXr8YtLRvz8rorKFDodNG33+JmZnivXhJvNXO5XAqFkpeX9+DBA19fX2NjYxQj7xOB41hJidK+fYp//gkf5gIF79ev+fBh0bBh5EuOH4N2vXfROpwsALTWQCUlJSUlpZycHIlHuNra2urq6mjt2k7E9xgIWscNIXjz5g3xRhRiyJAhbemtioqKioqKRKJ4HaWWRvzWolCUTpxQOH8e+Pz/PwLH21RaAABQuHkTO3FCvbYWRCKcShVZWwvnzBFOnIjT6eJ2KarghQsXAgMDqVQq2kj/dGAYrqvL/+knUFBQOHYM43A6+U519+6CRYtElpak0n4kZM23ZWVlbDa7X79+nVtMfmWg3bL+/ftT7t/v5uaGVVZ2rE9TKP/58EBRUTRwIN/HRzhlCrR64t3Q0FBbW6uvr99WGLGPC45jNTUKf/2lFByMtXqL5v0oKAjt7XmhofiAAR9BOBIA2fPtFxdX6rPy786BUEjJz8eqqzs8EUk8txQIKAUFCpcvi4YNw1vZBWpqampqah8i7QeBYbiWlsDDA//mG6Vt2yh5eR2qLG5gIPjpJ/ybbz6egCTk90AdBMehsVE+X+RxOJTCQqy8/D0fAH4uevYUOjryg4JEw4a1f4cZ19Lir14tHDUKSBvtY0LqbQdRUMB1daHdj9TfA58PbTz0/iLo1k1ob88LCRE6O8N7X9jCMFxbm79+vXDmTOjgO3MkHYXU2w6CYaJvvxWZmsqhKAUFXE0N1NW/2O9yAQAUFUWWlrzff+evWoVrabUpKoUiMjXl+/sL587Fv47oil82pL+LDoN/841gwQKlggKsgy/rSpZDpYpsbPA+fb5ovQUACgXX0xP4+OADByqcOEFhMjEWC31PDxiGDBDhtGmC778XDRwIbWz4k8gXUm87TrduwmnTKM+fKx47JvvZjywoFNHgwcKpU7vG7IRhuIaGwMVFaGdHef4cy8/HqqpAJAIVFVxfX2Rigvfvj3/hhsPXBSav97z/txCJID8fAgLgr786uUA1NIR9+2DChK43QSG3HgLBvz8UFEBZmXxO+4kh59tOQaFA//6wZg3U18Ply9DG6+NtoqEBPj4wfnzXU1po8er4WR4sk7RA7kt1FkVFGDQIfv0VJk3q2GzTvTtMmwYzZ0Krd7ZISNoJqbcfgJISWFnBr7/ChAntDWmFYTBkCHh4gLQXOUlI2gmptx+GsjIMGQL+/jB5crtUV18fliwBa2tyQUjyIZCrlA9GSQkGD4bAQMBxuHy5zZefKBTo3RtWrIBJk8jXEkg+EHI/WU6IRFBYCAEBEB8v+XAIw0BJCczMwNMTXFygTx9ysiX5QEi9lSv19XDrFsTGwuPHgBy+dO8O+vowbhw4O4OFBdBoZHBnEjmAk8gXkQhvaMDz8vBr1/Bz5/CUFPz5c7y+HhcKP6UUx44dYzAYrUQTBQYGGhkZ6enprVq1isPh4DheXV29ePFiPT09S0vLY8eOoSNPnTrFaGHChAntvOjly5eHDh2qo6Mza9as169fo8QJEyYQRZ06dar9VZAtQ05OzqRJk3R0dBgMxq1btzot4fr16xliODs7t1/C1kht4UePHo0ZM0ZHR2fChAmPHz+WOCUrKwtVZNy4cenp6eJZzs7Oly9fRr8l5CT1tr0gJ3XTp09v7wkiES4U4gLBR5SpDQoLC2k0mrm5uUR6fHy8jo5Oenr6o0ePjIyMtm/fjuP4pEmTZs2alZeXh7xV37lzB8dxb2/vdevWPXr06NGjR627mlSQx7+YmJi8vDx3d3c0aiBngA8fPkRFlZeXt78WsmUwMjJat25dfn4+cifa0NDQOQmzs7NTW2AwGEuWLGm/hK1p3cIsFotGo0VEROTm5q5YsUJHR0ckEhHH19fX02i0nTt35ufn79q1i0ajNTQ0CASCvXv3Ojs7A8DZs2fRkRJyknrbXt69e/f7778nJSV9bkHeg0AgGDNmTGBgYGu9jYqKIibh6dOn//LLL8iN8NOnT1Givb19eHg4juMMBuPq1avI1RZx+v79+9FM4ufnJ2xlPvzxxx9Dhw5Fv4uLiwGgtLQ0JSXF3t4euUoS76/tQaoMCJFIRKPRrl+/juM48kxSU1PTOQmJ3AMHDtjb2/N4vA4JKUHrFn727NmUKVNQCnIh3NzcTBz/4MEDGo2GWga5Cs3MzBQKhYcPHz58+LCenh6htxJyknrbXiorK62trX18fAoLC62trTdu3DhlyhQAYDAYxcXFOI5XVVXNmzePRqPp6Oj4+vqi21NRUeHm5kaj0fT09Hx9fVEXdHFxcXZ2XrduHY1Gs7a2Tk9PR4PrmDFjUE/icrnr1q3T09PT09PbsGEDn89vv5whISGrVq1KS0trrbfl5eVGRkZLlixZunSpjo7Os2fPcByvra1FXRwZFNevX0cx6QjHPQcPHsRx/OzZszo6OlevXk1NTdXT04uIiJAofMuWLfPmzUO/kf/RzMzMffv2EeXY2dlVVFS0sxZSZRBn7969enp627ZtMzc3R3EnOich+rO2tpZGo6WlpbVTvLaQ2sI4jjc0NBw5csTc3NzPz0/8eIFAUF9fj36fP38eAIg/cRy3traW0FtCTlJv2wvytD59+vS8vDzUmSwtLVFYnZUrV+I4jlzvuri4WFtbA8DWrVtFIhHy/+jq6opyvby8cBw3MjICAHNzc5SLlN/c3BwAfH19cRxfsWIFKgodsG3btnYKmZOTY2Rk9O7dO6l6m5WVpaOjY25ubmdnR6PRbty4QWTl5uaamJg4OzuLRKKKigo7O7uUlBQejxceHg4AxcXFzs7Ov/32W3l5eXl5eXh4uLW1NZvNfthCTU2Nn5/fwoULiQKRebx79+7ly5fX1dUVFRWZm5uvWLGinRWRKoP4ASgqzaRJk2g02g8//NDY2Ng5CdFvf39/e3v7dsomg7Za+NWrV9OnT6fRaD///LPUURgpbUhIiHhia70l5CT1tr1I6K2enh6O41lZWQAwYcIE5At30qRJOI7X1NS4uLj88ssvT548IZbEfD4fqSuXy0U/SktLkeNLGo3W3NyMip0yZQpKtLa2zs/PR+EUaDRaO4W0trbevn3748ePDx06pKenR8wniKVLl86aNQvNrps3b7a2tkbpKACqj49PU1NT6zL19PROnz6NwlUR0Gg0VHdEfHz81q1bZ82ahU5BHvCZTKZ4ObGxsa2HknaCZCD+RFYuKr+srIxGo126dOlDJDQyMmo9pXeCtloY8ebNGwB48OCBeGJzczMKMXPkyBGJ0lrrLSEn+Uyik6AYAgYGBgDQ2NiIIv2h8Jy9evU6ffr07t27USgz5IpeUVERxYxANw8A+vbti9xK9unTR1lZGfm7FIlEKARGZmamsbHxwIEDAYDFYrUzKm9JSUlERMSkSZN8fHxKSkqcnJzEc588efLdd98h77MODg5oZXj06NHVq1enpKSEhISg0BhZWVloiiMQCAR6enpRUVHIDX95efm9e/csLS2J/vTDDz/07dsXjTLQEpWid+/ewcHBRGgPBQWFpnZ/9ihVBuL3ixcvAABZKL1797axsXnx4kXnJETN8urVq8mTJ7dTNhm0buH4+HgiaKuenp65uTkSnmDp0qW3b98uKCjw8PB4b+GEnKTedhIUdpDwwIwcL6MlIo/HW7FiRUBAAHJVS0T3QxrbVvAhAnSAubl5Xl5eXl7e7du3L1++3L179/ZIVdHClStXzM3NUTyUlJSU/fv3A8B333137NixwsJCHo938ODBMWPG8Hi8VatWxcTEjBgxggiNQaFQVq1adefOHQBITEwsKSmxs7ObMmXK/v373717h+O4r6/viRMnJC49atSoZ8+enT59msvlBgUFmZub0+n0R48e+fn5cbnchoaGsLAwtCPQHqTKkJ+fHxgYWFdXZ2ZmBgAHDhzg8XiPHz9OTU21srLqnIQAkJWVpaenhyJrfiCtW9jMzCwtLe3UqVNCofDatWvPnj0bPHhwc3NzYGDgkydPsrKyjh07du7cuT59+qD2x9t+n+I/cn64bfA/goSd7O7ujuM4CjLEYDA4HA6KKnTkyBEUsPTw4cPiiTt37gSAMWPG4C3rW1QsAJiYmOA4jnQMWdpoZt61a1d0dDRxVocQX9+uX78eWfVVVVVEUGxLS8vc3FxxSxIRExOD4/iePXsAgEaj0Wi0xMREHMc5HA4RkXTixInV1dWtL4rOAgAdHZ2MjAwcx4uKiojIN9OnT6+trW1/FVrLcOnSJQAoLCzEcTwmJobYtULr/85JiOP4r7/+ipr9w2ndwjiOE7GgAADtlqFonbGxsWg8FSc/P58obejQoeJ2sricpN62F0JvkZ0jobc4jj948AAtsWg02pIlS9Ba8eHDh8S6y97eHvU5qXqL4mujG5Ofn49UFymt+L38cFgsFnpqIhs+n19ZWSmRyOVyUQC0tuDxeK3Pqq2tbWxs7KicbclAIBAIysvLJZ7cdE5C+dK6hfl8fnl5uUB+D/PJ9xzlTHV1tZqamkSElNraWgUFhY66RK6rq1NUVJQaLILkfxxSb0lIuh7kvhQJSdeD1FsSkq4HqbckJF0PUm9JSLoepN6SkHQ9SL0lIel6kHpLQtL1+Nefo9R48yQkJF8m5HxLQtL1IPWWhKTrQeotCUnXg9RbEpKuB6m3JCRdD1JvSUi6HqTekpB0PUi9JQEWi7V69WoTExNbW1vk27E1bDbbwMCgvr5ePLGxsZHBYFy7dg39GRcXx2AwTExMVqxYQRwZGhpqZWVlYWGxcePG9vuFaw8vX750dXU1MDBwcHBIT09vT+7169dHjRplYGDg7u6OHJgQxMXFOTg4fKBIFRUVHh4eBgYGKKgKSszMzHRycjIwMHB2dmYymVJPlGjJ1nLiOB4cHGxhYWFiYrJu3bovS2/Dw8Pb6jefCysrK5XPFxj+0zTI5s2bnz59eubMGX9//5UrV966dUs89+3bt8HBwePHj6+qqkLuyAm2bdvGZDKRNl68eNHHxycgICA5ObmiomLVqlUAcOrUqcDAwPDw8BMnTty8eTMoKEheMotEou+//97Q0DAlJWXGjBkODg41NTWyc4uKimbMmLF06dLr169369Zt0aJFxPHFxcU+Pj7I7dOHsGLFiu7du1+/ft3T03PhwoUvX75ks9mTJk1ycXG5du3agAEDZsyYIdVThXhLSpXz1KlTERERR44cSUhIuHz58r/+Lr6Q96VUVFQMDQ2RL+IvBAsLi8LCQg6H81mu/mkahE6nx8XFjR8/HgBWrlypqKi4d+9eIresrOz69euNjY0///zz69evkZNKAEhNTQ0NDRUIBD/++OO0adN8fHxEIlFoaCgAPH78eNSoUQ0NDR4eHv3799+6dSsAnDx5cvv27fKqy9OnT4cNG1ZRUUGlUgHAxMRk165dM2fOlJFbU1MTGxuLfES+efPGzMzs5cuXffr0EQqFU6ZMcXBwiI+Pf/z4cadFwnF8ypQpkZGR/fr14/P56urqt27dUlNT8/PzS0pKAoAHDx6MHz++rq5OWVlZ/ESJljx48GBrOS9fvhwXF5eSkgIAc+bM+ejz7fbt2xkMRmBgoIGBgZWVVUxMDEqvqanx8PCg0+kGBgYbN27k8XijRo0CgMLCQgaDUVlZKVFOcXExg8HYvn07AHC5XAaD4eXlhRK3bNni6uqqoqLi4OCAfJ02NTVt2rTJxMTExMTkt99+EwgE1dXVDAZj48aNs2bNotPpc+fOvXfvHoPBoNPpXl5eAoEAFbV+/XpHR0c6nT5r1izCfypBYmKira0tnU53cHC4f/9+c3PzqFGjnJycRCIRAAQGBqLYcLIvBAAPHjxwcHCg0+mOjo7IcDp//jyDwQgLC7O1tVVRUfH09ESFy2gQecHlctlsNgq8AAD9+vWTqHifPn0WLFgwb9488cT6+vo1a9bs27eP8ERLpVIJIZGPu4qKirKysn79+qFEfX195GVaLlRVVVGpVGpLBHBTU1PxJpKaW1lZ2b9/f5SCvJmiU/74449BgwYRXo47DYZhly5d6tev35kzZ9zc3IYPH25tbW1qapqUlMRisWJjY3/66ad169ZJKG3rlpQq59SpU8vLy1esWLFy5cqMjIyPrrevX79mMpmhoaFTp04tKyvz9va+d+8eAMydO/fkyZPjxo3T1dUNCwsLDg4ePXo0OmXs2LHIO7E4PB6PyWQiT+JCoZDJZObn56PE4ODgkpISXV3d9PR0NN5v3Lhx3759NjY233zzzZ49e4KCgtCRYWFhTU1N2traycnJEydORJ31zz//vHr1KjogIiKiW7dudnZ2ly5dErejACA1NXXRokV1dXWzZ8/Ozs6eMGFCWVkZnU6/e/fukydPRCJRdHT0y5cv+/fvL/tCxcXF48ePf/ny5cyZM7OysiZOnFhfX19XV8dkMjdu3Ii8n8fHx1+6dEl2g3wI5eXlj1tobm4GAIWWUNpKSkrIl79s1q5d6+XlRegkADg5OSUnJ4eHh1+7dm3Dhg0AwOVyeTyeouK/78ATP+SCUCjs0aMH8aeioqK42FJz+Xw+4a8P1VcgEOTm5kZFRW3ZsqVzYoi3JAqSBgDJyclPnz7t0aMHYXjX1NRcuHChpKSkublZ3IE7SGtJqXJWVlay2eyMjIycnBwul/uJ1rcHDx6MiIiIjY0FgMTExNzc3LS0NEdHx/j4+AsXLkybNq2xsXHHjh0AYGho6O/vT4yU7UFXVzcjI+P06dMA8Pz5cw6HEx0dbWVl5e/vHxERAQD79u1DRxobGycnJyPddnJy+uuvv9CKq6CgAB1gYWFx7ty5s2fP2tjY3L17V3zr4vDhwwBw8uTJ8PDwsLAwAIiPj589ezYAXLlyJTs7u6qqys3NDfVOGRdCjbBr1661a9cuXbqUzWafOXMGXWLhwoVXrlwJCQlBR3a6Qd7LhQsXRrVQXV0NAEh7AYDL5UpMCK25cuVKSkqKjY0Nk8msq6srKCioqKhgMBhRUVEJCQnbtm3z8fEBgF69eikrKxN7UfLdlFJUVORyucSfXC4XBVuQkausrEwkIh1TVlb29PR0c3N7+fLl8+fPWSxWa4fSshFvScJOiYmJ+eeffxoaGuLi4lCKgYHBX3/9lZmZGRYW9s8//xCnS21JqXJGRUXZ29s/evTo1q1b3t7e8hwCZTB48GAAGDJkCAAUFRWhUBojRowAAA0Njfj4+A6VJrGyR5610YjF5XLRnMxkMi0tLdEBbDYbLVBNTU0xDEP+UFGIEHV1dQBAhi4A2NjYIHNl+PDhjx8/FtdbFJli0KBBAGBlZQUABQUFKLLL+fPn0QDp6uqKDpZxIeR+ecmSJUTJxcXFyKMyaiVDQ0MAaM+k12k8PT09PT2JP6lU6qtXr1BMk/z8fMJmbouysjIA+P777wGgqqrq+fPnPXv2nDRpkoWFBQoVef/+fSqVqq6urq+vT4yJr169MjY2llcVevfuzWazq6urtbS0cBzPy8tDQUNk5CorK6N1JgCgHkin00tLSw8cOHDgwAG0Xpg+fTrKaifiLVlYWEin03NycrS1tZWUlEaMGJGXl3fy5Mno6Gi0LtXV1UWLVcIXvNSW7NOnT2s5c3Jypk2bhjrn2LFjP9F8i+RDywk1NTVkDaKgxjweb82aNYRP97b6K5IYnYJUiABZRBiGoT/RxomZmRmTyWQymVevXj179qz4YIwgjhcHBQohLiTu8RhFHkBVQLEFevfuTaVSXV1dmUxmVFSUtrY2Chok+0IotkVERASTyczIyDhz5sz8+fPFK0KYrLIbRI44OzuHh4fX19dnZ2cnJCQg+/zYsWMoPkBrPDw8ilpwcHA4ePDgsmXL0JOM/Pz8pqamPXv2zJo1i0KhjB07NiYmpqCgoLq6OiwsTCJe0YdgbGysra29d+/epqamY8eOVVVV2dnZNTc3BwUF5eTkSM0dOXJkXl7euXPnuFzu3r17zczMdHR0iIqcO3fOzMysQ0orAVqLhoWFsdns4uLihIQEtL5NT08/c+aMUCi8ceNGXl7eoEGDCDmltqRUOa2trY8fP15UVMTj8WJiYj6R3m7cuPHcuXPr1q0DgPHjx1taWmpra0dHR8fGxv7yyy8HDx5EgyWVSi0tLd23b1/rMFba2toAcPv27ZCQkI0bN8q4lo6OzrBhw/Ly8s6fP3/v3r2JEyfu2bOnnYurlJSUsLCwqKio06dP6+rqDhgwgMhCsW3WrFmTnJyMBJg4cSIAIFO5tLTU3d1dQuWkgvruiRMnmEzmjh07XFxc0CggFRkNIke2bNnC4XB0dXXt7OwWL16MIkdFR0dfuHCh/YWYmZmtX79+5MiRmpqaLBZr27ZtADB37tyJEydaWlr269dPQ0Nj7dq18pKZQqEcPXr00KFDmpqaP/30U2RkJJ1Ob2pq2rp1a05OjtRcY2PjnTt3zps3T0tL69KlS1FRUfISBtGtW7e4uLjY2Fg6nW5ubj558uRFixYNHjw4ICDAzc1NVVV12rRpISEhFhYWhJxSy5Eqp5+fX9++fQcOHKihoZGTk/OJ7GR1dXW0ITl37ty5c+cqKysnJib++OOPy5cvp1KpCxcunDNnDgAsXrw4LCxs06ZNs2bNkljRqaqq+vr6hoSEbN68mTBHpc6ZAHDo0CFPT8/NmzcDgL29vUQUFnQW+p/YxENYWFiEhIRUVVVpa2vHxMSIhx1YtGjRy5cvw8PDUW/es2cPCmmLHp+AmJEs+0Ljx48PCAjYtGlTenq6trb29u3bR44c2dZGq4wGkSP6+vppaWl1dXXdu3cntnOQxUugqqoq9WHYuXPniN+//fabn59fY2OjqqoqSlFSUvrzzz8jIyMFAkFHwzW8lzFjxpSXl1dXV2toaKBxWU1NjRCydS4ArFq1ysvLq6GhoXV0NTs7uw95CIQYP358YWFhdXW1qqoqYeKtWbPG29u7trZWU1MTjezichKIt2RrOTU1Nc+fP89ms/l8voaGxkd/fuvp6RkfH89kMnv37o1hmETnq6mpUVVVFVcPPp/P5/NfvXolMcOYmZmpq6vX1dWJRCLiEaJs6uvrFRUV29nd8/Pzrays5s6dGx0dXVVVpaOjI3VQ4PF41dXVdDod3YCsrKwHDx74+voaGxu39SqMVAQCQU1NDbKZZYMaRL77yR8EjoNIBHw+iESAYaCkBAoK0MYASvKR+ETzLQBIjXPTWgOVlJSUlJRycnIknlhqa2urq6ujhXE7QVtBHYVCochQJ2VlZfGAixcuXAgMDKRSqX/88UeHrqKoqNgepYWWBulQ4R8LPh979w57/RorKsJKSzEuF1dVxQcMEJma4nQ6fCFC/m/w0efbsrIyNpvdr1+/9z5d+LzweLzi4mIqldqnT5/2n9XQ0FBbW6uvry/fh5NfHEIhVldHefRI4epVyv37WHEx1tiI5ltcVVVkZyeYP180ZgyuoQGUL+vN2a+VL+s9x68HoRD4fKy5GQQCoFDwHj1AWbmr9mkul/LsmWJCgsK5c1hZGbQ8MxMH79VLOG+eYNEi0YAB5MT7CSD1Vq4IhRiXCxUVlIIC7NUrrKoKcBxXV8dNTERmZriubhfr0yIRVlGhkJiodPgwVlAAskM3du8unDaNv2yZyMoKunf/VCL+j0LqrZxoaqK8eUPJyKDcu0f55x+svBxjs4HHAxwHCgXX0hKOGyecP1/03Xd4V4lnKxBQ8vKUAgMVLl+Gdj5DVlQUDh8uWLlSNHIkrqpKblZ9PEi9/WD4fOz1a8WzZxUSEykFBdDUJNWSBAoF79+f/9NPQldXXEPjS+/TIhH26pXy9u0KZ8/Cf7/dew8YJjIxEXh7C6dNw3v16qpLgy+er3o35RPQ2Kjw+LFiSIjC7dvQ8oqvdEQi7OVLpd9/BwwTzp6Ntzzh/ELhchVu3lS4erVjSgsAOE7Jz1cMC4MePYRTpnzp1eyykMPhB8DnU5hMxZ07FVJT36O0CBzHysqUwsMp6entOv7zgb17R3n4EDr31bFIRCkoUDhxAisqgv9++0IiL/6db7+gx/pdiPJyOH0aMjPbu/wDABzHXr3qdvIkmJmBqekXakai/af6+vdsRMlAJFLIzFQoKoJBg4DsWh+BL7LfdAlwHF6/hrS0Dk9KIhGkpcH9+/Al7ykoKMAHvlbJYkF5ObR8kkoiX0i97SwiEVRXQ3299F0o2VRUQE4O/NfH2hcEhoG6OgwaBK0+oupYIV33kfUXD9msnQXDoFs3oFA6szPM50NtLbDZnTdEPzYqKuDgALa2nd/31tWFfv3gy35JrutC6m1noVBAVxeMjDrzKgWai9rx0d9nQ0kJLCxg7VoYPfo/s247588ePWDaNPj2W1JvPxLkc6APQFcXXF3h5UsoKemYtaypCYaGoK7+RT/FpVJh7Fjo3Rvu3oWsLGCxgEaD/Hxo5an4P6DhbNIkmD8f+vYl7eSPBKm3HwCVClOmAIcDCQnw7Bk0NQGGvd/0VVaGESNgxAiQ9/eocgbDgEqFIUPg22+BwwE+HxoaIDCwTb1VVIQxY8DGBkxMYNgw6N+ffNvx40Hq7QeAYaCnB/Pmgbk5ZGdDURFkZICMr3AxDLS1wd4e3Nxg8OCuYUMqKgKV+u/eMoUia/6kUMDREebMAXV16NkTvu4PpD43ZON+GAoK0Lcv9OoFtrZQXw/h4W3qLYaBgwPMmAHW1mBqCp36NvgzI9uawDDo1Qt0db/odfvXAqm3HwyGQY8e//6T8c0AhsHIkTBnDtBooKT0Ra9s2wKt4d+rvSQfH1Jv5cd7O7S6Oqird+HpCD30kl3HL/bJ1tcFud0nV2T0WuSWqUtPR++VH8O6dgW7Dv/Ot+R3fHKAw1Hi85XannN4PJ6Aw+m68y3W1KTE58uw0HjNzV26gl0Icr6VKzJmm6/egCSN5E8IqbdyhUJps+9iWJf3Vyr73RJU8S5dwa4Dqbfyoz0TTpeekRQU3iO/jGGLRK6Qeis/cBxwXNaEIxR27elIdu0QXbqCXQdSb+UKUt3O5XYJ3lu7rl7BLgKpt3JF9md9Xf0le9k6Se5LfUK6eE/60pDddzv3sW7X4quv4JfBF6S3XC731atXgq7rSQzDqv4b00gcHMefP38urtVsNtvAwKD+v14vUBTZa9euSS0kNDTU39+f+DMuLo7BYJiYmKxYsQKVw+fzt2zZYmFhYWFhsWPHDmFHvTHKBsNkmQzopQtyypUrEnec4IvQ25KSkunTp2tpaQ0aNEhNTc3T07OhoaHTpZ05c0ZFRSUyMvLPP/9UUVGJjY2Vo6gyYL17d/zEibZyRTguEAjQdPT27dvg4ODx48dXVVVJqNa2bduYTGZTU5PE6WfOnPHx8fHz8yPCFF68eNHHxycgICA5ObmiomLVqlUAEBkZmZKScuLEicjIyIiIiPj4eHnV7tGjRwMGDHh/gW3Mt56enioqKi9fvhRPJO5UewSg0+kqKioo2ltcXJyKioqLiwvKcnFxUVFRuX//vtQTw8PDjxw50p5LfFG0vuPifH69ra2tHTlyZEpKyowZMzZt2mRvbx8fHz958uROT7zffvttQEAAEfod/1QzwIv8/EZos+PiAEJFRTQdYRimo6Pj4eEhcUxqauqLFy/GjRvX+vQePXoMGTLEwcGBSElJSfnhhx/GjBljZmbm5+eXlJQkEAgqKyv9/PwGDx5sb29vZ2dXVVUlt+oBAIBC9+5ttiaajWW2tsS9kLhTspk0aRIAoBC1t2/fBoBr164JBAKBQIDME2tra6knbtiwISQkpD2X+KJofcfF+fx6GxoaWlVVtXz58hMnTvj5+V26dGn06NFMJvPUqVM3btxgMBhBQUEMBoNOp+/cufP48eMmJiYGBgbR0dEAcP78eQaDERYWZmtrq6Ki4unp2dzcXFpampiYSAztDx8+RKevXr0azWPOzs7Ozs4od9myZQwGo7m5OTU1FRViYGCwc+dOHMcfPXo0f/78iIiIdlbEysbGJzhYuiWJYXwKpaklCGifPn0WLFiAAnkT1NfXr1mzZt++fRRpJTg5OS1YsMDS0pJIoVKpRKhRpJ8VFRX+/v5Tpkx58ODB5s2bs7KyiBlJLggAKJqabcUK51EouIrK33//zWAwDh06hG7NwYMHGQzGggUL0DF//vknun0ohjpxp6TeR4nyx4wZAwAZGRkAcPXqVZT45MmTvLw8AHBwcOjevXtUVJSFhQW6iXv27AGAUaNGAUBhYSGDwSgpKVm5ciWatx0dHXNzcwGAx+P5+/tbWFjQ6fS5c+eWlJRAG/2qo/3hA2l9xwlEItHn19snT54AANHDKBTK999/DwDZ2dl1dXVMJnPr1q1GRkZsNtvf33/ZsmW2trZVVVU+Pj7V1dXogI0bN6K4uPHx8ZcuXUKJxFQTExMzcODAPn36HDp0CJlk2dnZ2dnZKDc3N5fJZIpEojlz5hQVFe3bt2/AgAH+/v43b96srKw8e/bs33//3c6KKHTr1m3IEJGRkZQ8DCvq3p3Vu7eMd3fXrl3r5eXVr1+/dl7OyckpOTk5PDz82rVrGzZsAAAul4uy0tPTb9++raysXFdX187S2kMzwNtu3aQ6eRQBcKlUvHfvBhaLyWRWVlaiu7BmzRomk0mEDo+OjnZ0dORyub6+vg8fPiTulNT7KHEJBoMBAPfv33/16lVVVdXo0aNRTZEmjx8//u7du76+vioqKnv37qVSqb/99tuLFy/QYQAwduzYixcvHjlyxNnZedu2bWlpacuWLQOAoKCgnTt39urVa/LkycnJydOnTxcIBFLl6Wh/+HjU1NR8fr19+/Yt/DfGNGqsiooK9Ofy5ctjY2PXr18PADt37oyNjZ07dy4AFBcXowMWLlx45coVZAsVFBRIlL98+fLo6OgLFy4AwKlTp6TKwOPx2Gw2m83GcXzr1q1nz541NjY2MzMLDg6WmBVloaDQrKt7TkuLJZGOYbim5mlNTa6WVlunXrlyJSUlxcbGhslk1tXVFRQUENVvCwaDERUVlZCQsG3bNh8fHwDo1asXyvLx8blz587MmTO3b9/eXuHbARegQEVFZGEhaVNgGAeg3NAQ79dPYplgZWVVUVFBBPU+cOBARETE/v37AeD8+fMS5Uvcx4aGhmMtoNW1trZ2WlrarVu3AMDb2xsA7t27l56eDgAjR47U0tI6cODAsWPHnJ2dDQ0NAaCsrGzHjh0AYGho6O/vj0axuro6GxubM2fObNu2DQDCwsKoVOrVq1djYmI8PDzy8vLSWxzxSMjT4f7w0dDW1v78emtkZAQAr1+/JlKKiooAwMDAQPwApNj6+vrQEqWe2NEZPHgwAKBbxWsVOmDo0KEAoKenp62tLbEvQtCtW7ddu3YBgI+Pz8SJE3/99VeBQNC/f38vLy9HR8f218V78+YDtbWNPj4iU1NQVER+G0WWlvwNG1JVVQVtO6YpKysDgO+//37GjBmPHz/etWvX2bNnZV/r9evXFhYWaWlpaWlpxsbGVCpVXV3dxsbmypUr6AAGg/HPP/+0X/j3wgMo6tFD4OaGf/MNKCvjADgATqHwunXLBMi3tcW1tSW22ebMmUMV858+ZMgQALCysgKxYZdA4j6WlZWtaCE5ORnDMHQvdu/eDQAjRowYPXr0zZs3U1JSqFSqpaUlnU6/efPm6NGjBwwYcPPmzdbyu7u7W1lZXblyZerUqS4uLn///Xd9fT2bzbaxsenRowcA2NjYQEv3ay1PJ/rDRyI7O/vz6+3w4cMBICIigs/nA0B9ff3hw4cBYNiwYa0Plrq4Qo2u0IYJilYsjY2NVVVVRCz5qqqq5uZmoVCI9jkaGxttbW0zMzPPnj07evTonJycqKgogUDAYrEI4/O9ZGdnHz9xYl9SkvKyZezgYE5QEC8ggBcWxgsJEc6eXSvTx4WHh0dRCw4ODgcPHly2bFlzc3NQUFBOTo7UU9ATo/z8/Kampj179syaNYtCoQwYMGD//v1v375ls9lHjhyxt7dvp/DthKWkJJw6lb9+ff3gwa8BSgH4hoZFQ4bsAMhSUQElJQl7R+2/vu/Q8FReXg4AOjo6EoVL3Ec9Pb1TLcyfPx9alrilpaXGxsbq6ur29vZsNruqqmrcuHGKioqhoaEnT5709vZ+/vw52l0nQKM5h8PZu3fvP//8gzR/69at6OYii48QTFtbW6o8He0PcofoDxSpWyCfmCVLlmhra9++fXvkyJErV64cNmxYYWHh8OHDJ0yYIJfyg4KCEhMTkSU5duxYAPjmm28AwM/Pz9fXFx3D4/HGjh27ePFifX19JycnAFBTU7t69Wrv3r29vLzaeSFkX5lbWGiamWk7O/fx8Xk5erTAxUX03Xe4hgbW8a9Sm5qatm7dKq634sOWmZnZ+vXrR44cqampyWKxkNW3c+fOhoYGY2NjOp3O4XCQESFH7qalOXl4jI2MHPbo0RKA85MmiaKjC3/88Q7A4fj4sLAwZANLFRgA1q9fn5yc7OfnBwAjRoyQfS0qlTqpBTMzM2gZ4qFFgYmRHe3A19bWAkC/fv2qqqqSkpIAQCQSoXJKS0v37dt35MiRsWPH3rlzx9HREc2iampqrq6u+fn5W7du/euvv/bt20elUtva3+5of5AL4g1I9IeBAwd+fj81Kioqd+/e9fb2TklJQX10/vz5QUFBhMTiEylKJLLa2ttEZ6FcOzu7RYsWAYCNjc3GjRsBYN26dT/88MOBAwfMzMysrKyYTKa6unpoaOimTZuQpTRs2LAlS5Z0dAdi6dKlS5cubSv3zp07EimqqqocabGFzp07h36oqalJHCCxXv3tt9/8/PwaGxtVW8JVGhkZ3blzp6GhgUKh0OQaIJto6jt37wKAlZXVOFdXN29vkZLSsCFDvjt4MD09fePGjbNnzy4sLIQ2bo22tjbam3Bzc5sxY8bp06dB7E69F0NDQ11d3dLSUjs7O2gxawFg5MiRAODu7n7y5EmkV8OHDy8tLUXz5+LFi8PCwjZt2pSWlpaZmbl69Wp01t69e3v27BkUFMRisYKCggBAV1f30KFDGhoa7ZTnEyB+x8X7w/8Bvvb10hNQMngAAAAASUVORK5CYII=\">"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The coefficient for `pct_unemployment` is -0.15, which means **\"for every percentage point increase in unemployment, life expectancy drops 0.15 years.\"** While this might be _correct_ it unfortunately sounds very stupid.\n",
"\n",
"Income's coefficient is rather unfortunate, too. \"For every extra dollar in median income, life expectancy goes up 0.00004825 years.\" No thanks! Individual dollars don't really make sense here.\n",
"\n",
"Instead, it'd be nice to say something like **\"for every ten point increase of unemployment,\"** or **\"for every 10k increase in income.\"** If we were using the dataframes version of regression, we'd create new columns, maybe something like this:\n",
"\n",
"```python\n",
"merged['income_10k'] = merged.income / 10000\n",
"merged['unemployment_pct_10pts'] = merged.pct_unemployment / 10\n",
"```\n",
"\n",
"Since we're using the formulas method, though, we can do the division right in the regression!"
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>OLS Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>life_expectancy</td> <th> R-squared: </th> <td> 0.490</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>OLS</td> <th> Adj. R-squared: </th> <td> 0.490</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>Least Squares</td> <th> F-statistic: </th> <td> 8997.</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Thu, 12 Dec 2019</td> <th> Prob (F-statistic):</th> <td> 0.00</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>21:50:50</td> <th> Log-Likelihood: </th> <td>-1.6208e+05</td>\n",
"</tr>\n",
"<tr>\n",
" <th>No. Observations:</th> <td> 65656</td> <th> AIC: </th> <td>3.242e+05</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Df Residuals:</th> <td> 65648</td> <th> BIC: </th> <td>3.243e+05</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Df Model:</th> <td> 7</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>nonrobust</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>t</th> <th>P>|t|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> 81.2365</td> <td> 0.122</td> <td> 665.628</td> <td> 0.000</td> <td> 80.997</td> <td> 81.476</td>\n",
"</tr>\n",
"<tr>\n",
" <th>pct_black</th> <td> -0.0666</td> <td> 0.001</td> <td> -56.960</td> <td> 0.000</td> <td> -0.069</td> <td> -0.064</td>\n",
"</tr>\n",
"<tr>\n",
" <th>pct_white</th> <td> -0.0386</td> <td> 0.001</td> <td> -36.707</td> <td> 0.000</td> <td> -0.041</td> <td> -0.037</td>\n",
"</tr>\n",
"<tr>\n",
" <th>pct_hispanic</th> <td> 0.0131</td> <td> 0.001</td> <td> 10.298</td> <td> 0.000</td> <td> 0.011</td> <td> 0.016</td>\n",
"</tr>\n",
"<tr>\n",
" <th>pct_less_than_hs</th> <td> -0.0862</td> <td> 0.002</td> <td> -48.979</td> <td> 0.000</td> <td> -0.090</td> <td> -0.083</td>\n",
"</tr>\n",
"<tr>\n",
" <th>pct_under_150_poverty</th> <td> -0.0596</td> <td> 0.003</td> <td> -21.738</td> <td> 0.000</td> <td> -0.065</td> <td> -0.054</td>\n",
"</tr>\n",
"<tr>\n",
" <th>np.divide(income, 10000)</th> <td> 0.4825</td> <td> 0.006</td> <td> 83.217</td> <td> 0.000</td> <td> 0.471</td> <td> 0.494</td>\n",
"</tr>\n",
"<tr>\n",
" <th>np.divide(pct_unemployment, 10)</th> <td> -1.4900</td> <td> 0.045</td> <td> -33.408</td> <td> 0.000</td> <td> -1.577</td> <td> -1.403</td>\n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <th>Omnibus:</th> <td>2114.193</td> <th> Durbin-Watson: </th> <td> 1.520</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Prob(Omnibus):</th> <td> 0.000</td> <th> Jarque-Bera (JB): </th> <td>4788.035</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Skew:</th> <td> 0.183</td> <th> Prob(JB): </th> <td> 0.00</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Kurtosis:</th> <td> 4.271</td> <th> Cond. No. </th> <td> 792.</td>\n",
"</tr>\n",
"</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: life_expectancy R-squared: 0.490\n",
"Model: OLS Adj. R-squared: 0.490\n",
"Method: Least Squares F-statistic: 8997.\n",
"Date: Thu, 12 Dec 2019 Prob (F-statistic): 0.00\n",
"Time: 21:50:50 Log-Likelihood: -1.6208e+05\n",
"No. Observations: 65656 AIC: 3.242e+05\n",
"Df Residuals: 65648 BIC: 3.243e+05\n",
"Df Model: 7 \n",
"Covariance Type: nonrobust \n",
"===================================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"---------------------------------------------------------------------------------------------------\n",
"Intercept 81.2365 0.122 665.628 0.000 80.997 81.476\n",
"pct_black -0.0666 0.001 -56.960 0.000 -0.069 -0.064\n",
"pct_white -0.0386 0.001 -36.707 0.000 -0.041 -0.037\n",
"pct_hispanic 0.0131 0.001 10.298 0.000 0.011 0.016\n",
"pct_less_than_hs -0.0862 0.002 -48.979 0.000 -0.090 -0.083\n",
"pct_under_150_poverty -0.0596 0.003 -21.738 0.000 -0.065 -0.054\n",
"np.divide(income, 10000) 0.4825 0.006 83.217 0.000 0.471 0.494\n",
"np.divide(pct_unemployment, 10) -1.4900 0.045 -33.408 0.000 -1.577 -1.403\n",
"==============================================================================\n",
"Omnibus: 2114.193 Durbin-Watson: 1.520\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 4788.035\n",
"Skew: 0.183 Prob(JB): 0.00\n",
"Kurtosis: 4.271 Cond. No. 792.\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 75,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model = smf.ols(\"\"\"\n",
" life_expectancy ~ \n",
" pct_black \n",
" + pct_white \n",
" + pct_hispanic \n",
" + pct_less_than_hs \n",
" + pct_under_150_poverty \n",
" + np.divide(income, 10000)\n",
" + np.divide(pct_unemployment, 10)\n",
"\"\"\", data=merged)\n",
"\n",
"results = model.fit()\n",
"results.summary()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**That's much much nicer!** Unemployment now has a -1.49 coefficient, which (surprise!) leads us right back to the [original article](https://apnews.com/66ac44186b6249709501f07a7eab36da):\n",
"\n",
"> An increase of 10 percentage points in the unemployment rate in a neighborhood translated to a loss of roughly a year and a half of life expectancy, the AP found. \n",
"\n",
"Additionally, the negative coefficient for `pct_less_than_hs` is the source of this line:\n",
"\n",
"> A neighborhood where more adults failed to graduate high school had shorter predicted longevity.\n",
"\n",
"Even though we accomplished this pretty quickly in very few lines of code, the actual process at the Associated Press was somewhat more complicated. They didn't just throw a bunch of columns in a regression!\n",
"\n",
"Stats-wise, the AP did things like examined other columns that didn't meet a p-value threshold, and reasoned about [collinearity](http://www.stat.tamu.edu/~hart/652/collinear.pdf) between variables different. For domain expertise, they included a threshold for poverty that's actually _above_ the poverty line. This is due to the way that Medicaid provides health insurance to people in low-income households.\n",
"\n",
"## Review\n",
"\n",
"In this exercise, we reproduced [an article from the Associated Press](https://apnews.com/66ac44186b6249709501f07a7eab36da) that analyzed life expectancy and demographic information from the census. They used a linear regression to find the relationship between census tract qualities like unemployment, education, race, and income and how long people live.\n",
"\n",
"After running the regression once, we ran it a second time to get numbers that were **more human and easier to use in a story**, like a \"1.5 year decrease in life expectancy\" as opposed to a 0.15-year or 8-week decrease. We also used the **formula version** of a statsmodels linear regression to perform those calculations in the regression with `np.divide`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Discussion topics\n",
"\n",
"Let's say we have the sentence \"a 1 percentage point increase in unemployment translates to a 0.15 year decrease in life expectancy.\" How do we change this into something people can appreciate?\n",
"\n",
"Translate some of your coefficients into the form **\"every X percentage point change in unemployment translates to a Y change in life expectancy.\"** Do this with numbers that are meaningful, and in a way that is easily understandable to your reader.\n",
"\n",
"What is the difference between using percentage points (0 to 100) vs fractions (0 to 1) when doing a regression analysis?\n",
"\n",
"Income, education, and race are all related, which can cause problems related to [multicollinearity](https://www.statisticssolutions.com/multicollinearity/) in our analysis. How can we avoid this?\n",
"\n",
"We should consult someone who \"knows something\" before publishing this story. Where would we go to find someone who understands statistics and the subject matter?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}