{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Intro to R Notebook\n", "\n", "## Contents\n", "\n", "## Part I\n", "1. [Variables, if statements, functions, printing output](#basics)\n", "2. [Dataframes](#dataframes)\n", "3. [Exercise](#ex1)\n", "4. [An example of plotting](#plotting)\n", "5. [Logging output](#logging)\n", "\n", "## Part II\n", "1. [Regression](#regression)\n", " 1. [Robust standard errors](#robust)\n", "2. [Merging data](#merging)\n", "3. [Other data manipulation tasks](#datamanip)\n", "4. [Creating a regression table](#output)\n", "5. [Project](#project)\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part I" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Variables, if statements, functions, printing output \n", "***" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 3\n" ] } ], "source": [ "x<-3\n", "print(x)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] \"The value of x is less than 5\"\n" ] } ], "source": [ "if(x<5){\n", " print(\"The value of x is less than 5\")\n", "} else{\n", " print(\"The value of x is greater than or equal to 5\")\n", "}" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] \"42 plus five is equal to: 57\"\n" ] } ], "source": [ "add5<-function(x){\n", " if(is.numeric(x)){\n", " theanswer<-x+5\n", " return(theanswer)\n", " }\n", " else{\n", " print(\"Argument x is not a number\")\n", " return(0)\n", " }\n", "}\n", "print(paste0(\"42 plus five is equal to: \", add5(52)))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] \"Argument x is not a number\"\n", "[1] 0\n" ] } ], "source": [ "print(add5(\"sdafgdsfg\"))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 55\n" ] } ], "source": [ "running_sum<-0\n", "for (i in 1:10){\n", " running_sum = running_sum+i\n", "}\n", "print(running_sum)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 1
    2. \n", "\t
  2. 2
    3. \n", "\t
  3. 3
    4. \n", "\t
  4. 4
    5. \n", "\t
  5. 5
    6. \n", "\t
  6. 6
    7. \n", "\t
  7. 7
    8. \n", "\t
  8. 8
    9. \n", "\t
  9. 9
    10. \n", "\t
  10. 10
    11. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 1\n", "\\item 2\n", "\\item 3\n", "\\item 4\n", "\\item 5\n", "\\item 6\n", "\\item 7\n", "\\item 8\n", "\\item 9\n", "\\item 10\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 1\n", "2. 2\n", "3. 3\n", "4. 4\n", "5. 5\n", "6. 6\n", "7. 7\n", "8. 8\n", "9. 9\n", "10. 10\n", "\n", "\n" ], "text/plain": [ " [1] 1 2 3 4 5 6 7 8 9 10" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "1:10" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Data frames \n", "***\n", "### The basic useful way to manage data in R is with a so-called *data frame*. We can create one manually:" ] }, { "cell_type": "code", "execution_count": 94, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\n", "
variable1variable2
52
74
96
\n" ], "text/latex": [ "\\begin{tabular}{r|ll}\n", " variable1 & variable2\\\\\n", "\\hline\n", "\t 5 & 2\\\\\n", "\t 7 & 4\\\\\n", "\t 9 & 6\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| variable1 | variable2 |\n", "|---|---|\n", "| 5 | 2 |\n", "| 7 | 4 |\n", "| 9 | 6 |\n", "\n" ], "text/plain": [ " variable1 variable2\n", "1 5 2 \n", "2 7 4 \n", "3 9 6 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df<-data.frame(variable1 = c(5,7,9), variable2=c(2,4,6))\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### It's easy to add functions of variables as new variables. For renaming and dropping variables, see the next lesson notebook." ] }, { "cell_type": "code", "execution_count": 95, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\n", "
variable1variable2variable2_squared
5 2 4
7 4 16
9 6 36
\n" ], "text/latex": [ "\\begin{tabular}{r|lll}\n", " variable1 & variable2 & variable2\\_squared\\\\\n", "\\hline\n", "\t 5 & 2 & 4\\\\\n", "\t 7 & 4 & 16\\\\\n", "\t 9 & 6 & 36\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| variable1 | variable2 | variable2_squared |\n", "|---|---|---|\n", "| 5 | 2 | 4 |\n", "| 7 | 4 | 16 |\n", "| 9 | 6 | 36 |\n", "\n" ], "text/plain": [ " variable1 variable2 variable2_squared\n", "1 5 2 4 \n", "2 7 4 16 \n", "3 9 6 36 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df$variable2_squared <- df$variable2^2\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## To make things more interesting, let's loading in some real data.\n", "***\n", "### Download \"qcew_wages_industrybyyear_naics.csv\" from here and save it to your project `raw data/` folder.\n", "### This data gives the US average wage by industry and year, from the Quarterly Census of Employment and Wages" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "projectpath<-\"C:/Users/Len/Dropbox/Teaching/Data TA/R tutorial/sample project/\"\n", "df <- read.csv(paste0(projectpath,\"raw data/qcew_wages_industrybyyear_naics.csv\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exploring/working with a dataframe by \"*slicing*\" it, i.e. indexing it as a matrix\n", "***\n", "### We can look at the first 5 rows of the dataframe:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
naicswage
1111 14667.21
1112 12264.32
1113 10557.79
1114 14888.62
1119 13482.87
\n" ], "text/latex": [ "\\begin{tabular}{r|ll}\n", " naics & wage\\\\\n", "\\hline\n", "\t 1111 & 14667.21\\\\\n", "\t 1112 & 12264.32\\\\\n", "\t 1113 & 10557.79\\\\\n", "\t 1114 & 14888.62\\\\\n", "\t 1119 & 13482.87\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| naics | wage |\n", "|---|---|\n", "| 1111 | 14667.21 |\n", "| 1112 | 12264.32 |\n", "| 1113 | 10557.79 |\n", "| 1114 | 14888.62 |\n", "| 1119 | 13482.87 |\n", "\n" ], "text/plain": [ " naics wage \n", "1 1111 14667.21\n", "2 1112 12264.32\n", "3 1113 10557.79\n", "4 1114 14888.62\n", "5 1119 13482.87" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df[1:5,c('naics','wage')]" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 8008\n" ] } ], "source": [ "# A \"comment\" in R can be made with the \"#\" symbol\n", "# Let's ask R to print how many records are in df:\n", "print(length(df[,1]))" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 14667.2116085861
    2. \n", "\t
  2. 12264.3232139683
    3. \n", "\t
  3. 10557.7870162356
    4. \n", "\t
  4. 14888.6238310147
    5. \n", "\t
  5. 13482.8713409426
    6. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 14667.2116085861\n", "\\item 12264.3232139683\n", "\\item 10557.7870162356\n", "\\item 14888.6238310147\n", "\\item 13482.8713409426\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 14667.2116085861\n", "2. 12264.3232139683\n", "3. 10557.7870162356\n", "4. 14888.6238310147\n", "5. 13482.8713409426\n", "\n", "\n" ], "text/plain": [ "[1] 14667.21 12264.32 10557.79 14888.62 13482.87" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#First five values of wage only:\n", "df$wage[1:5]" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] \"Mean industry wage: 41390.2396774684\"\n", "[1] \"Max industry wage: 226189.599662489\"\n" ] } ], "source": [ "print(paste0(\"Mean industry wage: \", mean(df$wage)))\n", "print(paste0(\"Max industry wage: \", max(df$wage)))" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] \"Mean industry wage: $41390.24\"\n", "[1] \"Mean industry wage: $41k\"\n" ] } ], "source": [ "#What if I want to display these more pretty-like?\n", "print(paste0(\"Mean industry wage: $\", round(mean(df$wage),2)))\n", "print(paste0(\"Mean industry wage: $\", round(mean(df$wage),-3)/1000, \"k\"))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Suppose we wanted to look up the wage for the coal mining industry in year 2000. The 4-digit \"NAICS code\" for this industry is 2121. We can do this by slicing `df` with the logical vector `df$naics==2121 & df$year==1990`:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\n", "
yearnaicswage
211990 2121 39767.71
\n" ], "text/latex": [ "\\begin{tabular}{r|lll}\n", " & year & naics & wage\\\\\n", "\\hline\n", "\t21 & 1990 & 2121 & 39767.71\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| | year | naics | wage |\n", "|---|---|---|---|\n", "| 21 | 1990 | 2121 | 39767.71 |\n", "\n" ], "text/plain": [ " year naics wage \n", "21 1990 2121 39767.71" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df[df$naics==2121 & df$year==1990,]\n", "#Be careful, if you used \"&&\" here it would only apply the AND operation on the first row!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise \n", "***\n", "## What is the average wage for workers in the Logging industry in 2015?\n", "## Industry codes are listed here\n", "## Note, this dataset has records for the 4-digit version of industry codes. The longer codes (e.g. 5 digits) are more specific and the data is not dissagregated to that level in this dataset\n", "\n", "### Answer: $41,748.04" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# An example of plotting \n", "***\n", "\n", "## Suppose we wanted to plot average nominal wages for the coal industry over time:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
yearnaicswage
211990 2121 39767.71
3301991 2121 40668.72
6391992 2121 42295.11
9481993 2121 41716.56
12571994 2121 44389.39
15661995 2121 45734.90
18751996 2121 48031.03
21841997 2121 49145.38
24931998 2121 50243.96
28021999 2121 50638.49
31112000 2121 52560.33
34202001 2121 54335.28
37302002 2121 54639.26
40402003 2121 56627.55
43502004 2121 59325.00
46602005 2121 64058.03
49702006 2121 66638.51
52802007 2121 68861.89
55862008 2121 72240.77
58922009 2121 73500.63
61982010 2121 77466.31
65042011 2121 81257.74
68092012 2121 80450.46
71142013 2121 82067.37
74192014 2121 83749.85
77242015 2121 83594.60
\n" ], "text/latex": [ "\\begin{tabular}{r|lll}\n", " & year & naics & wage\\\\\n", "\\hline\n", "\t21 & 1990 & 2121 & 39767.71\\\\\n", "\t330 & 1991 & 2121 & 40668.72\\\\\n", "\t639 & 1992 & 2121 & 42295.11\\\\\n", "\t948 & 1993 & 2121 & 41716.56\\\\\n", "\t1257 & 1994 & 2121 & 44389.39\\\\\n", "\t1566 & 1995 & 2121 & 45734.90\\\\\n", "\t1875 & 1996 & 2121 & 48031.03\\\\\n", "\t2184 & 1997 & 2121 & 49145.38\\\\\n", "\t2493 & 1998 & 2121 & 50243.96\\\\\n", "\t2802 & 1999 & 2121 & 50638.49\\\\\n", "\t3111 & 2000 & 2121 & 52560.33\\\\\n", "\t3420 & 2001 & 2121 & 54335.28\\\\\n", "\t3730 & 2002 & 2121 & 54639.26\\\\\n", "\t4040 & 2003 & 2121 & 56627.55\\\\\n", "\t4350 & 2004 & 2121 & 59325.00\\\\\n", "\t4660 & 2005 & 2121 & 64058.03\\\\\n", "\t4970 & 2006 & 2121 & 66638.51\\\\\n", "\t5280 & 2007 & 2121 & 68861.89\\\\\n", "\t5586 & 2008 & 2121 & 72240.77\\\\\n", "\t5892 & 2009 & 2121 & 73500.63\\\\\n", "\t6198 & 2010 & 2121 & 77466.31\\\\\n", "\t6504 & 2011 & 2121 & 81257.74\\\\\n", "\t6809 & 2012 & 2121 & 80450.46\\\\\n", "\t7114 & 2013 & 2121 & 82067.37\\\\\n", "\t7419 & 2014 & 2121 & 83749.85\\\\\n", "\t7724 & 2015 & 2121 & 83594.60\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| | year | naics | wage |\n", "|---|---|---|---|\n", "| 21 | 1990 | 2121 | 39767.71 |\n", "| 330 | 1991 | 2121 | 40668.72 |\n", "| 639 | 1992 | 2121 | 42295.11 |\n", "| 948 | 1993 | 2121 | 41716.56 |\n", "| 1257 | 1994 | 2121 | 44389.39 |\n", "| 1566 | 1995 | 2121 | 45734.90 |\n", "| 1875 | 1996 | 2121 | 48031.03 |\n", "| 2184 | 1997 | 2121 | 49145.38 |\n", "| 2493 | 1998 | 2121 | 50243.96 |\n", "| 2802 | 1999 | 2121 | 50638.49 |\n", "| 3111 | 2000 | 2121 | 52560.33 |\n", "| 3420 | 2001 | 2121 | 54335.28 |\n", "| 3730 | 2002 | 2121 | 54639.26 |\n", "| 4040 | 2003 | 2121 | 56627.55 |\n", "| 4350 | 2004 | 2121 | 59325.00 |\n", "| 4660 | 2005 | 2121 | 64058.03 |\n", "| 4970 | 2006 | 2121 | 66638.51 |\n", "| 5280 | 2007 | 2121 | 68861.89 |\n", "| 5586 | 2008 | 2121 | 72240.77 |\n", "| 5892 | 2009 | 2121 | 73500.63 |\n", "| 6198 | 2010 | 2121 | 77466.31 |\n", "| 6504 | 2011 | 2121 | 81257.74 |\n", "| 6809 | 2012 | 2121 | 80450.46 |\n", "| 7114 | 2013 | 2121 | 82067.37 |\n", "| 7419 | 2014 | 2121 | 83749.85 |\n", "| 7724 | 2015 | 2121 | 83594.60 |\n", "\n" ], "text/plain": [ " year naics wage \n", "21 1990 2121 39767.71\n", "330 1991 2121 40668.72\n", "639 1992 2121 42295.11\n", "948 1993 2121 41716.56\n", "1257 1994 2121 44389.39\n", "1566 1995 2121 45734.90\n", "1875 1996 2121 48031.03\n", "2184 1997 2121 49145.38\n", "2493 1998 2121 50243.96\n", "2802 1999 2121 50638.49\n", "3111 2000 2121 52560.33\n", "3420 2001 2121 54335.28\n", "3730 2002 2121 54639.26\n", "4040 2003 2121 56627.55\n", "4350 2004 2121 59325.00\n", "4660 2005 2121 64058.03\n", "4970 2006 2121 66638.51\n", "5280 2007 2121 68861.89\n", "5586 2008 2121 72240.77\n", "5892 2009 2121 73500.63\n", "6198 2010 2121 77466.31\n", "6504 2011 2121 81257.74\n", "6809 2012 2121 80450.46\n", "7114 2013 2121 82067.37\n", "7419 2014 2121 83749.85\n", "7724 2015 2121 83594.60" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Step one: create a dataframe with average manufacturing wage, by year\n", "dfcoal<-df[df$naics==2121,]\n", "dfcoal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## We'll use the `ggplot` plotting library for this, which must be installed. `ggplot` has a formula based syntax that takes some getting used to." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAANlBMVEUAAAAzMzNNTU1oaGh8\nfHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enr6+vw8PD////agy6EAAAACXBIWXMA\nABJ0AAASdAHeZh94AAAgAElEQVR4nO2dAXuiyhJE56LGJBpX//+fvSKgiAwwTE9PNVP9vbfx\nZpNjlXgWQQR343A40eNyB+BwtjAUicMRGIrE4QgMReJwBIYicTgCQ5E4HIGhSByOwFAkDkdg\nKBKHIzAUicMRGIo0Os593qrn9FU5tztex37n9+CcO5yWggfk0R8ZzPVrgrxzTarrPUVz6+J2\n02Gmpmnzu/wXXG+mSmx0Suu7cDwiHbpnyqcuf1X7V/uF4BUiTT47j22o0zPerztOZpmYS9dm\nd1n6KxSJ8znjIv24qn6KXn6c+xv8wt/9n+/z/et57w7LwNIinVptju7Q3voa8X3Z3D3a1797\n2rtqsUmzCTc9xRafnnGRKtc+q37c8EXWzv20t/aTT990Il3bF3I7172kq9zoS9AFs3uuy45h\nLw8pEudtxkXybTnVq4PneuiveRL+3bemqq92xXXc3V/ynQe/+Lh5/+NYuV37quy+Avi+9V4Z\nPb78fbk76dy+eKq/XHf3tU77BL8+n+m7h+b1f+8fBrU+ve67uYOfFn5fdzarnd49fLa5U0/9\nu7j2f7FJMvq4deW+m3I/9y/NvzS9u93UUKTRGVdm793oOAxXQ6f+xlT1uv0p0v75d8fHra+B\nSOcn6SnSfVPteHxuCH23vKP7fdzx9+378XfNJlLvvu/fr+f4gP82369/5HUP7fRfE57uK98v\n10h2rlfEvV9skow+bm255h7PX48vP+93u62hSKMzLtKlXsn8DjeP6hm+irpvMn1fb9f7c/ZS\nP3/r5/r3YzfEp0jV+XY91H93/527CadqINLu4cdvvVZovnt37/raJbfvXm4+nvK3x5P+/Hii\nPnTo33f/Dv4ea6m/x+vQ3j00s+u1qVdHf+1q53Bn93+xSTL6uHXlTrXgVfNl93632xqKNDqe\nF3GXZrfd5z7u4Uu943N7//h8Xnavdt5+p1kPXOvbx2Yz61q9i/SxVeWa9cOh/c3n0//qqlut\ndP0ju1srd/++2zv4be7s2vzOYWTD5u0bjczX7r7ef/Hs+833cqfuW/3f3tZQpNHxbg1dT9/7\nkX3cw+firl1PdOuNv8ev3cZEet7u1gOHd5EO9Vrw0ruX9u+a1cTrlV1zp5dHtP3jVqvY8753\n3bP4cWevXdW9exhpU/9Hs0fwu15z9X/xc9fCUKRr75vDu93WbK6QzPh3K9Rz32D+fv/ObvDS\n7v0J9VM9nzwTInV/tX8XqXlL57Gp/ibS/fn/139l12wk/bYv5X67d5FG7ru5s9czuncPzVS9\nzu167vnH9JtFHzsbbrdBR4pU1IyK9LrZPLt68/aWzWXwhPq5r8GOv5elIg332t1Oj23176FI\n9Srp2t87/dovcH7cqiON3fenA897+GzTbHnVjp7aHQxjj9LHd7wiDX9jK7PZYnGz795x/eu9\niuttIg+fEL0dxpdqP3hpt3OvPd/zIl0/RLo9dlBXQ5HqtdHv26rx/jPtyuRxq15J9u+7/9Ku\nGr6l3NxD1+ZVevfcFGtWfv1fDBfp4243MxRpdI7dO67H3p7a136t34+NpMq1h6XdXyedXr//\n2NnQrlfmRDo0T7KfRqTr81dePz0Q6b6K6L+ye5jeBqtv7V4/fOrvbHjcwVcTsP/vRO/OPt+Q\nPbiv5t+K/i+GizRytxsZijQ693+1v+5P0cvR9Td+9q76vTbfHe63OzeHCF3vmyT18+3+WurY\n7P7+a496aPY6T4n0c6c3b7TU93S4tr/S7Jw+NnvJLm+IXfV+3MG327VrqPrWw4X+fTe7v5s7\n+Hu8l/NXvXZ/H3t70upDhNoDntpDhM7d+0z9XwwXqf/b2xqKND7dG6rtEm+eC5d9993v1zeH\nv3B4+892O6WZ86RIt4b+Xd8+v946bd8urZ/Qu8Em/sm97/Q4u26H9F+XvH/f7RuyzZb+6fXO\naO8eurl0u9eeB63uOmd7vxguUv+3tzUUyTN/9aE1u2P7POqeIKfD42MUf2/fbOb62C/+1f1T\n2z9E6Ke+eT4N3rP5fK79Pg8Rup13r5uPA3jqIH+715bS4y6dez+k1PWeyO2qtHffj0OE7ivO\n9l+F4+uwnec9vGb4MYrf58GEr19cIVL/bjc1FMnw/Kz6wNEG3w0FGIpkd+6bGsMjCybnsYlU\nb/dscIWQfSiS1Qnf1Og2kTa3xwxhKJLV2YVvapzq7Z59wMfHOYuHInE4AkOROByBoUgcjsBQ\nJA5HYCgShyMwFInDERiKxOEIDEXicASGInE4AiMt0r/ZWfAjKgyYICyTAqLDoEhCEBQGTJDC\nylAkIQgKAyZIYWUokhAEhQETpLAyFEkIgsKACVJYGYokBEFhwAQprAxFEoKgMGCCFFaGIglB\nUBgwQQorQ5GEICgMmCCFlaFIQhAUBkyQwspQJCEICgMmSGFlKJIQBIUBE6SwMhRJCILCgAlS\nWBmKJARBYcAEKaxMOpE4nCKHa6RNMGCCFFaGIglBUBgwQQorQ5GEICgMmCCFlaFIQhAUBkyQ\nwspQJCEICgMmSGFlKJIQBIUBE6SwMhRJCILCgAlSWBmKJARBYcAEKawMRRKCoDBgghRWhiIJ\nQVAYMEEKK0ORhCAoDJgghZWhSEIQFAZMkMLKUCQhCAoDJkhhZSiSEASFAROksDIUSQiCwoAJ\nUlgZiiQEQWHABCmsDEUSgqAwYIKYK+Oci2BQJCEICgMmiLUyzk2YRJHUICgMmCDGyjg3YdLU\nyqq7F4qEE4RlUkCCRBpVZnJl1d0LRcIJwjIpIKFrJPc2/2ZWVs97oUg4QVgmBSRuG2lqTfV2\nLxQJJwjLpIAsXCP5N4QokiYEhQETxFaZSU24jaQIQWHABDFVZkYT7rXTg6AwYIJYKjOvyfy9\nUCScICyTAjLPmPWIIulBUBgwQeyUmfeIIulBUBgwQcyUWeARRdKDoDBgglgps8QjiqQHQWHA\nBDFSZpFHFEkPgsKACWKjzDKPKJIeBIUBE8RCmfm3hxbnoEhCEBQGTBADZZZqRJEUISgMmCD4\nZZZ7RJH0ICgMmCDwZQI8okh6EBQGTBD0MiEeUSQ9CAoDJgh4mSCPKJIeBIUBEwS7TJhHFEkP\ngsKACQJdJtAjiqQHQWHABEEuE+oRRdKDoDBgggCXCfaIIulBUBgwQTDL1J91DfeIIulBUBgw\nQSDLLDn7wrocFEkIgsKACYJYZtH5gNbloEhCEBQGTBDEMhQpKQMmCMukgFAkmYZKEBQGTBDI\nMtxGSsmACcIyKSAvhps4nWpsDookBEFhwATBK7NKoaU5KJIQBIUBEwSuTIRHFEkPgsKACYJW\nJsYjiqQHQWHABAErE+URRdKDoDBggmCVifOIIulBUBgwQRKXWb777RZwtqDwHM+foEg4QVhm\nMSTgDaFb7OpoKsfrJygSThCWWQoJOUThFu8RRdKDoDBggsCIJOARRdKDoDBggqCI5HTKUCQh\nCAoDJojGNtIClZxSGYokBEFhwATR2Ws3o1L91xRJpqESBIUBEyTxGul1a8IlJxWEIqlBUBgw\nQbRE+udfLTmxIBRJDYLCgAmSdmfD8L/HVHJyQSiSGgSFARNEVaQxlbpvUCSZhkoQFAZMkJRl\nRl/KDTaWnv9BkWQaKkFQGDBB1EX6162WBjvGKZJMQyUICgMmSMIyE7u8m0+Tv62bKJJMQyUI\nCgMmSB6RRg56oEgyDZUgKAyYIOnKzL0HS5GSNFSCoDBgglAkimSZARMkWZm54+vwt5Gq+yz5\nSpEyMmCCZBNpuB8cTqSq/WPuK0XKyYAJkk+kBEEokhoEhQETJFWZBFfbk2BQJCEICgMmCEWK\nF+m/euZwHM6Kkd4nJj+LdzZwjYTMgAmSqEyK679KMAJE4ks7CwyYIGnKJLm2kQSDIglBUBgw\nQSgSRbLMgAmSpEyai4RJMCiSEASFAROEIvlN4pEN+AyYICnKJLrangQjRKSwQWmoBEFhwASh\nSBTJMgMmSIIy604+TJFkGipBUBgwQSgSRbLMgAkiX2bl2fApkkxDJQgKAyYIRaJIlhkwQSgS\nRbLMgAkiXmbtdY4okkxDJQgKAyYIRaJIlhkwQaTLrL7wHkWSaagEQWHABKFIFMkyAyaIcJn1\nV4KlSDINlSAoDJggFIkiWWbABJEtE3Fpcook01AJgsKACUKRKJJlBkwQ0TIRHlEkoYZKEBQG\nTBCKRJEsM2CCSJaJ8YgiCTVUgqAwYIJQJIpkmQEThCJRJMsMmCCCZaI8okhCDZUgKAyYIBSJ\nIllmwASRKxPnEUUSaqgEQWHABKFIFMkyAyaIWJlIjyiSUEMlCAoDJghFokiWGTBBpMrEekSR\nhBoqQVAYMEEoEkWyzIAJIlQm2iOKJNRQCYLCgAlCkSiSZQZMEAmGc/EeUSShhkoQFAZMEAGG\ncxImUSSZhkoQFAZMkHiGcyImUSSZhkoQFAZMEIpEkSwzYIJQJIpkmQEThNtIFMkyAyYI99pR\nJMsMmCASaySUIBRJD4LCgAlCkSiSZQZMEIGdDShBKJIiBIUBE4QiUSTLDJgg0QyHEmQRgyIJ\nQVAYMEEoEkWyzIAJEstwKEGWMSiSEASFAROEIlEkywyYIJEMhxJkIYMiCUFQGDBBKBJFssyA\nCUKRKJJlBkyQOIZDCbKUQZGEICgMmCAUiSJZZsAEiWJ0B31nD7KYQZGEICgMmCAUiSJZZsAE\niWE8P4WUO8hyBkUSgqAwYIJQJIpkmQETJILx+lisnTIUSQiCwoAJQpEokmUGTJD1jN55GuyU\noUhCEBQGTBCKRJEsM2CCUCSKZJkBE2Q1o38GLjtlKJIQBIUBE4QiUSTLDJggaxlvp4S0U4Yi\nCUFQGDBBKBJFssyACbKS8X6OYjtlKJIQBIUBE4QiUSTLDJgg6xiDk+bbKUORhCAoDJggFIki\nWWbABFnFGF7FxU4ZiiQEQWHABKFIFMkyAybIGsbHZcXslKFIQhAUBkwQisThZJpNPAe5RtoE\nAybICsbnBWPtlKFIQhAUBkwQikSRLDNggoQzRq5gbqcMRRKCoDBgglAkimSZARMkmDHikaEy\nFEkIgsKACUKRKJJlBkyQUMaYR4bKUCQhCAoDJghFokiWGTBBAhmjHhkqQ5GEICgMmCAUiSJZ\nZsAEoUgUyTIDJkgYY9wjQ2UokhAEhQEThCJRJMsMmCBBDI9HhspQJCEICgMmCEWiSJYZMEFC\nGD6PDJWhSEIQFAZMEIpEkSwzYIIEMLweGSpDkYQgKAyYIBSJIllmwARZynDO75GhMhRJCILC\ngAmykOHclEl2ylAkIQgKAybIMoZzkybZKUORhCAoDJggFIkiWWbABKFIFMkyAyZIiEgAQWIZ\nFEkIgsKACbKI4bjXjiJhMmCCLGBMKaQaRIBBkYQgKAyYIPOMWY0slaFIQhAUBkyQWcYCj+yU\noUhSEBQGTJAZxvzLOqUgQgyKJARBYcAEmWYs0shKmcdPUCScIOWUWeiRjTLNT1AknCCllFn2\nsk4hiCCDIglBUBgwQfyMxRpZKPP8CYqEE6SMMgEe4Zd5/QRFwglSQpnlL+sSBxFmUCQhCAoD\nJsg4I0gj9DJvP0GRcIJsv0ygR9hl3n+CIuEE2XqZsJd1CYMkYFAkIQgKAybIG+NxiHewRqBl\nxn+CIuEE2WyZmU8d6QVJx6BIQhAUBkyQHmPmY7B6QRIyKJIQBIUBE4QiUSTLDJggFIkiWWbA\nBPnYRkIIko5BkYQgKAyYIEORIIKkY1AkIQgKAybIm0goQdIxKJIQBIUBE6S/jYQSJCGDIglB\nUBgwQSgSRbLMgAnyYqz3CLCM9ycoEk6QjZahSBTJHgMmyJMR4RFeGf9PUCScINssQ5EokkEG\nTJCOEeMRXJmJn6BIOEE2WYYiUSSLDJggLSPKI7QyUz9BkXCCbLEMRaJIJhkwQRpGnEdgZSZ/\ngiLhBNlgGYpEkWwyYII8GJEeYZWZ/gmKhBNke2UoEkUyyoAJUjNiPYIqM/MTFAknyObKUCSK\nZJUBE+Qm4BFSmbmfoEg4QbZWhiJRJLMMmCA3AY+Aysz+BEXCCbKxMhSJItllwAS5CXiEU2b+\nJygSTpBNlZHwCKYMRVKDoDBgglAkimSZgRLEbakMRdKDoDBQglAk71SPaW/cJr5SpIwMkCBu\nS2VSrJGqx/8mv1KknAyQIBRp1iOKBM3ACOK2VGYRgyIJQVAYGEEo0rxHC0T6r54lXnK2OdI7\nrkyNoEiPQfmnQgmCwoAI4oRyQJRZxqBIQhAUBkIQJ5UDocxCRpBI1ZtNFAmRgRCEIlEk8wyA\nIE4sB0CZpQyKJARBYQAEoUjLROKRDciM/EGcXI78ZRYzgkQKGpSGShAURv4gFIkibYCRPUj3\n8YlNlFnOoEhCEBRG9iAUiSJtgZE7yPPzfFsoE8CgSEIQFEbuIBSJIm2CkTnI6wPmGygTwqBI\nQhAUBkXKw6BIQhAURt4gvTOe2C8TxKBIQhAURsYg7j6yOWAeVYqkBkFh5AviKBJF2g4jWxDn\n3k0yXSacQZGEICgMipSHQZGEICgMipSHQZGEICgMbiPlYVAkIQgKI69IwjlgHlWKpAZBYeR7\naZcgB8yjSpHUICiMXEE+Lj5hucwKBkUSgqAwMgX5vIiL4TJrGBRJCILCyBNk5GJIdsusYlAk\nIQgKI0uQsYuKmS2zjkGRhCAoDIqUh0GRhCAojBxBRq9yabXMSgZFEoKgMDIEGb9arNEyaxkU\nSQiCwtAP4rnqss0yqxkUSQiCwlAP4rt6ucky6xkUSQiCwtAO4vPIZJkIBkUSgqAwKFIeBkUS\ngqAwlIN4PbJYJoZBkYQgKAzdIH6PDJaJYlAkIQgKQzXIhEf2ysQxJkT6OTh32/9RJFMMzSBT\nHpkrE8nwinTd1R/Uujl3pkiWGBQpD8Mr0pc73i26/bo9RbLEUAwy6ZG1MrEMr0h3iZ7/p0h2\nGHpBpj0yViaaQZGEICgMtSAzHtkqE8/witS+tDu6L4pkiaEVZM4jU2UEGF6RrlVznrLqQpEs\nMZSCzHpkqYwEwyvS7fa9c253vK70iCLlYVCkPIwJkSIHpaESBIWhE2TeI0NlRBgUSQiCwlAJ\nssAjO2VkGF6R3HOqr1WbSSgNlSAoDI0gSzwyU0aIsUCk+6wxCaWhEgSFoRBkkUdWykgxvCKd\n3P6uz2XvflfuAkdpqARBYSQOMjjBd+IcMI9qhEg71+yvc7uVb8qiNFSCoDDSBnEUyfcTPpE6\neVYf3YDSUAmCwkgaZHgNpNQ5YB7VCJH23Uu7/e1cr5Uokg0GRcrD8Ip0eR3Z4NwPRbLCoEh5\nGF6RbtfnkQ3uuMIjipSHwW2kPAy/SLGD0lAJgsJIL5JiDphHlSKpQVAYaV/aKeeAeVRjRDp2\n78ZSJEsMipSH4RXp6RFFMsVIGSTEI/gywgyvSJX727vLdc+Tn9hiJAwS5BF6GWmGV6T7mujb\nnW5XnvzEFoMi5WFMiXSq3z/iSztbjHRBwjwCLyPO8Ip0cL8Xt7udKZItRrIggR5hl5FneEWq\nDdrX+xp48hNTDIqUh+EV6Xba1acSWndUA0XKxkgVJNQj6DIJGF6RTmsFokhZGYmCBHuEXCYF\nwyvS/TVdnEsoDZUgKAyKlIfhFelQbx8dflefjYsi5WGkCRLuEXCZJAyvSLfbX31dF7f/pUiW\nGEmCrPAIt0waxoRI9Rx5iJAxBkXKw5gS6XysnNt9UyRLjBRB1ngEWyYRwyvSqbaI20jmGAmC\nrPIItUwqhlek+2u69ef9pkjZGPJB1nkEWiYZwyvSuVkjrb0WBUXKxKBIeRhekTqXdmsPbUBp\nqARBYYgHWekRZpl0jCmR6hOgcK+dMYZ0kLUeQZZJyJgQ6e+7Pmh19fENKA2VICgMipSH4RXp\nq4qxiLORkT4lTgHDY+02wZANsnp9hFgmKcMrEo/+tsmgSHkYXpGiB6WhEgSFIRokwiO8MmkZ\nFEkIgsKQDBLjEVyZxAyKJARBYVCkPAyKJARBYQgGifIIrUxqBkUSgqAw5ILEeQRWJjljXKS3\nKzFTJEsMipSHQZGEICgMsSCRHmGVSc8YF0liUBoqQVAYUkFiPYIqo8CgSEIQFAZFysPwi8Tr\nI5lkSEACrsuXNAfOoxohEq+PZJMhAAm4UGzSHEKQzCLx+kg2GfGQkEuXp8whBcksEq+PZJNB\nkfIwpkTi9ZEMMihSHoZXJF4fySaD20h5GF6ReH0kmwyhvXYQOXAe1Zjd37w+kkmGiEhbKpNd\npNhBaagEQWFIvLTbUhmKJNVQCYLCoEh5GH6R6rND8g1ZcwyBvXZbKqPG8IrEIxtsMihSHoZX\npKp+EylmUBoqQVAYFCkPwyvS6jURRcrKiIY4oSAQZfQYXpEOLuqiLhQpE4Mi5WF4RbpU+/WX\ndKFI2RgUKQ/DKxI/am6TEQtxUkEQyigyKJIQBIVBkfIwvCJFD0pDJQgKIxLixIIAlNFkUCQh\nCAqDIuVhjIt0fz3Hl3Y2GRQpD4MiCUFQGHEQJxckfxlVxrhIEoPSUAmCwqBIeRgUSQiCwqBI\neRgUSQiCwoiCdB+M3UQZXYZXpOsXt5EsMihSHoZXpAN3NphkxECeZ2rYQhllhlck535XGkSR\ncjIoUh6GV6QdP0ZhkkGR8jC8Il12x7jPUaA0VIKgMCIgr3NwbaCMNsMr0u2X20gWGRQpD8Mr\nEnc22GRQpDwMr0jc2WCTsR7SO7uq/TLqDK9IB+5sMMmgSHkYXpFuhy9+1NwggyLlYXhF4tHf\nNhmrIf3z5psvo8+gSEIQFAZFysPwihQ9KA2VICiMtZC3C7lYL5OBQZGEICgMipSH4Rfpetw5\nF3F4A0pDJQgKgyLlYXhFurQXo6jW7rtDaagEQWGshLxfos94mRwMr0hfrj7T6mXPS1/aYlCk\nPAyvSN3eOu61s8WgSHkYFEkIgsJYBxlcfNl2mSwMr0h8aWeTQZHyMLwicWeDTcYqyMAj22Xy\nMLwicfe3TQZFysPwixQ7KA2VICgMipSHQZGEICiMNZChR6bLZGL4RTpWPGjVIIMi5WF4RTry\n6G+TDIqUh+EVqXI/Kw2iSDkZKyAfHlkuk4vhFWn1mogiZWVQpDwMr0gHF3daO4qUh0GR8jC8\nIl2qPc/ZYJARDvn0yHCZbAyvSPyouU0GRcrDoEhCEBRGMGTEI7tl8jG8IkUPSkMlCAqDIuVh\nUCQhCAqDIuVhUCQhCAojFDLmkdkyGRkUSQiCwqBIeRghIlX3WfKVImVkUKQ8jACRqvaPua8U\nKScjEDLqkdUyORkUSQiCwqBIeRgTIv0cnLvt/95FWiAURcrJCIOMe2S0TFaGV6Tr7vFmrHPn\nl0jNNtC8SP/VM7Fu46CM9B4mzshZhI71EeC/bt99p7OIayRkBtdIeRhekepDg7r/tyLNCESR\nEBhBEI9HNsvkZVAkIQgKgyLlYcy9tDu+ThBJkSwwKFIehlek68cJIimSBUYIxOeRyTKZGV6R\nbrfv4QkieWSDAQZFysOYEClyUBoqQVAYARCvRxbL5GZQJCEICmMx5P6yPWmQwpaMV6TXJ2Sr\nr1Unb0BpqARBYSyFPBZtyiCFLZkFIt1njUkoDZUgKIyFkHbBJgxS2JLxinTqro/0298FTpHQ\nGRQpD8Mr0q49r53brTxZJEpDJQgKgyLlYXhF6l/6kiLZYXAbKQ/DK9K+e2m3v53rtRJFssEI\nESlpkMKWjFek3qUv3arz6aM0VIKgMJa+tEsepLAl4xXpdn0e2eCOKzyiSHkYFCkPwy9S7KA0\nVIKgMJZBpj0yVgaCQZGEICiMRZAZj2yVwWD4RXpeso8iWWJQpDwMr0i89KVNxhLInEemyoAw\nvCJV7m/vLtf96+QnFMkCYwFk1iNLZVAYXpHua6Jvd7pdXyc/oUgWGPOQeY8MlYFhTIl0qt8/\n4ks7WwyKlIfhFengfi9udztTJFuMWcgCj+yUwWF4RaoN2tf7GtYc+U2RsjHmIEs8MlMGiOEV\n6Xba1acSWndUA0XKxqBIeRh+kWIHpaESBIUxA1nkkZUySAyvSPu1L+koUlbGNGSZR0bKQDG8\nIlWxayiUhkoQFAZFysPwivS3P6465wlFysuYhCz0yEYZLIZXJMdDhEwypiBLPTJRBoxBkYQg\nKAyKlIfhFSl6UBoqQVAYE5DFHlkog8agSEIQFIYfstwjA2XgGBMiDa8hS5EsMChSHoZXpM9r\nyFIkCwwvJMAj/DJ4DK9In9eQpUgWGD5IiEfwZQAZXpE+L31JkSwwKFIeBkUSgqAwPJAgj9DL\nIDK8In1eQ5YiWWCMQ8I8Ai8DyfCK9HkNWYpkgUGR8jC8Io1cQ5YiGWCMQgI9wi6DyfCKtPr9\nI4qUlTEGCfUIugwowyuS250okkHGCCTYI+QyqAyvSPfXddX36td1FCkXgyLlYXhFul2OlXOH\ntcc1UKRMjE9IuEfAZWAZfpHucz46t/ulSJYYH5AVHuGWwWVMinRfLfHzSMYYFCkPY1Kk89d9\njbTman0UKRvjDTJ5fcvEQQpbMn6RHttIX9xGMsboQyYvuJw6SGFLxitS/W7sD/famWP0IO2p\nAjIFKWzJeEVyB76PZJFBkfIwvCLFrIwoUj4GRcrD8IrUzvlYUSRLDG4j5WFMinT6qpyjSKYY\n/TXSv5UeIZZBZ/hFOn09LuqyeksJpaESBIXxJlLOIIUtGY9IjUXORWwpoTRUgqAwXpD1HgGW\ngWeMi9Sui1Yf1UCRsjGekAiP8MrgM3wiHa63iPM1UKRsjA4S4xFcGQMMn0hcIxlltJAoj9DK\nWGCMi8RtJLOMBhLnEVgZEwyPSE+X1h9sh9JQCYLCeEAiPcIqY4PhF+nG95EsMmpIrEdQZYww\nJkW68cgGcwyKlIcxJ9L6QWmoBEFh3CHRHiGVscKgSEIQFMa/W7xHQGXMMCiSEASFIbA+AiqD\nEoQiqUFQGA4lSGFLhiIJQUAYDiVIaUuGIglBMBgOJUhxS4YiCUEgGA4liBADJghFUoNAMChS\nNgZFEkpPEEwAABE0SURBVIIgMBxKECkGTBCKpAYBYDiUIGIMmCAUSQ2Sn+FQgsgxYIJQJDVI\ndkb3Rmz2IIIMmCA5ReKoDhcfynCNZJnxOjBoA2VkIdbXSCgNlSB5Gb0D7OyXEYZQJJmGSpB8\njPoskBQpN4MiCUGyMYbnJTZdJgWEIsk0VILkYnycKd9ymSQQiiTTUAlCkRIwYIJQJDVIBobr\nTdYgyRgwQSiSGiQxY3BdiZ483EaCYFAkIUhaxoQ4Q8fwyyhDKJJMQyVIUsb4i7gMQZQZMEEo\nkhpEQ6T8QZQZMEEokhqEIiVgwAShSGoQtW2kvEF0GTBBKJIaJO0aKeBqsPBltCEUSaahEiQl\nI+icj+hl1CEUSaahEiQhI+zcqeBl9CEUSaahEiQdI/AcxNhlMkAokkxDJUgyRui5vKHL5IBQ\nJJmGSpBUjOBz4iOXyQKhSDINlSCJGOHXlgAukwdCkWQaKkHSMFZcowW3TCYIRZJpqARJwlhz\nrSPYMrkgFEmmoRIkBWPVNcNQy2SDUCSZhkoQecbiYxlSB8nHgAlCkdQg4oy1l7CELJMTQpFk\nGipBpBmrLwWLWCYrhCLJNFSCCDPWX1IZsExeCEWSaagEkWVEXJocr0xmCEWSaagEEWVEeIRX\nJjeEIsk0VIJIMmI8giuTHUKRZBoqQQQZUR6hlckPoUgyDZUgcow4j8DKAEAokkxDJYgYI9Ij\nrDIIEIok01AJIsWI9QiqDASEIsk0VIIIMAYXOsoXBIYBE4QiqUHiGQGn3EobBIcBE4QiqUGi\nGSEngUwaBIgBE4QiqUEoUgIGTBCKpAahSAkYMEEokhqE20gJGDBBKJIaJH6NFHBa4qRBgBgw\nQSiSGiSW4YRyQJSRYsAEoUhqEBQGTJDCylAkIUgkw0nlQCgjxoAJQpHUIHEMJ5YDoIwcAyYI\nRVKDRDGcXI78ZQQZMEEokhokhtHtrNtEGUkGTBCKpAaJYDx3em+hjCgDJghFUoOsZ7zePNpA\nGVkGTBCKpAZZzei9CWu/jDADJghFUoOsZfQPZjBfRpoBE4QiqUFWMt4OCrJeRpwBE4QiqUFQ\nGDBBCitDkYQg6xjvR6kaLyPPgAlCkdQgqxiDo71tl0nAgAlCkdQgaxjDT02YLpOCAROEIqlB\nVjA+Pn1kuUwSBkwQiqQGCWd8forPcJk0DJggFEkNEswY+TSs3TKJGDBBKJIaJJQx9qlys2VS\nMWCCUCQ1SCBj9OwMVsskY8AEoUhqkDDG+FlOjJZJx4AJQpHUICgMmCCFlaFIQpAghue0WzbL\nJGTABKFIapCljKlLTpgrk5oBE4QiqUEWMiZPp2qtTHIGTBCKpAZZxpg+wbexMukZMEEokhqE\nIiVgwAShSGqQBQz3nIQ5YB7VwspQJCHINOOpD7eR1CEUSaahEuSN8SaLG/5X0hwwj2phZSiS\nEKTPWLb+SZMD5lEtrEyISFU97dfbxNfSRZrdFkqZA+ZRLaxMkEi9L5X/K0VaexFLxDJZGTBB\nKJIahCIlYMAEERWp6n+lSF7G+qsqA5bJy4AJIitSt4k0K9J/9czhtjr143j3KHcMTq5Ztkaa\nEIhrpH/16ihzDphHtbAyASJ1NlEkLyPysuRYZQAYMEEokhqkZsRpBFYGgQETRH5nA0XyMiJX\nR0I5YB7VwsoEirRsZ0OJIsVrBFQGhQETRP7IhiVfSxTJbakMDAMmiPw20vJBaagBua+OtlMG\niAEThCKpQJxQEIgySAyYIBRJAdJsHW2kDBYDJghFSg7pdjJsogwaAyYIRUoHaQ6ne+6rs10G\nlAEThCIlgwwPTTVdBpUBE4QipYJ8HOJtuQwsAyYIRUoFoUgaDJggFCkVhCJpMGCCUKRUkI8P\n71kuA8uACUKR0kBqhQYfgrVbBpgBE4QipYCMHp9qtQw0AyYIRUoASXe1vcKee0oQiiTTUBji\n+7iEyTLoDJggFEkakvS83YU995QgFEmmoSQk8Xm7C3vuKUEokkxDOcjkh2CtlTHBgAlCkQQh\n058lN1bGBgMmCEUSg8ydk8FUGSsMmCAUSQoye2oTS2XMMGCCUCQZyIJTBNkpY4gBE4QiRUJc\nczCQUpDCnntKEIok0zAGEnBdCfwyBhkwQShSFCTkAi3wZSwyYIJQpCgIRcrMgAlCkaIgFCkz\nAyYIRYqDBFx6D7+MQQZMEIoUBwm4hCV+GYMMmCAUKQoScoEJ+DIWGTBBKFIMJOhCLehlTDJg\nglCkCEjYBY/Ay9hkwAShSBEQipSdAROEIq2HBF6BD7uMUQZMEIq0GhJ6JUvoMlYZMEEo0lpI\n8BVhkcuYZcAEoUhrIRQJgQEThCKthIRfohy4jF0GTBCKtA4S7hFwGcMMmCAUaRVkhUe4ZSwz\nYIJQpDWQNR7BljHNgAlCkdZAKBIKAyYIRVoBWeURahnbDJggFCkcss4j0DLGGTBBKFIwZKVH\nmGWsM2CCUKRgCEUCYsAEoUihkLUeQZYxz4AJQpECIas9QixjnwEThCKFQdZ7BFhmAwyYIBQp\nCBLhEV6ZLTBgglCkIAhFAmPABKFIIZAYj+DKbIIBE4QiBUCiPEIrsw0GTBCKtBwS5xFYmY0w\nYIJQpOUQioTHgAlCkRZDIj3CKrMVBkwQirQUEusRVJnNMGCCUKQlE3Cq/LRBCnvuKUEokkzD\n2Qm4dkvaIKU995QgFEmm4dyEXE0saZDinntKEIok03BuKBIuAyYIRZofioTLgAlCkeZHxiOQ\nMkIQFAZMEIo0O+6hEkAQIQZMkMLKFC5SoxBAEDEGTJDCypQtkkMJIseACVJYmaJF6l7RZQ8i\nyIAJUliZgkV6bRnZWVxKEBQGTBCK5J/eDgY7i0sJgsKACUKRvNPfUWdncSlBUBgwQSiSZ953\neNtZXEoQFAZMkJwiIU9JXTnKU9Aaafj+q51/95QgKAyYIHxpNzKfxzHYWVxKEBQGTBCK9Dkj\nhwPZWVxKEBQGTBCK9DFjh9XZWVxKEBQGTBCKNJjxw1PtLC4lCAoDJghFeh/PUd52FpcSBIUB\nE4QidfP4qITv0xJ2FpcSBIUBE4QitTP96T07i0sJgsKACUKRmpn5OLmdxaUEQWHABKFIzVAk\nmwyYIBSpGYpkkwEThCK1M32CEzuLSwmCwoAJQpHamT5RkJ3FpQRBYcAEoUjNzJwlyM7iUoKg\nMGCCUKRmKJJNBkwQivSYudPW2VlcShAUBkwQilTP7Okf7SwuJQgKAyYIRaqHIlllwAShSP+W\nXIzPzuJSgqAwYIJQpEUXtbSzuJQgKAyYIBSJIllmwAShSIuusmxncSlBUBgwQSjSogu22Flc\nShAUBkwQikSRLDNgghQv0rIriNlZXEoQFAZMkNJFWnglPjuLSwmCwoAJQpEEILYYMEEKK7Np\nkZZeGtbO4lKCoDBggpQt0uJLLNtZXEoQFAZMEIoUDbHGgAlSWJkNi7TYI0OLSwmCwoAJQpFi\nIeYYMEEKK7NdkZZ7ZGhxKUFQGDBBChYpwCNDi0sJgsKACUKR4iAGGTBBCiuzVZFCPDK0uJQg\nKAyYIMWKFOSRocWlBEFhwAShSDEQkwyYIIWVsSfS+zlTxxlhHhlaXEoQFAZMkA2KNDiL9ygj\n0CNDi0sJgsKACbI9kYbXlaBIKSAoDJggZYoU6pGhxaUEQWHABNmsSE9XRhjBHhlaXEoQFAZM\nkO2J1G0jdS5RpBQQFAZMkO2J5F577R43PhnhHhlaXEoQFAZMkK2JNLxa2Mjlw1Z4ZGhxKUFQ\nGDBBNibSiCS3oUsUaUMMmCDbEmnMkZrRf7G3xiNDi0sJgsKACbIpkUYdaRndDgiKBBOksDJ2\nRBpX5MUY7hgPGDuLSwmCwoAJsiGRPIL0GBRJDILCgAmyHZF8flCkFBAUBkyQzYjk1aPP4DaS\nFASFARNkIyJN6PHGWOmRocWlBEFhwATZhkhTdhS2uJQgKAyYIJsQaXItU9jiUoKgMGCCbEGk\n6VdrhS0uJQgKAybIBkSa2eopbHEpQVAYMEHsizS396CwxaUEQWHABDEv0uxeuMIWlxIEhQET\nxLpI83uzC1tcShAUBkwQ4yIteFeosMWlBEFhwAQxLdKid1cLW1xKEBQGTBCrItWHKCw7SKGw\nxaUEQWHABDEqUsBBc4UtLiUICgMmiE2RQg7jLmxxKUFQGDBBKJIaBIUBE6SwMhRJCILCgAlS\nWBlAkYI+oFfY4lKCoDBgghgU6WHQ8g8WFba4lCAoDJgg5kQK/mReYYtLCYLCgAliTKQVH3At\nbHEpQVAYMEHERaoef9xn6utKkbKdt8TQ4lKCoDBggkiL9BClao3yfV0nUsbTLRhaXEoQFAZM\nEGGRqlsqkVZqVNziUoKgMGCCyIrUyiIv0mqNiltcShAUBkyQTCL9V88s7nnH0rs5OJyMM/t0\nrm6yayTXvmWU+J8KJQgKAyZIYWWWi/T0REik9WcYDmuoBEFhwAQprEyASM1IiSTjUWmLSwmC\nwoAJkuR9JIqEzIAJUlgZiiQEQWHABCmszAqRhI5sEPGotMWlBEFhwASBPtZOwqPSFpcSBIUB\nEwRaJJhHCSYIy6SAUCSZhkoQFAZMkMLKUCQhCAoDJkhhZSiSEASFAROksDIUSQiCwoAJUlgZ\niiQEQWHABCmsDEUSgqAwYIIUVoYiCUFQGDBBCitDkYQgKAyYIIWVoUhCEBQGTJDCylAkIQgK\nAyZIYWUokhAEhQETpLAyFEkIgsKACVJYGYokBEFhwAQprAxFEoKgMGCCFFaGIglBUBgwQQor\nQ5GEICgMmCCFlaFIQhAUBkyQwspQJCEICgMmSGFlKJIQBIUBE6SwMhRJCILCgAlSWBmKJARB\nYcAEKawMRRKCoDBgghRWhiIJQVAYMEEKK0ORhCAoDJgghZWhSEIQFAZMkMLKUCQhCAoDJkhh\nZSiSEASFAROksDIUSQiCwoAJUlgZiiQEQWHABCmsTDqR5mf5hc8TD0oQlBwMMpyQHBQp+6Dk\nYJDhUKRFgxIEJQeDDIciLRqUICg5GGQ42CJxOBscisThCAxF4nAEhiJxOAJDkTgcgaFIHI7A\n6IhUNX/eZ+yr4viDaCbxPRDqD8hUENUkn4/A+5JCyDEdREWkNlz7x/Cr4owHUQ7hfSDUHxBv\ngCz/vr09Au9LKneOBRk0RKpub49KPpE8QSjSMEB2kQZLKncOEJFu749KxjXSRBDtQRBpMoj6\nvD8CWUSazDE9FOmmv0UwvP/B1/xBcjwg0CLNPyDKIlXV+9dsIiEEgRAJMsjbzYxLJiCH8l67\n9slbZV0jDYP0/kY1CNzz9y0ARXrlGNwaG02R+reyijQIMvgbtRwIIo0GufW/5ggy8t9Zcwxu\njQ23kXIEeX+uZHxAQIPcBkvKQA6K1FuZq6ZAeEAmguTwKL9IEzlmYvDIBv0gVbcTKPcDghtk\nuKTwc/BYOw5HYCgShyMwFInDERiKxOEIDEXicASGInE4AkOROByBoUgcjsBQJA5HYCiSwfly\nh+bGwX3lTcLphiJZnMr91l9+nfYhVhzfUCSLc3bV9Xa7Vu6cOwmnHYpkch4v7toXdtcv576u\n9a3zwbnqeL/h3F+1zxqwuKFINqdyf3/tC7vK3Wd3v3FyjznWIu259aQ7FMnmnN3h0Lyw+67N\nObqf221Xbzn9OVeLdMwdsLShSEanfj33uLF7LMJmP97l9L1vRLpkjFbkUCSr41z3tZn7zX13\ny3Gxag8fcavzKdKX2/2cLhQpy/ARtzqdLDv39p0rRcoyfMStTifLsd6x8Ov29XfOt+ueImUZ\nPuJWp5Pl+tj97f5qpbiNlG34iFudpyyXL+f2jz3hjxsUKcvwEedwBIYicTgCQ5E4HIGhSByO\nwFAkDkdgKBKHIzAUicMRGIrE4QgMReJwBIYicTgCQ5E4HIGhSByOwPwP3C8Zb7A0URIAAAAA\nSUVORK5CYII=", "text/plain": [ "plot without title" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#install.packages(\"ggplot2\")\n", "library(ggplot2)\n", "\n", "#ggplot titles by default show up on the left. I like them in the middle:\n", "theme_update(plot.title = element_text(hjust = 0.5))\n", "\n", "ggplot(data=dfcoal, aes(x=year, y=wage, group=1)) +\n", " geom_line()+\n", " geom_point()+\n", " labs(title=\"U.S. Coal Industry Wages Over Time\", x=\"Year\", y = \"Average nominal wage\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Follow up question (\"extra credit\"): have these nominal wages kept up with inflation? You could find out by integrating data from here ." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## We can also save this plot to a file (e.g. to put it in our paper). Remember the project folder!" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "ggsave(filename = paste0(projectpath,'figures/coal_wagesbyyear.png'), width = 10, height = 5, dpi=500)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Logging output \n", "***\n", "\n", "## If you have experience with Stata, you might be used to creating a \"log file\" from your code. We can do a similar thing in R, using the `sink` command:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "projectpath<-\"C:/Users/Len/Dropbox/Teaching/Data TA/R tutorial/sample project/\"\n", "\n", "#split=TRUE allows output to also go to the screen (otherwise you would only see it in the log file)\n", "#append=FALSE tells R to overwrite \"mylog.txt\" rather than appending to itself, if it already exists\n", "sink(file=paste0(projectpath,'/log files/mylog.txt'),split=TRUE, append=FALSE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## This doesn't work in a Jupyter notebook, so copy this code into R Studio\n", "\n", "### Note: it's good practice to close the connection to this file at the end of your code (this makes more sense when working in RStudio in which your have a \".R\" script file). The following code will close any open file connections:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "sink.reset <- function(){\n", " for(i in seq_len(sink.number())){\n", " sink(NULL)\n", " }\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part II" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Regression analysis \n", "***\n", "## Let's move on to regression analysis. For this, download this data set to your `raw data/` folder. Or, you can read it right of the web (and save a local copy):" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "df<- read.csv(\"http://www.columbia.edu/~ltg2111/teaching/wages_bycounty_yr2000.csv\")\n", "write.csv(df, paste0(projectpath,\"raw data/wages_bycounty_yr2000.csv\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## This data set contains the average wage and employment count by US county (`fips') in the year 2000:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
Xfipswageemployment
1 1001 24803.96 8844
2 1003 21926.4540999
3 1005 22849.07 9890
4 1007 20039.41 2686
5 1009 21934.78 7372
\n" ], "text/latex": [ "\\begin{tabular}{r|llll}\n", " X & fips & wage & employment\\\\\n", "\\hline\n", "\t 1 & 1001 & 24803.96 & 8844 \\\\\n", "\t 2 & 1003 & 21926.45 & 40999 \\\\\n", "\t 3 & 1005 & 22849.07 & 9890 \\\\\n", "\t 4 & 1007 & 20039.41 & 2686 \\\\\n", "\t 5 & 1009 & 21934.78 & 7372 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| X | fips | wage | employment |\n", "|---|---|---|---|\n", "| 1 | 1001 | 24803.96 | 8844 |\n", "| 2 | 1003 | 21926.45 | 40999 |\n", "| 3 | 1005 | 22849.07 | 9890 |\n", "| 4 | 1007 | 20039.41 | 2686 |\n", "| 5 | 1009 | 21934.78 | 7372 |\n", "\n" ], "text/plain": [ " X fips wage employment\n", "1 1 1001 24803.96 8844 \n", "2 2 1003 21926.45 40999 \n", "3 3 1005 22849.07 9890 \n", "4 4 1007 20039.41 2686 \n", "5 5 1009 21934.78 7372 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df <- read.csv(paste0(projectpath,\"raw data/wages_bycounty_yr2000.csv\"))\n", "df[1:5,]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Suppose we want to run a regression of wages on employment in this data. If we just wanted to visualize the regression line, we could let `ggplot` do it all automatically:" ] }, { "cell_type": "code", "execution_count": 131, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAPFBMVEUAAAAA//8zMzMzZv9N\nTU1oaGh8fHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enr6+vw8PD///+/r5BdAAAA\nCXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2diZabuhZE8Yu7M9xMjv//X59ng5BAI6ri\n1FkrPdh4d3FgB5AxDGeVSlVcQ+8AKtUeSiKpVBVKIqlUFUoiqVQVSiKpVBVKIqlUFUoiqVQV\nSiKpVBWqlUin1YqYpKSE78g3g5dItvHk8XHwEsk2njw+Dl4i2caTx8fBSyTbePL4OHiJZBtP\nHh8HL5Fs48nj4+Alkm08eXwcvESyjSePj4OXSLbx5PFx8BLJNp48Pg5eItnGk8fHwUsk23jy\n+Dh4iWQbTx4fBy+RbOPJ4+PgJZJtPHl8HLxEso0nj4+Dl0i28eTxcfASyTaePD4OXiLZxpPH\nx8FLJNt48vg4eIlkG08eHwcvkWzjyePj4CWSbTx5fBy8RLKNJ4+Pg5dItvHk8XHwEsk2njw+\nDl4i2caTx8fBSyTbePL4OHiJZBtPHh8HL5Fs48nj4+Alkm08eXwcvESyjSePj4OXSLbx5PFx\n8BLJNp48/lb4L+uTSiTTePL4G+G/SCTh+/J3gf/yZd0kiWQbTx5/C/yXGI8kknE8efwN8HEe\nSSTjePL47fGRHkkk43jy+M3xkRpJJOt48viN8bGbo5NEso4nj98Wn+CRRDKOJ4/fFJ/ikUQy\njieP3xKfopFEso4nj98Qf/VI59oJD8Knxd+2RxJJeBA+K/6+XyeRhAfhk+Ifx0cSSXgQPiX+\nNVwnkYQH4TPi38PeEkl4ED4hfjTsLZGEB+Hz4cdvH0kk4UH4dPjJ27ASSXgQPht+ejqDRBIe\nhE+Gd04LkkjCg/Cp8LOzVCWS8CB8Jvz8bG+JJDwInwjvOdtbIgkPwufB+z41IZGEB+HT4L2f\nPpJIwoPwWfD+T/FJJOFB+CT4wKdhJZLwIHwOfOhT5RJJeBA+Az58kROJJDwInwC/cLEgiSQ8\nCB8fv3SxIIkkPAgfHr940S2JJDwIHx2/fPE6iSQ8CB8cv3IRyGYiHe9fLxXzXSLh48njF+LX\nLqbaSqSHJ48va98lEgGePH4ZfvWixI1EOp4l0t7w5PGL8OsX9267ayeRdoQnj1+Cj7hIPopI\n/7tWBFal2rxuF8lvUdoiWcSTx8/GR96qPCGJRDKNJ4+fi4+9VXlCEolkGk8ePxMffavyhCQS\nyTSePH4ePv5W5QlJJJJpPHn8HHzKrcoTkqSLpDMbdoQnj5+BT7pVeUKSJJFyqmrcnBK+Ix8O\nn3ar8oQkEsk0njx+Mj7xVuUJk0ok03jy+Kn4NI8kkvAofCx8okcSSXgUPhQ+1SOJJDwKHwmf\n7JFEEh6FD4RP90giCY/Cx8FneCSRhEfhw+BzPJJIwqPwQfBJb8Om408SyTqePH4kPtMjiSQ8\nCh8Cn+uRRBIehY+Az9VIIgkPw++Pz94cxeFfk0ok03jy+Ov4Eo8kkvAo/N74Io8kkvAo/M74\nIo0kkvAw/L74Qo8kkvAo/K74Uo8kkvAo/J74Yo8kkvAo/I74co8kkvAo/H74Ch5JJOFR+N3w\nNTySSMKj8Hvhq3gkkYRH4ffBl70Nu4r3TiqRTOPJ4/vxtTySSMKj8Hvgq3kkkYRH4XfAV9NI\nIgkPw98cX29z5MWHJ5VIpvHk8Wf4qh5JJOFR+Bvj63okkYRH4W+Lr6uRRBIehr8pvrZHEkl4\nFP6W+OoeSSThUfgb4ut7JJGER+Fvh2/gkUQSHoW/Fb7ycJ2Lj5hUIpnGk8d/4tt4JJGER+Fv\ng2+jkUQSHoa/Cb6VRxJJeBT+FvhmHkkk4VH4G+DbeSSRhEfht8c39EgiCY/Cb41vNFz3xCdM\nKpFM48njt/VIIgmPwm+Lb6uRRBIeht8U39ojiSQ8Cr8lvrlHEkl4FH5D/NUjnPQSyTaeN/5t\ne4STXiLZxtPGv+/X4aSXSLbxpPGfw9446SWSbTxn/NfbRzjpJZJtPGX893AdTnqJZBvPGH80\n7I2TXiLZxhPGH799hJNeItnG88WfvA2Lk14i2cbTxZ+ezoCTXiLZxpPFd8/2xkkvkWzjueLP\nPjWBk14i2cZTxZ+fpYqTXiLZxjPF95ztjZNeItnGE8X3fWoCJ71Eso3nie/99BFOeolkG08T\n3/8pPpz0Esk2niV+4NOwOOklkm08R/zgxYJw0ksk23iK+OGLbuGkl0i28QzxFy5ygpNeItnG\nE8RfulgQTnqJZBuPH3/xols46SWSbTx8/OWL1+Gkl0i28ejxVy4CiZNeItnGY8dfvUY+TnqJ\nZBsPHX/9XhM46SWSbTxy/Ihre+Okl0i28cDxY66Rj5NeItnG48aPutcETnqJZBsPGz/uni04\n6SWSbTxq/Mh7H+Gkl0i28aDxY+8hhpNeItnGQ8aPv8UyTnqJZBuPGD/hVuU46SWSbTxg/JRb\nw+Kkby+SSpVSt1vDMpe2SBbxcPHTblWOk14i2cajxU/zCCi9RLKNB4uf6BFQeolkGw8VP2G4\nLgefXhJJeBB+Cj7dI6D0Esk2Hih+ukZI6SWSbTxO/ByPcNJLJON4mPhZHsGkl0jW8Sjx8zxC\nSX+SSNbxIPEzPQJJf5tUIpnGY8TP9Qgj/X1SiWQajxA/Y9g7BV9QEkl4EH4EvsAjgPSvSSWS\naXz/+AUaAaR/TyqRTON7xy/ZHEXgC0siCQ/CX8EXegTUfIlkG983fqlHQM2XSLbxXeOXaoTU\nfIlkG98xfvHmaBlfoySS8CD8ML6GR0DNl0i28d3iV/EIqPkSyTa+V/wqGiE1XyLZxneKX8kj\noOZLJNv4PvFreQTUfIlkG98lfjWPgJovkWzje8Sv5xFQ8yWSbXyH+BU9Amq+RLKN3z5+TY+A\nmi+RbOM3j1/VI6DmSyTb+I3j13kbNoivXRJJeBD+FF/bI6DmSyTb+E3jV/cIqPkSyTZ+y/jV\nNUJqvkSyjd8ufv3N0Qmp+RLJNn6z+E08Amq+RLKN3yp+G4+Ami+RbOM3it9GI6TmSyTb+G3i\nt/IIqPkSyTZ+k/jNPAJqvkSyjd8ifjuPgJovkWzjN4jf0COg5ksk2/j28Vt6BNR8iWQb3zx+\nU4+Ami+RbONb89t6BNR8iWQb35bf6G3Yd+E0XyLZxjflN/cIqPkSyTa+Jb+9R0DNl0i28Q35\nN424uyORhO/Nf2yOuLsjkYTvzH/u1nF3RyIJ35f/Ojzi7o5EEr4r/z3KwN0diSR8T/5otI67\nOxJJ+I788ag3d3ckUmP84XBoiR8Xn0iTd4+4uyOR2uIPh5FJdOkb86fvwnJ3RyI1xR8OY5PY\n0jfmO2czcHdHIjXFS6RwuWcFcXdHIjXFS6Rgzc6u4+6ORGqL1zGSvzxnqXJ3RyI1xmvUzle+\ns725uyORhN+e7/3UBHd3JJLwm/P9Hz7i7o5EEn5jfuhDfNzdkUjCb8sPfhiWuzsSSfhN+eEP\nlXN3RyIJvyV/4doM3N2RSMJvyF+6xgl3dySS8NvxF68VxN0diST8Zvzla25xd0ciCb8Vf+Xa\nddzdkUjCh2t8flMxf+0akHTdycVLJGv4yRm3pfzVa6mydScbL5GM4aefASnkr1+TmKw7+XiJ\nZAxfUaSYa3uTdScfL5GM4euJFHWNfLLu5OMlkjV8rWOkuHtNsHUnGy+RzOHrjNpF3rKFrju5\neIlkG5/Jj771EXd3JJLwTfnxtxDj7o5EEr4lP+FWfNzdkUjCN+Sn3NGSuzsSSfh2/KQ7w3J3\nRyIJ34yfdodl7u5IJOFb8RPvVM7dHYkkfCN+okfk3Wkr0vFWzx8e38+e7xIJH5/GTxiuy8Gn\nFw4+Q6S7TY9/z2/H+XeJRIBP4qd7RN6d9iKNZZFIxPgUfrpG7N3ZRqTj6GeJRIpP4Od4RN6d\n5iK9DoVeQnlF+t+1UvxU4dbVo94ZKCpdpMcXbZGI8dH8rO0Re3e22SI9f5JIxPhYfqZH5N1p\nLdJx8qNEIsZH8nM9Iu/ONiJp124H+Dh+tkfk3dlOpOXBBokEj4/i53tE3p2Ndu1CZzTozAYe\nfAy/wCPy7uhcO+Gr8TNOZ0jBlxUOXiLZxq/yyzwi745EEr4Sv9Aj8u5IJOHr8As1Yu+ORBK+\nBr90c7SCr1A4eIlkG7/Er+AReXckkvDl/BoekXdHIglfzK+hEXt3JJLwpfw6HpF3RyIJX8iv\n5BF5dySS8GX8Wh6Rd0ciCV/Er+YReXckkvAl/HoekXdHIglfwK/oEXl3JJLw+fyaHpF3RyIJ\nn8uv8jZsGF+5cPAIIk1valq9cJoNiHf5lT0i7w6XSM5ttqsXTrMB8Q6/tkfk3aES6XBobBJO\nswHxU35tjdi7I5Ems9gOzY8f86tvjk7s3ZFIk1lsh+bHj/gtPCLvDpVIOkbqiX/zm3hE3h0u\nkTRq1xH/4jfRiL07ZCLhdMMe/slv5BF5dySS8Gn8Vh6Rd0ciCZ/Eb+YReXckkvAp/HYekXdH\nIgmfwG/oEXl3JJLw8fyWHpF3RyIJH81v6hF5dySS8JHV5m3Yd3F3RyIJH1etPeLujkQSPq6a\ne0TdHXsiLZ9jhNNsNPxVI+L4UPg9iLRy1itOs7Hw980RbXww/A5EWvscBk6zofCP3TrW+Gh4\niVRapPjn4RFpfDi8RCotTvxrlIEzPh5+ByLpGCm93qN1lPEB8XsQSaN2qTUa9WaMj4jfhUjC\np9X43SPC+JB4iWQPP3kXli8+Jl4imcNPz2agiw+Kl0jW8M5ZQWzxUfESyRjePbuOLD4sXiKZ\nws/PUqWKD4yXSJbwnrO9meIj4yWSIbzvUxNE8aHxEskO3vvhI5742HiJZAUf+BAfS3x0vEQy\ngg99GJYkPjxeItnABz9UzhEfHy+RTODD12agiE+Al0gW8AvXOGGIz4CXSAbwS9cKIohPgZdI\n+8cvXnMLPz4HXiLtHr987Tr4+CR4ibRz/No1IMHj0+Al0r7xq9dSxY7Pg5dIu8avX5IYOj4R\nXiLtGR9xaW/k+Ex4ibRjfMwl8oHjU+El0n7xUbeawI3PhZdIu8XH3bIFNj4ZXiLtFB976yPQ\n+HR4ibRPfPQtxDDj8+El0i7x8Xfig4xPiJdIe8Qn3NESMT4jXiLtEJ9yZ1jA+JR4ibQ/fNId\nlvHic+Il0u7waXcqh4tPipdIe8OneQQXnxUvkfaFjx72zuSnlhm8RNoVPtkjrPjE+D2JFLgD\nJk6zm+OTNcKKz4zfkUihezLjNLs1PsMjpPjU+P2IdDgETMJpdmN8jkdA8bnxEqm0YPBZHuHE\nJ8dLpNJCwed5BBOfHb8fkWwfI6UP16Xxc8sMfkciWR61y/YII/4O8HwiBXSphU8uBHy2Rhjx\n94CnEym0A1cJn14A+AKPEOLvAt9epLr1GFLoHQOqrh71zqB6FccWKTg2VwefUd3xJdsjgPg7\nwUuk0uqNL/Ooe/y94NlE0jGSU4Ue9Y6/GzydSBq1G1f+sHccv7jM4PlESq4d48s92nN3NsVL\nJGJ8uUZ77s62eInEi6/h0X67szFeItHiq3i02+5sjZdIrPg6Hu21O5vjJRIpvpJHO+3O9niJ\nRImvMFy3yK9WZvASiRFfz6M9dqcLXiIR4utptMfu9MFLJD58TY/2151OeIlEh6/q0e660wsv\nkdjwdT3aW3e64SUSGb6yRzvrTj+8ROLC1/ZoX93piJdITPiKw95efv0yg5dIRPgGHu2oO33x\nEokH30CjHXWnM35BpP++DsP5849EwsC32ByddtOd7vigSP8+hkudh+E3v0jJn05PwzetJ76R\nR0CrIjc+KNK34fvFovPP4ZNepPTrpSThm5HH+FYeAa2K3PigSBeJXv+4Rcq4gldKbbIsW2mE\ntCpy4yVSaW2wLJttjk5IqyI3PijSY9fu+/BNIi1W+2XZ0iOgVZEbHxTp33G41fEvu0jntiY1\nX5ZNPQJaFbnxQZHO5x8fw/Dx/V+hRxKpEN9UI6RVkRu/IFKlqho3p7hFauwR0KrIjbcgUsb1\nwpPwLau1R0CrIjc+KNLwquO3osOkqnFz6pxxvfA0fLtq7hHQqsiNjxDpUiUmVY2bU8T49h4x\ndwcKHxTp1/B50efv5/CzcAi8atyc4sXfbsXXDn8r3u5g4YMifQz38brho/BN2apxc4oWf9se\nSSQOfFCkpzzFZzdUjZtTrPj7fp1E4sAHRfp87tp9nn9ft0oSaVv8821YicSBD4r0931mwzD8\nJ5E2xr9OZ5BIHPigSOd/rzMbhu8FHkmknHqfFiSROPBhkWpV1bg5RYgfDXtLJA68RMLDT85S\nlUgc+LBI35/vxkqkbfHTs70lEgc+KNLLo84ilZ/dg9PsqHI+NSGROPBBkY7Dn8/h77/PrS9+\n4ohT4XxTnGbHlHtWkETiwAdFumyJfgy/zv/aX/xkYoojTo1PQOA0O6JmZ9dJJA78kki/ru8f\nNd+1m5jiimNNpPlZqhKJAx8U6evw8+/wcf7dWqSpKsZF8pztLZE48EGRrgZ9XscaGl/8ZFkk\nW8dIvk9NSCQOfFCk86+P66WEys5qSBZpLo6hUTvvp48kEgc+KNKvUoEiRfIML+TM8tIsVua1\nwvs/xSeROPBBkS77dHVcWo3Q9MrcJ6RmL1bg07ASiQMfFOnr9fjo68/iq3HpzIaoCl68TiJx\n4IMinc9/rvd1GT5/SqT2+PBFICUSB35BpGt93/gUoRa7eTjNDtbCxVQlEgd+SaTf34/D8PFj\nQ5GaXIAOp9mhWrpYkETiwAdF+nW1aONjpDaXRMVptr+Wr+0tkTjwQZEu+3Tl1/2WSOu1co18\nicSBD4r0+75FKr0XhURaqbV7TUgkDnxQpKdLH6WnNqTEtXeMtHotVYnEgV8S6XoBFI3aNcWv\nX5NYInHgF0T68+N60mrx+Q1V4+ZUOj5J55L0Edf2lkgc+KBI3441LKIUKW0HsyB9zDXyJRIH\nPijSdufaAXXjVolDHvnpo+41IZE48EGRNjv7G+6k1a1Eirtni0TiwAdFqlZrCZreTu82i4nT\nbyRS5L2PJBIHvrtIbW/wepvF1BdscowUew8xicSBl0ieaj9qt/Y2bCE+oXBWRW68ROqBj/dI\nIpHgu4vU8hjpDsZp9rMSPJJIJPj+IrUbtXsoitPsRyXdYlkiceABRGpl0nOnEafZt0rZHGXg\nUwusO7R4AJFa7dthipTokUQiwfcXqdloA6RIqR5JJBL8jkVCPEZK1UgiseD3LBLeqF26RxKJ\nBN9fpObnCOE0O8MjiUSCBxCp9VmrMM3O8UgikeARRMLpRlN8lkcSiQQvkbbC53kkkUjwEmkj\nfKZHEokEL5G2wed6JJFI8BJpE3y2RxKJBC+RtsDneySRSPASaQN8gUcSiQQvkdrjSzySSCR4\nidQcX+SRRCLB54h0vNbj+3nhe5pI3ruZ1zjpoW+zk0/3TsMXF86qyI3PEmn07Rj+nibS/Obm\nl19zTsObvaJrs0s9kkgkeBSRnHPAD+NKmvf5K3o2u1QjicSCzxDpOP6OJpLnJf2aXbw5WsZX\nKZxVkRufI9LzEGlVpP9dKw76EMD5dfpgBqdr3TzqHUK1daVtkRYE6nqMBLRFqrE90haJBZ8h\n0tOmuiLVGrWDOUaqopFEYsHjiFSrMEbt6myOgvh6hbMqcuMzRGqza9eueuCreSSRSPCZIsUN\nNsSJVO+t18AsNuIu4Ot5JJFI8BkirZ7RkHZmQ/5br7Gz2Aa7gK+nkURiweeIlFYrAfLfeo2e\nxSbUBXxNjyQSCV4iVcBPk1f1SCKR4K2LVP5Hz85uaV2PJBIJvrtIfY+RKvzZ8/Q/gcoeSSQS\nfH+RIkft8lf4hQQ1NoRTkWp7JJFI8P1FiluVC1b4LUWq7pFEIsF3FyluXS5Z41uLdGrpkUQi\nwROKlLrmtz5GeiWq+DbsBN+0cFZFbjyfSMnrfvNRu0c18UgikeC7i5R6jJS+N7ZRs9t4JJFI\n8P1FitwojHbsIEVqo5FEYsEDiJTWDUyRGm2OThKJBU8nUuVjpPK64tt5JJFI8OgieaSZPbQi\nVvtmN/RIIpHgwUWKHtEb/ZKAr1HnZodHD3zbwlkVufHYIsUcEE2m8U3eutlNPZJIJPh9ieSd\nvnGz23okkUjwEqmwGnskkUjw2CKlHiNtL1JrjyQSCR5cpIh3a7seIzX3SCKR4NFFWq/pVmjb\nUbvbNYnb4a8lkTjwPCKFtk2rx1Htmn3bHuEsS0S+GTyNSGFb1o6jWjX78TYszrJE5JvBs4i0\ntN3pc2bD83QGnGWJyDeDBxAp6tS59HNV37OY86LVep0WVBnvzqRE4sD3FylOkPYipdHfw3V1\nl+VsLiUSB767SIdDlCPNRUrCj89Srbos57MpkTjwOCIlnFGXVvFjGbH8ydneEkn4E5RIi2tx\n5Irum6S6SNNPTUgk4U8AIiVtkqJYs1mcTRNOsfYHruWczaBjJOFPCCIlHCVFktxZnE8Teu3a\nH7iWe1aQRu2EPyGIFL1BWn0kRqTw38r0CGhZIvLN4LuLlLtj53tNkUhzlu/h+VmqOMsSkW8G\njyJS6kT+l70fGz2XJVJgKs/Z3jjLEpFvBg8j0uLaHSvS9Op33m4keDSfzvepCZxlicg3g+8u\nUtT9kaJF8j0bNWq3+gev5f30Ec6yROSbwfcX6fTyaM2klUfcJwMiRZU3j/9TfDjLEpFvBg8g\n0jXu6qGL93hoYeJCkTyehi5eh7MsEflm8AgiHQ6HgEhBWRbFc2h5zY71CGhZIvLN4AFEmnvk\nGzLwv2TpyYxuhCt8MVWcZYnIN4PvL5LXo9svQVnCm7DRs043osYYgrVwkROcZdmPH24ud3c4\nRXIeyBVp9tT59WhMQ3yTLV7bu8eyLPtfYZ2fWAvNxVnTG+N3K9K0GwsbMM+rnceWr5HfYVlG\n/6+QyU+spebirOmN8TgiuQ+cMo+RikTyTbdyr4ntl2X0/wqZ/NSSSCcgkUa7K+/FElxdlpad\n+0yhSGvXgJRIEumEIJLn1IbYnbC4Z1KOkebrxOq1VCWSjpFOMCLVXTmcWXz+lZiJkz3SMdJJ\no3YnBJEOYZFSVpdayzLVI43aCX+btLdIh0nNn4qcjzZ7FzHXyMdZloh8M3gskXzPRM1Gm+Pd\nqHtN4CxLRL4ZPJRI3meiZqOJSHH3bMFZloh8M/juIoXfXQ08vATxz2IMwFOR9z7CWZaIfDP4\n/iKF311N3iQFZjEK4NbK27Cl+OiSSBx4AJEucf0aJA2Jh4f8spod7RHQskTkm8FDiBTypeT9\nkvdrc5od7xHQskTkm8EjiBT2ZdmjpWdHW7OMZqfcGhZnWSLyzeABRJruwa1vg55TLG6vSkRK\n2Bzl4NNKInHg+4s0HWpY35t7TrF8BFUgUppHQMsSkW8G312k6Zjd+vjCa4qVSbOPkRI9AlqW\niHwzeBiRJr+tT78+aeaoXaJGSMsSkW8GDyNS8hYpekwvqdnJHgEtS0S+GXx3kbKPkaLPgU5p\ndrpHQMsSkW8G31+k7FG76FmMnzTDI6Blicg3g4cTqX7FdyPHI6Blicg3g+8v0vQgqUHFdiN1\nuC4Rn1kSiQPfXaTDukjuc4nWRXYj0yOgZYnIN4MnEMl9MnUDFteNTI2QliUi3wweR6SVCeJf\nMJvFmImyPUpdlg2HSrIKZ1XkxvOI5LzVFL9CxnQj36PEZZl8NCiREPDrC41GpNc0B/eB2G4s\nTF7gUdqyTB9XkUgA+IiFhiNSM5POd4UWJi/xSCLtHx+z1NBEeued/OROk2LSeWXy3OG6Jz5l\nYonEiCcUybFlNpH/VbF/4fA8R2/8dKFHOkbaP55LpGniafawR+kiOa8q1EijdhbwEWsajkhT\nfbzmzF8Ts1J6PHq/rtgj8lWFPL5G7d7lX8tdT5yfEzx6XVvFBy73iHxVIY+Pg0cSybtpWnhF\ndDc8Ft5/qeAR0LJE5JvBdxcp3aPxszE6OQkqewS0LBH5ZvCMIp3GHq2b5Cao6xHQskTkm8FT\ninQae7R+GOh/uHTYewVfqyQSBx5IJHdLs/65ihKRankEtCwR+WbwMCKdkjxK2iL5JqmlEdKy\nROSbwQOJ5PXDI8F0rGH9GMk3UT2PgJYlIt8Mvr1IK/UU6X1G3POn17P+Vzx/jvwDk8euHlXK\nr1JNq+8WyTPqcArsu8Xt0C1NXnF7hPSfIiLfDB5FpMNMpNGTzg5eqUhVPUo+eS6xJBIHHkek\n8fjdYTos55qQYdL712rDdV569ZJIHHgYkaYmOc85a2v0zt1toumo3cijGgYk7meml0TiwOOI\n5Lh0mOzYRYk0e+g+1Xn85GhzVMUAiST8fVIkkXxDDu4WyWvWm+RBn0dPuh45esb3zfkb6S+M\nLYnEgccSaalGk3s3JvNVeiTS48fx4ZE7faYQjT2SSCR4CpHeG4vZA0vnC81EmgwzONPPXx5Z\nbT2SSCR4CpEOr8/lBWwZcWbs00skZ7huOnm2SDjLEpFvBs8i0mM1n63uc+Ec+KMbc4+cbYlE\nEr4EDy2SxxN3bV/0aNwN5+2j+bSZHgEtS0S+GTyqSM5HjiZqhV8e7sbcI59J8X1z8O1KInHg\nu4sUMCnw5MrLn4/MuuE7PKo0SICzLBH5ZvDdRQptkbxPrrz6/ZDTDffwSCJtxTeDhxNpIsLY\nkaBHsw8wOY6c52epoom0kEUiceDhRHLfI1pa4X2SeUTynKVaz6May3IpjUTiwAOK5Hty8bUr\nj3nP9q7mUYVluTiPEokDjyeSd8/Ov5Z5n3MeqvypiXlJJOFPqCK999giTPI8+P7F51G9rdG1\nJJLwJ2iRDh6RDq93lya/T3DjX28XZ/D9xfgOrZWOkYQ/oYm04NV74qUVbyzg6XFW0NkzSUWT\nNGon/AlNpIWt03zaddh9t45ApH548vg4eBKRfNOGzpebeHSQSF35ZvAcIoWnXYDdPJqLBHiM\n1BFPHh8HzyHSxKSlXbyZRx6R4EbteuLJ4+PguUQaj9oFTTo9Nbp7BNRsQDx5fBw8lUjzrZNX\ns9PEI6BmA+LJ4+Pg0UUaH9E43vg9u/00eRcWp9mAePL4OHh0kUam+DZA40feP0/PZihs9mHl\niApnWSLyzeBhRfJXwy8AAA7hSURBVFp4p3b+6unPzllBZc12/+KscJYlIt8MHlUkz3tGHo9O\nnj2/2dneRc32/c1p4SxLRL4ZPKpIkR7NP4o0/9TE2TdtbEkk4eMmRRfpvvYHPXJZnrO9z6On\nk02SSMLHTQoukmeqJZrvUxPPBBEvD+VbmgBnWSLyzeAJRHI2SUvrtfdDfGUiadRO+KhJwUV6\nvZEUYdLs2t7TbmSKtFY4yxKRbwYPL9Jr9Z+Z5ErhuUb+tBtNPAJaloh8M3gekdzjJFeL19ne\nY+y0Gy08AlqWiHwzeCKRxlN7dtRGZ9edAiI1KW48eXwcPIlIjymnl7KbiPTwSCJh8c3gSURy\ntkQnV6QvjkfOMdKhxT7do3CWJSLfDB5OpBWzxp+VeLtyf/toqsvrt3OjUYbTE9+0JBIHHk6k\n1U2S+8pT4O2j9yzOX/p4eXybFgpnWSLyzeB3INLatVT9IlXbSuEsS0S+GTy/SKvXJPaK5N9K\n5RTOskTkm8ETi3Sae+RV43mMNHlSIm3DN4PHF2l6pRP3tY5GPjfOo4uluH82vlHBwlmWiHwz\neHyRfBc+eZa7OfLKcR4/O/278X0KF86yROSbwROINLns96S+PK6lOgX5uzF7spJHQMsSkW8G\nzyJS2CP3GkOxItUqnGWJyDeDJxHJp8HzZAb/xbpm3WjkEdCyROSbwROIdPJL8vJoum8X7kYb\nj4CWJSLfDB5bpOfK7xHpeVLQdPygsBs5xY0nj4+DxxbJ+RzS6GWP4bqYrQxOswHx5PFx8Ogi\nBa7ElXKHZZxmA+LJ4+Pg4UVyNk/3SrpTOU6zAfHk8XHwTCK9XpPkEVCzAfHk8XHwRCK9XrHk\n0cH/PpLn4VqFsywR+WbwJCKNTFg823u65Xp1w/dwrcJZloh8M3gSkd4mfJlcLMgPc7vhfbhW\n4SxLRL4ZPJNIzw/Dji4W5Ie53ZBI3fhm8GQiOacFBWBuNyRSN74ZPJVIX1Y80jESHN8MHlek\n+ZOPYYbFzYvnGY3adeSbweOKdHAvzfUerlvbwEyfxWk2IJ48Pg4eWqSZR5OXLCOzupFT3Hjy\n+Dh4CpGcizOsiOQ+jdNsQDx5fBw8hUgH30VOJBID3wyeRKTxxRncfTfXKYkExDeD5xDJPZvB\n8eh9fRTnkdRu5BQ3njw+Dp5CpLizguZGJXcjp7jx5PFx8AwiLXkUvmFSTjdyihtPHh8HTyDS\n8qePJBI03wweXqTVa+RP9uwkEhjfDB5dpFWPgmMMOd3IKW48eXwcPLhIER5NUIXdyCluPHl8\nHDy2SGkXZwjNYjFhx3jy+Dh4ZJG+TN+GzS2cZgPiyePj4IFFWvxQeULhNBsQTx4fB48q0un1\nKb6iTtxmsZiwYzx5fBw8qEiea+RnF06zAfHk8XHwoCId3h5JJGa+GTyoSCOPJBIz3wweU6Sa\nHgE1GxBPHh8HDynSa7ju5VGJUDjNBsSTx8fBI4o0H64r2jThNBsQTx4fB58j0vFSz+9H5/fx\n90yRAh5lm4TTbEA8eXwcfIZIx+eX4/R393umSM/jo/k0Oa04ITUbEE8eHwcPJ5Lv/SOJRMs3\ng88Q6WnP8f1jPZH843U6RmLlm8GXiPQ8RAqK9L9rrYDCHh2uz72nSoqnUvWqJJGWBCrYIk08\nKtqfG/9fUY7YL548Pg6+QKTnD/VE8nikN2Tb4snj4+DzRDqOf6ol0hevRz6R0vTCaTYgnjw+\nDj5LpOP7az2RAh55lEncUOE0GxBPHh8HnyPSaNh7ebAhRaSARiGPEkzCaTYgnjw+Dj5DpOPa\nGQ05ZzY4Hp3CHkkkJr4ZfM4WKa1WAvg9WrJFIhHxzeAxRJp5NDsQGv+iYyQevhk8hEju8dHz\nielU4Ru5VOtGTnHjyePj4AFEmg3X+SZK3ApNZjHvZTbw5PFx8AAiHVY9kki0fDN4BJGcCk+U\n0oPRLOa9zAaePD4OnkOkotO/cZoNiCePj4PHEyk4WfxMObOY+0ILePL4OHg8kSqc7+3MYm3g\nnvDk8XHwEsk2njw+Dl4i2caTx8fB44lUNOfeWaxO3BGePD4OHk6ksotBemexKm1nePL4OHg0\nkQovdOKdxZqwveHJ4+PgwUQa/V7UgcksViPtEE8eHwcPJtLhdKhuEk6zAfHk8XHwaCI1GHbA\naTYgnjw+Dl4i2caTx8fBSyTbePL4OHhgkYpaMJ7FWqA94snj4+DhRLq/jVTRI6BmA+LJ4+Pg\n0UQqmu/ALDZg7gZPHh8HjyVS0VwHZ7EJdSd48vg4eCaRMk3DaTYgnjw+Dp5IpNxtFk6zAfHk\n8XHwUCLFTBg/a69ZTH+JHTx5fBw8kkhRE8bP2msW019iB08eHwePJNKyJRKJkW8GjyXSuknx\nc/aexYzXmMGTx8fBg4kUuNp36JG4Wcx5kRU8eXwcPKBIr6+V3lfCaTYgnjw+Dh5PpPu37CMi\nzyzWgOwVTx4fBw8mUsoxU/Qs1oDsFU8eHwePJZL7u0SSSCR4JJHGv+sYaRs8eXwcPJBI4wdO\nFT9JgdNsQDx5fBw8nki6rt2WePL4OHhAkWoXTrMB8eTxcfASyTaePD4OHkikVibhNBsQTx4f\nB48kUiOTcJoNiCePj4OXSLbx5PFx8BLJNp48Pg4eSaSieV6YxUbcXeDJ4+PggUQqmuWlWWwF\n3gOePD4OHkakohlensV2aH48eXwcvESyjSePj4OXSLbx5PFx8BLJNp48Pg4eRqR2JuE0GxBP\nHh8HL5Fs48nj4+Alkm08eXwcPI5IekO2B548Pg4eSSSdIrQ9njw+Dh5KJN0faXM8eXwcPJZI\nTUzCaTYgnjw+Dl4i2caTx8fBSyTbePL4OHgUkU4SqQuePD4OHkakdibhNBsQTx4fB48mUgOT\ncJoNiCePj4OXSLbx5PFx8BLJNp48Pg4eR6SKl813ZrEBczd48vg4eCCRql/0+zmLLaB7wZPH\nx8EjidSocJoNiCePj4OXSLbx5PFx8Cgi6aPmffDk8XHwMCK1cwmn2YB48vg4eCSR9Hmk7fHk\n8XHwUCLp80ib48nj4+Alkm08eXwcvESyjSePj4OHEqlopsOz2Aa7Dzx5fBw8kkhF87wwi424\nu8CTx8fBo4hUNL8rs9iQTY8nj4+Dl0i28eTxcfASyTaePD4OXiLZxpPHx8FLJNt48vg4eIlk\nG08eHwePIpLO/u6DJ4+Pg5dItvHk8XHw7UVaqadIff66SlW7dIxkEU8eHwePIpJ27frgyePj\n4CWSbTx5fBy8RLKNJ4+Pg0cRqWh+V2axIZseTx4fBy+RbOPJ4+PgJZJtPHl8HLxEso0nj4+D\nl0i28eTxcfASyTaePD4OXiLZxpPHx8FLJNt48vg4eIlkG08eHwcvkWzjyePj4CWSbTx5fBy8\nRLKNJ4+Pg5dItvHk8XHwEsk2njw+Dl4i2caTx8fBSyTbePL4OPjuIp1aewTUbEA8eXwcvESy\njSePj4PvLpJ27briyePj4CWSbTx5fBy8RLKNJ4+Pg5dItvHk8XHw3UXSYENXPHl8HHx/kU6N\nPQJqNiCePD4OHkAkoG7Yw5PHx8EDiKQtUkc8eXwcfH+RdIzUE08eHwffXSSN2nXFk8fHwUsk\n23jy+Dh4iWQbTx4fB99dJB0jdcWTx8fB9xdJo3Y98eTxcfAAIgF1wx6ePD4OXiLZxpPHx8FL\nJNt48vg4eIlkG08eHwcvkWzjyePj4CWSbTx5fBy8RLKNJ4+Pg5dItvHk8XHwEsk2njw+Dl4i\n2caTx8fBSyTbePL4OHiJZBtPHh8HL5Fs48nj4+Alkm08eXwcvESyjSePj4OXSLbx5PFx8BLJ\nNp48Pg5eItnGk8fHwUsk23jy+Dh4iWQbTx4fBy+RbOPJ4+PgJZJtPHl8HLxEso0nj4+Dl0i2\n8eTxcfASyTaePD4OXiLZxpPHx8FLJNt48vg4eIlkG08eHwcvkWzjyePj4CWSbTx5fBy8RLKN\nJ4+Pg28v0nr9r9tfrlHc6cnjI6aXSHnFnZ48PmJ6iZRX3OnJ4yOml0h5xZ2ePD5i+n4iqVQ7\nKomkUlUoiaRSVSiJpFJVKImkUlUoiaRSVaiNRDpeavGJ4AQYFcjGkZ67+STptxHp+PoSeCI4\nAUaFFuXjC3Z67uazpO8m0nHyBEg3AnWcReNKfyZuPkv6jUV6b4YhuxGsR9ql9Ni12nzoIki/\nrUgjXwhFWkgPsZu+VMvNR49PkL6PSMd7kYo0T388Q6d//Jf++IEt/vHIkH5jkY6vDTTjFok1\n/XvfiDc+evoOW6TRQ2jdCJabjys9e3yGVUciRRX3mjj+j9x9AD0+y6rTbbBh5QmwOk6/OI+C\npz+OvtHFZ1l1Nj6zYfYuNNbb08F6/CdImf51dMEbf/x9/YkupXPtVKoKJZFUqgolkVSqCiWR\nVKoKJZFUqgolkVSqCiWRVKoKJZFUqgolkVSqCiWRqGvIXX7/AZwMsKuSSNSVLVL2C1X+Uj+p\nSyKhlPqJXP++DcO3f+fbev91+Hr++zF8/ff47fPv+eHD3+tUf8//ho/bay7fppOPMX+/Dsfv\n1x9kUt1SO5HreF3hr3oMw9fLTz8/Ll++XX+7mDEc70qd/92muvz2ffh1mfTn8MOZfIy5/fhd\nIlUvtRO4flxW+Yse/93UuRhy+fXnVYBh+Px3/hy+30X6Pnyeb7/9uf5w2RT9diYfYy4v/G84\nateueqmfwPVxWzqXfbTrPtn1y2MbNAx/Ljt0123M9beP63O3375eH79rMp7cwTwQXedsf6V+\nAtfwqMd6//pyt+D50/u3PxdZft33/aaTTzASqUWpn8CVKtJ143Q7UJJIm5f6CVwfr6XjmnHb\nmft0d+0um6Pvx2E2uYORSC1K/QSu79dRgp9PYcYifZ7/fQ4/nMGG89WZ24DDdHIHI5FalPoJ\nXPeB7dsIgivSdcD77Ax/n6+bpOHn2Z3cwTwQOkeoakkk5Lq+1fr5++zZtfu8vgU7fUP2Vu9x\nudGXKeb29T+JVLckEmGF98t+389uUG1eEomwwiJ93k5uUG1fEomwQiIN96EGVYeSSIQVEul4\nPXtB1aUkkkpVoSSSSlWhJJJKVaEkkkpVoSSSSlWhJJJKVaEkkkpVoSSSSlWh/g8D9POnU4xk\nqgAAAABJRU5ErkJggg==", "text/plain": [ "plot without title" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Scatterplot of wage vs employment with regression line\n", "\n", "library(ggplot2)\n", "ggplot(data = df, aes(x = employment,y = wage)) + \n", " geom_point(color='cyan') +\n", " geom_smooth(method = \"lm\", se = FALSE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Note: the above relationship does not look particularly linear overall. We might want to add employment squared as a regresssor as well, or use so-called *non-parametric* regression techniques that let the curve take an arbitrary shape.\n", "\n", "### By comparing the linear regression line to the `loess` local smoother, we can assess the linearity if the relationship. Here we zoom in on counties with average wage between 20-30k and employment below 100k, where the bulk of the data is." ] }, { "cell_type": "code", "execution_count": 128, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAS1BMVEUAAAAA/wAA//8zMzMz\nZv891tZNTU1cetZoaGh8fHyMjIyampqnp6eysrK9vb3Hx8fKysrQ0NDW1tbZ2dnh4eHp6enr\n6+vw8PD///8gDbNwAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2dDZujts6Gp7w7\ns92253S6PaX5/7/0nSSAJVmyJdsEG6Tr2p0A5vGX7ggbQ95ubm5u1fZ2dAHc3M5gDpKbWwNz\nkNzcGpiD5ObWwBwkN7cG5iC5uTUwB8nNrYE5SG5uDaw1SHPWFEnstouoF3UY0aOK6iAdp3r1\nop6q/g7ScapXL+qp6u8gHad69aKeqv4O0nGqVy/qqervIB2nevWinqr+DtJxqlcv6qnq7yAd\np3r1op6q/g7ScapXL+qp6u8gHad69aKeqv4O0nGqVy/qqervIB2nevWinqr+DtJxqlcv6qnq\n7yAdp3r1op6q/g7ScapXL+qp6u8gHad69aKeqv4O0nGqVy/qqervIB2nevWinqr+DtJxqlcv\n6qnq7yAdp3r1op6q/g7ScapXL+qp6u8gHad69aKeqv4O0nGqVy/qqervIB2nevWinqr+DtJx\nqlcv6qnq7yAdp3r1op6q/g7ScapXL+qp6u8gHad69aKeqv4O0nGqVy/qqervIB2nevWinqr+\nDtJxqlcv6qnqbwHp48s0fx2kA0UHKuqp6m8A6WP5L/f3zCBN09RetK2N06qnqr+DZFGdppYk\nncqRehEdAKSVpuuCNE1GktKJh3GkpnE42DD1Pwyk/7ubSq7Qvnp2T/lEtpMla1vqLu1egRdW\nY/j2uqlAek4mHB+RUnGho4iUSz7AN/IULD5SKx6XtMWl85ki0s4gJR20ozHS+CBNIki7uLz5\n0lml2sIcpEZWNGt3YpD2cfmLgNTNrJ0JpDYD5bKiZtxiKJDYA3XqDtLBIFnGSC0uQWJVraXz\nvmkSma39GIkp4U4uf5ExUj8rGxLNjUWbdHik2sjADEZz1Va2lI2J87vEjgZNMQBINmtROr1t\nPcCC1PHXXCvUsWprexVIDcxBultp54RuFUDatc8L5R2kPcxBmlWXOSjBugH6lRsj1Xd6sjal\n8laQdMle40jdhnkHaVY5FUqwbcggNbq4S9WmWN84RlImfJEjkbKUNICDdBxIKEXYSIDU5jLE\nApI6L9usnbYahzhSUQs7SN2BJI6RZnyw2Awg6XOzdXnPIJV9VzlIe4GUd0IJJGnWbsYHi00/\nRjL4lIO0hzlId8t2R/bLX9uOtn7Xz9rtBlJnYyRkDlKPIGVJEja0VVZmVKI67wlSV7N2xHyM\n1BdI9m+2KLWuHa0ZGXpnrzGS1o5xpJKr5zPVf3SQ4uSHg2SetWtsBzlSL6JXBQl7ndW/mfTH\ng3Ss6EBF7br+34yiB4NE3bkkIJWAtN8YyWBdO9Kooo1Uv40FUgyCeYRUCFLTWbtC69mRhhVt\novrtbjbR3kAyWukYyWoXc6SBRRuofos4Oj9IpbN2VruUIw0tWq/KcNQ9SOxQpY6sTntnJ9Hy\nG2k2G0e0VpXDaACQpKeYlbUWRNtbXvWQG5JcW43j8z2CxGI0AkiR1V7tHdQ7hyyRYdtqHJ/v\nDyQ+HKlEHaQmqkWldpC6UA0tKGHkIDUzB2kU0ZIVjM8mFMORSrQ7kMYcIx0Dko+RSlVDs209\nl8BoTJB0w3YxERWtmwOUVLnyHPEYQctZu1QFTgYS6K0FpFQ4Uol2CJLGZL+9aROarNNZu5aq\nyZY6F0jw+gFwVCU6JkiJK6mbNqHJzuVInKVb6lz1R3XVcOQgSQmtcJ3LkTi7LEhzHiMHSUpo\nDlPnciTOrgQS7n8FR6cFSTtGSnHkq7+xXWiMhK5INBydFyTlrJ2DZLALzdoFU2F0ZpB0oi1A\neiQ8rSOdTtSmquTo8iBJ1ys2jr6SntWRzidqUc3cPDKJDgqSelgsJJTPJ0eW4HVORzqjqEFV\njdGgICmCRSKiVPYOVXaQakVfPBWoVtWHI5VofyApLrtSY5y63omUHaRK0TbfeYapIaWqBaMh\nQdJMBLwQJDxGMt59ytgVQGrTVZZJVpWqKRypRB2knDKYtbNOmmfMQaoUYYU1qkaMTgvS68ZI\n68ZNXTSDOUh1IryyQtXM0Ygg6b729xvBxiOkTfVAkNqPEWx27BiJb3ihO7KqdoyGBKnWW1v6\nUShLGqTC8u4wRpj6mxdp8Z0nByQ7SCUcXRmkFqGDgNTg9i6y6jFCs5Jk7Wg6hRFSAUhFHF0Y\npCYuRUFi+nPazC7fHKTW157BjgaJtYIxUhlGQ4LU5p0NjVwKj5GkBA7SQaJsXVOqpRwNABJ/\n7VtSVyDazqVWDXkRxitA0n65FJVEdUKfIBlVrTePVKJbimNBsnS8LmVbkLAqV6RyjnaYtSvj\nSHHKGUAqx6h/kFJOSPcrnaTlGClSpVbF0R6ONFln7ZTFHx+kinAki8IUvYJED2gdtuWsXaxK\nrYrXHm7IXgWkKoxGBik6YgWprYmqNbxWr2ZqoHoNkOrCkSBKUhwKkvyNHl0y9RmRjhMVW8Oo\neoUxUi1GZwCJkKSu8mvGSMeJyt8rVlV1q7b9appfB1I9R92DlBsj0Ys7dZUrpgASqo2tD5BU\nZpxgVYvuYFS1AUZDg1TKwp4gNfakoUBq3aLzi0BqwtGwID32lXXcjiC19qQuxkg6GxakNhx1\nD5LgDsveon6rGCMlzrjN7eHsYtZOaYOC1AijAUCyLOvVGZ61CzJ5wVSmHYL0UtUxx0jNOBoA\npGDI5+tBCkrkk2TJXLdv5OuC9MJZu6qcNtXqm0ecqJyiF5CAkzYDKQgpJNMg8RzV+dZgIL1M\ntO77alVtiNFIICE3rWnJfUASOWry3dnUhgepMvI/VVuGo3lYkGq+6YtBSmLBnl97recgsdYC\npLYYjQtSsRGJoKlRT6RwkF4nWg9S43A0jwRSm0mhSKRJkJvXMRKXm4PUXLS2VZtjNBRILSaF\nGoU1xm7yPH2NaNYK9Gt9ns3ytXTWteoOHA0FUgPbE6TdnopPWUl1TK0qhNk6Ua3tIroHRg5S\nM9vzhqRxiiNl98SWVo3lhSyHAWkfjq4G0g7OvtiOS2QSstY8H6kNrcroHwtSdQPvxNHlQGp/\n932x/UBK6RrzLAlgXYFU28IPjA4KnicDabepsIl4WAumsiAZPasFSEeOkWq/q57hyEEChpvT\n1Lj7gMRwpCtiqux5kGxVtzsil55VGACk9eaRgxQMt6etdXdpR5YjTRGTZc+OkYqKWTlrx1v/\nIG2jIwdpM9ygxuZVt2PNRVOGK/E8tqjNOLLP2qntlWOkkvYAaxkcpM1eApJJdFeQGtvAIM0L\nR2aS4GTdtUFCbWcBKTqibEcjnjRxkittLnuC1DDMBdHGxooWXN/hpXWXBom0HeOmgkR8bB+Q\nbjRtiivNgYeoNneThSmM1qKtrRFI5N7RlUGKGi9yU0GBafWdQMrLGQ8oRMtsWc3UlqSOQaL3\nYB2kgn4vB6mHuUCNqL1ZBgfJ2DHxUgYH6bUgdXB3Srf629ouo4Nk6hhmSdCVQSq/pC8eIxnt\noN4pAWLwMZLJuKV1lwapfJIpOpF/cqjWRgNp3Fk7g/ErVK8NUjtjn2WtV20tmBTdKlAOUmMr\nFt1xskVY6e0gbVYCAvW9fR4damvJG9LRx1rVKisVTZa/qqTiexkcpNVKMAjnTLuQ9OLBNiKp\njWqdFYqme6KmpPKDR9cDSWjhEgzAOVNMUgOoDgOpmarShIw7Ayn1mqDTgZSzRxNL+9kjGa3l\nnABSJqOXm7IYJQ3QyhrnvFNVng/w9Wovjkjil1VlRIre0t3kQq/F11xUDMUYqcBqiiq2VU9j\npMxb604XkTIZy/5d4kfhnOg2Si8gxeXIz9qVWFcg7TBrl3stg4MED5lrup1Db6NMA4JUZX2B\nlLIi0ezrTa4GUvOb72uFmFzqM7oISGKn9AKS4i1BlwPJFnjUiXmQ9BmpVO3GAZ0SlYaP+Zx6\nmrVLm11U87at64FkqbKeBiQ6cf5bYrW9w5YjIcoWWleTTnx+D1HVW+scpISlcUBHbvhAG45q\ne4cvhizKplfWpQ+f30FU+fJHB0m2NA/40I0caDMSc5AOF9W+RNVBki0JBDkWg6TNJWUO0gtF\nUSWXDf1PHjlIsuUD0nYUJusHJH5409cYKan98iVSZMPwTu/LgkRcn02ZGSGFwzhhM46azNqZ\nRNmCq2pTsyyutWjaeFH6vfhlpl/guypIsAOT12+iDGh5GoNacTTQ9VKhaiZ6HwqS7ScmLgpS\nFE5Srk8vnqMwVHIxt+fXvE20DfbjgMTnR0Gy/iCsg5TFIL54DifGO5UGo1m6qM2Nu3fcXFVp\nB4AkZYi72fyLRw4S+Mz6dhzz454o4+h+SvrMF/ROSTTNq6rt5fWX6wv32n857KIgsWMk3rc1\nIJkvj1aVA76R+wKpNiKbi66pb8kP8F0VJGbWbkIGD+ZBsrajgzRrINhhCbKivkU/ZHlZkIDB\nSMSQgjalfrC2I85TXVSrKYrahqOCoioyzomWfAvsw9E1QcItSaJD3DkyVpYqE5vAzHmiqHWm\nCp5NOCpqgCwEe4CUqW/p7ypfESQu4gCAir6ii9sxmdeLlwgV5VCqehhISdHi3ye/IEggAsFt\nRBI9wV7lg66XaBn2BQnIHwJS0ZdeQtR680inWm5DgBT6QLymA+mtVS4Ka1lVo+0NEtQ/ZIxU\n9H0li5ZjdHmQGJKE5LyQVOWia47YSnoHZssWoheQWszalZgkWhGOEqpVNhpIcKeUXNBZvaEb\nkHC+3OduQMrbK0WrMLokSOBdjvHtIjEiRUewSC8gJeMny3yNgbwGB6kuHEmqtdY5SM/vZupw\nuZBED2Ec+xgjyXEVENbJrJ3GXiZai9FFQXqaQEZ8Ak9SGiR+AGCEy8qiFDzBsbl9l++j+rBX\nidZzdGWQIktEkUxIUrajMUyh5IozUxztB5KlAaz2GtEGGDlIyEQv5KPVtC6HUFXZPHBCyTVn\nYo5Y7mdlUfW2ZrnLeomXeGcTjhwkbGmShAO6KleBpDo1Q966/RKQjMFXsFd4ZxuOHCRkiW90\nuMn4iGaMtDdI0HsTJ7wCJGNVJdvfOxthdG2Q4p4G/Z9wBPnLPpUmrckJLcmnKfu4RXz6y0Di\nx0ijgNSMoyuDxHR16P+EJ8SHWLIEkpJFIyetHIGxmNaKQCpz/YkZJI4BUvXNI1a1oQ0BEtfX\nZSAxacscKTrrFvZZ1eT8hSaq9P0Bx0gNMWpV1M9Po2inIIX+7w4ks4kn8U00TcU5saotONoV\npJbhaG5S1M+7GUV7BSn0f8Kp6CExtuXLwgi3Akk0tok2jlqB1MR2FG2LUX1RPxczinYLUrDE\nMQ01epcEKanQTRCvsheA1KrEu4HUOBzNtUX9/BwXpJYeWiOEikGEbrXinKVBaqHarGUVvWnP\n6dY+HM11IH1+Dg1S0+/64nZM+u8Lr5dqLyGBarur0Xz9C3K67cFRRVd9fg4OktzfBV4wPEi1\nge8YkEqy2gOj4q76pGYU7RqkZN8Ih8YHqZ1q3yDtw1GuqHwhI4xGBwlVM9k50qFy70x5wpgg\nvXCMZAdpJ44yReVKyVA0JkjkphGutNA54rEK70w4wqAgvXDWzsjRXhhlihq7DU/RoCDNkCMy\nB/1CkBI2KkgvFLVzdED9iduIFI0KEqolCUmqFtl26vzI/EXtIDW1Zzw6GKQURUODNEUgZe7F\nssuKNL1jHzqM4J17qjYVXS/rjqj/2vUZik4BUv7ENTl7vupMK0n9eyer2mqI1LSo2/DokC+S\nR5PkMToBSPnzBLOCZMlpTJDq2lMQrbUwzXBQRNZQNDRItf1uBsmQl1SbqgLf6iVE1UW5FUnN\nfB7O1h0BkpKisUGq9SnbGMnkY/IihNrFPO2iBlSdN+lYvCS/Vj6PZr1fDpKeorFBqvWpyTJr\n1wKkyq/8W72EoPo0XrsovzY+T5Z6vxYkE0VDg9RkjKTsnUuAJD1VUpBhE5+nN2Hj/m+RC1tU\nK0VXBmk5XdnlTcZIYol12ruDxJXjMJDiJ4+waKt2YIpagJGDpO1yW07ZMRKSU9Zj7zESa0eB\nxKwJoouZ2rREVNQijEYGqdKnjCDZLDdrh4qu9Ym9Z+14O2SMxD4I+xKQCjEaGiSNT4lJ1o5g\n863uIdPXvA2k1pZVLWmMyqLyS1RfAFIxRqOClGhC3UXT1g9cvvVddCaQXi8qLPXee4xUQdGo\nIAljDXxoTrloCqTaL7v8BSPJwTJGam7dgSQ+MbHrrF0dRYOCFPwwckHtd30pSLqryWwikkbn\nE9cASX7yaL/6V1M0Nkjcm7XVF03rARtICkh08azk6/QSICWe4Nur/i0wugRIvNpywDRG0kDC\nlil5htYuAFLyQdhdSvr5sw1HQ4K0OSvj2FqOklVOXg6CY1y6qEyaAKWyPUFqBTsSNVv6gfL2\n9b8DcGmQVufkIoTqwi6k1PfOxlFmkoDjqM38UmNvX1Tv/zWDHYpaLfcW1cYgLQBUg/T+/n7/\nH7XfMCAtuEQhAlns9vFh9aJViSOeJMp8Ey9t7e1P22XhUYnPZ99v0hKkwEEtSO93+/oftV/3\nIHHXcSgEMUfhVnxU0zsTMkGfFlWRymbNvf1pnYCkeKl3M5AQCJUgva+GGrB3kKLuRh7AH+VS\ngh2KfHmMBPebornApgHpdSBVZWX2ec3rthqBREi4Ikhrf0/k8gp8FEf5jUCao0HYzO2iq1ny\nmaiL0UIKmjBGqsvL6PO635hoAVJMwoVBooCET6gyzFYcUDS9g/mFMoGjsNBCq2q3ao7Y0/lZ\nu0pqbfVXvv2xulFZEq44RoIgZead42Tsl66udwg7VAjs2hWk2lk7vtn4or4SJO1bVKsaVQTh\nkrN2SZCYyzyUjPcL9awdzOQgkBTFS5jQbkeDpH8ZcWmjpjm45n2kib8LuxjDUdYVrE8RxLIo\np+XTS0FSurwJpJeNkQwv9S5r1BwA1wTpYfk+FjmK9+Ty5eMPkxnK4ZUgaYOHDaQXzdpZXo5f\n0KgKAC4MUr6PExzhST/jk0Ns7lxWxj5nKmSMHeqQFO09cq2d7TcmrCXVAXBlkJCxPsSTxF3x\n2UHi02iKahDgcq0FKTVr19h0osbfajGVVA2Ag/Q0Fhdx1i4mqQFIBReMuTzYbCvHSIIdBpL5\nB8r1JbUA4CA9LPa4ZZsCQzjSgmT20sk+RqoF6cWLEBqJ2n86TFlSIwAO0sMij9t24COUIzVI\nRi9dtF8LUpUdA1I+HBWGeTMADtLDZJCEG7QTiTFt/WjN45VjJI184ughIOXDUdEXSQkADtLT\n2Cs7fHeU+WypssXKQKqbtVOVST58AEiK0REXkjMlLQTAQVqMNjcPz8xdK+mqbLFCkHRWJspe\nJVarZiwpqhkdWUHCTv1cteMgSXnnk9yNXs4xG3ZRrZWMkbR2EpBUswwmkKhPP9eROkhi3vkk\n2GBfhDFSpWg+zwLV9DBmsVOApJ2sU4+RYpdenmtwkMS880mIgb6YBJLqRzOc2VUVumcYI+kn\nvTXfebxLnw6kjy9b/oQPYf/2twQkpUPjMwpGsJyIIp1GNYqWLUTTGbVUTZokar95JIrKLn02\nkD62/5a/H3g/Om4AaYqmsHVWBBJOrvR4lXcqxm920QJ7IUhVGCHRtE+fbIwEQeE+20Ga0OIf\n++18O0gkvSrfSTVrB6UuAlIlR5to3qlPOGsH4IHBqQikiZqq+WE66xiJ5qPJ95HCBtK+Y6QX\nqaIaMKLmpXWRPUQbuX2woUDahkQbUAik/7tbViri6Mtf+WTLf+EsKACTKbMkO2ynqNKpStO1\n5Wr9wKg6l5/9WnGdVCBFV3Yf5RGJASmZCm5uH1Gy7HdHnE02dqiDpT6mrtZzRCLVpqL14ejL\nWsUObGNEpA/05wUg0aPhIzolbJjGSOwetgDZ2hRMPGpd3ib8CpAaYNTQ5bENAVLEUQ1I4AEJ\n2Z8bg8QvQMqTdOTKBmOo2x+k6nDU2OWxjQASml6ovrTDr6MSnCWKVwxSswEkQT+TRFqobsuL\nmq6oUgml7BsxjzOFonUY7eDy2AYA6QN/+AD/au4jpY2CRJbe4VRYVOPpyis3pqjGQKETZUwo\noZj93rN2NRzt4/L7qKKCt72P9LEtYQg7uL9NQRLXAnELhpCoztN1qeKiKgm0iXLGZyRnv+9V\naAVGu7n8Pqp8/eUWMkQkm+VbVv+m0bzPxiBpPH2aJt0V2oEgMTeRk9k3WHgUC6+ixRx9MndV\nHaRGIBk8MZ9wdSwLSPoCHAnSTDkC60EqVONMEreTF9Eyjp6OGa/zcZDagNTAFRk1A0hiCvkb\nmZ5dVlZZNGtbmcXsqx/O4JrlIVqC0eaXzMpTB+lwkJjTIpDK77XK38g0lba8vFWB1HjWTgGS\nnSPolw5ShyBx58UgZa8HhQKI38jZMmgzFkXzlm+0nUAy3zwifukg7QaSKmSI5/EkmfyoBqSU\nQ2u/H4pcPiu+zxjJiBHjmD5G2g2k3Kyd5DMyAKoqK3SqQMoHjeVocexIHt9j1s4UjiTP9Fm7\n3UBKHxYdkh5AicRLO/lqkNXPFDV5bgak9XDPi1aR6TE6xOX3UTU26oggEU9ftqKveTZVQgdl\nmylqniMZpO34ICCpw9FRLl+jKj8yaGzUg0FSORx7iCaLvubx6byYxBFXVE3RNByNBtLzwaNs\nst1dfhfVxEPsxkY9FiTFV3dWBYIk35DNunc63U2TSHdhNxhIz3CUEX2By++imnqtirFRDwVJ\ncLptl8bzjSDlJJuApCrwECCtV3UJ0de4/D6q5wZJFzpIesY7sc7GUdFNWhYkiSRVgSPRVtZS\ndRsdCaKvc3mlfVHhIC2mvQiLzoi8E6tsHKXvo/IJ+DGSanWRlNEAIIVZBka03uvbg2R6b1fm\nBGOjHgpScomCyViQUsmkgwIKwqxd6RqhISISnKzDoi18/nMHkGxvklxPEY4YG/VYkFKL5my2\nnVMMUjJf08qGrEVXoS2tlSqa9IaiJV7JWg8gyWZs1INB4qzMQ9dz0vnmAtJLQFrP7RkkcvNo\nE836scFVHaR9QdKu+bSJZrVtIBXijnIqLmrSmoBE78E+RRu7cQdjpIQZG7UDkILDVBGERQsU\nU2jwj1EYihVllAWplNQGIMVrGW669wofDRKatfsla0kpY6MeD1JwGI3rZObcSJVtzhilDTtq\nvFOYT6mZF0lZPUgUI6XPHwNSHpcCuwsbG/VwkILDaFwnM1dAqpyZW0jmgjOr8E6uCHFR+dOO\nAAmHI4vPv2KMtAs3rBkbdSyQ5CTbESVIqayWs+DZ9ddxjMU3p3Tnpa0SJICR2ed3nLUrIKFS\n1NioQ4C07W0GUiZWQdt25SuUUOMOMcslEttaqwNp46jc59WWEK3ApUVRRwQJ3AdNen34pAdJ\ndkYbSMXRIX1uXNRoKYY6D0HVaCtGu3inQrRFmLnu80gqjiBurNh6RJ61I58f6Tm1CRyERau7\n+8pYBqTCHCpAenK0p3cKoi0AilXrDLXMECDdTRkhUnPl5DCvgXcITg4jEVnJV2LimY1AImeW\ng/RF0ffv+3onFtVxc9CkOmqaMUCKXDW6xGF8mXU7TVADeUjOC3dWcyRbeoykNg6kAqkviCSM\n2oNkiDvWWfXrghRxQl2KA4l17hQY3DENIftxlJ610xsDkr3An99THDUFyXjh5iDl84Z+QDhC\n12kdgJSvjd1arbXDBbzZy/z5xEjmqJF3Fg18HKR83sENYMevWww5M960g8QdyJTykWaX9aXM\nsxllQuhEK0h3p8lgVOmdUghSivoYKZv36gUsI0wMgifjEBbvnKPEAmKMNj56/98EElaTPZp7\nWtCSj6iqBunpMrlwVO6d6YGQVvSYNeWpruKa/UiQNkdlGBF9H6QCyelOKbF0lDsOdltASkZS\naOzz65ycIXP9GGnzmDxGRd6ZHwod/86GlCW6im32A0Fau1u6RNN9rZJ0BRdhQk7bbtulHVZL\n1UIFkjlOqWbtgr8owpHRO9VTCQ5SG5BYzwFBKu9Cq5vvCpL2MolXqwXJnL2iAaC7qDDSeqdp\nNs5B2hOkEAN0HMU3g8SVDbKaBiS9K5eClLgz1hIk7C1KjhTeaYRIJ1piDtJs9N01FU4rrbVL\nKQrHakJSTjsq6sxT3xQk6itajDLeWTKpnRUttuuBxPiYyXW3ZCitsPo7rSkcKQOpdNaOZhx9\n1JmgyriKniPZO8sQyohW2QVBin/WxeS4bEK0Iw2SjtVm09KxiU2ECmvNnFNlPcXAEe+dNRCJ\notV2RZDiXaYIwHPEL3+IVfV8TC+6IQszLIc3UgXOAe7IWDBivLOSIV60iTlIszEizWyQgbvg\nJsuR5jaVUNQG9gqQkG+ANQI2juQnHmpc1UHaC6QSjlJ0UJCYs5PiiaK2MNUYqUoVu0ZYtaa7\neUS8s2BeLi/a3hKqpiUSyq7aUgwBEutWTFq6Qw5I+dkHdPTFIJkHRpwq50dPkKwY7QDR3V4N\nkm3Rnrar1hSHgxR8JsURv3aGD0l4G6aFu/MB6UCQqlVFR3p/N4ajHQha7MUgGZeRG7vqcJBI\n0BA5SpBE9nFZ4ZTPz6lv/W5ByoeqxLscA0dKf9oLoac5SA1BCh4rcpQDaf2cqjKSSF/VgTS4\nqM2tdDVT4njKkR6+9H3hKO9QCKJXv0Wo3OSaXQ8k6DbzrAIpJODyJcJsNjQDvN0JSJmvgIcD\nJN1zwyjjUTQQDQNSqmanHiOlQGIwgS4VnytUmdWNNpM+2gQkKt8WpNUBUu75HB3lQGKu5iTR\nqnfW7/Pu7zRJeiljVx0NEjdGir2G4wjOFyRBir1P5ki7LK7EIvmGIAEHkN1znWRIehs/JhJE\n63794eUgmQy17wAgMbN2Oa8JhwiDfJUZHbhJQOJJagBSLN9sjIQcQHTPMFknOZs8s8CLVrqt\ng9QUJNFPFCAx2yqQmMNjgMTM2lEHkNwTTnozvpaenxsFpIa/kIQaeTSQmDWn7JkZkOL4k5rq\nGgokYowD8O6ZvHeUn+QeBqS6URs01M7DgYRN9n/ikxik5JWclE2So17GSNh4B+DcU74Hq7xT\nZBojaX3Z19rtAVLsw5N003Q9QC8DF9HMpe1RpKgAACAASURBVFwiS/msTmbtgBkc6bvAkeF2\nq2XWTh2lHKTmIHG+H19rTXgyHB3cPupBsgSuTu4jrWZyJBYj45oFg3fqr/ccpNYgceOTeNSy\nfF53Crzc0tdoUQZkR66oLY3OWmrN5EgcRvZ1Pw5StjM7AIkd6EcgUYAyIGVzZ+KdfN5+ICkv\nQzezORKDUdHiOQcp25mHgwSIuSBI2tI+zepIEUelS1At3uljpENAEjhKLHRYd/IeeFN/x5N0\nO4IEZGEORpDsjoQxsl/QsaJ581m7A0BiOFo/TWTWbksDE3Ci2m94ki7l0VUg0bAKRbUgFTgS\nxMg2tZAQbWkOUhuQQniJL7HkE1SqZZY4tQYkUEFcV/UYqciRAkeVEEHRpuYgNQHp6UERSKqv\naDmBZcRhMD1IcfZpkPLklznSylE9REC03LhaOUgtQFp9irmwy7LAp4gvEtuZGiQm+wxIGStz\npAWjNhR92r2zoFErS1hcVMkSRWUrczBIgSS8X3MqTjbxu9uYFiQ2f7ALHVWIFjrSHaMGF3RY\nVGPKZkrV/1VFzZmiqCjF0SBF89jKgERddmL3trIqkFKzdkkrdCQEUZtXLWS8U9k6xDL136eo\narMUdT58jCQgYTgzsa+h1YFkFH0KlDrAz9YQPUT53cpWMdafWJuims1Y1ONAUj2+kD9T3tXS\nasZIVtG7QvGjAFqIjE8bYO/UNkZZ/UUrKWq5GYt6IEiEJEuTyiDt+rI4ZdG0U/DSC4/KH/JR\nhyJrFj9j36q3kq7SFrXejEU9EqTbQgIAQumC/AXhY/NokGZ9VOJFi5+WM1zQmfP42Roiuf55\nyxW1jRmLeihI4MGI555AR8YVmTC2bBW9UURRVH16dYRlivq5ObmRpI2i7xpHsmXBF7XeXvU8\nVpkZi3osSMRCnMm6onQ4eRp7UOX1oKiaopWBtPSgnaQwMPr+XeVIuhwSRW1i1aJsqR2kxERe\nznBIS6SKD+qySj+DWKI40/pb3XyxXyBF37WOlM0gWdRG1kQ0Kvk1QYquzopA2hInz2NGVblT\nYFEFlaAUFydZYtw7kZdrQPoFU/RcWqdzpIR+uv4NrZkoKv4lQWLGOQUgUf50IC0biVPW3RP8\nxT4mPbdDUWLQO7GXm97N/QmXqFY5klTgvkF62FqFK4LExIhljz0gWUHatlIchZgV7xZ0DQVe\nmojpxCxH5HYRfGKi2JFSJR4ApLtV1T/VHAOBtOEjcST5aoSRbowU0qY4ilcwxenrQCrpYzrR\njR7gK3KkXIkHAeluFwcJgSBwxDqrniOMgCYxD5KQvawknXQrWsEQ3y7CD8K+YKF2I9vvlp+9\nVTPt0jlIZJyScl3JWyOKmi3mUYFUxNH2LFbF7aKwi76WwfKektrVTHW2671zW7syZizqsSDN\nhCO0jdIy+4CXA3fP5As0ct6/Ztieo2feFbeLwM7oNUH6Jx4MBR8PpLmWJWNRDwbpbhN1Rx1I\nOBDJIKGzTE6/Jp3grB1XqoJ1sp/W5QXCatT4dVsqkEDhVcUdEqS5iiVjUbsCaT3AdfC2L/yB\nJ26faL6xcIHbpxZqFyg+e8oAkrSmm3uHah6kgsYYFaS5nCVjUQ8GCQST1MQY2LemxByF7Vt0\nFpBuDdI82RVDV2k5Elejsq/0Vj+DdxGQ5kKWjEU9FiRwWaZspigETWRkRcZADGqVIE1MYNQL\nIZfXcCSv6RZejZ8AiWtLc/3b2QsX6ivASTbWkSApbEJmOWX9FCRuVAjtRZnVl3n789yhPvmn\n0TaKmGMPjPRSXDFr22Ios7Z9cUbH3UeyBAocYWYYiLAQ2jvXxSPwjWQvcTBD7LibNDB6mPyL\nR4zqbIo+yfq3tJc/OqaKRKDRLEXtCSQ9SWRz2ceCFE9OaHKJihrnZ5VSufxmSYqSP8DHPhVe\nUXVc/5Z2xDOYSoo+xwKJ3k6FB8Vuxwf4MDTHrnMwSHmXB5ahKBGOiCqtfUnVcf1b2jEPM58S\nJNkvtf1OzuU/ciltRlZ/Wy+V8i7PYSR1cfL3YIMq00wlVSf1b2hHvRXgjCBFs9iLqTuenEtA\n4pLmNYWiBhnwR2NZlw/2nqUox5HwnpI6js4Fkg4lo+jBIG3OPUUXdqqe5zjC95toclEnV1RN\nMsYSXUVAetdQlMNIfE9JFUfHgmQruk71tCARU4JEg1nYigTkjZlSx2R8Y5JpLNlVGKTcwEjF\n0TxS8FCKGltdXdRLgKRoPe6iUAYpThZLwQl1ktMNq+ssdAp/5xWAFChK36RNc7S2ans7UNTa\n6paingckGZgSjvCKPJ6w+BjRUmumDfaJsBZoBQlRxKZcdiYxgq3a3E4K0pxgySh6NEil1+48\nR0COCUgqkCLRsC1kyRnFgOXjJ8ZIfg3XsjvBEcjaQTLaSUAqNAkkdMcnSs5sRGoySOqbSBwH\nAkhhXJTh6F28eUQyPxtIu42Rgp0cJEoCPZgBKRKTx0hP+VWND1c0U7lzBRA4kNDsAkkSNgNH\nOYgSDVBnp5u1I3ZmkPKuz/rzNr8Wywna0V5yLI5lCZLYLzch0iQn6cA5Ujji283gSHoPHYfO\nUtXTgoR8lXPciV1gsKGQ8hGJA+GkKTzllAGJZ4KdtcvMdaMoxoUjsXJ6RzJcM50fpJn2nlF0\nXJDAEbqZvvQqWeMAiyqCJGEUW/6WEQQpDkfaoibN0gqXAAn3oFG0N5B4dsQ+JweyV14ZOV1R\neX07Rak7rwAkE0YOUp3q2CCF7gxdawhIZpDsSxSi2tDz7RhlnkeSOLIXVbK2IE3c8LJWtMRq\nVQcGidITrUwoA4lLzHxUWaOHXHAsSj/Y9+QIY1RfVGiGr5Os6FPL9gXVJ0hrdxpFOwAJ+L0Z\nJGmuWuAo18tigkRtCinKgvQwyFG66JqiqqtrFlVeCdhEy6yF6llBSjHARiQxYbqX5QRSbdQU\ncfMLeZACR6li64paZXuCZL060KnW2eggra4cdUmaI3qbJ5/SmoKvjR0jtDMH0oZRoszaolba\njiCZTlGrVtuAIKGm5DliLIphudNw4MukkIoKrY6iPEhPjuISplpFKGrKVG683xjJCJ9StbWN\nARKa855B24IGpm29tj4DEk5FTmHv4tIEiaJuVo1RBqQ7Rd++qQuYKmrSdG6836ydg9QQpNWg\no+PvN9LYIPbQRPJKICSrubsbF3U73ACjNEh3ighG4fop0aZWR1L68X7e6SC1B4l4Omhi0trh\nGLN/+yiclZkDTBV1Pc9OUfJ5JAajb99ijs4JUq9jJKtoRyAFN4kcHnc3Bo1cyt2iJFi/7DsQ\nfnkaKboHo9TzSLHND4ziIp4TpC5n7eyi/YAE3YQ6fBKkKRblDhO1dV++mFhV+TssaLZbeI6C\nB+mJ0TfOvVVef9AYqcD6vY9UINoNSAwWYJMJLRN9STEQFTnC/OicCKhafz9i2dSD9JXLhlHZ\nPeWDZu1KzEHaE6Ro3/Yx3s9FJTSCjemDaukv+BhP5e+wxPMLSpAeOT1HR3rEIxvH5x2kHUHa\nNuiumSSGJ8FU0giWRKTsYCnCcyHCjtGnaoz0zGebZCjmaCCfd5D2AIm5kMs5+gxJIqLxeYSj\n9PQdPXbL8ZPCSCIwgLRmw07WGW0cn3eQdgFpm2HQO/oMSMpXGXOaujM7U5B0P8yqeMyIAylk\nauaIK/44Pu8g7QLSRC3l6IQ1HiTmZBLp5AsoqPqpWV+qeVqPAylkaQ9HbPtErWqSFMxByqY4\nGCQaJhBIshvEIE1AFKRhTtNNfm1pYPBoiRF9S3chR1FF4pGnSZQ3Bymb4liQQkdzHMnGgDSR\nKjMyMNRlSzfBRQxpkEoo+opFsIm+GTkClSdHbjhVG5IcpGyKQ0ECFCAk9AGjDKR1K4MUDB6t\nKZpR71jDUeIbhwGpniUHKZuiE5DICrvMmSCREqRpu6SLryWlXIDfZ39cT4fPMnm31X8xazhK\ncSSApCEpkeqWPlxm+4LUtLgDgWQhiaRAZ/BjpOUzy5GQGSIg8+N6Ooyet5Ng/Z9WPjriys6v\nF9Fdy4qpbunDZbYrSG2L2ztIckhKtQFJgbfYWTtGOJcXRoADKUsRvXl0xwjm9iyqORzFLYBM\nmgut0QzfebZypm1PkBoXt3uQUEeXgUTSs/kyyum8KBMRSIppumg5wzvJ7VFUC0YTjbCcRQ3g\nILUSTaY4GKS7bVWmzi20hBGkiVi0OzojhkL6cb3EkiG6wC7u3HtR9RxFpadH0w1g1CfmIGVT\ndAASdw22bXAnogM5kChH8ESeVYYK6cf1tCAx5b4X1RCOZOyJMtPwOo9KeZ6PkbIpegAJWnCX\nhOOg3fHXPKO3UpNvXJaK+Mf1xB9riUESyj3fSjii8/nwKNsAeks0js/aZVP0BtLi7bzj8Ea+\n5hk5MDBKnPowngry43qBFJGj9Xii3MbhUdweOHZXg5SwcUSveR+JNcyR9WuFfncmhaL9EhT4\nx/VCwElw9EiSKKhtsk7kKAre4/i8g7QvSLP1RiIVxeclpOh+mQny43r5cHS3ZDkL78EKu8Cx\ncXzeQWoPEu8iYhhJiMZnymJkd4IJOtedGyDl3o5acPMoLj8sfjg6js87SM1BEi5a8lN2/NEl\nCf6bA8mAUR6kTM2fGOm73NIU4/h8PyBlr3wGASl29AxHYr0DR2GcBUHiIt9zpwGj981wQjLZ\nLdoajtRdnmqMaN+hPp/1yBJRmxWopr+bdaJdgURIkiutAIn5yLTYelSMLcwChgRHZLabt+2q\nzuCd2Z4OdiRIlnL2A5KifUcDSTfXnUiKgw8FiD3zwZFwlbbdev1JeJE4es/WAIyOTgeSraAO\nUnOQMnRgk9NuwQenY1dNwLOk4Q5cwRCBFKWmS1J5g5MMDtIe5iBp+0BK/NwXg0NPoyAJYMB1\nQGmQnlt0SSpraLKufozEmYNktnypBwFpkpycqTFIHx27243umektXqIJbwkFQLahEctNvD1r\nuoPMeWebaNMzeKePkQrOyZZ6PJCECm1ATPFYB9MVPY6D1BmOAD4bItFKILT6O+Yo8Joyeu8o\n10Q2t9SqFlnbWbslVT8gNRA9GiQIhuw4zyMhCUgJw04E0hQZVaajnDsk8SKG1Dsb3t+VDh/d\ng002UeabRbaXO5K5jGu1HKR2IE0YEaFPYhq2pNGhW+K8NEcrSNwDewQkGJOU/s4sZch4JxVW\nOuyrHclM+1YvB6kZSHlPf1iUCAOFThdB4hwRQyFiREEK898hl0xduSVBOe/EwlqHfbEjKb9H\nuDMcpP1A4hPGvEkgLRFpisMVIxug2PiQHh+nYyQ4161wJX5pnQoksieVS1a13BykbIqOQBIT\nGkCabjN/4RerQigyHEUgQVW59Ote4Rf4NCDRPfIpCtVyawmSj5Hag4TWw4npIA4TmnRgSLqh\nvk2AhKLL42PibSYIJIInpDou9xqOmEJkx0hMI6TaM69abC3HSD5rtwNIuQHGRMLOuhNtoBQ3\n6nNoK+zGl2n3j6m3meBfYOHKE3vUkgT8chipaG7WLtrR4xipYNZOI1puNapiVfoHiS88dtEo\nrMRRJqSIQIoC1GMnHvCsGDEcPY+SX2CJy8OEjIDRNz6gGLs877DTbiOPcURrVOXvqhFASlRo\n9b4IHM4tlz1ojERSh9MiUtblQAxH953kF1gW1Th/mi38Bb5akLK2fZO0t3FEK1S5PlKLdgIS\n6/iBn6zThl23SI7Ro7iIo6P1yu9n9uEIpkyAI/bbrrEjLQXYzefzEbFAtL1dGiRcg4mYVOe4\nDZQgsRyxY6MFJFXvRGXFc95xRcYCSXSyGtEd7MogYWgoR2BSgVSaNsEKEiVw2RA4SkzWbQvp\nVL2T5ChV/za2L0gJLysX3cP6GyP997e3t9uP/70WpGiN6USnCZhKw6Q3vLklCOeow9HdgqrV\n8q83GWqMdAWQdpm1+/fXty+7vb39/VqQEhxJEwgoot1ogIsaSs9RaFvmgjFtmrcENXekUubz\ndhGQKkRFkH5/++OLotufbz/2BylJEkww0zNSN2QDaMsJzzNjjqRoRIoqO1K8X/W2rbEcqTlH\ng9U/l0IC6Qui7d/uICHf50iSQOLBo2IzcAOOI/ZJc1rUqAi4LHCH8q11gzlSa45Gq38mRR8g\noQmCNZJA142cFUEHL/CQGP0rcRSTFBVVBCk6oH3546kcqRfR7kBaLu3+ePv9NSDdDbgjiUjS\nrF3EGplgoCCJHBGSmKIKIJFiWl6ieipH6kW0O5D+/Xh72Mc/+4MEaQGf5IupEKBgisfnG5JM\ngrSMj94jkviisqUBxXweNLyLuPkSoRJVnY0j2h1It9t/fn17+/WPfws5MoBErsjCvjRH2zlw\nd1jAh0EDHK3MrPMM75QkqahMeTBHYWmdymxdLjZGlarSxhHtEKRKU5cOuyIliT0zcexGjgM9\n/CQffrtJIEkuKpMrx1G+4rj+Oks1R7mq1sYRdZBoDEp5jQEkYCj6MBy9ixzlQFo+GX9iwkEa\nRrVq1m61j9+Lhkna0kUcUU/lTpWdSgYJQcNx9C6/t1tUBZHP/EstDtIwqk1A+rISkpSly3KU\nJwkluUWHV2M5QgEpU1RWddtj/8UjHyMNo1oB0l9vP77w+efH25+FU+DK0iGE5CunyMIBnATP\n2gGD0ID1DCqOZNXN7Bz5rN0eot2B9Ovbc77u7dfCm7K60kUc8SDJR+j8Hlp3BE6A0MQ/vicO\njkBRGasIR7JopY3j86eqvwjSCk/x6gZV6WKO0C0l9sh2mIIUlpduKcIJIPw8Hob9hCTlOBIa\ncstCxVEUUU7lSL2IdgfSj/XS7sft73tU2h0kbpkCXgjELVolgyxycyqcEaBh381Q1JBrFrqb\nR1FQzfVOBJ7OxvH5a4D0T1jZ8Pb238eujy9b/36QbfhXDxJmAPg+jFJkGQ7GZ6YkRSBNZIkd\n93KGsoZcstBd1kVfBbneicHT2Tg+fw2Qbv9uKxve/nhytP630LJu079lIN130IAUrRcKaaJd\nEkjP/xFHdGVdaUMGjBTDIytITHqdjePz0nC2lWpDa3pDdleQnnsidNAGgCZa+RCBhK4K3yFH\nKB69P37aqLAh1eFodpBSoqXBN63a0tqvbPgIrDQASVwIxHMEI088NfHcIl9zT47eIUcwGD1e\nyVDekMtLVBUC5jHShUAqrmpStanVgPTHejc2AmkdIokg/d/dVFzeVgKWz3gvSTJBW/YzMrH+\n+/vPn+/3/+4c/fy5bDz/smcpbX0XsbqqFnG+Pi1z6MWkvhvUCEgbR2h/CqCSiMSvtMNfUVEK\nkGzGaZjvjiUgbQvs1s/Le/OV34VsbbZ3eseHykWNGuSEacQXRJ47In28/e/H2z///sAvP/kA\nH3YCiV610ZtJfLNv+6J8A0fv4beVF5JUP6DMqgKMeI4Uqo27fP1ib6v6NB8jZVNIIH1Fov+8\n/XX7F7385AN+agMSR5ICpEhs3RcdC3dif4lASnchOBTXJstR3jUcpNUactQjSH/d7x/BS7uP\n8H87kLiLO+SHlK+kj8YJNpAIR+/gbVuyFC5qsMRlnQ6kQpfPlXdMkLpXrQDpt7c//3n79fY3\nAAlMe6cnG2wgzZSlmfhLRJIsGaO2Tc+xHCUM6ZAvzxCOSkFShSzjac+D4/j8NUC6E/TjPtew\nrfz+yK1osK9s2IxwRAwBhpOQ9LHMNtF9n/g2cCSA9NgHfvIoca5aPDpUctp26jg+fw2Qbn/9\nen+V0LKqocBspZvg+h9qMUhrOupXEUjwTuy6QDX5CB+T6VrUZc+2liHn0lrx+Ij9tGDj+Pw1\nQPqrFKAqkHgL7kMiV+xYEkj0TqyUD5UCRYUc5SuXNpGINCoO0kGqVZMNv9exZCod8I/YUSZ8\nFKCCHWvCB+9m5YiSBIpKwlGlJQOSbgJEsHF8/hog/XYfH/32Z/HbuMomG+CiVHC5h2mJFtwR\nkZijTzVHa95xUQNHTaZqJ35+LRdzshwN5PPXAOl2+9/9d13efvy5N0hTymbOt8gxRmZJSDh6\nT71xK3mf6rm+1rK0TmHykvLmqrU2jmiHIN3tD7rWrj1ISY7i67dwzvoJ7hM52l6pKk8zRASj\nok74B2HztdPWny1Gc9VKG0e0R5D+/uPj7e3X/+wLUoYjEnaexvo6gwHDUWreOwMS/mHlfPWy\nJjRRpfw4Pn8NkP66U7T/GCnHETspzvp6HJLiC7v39ONHKZDCVR2fe4ElIlKF/Dg+fw2Qvq7p\nyt/7rQZJgkfw58WEYxNZpkcmGp4cqe7vxGMkMDrqECSYXn1TwWK3tuvhNtEdrDeQ/n5GpNLf\norC/swFHgwncdV0MbCZcjeMIvAsy5wwCSHDS+yiQ5DzRCfxNhUq7NdQContYbyCtLP1aurRB\nUzoZpDm4wQRmwddT5U69QY7IO+zyHBGSlr/k3lErn7KNkeRcMXrRTYXKUj5FG2oF0aZqu6rW\nztr9+599Z+1SIK2flz18T8Z7tnzDmu/cDSSsBTJaPkT3YEOuVb4lNZHMEZ+bg7S7at19pP/c\nF60Wr29QlQ7SwwWk6ChWYHat+QaOfsEcJb/XZ+h+zw90LQM8X+NccgpTlztIR6pWgPT7Rw1F\nWpCQKyIX1oDE7cMgsRxl1uYwHFHvTOTP6sr1V1sqL3TEx0iHiIogvWStHembdROgQ6/+ZpJe\nAimMkMJb7KRTqBbgaA1HNyaZCqRUEluXp7KCR3zW7hBREaRXrP5O+/RyNEGSDNIyubA9zLeN\njzQghe1t0vsFIGUdVenJ49zyuQZI1ZYvHWQlmtzWk8RVeZmkCw/FZkGiWls4mmSQFNc7apDy\nUkobx+cdpMYgUUZibjiQNgSZKgeOAkhQXCoO2lqfJ8cNGcUtTSX5XG402Z6T6lcRzauWNPMQ\nIMWQ0F1iGqHKMUdghZ22HTeO7v/wCD5RodxOWGwH6QjVonbuG6QbCUhxtMlwxDUIBkn9hgZq\nZGmdqs8VXYSK7SAdoFrW0J2DFMWkGX2Brx8ZjpIg0ZmGUo7SGUWmSSmC5GOkF6meE6Q42tBa\nopC17Ylnz7Y/N7jomw6QtAbuwWobPkI9leq5YZy1U9o4Pu8gtQOJ5YjMLuObsQxxy8fnn9u8\nzHaXByQ4y6BteFSLXLpQ/x1sHJ/3MVIzkAhGaCKPn4tAn2Z4zmK3JQb9At5jZ2s08jLip2rm\nHC1H4qxdOxvH533WrhVIlBJuIo+kYdwVJ1mDEHgq1tZ84C2q4VpSDZIuj6X+lsR21RKHyYq2\ntFN9kXQFUgajacIvEFpkWJB+CY/FhkRMGejeZXREsVCClK8zqr8tuVG1oEB50abmIO0B0hyP\niTiQ8IDoYejIJ+Yot6CB7gUvUTWBpL7Pi+qfEy2yGyhQM5IcpGyKI0ECVMxyPCI3miZyhYc+\nzZ8/V5A+FSDR3WTW29KQCY4kd3aQhlHtHCTc1QJA8EiclszovX/+/GQCkg4k9OARgFPVkIKl\n/NlBGka1b5BkjoAbo8UNbGJ40vv7Ty4gSYEByXIvI14PjwfSS8dIZTk5SG1Awj5GoCApNCBN\n0UwDfLac7+ugyr7Ue8t2QJBeOGtXyKyD1BwkJrzA3ZE/ChxNMkfLWVwhHn/5d+PD/LP1SdSS\nP3QWR0p9WRSLFtuVQeKYmLn9mXmJEJDeGZASvR2FI1yAzFdu2onko626HOfgIB0ieiBIeLWP\nDqSAkwwSDkhbD3PdvWxzrwmCg7S0p5R50dysy0n+DtIhokeCdJsoL2SpDc9RuI2Ez5jQmoYF\npNDFTHc/d8SjI1yEjKeUulGrLqf5H+BIhQ3gILUCafkL4UmRlFtp9+DoJwpIsRbI/7mDGR3R\npNGpzMHDvpE7AMln7foAKUwqEG4QJmBz+YhPnrZ4hC/sNi+LTpU4YsERN3EJLXYekDoSvTRI\nwpgHHQMb6MO8XeC9r08hffIgcTenhB+yjLiIOWpA0lnGSF2JXhskE0kwPRAITyFFa4NgNmic\nJf4gLKXiRg5S1SNBOnzWrivRi4OUIonz0sh3p+WdW9EtJNbjF6D0P6ycBEkx4cAePJUj9SJ6\ndZCkGe05A9IWqLgREjEsbPphZVybEo6Ywy9zJLls2jCqLKotKjtILwQpCj00+bIZTX0z+WG5\n8PZHhZHalHAUJ6jrcinHWFUuXabcKdEqOYuo1Rwk5OZpjqID9F7s+ktI6MwtFm03j9TdnqpN\nlqNdQBLzjFTl8uVLLoqKRTKQ5CDtDRLdpt2z8bDt/twC0s+Fo3BVSM+7G3q/yXYQFw/u6A4k\nOVMH6RDRDkGi21H3gMiyHng8U/4YIf1cOZqYpT3LBveTRwKtioZU+A+fwEFKSall9Krl1jdI\nLB5bXzAg4SNzcNDwJrvPJEhPY1aoxrjiHfm1ZpaqakQzZgBpwDGSUUipWmFdg0Rbi2DDcoSn\n8bYPG0cLSFMCJG6ldx1Ipd+fLxojjTdrJ39L1KjWWM8gSa5LKCEcBTwASICjB0gYOZxtfO+o\nHqRCe9WsXQN7raiDVAASvusag8ScIIP0uYC0eZmCI/7KzjBGKrZTDbZbijpIRSDF13EqkLZd\nKCDdQZLzTDwHWz5rV2wOkmQVHF0PJEgSB1D8tRRxtDj7BtJ9kdDPT3CEmLQkKN9r43jnQEX1\nWbs2IDGLGcBe7msJMAb2hoB0nwL/+QlxnGCfMJPeqbKhZsrXxm4O0jCqfYMU/z4SmpNjPD1i\n427PpQzr75f/jPEEN49iLd7oQXekUUSvCZJqbASNCVUbR0+QOI7uqemSoFw2+PjRjpSMnsWq\nehtH9KIgcQsYYNjhHAi4+LRONWwB6Z0HaV1aF+QoKlwW4fANHMlXTGn6Lk9Gz2JVg40jelWQ\nmIHSxpHg6Rg59Ktin0JAAj/Vsigk9EmhUENaPFpXf40lkK9Qtdg4opcHiUSByJWDYda25/me\nH0SMwE8e0fyYogkgmTxaV3+NOUhHqw4A0syxkwQJrVsIz/Mtv3EpcBQzqr24gw3pIHUvel2Q\nwOWSEqR1oIJAWn/kEv96GRodYV7ng+Fp7AAAIABJREFUrIOCg68BKQ11qWoLG0f0uiBFfk1N\nFlhB2n53ebr/2mvMEdRBsrL+SjNpyJYcRb0jiBvpHcfnHaR9QDJzNM80IEVCEUczlk1wxN9H\nUnm00u1JE0n1dZCOVx0JJIGjbZqNtc/tAfOVI0jS9noTcL7KLXGiKf+r5vHZinQ1IMk5jOPz\nDlJDkAIuCZBEzwwcPR6RJScwz5ODuY2UQYjNl3PqCKIEiStAIodxfN5BagcSpYUFSfbMFaTn\nTVl8yhwu68AZOi+PKLaQVAqSTAfPkZDFOD7vIDUDCTqseGs2cc2D3xyETiL3jsI5Knc9BCRp\nHYcti3F83kFqDxJHEprmZtwmTNk913+jc59LVC0g4f0RSFaSFOk0XS59hThIL1UdB6QVHLS9\nmsARDkif4PTtvQzxmUmOKEnskbwpU+ufNIk0EyUax+cdpGYgEZIkkGIvmuOANJPZOu5Mfs+6\nO6JXOJKzhhEJf8fgA+WqdqsW5crrIO0EEiaJPQUlRQEpHGZ/nzxriVxLOCoaI8laBpi7BIkt\nvYO0G0hk0ARsjULAtoD0WK66iYJJb53rTblVDtpLtTVtQ5DWUg0OEl98B2k/kOSJBRqwQkBa\nFtk9Da30VvneliqV2Lq+tBlI+BtEc4KDdIhoXyAl4hG1eQVpWRz0TEgxyjufLpW2dywcmbpc\nzZGDdIxoVyCJbsgl+MSrVR/pGI4OAkmT2NDleo66BMnHSAeBRFNuUQgkAA+YLyBxGB0Dki6x\nvssNon2C5LN2+4Jk5AjeHIVvanhwNBVxpHRS2xhJmdY68GqsarFxRC8JUhRoeO8PO8Hh8KaG\n56+KTWgtwwadopE0qbaF6k3UkKjCHKSjVccAKReTwD7E0fpy1fvU9/q2LSAQux6HaL6Yj6Ki\nzBuZgzSM6mAgCWN1uG/9C972fQcp3DxK+Fx0RO+e4IkPVXq1qM4sGY/j8w5SK5CYNUGMt3Jw\ngZ+f+OIIvrXuluYoWsKnn2A7EKRdLhgtNo7oNUGaJ3w1NovLTMluFJDgg0cJj0uDlHbVo0E6\nWnUc0YuC9LAJQQB9FVysISeGAen7srSO4w1ZEqQMI8eOkY5XHUf0fCAZ7Ms/n/8//qL9912A\npOXAz7vdf3n558/vd5Bg8uUMTQZhB5ZPlzNxIC/hdno7+FfN8SQ33BMvCg8/LPb9wRESIHFp\nO/+ZYKZ50szJgeVzujag4JaQdapv5F5EzxeR1KWb0DTDFNwRMIEB2a7svt9BWuUISJO8MzZy\nlJCYeYsQ4T1fcVz/xjaOz5+q/seDhL7OecMHQ0BCHLHMcFqsxRxhElP1cJA6Er0sSCwr1PnR\n1vYLfXeOoB7j0GqQ6MXcWCDBXMfxeQfpxSDNDEhLQMI/vDwto5kIAwVIRMcC0uFjJJTvOD7v\nIO0CkkgSx9H3B0efzC+Y35bYsGyiWJF1c0QPVMjUBpyXr3YoqiFtJndQ3nF83kFqCVIyfDDX\nZ48R0vfvyzK7WJTigiFNFwgVBSn05kioKg7S7qojgARm7YC/Q+dAXD0w+v6LBBJyKopDxhKw\ndeZIuJwO0u6qQ4D0NBqIwn/4wuzB0fflF/pmomBarRDZMCDRgqKtvor6ctGrg4Q5IuvuwO73\n799DQHqPLnCYGTt9aw0Lks/a7a06Nkjh87abchTf++FO15qcvi9HStasr6K+XPSSIMUcxAEJ\nRZdvD47e199xge60Jr5R+iztVeqdxmx0oun8xBxvZYVJm4OUTXEgSBMDAsfG9vnBUfhBJDbx\nrfi2jpD4sVu31s4o32rWrqAwVnOQsimOA4nHBnoBOfzt24Oj92XG7h2nXz7eAH9WjrjkG57Z\niqjl1w+7dHlBHFaYg5RNcTxI0jzuDGfGH6/0fmC0vswO/bBYSHwrdCXhrGV3NUggzfbBQRpG\ndXSQwO7703vf7wMj+BPmkUKvIIGChk8O0jCqXYPELDdQcARfCsmDVOZLFSDlhyWwoDuD5GOk\nY0SPBGmdYAMdz3jzitG3x7iIDUgtQKoYI+UHY7Cce4Pks3aHiB4K0rxyJN2lX3Y8326ygfSO\nR0g4eTFI6lm7MulQou3jqRypF9FLgrRhFHwzhmD54bAHR8vLVcEPXUZujR+jaGWgNqWQorBL\nRRvaOD5/qvofCRK8OIsmwbdU68u2MEef7GXgWuXmVzehNg3D3akcqRfRC4I0JUFa9izhaHlP\nw+e6puF95l26kiDx9GqQ0qItbRyfP1X9ewGJ2zsvHD32P37f8rmmYflhsejU6ms6+XQHaRTR\ni4PE7/6GOPoiKXCUglDXNnHCxOm1YyTWbkI5YFalqo1tHFEHKd4Nfp38ff3l5cdPL7cAiUmp\nA6ndAGwdzkl6ZciO4/MOUhuQeJJYjgJIDEfRCnJNy3BJmX3r5m69kyiyKb5i1eY2jugVQYLv\nIl73sRjNn5qAZOGIWyLL38Pa8ZaPgzSM6rggAY7uM98bR/EIiZ4KNhN5M+fH54QEJwcpm42D\nlE3RB0jbPpajx89bvm8cvU8TBwLezDkhTxKbZs43ZIHD7z9GUp+dz8hByqY4EiQJh/Xe0bId\nZhreGY7CyY+PN/BZSVIuyZxtyBqXT5xZjqepUNl2cJB6B4mQ9Ph/5Wij430DaeOIghRImtUg\naRJthxWrvzsZzWyq+kLlv1EcpHyKQ0GiPNw7c/sdy4Wr9zggsYarLLkG3pd1tfVwevr7LCDJ\niR2kbIreQAo/B/vs2s84ICF6qBcgP4qyJHvV/k9vyEqTEibrBSRFSHKQsin6Aukb4ej+qu9k\nQIqGWXCIEOVI/EXva9ESIWGa3GQ7O5KhUA5SteixIFGSws8qr4Q8lnpzU98z+RuBxFlLkNgb\ntybb25Eshco0hIOUTXEsSIikb4CjxR4PIcELu3cKUB8gFVlXjpSukIOUTXEwSMxaBsyRMPUd\nT9yRMRJvhAA1EPEY6WwgnUX0qiCRePTYH2LNOzPTsH150oAEIpLs5lPZJRmZtUuDVCLazsbx\n+VPV/2CQZIzujEQcvdN4QgLTDOgsazG5qCjfVAYlYa6hjePzp6r/sSCpOfoFcwQ9lZB0A47e\nDicDSPoLv1M5Ui+i1wYprAkCUFCOPjeQsKdikuicOs20DC5cm3xAcpCOEr0iSD9niNFEQIIc\nfX6GX5+IKcEg8WlW0zo5baZIRUzqIB0qekmQPsES1cX7AEifYYD0BCniKFzohR0CbIvpvdza\nkCSP1qJ6G8fnT1X/Y0H6nOnNow2E5+sgt5+eEECil4PirPhi/N4WDYkzaS+qtnF8/lT1Pxik\nz+ge7LqiYYHnlwAStzKPhCT4mb7CFac0ktSXIyWL31dRXy56TZC+f/9O1zI8/X9h5xNw9MlE\nGwLSo0IYKppnIUldOVK6+F0V9fWilwTp8UOW8bNFywq7iKO7EXrAzmVznbXDKWY2bcuGLLDy\nxTyp8vdU1ANErwjSg6PtFywBSZ+Ao8dIaeNoDpjMiJLNu6JVceEzTWqwnhzJQXq9avcgPTCB\nAWQKl3XgvarvASTuyiYL0vaxkKOuHMlBer1q3yD9/P4ZLViYnth8orcGvT/e9b1Y7EcTxOMW\n7QXnlHLUlyOlK9BVUV8vekmQVkru7Dz2LHdhAUcrauGsiATE0bZECEFDQDK24tybIyWr0FdR\nXy56SZAAR4tFHKVAmrhhD/PY7bx9iZ8EpNerjiN6bZA+PwFJn+EXLj/Di+w2w2FmpnRES4TQ\nzEQhR+5Iw4heGyTEUrgPuz7Xx42IcGQKKQJIbPwp48gdaRjRS4L0KZEUHpyIOOIu0yAdGZCg\nSr6EQdWQ9ljRgYp6qvofCxIm6R2vZwgcsVFl4jm5SaEKm/Ear7Z32MxO5Ui9iF4UJETSGpfC\nr06sHHEuP0UkTfDBvi0Jm7V11qGyd/jMTuVIvYheFaSYpF/Q/aN32eMpSM+NW5QgcW6+jGsz\nqVMacjuVI/UielmQKEm/MPFo4h8bx/65bDlI1xa9LkifEUYCRwJJaEMLkpUkB2kU0euCFIcj\niSPsitM0oTAVg5ScUVgvDF/yDJ6PkV4m6iAt0eiX9eYRGh/FTh+DEI2RknPck4kkn7UbRfSK\nIIHLuF9+2TB6j+NRdoHdtk/fjpxAopmUqiY7lSP1InpZkLAxHMGHYMO5lIPt4E17tYYEcidd\n0pFAo1yy/ibRjkAi03eQk4lyQ0FCSFC2BAMC2cBEaqOFNc4xIdrImqnCRnGQsikOBOnnevc1\ngugdvcMufhhd4EhegScYVlM3pEZazk0QbWWtVFGjOEjZFEeDhOYb6H1YhT2UmL3IDyZyVUg+\n20Ca8slZi05TNFFBNg7SIaJdgPS5cbQ9NtEUJOTADAQmkHDOBisAqSQfB+kQ0d5AsnJEh0zc\nU0pcSlySrMOG2tCc9WYHqSgjHyMdItoZSFaMtr7ePt/wHnQs2tgs564MSPn6UaOndQ6Sz9pZ\nRLsAieUIfverSYpdfk6DpPTTGKQoiUKJJOkdpCFFrwxSmGvgw1F61g7BMsUuP5NQIG/oGlI4\nZReXLxAdyOcdpFYgvYPlQPrJuvupIkgTBxIMBWCv4SufzNrFCdRKgqhgdo4G8nkHqRFIqXtH\nVpLSIAGrB4m1nUAqsHF8/lT1PxCkSeCImc02kJS7CHOQzi3qIME1QQqQZrpn+aC9CGPuMFU1\npJ2jczlSL6KXB4mjYtvBgcSOk/JVjnHMFnSOVeFpUzQOU9qpHKkX0QuCNAscMSFJJCm+r6MH\nyeL38aLV7fSCUMSLNrJxfP5U9T8SpJ88R9y1HU/Sw8i+G9zPZNoCJBLSSkk6lSP1InpJkD45\njNjJBn7f3ej2mi9kEuXqIJ1a9JIgqTkC13JkZ5TmtiV97ozdvIAjJUhVoq1sHJ8/Vf07AEnN\n0cPEtLjKssB6WNF8oJmi05kxkhXPUzlSL6LXBSlyeQakCbwwqAikaTKzg5uJbEOxqfQi71SO\n1IvoCCB9fBn9+8HsV4P0GXOUvB17P0lMS6oc8InT2E3TO21AquJdVK21cUQHAOlj+W/9e/sQ\n9lsiEgdLkiSJq6jKTUl6GUh1vEuq1TaO6BVBYlfX3Q/oQHqKTOSqLVoi9DKQWoyRKospqNbb\nOKIDgBSoAeBUgDQJD03Q1T/oWJYKZtHqq0BqMGvnIPWpui9I6xApAun/7paVCu6Nf66S/Hgl\nSIhOemwq9G9YoG8bpZxusqlAglGIAmSOSEIMEh82mmlCKoe/O/BZbWftmpiPkYZRbR6RPsjf\nGpC4+5oxSCQdJQnvn5kxEtrWNBrbTC1EZFFg1frj+PyFQYo4qgLpxgaaWRojLWfFx8DeGb2y\nOOKo2EnxXGChiCja1Mbx+VPV3wQS5qj20u5euikyug4IIybcrU0uCVoM42g0fHeqEUmncqRe\nREcAicSjjxs/VjKApCVlOy4SxigRQ0esMBSDlEp8KkfqRXQAkD6WaboPaUWDdWUDA1IKIwtg\nTG7wiJGGcpCSqU/lSL2IDgCS0TSlS/PA3VHFSWe4ARNz2YUjChzI4cIxUjqjUzlSL6LXBElc\nObc48sRdtqGoglPMKU8naZI80OOFs3YO0stFLwqScAOJppQ4InPiD1HRcQ0gRQkKe6cKJBOy\natUyG0f0qiAJjx9hJ2I4Ij4aPsv5UpASJWsFUs0YKVfGMtVCG0fUQQIhCV22CasZyPb2MZGv\nmqN2IJXP2mWjZpFqqY0jelmQECVccJJWBQl+hvIlKSbdhR0j/vrecZB6Uh0LpHVLnoJARndM\ndImQwJrKRUkCWpsiF6fmIA2jOhZI2zydBBK/8nQCc9rxWjtx5iJfPlpUImESUIgSK81jHJ93\nkFqCxK5JIOwwKUJQwLgpQGJdNOO0uDbF0SIlGllhDuP4vIPUFKQcSbGF5XghJQMSAQ8Yy1HS\nbY8AqSfVcUQvDJKdJMyUAJLMUWzZpA7SKKIOEgIpTQ13LBojWbzdCNJLxkhdqY4j6iDFVKj2\nIvjmF4H0glm7vlTHEXWQsrFnlo+xVbaEjVxad6RRRK8Mkm5AFMJGiiMya6cOHZl07kijiF4a\nJA1JcPaAQUyosiUqrSeki9rQTuVIvYg6SDmOwNrwmKPg/nSJkI0kKbk70iiiDpLM0OrbAmGY\nliqQxPTuSKOIXhskNsZEXs0ed5BeoTqO6MVBwiRJC1WZ43MSJOsYyUEaXvTiINHLNX5hAsvR\nnBgjmW/4SOC5I40iemmQJmIzz1EEEpeush0F8NyRRhG9MkiUI3y9Bl07SoejmK7KJeaONIro\nhUGKOGLvGUlpuwBJiGN1ooU2js+fqv6dgxST8nqQJvyUE5+im2dZB/L5U9X/eJDkIMOiIrCW\nrXJZ0FgzTNaGFkFrp3KkXkQvC5KaoyRI2SpjhXzZ8InpMxykjkTnW5ul+VQ0m2I0kOijsTqQ\nQDKb08cZJMVN5iDtolrWGRnRfIphQNpOWz8nD0JDyYxerwEpxWbq1LiJWnjAOD6/D0ilX2tJ\nGwAk1crvNB50flxOOJeCVDprl8wrEm3iAeP4vIPUFCTuoYiIpC1l+MBBxjUikUm3s4BhaZ+n\nM6OibVxgHJ93kJqCFEPD0BRS0lM2NW4f3j9jFca4Y489idqkOq0IpEovGMfnfYy0L0grTBQX\ncJyBhoULngZ3CAWS3ViuTbLbSkGq8YNxfN5n7XYFaXEkESQGsbkNRyUgZRw/eZAfI9WSNI7P\nn2rW8niQNCTJ6RiQlk20f57ZTWKyF5eClJ21w8cdpE5VRwCJ54PjIMERcUH8gUkllEg8XAxS\nym5cxHKQelTtHiSBIh4kZhZinmm6G9qKzs/4qRmkminrG5tdJUcD+byD1AwkmSN8NQfPookA\nXWuVXwlSxdCWB6l2rDyOzztIrUDSciSRhLemrcr0MH8ma8LRSufmTQCpWnUHG0f04iDRy7bH\n8QiszcLOKE1Yd7Qljc5Lea7E0Q4ksWOkFqrtbRzRi4PEMiGDlLg3e5tT7lkWAfYIHDM3a9dK\ntbmNI3pJkOQx0v0wt28zBrZFVGKvwnYEqb2N4/Onqv+xIIkkzdwhcGrYfn4KRyXyasxBGkf0\nmiCJEUnE62lwmyFsn5DUTm6xUzlSL6IOUilI2Mwg6dLtwdG5HKkXUQfJQpIKJFULKVOW9U5G\n+lSO1IvoNUHiBkJ6khjJIo4UaYt6Jyd9KkfqRfSiIEWrfubUsqEoJkWSJopmJUhT9glZklqp\nfSpH2knU0JcGVbP1D9LMPHqkAUn4vt8DJJMgnU90kKpEbW2vVbVb9yCxrKTIep7JgYUOqFtI\nyZFaEaR2kKpFrb2pUy2w3kESAw9FR9gKy1ZnJpnKsomLQfIxUrWog6QDSckR3cEdfioWgJS1\ncpB81q5W1EFqCVK0h8XsobgHSMVjpKw5SFkr6EwHqcDuKgidHTgqnrXLmoOUN3tnXhEk8w3Z\nBRNEy/4gzSN550BFPVX9jwUpJkk3/4Bg2XaFjebmjjSKqINktk1n42/daG/uSKOIXhMkEzQT\njVdAKuy57cKSO9Ioog6SKvgIq7vBntsuV3fuSKOIOkg8QZQZBiS850bPaGLDO1Jliwxf/91F\njwUpS1IKJH5PCqRybxrdkWq/W0av//6iB4PErPZJgpTlKAVShTcN7kjVUXrw+r9A9GiQEqQ8\nd/NHgwglSx4j1XjT4I7kIO0uejhI1O8JR/h6LA5IEVvyrJ2D1Fa02hykViBN4qLTCe9cz4gP\nRk4i53thkHyMtLvokSA9excHH0qSFH/CPuIjiXwrvGl4R6rjaPz67y56IEgSNHAfSw2NVxN8\nr10q33JvckcaRdRBSi+zY05it9u340PZHWkUUQeJowdRMoUR1XMPhU9ZZZs9lQXV8hB3t1M5\nUi+iDhIXhQhH0wTnH5j0qMp1Xo4LKQ48upsKG8jnT1X/40CKo05MBZhrwPtZDHGVK70cZyNP\nhVXkcSpH6kX0uiAJJIVUKHm8GWBEVa71cizjII0i6iARjgIf3OsYKHDhc2OQUmMkPgt9pqdy\npF5ErweS4aG+xGQd1prbg5SatRM40uZ6KkfqRfSCIJU9HotPJVKwyo04WppJrAFfpxrRShvH\n509V/yNBKnr5yXIucVcGpKYPyqp7x0E6WNRBspMEpDiQWpqDNIroJUEq4QhezQWpbkDyMdLB\nolcEqYgjPL2w8bQdOrx31BwdX9Qzil4QJB02SY7AvvXQqXqnG9VxRB2kDEip55ZgCDhV73Sj\nOo7o+UDKmgqkkCp15quL7uYmWp9jJO49DVyqsJC1/Gsn9X0zjOhART1V/Y8E6YZB4UHCtDxM\n4oima2juSKOIXhKk5x8BH3kglOZtF5LckUYRvShIfEx6HmYDVDiSeg6josnEora3UzlSL6JX\nBEkiAR9mjsHzEskamjvSKKIXBIklIIaDpAMnCvFsD5LckUYRvR5IMkdJk55Mmplj7cwdaRRR\nB6kMpAgcB+naog6SlqPo1UKYm9seHLkjDSN6PZDKn6GAmzFIu9yUdUcaRfR6IJVyRLZpADpV\n73SjOo6og1QBlq3KJVaiOuUiYz9FPZGog1TIkb3KJVagKpawRlRh4/j8qep/HEj1JOWqnHNk\ng9l7J1XGYlGNjePzp6r/gSDVkrRIiFXOOrLBHKRRRK8IkunlJzMFb55ZVrZ8FZ6sNwdpFNFL\ngiRTk96LXtRAnJUuKV/31zHlY6RRRB0kJUlQQAYJnAOS6hqNLWrBOdkMT+VIvYg6SAqQKEki\nSDxHFSS5I40iekmQ5DFSavQUBDg6MEg4pYN0flEHSW/rVRMLBwQp7HWQriJ6SZDKOCLhiq1y\ndKSOI3ekYUSvCJIm8CiwYqscHajiyB1pGFEHqRQroniq3ulGdRxRB6kNRx2+104s06kcqRfR\nK4JkGSNp3s3A7MqaMn1F78hlOpUj9SJ6SZD0MWnm32aHTQfSFM3m5Uta0TuJQp3KkXoRdZAy\nHHFv7yJ6aK8ECDpRh97sII0j6iDlOAouyZEUDoDETJ74RAfpdKIOUhakiBP67iAwaye678tB\n8jHSa0WvCZLxThFPAE5IQSJeTCSVHPms3TCiFwVJJAke2E5gCSAJCUgRKWSHjiN3pGFErwqS\nIiZt6Wksmen5oMqQo4ikfNHYora2UzlSL6JXBSnPkXQ5h/duR25h/2wYBCmK2txO5Ui9iF4U\nJA1HJPRsn2e0e93G+TpIVxO9Jkg6jtBz5VGcSYLUatWQO9IoopcEKc9PDBL5ONNREs23CUfu\nSMOIXhEkRSDioxAGCY+TTtU73aiOI3pBkLIcbfzMbBDilU7VO92ojiPqIMWRaE01szFIkDpV\n73SjOo7o9UBKBiKJFE4Fi52qd7pRHUfUQaIc8ZucCPp8qt7pRnUcUQcJskK2Y9CwRtg6We90\nozqO6PVAkpfZ0SPRNj3dVuUSc0caRfSCIMkPj2dAohd/DtLuquOIXhEkdonQ/QBPmMiascol\n5o40iuglQQLRJkwbcCEoCZK1yiXmjjSK6BVBwgRlAQqnzAAke5VLzB1pFNELggRiEd6TCjwT\nmWsoqHKJuSONInphkOiOOASxxPB7z9U73aiOI+ogxdd1617TsxCn6p1uVMcRvSBIiachlv3h\n2k9P0ql6pxvVcUSvDBLCBYDEBagGVS4xd6RRRC8MktbwmdIje6fqnW5UxxG9Iki2HxoL5yWj\n1Kl6pxvVcUQvCZIpJm2nSWFKW2UopE3qjjSKqIOkJUmKU1BU1ThCSBOKuoOdypF6Eb0kSCaO\nVr/PgaTlQxLgi6pLZrNTOVIvog5SYANQwpCU5kjPhyrhmsAdaRRRB4njiFmBl5lraAvSlsId\naRTRS4LEjZHY1zOwnPEz4GqQFGOkoOWONIroJUFKxiO0SpVJKSwdUo+RFLN2DtJ4olcEiUGD\nRwauXWVPQlXWcqQu4OyONI7oBUFiiUnFnmniph/wndq27bjpuyONIuogiaAAkNjTsGDbdlzl\n3ZFGEXWQUkZ/yIU8XIv0ihsr2UzDiA5U1FPV/ziQ0uBEt4sIInE8cpB2VB1H1EFiIhBDCzmf\nEVQ3jsXckUYRdZAyEUo4Pdp1qt7pRnUc0euBpCFpfUZWPpXuPlfvdKM6jugFQVKTlDoxumt0\nqt7pRnUc0SuCFC0Ryt5tnejwKR4pnap3ulEdR9RBEkISND5g4VSn6p1uVMcRvSRIeY7QxDd7\nyEF6heo4olcESc3RxK1pmNmbSefqnW5UxxG9IEgWjth7tOwTF6fqnW5UxxF1kJIYcVvRvaTH\nx1P1Tjeq44g6SBxITAhCi4UoRzQ8NTN3pFFELwiS8QWRkkSkpm8evbkjjSLqIJVwxKtpWsZo\n7kijiDpIiBqWozQjDtKOquOIOkgxNWRfDhIHaT/VcUQdpPCwHnOQXcQg6CnbxmTuSKOIXhAk\nORqFwyGZ7YWOrc0daRTREUD6+DLN3yqQUEwK6bifnC2scom5I40iOgBIH8t/ub+VIMVxZwKj\npuxI6VS9043qOKIXBMnwsrotff4C71S9043qOKIDgLTS1Ayk23alpgPpbg7SIarjiJ4KpP+7\nm0putQcdt+Xf43MqoXzYze1oU4H0cWsZkcDnEGUyEcfHSEeojiM6SETaDSRgmWk5n7U7QHUc\n0TFA+oD/7QVSnZ2qd7pRHUd0CJA+EE0OUo+iAxX1VPW3gPSBw5KD1KPoQEU9Vf0NIH18LEsX\nGq1s8N7ZQ3Sgop6q/paIZLMWpSuwU/VON6rjiDpIjexUvdON6jiiDlIjO1XvdKM6jqiD1MhO\n1TvdqI4j6iA1slP1Tjeq44g6SI3sVL3Tjeo4og5SIztV73SjOo6og9TITtU73aiOI+ogNbJT\n9U43quOIOkiN7FS9043qOKIOUiM7Ve90ozqOqIPUyE7VO92ojiPqIDWyU/VON6rjiDpIjexU\nvdON6jiiDlIjO1XvdKM6jqiD1MhO1TvdqI4j6iA1slP1Tjeq44g6SI3sVL3Tjeo4og5SIztV\n73SjOo6og9TITtU73aiOI+rtiZsOAAAD0UlEQVQgNbJT9U43quOIOkiN7FS9043qOKIOUiM7\nVe90ozqOqIPUyE7VO92ojiPqIDWyU/VON6rjiDpIjexUvdON6jiiDlIjO1XvdKM6jqiD1MhO\n1TvdqI4jej6Q8mb7dcxDzYu6g41TUlNRHaSEeVF3sHFK6iC1Mi/qDjZOSR2kVuZF3cHGKWnn\nILm5ndAcJDe3BuYgubk1MAfJza2BOUhubg3MQXJza2CvBunjy16cpdbWon0sH+A2/NuB5YrY\nTVE/YEm7btRnKXINKhf3xSB9bP91Z1vRPvA2/duDZYrYU1FvYzTqRyhJokETxXWQVnOQdjJY\nml5L+gFK4iA1sI9Qul77/JYvYkdFvaEvp45L6iA1tY9buJpft8nfDixXxI6Kuham+0Z1kFra\nIN5Z0d+vN5WHHm8OUkv7AB/67fOnDQXS+qnfkjpIDW2QPn/aICB9oI/9ltRBamcf4f+u+3yo\nS7tRGtVBamYf4c/6r88+zxaxo6ICkPpu1MFA6udOdmQfuZvvHRW9/Ab8y239euq9pAvwo6xs\ncHM7pzlIbm4NzEFyc2tgDpKbWwNzkNzcGpiD5ObWwBwkN7cG5iC5uTUwB8nNrYE5SMPaW2nf\n/bePtQTnMgdpWCsGqfhEN9m8TYc1B6kn8zbt1f79/e3t939vD7//7e232z+/vv3277L145/b\nwsM/91T/3P59+/VxztcfnBzK/PPb28cf9w9OUnvzJu3VPu4Of8fj7e23r09//vr13+/3rS8y\n3j6eSN3+faT62vrj7a+vpH++/YckhzKPj384SLuYN2mn9p8vl//C478PdL4I+dr88w7A29uP\nf28/3v54gvTH24/bY+t/9w9foehvkhzKfJ3437cPv7TbxbxNO7VfHz3zdY12vya7/7fEoLe3\n/31d0N1jzH3r1/uxx9Zv9/1PTGByIrNIHFqzc5q3aaf2ttji99t/TwrWT2Hrf1+w/PW89sPJ\nkYyDtJd5m3ZqVpDuwekxUHKQDjFv007t161nKBmPi7kf9NLuKxz98fEWJScyDtJe5m3aqf1x\nnyX4cwUGgvTj9u+Pt/+QyYbbnZnHhANOTmQcpL3M27RTe05sP2YQKEj3Ce8bmf6+3UPS2583\nmpzILBK+Rqi5OUi92v1W64+/b8yl3Y/7LVh8Q/ZhYV4O/IdlHv//10Fqbw7SYCZfl/39XN3g\ndog5SIOZDNKPx+IGt2PMQRrMJJDenlMNbgeZgzSYSSB93FcvuB1mDpKbWwNzkNzcGpiD5ObW\nwBwkN7cG5iC5uTUwB8nNrYE5SG5uDcxBcnNrYP8P1F3jEefAAeIAAAAASUVORK5CYII=", "text/plain": [ "plot without title" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ggplot(data = df, aes(x = employment,y = wage)) +\n", " ylim(20000, 30000) + xlim(0,10^5) +\n", " geom_point(color='cyan') +\n", " geom_smooth(method = \"lm\", se = FALSE) +\n", " geom_smooth(method = \"loess\", color=\"green\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Now let's perform the linear regression using the \"linear model\" or `lm` command to see the estimated coefficients:" ] }, { "cell_type": "code", "execution_count": 122, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "\n", "Call:\n", "lm(formula = wage ~ employment, data = df)\n", "\n", "Coefficients:\n", "(Intercept) employment \n", " 2.345e+04 2.311e-02 \n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "lm(wage~employment,data=df)" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "Call:\n", "lm(formula = wage ~ employment, data = df)\n", "\n", "Residuals:\n", " Min 1Q Median 3Q Max \n", "-66231 -3382 -681 2471 47744 \n", "\n", "Coefficients:\n", " Estimate Std. Error t value Pr(>|t|) \n", "(Intercept) 2.345e+04 1.041e+02 225.33 <2e-16 ***\n", "employment 2.311e-02 8.014e-04 28.84 <2e-16 ***\n", "---\n", "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n", "\n", "Residual standard error: 5701 on 3214 degrees of freedom\n", "Multiple R-squared: 0.2056,\tAdjusted R-squared: 0.2053 \n", "F-statistic: 831.7 on 1 and 3214 DF, p-value: < 2.2e-16\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#It's usually useful to save the estimated model as an object:\n", "model1<-lm(wage~employment,data=df)\n", "summary(model1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Important! The standard errors that `lm` computes by default assume homoskedasticity, which is not reasonable in most social science contexts. Based on the above scatterplot, would it be reasonable to assume homoskedastic errors for this regression? Doesn't look like it.\n", "\n", "## To get heteroskedasticity-robust standard errors, you can use the `sandwich` library: " ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\t
(Intercept)
\n", "\t\t
171.42093890755
\n", "\t
employment
\n", "\t\t
0.00506838094830457
\n", "
\n" ], "text/latex": [ "\\begin{description*}\n", "\\item[(Intercept)] 171.42093890755\n", "\\item[employment] 0.00506838094830457\n", "\\end{description*}\n" ], "text/markdown": [ "(Intercept)\n", ": 171.42093890755employment\n", ": 0.00506838094830457\n", "\n" ], "text/plain": [ " (Intercept) employment \n", "1.714209e+02 5.068381e-03 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#install.packages(\"sandwich\")\n", "library(sandwich)\n", "cov <- vcovHC(model1, type = \"HC\")\n", "model1.robustses <- sqrt(diag(cov))\n", "model1.robustses" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
\n", "\t
(Intercept)
\n", "\t\t
23452.3821454205
\n", "\t
employment
\n", "\t\t
0.0231118203981946
\n", "
\n" ], "text/latex": [ "\\begin{description*}\n", "\\item[(Intercept)] 23452.3821454205\n", "\\item[employment] 0.0231118203981946\n", "\\end{description*}\n" ], "text/markdown": [ "(Intercept)\n", ": 23452.3821454205employment\n", ": 0.0231118203981946\n", "\n" ], "text/plain": [ " (Intercept) employment \n", "2.345238e+04 2.311182e-02 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "model1$coefficients" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "Call:\n", "lm(formula = wage ~ employment + 0, data = df)\n", "\n", "Coefficients:\n", "employment \n", " 0.06989 \n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#If we wanted to supress the constant:\n", "lm(wage~employment+0,data=df)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAMFBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAgAElEQVR4nO3diXqiyhZA4UIRbQd8/7dtAUVwJOWuYe9a/3fviUkHC6la7UQn\n7gzgZy71DgAWEBIggJAAAYQECCAkQAAhAQIICRBASIAAQgIEEBIggJAAAYQECCAkQAAhAQII\nCRBASIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBASIAAQgIEEBIggJAAAYQECCAkQAAhAQII\nCRBASIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBASIAAQiqZY/qlcCRHblxW10vu5TrbRNmZ\nbeWc30iv9nrytfp+sRtjtfO83tcH5+uemFXQTf1mUUiHKsoR217GDhLS/n5x7XqN3/US0lxB\nN/WbRSFFWhwr546em37c666j68Wdu1o4ECF9VNBN/eYppM/fFGlfRDa9fu1Uu3tIl1b/9f+T\nu+P7+7fYUdBN/ebdPVK77R4E1f+Gr4wrcb/pHn3tr5ucLp+td5MtT6v+QdO/bu2umtPt+nYr\ntzpc7g4qtz7Mh59d34u//l9s+bgLVfeU57Zp21SuGgae7NZ6vObuQv+/arbD080mt/y+S/Nh\nbl+9fXy8vY/XYhYhjd6EdKqu+axnIV2fYbi63+Jw/Zb7lqt+g9t3ucN5XMeXNduMXxtNr2/a\n622HXmz5cRduu32Y3aD6fiurIaRxgGGHp5tNb/l4Hc+3dHrIHm/v07WYRUijNyFd/ta//GXa\nXpbIbrLC69uKGZZxNX5629J1m12eh6zb87mZ9XG5D5gGMJhd36uQnrf8vAu3T6vJDVrvJ7fy\ncrs201vc7/Bss+ktHzd8vqWTQ/Z0e5+uxSxCGrmp6xeG/3YPVNrL39nj1/on7bv28qjl8vGy\nOi/PNKruQ3XfsltQ3fOQ0+yaLl/tnuSvjv2H+9gP1/fiod3Tlq93YV9Nl3Q7rOLp1Y0XT0OD\n+8kAD5u9uuUPwzyE9HR7n46fWYQ0ehNSt2TG5yG3hbO5/f3a9E/W62H5d8vstuX+4aqH/x5m\nH+7f8HB9L0J62vLNLlxf37582g5b1ufXIXWv5Hduj9v6rWebvbrlD8M8PUd6uL1Px88sQhq9\nCWk7fGEz/tV9/dAvuP7v9eH5xss/vnzDv2btpgvr5YvrD9f35rXmx3vK6SZuvgv3W1Kd34Q0\ndD+8keTGgCabvbrl7vmWzq/26fbOr8UsQhrdl9h8kTS3pXV6+qPbpXfL6/xvNSnzc0izS8tC\n+rAL878T3oR0+eQ4e8T6sNmLW/4lpFe3d3YtZhHS6F1I5/bf8GLUevZH491B9fIeqf+0+yt/\ntdkdl4Q0vb6/3iNVTyu8ml3Bi5CqbnPnHvqYb/biln8O6eXtnV2LWYQ0ehtSp3/L5v61+utz\npP5PV9evfw2p/voc6enD512oZ8/SXoR0uT3b55Dqhyd3T7f8xS3tcj4Mn768vbNrMcv0jfub\nNyGtxqcPt7uK9u2rdu71A57v90jfX7V7+vCwyW7+wmG3R4f+w/r8MqRu88v/T/NvmG32fMuf\nhunuwprr6xZvbu/8WswipNGbkC6rbX3qnzN3z8q7FdN9HN95HM6wefPuyrr/5v18mb0K6fH6\nFoT0ZRfGTx9eIRwvjm9DzV4fn242v+WP1zt8upl9+ur2zq/FLEIavQlpfLLcP8Tf3C7clvFw\nptr++i0PkVxPA+jeRz18Dunh+paE9LDJdaz6+qfXPRpW76uQbpuvZ1+dbTa75ddveRjmdP32\n4dOXt3d2LWYR0uhdSMPj+/X1ffl6bGdTTV7UPXbn2u2fIum+XG2Op9sJCy+u/WZ2fYtCetiF\n+el+57ZZTd9vfb6V3dl7t3+PdP/qdLPZLb99y8Mw/Q3/Nz6oe3F758fPKkKS1Cp7ImD7+X9U\nHEkJbnimcVzPT6DLHiGJ4UhKuD/vf3r1GGUgJAnjPxWw/coU3iMkEe22exmrsn5CGd4iJEAA\nIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAAIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAA\nIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAAIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAA\nIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAAIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAA\nIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAAIWXPvZV6z3DHZGTv7RQxdxlhMrJHSBowGdkj\nJA2YjOwRkgZMRvYISQMmI3uEpAGTkT1C0oDJyB4hacBkZI+QNGAyskdIGjAZ2SMkDZiM7BGS\nBkxG9ghJAyYje4SkAZORiff/WIKQNGAyMvF+IghJAyYjE4SkG5ORCULSjcnIBCHpxmRkgpB0\nYzIyQUi6MRmZICTdmIxMEJJuTEYmCEk3JiMThKQbk5EJQtKNycgEIenGZGSCkHRjMjJBSLox\nGZkgJN2YjEwQkm5MRiYISTcmIxOEpBuTkQlC0o3JyAQh6cZkZIKQdGMyMkFIujEZmSAk3ZiM\nTBCSbkxGJghJNyYjE4SkG5ORCULSjcnIBCHpxmRkgpB0YzIyQUi6MRmZICTdmIxMEJJuTEYm\nCEk3JiMThKQbk5EJQtKNycgEIenGZGSCkHRjMjJBSLoxGZkgJN2YjEwQkm5MRiYISTcmIxOE\npBuTkQlC0o3JyAQh6eY/GYdt7Tp1cxDcn2IRkm6+k9Gu3N1adJfK5BPSe0F2ER/4HvLGVf+O\n/aXTvnKN3A6VSvQeiZCi8z3klTuOl4+uktmZkhGSbr6HfPbogYcSvyMk3bhHygQh6fbDc6T9\nqb/EcyQRhKSb9yFfT14jWrWSu1QmQtLth/eRmv59pKre8j6SAELSjUOeCULSjUOeCULSjVOE\nMkFIunGKUCYISTdOEcoEIenGG7KZICTdAp0ixKnIf0VIukW4R2JWlyAk3SKcIsSsLkFIukU4\nRYhZXYKQdItwihCzugQh6RbhkDOrSxCSboSUCULS7cdDvls5V++DDlEIQtLtt/eRrq84fD6x\ngVldgpB0+ymkxjXt+Xxq3C7EEGUhJN1+Cqly/everVuFGKIshKTbTyHdzv75fBYQs7oEIen2\nU0ibW0icIvQzQtLNP6R6u9u7f5eLbcMpQr8jJN38QxrP7Hau4hShnxGSbt6H/Hjc7eq6f8mh\n+fzTuJjVJQhJN85syAQh6UZImSAk3QgpE4SkGyFlgpB0I6RMEJJuhJQJQtKNkDJBSLoRUiYI\nSTdCygQh6UZImSAk3QgpE4SkGyFlgpB0I6RMEJJuhJQJQtKNkDJBSLoRUiYISTdCygQh6UZI\nmSAk3QgpE4SkGyFlgpB0I6RMEJJuhJQJQtKNkDJBSLoRUiYISTdCygQh6UZImSAk3QgpE4Sk\nGyFlgpB0I6RMEJJuhJQJQtKNkDJBSLoRUiYISTdCygQh6UZImSAk3QgpE4SkGyFlgpB0I6RM\nEJJuhJQJQtKNkDJBSLoRUiYISTdCygQh6UZImSAk3QgpE4SkGyFlgpB0I6RMEJJuhJQJQtKN\nkDJBSLoRUiYISTdCygQh6UZImSAk3QgpE4SkGyFlgpB0I6RMEJJuhJQJQtKNkDJBSLoRUiYI\nSTdCygQh6UZImSAk3QgpE4SkGyFlgpB0I6RMEJJuhJQJQtKNkDJBSLoRUiYISTdCygQh6UZI\nmSAk3QgpE4SkGyFlgpB08z/kh23tOnVzCDVESQhJN99D3q7c3TrIEGUhJN18D3njqn/H/tJp\nX7kmxBBlISTdfA955Y7j5aOrQgxRFkLSzfeQO/fuE7EhykJIunGPlAlC0u2H50j7U3+J50gi\nCEk370O+nrxqt2qDDFEUQtLth/eRmv59pKre8j6SAELSjTMbMkFIuhFSJghJN04RygQh6cYp\nQpkgJN04RSgThKQbb8hmgpB0C3SKkJvyHKIshKQb90iZICTdOEUoE4SkG6cIZYKQdOMUoUwQ\nkm6c2ZAJQtKNkDJBSLr9fsi/vrzNrC5BSLoRUiYISTf/N2QXv+fKrC5BSLr5HvJDRUiiCEk3\n70Pe1m7dvyPLQzsRhKTbD4f8n3P/zoQkhJB0++WQn9aubglJBiHp9tsh37pqT0giCEm3Hw/5\ncfX9n0kwq0sQkm4/H/INIYkgJN04RSgThKQbIWWCkHQjpEwQkm6ElAlC0o2QMkFIuhFSJghJ\nN0LKBCHpRkiZICTdCCkThKQbIWWCkHQjpEwQkm6ElAlC0o2Q4nJvvd/kz3/AIY+PkOKSrIKQ\nMkJIcRGSUYQUFyEZRUhxEZJRhBQXIRlFSHERklGEFBchGUVIcRGSUYQUFyEZRUhxEZJRhBQX\nIRlFSHERklGEFBchGUVIcRGSUYQUFyEZRUhxEZJRhBQXIRlFSHERklGEFBchGUVIcRGSUYQU\nFyEZRUhxRQrJ44cV4SeEFFfqeyQmIxBCiouQjCKkuAjJKEKKi5CMIqS4CMkoQoqLkIwipLgI\nyShCiouQjCKkuAjJKEKKi5CMIqS4CMkoQoqLkIwipLgIyajpgV1tT6GHKB4hGTU9sM65EC0x\ndxOEZNT0wLb/NiFaYu4mCMmoxwN72K6kW2LuJgjJqBcH9lhd7pd2QYcoFyEZ9Xxg9+v+nySv\nAw5RMEIy6uHAttvL3dFq315qqgMNUTZCMmp2YA/diw3NcfgDsUPO3E0QklGz95Eud0a79vYH\nVYghikdIRs3eR6r3oYcoHiEZNXsfKfwQxSMko2YHtm26x3NVI1sUczdBSEZND+yp6l9hcK4S\nPbeBuZsgJKOmB3btNt19UdvIvfT9OETxCMmo+UmrjxfEhygeIRk1PbCVG54ctYQUDCEZNT2w\njVsfLh8Oa9eEGqJ4hGTU7MCur7/6Q+48u6chSkdIRs0P7L+6y0jwzO/nIQpHSEZFOLDM3QQh\nGeV/YA/bun8cWDeHUEMYREhG+R7YdjX5dYqfn1MxdxOEZNTswG5Xi3/TaOOqf8O/tzjtq8+v\n8jF3E4Rk1PTAbv/wK3srdxwvHz//kwvmboKQjJq/Ibv89bpZa5/DY+4mCMmo5T3McY/kh5CM\nmh7Y2i3/9xOX50j74RxxniP9BSEZNT2wp2r95ZXsifXkVbvVxwCZuwlCMmr+0G75iw3n86Hp\n30eq6i3vIy1HSEb5h+Q1RPEIyShOEYqLkIziFKG4CMmo+YHd192junrBj2zgFCE/hGTU7MCu\nh6dHS374CacI+SEko6YHdufW/b8y37nN1+14Q9YPIRk1P0WoPV9/INf37T6fEuFcoJcAtSMk\nox57WBoS90h+CMmo6YFdXe+Rjm71dTtOEfJDSEa9eI60X3QWOKcIeSEko2YHtv7LTxHiFCEf\nhGTU8/tIrv4XcojCEZJRnCIUFyEZ9euB3VVu9eUZFXM3QUhGeR/YY+2q3fXHPHCK0GKEZNT8\nfaTl76Ee+29rul8Ec6o/v8rH3E0QklG+IW26946a4Z3Y9vP7TszdBCEZ9eLAHtYLfs/Y0Nr1\nN5LxU4QWIySjXh3YdsFJq0M7/4bHdJwitBghGfXywC56aLe5nc7QbjhFaDFCMurVgd19vofp\ntdXkF2V+/nbmboKQjHr9YsN2wZbNLZ/qy+/3Y+4mCMmoVyF9e4f1hyGKR0hGcYpQXIRkFCHF\nRUhGvXlDVvJfiDN3E4RkFCHFRUhGzQ7sttpf/nuoFv3DPr8hSkdIRk0P7Pb6A02ObsE5Qn5D\nFI+QjJo/tHu8ID5E8QjJqOmBrcZ7pO8/RchziOIRklHTA9v9iK3Lh2U/RchviOIRklGzA3v7\nEVtfzvn5ZYjSEZJR8wP7r/8pQvuQQxSOkIzizIa4CMkoQoqLkIyaH9jlv2jMe4jCEZJRzy82\nnBf9ojHfIUpHSEZND+xfftGY5xDFIySj5m/ILv9FY55DFI+QjHo8RYiQwiIko6YH9i+/aMxz\niOIRklEvniNxilBAhGTU7MD+6ReN+Q1ROkIy6vl9JH7RWEiEZBRnNsRFSEZND2wte9b3qyGK\nR0hGPb78HXiI4hGSUY8vfwceoniEZNT0wLb1+hB4iOIRklHzh3biP9PucYjiEZJRhBQXIRnF\ny99xEZJRhBQXIRl1O7CBXvqeDoEzIZk1DylITszdBCEZRUhxEZJRhBQXIRlFSHERklGEFBch\nGUVIcRGSUfeQgvzay+kQOBOSWYQUFyEZxZkNcRGSUYQUFyEZRUhxEZJRhBQXIRlFSHERklGE\nFBchGUVIcRGSUYQUFyEZRUhxEZJRhBQXIRlFSHERklGEFBchGUVIcRGSUYQUFyEZRUhxEZJR\nhBQXIRlFSHERklGEFBchGUVIcRGSUYQUFyEZRUhxEZJRhBQXIRlFSHERklGEFBchGUVIcRGS\nUYQUFyEZRUhxEZJRhBQXIRnlf2AP27r/eft1cwg1hEGEZJTvgW1Xk99dsQ4yhEmEZJTvgW1c\n9e/YXzrtK9eEGMKk5CG99XYTLOF7/Cp3HC8fXRViCJOSh/T3TbCE7/Gb/Q32+a8zpmiCkIzi\nHikuQjLqh+dI+1N/iedIf0FIRnkfv/XkeeqqDTKERYRk1A/vIzX9+0hVveV9pOUIySjObIiL\nkIwipLgIyShOEYqLkIziFKG4CMkoThGKi5CM4g3ZuAjJqECnCHE65BuEZBT3SHERklGcIhQX\nIRnFKUJxEZJRnCIUFyEZxZkNcRGSUYQUFyEZ5X382o1z6/31Sj5eC1M0QUhGeZ8iVA0n2g1X\nQkhLEZJR/i9/7y417ar+NDtCWoyQjPJ/Q7b/cKpWJ0L6A0Iy6tdThNr1mpD+gJCM8j1+K3d7\nE3a1JqTlCMko3+O3c5vrpZNbE9JihGSU9/Frxnr2X07wZoomCMko/+N3rG+XThtCWoqQjOLM\nhrgIyShCiouQjCKkuAjJKEKKi5CMIqS4CMkoQoqLkIwipBDe/4JJQjKKkEKIUwUhZYSQQiCk\n4hBSCIRUHEIKgZCKQ0ghEFJxCCkEQioOIYVASMUhpBAIqTiEFAIhFYeQQiCk4hBSCIRUHEIK\ngZCKQ0ghEFJxCCkEQioOIYVASMUhpBAIqTiEFAIhFYeQQiCk4hBSCIRUHEIKgZCKQ0ghEFJx\nCCkEQioOIYVASMUhpBAIqTiEFAIhFYeQQiCk4hBSCIRUHEIKgZCKQ0ghEFJxCCkEQioOIYVA\nSMUhpBAIqTiEFAIhFYeQQiCk4hBSCIRUHEIKgZCKQ0ghEFJxCCkEQioOIYVASMUhpBAIqTiE\nFAIhFYeQQiCk4hBSCIRUHEIKgZCKQ0ghEFJxCCkEQioOIYVASMUhpBAIqTiEFAIhFYeQQiCk\n4hBSCIRUHEIKgZCKQ0ghEFJxCCkEQioOIYVASMUhpBAIqTiEFAIhFYeQQiCk4hBSCIRUHEIK\ngZCKQ0ghEFJx/I/fYVu7Tt0cQg2hFiEVx/f4tSt3tw4yhGKEVBzf49e46t+xv3TaV64JMYRi\nGkN67+2V4c73KFXuOF4+uirEEIppDMnjynDne5Rmf099/kurwIkgpOJwjxQCIRXnh+dI+1N/\niedIzwipON5HaT15NrpqgwyhFyEV54f3kZr+faSq3vI+0iNCKg5nNoRASMUhpBAIqTicIhQC\nIRWHU4RCIKTicIpQCIRUHN6QDYGQihPoFKHCT3okpOJwjxQCIRWHU4RCIKTicIpQCIRUHE4R\nCoGQisOZDSEQUnEIKQRCKg4hhUBIxSGkEAipOIQUAiEVx//MhsUnLxQ4EYRUHN+jtCOkDwip\nON5H6Vh9/scTAkPoRUjF8T9Kx88nBkkMoRYhFeeHo7SbnLcaaAitCKk4vGoXAiEVh5BCIKTi\nEFIIhFQcQgqBkIpDSCEQUnEIKQRCKg4hhUBIxSGkEGyFxO/EXICQQrAV0t83KRAhhZB6Iace\nv0CEFELqhZx6/AIRUgipF3Lq8QtESCGkXsipxy8QIYWQeiGnHr9AhBRC6oWcevwCEVIIqRdy\n6vELREghpF7IqccvECGFkHohpx6/QITk7/25M6kXcurxC0RI/vJdyKnHLxAh+ct3Iacev0CE\n5C/fhZx6/AIRkr98F3Lq8QtESP7yXcipxy8QIfnLdyGnHr9AhOQv34WcevwCEZK/fBdy6vEL\nREj+8l3IqccvECH5y3chpx6/QITkL9+FnHr8AhGSv3wXcurxC0RI/vJdyKnHLxAh+ct3Iace\nv0CE5C/fhZx6/AIRkr98F3Kk8flhxneE5C/5QlY4vlmE5E/jQk49vlmE5E/jQk49vlmE5E/j\nQk49vlmE5E/jQk49vlmE5E/jQk49vlmE5E/jQk49vlmE5E/jQk49vlmE5E/jQk49vlmE5E/j\nQk49vlmE5E/jQk49vlmE5E/jQk49vlmE5E/jQk49vlmE5E/jQk49vlmE5E/jQk49vlmE5E/j\nQk49vlmE5E/jQk49vlmE5E/jQk49vlmE5E/jQk49vlmE5E/jQk49vlmE5E/jQk49vlmE5E/j\nQk49vlmE9JXHD53KdyFHGr+8H9RFSF9pXMj5jm+1MEL6ytZC1je+DoT0VeqFVPr4OhDSV6kX\nUunj60BIX6VeSKWPrwMhfZV6IZU+vg6E9FXqhVT6+DoQ0lepF1Lp4+tASF+lXkilj68DIX2V\neiGVPr4OhPRV6oVUyPjKz3kgpK8KWcj6xs8KIX2V7UIqffysENJX2S6k0sfPCiF9le1CKn38\nrBDSV9kupNLHzwohfZXtQip9/KwQ0lfZLqRSxlfxujghfZV8ITH+XzdJgJCuPN4QzHchlT5+\nAoR0ZWohlT5+AoR0ZWohlT5+AoR0ZWohlT5+Av47c9jW/TOIujmEGiImUwup+PE9vL2yZXy3\nb1eTfVgHGSIIjwOpcSGVMX6cyVzGd/vGVf+O/aXTvnJNiCGCMLWQGD/GlS3ju33ljuPlo6tC\nDPGVzx24wolk/LTjL+O7/WxFPi9PrzUOZCJeSH+4RwLs++E50v7UX/r6HAmwz/uh4XpyR7hq\nJXcJ0OeH95Ga/n2kqt5+eR8JsC+r16YBrQgJEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJ\nEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJEKA5pEQ/qgm5SL0Ap7LamT9Kve+MX/b4M1nt\nzB+l3nfGL3v8max25o9S7zvjlz3+TFY780ep953xyx5/Jqud+aPU+874ZY8/k9XO/FHqfWf8\nssefyWpn/ij1vjN+2ePPZLUzf5R63xm/7PFnstqZP0q974xf9vgzWe3MH6Xed8Yve/yZrHbm\nj1LvO+OXPf5MVjvzR6n3nfHLHn8mq50BtCIkQAAhAQIICRBASIAAQgIEEBIggJAAAYQECCAk\nQAAhAQIICRBASIAAQgIEEBIggJAAAUpC2q1c1bT9xaZaflFyD25HKtH4j+KMMkh729PP/RI6\nQmr6Xz5QdUdo3V9cnZdcFHS8/e6DROM/ijPKIO1tTz/3i6gI6eg2bff34uZ8PrjqeD5W7rDg\nouQeVNfFlGj8R3FGGaS97ennfhkVIdXDXnbz2bj95dI/t11wUc7Ora+LKc34T+KM0kt825PP\n/UIqQrrqDmbtTufur6l6wUXBkZvzdTGlGf9JnFF6edz2dHO/dAcTjOmpdevzbU67D18vyjk+\nXm/k8Z/EGaWXxW1POPcLKQpp191xJzuY6RfT+70JLf1tTzv3S+gJ6VR199iE9GJvQkt+2xPP\n/RJqQmqrdfeBkF7sTWipb3vquV8i65Cmv7t6Pbw7UN2P1deLouNfP0Yd/704o9ykvu0p5v6v\nlIR0Wq1P/ZeGF2ZO99doPlwUHH+cnajjvxdnlJu0tz3N3P9V1iGN9m59vbTt3yrYu2bBRVHX\nxZRs/Lk4o9wkve0ZzP0SKkI6jccy3bvb18WUybvrcd/DT3nbc5j7JVSEtHFufJS16i/0B/fr\nRUm3h3ipxn8QZ5SrlLc9h7lfQkVIbnIw2/4E3/7LXy/K7sN52aCBxn8QZ5SrlLc9h7lftJ9J\nRgWMISRAACEBAggJEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJEEBIgABCAgQQEiCAkAAB\nhAQIICRAACEBAggJEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJEEBIgABCAgQQEiCAkAAB\nhKQPc5YhJiU3D7/c/vbp7vrxtOl+KV27/Ercuyl++wfwwMHMzeuQjvePnWr5lRBSFBzM3LwM\n6Vhdv7x2Tevay399r2/BH8ADBzM3r0LaufX1y90Hd26/3iW9vb4FfwAPHMy4Lqt366rt+dy4\n673KbuVWu+EPm+rytWF9X75a7a4bXP5z+/K5cu1tzpxrV66efu95f+ltvb9tdb++4dP+v/va\nXX/x9/DpbQv8hpDicm7bPcfp1u9QUn/Brc+3i3W/vuvxq/2nx/H+o3Gr/RhS3V/F/Xt3wxOo\n3fW779d3D2k7fEtz/fS+BX5DSHFdVnzbLd/+v5fHZ/9cdeyeAf27X7xMyb7788szof39Adjt\n4+ay7jeH8apm31u5Y3c1q+G7J9d3D8kNI7nrp/ct8BtCisu5Q//f03lYyXUXQBdDd/HQX+y/\n2iXSdg/cHkM6H5vufma8qvn37s/3755c3+Sh3f3Khj/gYZ0MQopruqTvS/z54tWLkC4X96vx\n4dvsey+J1cfj26se/nvab9djSPct8BtCikskpMv9z+pVSOdt1b3HdPoU0nr85vkW+A0hxbU0\npIcNHkJ6/P67fbO6RfYypI1b7fane0jjFvgNIcX1GNLtOVJ9u3i4f3WywfhxePm7fx/p+pX6\n8WnOrZzJ9fXfer80C+n82CJ8cAjjegxp8qrd/v4qW//V8+7Fiw0bV9/ObLh+ZfK9q+Elues9\n0uT6uudU7XoI6XA+3p8j3bfAbwgprseQpu8j9W8IbSZfvT3buW9ybqvxXLvbn9y/99/w/Odw\n/bP79e3Gd5QaN/ue+xb4DSHF9RTSeVeNZzZsZ2c2uM3p/BTS+dTczv4eH5CN3zucp3AY/zUs\nB1sAAACGSURBVOx+fZdLm+HSpvuO/XhfN26B3xCSPsxZhpgUfZizDDEpgABCAgQQEiCAkAAB\nhAQIICRAACEBAggJEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJEEBIgABCAgQQEiCAkAAB\nhAQIICRAACEBAggJEEBIgABCAgQQEiCAkAABhAQI+A9L3MO8zQW9GQAAAABJRU5ErkJggg==", "text/plain": [ "Plot with title \"Histogram of model1$residuals\"" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Suppose we wanted to inspect the residuals:\n", "hist(model1$residuals, breaks=50,xlim=c(-25000,25000))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Merging Data \n", "***\n", "## Now we've found another data source that includes unemployment by county in 2016. How can we combine this with our existing dataset, to perform a multiple linear regression?\n", "\n", "## This data is available here. Download it and save to your `raw data/` folder." ] }, { "cell_type": "code", "execution_count": 132, "metadata": {}, "outputs": [], "source": [ "temp<- read.csv(\"https://raw.githubusercontent.com/plotly/datasets/master/fips-unemp-16.csv\")\n", "write.csv(temp, paste0(projectpath,\"raw data/fips-unemp-16.csv\"))" ] }, { "cell_type": "code", "execution_count": 133, "metadata": {}, "outputs": [], "source": [ "df2<-read.csv(paste0(projectpath,\"raw data/fips-unemp-16.csv\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## We can perform the merge using the `merge` command:" ] }, { "cell_type": "code", "execution_count": 134, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
fipsX.xwageemploymentX.yunemp
1001 1 24803.96 8844 1 5.3
1003 2 21926.4540999 2 5.4
1005 3 22849.07 9890 3 8.6
1007 4 20039.41 2686 4 6.6
1009 5 21934.78 7372 5 5.5
\n" ], "text/latex": [ "\\begin{tabular}{r|llllll}\n", " fips & X.x & wage & employment & X.y & unemp\\\\\n", "\\hline\n", "\t 1001 & 1 & 24803.96 & 8844 & 1 & 5.3 \\\\\n", "\t 1003 & 2 & 21926.45 & 40999 & 2 & 5.4 \\\\\n", "\t 1005 & 3 & 22849.07 & 9890 & 3 & 8.6 \\\\\n", "\t 1007 & 4 & 20039.41 & 2686 & 4 & 6.6 \\\\\n", "\t 1009 & 5 & 21934.78 & 7372 & 5 & 5.5 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| fips | X.x | wage | employment | X.y | unemp |\n", "|---|---|---|---|---|---|\n", "| 1001 | 1 | 24803.96 | 8844 | 1 | 5.3 |\n", "| 1003 | 2 | 21926.45 | 40999 | 2 | 5.4 |\n", "| 1005 | 3 | 22849.07 | 9890 | 3 | 8.6 |\n", "| 1007 | 4 | 20039.41 | 2686 | 4 | 6.6 |\n", "| 1009 | 5 | 21934.78 | 7372 | 5 | 5.5 |\n", "\n" ], "text/plain": [ " fips X.x wage employment X.y unemp\n", "1 1001 1 24803.96 8844 1 5.3 \n", "2 1003 2 21926.45 40999 2 5.4 \n", "3 1005 3 22849.07 9890 3 8.6 \n", "4 1007 4 20039.41 2686 4 6.6 \n", "5 1009 5 21934.78 7372 5 5.5 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "merged_data <- merge(df, df2, by=\"fips\", all.x = TRUE, all.y=FALSE)\n", "merged_data[1:5,]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Note that `by` here can be a list of variables to merge on combinations of, e.g. `by=c(var1, var2)`.\n", "\n", "### The options `all.x = TRUE, all.y=FALSE` makes sure to keep all of our original county observations, even if we can't find an unemployment number for them. Counties that show up only in the unemployment dataset are ignored.\n", "\n", "### We can see that 9 counties were unmatched, because they have NA's for `unemp`:" ] }, { "cell_type": "code", "execution_count": 135, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " fips X.x wage employment \n", " Min. : 1001 Min. : 1.0 Min. : 8053 Min. : 19 \n", " 1st Qu.:19039 1st Qu.: 804.8 1st Qu.:20264 1st Qu.: 2265 \n", " Median :30028 Median :1608.5 Median :23167 Median : 6308 \n", " Mean :31464 Mean :1608.5 Mean :24230 Mean : 33644 \n", " 3rd Qu.:46114 3rd Qu.:2412.2 3rd Qu.:26846 3rd Qu.: 18226 \n", " Max. :78030 Max. :3216.0 Max. :79544 Max. :3548191 \n", " \n", " X.y unemp \n", " Min. : 1.0 Min. : 1.700 \n", " 1st Qu.: 809.5 1st Qu.: 4.000 \n", " Median :1611.0 Median : 5.000 \n", " Mean :1612.3 Mean : 5.456 \n", " 3rd Qu.:2415.5 3rd Qu.: 6.300 \n", " Max. :3219.0 Max. :23.500 \n", " NA's :9 NA's :9 " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
fipsX.xwageemploymentX.yunemp
86 2201 86 27511.05 1165 NA NA
88 2232 88 28550.83 1068 NA NA
91 2270 91 16724.99 878 NA NA
92 2280 92 25470.78 1736 NA NA
241246113 2412 22689.80 2244 NA NA
290951515 2909 25314.70 3970 NA NA
321478010 3214 28077.1212053 NA NA
321578020 3215 20178.25 1991 NA NA
321678030 3216 22985.7114565 NA NA
\n" ], "text/latex": [ "\\begin{tabular}{r|llllll}\n", " & fips & X.x & wage & employment & X.y & unemp\\\\\n", "\\hline\n", "\t86 & 2201 & 86 & 27511.05 & 1165 & NA & NA \\\\\n", "\t88 & 2232 & 88 & 28550.83 & 1068 & NA & NA \\\\\n", "\t91 & 2270 & 91 & 16724.99 & 878 & NA & NA \\\\\n", "\t92 & 2280 & 92 & 25470.78 & 1736 & NA & NA \\\\\n", "\t2412 & 46113 & 2412 & 22689.80 & 2244 & NA & NA \\\\\n", "\t2909 & 51515 & 2909 & 25314.70 & 3970 & NA & NA \\\\\n", "\t3214 & 78010 & 3214 & 28077.12 & 12053 & NA & NA \\\\\n", "\t3215 & 78020 & 3215 & 20178.25 & 1991 & NA & NA \\\\\n", "\t3216 & 78030 & 3216 & 22985.71 & 14565 & NA & NA \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| | fips | X.x | wage | employment | X.y | unemp |\n", "|---|---|---|---|---|---|---|\n", "| 86 | 2201 | 86 | 27511.05 | 1165 | NA | NA |\n", "| 88 | 2232 | 88 | 28550.83 | 1068 | NA | NA |\n", "| 91 | 2270 | 91 | 16724.99 | 878 | NA | NA |\n", "| 92 | 2280 | 92 | 25470.78 | 1736 | NA | NA |\n", "| 2412 | 46113 | 2412 | 22689.80 | 2244 | NA | NA |\n", "| 2909 | 51515 | 2909 | 25314.70 | 3970 | NA | NA |\n", "| 3214 | 78010 | 3214 | 28077.12 | 12053 | NA | NA |\n", "| 3215 | 78020 | 3215 | 20178.25 | 1991 | NA | NA |\n", "| 3216 | 78030 | 3216 | 22985.71 | 14565 | NA | NA |\n", "\n" ], "text/plain": [ " fips X.x wage employment X.y unemp\n", "86 2201 86 27511.05 1165 NA NA \n", "88 2232 88 28550.83 1068 NA NA \n", "91 2270 91 16724.99 878 NA NA \n", "92 2280 92 25470.78 1736 NA NA \n", "2412 46113 2412 22689.80 2244 NA NA \n", "2909 51515 2909 25314.70 3970 NA NA \n", "3214 78010 3214 28077.12 12053 NA NA \n", "3215 78020 3215 20178.25 1991 NA NA \n", "3216 78030 3216 22985.71 14565 NA NA " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "summary(merged_data)\n", "merged_data[is.na(merged_data$unemp),]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### On the rest of the data, let's do an OLS regression of wage on employment count and unemployment rate (counties with missing `unemp` will be ignored by the `lm` command):" ] }, { "cell_type": "code", "execution_count": 136, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "Call:\n", "lm(formula = wage ~ employment + unemp, data = merged_data)\n", "\n", "Residuals:\n", " Min 1Q Median 3Q Max \n", "-64582 -3441 -673 2457 48149 \n", "\n", "Coefficients:\n", " Estimate Std. Error t value Pr(>|t|) \n", "(Intercept) 2.565e+04 2.575e+02 99.590 <2e-16 ***\n", "employment 2.261e-02 7.934e-04 28.501 <2e-16 ***\n", "unemp -3.997e+02 4.294e+01 -9.308 <2e-16 ***\n", "---\n", "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n", "\n", "Residual standard error: 5631 on 3204 degrees of freedom\n", " (9 observations deleted due to missingness)\n", "Multiple R-squared: 0.2267,\tAdjusted R-squared: 0.2262 \n", "F-statistic: 469.6 on 2 and 3204 DF, p-value: < 2.2e-16\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "model2<-lm(wage~employment+unemp,data=merged_data)\n", "summary(model2)" ] }, { "cell_type": "code", "execution_count": 137, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\t
(Intercept)
\n", "\t\t
359.606291243014
\n", "\t
employment
\n", "\t\t
0.00495834913777904
\n", "\t
unemp
\n", "\t\t
50.0528772565186
\n", "
\n" ], "text/latex": [ "\\begin{description*}\n", "\\item[(Intercept)] 359.606291243014\n", "\\item[employment] 0.00495834913777904\n", "\\item[unemp] 50.0528772565186\n", "\\end{description*}\n" ], "text/markdown": [ "(Intercept)\n", ": 359.606291243014employment\n", ": 0.00495834913777904unemp\n", ": 50.0528772565186\n", "\n" ], "text/plain": [ " (Intercept) employment unemp \n", "3.596063e+02 4.958349e-03 5.005288e+01 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Again, if we are interested in statistical inference on our coefficients, we should generally use\n", "# standard errors that are valid under heteroskedasticity:\n", "cov <- vcovHC(model2, type = \"HC\")\n", "model2.robustses <- sqrt(diag(cov))\n", "model2.robustses" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aside: deleting, renaming, ordering variables in a dataframe \n", "\n", "### First, let's drop these meaningless variables X.x and X.y. These come from the fact that there was a variable named \"X\" in both datasets before the merge, which simply counted the row number." ] }, { "cell_type": "code", "execution_count": 138, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
fipsX.xwageemploymentX.yunemp
1001 1 24803.96 8844 1 5.3
1003 2 21926.4540999 2 5.4
1005 3 22849.07 9890 3 8.6
1007 4 20039.41 2686 4 6.6
1009 5 21934.78 7372 5 5.5
1011 6 20437.67 2774 6 7.2
\n" ], "text/latex": [ "\\begin{tabular}{r|llllll}\n", " fips & X.x & wage & employment & X.y & unemp\\\\\n", "\\hline\n", "\t 1001 & 1 & 24803.96 & 8844 & 1 & 5.3 \\\\\n", "\t 1003 & 2 & 21926.45 & 40999 & 2 & 5.4 \\\\\n", "\t 1005 & 3 & 22849.07 & 9890 & 3 & 8.6 \\\\\n", "\t 1007 & 4 & 20039.41 & 2686 & 4 & 6.6 \\\\\n", "\t 1009 & 5 & 21934.78 & 7372 & 5 & 5.5 \\\\\n", "\t 1011 & 6 & 20437.67 & 2774 & 6 & 7.2 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| fips | X.x | wage | employment | X.y | unemp |\n", "|---|---|---|---|---|---|\n", "| 1001 | 1 | 24803.96 | 8844 | 1 | 5.3 |\n", "| 1003 | 2 | 21926.45 | 40999 | 2 | 5.4 |\n", "| 1005 | 3 | 22849.07 | 9890 | 3 | 8.6 |\n", "| 1007 | 4 | 20039.41 | 2686 | 4 | 6.6 |\n", "| 1009 | 5 | 21934.78 | 7372 | 5 | 5.5 |\n", "| 1011 | 6 | 20437.67 | 2774 | 6 | 7.2 |\n", "\n" ], "text/plain": [ " fips X.x wage employment X.y unemp\n", "1 1001 1 24803.96 8844 1 5.3 \n", "2 1003 2 21926.45 40999 2 5.4 \n", "3 1005 3 22849.07 9890 3 8.6 \n", "4 1007 4 20039.41 2686 4 6.6 \n", "5 1009 5 21934.78 7372 5 5.5 \n", "6 1011 6 20437.67 2774 6 7.2 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "head(merged_data)" ] }, { "cell_type": "code", "execution_count": 141, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
fipswageemploymentunemp
1001 24803.96 8844 5.3
1003 21926.4540999 5.4
1005 22849.07 9890 8.6
1007 20039.41 2686 6.6
1009 21934.78 7372 5.5
1011 20437.67 2774 7.2
\n" ], "text/latex": [ "\\begin{tabular}{r|llll}\n", " fips & wage & employment & unemp\\\\\n", "\\hline\n", "\t 1001 & 24803.96 & 8844 & 5.3 \\\\\n", "\t 1003 & 21926.45 & 40999 & 5.4 \\\\\n", "\t 1005 & 22849.07 & 9890 & 8.6 \\\\\n", "\t 1007 & 20039.41 & 2686 & 6.6 \\\\\n", "\t 1009 & 21934.78 & 7372 & 5.5 \\\\\n", "\t 1011 & 20437.67 & 2774 & 7.2 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| fips | wage | employment | unemp |\n", "|---|---|---|---|\n", "| 1001 | 24803.96 | 8844 | 5.3 |\n", "| 1003 | 21926.45 | 40999 | 5.4 |\n", "| 1005 | 22849.07 | 9890 | 8.6 |\n", "| 1007 | 20039.41 | 2686 | 6.6 |\n", "| 1009 | 21934.78 | 7372 | 5.5 |\n", "| 1011 | 20437.67 | 2774 | 7.2 |\n", "\n" ], "text/plain": [ " fips wage employment unemp\n", "1 1001 24803.96 8844 5.3 \n", "2 1003 21926.45 40999 5.4 \n", "3 1005 22849.07 9890 8.6 \n", "4 1007 20039.41 2686 6.6 \n", "5 1009 21934.78 7372 5.5 \n", "6 1011 20437.67 2774 7.2 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "somevars <- names(merged_data) %in% c(\"X.x\", \"X.y\")\n", "cleaned_data <- merged_data[!somevars]\n", "head(cleaned_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Some other data frame manipulation tasks that may be useful:" ] }, { "cell_type": "code", "execution_count": 142, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
employmentemployment_quartilecountywageunemp
202819 [19,2.26e+03]38087 13620.74 2.5
162932 [19,2.26e+03]30069 16411.38 4.8
173150 [19,2.26e+03]31165 15092.88 2.8
161361 [19,2.26e+03]30037 15563.97 4.4
165361 [19,2.26e+03]31005 13003.02 6.2
165461 [19,2.26e+03]31007 19724.15 3.8
\n" ], "text/latex": [ "\\begin{tabular}{r|lllll}\n", " & employment & employment\\_quartile & county & wage & unemp\\\\\n", "\\hline\n", "\t2028 & 19 & {[}19,2.26e+03{]} & 38087 & 13620.74 & 2.5 \\\\\n", "\t1629 & 32 & {[}19,2.26e+03{]} & 30069 & 16411.38 & 4.8 \\\\\n", "\t1731 & 50 & {[}19,2.26e+03{]} & 31165 & 15092.88 & 2.8 \\\\\n", "\t1613 & 61 & {[}19,2.26e+03{]} & 30037 & 15563.97 & 4.4 \\\\\n", "\t1653 & 61 & {[}19,2.26e+03{]} & 31005 & 13003.02 & 6.2 \\\\\n", "\t1654 & 61 & {[}19,2.26e+03{]} & 31007 & 19724.15 & 3.8 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| | employment | employment_quartile | county | wage | unemp |\n", "|---|---|---|---|---|---|\n", "| 2028 | 19 | [19,2.26e+03] | 38087 | 13620.74 | 2.5 |\n", "| 1629 | 32 | [19,2.26e+03] | 30069 | 16411.38 | 4.8 |\n", "| 1731 | 50 | [19,2.26e+03] | 31165 | 15092.88 | 2.8 |\n", "| 1613 | 61 | [19,2.26e+03] | 30037 | 15563.97 | 4.4 |\n", "| 1653 | 61 | [19,2.26e+03] | 31005 | 13003.02 | 6.2 |\n", "| 1654 | 61 | [19,2.26e+03] | 31007 | 19724.15 | 3.8 |\n", "\n" ], "text/plain": [ " employment employment_quartile county wage unemp\n", "2028 19 [19,2.26e+03] 38087 13620.74 2.5 \n", "1629 32 [19,2.26e+03] 30069 16411.38 4.8 \n", "1731 50 [19,2.26e+03] 31165 15092.88 2.8 \n", "1613 61 [19,2.26e+03] 30037 15563.97 4.4 \n", "1653 61 [19,2.26e+03] 31005 13003.02 6.2 \n", "1654 61 [19,2.26e+03] 31007 19724.15 3.8 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Create a new variable that groups employment into 4 quartiles\n", "qnt <- quantile(cleaned_data$employment)\n", "cleaned_data$employment_quartile<-cut(cleaned_data$employment,unique(qnt),include.lowest=TRUE)\n", "\n", "#Sort the data frame by employment\n", "cleaned_data <- cleaned_data[order(cleaned_data$employment),]\n", "\n", "#Re-order first few columns in data frame.\n", "#In this exmaple, make \"employment\" and \"employment_quartile\" first and retain ordering of the rest\n", "somecols <- c(\"employment\", \"employment_quartile\");\n", "cols <- c(somecols, names(cleaned_data)[-which(names(cleaned_data) %in% somecols)]);\n", "cleaned_data <- cleaned_data[cols]\n", "\n", "#rename the \"fips\" variable to \"county\"\n", "colnames(cleaned_data)[which(colnames(cleaned_data) == 'fips')[[1]]] <- \"county\"\n", "\n", "head(cleaned_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### If you have experience with Stata, you'll notice that doing some of these things was more tedious than in Stata.\n", "\n", "### There are some tools that make data manipulation easier in R. One is the `dplyr` package, another is `data.table`\n", "\n", "### For example, with `data.table`, we can quickly extract the average employment across counties within each quartile of employment" ] }, { "cell_type": "code", "execution_count": 117, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
employment_quartilemean_emp_byqt
[19,2.26e+03] 1159.604
(2.26e+03,6.31e+03] 4009.684
(6.31e+03,1.82e+04] 10953.877
(1.82e+04,3.55e+06]118451.797
\n" ], "text/latex": [ "\\begin{tabular}{r|ll}\n", " employment\\_quartile & mean\\_emp\\_byqt\\\\\n", "\\hline\n", "\t {[}19,2.26e+03{]} & 1159.604 \\\\\n", "\t (2.26e+03,6.31e+03{]} & 4009.684 \\\\\n", "\t (6.31e+03,1.82e+04{]} & 10953.877 \\\\\n", "\t (1.82e+04,3.55e+06{]} & 118451.797 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| employment_quartile | mean_emp_byqt |\n", "|---|---|\n", "| [19,2.26e+03] | 1159.604 |\n", "| (2.26e+03,6.31e+03] | 4009.684 |\n", "| (6.31e+03,1.82e+04] | 10953.877 |\n", "| (1.82e+04,3.55e+06] | 118451.797 |\n", "\n" ], "text/plain": [ " employment_quartile mean_emp_byqt\n", "1 [19,2.26e+03] 1159.604 \n", "2 (2.26e+03,6.31e+03] 4009.684 \n", "3 (6.31e+03,1.82e+04] 10953.877 \n", "4 (1.82e+04,3.55e+06] 118451.797 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "library(data.table)\n", "dt = data.table(cleaned_data)\n", "dt <- dt[, .(mean(employment)), keyby = .(employment_quartile)]\n", "#the column with mean employment ends up with the non-descriptive name \"V1\", so let's change it\n", "colnames(dt)[[2]]<-\"mean_emp_byqt\"\n", "dt" ] }, { "cell_type": "code", "execution_count": 118, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
employment_quartileemploymentcountywageunempmean_emp_byqt
(1.82e+04,3.55e+06]18335 17177 30041.86 5.8 118451.8
(1.82e+04,3.55e+06]18434 48203 29760.17 6.1 118451.8
(1.82e+04,3.55e+06]18435 27131 28022.64 3.5 118451.8
(1.82e+04,3.55e+06]18449 36105 22923.15 4.8 118451.8
(1.82e+04,3.55e+06]18460 5093 27086.02 6.8 118451.8
(1.82e+04,3.55e+06]18534 48231 28889.91 4.3 118451.8
\n" ], "text/latex": [ "\\begin{tabular}{r|llllll}\n", " employment\\_quartile & employment & county & wage & unemp & mean\\_emp\\_byqt\\\\\n", "\\hline\n", "\t (1.82e+04,3.55e+06{]} & 18335 & 17177 & 30041.86 & 5.8 & 118451.8 \\\\\n", "\t (1.82e+04,3.55e+06{]} & 18434 & 48203 & 29760.17 & 6.1 & 118451.8 \\\\\n", "\t (1.82e+04,3.55e+06{]} & 18435 & 27131 & 28022.64 & 3.5 & 118451.8 \\\\\n", "\t (1.82e+04,3.55e+06{]} & 18449 & 36105 & 22923.15 & 4.8 & 118451.8 \\\\\n", "\t (1.82e+04,3.55e+06{]} & 18460 & 5093 & 27086.02 & 6.8 & 118451.8 \\\\\n", "\t (1.82e+04,3.55e+06{]} & 18534 & 48231 & 28889.91 & 4.3 & 118451.8 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| employment_quartile | employment | county | wage | unemp | mean_emp_byqt |\n", "|---|---|---|---|---|---|\n", "| (1.82e+04,3.55e+06] | 18335 | 17177 | 30041.86 | 5.8 | 118451.8 |\n", "| (1.82e+04,3.55e+06] | 18434 | 48203 | 29760.17 | 6.1 | 118451.8 |\n", "| (1.82e+04,3.55e+06] | 18435 | 27131 | 28022.64 | 3.5 | 118451.8 |\n", "| (1.82e+04,3.55e+06] | 18449 | 36105 | 22923.15 | 4.8 | 118451.8 |\n", "| (1.82e+04,3.55e+06] | 18460 | 5093 | 27086.02 | 6.8 | 118451.8 |\n", "| (1.82e+04,3.55e+06] | 18534 | 48231 | 28889.91 | 4.3 | 118451.8 |\n", "\n" ], "text/plain": [ " employment_quartile employment county wage unemp mean_emp_byqt\n", "1 (1.82e+04,3.55e+06] 18335 17177 30041.86 5.8 118451.8 \n", "2 (1.82e+04,3.55e+06] 18434 48203 29760.17 6.1 118451.8 \n", "3 (1.82e+04,3.55e+06] 18435 27131 28022.64 3.5 118451.8 \n", "4 (1.82e+04,3.55e+06] 18449 36105 22923.15 4.8 118451.8 \n", "5 (1.82e+04,3.55e+06] 18460 5093 27086.02 6.8 118451.8 \n", "6 (1.82e+04,3.55e+06] 18534 48231 28889.91 4.3 118451.8 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#We may want to merge this information back to the original dataset:\n", "cleaned_data <- merge(cleaned_data, dt, by=\"employment_quartile\")\n", "head(cleaned_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Suppose we wanted the average of the variables `employment` and `wage` across all counties. We could do this with the lapply funciton, which would allow us to loop through any set of variable names and ask for its mean:" ] }, { "cell_type": "code", "execution_count": 119, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] \"Mean of variable employment: 33643.7406716418\"\n", "[1] \"Mean of variable wage: 24229.9502373469\"\n" ] } ], "source": [ "nothing<-lapply(c(\"employment\", \"wage\"), function(thisvar){\n", " print(paste0(\"Mean of variable \", thisvar, \": \",mean(cleaned_data[,thisvar])))\n", "})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# Now that we have several models estimated, how can we neatly display all the output?\n", "***\n", "## This can be done with the `stargazer` libary:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "=================================================================================================\n", " Dependent variable: \n", " -----------------------------------------------------------------------------\n", " wage \n", " not robust robust robust \n", " (1) (2) (3) \n", "-------------------------------------------------------------------------------------------------\n", "employment 0.023*** 0.023*** 0.023*** \n", " (0.001) (0.005) (0.005) \n", " \n", "unemp -399.665*** \n", " (50.053) \n", " \n", "Constant 23,452.380*** 23,452.380*** 25,647.870*** \n", " (104.080) (171.421) (359.606) \n", " \n", "-------------------------------------------------------------------------------------------------\n", "Observations 3,216 3,216 3,207 \n", "R2 0.206 0.206 0.227 \n", "Adjusted R2 0.205 0.205 0.226 \n", "Residual Std. Error 5,700.850 (df = 3214) 5,700.850 (df = 3214) 5,630.677 (df = 3204) \n", "F Statistic 831.685*** (df = 1; 3214) 831.685*** (df = 1; 3214) 469.627*** (df = 2; 3204)\n", "=================================================================================================\n", "Note: *p<0.1; **p<0.05; ***p<0.01\n" ] } ], "source": [ "#install.packages(\"stargazer\")\n", "library(\"stargazer\")\n", "stargazer(model1, model1, model2, se=list(NULL, model1.robustses, model2.robustses), column.labels=c(\"not robust\",\"robust\", \"robust\"), align=TRUE, type=\"text\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### We'll also save a LaTeX version of this table into a file, to use in the next lecture on LaTeX:" ] }, { "cell_type": "code", "execution_count": 125, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "% Table created by stargazer v.5.2.2 by Marek Hlavac, Harvard University. E-mail: hlavac at fas.harvard.edu\n", "% Date and time: Wed, Oct 30, 2019 - 6:00:38 PM\n", "% Requires LaTeX packages: dcolumn \n", "\\begin{tabular}{@{\\extracolsep{5pt}}lD{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} } \n", "\\\\[-1.8ex]\\hline \n", "\\hline \\\\[-1.8ex] \n", " & \\multicolumn{3}{c}{\\textit{Dependent variable:}} \\\\ \n", "\\cline{2-4} \n", "\\\\[-1.8ex] & \\multicolumn{3}{c}{wage} \\\\ \n", " & \\multicolumn{1}{c}{not robust} & \\multicolumn{1}{c}{robust} & \\multicolumn{1}{c}{robust} \\\\ \n", "\\\\[-1.8ex] & \\multicolumn{1}{c}{(1)} & \\multicolumn{1}{c}{(2)} & \\multicolumn{1}{c}{(3)}\\\\ \n", "\\hline \\\\[-1.8ex] \n", " employment & 0.023^{***} & 0.023^{***} & 0.023^{***} \\\\ \n", " & (0.001) & (0.005) & (0.005) \\\\ \n", " & & & \\\\ \n", " unemp & & & -399.665^{***} \\\\ \n", " & & & (50.053) \\\\ \n", " & & & \\\\ \n", " Constant & 23,452.380^{***} & 23,452.380^{***} & 25,647.870^{***} \\\\ \n", " & (104.080) & (171.421) & (359.606) \\\\ \n", " & & & \\\\ \n", "\\hline \\\\[-1.8ex] \n", "Observations & \\multicolumn{1}{c}{3,216} & \\multicolumn{1}{c}{3,216} & \\multicolumn{1}{c}{3,207} \\\\ \n", "R$^{2}$ & \\multicolumn{1}{c}{0.206} & \\multicolumn{1}{c}{0.206} & \\multicolumn{1}{c}{0.227} \\\\ \n", "Adjusted R$^{2}$ & \\multicolumn{1}{c}{0.205} & \\multicolumn{1}{c}{0.205} & \\multicolumn{1}{c}{0.226} \\\\ \n", "Residual Std. Error & \\multicolumn{1}{c}{5,700.850 (df = 3214)} & \\multicolumn{1}{c}{5,700.850 (df = 3214)} & \\multicolumn{1}{c}{5,630.677 (df = 3204)} \\\\ \n", "F Statistic & \\multicolumn{1}{c}{831.685$^{***}$ (df = 1; 3214)} & \\multicolumn{1}{c}{831.685$^{***}$ (df = 1; 3214)} & \\multicolumn{1}{c}{469.627$^{***}$ (df = 2; 3204)} \\\\ \n", "\\hline \n", "\\hline \\\\[-1.8ex] \n", "\\textit{Note:} & \\multicolumn{3}{r}{$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01} \\\\ \n", "\\end{tabular} \n" ] } ], "source": [ "stargazer(model1, model1, model2, se=list(NULL, model1.robustses, model2.robustses), column.labels=c(\"not robust\",\"robust\", \"robust\"), align=TRUE, float=F, type=\"latex\", out=paste0(projectpath,\"/tables/table1.tex\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "

Project!

\n", "\n", "## Suppose we want to add to this regression a \"fixed effect\" for each state. See if you can figure out how to do this.\n", "\n", "## Some hints:\n", "- the first two digits of a FIPS code are constant within a state\n", "- dummies for each value of a variable can be added to the regression by creating a \"factor variable\", e.g. `df$state.f <- factor(df$state)`\n", "\n", "\n", "

\"Extra credit\"!

\n", "\n", "## If the first task was too easy, now try adding to the regression a dummy variable for the first letter of the county's name beginning with a \"W\". County names can be downloaded here: https://www2.census.gov/programs-surveys/popest/geographies/2017/all-geocodes-v2017.xlsx" ] } ], "metadata": { "kernelspec": { "display_name": "R", "language": "R", "name": "ir" }, "language_info": { "codemirror_mode": "r", "file_extension": ".r", "mimetype": "text/x-r-source", "name": "R", "pygments_lexer": "r", "version": "3.6.1" } }, "nbformat": 4, "nbformat_minor": 2 }