---
title: "IPUMS"
author: "Robby"
date: "March 27, 2020"
output: word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(ipumsr)
#library(survey) # Some easy weighting functions.
library(srvyr) # Like survey but with syntax like dplyr.
library(labelled) # for working with labelled data
```

## Objective 

To get up and running with IPUMS.

## Get the data

You need to go to the ipums website to get the data. Login with the credentials and click on My Data. Follow the instructions [here](https://cran.r-project.org/web/packages/ipumsr/vignettes/ipums.html).

### About the data extract

As of 3/27/2020, the data extract is as follows: 

- ACS from 2012-2018. That's seven years. Usually, we do 5-year samples but I grabbed two more in case we want to increase the sample size for some breakdowns. 
- Entire state of Louisiana. We will use PUMA codes to filter down to the region or its parishes as necessary/possible. Some PUMAs cover multiple parishes.
- A wide assortment of housing, demographic, economic, and workforce variables. 


### Load and clean

Basically, you load a "ddi" xml file that has the codes and metadata, and then load in the .dat file that actually has the data. If you open the DDI file in a browser, it's a handy codebook. For some reason, it doesn't work when you open it locally, but if you open it from IPUMS website, you can just keep it open in your browser for reference. 

```{r}

# Note that you can pass in the loaded DDI into the `read_ipums_micro()`
# NOTE: To load data, you must download both the extract's data and the DDI
# and also set the working directory to the folder with these files (or change the path below).

if (!require("ipumsr")) stop("Reading IPUMS data into R requires the ipumsr package. It can be installed using the following command: install.packages('ipumsr')")

ddi <- read_ipums_ddi("inputs/ipums/usa_00007.xml")
ipums.la <- read_ipums_micro(ddi)
```

This pulls in *labelled* data. The raw data is just numeric codes for different characteristics, but the labels make it easier because you don't have to go back. It's a bit of a trade-off though. If memory serves, the labelled fields are a bit weird and sometimes don't play nicely with dplyr or other common packages. 

Let's clean up the labels to make them easier to use. The default is a "labelled" -- it includes both the original codes (which you can find in the DDI file data dictionary) and the labelles attached the code. For example, for HHTYPE, it will include the code "1" and the label "Married-couple family household" and all of the other code/label pairs. As a convenience I use `labelled::to_factor` to convert the labels to factors. This way, it just shows the verbose labels like "Married-couple family household". More importantly, I've had issues with labelled data causing conflicts with dplyr or srvyr.

Also, be aware that some of the continuous variables have special codes. For example, the income variables have top coding, so they're not pure double/integer variables. 

```{r}
# Note this list of variables converted to factors may be incomplete.

# to_factor pulls the labels out and makes them factors.
# This is an alternative to the hybrid int+lbl default. 
ipums.la.fac <- ipums.la %>% 
  mutate_at(vars(NUMPREC:SUBSAMP,
                 HHTYPE,
                 STATEFIP:PUMA,
                 GQ:MORTGAGE,
                 TAXINCL:INSINCL,
                 FOODSTMP,
                 SSMC:MULTGEND,
                 FAMUNIT:SCHOOL,
                 EDUC:LABFORCE,
                 WKSWORK2:WORKEDYR,
                 MIGRATE1:MIGPUMA1,
                 DIFFREM:GCRESPON), 
            to_factor) %>% 
  mutate(OCC2010.fac = to_factor(OCC2010), # to have factor and numeric versions of occ and ind codes
         OCCSOC.fac = to_factor(OCCSOC),
         IND.fac = to_factor(IND)) #%>% 
  #mutate_at(vars(INCTOT:INCBUS00), to_character) # this takes the labels for the income variables but makes them character

# clear the original object to save memory.
rm(ipums.la)

```



For convenience, define PUMA codes for querying by geography. Maps with the exact definitions for 2010 PUMAs are on the [IPUMS website](https://usa.ipums.org/usa/volii/pumas10.shtml).

```{r}
## Create filters for PUMA codes -- these probably changed in 1990 and 2000.
# ipums provides the best documentation. The fed websites are becoming unusable. 

