# ipums EDA

# load libraries
library(tidyverse)
library(ipumsr)


# Step 1: Read the IPUMS extract
ddi <- read_ipums_ddi("../data/IPUMS/usa_00003.xml")
ipums <- read_ipums_micro(ddi)

ipums <- ipums |>
    janitor::clean_names() |>
    # Only convert genuinely categorical columns to factors; leave
    # numeric variables (inctot, uhrswork, trantime, rentgrs, age) as
    # labelled numerics so arithmetic (quantile, mean, etc.) still works
    mutate(across(
        c(degfield, tranwork, sex, race),
        haven::as_factor
    )) |>
    select(-gq, strata, cluster, cbserial, sample)

# load and save smaller versions of the data
# for each year of interest, load the data and save it as a separate RDS file

years <- c(2006, 2016, 2024)
for (year in years) {
    ipums_year <- filter(ipums, year == !!year)
    saveRDS(ipums_year, paste0("../data/IPUMS/ipums_", year, ".rds"))
}

# remove NA from inctot and degfield
# don't remove NA from degfield, all of 2006 is that way!
ipums <- ipums |> filter(!is.na(inctot))

# filter out the top 1% of inctot to avoid extreme outliers
ipums <- ipums |> filter(inctot < quantile(inctot, 0.99))


# Step 2: Clean and filter the sample
df <- ipums |>
    filter(
        degfield != "N/A", # only include those with a degree field
        age >= 25 & age <= 64, # prime working age
        inctot > 0 & inctot < 9999998, # valid non-zero earnings
        uhrswork >= 35 # full-time workers
    ) |>
    glimpse()

df |>
    filter(inctot < 150000) |>
    ggplot(aes(
        y = reorder(degfield, inctot),
        x = inctot,
        group = as_factor(degfield)
    )) +
    geom_boxplot()


# use group_by * count as part of EDA
group_by(ipums, race) |> summarize(n = n())

# is equal to:
count(ipums, race)

group_by(ipums, race) |>
    summarize(
        n = n(),
        mean_income = mean(inctot, na.rm = TRUE),
        median_income = median(inctot, na.rm = TRUE)
    ) |>
    arrange(desc(mean_income))