Code
# load libraries
library(tidyverse)
library(tidycensus)
library(readr)
library(modelsummary)
library(fixest)
# We also set the census API key but keep that hidden here.Estimated Food Insecurity in Western North Carolina
Technical Report
Introduction
Nourish WNC, is an arm of The Valley Hope Foundation that aims to provide nutritious food to individuals and families facing food insecurity in Western North Carolina. Nourish WNC is also especially concerned with reaching underserved rural communities.
To enhance Nourish WNC’s ability to locate and serve communities most in need we estimated food insecurity rates and population density in census tracts in Western NC. A visualization of the results of our estimation can be found on the Valley Hope Foundation website.
Below is a detailed explanation how we generated our food insecurity estimates.
We used the R programming language to generate our estimates and to access some of the necessary data. Examples of our R code are provided throughout this document.
Geography
To begin we downloaded the NC places, NC primary and secondary roads, NC census tracts and US counties shapefiles from the US Census Bureau website. We limited the data to just the 26 counties relevant to Nourish WNC: Alexander, Alleghany, Ashe, Avery, Buncombe, Burke, Caldwell, Catawba, Cherokee, Clay, Cleveland, Graham, Haywood, Henderson, Jackson, Macon, Madison, McDowell, Mitchell, Polk, Rutherford, Swain, Transylvania, Watauga, Wilkes and Yancey.
Population Density
We also downloaded each census tract’s total population estimate from the American Community Survey using the tidycensus package in R. The 2024 data is the most recent available. We then divided each census tract’s total population (variable B01003_001) by its land area (reported in the shapefiles) to find its population density.
Code
# Pull 2025 ACS 5-year total population for NC census tracts
nc_pop <- get_acs(
geography = "tract",
state = "NC",
variables = "B01003_001",
year = 2024,
survey = "acs5",
geometry = FALSE,
cache_table = TRUE
)
# Export data for use in QGIS or elsewhere
write_csv(nc_pop, "Data/nc_pop.csv")
Estimating Food Insecurity
To estimate food insecurity rates for census tracts in Western NC we used a modified version of the method employed by Feeding America’s Map the Meal Gap study. Map the Meal Gap estimates food insecurity rates by US county. An overview of the study’s methodology can be found here.
Our method relies on using the food insecurity rates per state (provided by the USDA) to fit a regression model that estimates the extent to which population-level attributes predict food insecurity rates in a community. The coefficients derived from that model are then applied to data at the census tract level to estimate food insecurity rates for census tracts in Western NC.
Collecting food insecurity rates per state
To find food insecurity rates per state we consulted the Household Food Insecurity in the United States 2024 report, accessible on the USDA Economic Research Service website. This gave us the average food insecurity rate per state from the 2022 to 2024 samples and the 2019 to 2021 samples. We got the 2016 to 2018 average and the 2013 to 2015 average from the Household Food Insecurity in the United States in 2018 report found here.
Code
# import the USDA data stored on a local drive:
Food_Insecurity_Rate_by_State_csv <-
read_csv("Data/Food Insecurity Rate by State from HFSinUS2024and2018 - Data.csv")
# format data
FI_Rates_by_State <- Food_Insecurity_Rate_by_State_csv |>
pivot_longer(cols = starts_with("20"),
names_to = "Time_Period",
values_to = "FI_Rate") |>
mutate("observation" = paste(State, Time_Period, sep = "-"),
"FI_Rate" = FI_Rate * .01)Households are considered food insecure if they report three or more “food-insecure conditions” during the previous 12 months.1 Some examples of food insecure conditions are worrying about running out of food before being able to buy more, eating less than one should because there is not enough money for food, or failing to eat for a whole day due to not having enough money for food. A detailed treatment of how the data are collected can be found in the Household Food Insecurity in the United States reports.
Getting state-level food insecurity predictors
We then gathered state-level census data for each state using the tidycensus R package to access US Census Bureau APIs. We collected data on the same or similar measures to those used in the Map the Meal Gap analysis. The following data points were collected:2
total_pop= Total population from the ACS Demographic and Housing Estimates table (tableDP05).black_alone= Population identifying as Black alone from the ACS Demographic and Housing Estimates table (tableDP05).hispanic= Population identifying as Hispanic or Latino from the ACS Demographic and Housing Estimates table (tableDP05).total_diability= Total non-institutionalized population from the Disability Characteristics table (tableS1810).with_disability= Total non-institutionalized population with a disability from the Disability Characteristics table (tableS1810).occupied_housing= Number of housing units that are occupied from the Selected Housing Characteristics table (tableDP04).owner_occupied= Number of housing units that are owner occupied from the Selected Housing Characteristics table (tableDP04).median_income= Median household income from the Median Household Income in the Past 12 Months table (tableB19013).poverty_total= Total population from the Poverty Status in the Past 12 Months by School Enrollment by Level of School for the Population 3 Years and Over table (tableB14006).poverty_below= Population with income in the past 12 months below the poverty level from the Poverty Status in the Past 12 Months by School Enrollment by Level of School for the Population 3 Years and Over table (tableB14006).unemployment_rate= Unemployment rate from the Employment Status table (tableS2301).
Code
# Get ACS data for all states in 2024
# check each variable for each year: # valid for:
vars2024 <- c(
total_pop = "DP05_0001", # 2024
black_alone = "DP05_0045", # 2024
hispanic = "DP05_0090", # 2024
total_disability = "S1810_C01_001", # 2024
with_disability = "S1810_C02_001", # 2024
occupied_housing = "DP04_0002", # 2024
owner_occupied = "DP04_0046", # 2024
median_income = "B19013_001", # 2024
poverty_total = "B14006_001", # 2024
poverty_below = "B14006_002", # 2024
unemployment_rate = "S2301_C04_001" # 2024
)
State_FI_Predictors_2024 <- get_acs(
geography = "state",
variables = vars2024,
year = 2024,
survey = "acs1",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
State_FI_Predictors_2024 <- State_FI_Predictors_2024 %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
State_FI_Predictors_2024 <- State_FI_Predictors_2024 %>%
mutate(
disability_rate_2024 = with_disability / total_disability,
homeownership_rate_2024 = owner_occupied / occupied_housing,
percent_black_2024 = black_alone / total_pop,
percent_hispanic_2024 = hispanic / total_pop,
poverty_rate_2024 = poverty_below / poverty_total,
unemployment_rate_2024 = unemployment_rate * 0.01,
median_income_2024 = median_income,
total_pop_2024 = total_pop
)
# Get ACS data for all states in 2023
# check each variable for each year: # valid for:
vars2023 <- c(
total_pop = "DP05_0001", # 2023
black_alone = "DP05_0038", # 2023
hispanic = "DP05_0076", # 2023
total_disability = "S1810_C01_001", # 2023
with_disability = "S1810_C02_001", # 2023
occupied_housing = "DP04_0002", # 2023
owner_occupied = "DP04_0046", # 2023
median_income = "B19013_001", # 2023
poverty_total = "B14006_001", # 2023
poverty_below = "B14006_002", # 2023
unemployment_rate = "S2301_C04_001" # 2023
)
State_FI_Predictors_2023 <- get_acs(
geography = "state",
variables = vars2023,
year = 2023,
survey = "acs1",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
State_FI_Predictors_2023 <- State_FI_Predictors_2023 %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
State_FI_Predictors_2023 <- State_FI_Predictors_2023 %>%
mutate(
disability_rate_2023 = with_disability / total_disability,
homeownership_rate_2023 = owner_occupied / occupied_housing,
percent_black_2023 = black_alone / total_pop,
percent_hispanic_2023 = hispanic / total_pop,
poverty_rate_2023 = poverty_below / poverty_total,
unemployment_rate_2023 = unemployment_rate * 0.01,
median_income_2023 = median_income,
total_pop_2023 = total_pop
)
# Get ACS data for all states in 2022
# check each variable for each year: # valid for:
vars2022 <- c(
total_pop = "DP05_0001", # 2022
black_alone = "DP05_0038", # 2022
hispanic = "DP05_0073", # 2022
total_disability = "S1810_C01_001", # 2022
with_disability = "S1810_C02_001", # 2022
occupied_housing = "DP04_0002", # 2022
owner_occupied = "DP04_0046", # 2022
median_income = "B19013_001", # 2022
poverty_total = "B14006_001", # 2022
poverty_below = "B14006_002", # 2022
unemployment_rate = "S2301_C04_001" # 2022
)
State_FI_Predictors_2022 <- get_acs(
geography = "state",
variables = vars2022,
year = 2022,
survey = "acs1",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
State_FI_Predictors_2022 <- State_FI_Predictors_2022 %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
State_FI_Predictors_2022 <- State_FI_Predictors_2022 %>%
mutate(
disability_rate_2022 = with_disability / total_disability,
homeownership_rate_2022 = owner_occupied / occupied_housing,
percent_black_2022 = black_alone / total_pop,
percent_hispanic_2022 = hispanic / total_pop,
poverty_rate_2022 = poverty_below / poverty_total,
unemployment_rate_2022 = unemployment_rate * 0.01,
median_income_2022 = median_income,
total_pop_2022 = total_pop
)
# Get ACS data for all states in 2021
# check each variable for each year: # valid for:
vars2017_to_2021 <- c(
total_pop = "DP05_0001", # 2021, 2020, 2019, 2018, 2017
black_alone = "DP05_0038", # 2021, 2020, 2019, 2018, 2017
hispanic = "DP05_0071", # 2021, 2020, 2019, 2018, 2017
total_disability = "S1810_C01_001", # 2021, 2020, 2019, 2018, 2017
with_disability = "S1810_C02_001", # 2021, 2020, 2019, 2018, 2017
occupied_housing = "DP04_0002", # 2021, 2020, 2019, 2018, 2017
owner_occupied = "DP04_0046", # 2021, 2020, 2019, 2018, 2017
median_income = "B19013_001", # 2021, 2020, 2019, 2018, 2017
poverty_total = "B14006_001", # 2021, 2020, 2019, 2018, 2017
poverty_below = "B14006_002", # 2021, 2020, 2019, 2018, 2017
unemployment_rate = "S2301_C04_001" # 2021, 2020, 2019, 2018, 2017
)
State_FI_Predictors_2021 <- get_acs(
geography = "state",
variables = vars2017_to_2021,
year = 2021,
survey = "acs1",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
State_FI_Predictors_2021 <- State_FI_Predictors_2021 %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
State_FI_Predictors_2021 <- State_FI_Predictors_2021 %>%
mutate(
disability_rate_2021 = with_disability / total_disability,
homeownership_rate_2021 = owner_occupied / occupied_housing,
percent_black_2021 = black_alone / total_pop,
percent_hispanic_2021 = hispanic / total_pop,
poverty_rate_2021 = poverty_below / poverty_total,
unemployment_rate_2021 = unemployment_rate * 0.01,
median_income_2021 = median_income,
total_pop_2021 = total_pop
)
# Get ACS data for all states in 2020
State_FI_Predictors_2020 <- get_acs(
geography = "state",
variables = vars2017_to_2021,
year = 2020,
survey = "acs5",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
State_FI_Predictors_2020 <- State_FI_Predictors_2020 %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
State_FI_Predictors_2020 <- State_FI_Predictors_2020 %>%
mutate(
disability_rate_2020 = with_disability / total_disability,
homeownership_rate_2020 = owner_occupied / occupied_housing,
percent_black_2020 = black_alone / total_pop,
percent_hispanic_2020 = hispanic / total_pop,
poverty_rate_2020 = poverty_below / poverty_total,
unemployment_rate_2020 = unemployment_rate * 0.01,
median_income_2020 = median_income,
total_pop_2020 = total_pop
)
# Get ACS data for all states in 2019
State_FI_Predictors_2019 <- get_acs(
geography = "state",
variables = vars2017_to_2021,
year = 2019,
survey = "acs1",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
State_FI_Predictors_2019 <- State_FI_Predictors_2019 %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
State_FI_Predictors_2019 <- State_FI_Predictors_2019 %>%
mutate(
disability_rate_2019 = with_disability / total_disability,
homeownership_rate_2019 = owner_occupied / occupied_housing,
percent_black_2019 = black_alone / total_pop,
percent_hispanic_2019 = hispanic / total_pop,
poverty_rate_2019 = poverty_below / poverty_total,
unemployment_rate_2019 = unemployment_rate * 0.01,
median_income_2019 = median_income,
total_pop_2019 = total_pop
)
# Get ACS data for all states in 2018
State_FI_Predictors_2018 <- get_acs(
geography = "state",
variables = vars2017_to_2021,
year = 2018,
survey = "acs1",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
State_FI_Predictors_2018 <- State_FI_Predictors_2018 %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
State_FI_Predictors_2018 <- State_FI_Predictors_2018 %>%
mutate(
disability_rate_2018 = with_disability / total_disability,
homeownership_rate_2018 = owner_occupied / occupied_housing,
percent_black_2018 = black_alone / total_pop,
percent_hispanic_2018 = hispanic / total_pop,
poverty_rate_2018 = poverty_below / poverty_total,
unemployment_rate_2018 = unemployment_rate * 0.01,
median_income_2018 = median_income,
total_pop_2018 = total_pop
)
# Get ACS data for all states in 2017
State_FI_Predictors_2017 <- get_acs(
geography = "state",
variables = vars2017_to_2021,
year = 2017,
survey = "acs1",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
State_FI_Predictors_2017 <- State_FI_Predictors_2017 %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
State_FI_Predictors_2017 <- State_FI_Predictors_2017 %>%
mutate(
disability_rate_2017 = with_disability / total_disability,
homeownership_rate_2017 = owner_occupied / occupied_housing,
percent_black_2017 = black_alone / total_pop,
percent_hispanic_2017 = hispanic / total_pop,
poverty_rate_2017 = poverty_below / poverty_total,
unemployment_rate_2017 = unemployment_rate * 0.01,
median_income_2017 = median_income,
total_pop_2017 = total_pop
)
# Get ACS data for all states in 2016
# check each variable for each year: # valid for:
vars2015_to_2016 <- c(
total_pop = "DP05_0001", # 2016, 2015
black_alone = "DP05_0033", # 2016, 2015
hispanic = "DP05_0066", # 2016, 2015
total_disability = "S1810_C01_001", # 2016, 2015
with_disability = "S1810_C02_001", # 2016, 2015
occupied_housing = "DP04_0002", # 2016, 2015
owner_occupied = "DP04_0046", # 2016, 2015
median_income = "B19013_001", # 2016, 2015
poverty_total = "B14006_001", # 2016, 2015
poverty_below = "B14006_002", # 2016, 2015
unemployment_rate = "S2301_C04_001" # 2016, 2015
)
State_FI_Predictors_2016 <- get_acs(
geography = "state",
variables = vars2015_to_2016,
year = 2016,
survey = "acs1",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
State_FI_Predictors_2016 <- State_FI_Predictors_2016 %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
State_FI_Predictors_2016 <- State_FI_Predictors_2016 %>%
mutate(
disability_rate_2016 = with_disability / total_disability,
homeownership_rate_2016 = owner_occupied / occupied_housing,
percent_black_2016 = black_alone / total_pop,
percent_hispanic_2016 = hispanic / total_pop,
poverty_rate_2016 = poverty_below / poverty_total,
unemployment_rate_2016 = unemployment_rate * 0.01,
median_income_2016 = median_income,
total_pop_2016 = total_pop
)
# Get ACS data for all states in 2015
State_FI_Predictors_2015 <- get_acs(
geography = "state",
variables = vars2015_to_2016,
year = 2015,
survey = "acs1",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
State_FI_Predictors_2015 <- State_FI_Predictors_2015 %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
State_FI_Predictors_2015 <- State_FI_Predictors_2015 %>%
mutate(
disability_rate_2015 = with_disability / total_disability,
homeownership_rate_2015 = owner_occupied / occupied_housing,
percent_black_2015 = black_alone / total_pop,
percent_hispanic_2015 = hispanic / total_pop,
poverty_rate_2015 = poverty_below / poverty_total,
unemployment_rate_2015 = unemployment_rate * 0.01,
median_income_2015 = median_income,
total_pop_2015 = total_pop
)
# Get ACS data for all states in 2014
# check each variable for each year: # valid for:
vars <- c(
total_pop = "DP05_0001", # 2014, 2013
black_alone = "DP05_0033", # 2014, 2013
hispanic = "DP05_0066", # 2014, 2013
total_disability = "S1810_C01_001", # 2014, 2013
with_disability = "S1810_C02_001", # 2014, 2013
occupied_housing = "DP04_0002", # 2014, 2013
owner_occupied = "DP04_0045", # 2014, 2013
median_income = "B19013_001", # 2014, 2013
poverty_total = "B14006_001", # 2014, 2013
poverty_below = "B14006_002", # 2014, 2013
unemployment_rate = "S2301_C04_001" # 2014, 2013
)
State_FI_Predictors_2014 <- get_acs(
geography = "state",
variables = vars,
year = 2014,
survey = "acs1",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
State_FI_Predictors_2014 <- State_FI_Predictors_2014 %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
State_FI_Predictors_2014 <- State_FI_Predictors_2014 %>%
mutate(
disability_rate_2014 = with_disability / total_disability,
homeownership_rate_2014 = owner_occupied / occupied_housing,
percent_black_2014 = black_alone / total_pop,
percent_hispanic_2014 = hispanic / total_pop,
poverty_rate_2014 = poverty_below / poverty_total,
unemployment_rate_2014 = unemployment_rate * 0.01,
median_income_2014 = median_income,
total_pop_2014 = total_pop
)
# Get ACS data for all states in 2013
State_FI_Predictors_2013 <- get_acs(
geography = "state",
variables = vars,
year = 2013,
survey = "acs1",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
State_FI_Predictors_2013 <- State_FI_Predictors_2013 %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
State_FI_Predictors_2013 <- State_FI_Predictors_2013 %>%
mutate(
disability_rate_2013 = with_disability / total_disability,
homeownership_rate_2013 = owner_occupied / occupied_housing,
percent_black_2013 = black_alone / total_pop,
percent_hispanic_2013 = hispanic / total_pop,
poverty_rate_2013 = poverty_below / poverty_total,
unemployment_rate_2013 = unemployment_rate * 0.01,
median_income_2013 = median_income,
total_pop_2013 = total_pop
)Because the food insecurity estimates we have from the USDA are three year averages, we generated corresponding three year averages for the census data we collected.
Code
State_FI_Predictors_2022_to_2024 <- State_FI_Predictors_2024 |>
select(GEOID,
NAME,
total_pop_2024,
disability_rate_2024,
homeownership_rate_2024,
percent_black_2024,
percent_hispanic_2024,
poverty_rate_2024,
unemployment_rate_2024,
median_income_2024) |>
left_join(State_FI_Predictors_2023 |>
select(GEOID,
total_pop_2023,
disability_rate_2023,
homeownership_rate_2023,
percent_black_2023,
percent_hispanic_2023,
poverty_rate_2023,
unemployment_rate_2023,
median_income_2023)) |>
left_join(State_FI_Predictors_2022 |>
select(GEOID,
disability_rate_2022,
total_pop_2022,
homeownership_rate_2022,
percent_black_2022,
percent_hispanic_2022,
poverty_rate_2022,
unemployment_rate_2022,
median_income_2022))
State_FI_Predictors_2022_to_2024 <- State_FI_Predictors_2022_to_2024 |>
mutate(Year_Range = "2022 to 2024",
Avg_pop =
(total_pop_2022 +
total_pop_2023 +
total_pop_2024) / 3,
Avg_disability_rate =
(disability_rate_2022 +
disability_rate_2023 +
disability_rate_2024) / 3,
Avg_homeownership_rate =
(homeownership_rate_2022 +
homeownership_rate_2023 +
homeownership_rate_2024) / 3,
Avg_percent_black =
(percent_black_2022 +
percent_black_2023 +
percent_black_2024) / 3,
Avg_percent_hispanic =
(percent_hispanic_2022 +
percent_hispanic_2023 +
percent_hispanic_2024) / 3,
Avg_poverty_rate =
(poverty_rate_2022 +
poverty_rate_2023 +
poverty_rate_2024) / 3,
Avg_unemployment_rate =
(unemployment_rate_2022 +
unemployment_rate_2023 +
unemployment_rate_2024) / 3,
Avg_median_income =
(median_income_2022 +
median_income_2023 +
median_income_2024) / 3)
State_FI_Predictors_2019_to_2021 <- State_FI_Predictors_2021 |>
select(GEOID,
NAME,
total_pop_2021,
disability_rate_2021,
homeownership_rate_2021,
percent_black_2021,
percent_hispanic_2021,
poverty_rate_2021,
unemployment_rate_2021,
median_income_2021) |>
left_join(State_FI_Predictors_2020 |>
select(GEOID,
total_pop_2020,
disability_rate_2020,
homeownership_rate_2020,
percent_black_2020,
percent_hispanic_2020,
poverty_rate_2020,
unemployment_rate_2020,
median_income_2020)) |>
left_join(State_FI_Predictors_2019 |>
select(GEOID,
total_pop_2019,
disability_rate_2019,
homeownership_rate_2019,
percent_black_2019,
percent_hispanic_2019,
poverty_rate_2019,
unemployment_rate_2019,
median_income_2019))
State_FI_Predictors_2019_to_2021 <- State_FI_Predictors_2019_to_2021 |>
mutate(Year_Range = "2019 to 2021",
Avg_pop =
(total_pop_2019 +
total_pop_2020 +
total_pop_2021) / 3,
Avg_disability_rate =
(disability_rate_2019 +
disability_rate_2020 +
disability_rate_2021) / 3,
Avg_homeownership_rate =
(homeownership_rate_2019 +
homeownership_rate_2020 +
homeownership_rate_2021) / 3,
Avg_percent_black =
(percent_black_2019 +
percent_black_2020 +
percent_black_2021) / 3,
Avg_percent_hispanic =
(percent_hispanic_2019 +
percent_hispanic_2020 +
percent_hispanic_2021) / 3,
Avg_poverty_rate =
(poverty_rate_2019 +
poverty_rate_2020 +
poverty_rate_2021) / 3,
Avg_unemployment_rate =
(unemployment_rate_2019 +
unemployment_rate_2020 +
unemployment_rate_2021) / 3,
Avg_median_income =
(median_income_2019 +
median_income_2020 +
median_income_2021) / 3)
State_FI_Predictors_2016_to_2018 <- State_FI_Predictors_2016 |>
select(GEOID,
NAME,
total_pop_2016,
disability_rate_2016,
homeownership_rate_2016,
percent_black_2016,
percent_hispanic_2016,
poverty_rate_2016,
unemployment_rate_2016,
median_income_2016) |>
left_join(State_FI_Predictors_2017 |>
select(GEOID,
total_pop_2017,
disability_rate_2017,
homeownership_rate_2017,
percent_black_2017,
percent_hispanic_2017,
poverty_rate_2017,
unemployment_rate_2017,
median_income_2017)) |>
left_join(State_FI_Predictors_2018 |>
select(GEOID,
total_pop_2018,
disability_rate_2018,
homeownership_rate_2018,
percent_black_2018,
percent_hispanic_2018,
poverty_rate_2018,
unemployment_rate_2018,
median_income_2018))
State_FI_Predictors_2016_to_2018 <- State_FI_Predictors_2016_to_2018 |>
mutate(Year_Range = "2016 to 2018",
Avg_pop =
(total_pop_2016 +
total_pop_2017 +
total_pop_2018) / 3,
Avg_disability_rate =
(disability_rate_2016 +
disability_rate_2017 +
disability_rate_2018) / 3,
Avg_homeownership_rate =
(homeownership_rate_2016 +
homeownership_rate_2017 +
homeownership_rate_2018) / 3,
Avg_percent_black =
(percent_black_2016 +
percent_black_2017 +
percent_black_2018) / 3,
Avg_percent_hispanic =
(percent_hispanic_2016 +
percent_hispanic_2017 +
percent_hispanic_2018) / 3,
Avg_poverty_rate =
(poverty_rate_2016 +
poverty_rate_2017 +
poverty_rate_2018) / 3,
Avg_unemployment_rate =
(unemployment_rate_2016 +
unemployment_rate_2017 +
unemployment_rate_2018) / 3,
Avg_median_income =
(median_income_2016 +
median_income_2017 +
median_income_2018) / 3)
State_FI_Predictors_2013_to_2015 <- State_FI_Predictors_2013 |>
select(GEOID,
NAME,
total_pop_2013,
disability_rate_2013,
homeownership_rate_2013,
percent_black_2013,
percent_hispanic_2013,
poverty_rate_2013,
unemployment_rate_2013,
median_income_2013) |>
left_join(State_FI_Predictors_2014 |>
select(GEOID,
total_pop_2014,
disability_rate_2014,
homeownership_rate_2014,
percent_black_2014,
percent_hispanic_2014,
poverty_rate_2014,
unemployment_rate_2014,
median_income_2014)) |>
left_join(State_FI_Predictors_2015 |>
select(GEOID,
total_pop_2015,
disability_rate_2015,
homeownership_rate_2015,
percent_black_2015,
percent_hispanic_2015,
poverty_rate_2015,
unemployment_rate_2015,
median_income_2015))
State_FI_Predictors_2013_to_2015 <- State_FI_Predictors_2013_to_2015 |>
mutate(Year_Range = "2013 to 2015",
Avg_pop =
(total_pop_2013 +
total_pop_2014 +
total_pop_2015) / 3,
Avg_disability_rate =
(disability_rate_2013 +
disability_rate_2014 +
disability_rate_2015) / 3,
Avg_homeownership_rate =
(homeownership_rate_2013 +
homeownership_rate_2014 +
homeownership_rate_2014) / 3,
Avg_percent_black =
(percent_black_2013 +
percent_black_2014 +
percent_black_2015) / 3,
Avg_percent_hispanic =
(percent_hispanic_2013 +
percent_hispanic_2014 +
percent_hispanic_2015) / 3,
Avg_poverty_rate =
(poverty_rate_2013 +
poverty_rate_2014 +
poverty_rate_2015) / 3,
Avg_unemployment_rate =
(unemployment_rate_2013 +
unemployment_rate_2014 +
unemployment_rate_2015) / 3,
Avg_median_income =
(median_income_2013 +
median_income_2014 +
median_income_2015) / 3)
State_FI_Predictors <- bind_rows(State_FI_Predictors_2013_to_2015,
State_FI_Predictors_2016_to_2018,
State_FI_Predictors_2019_to_2021,
State_FI_Predictors_2022_to_2024) |>
select(GEOID,
NAME,
Year_Range,
Avg_pop,
Avg_disability_rate,
Avg_homeownership_rate,
Avg_percent_black,
Avg_percent_hispanic,
Avg_poverty_rate,
Avg_unemployment_rate,
Avg_median_income) |>
mutate("observation" = paste(NAME, Year_Range, sep = "-")) |>
filter(NAME != "District of Columbia", NAME != "Puerto Rico")Fitting predictor models
We then used OLS regression to estimate the extent to which each of the variables we collected can be used to predict food insecurity rates at the state level.
We fit three models. The first model uses all seven predictors plus fixed effects for state and time period. It is as follows:
\[ \begin{gathered} FoodInsecurityRate_{st} = \beta_0 + \beta_1 PovertyRate_{st} + \beta_2 UnemploymentRate_{st} \\ + \beta_3 MedianIncome_{st} + \beta_4 PercentHispanic_{st} + \beta_5 PercentBlack_{st} + \\ \beta_6 HomeOwnershipRate_{st} + \beta_7 DisabilityRate_{st} + \mu_s + \mu_t + \epsilon_{st} \end{gathered} \]
The second model only includes economic determinants of food insecurity and the state and time fixed effects. The third model only includes predictors that were statically significant in the first model in addition to the state and time fixed effects.
Code
# create panel data
Regression_Data <- left_join(State_FI_Predictors, FI_Rates_by_State |> select(FI_Rate, observation), by = "observation") |>
filter(NAME != "District of Columbia", NAME != "Puerto Rico") |>
mutate(State = NAME)
# fit models
OLS_FI <- feols(FI_Rate ~
Avg_poverty_rate +
Avg_unemployment_rate +
Avg_median_income +
Avg_percent_hispanic +
Avg_percent_black +
Avg_homeownership_rate +
Avg_disability_rate |
State +
Year_Range, data = Regression_Data)
OLS_FI_Eco_Determinantes <- feols(FI_Rate ~
Avg_poverty_rate +
Avg_unemployment_rate +
Avg_median_income +
Avg_homeownership_rate +
Avg_disability_rate |
State +
Year_Range, data = Regression_Data)
OLS_FI_Only_Sig_Predictors <- feols(FI_Rate ~
Avg_poverty_rate +
Avg_unemployment_rate +
Avg_disability_rate |
State +
Year_Range, data = Regression_Data)
# display model summaries
modelsummary(list("All Variables" = OLS_FI,
"Just Eco Determinants" = OLS_FI_Eco_Determinantes,
"Just Sig Predictors" = OLS_FI_Only_Sig_Predictors),
stars = TRUE)| All Variables | Just Eco Determinants | Just Sig Predictors | |
|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||
| Avg_poverty_rate | 0.335+ | 0.342+ | 0.278+ |
| (0.189) | (0.189) | (0.166) | |
| Avg_unemployment_rate | 0.511* | 0.397+ | 0.389+ |
| (0.235) | (0.226) | (0.214) | |
| Avg_median_income | 0.000 | 0.000 | |
| (0.000) | (0.000) | ||
| Avg_percent_hispanic | 0.065 | ||
| (0.253) | |||
| Avg_percent_black | -0.421 | ||
| (0.256) | |||
| Avg_homeownership_rate | 0.087 | 0.098 | |
| (0.130) | (0.127) | ||
| Avg_disability_rate | 0.451+ | 0.387 | 0.414 |
| (0.258) | (0.256) | (0.250) | |
| Num.Obs. | 200 | 200 | 200 |
| R2 | 0.895 | 0.892 | 0.892 |
| R2 Adj. | 0.850 | 0.849 | 0.850 |
| R2 Within | 0.153 | 0.134 | 0.130 |
| R2 Within Adj. | 0.111 | 0.104 | 0.112 |
| AIC | -1212.0 | -1211.6 | -1214.6 |
| BIC | -1014.1 | -1020.3 | -1029.8 |
| RMSE | 0.01 | 0.01 | 0.01 |
| Std.Errors | IID | IID | IID |
| FE: State | X | X | X |
| FE: Year_Range | X | X | X |
We assessed the first model (which includes all of our variables) to be the best for generating the coefficients we needed for our analysis. This model predicted the largest portion of the state-to-state variation in food insecurity rates with an R2 of .895, had the lowest BIC, and the second lowest AIC. The range of residuals for our model is +/- 2.7 percentage points. The average food insecurity rate reported across all our observations is 12%, the average for NC is 13.1%. The residual range for North Carolina’s observations in our model is + 1.2 to - 1.5 percentage points.3
Generating food insecurity estimates for Western NC census tracts
We then downloaded census tract-level data for Western NC. Here we relied on ASC five-year estimates since ASC one-year estimates are not available for such small geographies.
Code
# Define variables
vars2024_acs5 <- c(
total_pop = "DP05_0001",
black_alone = "DP05_0045",
hispanic = "DP05_0090",
total_disability = "S1810_C01_001",
with_disability = "S1810_C02_001",
occupied_housing = "DP04_0002",
owner_occupied = "DP04_0046",
median_income = "B19013_001",
poverty_total = "B14006_001",
poverty_below = "B14006_002",
unemployment_rate = "S2301_C04_001"
)
# Get ACS data for all tracts in the state
FI_Predictors <- get_acs(
geography = "tract",
state = "NC",
variables = vars2024_acs5,
year = 2024,
survey = "acs5",
geometry = FALSE,
cache_table = TRUE
)
# Reshape to wide format
FI_Predictors <- FI_Predictors %>%
select(GEOID, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate)
# Calculate rates
FI_Predictors <- FI_Predictors %>%
mutate(
disability_rate = with_disability / total_disability,
homeownership_rate = owner_occupied / occupied_housing,
percent_black = black_alone / total_pop,
percent_hispanic = hispanic / total_pop,
poverty_rate = poverty_below / poverty_total,
unemployment_rate = unemployment_rate *0.01
)Finally, we were able to use the coefficients from our state-level model to predict the likely food insecurity rate for census tracts in Western North Carolina.4
Code
FI_Predictors <- FI_Predictors |>
mutate("FI_Rate" =
(poverty_rate * 0.335) +
(unemployment_rate * 0.511) +
(median_income * 0.00) +
(percent_hispanic * 0.065) +
(percent_black * -0.421) +
(homeownership_rate * 0.087) +
(disability_rate * 0.451) +
0.024 + # the intercept
-0.016 + # North Carolina fixed effect
-0.0098 # time period fixed effect
)
# Export data for use elsewhere
write_csv(FI_Predictors, "Data/FI_Rates.csv")
Footnotes
Rabbitt, M. P., Reed-Jones, M., Hales, L. J., Suttles, S., Burke, M. P. (2025). Household food security in the United States in 2024. (pages 3 and 4)↩︎
The specific variable codes change from year to year and can be found in the R code below. We used American Community Survey one-year estimates for each year except 2020 during which only the ACS five-year estimate is available.↩︎
The Map the Meal Gap method weighted the model using state population. We chose not to do this in part because this would give more weight to sates much larger than North Carolina. We found that when we used weights in our model that 3 out of 4 observations for North Carolina had larger residuals than in our unweighted model.↩︎
The intercept is 0.024, the state fixed effect for North Carolina is -0.016. The data we have for Western NC census tracts is from the ASC five-year estimates covering the years between 2020 and 2024. Therefore we add the weighted average of the time period fixed effect for the 2019 to 2021 (-0.023) and 2022 to 2024 (-0.001) time periods.↩︎