super.region.pumacodes.2010 <- c(
  "1400",	# Capital Region Planning Commission 1--Baton Rouge Metropolitan Area (West) PUMA
  "1500",	# Capital Region Planning Commission 2--East Baton Rouge Parish (South) PUMA
  "1501",	# Capital Region Planning Commission 3--Baton Rouge City PUMA
  "1502",	# Capital Region Planning Commission 4--East Baton Rouge Parish (North) PUMA
  "1600",	# Capital Region Planning Commission 5--Ascension Parish PUMA
  "1700",	# Capital Region Planning Commission 6--Livingston & St. Helena Parishes PUMA
  "1800",	# Capital Region Planning Commission 7--Tangipahoa & Washington Parishes PUMA
  "1900",	# South Central Planning & Development Commission 1--River Parishes PUMA
  "2000",	# South Central Planning & Development Commission 2--Lafourche & Assumption Parishes PUMA
  "2100",	# South Central Planning & Development Commission 3--Terrebonne Parish PUMA
  "2200",	# Regional Planning Commission 1--St. Tammany Parish (Northwest) PUMA
  "2201",	# Regional Planning Commission 2--St. Tammany Parish (Southeast) PUMA
  "2300",	# Regional Planning Commission 3--Jefferson Parish (North) PUMA
  "2301",	# Regional Planning Commission 4--Jefferson Parish (Central) PUMA
  "2302",	# Regional Planning Commission 5--Jefferson Parish (West Bank) PUMA
  "2400",	# Regional Planning Commission 6--New Orleans City (Northeast) PUMA
  "2401",	# Regional Planning Commission 7--New Orleans City (Central) PUMA
  "2402",	# Regional Planning Commission 8--New Orleans City (South) PUMA
  "2500"	# Regional Planning Commission 9--St. Bernard, Jefferson (South) & Plaquemines Parishes PUMA
)
  
orleans.pumacodes <- c("2400","2401","2402")

# This should be double checked.
NOmsa.pumacodes <- c("2400","2401","2402", # Orleans
                     "2500", # St Bernard, Plaquemines, south Jefferson
                     "2300","2301","2302", # Jefferson
                     "2200","2201", # St Tammany
                     "1900") # river parishes (Includes St. James, I think)

```



## Examples

Here are some examples using srvyr. `as_survey_design` and `survey_mean` are from the [srvyr package](https://cran.r-project.org/web/packages/srvyr/srvyr.pdf). The command `as_survey_design` creates a tbl_svy survey object. 

[This link](https://usa.ipums.org/usa/complex_survey_vars/userNotes_variance.shtml) describes the STRATA and CLUSTER variables in IPUMS. Here is the important part:

> Three Technical Variables Are Needed for Analysis of the IPUMS Data:

> 1. A sampling weight must be chosen (PERWT, HHWT, or SLWT), based on the type of question a researcher is trying to answer. Analysts should review variable descriptions of the variables of interest and the "Sampling Weights" note for more information about which weight to use.
2. STRATA is an integrated variable that represents the impact of the sample design stratification on the estimates of variance and standard errors.
3. CLUSTER is an integrated variable which uniquely identifies each household record in a given sample.

As an alternative, replicate weights can be used as described [here](https://usa.ipums.org/usa/repwt.shtml). The current extract does not include the replicate weight variables, but these can easily be added, and the command `srvyr::as_survey_rep()` could be used to prep the data frame. 



```{r}
## Median income
# Find codes for employment status
ipums.la.fac %>% count(EMPSTAT)

# Total income for employed workers.
ipums.la.fac %>% 
  filter(PUMA %in% orleans.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  summarise(mean.income = survey_mean(INCTOT, vartype = "se"))

# vartype can be "se", "ci", "var", or "cv".

## Race
# Estimate race breakdown using srvyr
ipums.la.fac %>% 
  filter(PUMA %in% orleans.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  survey_count(RACE, vartype = "se") %>% 
  mutate(pct = n/sum(n))

# Same thing with normal dplyr functions -- same results, no variance.
ipums.la.fac %>% 
  filter(PUMA %in% orleans.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  count(RACE, wt=PERWT) %>% 
  mutate(pct = n/sum(n))

## Income by race
ipums.la.fac %>% 
  filter(PUMA %in% orleans.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  group_by(RACE) %>% 
  summarise(mean.income = survey_mean(INCTOT, vartype = "se"))

```

### Additional concerns

- Use the appropriate weight: PERWT for individuals and HHWT for households. 
- The difference between household and individual data is very important. I'm not exactly sure how this works with CLUSTER, but you wouldn't want to count a household with 5 people 5 times when calculating the median household income. 
- You might need a weighted median/weighted quantile function when calculating medians (e.g., Hmisc). survey_median and survey_quantile are probably preferable. 
- When you use weights on a five year sample, you will get a five year estimate, so you will have to divide by 5. If you only want a median income or a percentage, you can ignore this issue. 

