--- title: "R Basics for Lecture 1" output: html_document --- ## Econ B2000, Statistics and Introduction to Econometrics ## Kevin R Foster, Colin Powell School, the City College of New York, CUNY For this class we'll be using the statistical analysis program, R. The computer labs here on campus have the necessary software. If you have your own computer, download and install the program from [R Project](http://www.r-project.org/) along with [R-Studio](http://www.rstudio.com/). Depending on your particular device these details will be different so take it patiently and figure it out - there are sufficient online sources for help, to address any problem that arises. The best way to learn it is to just do it. How many times have you got a new game, skipped reading the instructions, and learned by crashing it a few times? This isn't quite so simple but the basic gist remains - try it. Below I give you some pointers about how to just do it, but they're meant to be read then immediately done in real life, I'll give you some commands to just copy-and-paste into the program. The program, R, is the machine underneath that does the work, while R-Studio is a skin on top that makes it easier to use. Install both then run R-Studio, and you'll get something that looks like this: ![image of R Studio](screenshot_RStudio1.png) The screen has 4 parts: the Console at lower left (where I drew the green arrow) is most important since that's where you type in the commands to R. You can start by just copying and pasting commands from this help into the "Console" and seeing the output from the program. The guide, *An Introduction to R*, suggests these commands to give a basic flavor of what's done and hint at the power. Don't worry too much about each step for now, this is just a glimpse to show that with a few commands you can do some relatively sophisticated estimation -- i.e. for a small cost you can get a large benefit. Copy them and paste into the "Console". ```{r message=FALSE} x <- 1:50 w <- 1 + sqrt(x)/2 example1 <- data.frame(x=x, y= x + rnorm(x)*w) attach(example1) ``` This creates x and y variables (where the rnorm command creates random numbers from a normal distribution), puts them into a data frame, then attaches that data frame so that R can use it in later calculations. (What's a data frame? A way to collect together related data items.) Next some stats - create a linear model (that's the "lm") then a "lowess" nonparametric local regression, and plot the two estimations to compare. (Just copy and paste these, don't worry about understanding them for now!) ```{r} fm <- lm(y ~ x) summary(fm) lrf <- lowess(x, y) plot(x, y) lines(x, lrf$y) abline(0, 1, lty=3) abline(coef(fm)) detach() ``` Your own graph should look similar although with the random numbers, not exactly. The final "detach" command just cleans up, it is the opposite of "attach". Later you will find that while it is possible to go back and fix up code that you previously ran in R (the up arrow brings most recent commands), it is a bit of a pain. It is easier to write out the code in the upper-left panel, then R-Studio can obediently run all or part of it in R for you. Then if you make a mistake or just want to do it again, it's a bit easier. You can save your code for next time, as well. One big thing to learn is that while pushing buttons is easier at first, eventually you want to be able to write code. I won't insist that you do that from week #1, but keep it in mind as a goal to work towards. In R, you can usually see the code generated by the button-pushes so you can learn. For all of these commands, you can use the help box on R Studio. But as I said, don't worry much about those commands for now, I'm not expecting you to become an R-ninja overnight. # With Some Data Next we will go through some basic stats with R using the Census Bureau's PUMS (Public Use Microdata Sample, from the American Community Survey, accessed from IPUMS). I have restricted the sample to contain only people living in the state of New York, a sample of 196,585. The full sample, with 3,190,040 observations, is also on the class website if you'd like -- although you might find that it slows your computer quite a bit! You have to learn a bit about your computer's filing system. When you download stuff it probably goes into a Download folder. When you installed R and R-Studio, those went into their own folders (probably within some Application folder) and then the program might have created a R folder for your work. When you do homework or other projects, those might (well, really should) get their own folder. But joining those up is a bit of housekeeping. If you read the pros, you'll get more tips on organization. Go and download the PUMS data from the [class page](http://kfoster.ccny.cuny.edu/classes/fall2019/), that will likely put that zip file into your Downloads folder. Within that zip file is one particular file, acs2017_ny_data.RData - move that into some appropriately-titled folder for this project. (Personally I might use the default "R Projects" folder then create "PUMS Data" as a folder within that, but you can choose. Young people with supple minds can get by even if they have zero organization like a hippie, but later you might appreciate being tidier.) Note that you can't just download the zip file and run the program, you have to extract the particular .RData file -- that regularly trips up people as they begin. R is looking for files in a particular directory. Type the command, "getwd()" to see where it's currently looking. Then use the command, "setwd," to tell it where it ought to be looking (the R Projects/PUMS Data folder that you created just now, or wherever you put that data file). Alt you can click "Session" then "Set Working Directory" then "Choose Directory" and click the folder, and that will insert the "setwd" command onto the Console line. Tell R to look for the data in the folder where you put the data. If you have a Dropbox (or Git or whatever) account then that can be worked in the same way, if you set that up as a folder on your computer. Then run these commands (output from those commands is below), ```{r} rm(list = ls(all = TRUE)) # clear workspace # change this command to your own directory where you put the data file # setwd("C:\\Users\\kevin\\Documents\\R") # ^^^^^^ load("acs2017_ny_data.RData") #glimpse(acs2017_ny) try this later acs2017_ny[1:10,1:7] attach(acs2017_ny) ``` In the next section I'll explain more about the data and what the lines mean but for now AGE is the person's age in years and female is a 0/1 variable (takes value 1 for true and 0 for false). So if you look at the output, see that the first person on line 1 is a 72-year-old female, next is a 72-year-old male, then a 31-year-old male, etc. The set of variables beginning "educ_" are also 0/1 so the first 2 people have advanced degrees and next 2 have college degrees. You can also use the command, summary, to find out about data. ```{r} summary(acs2017_ny) print(NN_obs <- length(AGE)) ``` So this shows that there are `r prettyNum(NN_obs, big.mark = ",")` people in this dataset. ###Simple Stats We compare the average age of the men and the women in the data, ```{r} summary(AGE[female == 1]) summary(AGE[!female]) ``` This uses the female dummy variable. The comparison is between those who have the variable female=1 (i.e. women) and those not female=1 (so logical not, denoted with the "!" symbol, i.e. men). *I know, you can -- and people do! -- worry that this binary classification for gender misses some people; government statistics are just not there yet. Progress is slower than we'd like.* Women in this dataset are, on average, a bit older, with an average age of `r round(mean(AGE[female == 1]), digits = 1)` compared with `r round(mean(AGE[female != 1]), digits = 1)` for men. You might wonder (if you were to begin to think like a statistician) whether that is a big difference - hold onto that thought! Alternately you can use the commands to calculate the average, mean(), and the standard deviation, sd(), to get those statistics: ```{r} # here i want to find average ages of men and women mean(AGE[female == 1]) sd(AGE[female == 1]) mean(AGE[!female]) sd(AGE[!female]) ``` Later you might encounter cases where you want more complicated dummy variables and want to use logical relations "and" "or" "not" (the symbols "&", "|", "!") or the ">=" or multiplication or division. As you're going along, if you're copying-and-pasting then you might not have had trouble, but if you're typing then you've probably realized that R is persnickety - the variable AGE is not equivalent to Age nor age nor aGe... You might be wondering if there's an easier way than copy-paste; *there is*. If you highlight text in the Source pane and then hit CTRL-Enter, that will run the block of code. Alt if you click the "Knit HTML" button from the Rmd file then that will recreate the webpage as well as run all of the commands. ###More Details At the beginning you'll be just copying commands that I give you but as you start making changes, you have to understand what's going on. For better or worse there is rarely just a single way of doing things, when people talk about the R language it does indeed have many features of a language. Just like a spoken language has many ways of saying "hi" so too there are lots of different ways to ask R to produce means of different groups. I am showing you a few particular ways but if you look around online you'll find others. Remember to save your work. The computers in the lab wipe the memory clean when you log off so back up your files. Either online (email it to yourself or upload to Google Drive or iCloud or Dropbox or whatever) and/or use a USB drive. But be careful - you have the option of saving your code and/or your workspace; it is usually better to save the code. For example, your code might load some data (at this stage, usually data that I've provided). If you make changes to this data during your session, you might not want to restart with those changes. If you just save your code and re-run the program, you will start fresh. If instead you save the workspace, then you end up saving all of your scraps and changes - which eventually start to build up. So only save your workspace if you've run a large chunk of code that you don't want to have to later re-run. (This is even more true if you "attach" data frames without later "detach"-ing them!) There are 2 main file types: R script (lines of instructions to R on how to calculate some statistics) and R Markdown (which combines the instructions to R along with text around it, which is how I created this file). There are more types of course but those are the ones for now. To create files like the one you're reading now, use "File \\ New File \\ R Markdown..." Give it a title and save the file with some name that you'll remember and the extension, ".Rmd". Remember the subdirectory or folder. There is a bigger philosophy behind this: the idea of "reproducible research," that if you share your .Rmd file then anybody else can run the exact same program and get the exact same results. It's a way of convincing an audience that there was no skulduggery in how the data relates to the conclusion. It's useful in class since it's an easy way to submit homework or share work with your study group. ## Variable Coding Some of the PUMS variables here have a natural interpretation, for instance Age is measured in years. Actually even this has a bit of a twist, look at the histogram. ```{r} hist(AGE[(AGE > 90)]) ``` There is a bit of weirdness in the right, where it looks like there are suddenly a bunch of people who are 95 but nobody is 94 or 96. This is due to a coding choice by the Census, where really old people are just labeled as "95" (top-coding) so it actually should be interpreted as meaning "92 or older". So if you were to get finicky (and every good statistician is!) you might go back to the calculations of averages previously and modify them all like this, to select just those who are female and who are coded as having age less than 90. Many variables are topcoded! *And recall that topcoding wouldn't change the median values calculated before, which is a point in favor of that statistic.* ```{r} mean(AGE[ (female == 1) & (AGE<90) ]) ``` You go make those other changes, figure out how top-coding changes the calculations of average age by gender -- I'll wait right here... ## Variable Coding Again So we were saying that some variables, like Age - ahem! -- have a natural interpretation as a number. Others are logical variables (called dummies) like female, Hispanic, or married - there is a yes/no answer that is coded 1/0. *Note that if you're creating these on your own it's good to give names that have that sort of yes/no answer, so a variable named 'female' is better than one named 'gender' where you'd have to remember who are coded adtrue and who are false.* Many variables, like PUMA, have no natural explanation at all. Here are the first codes, ```{r} str(as.numeric(PUMA)) ``` You have to go to the codebook (or, in this case, the file PUMA_levels.csv or acs2017_codebook.txt from the zip file) to find out that 3801 codes for Washington Heights/Inwood, 3802 is Hamilton Heights/Manhattanville/West Harlem, etc. The program will happily calculate the average value for PUMA (type in *mean(PUMA)* and see for yourself!) but this is a meaningless value -- wtf is the average neighborhood code value!? If you want to select just people living in a particular neighborhood then you'd have to look at the list below. PUMA | Neighborhood ---- | ------------ 3701 | NYC-Bronx CD 8--Riverdale, Fieldston & Kingsbridge 3702 | NYC-Bronx CD 12--Wakefield, Williamsbridge & Woodlawn 3703 | NYC-Bronx CD 10--Co-op City, Pelham Bay & Schuylerville 3704 | NYC-Bronx CD 11--Pelham Parkway, Morris Park & Laconia 3705 | NYC-Bronx CD 3 & 6--Belmont, Crotona Park East & East Tremont 3706 | NYC-Bronx CD 7--Bedford Park, Fordham North & Norwood 3707 | NYC-Bronx CD 5--Morris Heights, Fordham South & Mount Hope 3708 | NYC-Bronx CD 4--Concourse, Highbridge & Mount Eden 3709 | NYC-Bronx CD 9--Castle Hill, Clason Point & Parkchester 3710 | NYC-Bronx CD 1 & 2--Hunts Point, Longwood & Melrose 3801 | NYC-Manhattan CD 12--Washington Heights, Inwood & Marble Hill 3802 | NYC-Manhattan CD 9--Hamilton Heights, Manhattanville & West Harlem 3803 | NYC-Manhattan CD 10--Central Harlem 3804 | NYC-Manhattan CD 11--East Harlem 3805 | NYC-Manhattan CD 8--Upper East Side 3806 | NYC-Manhattan CD 7--Upper West Side & West Side 3807 | NYC-Manhattan CD 4 & 5--Chelsea, Clinton & Midtown Business District 3808 | NYC-Manhattan CD 6--Murray Hill, Gramercy & Stuyvesant Town 3809 | NYC-Manhattan CD 3--Chinatown & Lower East Side 3810 | NYC-Manhattan CD 1 & 2--Battery Park City, Greenwich Village & Soho 3901 | NYC-Staten Island CD 3--Tottenville, Great Kills & Annadale 3902 | NYC-Staten Island CD 2--New Springville & South Beach 3903 | NYC-Staten Island CD 1--Port Richmond, Stapleton & Mariner's Harbor 4001 | NYC-Brooklyn CD 1--Greenpoint & Williamsburg 4002 | NYC-Brooklyn CD 4—Bushwick 4003 | NYC-Brooklyn CD 3--Bedford-Stuyvesant 4004 | NYC-Brooklyn CD 2--Brooklyn Heights & Fort Greene 4005 | NYC-Brooklyn CD 6--Park Slope, Carroll Gardens & Red Hook 4006 | NYC-Brooklyn CD 8--Crown Heights North & Prospect Heights 4007 | NYC-Brooklyn CD 16--Brownsville & Ocean Hill 4008 | NYC-Brooklyn CD 5--East New York & Starrett City 4009 | NYC-Brooklyn CD 18--Canarsie & Flatlands 4010 | NYC-Brooklyn CD 17--East Flatbush, Farragut & Rugby 4011 | NYC-Brooklyn CD 9--Crown Heights South, Prospect Lefferts & Wingate 4012 | NYC-Brooklyn CD 7--Sunset Park & Windsor Terrace 4013 | NYC-Brooklyn CD 10--Bay Ridge & Dyker Heights 4014 | NYC-Brooklyn CD 12--Borough Park, Kensington & Ocean Parkway 4015 | NYC-Brooklyn CD 14--Flatbush & Midwood 4016 | NYC-Brooklyn CD 15--Sheepshead Bay, Gerritsen Beach & Homecrest 4017 | NYC-Brooklyn CD 11--Bensonhurst & Bath Beach 4018 | NYC-Brooklyn CD 13--Brighton Beach & Coney Island 4101 | NYC-Queens CD 1--Astoria & Long Island City 4102 | NYC-Queens CD 3--Jackson Heights & North Corona 4103 | NYC-Queens CD 7--Flushing, Murray Hill & Whitestone 4104 | NYC-Queens CD 11--Bayside, Douglaston & Little Neck 4105 | NYC-Queens CD 13--Queens Village, Cambria Heights & Rosedale 4106 | NYC-Queens CD 8--Briarwood, Fresh Meadows & Hillcrest 4107 | NYC-Queens CD 4--Elmhurst & South Corona 4108 | NYC-Queens CD 6--Forest Hills & Rego Park 4109 | NYC-Queens CD 2--Sunnyside & Woodside 4110 | NYC-Queens CD 5--Ridgewood, Glendale & Middle Village 4111 | NYC-Queens CD 9--Richmond Hill & Woodhaven 4112 | NYC-Queens CD 12--Jamaica, Hollis & St. Albans 4113 | NYC-Queens CD 10--Howard Beach & Ozone Park 4114 | NYC-Queens CD 14--Far Rockaway, Breezy Point & Broad Channel Now you're probably thinking, isn't there some easier way? Yes there is. R has variables called "factors" that join together the long list of codes with a separate file telling what those codes mean. Later when we do further statistics, R will know how to appropriately treat these factors. (Also it will then give an error if you calculate mean(PUMA), which is proper.) ```{r} PUMA <- as.factor(PUMA) female <- as.factor(female) ``` I will leave you to worry over the recoding of the other variables, because it's good for the soul. I will show you 2 ways -- the quick and dirty way, and the fancy correct way. First the quick and dirty way. ```{r} print(levels(female)) levels(female) <- c("male","female") ``` Well, ways, ```{r} educ_indx <- factor((educ_nohs + 2*educ_hs + 3*educ_somecoll + 4*educ_college + 5*educ_advdeg), levels=c(1,2,3,4,5),labels = c("No HS","HS","SmColl","Bach","Adv")) ``` (If you can figure out how that bit of code works, that would be good) These just type in the levels. But for things like PUMA, it could be a long list and might not even match every one. To do it better, we need help from an R package. ###Detour on Packages But first a bit of a detour, to mention how to use packages. R depends crucially on "packages" - that's the whole reason that the open-source works. Some statistician invents a cool new technique, then writes up the code in R and makes it available. If you used a commercial program you'd have to wait a decade for them to update it; in R it's here now. Also if somebody hacks a nicer or easier way to do stuff, they write it up. So enter this into the Console, ``` install.packages("tidyverse") install.packages("plyr") ``` then ```{r message=FALSE} library(tidyverse) library(plyr) levels_n <- read.csv("PUMA_levels.csv") levels_orig <- levels(PUMA) levels_new <- join(data.frame(levels_orig),data.frame(levels_n)) levels(PUMA) <- levels_new$New_Level ``` Alt, from R-Studio, click "Tools" then "Install Packages..." and tell it to install the packages, "plyr" and "tidyverse". That is nice if you want to see some of the packages or if you don't quite remember the name. Then the next piece of code, library, tells the program that you want to use commands from this package. Those commands read in a little csv file that I had made, with the PUMA codes, then matches the old codes with the new complete text. Note that I'm lazy so codes in NY state outside of NYC are coded NA. ###Back from Detour R will do the summary differently when it knows the variable is a factor, ```{r} summary(female) summary(PUMA) summary(educ_indx) ``` To find mean and standard deviation by neighborhood, you could use something like this, ```{r} ddply(acs2017_ny, .(PUMA), summarize, mean = round(mean(AGE), 2), sd = round(sd(AGE), 2)) ``` Although tapply would also work fine. Here's the 90th and 10th percentiles of wages by neighborhood, ```{r} dat_use1 <- subset(acs2017_ny,((INCWAGE > 0) & in_NYC)) ddply(dat_use1, .(PUMA), summarize, inc90 = quantile(INCWAGE,probs = 0.9), inc10 = quantile(INCWAGE,probs = 0.1), n_obs = length(INCWAGE)) ``` You could also use table (or crosstabs) for factors with fewer items, ```{r} table(educ_indx,female) xtabs(~educ_indx + female) ``` Want proportions instead of counts? ```{r} prop.table(table(educ_indx,female)) ``` *Remember prop.table later when we do marginals.* Try it and see what happens if you use table with PUMA... This data includes not just whether a person has a college degree but also what field was the degree in: Economics or Psychology, for instance. Look over the codebook about DEGFIELD and DEGFIELDD (that second D means more detail) to see the codes. Maybe look at 10th and 90th percentiles by degree field? In general, R is very flexible so there are often many different ways to get the same answer. There are some people who love to debate which is best. (Often, tradeoff between speed and intelligibility.) For now just worry about learning at least one way. Later on you can go back and refine your techniques. Sometimes attaching a dataset makes things easier. But as you get more advanced you might find it better to include the dataset name inside the function. There are advantages and disadvantages each way and some of the intro texts suggest one or the other. If you do a lot of analysis on a particular subgroup, it might be worthwhile to create a subset of that group, so that you don't have to always add on logical conditions. These two sets of expressions, looking at "prime-age" people, get the same results: ```{r message=FALSE} mean(educ_nohs[(AGE >= 25)&(AGE <= 55)]) mean(educ_hs[(AGE >= 25)&(AGE <= 55)]) mean(educ_somecoll[(AGE >= 25)&(AGE <= 55)]) mean(educ_college[(AGE >= 25)&(AGE <= 55)]) mean(educ_advdeg[(AGE >= 25)&(AGE <= 55)]) # alternatively restrict1 <- as.logical((AGE >= 25)&(AGE <= 55)) dat_age_primeage <- subset(acs2017_ny, restrict1) detach() attach(dat_age_primeage) mean(educ_nohs) mean(educ_hs) mean(educ_somecoll) mean(educ_college) mean(educ_advdeg) detach() ``` So you detach the original data frame and instead attach the restricted version. Then any subsequent analysis would be just done on that subset. Just remember that you've done this (again, this is a good reason to save the commands in a program so you can look back) otherwise you'll wonder why you suddenly don't have any kids in the sample! ## Why All These Details? You might be tired and bored by these details, but note that there are actually important choices to be made here, even in simply defining variables. Take the fraught American category of "race". This data has a variable, RACED, showing how people chose to classify themselves, as 'White,' 'Black,' 'American Indian,' 'Asian,' various combinations, and many more codes. Suppose you wanted to find out how many Asians are in a particular population. You could count how many people identify themselves as Asian only; you could count how many people identify as Asian in any combination. Sometimes the choice is irrelevant; sometimes it can skew the final results (e.g. the question in some areas, are there more African-Americans or more Hispanics?). Again, there's no "right" way to do it because there's no science in this peculiar-but-popular concept of "race". People's conceptions of themselves are fuzzy and complicated; these measures are approximations. ## Basics of government race/ethnicity classification The US government asks questions about people's race and ethnicity. These categories are social constructs, which is a fancy way of pointing out that they are based on people's own views of themselves (influenced by how we think that other people think of us...). Currently the standard classification asks people separately about their "race" and "ethnicity" where people can pick labels from each category in any combination. The "race" categories include: "White," "African-American," "American Indian," "Asian," and others. Then the supplemental race categories offer more detail. These are a peculiar combination of very general (well over 40% of the world's population is "Asian") and very specific ("American Indian") representing a peculiar history of popular attitudes in the US. Only in the 2000 Census did they start to classify people in mixed races. If you were to go back to historical US Censuses from more than a century ago, you would find that the category "race" included separate entries for Irish and French. Stephen J Gould has a fascinating book, The Mismeasure of Man, discussing how early scientific classifications of humans tried to "prove" which nationalities/races/groups were the smartest. Ta-Nehisi Coates notes, "racism invented race in America." Note that "Hispanic" is not "race" but rather ethnicity (includes various other labels such as Spanish, Latino, etc.). So a respondent could choose "Hispanic" and any race category -- some choose "White," some choose "African American," some might be combined with any other of those complicated racial categories. If you wanted to create a variable for those who report themselves as African-American and Hispanic, you'd use the expression (AfAm == 1) & (Hispanic == 1); sometimes stats report for non-Hispanic whites so (white == 1) & (Hispanic != 1). You can create your own classifications depending on what questions you're investigating. This data includes items on birthplace and ancestry (more detail on relatives!). The Census Bureau gives more information [here](http://www.census.gov/newsroom/minority_links/minority_links.html). All of these racial categories make some people uneasy: is the government encouraging racism by recognizing these classifications? Some other governments choose not to collect race data. But that doesn't mean that there are no differences, only that the government doesn't choose to measure any of these differences. In the US, government agencies such as the Census and BLS don't generally collect data on religion. ### Re-Coding complicated variables from initial data If we want more combinations of variables then we create those. Usually a statistical analysis spends a lot of time doing this sort of housekeeping - dull but necessary. It has a variety of names: data carpentry, data munging... Educational attainment is also classified with complicated codes in this data: the original data has code 63 to mean high school diploma, 64 for a GED, 65 for less than a year of college, etc. I have transformed them into a series of dummy variables, zero/one variables for whether a person has no high school diploma, just a high school diploma, some college, a bachelor's degree, or an advanced degree. As with race, there is no rule that you must always do it thus. Your own analysis might be different. (Is it valid to lump together everybody who lacks a high school diploma, no matter whether they completed just first grade or up to 12th?) That's the whole point of learning to do the data work for yourself: you can see all of the little decisions that go into creating a conclusion. Some conclusions might be fragile so a tiny decision about coding could change everything; other conclusions are robust to deviations. You must find out. # De-bugging Without a doubt, programming is tough. In R or with any other program, it is frustrating and complicated and difficult to do it the first few times. Some days it feels like a continuous battle just to do the simplest thing! Keep going despite that, keep working on it. Your study group will be very helpful of course. Often a google search of the error message helps. If you've isolated the error and read the help documentation on that command, then you're on your way to solving the problem on your own. If you have troubles that you can't solve, email me for help. But try to narrow down your question: if you run 20 lines of code that produce an error, is there a way to reproduce the error after just 5 lines? What if you did the same command on much simpler data, would it still cause an error? Sending emails like "I have a problem with errors" might be cathartic but is not actually useful to anyone. If you email me with the minimal code that recreates the error, along with the text of the error and/or a screenshot, then that will help more. ##Do it The first homework assignment asks you to start working on these questions. Begin by running the code that I give here, just to see if you can replicate my results. Then start asking more questions about the data - there is so much info there! Immigration/ancestry status, how they commute, health insurance, poverty, income, owner or renter (and how much they pay!), all sorts of info. Have some fun.