What Predicts Household Income Across Mid-Atlantic Counties?

Quarto

Quarto enables you to weave together content and executable code into a finished document. To learn more about Quarto see https://quarto.org.

Running Code

When you click the Render button a document will be generated that includes both content and the output of embedded code. You can embed code like this:

1 + 1
[1] 2

You can add options to executable code like this

[1] 4

The echo: false option disables the printing of code (only output is displayed).

[![Map of U.S. counties shaded by per capita income; darker areas have higher income. Image: Jajhill, “United States Counties Per Capita Income,” Wikimedia Commons, licensed CC BY-SA 3.0 (https://creativecommons.org/licenses/by-sa/3.0/). Source: https://commons.wikimedia.org/wiki/File:United_States_Counties_Per_Capita_Income.png](https://commons.wikimedia.org/wiki/Special:FilePath/United_States_Counties_Per_Capita_Income.png){width=80%}](Map of U.S. counties shaded by per capita income; darker areas have higher income. Image: Jajhill, “United States Counties Per Capita Income,” Wikimedia Commons, licensed CC BY-SA 3.0 (https://creativecommons.org/licenses/by-sa/3.0/). Source: https://commons.wikimedia.org/wiki/File:United_States_Counties_Per_Capita_Income.png)

Introduction

Household income is one of the most common ways that governments and economists compare communities. It varies widely, even between neighboring counties and neighborhoods, and this project explores why. I look at what is associated with median household income across counties in seven Mid-Atlantic jurisdictions: Maryland, Virginia, the District of Columbia, Delaware, Pennsylvania, West Virginia and New Jersey.

I plan to do four things. First, I explore the distribution of income and how it differs by state. Second, I use multiple linear regression to test which county characteristics are associated with income: education, home values, commute time, homeownership rate and population density. Third, I use diagnostic plots to check whether the regression model meets its assumptions. Fourth, I present the results with an interactive chart, a static chart and a county map, which make the patterns easier to see.

Where the data came from. The data is the County Demographics dataset from the CORGIS Dataset Project (Bart et al., version 3.0.0, 2021), which assembled it from the U.S. Census Bureau’s QuickFacts tables. The original source is the U.S. Census Bureau (https://www.census.gov/quickfacts). The Census Bureau collects these figures through surveys and censuses. The income, education and housing estimates come mainly from the American Community Survey for 2015 to 2019, an ongoing survey of a sample of U.S. households. The population counts come from the decennial census, which is taken every 10 years, and the business figures come from the 2012 Economic Census. Because the American Community Survey is a sample survey, each county value is an estimate with a margin of error, so I interpret small differences between counties with caution.

The variables. Each row is one county. The full dataset has 43 variables: 2 categorical variables (County and State) and 41 quantitative variables. I use County and State as categorical variables, because income and home values can be very different from one state to another. I use six quantitative variables: median household income (the response I want to predict), the percent of adults with a bachelor’s degree or higher, median home value, mean travel time to work, homeownership rate, and population per square mile.

How I cleaned the data. The source codes missing values as -1, so I recoded every -1 as NA to keep the analysis accurate. I selected and renamed the columns I need, then used dplyr’s filter() to keep only the seven jurisdictions and only counties with no missing values on my variables, which brought the data well under 800 observations. Population density is strongly right-skewed, so I used mutate() to create a new variable, log_density, a base-10 log of density.

Why I chose this topic. I wanted to see how demographic and housing factors like income, education and homeownership relate to each other. The seven jurisdictions I chose are very diverse, from dense urban areas like Washington, D.C. to rural counties in West Virginia and Pennsylvania, so they show a wide range of conditions in one region.

library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.3     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(RColorBrewer)
library(highcharter)
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 
library(maps)
Warning: package 'maps' was built under R version 4.6.2

Attaching package: 'maps'

The following object is masked from 'package:purrr':

    map
library(scales)

Attaching package: 'scales'

The following object is masked from 'package:purrr':

    discard

The following object is masked from 'package:readr':

    col_factor
library(ggfortify)
library(DataExplorer)
Warning: package 'DataExplorer' was built under R version 4.6.2
downloads <- file.path(Sys.getenv("USERPROFILE"), "Downloads")
path <- file.path(downloads, "county_demographics.csv")
if (!file.exists(path)) path <- paste0(path, ".csv")
counties_raw <- read_csv(path)
Rows: 3139 Columns: 43
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): County, State
dbl (41): Age.Percent 65 and Older, Age.Percent Under 18 Years, Age.Percent ...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

This chunk loads the packages I need and reads the county data from my Downloads folder with read_csv.

states_keep <- c("MD", "VA", "DC", "DE", "PA", "WV", "NJ")

counties <- counties_raw |>
  mutate(across(where(is.numeric), ~ na_if(.x, -1))) |>
  select(county = County,
         state = State,
         income = `Income.Median Houseold Income`,
         bachelors = `Education.Bachelor's Degree or Higher`,
         home_value = `Housing.Median Value of Owner-Occupied Units`,
         commute = `Miscellaneous.Mean Travel Time to Work`,
         homeownership = `Housing.Homeownership Rate`,
         density = `Population.Population per Square Mile`) |>
  filter(state %in% states_keep,
         !is.na(income), !is.na(bachelors), !is.na(home_value),
         !is.na(commute), !is.na(homeownership), !is.na(density)) |>
  mutate(log_density = log10(density))

nrow(counties)
[1] 304

Here I cleaned the data. I changed the source’s -1 missing codes to NA, kept and renamed the columns I need, and used filter() to keep seven jurisdictions with complete data. That left 304 counties, under the 800 limit. I also created log_density because population density is very skewed.

counties |>
  count(state) |>
  ggplot(aes(x = reorder(state, n), y = n, fill = state)) +
  geom_col(show.legend = FALSE) +
  scale_fill_brewer(palette = "Set2") +
  coord_flip() +
  labs(title = "Number of Counties per State in the Sample",
       x = "State", y = "Number of counties",
       caption = "Source: CORGIS County Demographics (U.S. Census Bureau)") +
  theme_minimal(base_size = 13)

ggplot(counties, aes(x = income)) +
  geom_histogram(bins = 25, fill = "steelblue", color = "white") +
  scale_x_continuous(labels = dollar) +
  labs(title = "Distribution of Median Household Income",
       x = "Median household income", y = "Number of counties",
       caption = "Source: CORGIS County Demographics (U.S. Census Bureau)") +
  theme_minimal(base_size = 13)

ggplot(counties, aes(x = reorder(state, income, median), y = income, fill = state)) +
  geom_boxplot(show.legend = FALSE) +
  scale_fill_brewer(palette = "Set2") +
  scale_y_continuous(labels = dollar) +
  labs(title = "Median Household Income by State",
       x = "State", y = "Median household income",
       caption = "Source: CORGIS County Demographics (U.S. Census Bureau)") +
  theme_minimal(base_size = 13)

These three plots are my first look at the data: a bar chart of how many counties each state has, a histogram of income, and boxplots comparing income by state. [One sentence on what you notice, for example which state has the highest typical income.]

counties_num <- counties |>
  select(income, bachelors, home_value, commute, homeownership, log_density)

plot_correlation(counties_num)

fit_full <- lm(income ~ bachelors + home_value + commute + homeownership + log_density,
               data = counties)

final_model <- step(fit_full, direction = "backward", trace = 0)
summary(final_model)

Call:
lm(formula = income ~ bachelors + home_value + commute + homeownership + 
    log_density, data = counties)

Residuals:
     Min       1Q   Median       3Q      Max 
-22124.4  -3703.6    832.2   3838.4  22061.4 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -5.624e+04  4.972e+03 -11.312  < 2e-16 ***
bachelors      4.516e+02  7.923e+01   5.700 2.88e-08 ***
home_value     1.062e-01  8.895e-03  11.944  < 2e-16 ***
commute        3.014e+02  8.365e+01   3.603 0.000368 ***
homeownership  8.219e+02  5.413e+01  15.185  < 2e-16 ***
log_density    7.296e+03  9.355e+02   7.799 1.05e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6340 on 298 degrees of freedom
Multiple R-squared:  0.9062,    Adjusted R-squared:  0.9046 
F-statistic: 575.9 on 5 and 298 DF,  p-value: < 2.2e-16

The correlation plot shows which variables move together with income. I fit a model with all five predictors, then used backward elimination to remove the weakest ones and keep a simpler model. [One sentence on which predictor has the strongest link to income.]

autoplot(final_model, 1:4, nrow = 2, ncol = 2)

coefs <- coef(final_model)
equation <- paste0("income = ", round(coefs[1], 1),
                   paste0(sprintf(" %+.2f x %s", coefs[-1], names(coefs)[-1]),
                          collapse = ""))

model_summary <- summary(final_model)
adj_r2 <- model_summary$adj.r.squared
f <- model_summary$fstatistic
model_p <- pf(f[1], f[2], f[3], lower.tail = FALSE)
coef_table <- as.data.frame(model_summary$coefficients)
sig_vars <- rownames(coef_table)[coef_table[, 4] < 0.05 & rownames(coef_table) != "(Intercept)"]

equation
[1] "income = -56240.8 +451.64 x bachelors +0.11 x home_value +301.43 x commute +821.92 x homeownership +7296.39 x log_density"

My final model is: income = -56240.8 +451.64 x bachelors +0.11 x home_value +301.43 x commute +821.92 x homeownership +7296.39 x log_density

The adjusted R-squared is 0.905, so the model explains about 90.5% of the variation in income. The model p-value is <2e-16, and the predictors below 0.05, which are the significant ones, are bachelors, home_value, commute, homeownership, log_density. The diagnostic plots check the model: Residuals vs. Fitted checks that a line fits, Q-Q checks that the errors are roughly normal, Scale-Location checks for even spread, and Cook’s Distance flags influential counties. [One sentence on what yours show.]

cols <- brewer.pal(7, "Set2")

hchart(counties, "scatter",
       hcaes(x = bachelors, y = income, group = state, name = county)) |>
  hc_colors(cols) |>
  hc_title(text = "Counties With More College Graduates Have Higher Household Incomes") |>
  hc_xAxis(title = list(text = "Adults with a bachelor's degree or higher (%)")) |>
  hc_yAxis(title = list(text = "Median household income ($)")) |>
  hc_legend(align = "right", verticalAlign = "top", layout = "vertical") |>
  hc_caption(text = "Source: CORGIS County Demographics (U.S. Census Bureau)") |>
  hc_tooltip(pointFormat = "<b>{point.name}</b><br>Education: {point.x}%<br>Income: ${point.y:,.0f}")

This interactive chart shows that counties with more college graduates tend to have higher incomes. Hover over a point to see the county, its education rate and its income.

ggplot(counties, aes(x = home_value, y = income, color = state, size = commute)) +
  geom_point(alpha = 0.7) +
  geom_smooth(aes(group = 1), method = "lm", se = FALSE,
              color = "black", linetype = 2, linewidth = 0.5, show.legend = FALSE) +
  annotate("text", x = max(counties$home_value) * 0.75, y = min(counties$income) * 1.1,
           label = "Dashed line = overall linear trend", size = 3.5) +
  scale_color_brewer(palette = "Set2", name = "State") +
  scale_size_continuous(range = c(1.5, 8), name = "Mean commute (minutes)") +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = dollar) +
  labs(title = "Higher Home Values Go Along With Higher Household Incomes",
       subtitle = "Mid-Atlantic counties; point size shows mean travel time to work",
       x = "Median value of owner-occupied homes",
       y = "Median household income",
       caption = "Source: CORGIS County Demographics (U.S. Census Bureau)") +
  theme_bw(base_size = 13)
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
`geom_smooth()` using formula = 'y ~ x'
Warning: The following aesthetics were dropped during statistical transformation: size.
ℹ This can happen when ggplot fails to infer the correct grouping structure in
  the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?

This chart shows that counties with higher home values also tend to have higher incomes. Color shows the state, point size shows commute time, and the dashed line is the overall trend.

state_lookup <- tibble(state = c(state.abb, "DC"),
                       state_name = c(tolower(state.name), "district of columbia"))

income_keys <- counties |>
  left_join(state_lookup, by = "state") |>
  mutate(key = paste(state_name,
                     county |> tolower() |> str_remove(" county$") |> str_remove_all("[.']"),
                     sep = "|"))

county_map <- map_data("county") |>
  filter(region %in% state_lookup$state_name[state_lookup$state %in% states_keep]) |>
  mutate(key = paste(region, subregion, sep = "|")) |>
  left_join(select(income_keys, key, income), by = "key")

ggplot(county_map, aes(x = long, y = lat, group = group, fill = income)) +
  geom_polygon(color = "white", linewidth = 0.1) +
  annotate("text", x = -77.04, y = 38.55, label = "Washington, D.C. area", size = 3.5) +
  scale_fill_viridis_c(option = "plasma", labels = dollar,
                       name = "Median household\nincome", na.value = "grey85") +
  coord_quickmap() +
  labs(title = "Median Household Income by County Across the Mid-Atlantic",
       subtitle = "Grey areas had no match in the map data",
       caption = "Source: CORGIS County Demographics (U.S. Census Bureau); map: R maps package") +
  theme_void(base_size = 13)

The map shows median household income for each county, with darker and lighter colors showing high and low income. Grey counties did not match the map data. [One sentence on a pattern, for example whether the Washington, D.C. area stands out.]

Conclusion

In conclusion, [the strongest predictor] has the strongest link to county income, and my model explains about 90.5% of the variation. [One limit from the diagnostic plots.] Next I would add more states or more predictors.

Sources and Bibliography

I based my code on our class notes, and I also checked the online documentation for ggplot2, the step function, ggfortify and highcharter to make sure I was using the functions correctly.