County-Level Geographic Differences in Adult Obesity Prevalence After Controlling for Known Risk Factors

1 Abstract

Adult obesity prevalence varies across U.S. counties. This study examines whether county-level geographic differences remain after accounting for measured behavioral, socioeconomic, food-related, and environmental characteristics. County-level data are assembled from the IHME datasets and supporting public data files maintained in the project’s GitHub repository. The primary outcome is an approximate county-level adult obesity prevalence measure calculated as the arithmetic mean of the male and female 2009 IHME county estimates. Predictors are selected from the available county-level measures of physical inactivity, alcohol use, educational attainment, unemployment, household income or wealth, food insecurity, sunlight exposure, environmental quality, and population.

The analysis uses descriptive statistics, exploratory visualizations, multiple ordinary least squares regression, and a state fixed-effects model. Model fit is assessed using adjusted R-squared, root mean squared error, a nested-model comparison, and regression diagnostics. Because the data are observational and aggregated at the county level, all findings are interpreted as associations rather than causal effects. The final interpretation will focus on the amount of variation explained by measured county characteristics and whether adding state indicators provides additional explanatory value.

The numerical findings in this abstract are generated from the current data when the document is knitted. The outcome definition and the years represented by each source must be considered when interpreting the results.

Keywords: adult obesity, county-level analysis, geographic variation, physical inactivity, socioeconomic factors, state fixed effects, spatial disparities

2 1. Introduction

Adult obesity is an important public-health issue in the United States, and prevalence differs across geographic areas. National and state summaries can conceal substantial variation among counties. Counties differ in economic resources, educational attainment, employment conditions, food security, behavioral patterns, and environmental characteristics. These differences provide an important context for examining county-level variation in adult obesity.

Many known characteristics have been associated with obesity, including physical activity, socioeconomic conditions, and food-related circumstances. The purpose of this study is not simply to identify whether these established factors are associated with obesity. Instead, the central question is whether geographic differences among U.S. counties remain after accounting for measured county-level characteristics.

This study combines publicly available county-level data from multiple sources and uses regression modeling to estimate the relationship between adult obesity prevalence and selected explanatory variables. A baseline ordinary least squares (OLS) model is compared with a model that includes state fixed effects. The comparison evaluates whether state-level geographic indicators provide additional explanatory information after adjustment for the measured risk factors.

The analysis is observational. It does not establish that county location or any individual predictor causes obesity. Because the unit of analysis is the county, results describe county-level associations and should not be interpreted as individual-level relationships.

2.1 1.1 Importance of the Study

Understanding geographic differences in obesity may help identify communities where observed prevalence is higher or lower than expected given measured county characteristics. If geographic differences remain after adjustment, this may indicate that additional contextual, policy-related, healthcare, demographic, or spatially patterned factors warrant investigation.

The study also demonstrates how data from multiple public sources can be integrated and analyzed using reproducible statistical methods in R. The emphasis is on model comparison, effect estimates, uncertainty, and diagnostic evaluation rather than statistical significance alone.

3 2. Literature Review

3.1 2.1 Obesity and Geographic Differences

Public-health research has examined obesity in relation to behavioral, socioeconomic, and environmental characteristics. A county-level perspective is useful because counties within the same state can differ substantially in population composition, economic conditions, access to resources, and health-related environments. The present study uses this geographic variation as the focus of its analysis.

3.2 2.2 Behavioral and Socioeconomic Characteristics

Physical activity is an important behavioral characteristic in the study of obesity. County-level differences in physical activity may be associated with differences in obesity prevalence. Economic and educational characteristics may also be related to community resources, opportunities, and health conditions. This study therefore considers physical inactivity, alcohol use, education, unemployment, and household income or wealth as potential explanatory or control variables.

3.3 2.3 Food-Related and Environmental Characteristics

Food insecurity is included as a county-level measure of food-related economic hardship. Environmental characteristics, including sunlight exposure and the Environmental Quality Index where available, are considered as additional contextual measures. These variables are not assumed to represent every aspect of the local food or physical environment; they are included because they are available in the project’s compiled public datasets.

3.4 2.4 Geographic and Spatial Variation

State fixed effects represent systematic differences among states that are not captured by the measured explanatory variables. If the fixed-effects model improves model fit, this suggests that state membership provides additional explanatory information. State fixed effects do not, by themselves, identify local spatial dependence between neighboring counties. A separate spatial-autocorrelation analysis would be required to assess that question.

3.5 2.5 Research Contribution

The contribution of this study is its focus on the amount of county-level geographic variation that remains after adjustment for measured behavioral, socioeconomic, food-related, and environmental characteristics. The study does not claim that the selected predictors or the geographic approach are unprecedented. Rather, it applies a consistent county-level modeling framework to assess the additional explanatory value of state-level geographic indicators.

4 3. Research Questions and Hypotheses

4.1 3.1 Primary Research Question

To what extent do county-level geographic differences in adult obesity prevalence remain after controlling for known behavioral, socioeconomic, food-related, and environmental factors?

4.2 3.2 Supporting Research Questions

1. How much variation in adult obesity prevalence exists across U.S. counties before accounting for other county-level characteristics?
2. To what extent are measured behavioral, socioeconomic, food-related, and environmental characteristics associated with county-level adult obesity prevalence?
3. Does adding state fixed effects improve the explanatory performance of the model after adjustment for the measured characteristics?
4. Do residual diagnostics indicate important departures from the assumptions of the OLS regression model?

4.3 3.3 Hypotheses

• H1: County-level adult obesity prevalence varies across U.S. counties.
• H2: The selected behavioral, socioeconomic, food-related, and environmental characteristics jointly explain a statistically meaningful portion of the variation in county-level adult obesity prevalence.
• H3: Adding state fixed effects improves model fit beyond the baseline model containing measured county-level characteristics.
• H4: After adjustment, residuals may still exhibit patterns that indicate unexplained county-level variation. This hypothesis will be evaluated descriptively; state fixed effects alone are not a test of neighboring-county spatial autocorrelation.

5 4. Data and Variables

5.1 4.1 Data Sources

The project repository contains the county-level source files used for this analysis. The files are read directly from the public GitHub repository so the analysis can be reproduced. The principal sources identified in the approved proposal include the Institute for Health Metrics and Evaluation (IHME), U.S. government datasets, the U.S. Environmental Protection Agency, Feeding America, and CDC environmental data.

The data files include:

• IHME_Obesity.csv
• IHME_PhysicalActivity.csv
• AlcoholConsumption.csv
• Education.csv
• Unemployment.csv
• Wealth.csv
• FoodInsecurity.csv
• Sunlight.csv
• EnvironmentalQualityIndex.csv
• population.csv

The analysis uses only variables that can be identified and merged reliably from the available source files. The code below prints the selected source columns and reports any optional measure that is not found. Missing variables are not silently replaced with unrelated measures.

5.2 4.2 Outcome Variable

The IHME obesity file contains separate male and female county prevalence estimates for 2009. For this analysis, the primary outcome is defined as the arithmetic mean of those two estimates:

\[ \text{Adult Obesity}_{i} = \frac{\text{Male Obesity}_{i,2009}+ \text{Female Obesity}_{i,2009}}{2} \]

This is an approximation, not a population-weighted combined prevalence estimate. A population-weighted overall estimate would require county-level male and female population weights that are aligned to the same year and definition. The outcome is expressed in percentage points.

5.3 4.3 Explanatory Variables

Variable	Role	Measurement / interpretation
Physical inactivity	Behavioral	100 minus the average male/female sufficient physical activity prevalence, where available
Alcohol use	Behavioral	County-level alcohol-use measure supplied in the project repository
Educational attainment	Socioeconomic	County-level college or bachelor’s attainment measure
Unemployment	Socioeconomic	County unemployment rate
Household income / wealth	Socioeconomic	County-level income or wealth measure
Food insecurity	Food-related	County-level food-insecurity measure
Sunlight exposure	Environmental	Average daily sunlight, where available
Environmental quality	Environmental	Environmental Quality Index, where available
Population	Contextual	County population; considered as log population in a sensitivity specification
State	Geographic	State indicator used for fixed effects

5.4 4.4 Geographic Data Integration

County records are matched using standardized state and county names. The original proposal described the use of five-digit FIPS codes. The current repository files do not all expose a consistently named FIPS field, so this version uses a normalized state-plus-county key. This approach avoids matching counties with the same name in different states, but it is not equivalent to a verified FIPS join. The merge diagnostics below report unmatched records and should be reviewed before the results are finalized.

repo_base <- "https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/"

read_project_csv <- function(filename) {
  url <- paste0(repo_base, filename)
  tryCatch(
    readr::read_csv(url, show_col_types = FALSE, name_repair = "unique"),
    error = function(e) {
      stop("Could not read ", filename, ". Check the repository filename and internet connection. Details: ",
           conditionMessage(e))
    }
  )
}

normalize_state <- function(x) {
  x <- stringr::str_squish(as.character(x))
  x <- stringr::str_replace_all(x, "\\.", "")
  abbr <- match(toupper(x), state.abb)
  x[!is.na(abbr)] <- state.name[abbr[!is.na(abbr)]]
  dplyr::case_when(
    toupper(x) %in% c("DC", "DISTRICT OF COLUMBIA") ~ "District of Columbia",
    TRUE ~ stringr::str_to_title(x)
  )
}

normalize_county <- function(x) {
  x <- as.character(x)
  x <- stringi::stri_trans_general(x, "Latin-ASCII")
  x <- stringr::str_to_lower(x)
  x <- stringr::str_replace(x, ",.*$", "")
  x <- stringr::str_remove(x, "\\s+(county|parish|borough|census area|municipality)$")
  x <- stringr::str_replace_all(x, "[^a-z0-9]", "")
  x
}

find_column <- function(data, patterns, required = TRUE, label = "variable") {
  nm <- names(data)
  hits <- unique(unlist(lapply(patterns, function(p) {
    which(stringr::str_detect(nm, stringr::regex(p, ignore_case = TRUE)))
  })))
  if (length(hits) == 0) {
    if (required) {
      stop("Could not identify ", label, ". Available columns: ",
           paste(nm, collapse = " | "))
    }
    return(NA_character_)
  }
  nm[hits[1]]
}

standardize_geo <- function(data, source_name = "dataset") {

  original_names <- names(data)

  # Normalize column names: remove spaces, underscores,
  # punctuation, and convert to lowercase
  normalized_names <- tolower(
    gsub("[^a-zA-Z0-9]", "", original_names)
  )

  # Identify state column
  state_idx <- which(
    normalized_names %in% c("state", "stabr", "statename", "stname")
  )[1]

  # Identify county column
  county_idx <- which(
    normalized_names %in% c(
      "county",
      "location",
      "areaname",
      "countyname",
      "countystate",
      "ctyname"
    )
  )[1]

  if (is.na(state_idx)) {
    stop(
      paste0(
        "Could not identify state column in ",
        source_name,
        ". Available columns: ",
        paste(original_names, collapse = " | ")
      )
    )
  }

  if (is.na(county_idx)) {
    stop(
      paste0(
        "Could not identify county column in ",
        source_name,
        ". Available columns: ",
        paste(original_names, collapse = " | ")
      )
    )
  }

  state_original <- original_names[state_idx]
  county_original <- original_names[county_idx]

  data %>%
    mutate(
      State = stringr::str_squish(as.character(.data[[state_original]])),
      County = stringr::str_squish(as.character(.data[[county_original]]))
    ) %>%
    mutate(
      County_key = paste(
        normalize_state(State),
        normalize_county(County),
        sep = "|"
      )
    )
}

numeric_value <- function(x) {
  readr::parse_number(as.character(x), na = c("", "NA", "N/A", "-", "—"))
}

make_one_variable <- function(filename, output_name, patterns,
                              required = FALSE, transform = identity) {
  raw <- read_project_csv(filename)
  geo <- standardize_geo(raw, filename)
  value_col <- find_column(
    geo, patterns, required = required,
    label = paste(output_name, "in", filename)
  )

  if (is.na(value_col)) {
    message("Optional variable not found; omitted: ", output_name,
            " (", filename, ")")
    return(NULL)
  }

  message("Mapped ", output_name, " from ", filename, " column: ", value_col)

  geo %>%
    transmute(
      State, County_key,
      !!output_name := transform(numeric_value(.data[[value_col]]))
    ) %>%
    distinct(State, County_key, .keep_all = TRUE)
}

5.5 4.5 Constructing the County-Level Analytical Dataset

# Outcome: average of the male and female 2009 obesity estimates.
obesity_raw <- read_project_csv("IHME_Obesity.csv")
obesity_geo <- standardize_geo(obesity_raw, "IHME_Obesity.csv")

male_obesity_col <- find_column(
  obesity_geo, "^Male obesity prevalence, 2009",
  label = "male 2009 obesity prevalence"
)
female_obesity_col <- find_column(
  obesity_geo, "^Female\\s+obesity\\s+prevalence,\\s*2009",
  label = "female 2009 obesity prevalence"
)

obesity_data <- obesity_geo %>%
  transmute(
    State, County_key,
    male_obesity_2009 = numeric_value(.data[[male_obesity_col]]),
    female_obesity_2009 = numeric_value(.data[[female_obesity_col]]),
    adult_obesity = rowMeans(
      cbind(
        numeric_value(.data[[male_obesity_col]]),
        numeric_value(.data[[female_obesity_col]])
      ),
      na.rm = TRUE
    )
  ) %>%
  mutate(
    adult_obesity = ifelse(is.nan(adult_obesity), NA_real_, adult_obesity)
  ) %>%
  distinct(State, County_key, .keep_all = TRUE)

# Physical inactivity: complement of sufficient physical activity prevalence.
activity_raw <- read_project_csv("IHME_PhysicalActivity.csv")
activity_geo <- standardize_geo(activity_raw, "IHME_PhysicalActivity.csv")

male_activity_col <- find_column(
  activity_geo,
  "^Male\\s+sufficient\\s+physical\\s+activity\\s+prevalence,\\s*2009",
  label = "male 2009 sufficient physical activity"
)
female_activity_col <- find_column(
  activity_geo,
  "^Female\\s+sufficient\\s+physical\\s+activity\\s+prevalence,\\s*2009",
  label = "female 2009 sufficient physical activity"
)

activity_data <- activity_geo %>%
  transmute(
    State, County_key,
    physical_inactivity = 100 - rowMeans(
      cbind(
        numeric_value(.data[[male_activity_col]]),
        numeric_value(.data[[female_activity_col]])
      ),
      na.rm = TRUE
    )
  ) %>%
  mutate(
    physical_inactivity = ifelse(
      is.nan(physical_inactivity), NA_real_, physical_inactivity
    )
  ) %>%
  distinct(State, County_key, .keep_all = TRUE)

# Additional candidate measures. A source variable is included only when a
# corresponding county-level column can be identified.
optional_sources <- list(
  make_one_variable(
    "AlcoholConsumption.csv", "alcohol_use",
    c("^Hvy$", "heavy.*(prevalence|rate)", "heavy drinking")
  ),
  make_one_variable(
    "Education.csv", "college_education",
    c("^College$", "bachelor", "college.*percent", "college.*rate")
  ),
  make_one_variable(
    "Unemployment.csv", "unemployment",
    c("^Unemployment_rate_2019$", "^Unemployment$", "unemployment.*2019")
  ),
  make_one_variable(
    "Wealth.csv", "median_income",
    c("^Income$", "median household income", "median.*income", "wealth")
  ),
  make_one_variable(
    "FoodInsecurity.csv", "food_insecurity",
    c("^FoodInsecurity$", "food insecurity.*(rate|percent|pct)?",
      "food insecurity")
  ),
  make_one_variable(
    "Sunlight.csv", "sunlight",
    c("^Sun$", "Avg Daily Sunlight", "average_daily_sunlight",
      "average daily sunlight")
  ),
  make_one_variable(
    "EnvironmentalQualityIndex.csv", "environmental_quality",
    c("^EQI$", "environmental quality index")
  ),
  make_one_variable(
    "population.csv", "population",
    c("^Population$", "POPESTIMATE2019", "population.*2019")
  )
)

optional_sources <- optional_sources[!vapply(optional_sources, is.null, logical(1))]

analysis_data <- obesity_data %>%
  left_join(activity_data, by = c("State", "County_key"))

for (source_data in optional_sources) {
  analysis_data <- analysis_data %>%
    left_join(source_data, by = c("State", "County_key"))
}

analysis_data <- analysis_data %>%
  mutate(
    log_population = if ("population" %in% names(.)) {
      ifelse(population > 0, log(population), NA_real_)
    } else {
      NA_real_
    }
  )

# Keep counties with a valid outcome.
analysis_data <- analysis_data %>%
  filter(is.finite(adult_obesity))

# Merge diagnostics.
merge_summary <- tibble(
  Measure = c(
    "Counties with obesity outcome",
    "Counties with physical inactivity",
    "Counties with missing physical inactivity"
  ),
  Count = c(
    nrow(analysis_data),
    sum(!is.na(analysis_data$physical_inactivity)),
    sum(is.na(analysis_data$physical_inactivity))
  )
)

knitr::kable(
  merge_summary,
  caption = "County-level data availability after merging"
) %>%
  kableExtra::kable_styling(full_width = FALSE)

County-level data availability after merging

Measure	Count
Counties with obesity outcome	3194
Counties with physical inactivity	3194
Counties with missing physical inactivity	0

# Show available variables and their missingness.
missingness_table <- tibble(
  Variable = names(analysis_data),
  Non_Missing = vapply(analysis_data, function(x) sum(!is.na(x)), integer(1)),
  Missing = vapply(analysis_data, function(x) sum(is.na(x)), integer(1))
)

knitr::kable(
  missingness_table,
  caption = "Missing-value summary for the merged analytical dataset"
) %>%
  kableExtra::kable_styling(full_width = FALSE)

Missing-value summary for the merged analytical dataset

Variable	Non_Missing	Missing
State	3194	0
County_key	3194	0
male_obesity_2009	3194	0
female_obesity_2009	3194	0
adult_obesity	3194	0
physical_inactivity	3194	0
college_education	0	3194
unemployment	0	3194
median_income	0	3194
food_insecurity	0	3194
sunlight	3153	41
population	3135	59
log_population	3135	59

5.6 4.6 Data Quality and Analysis Sample

The number of counties available for each measure may differ because the source datasets may have different geographic coverage and missing values. The regression sample is therefore defined by the observations with complete data for the outcome and the predictors included in the final model. The final paper should report this sample size and should not imply that every model uses all U.S. counties.

summary(analysis_data$adult_obesity)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.00   34.70   36.85   36.94   39.40   52.25

county_counts <- analysis_data %>%
  summarise(
    counties = n(),
    states = n_distinct(State),
    mean_obesity = mean(adult_obesity, na.rm = TRUE),
    sd_obesity = sd(adult_obesity, na.rm = TRUE),
    min_obesity = min(adult_obesity, na.rm = TRUE),
    median_obesity = median(adult_obesity, na.rm = TRUE),
    max_obesity = max(adult_obesity, na.rm = TRUE)
  )

knitr::kable(
  county_counts,
  digits = 2,
  caption = "Summary of the county-level obesity outcome"
) %>%
  kableExtra::kable_styling(full_width = FALSE)

Summary of the county-level obesity outcome

counties	states	mean_obesity	sd_obesity	min_obesity	median_obesity	max_obesity
3194	52	36.94	4.28	18	36.85	52.25

6 5. Methodology

6.1 5.1 Exploratory Data Analysis

Exploratory data analysis is used to describe the distribution of adult obesity prevalence and examine the relationships among the outcome and available explanatory variables. Descriptive statistics and visualizations are reviewed before fitting regression models.

ggplot(analysis_data, aes(x = adult_obesity)) +
  geom_histogram(bins = 30, color = "white") +
  labs(
    title = "Distribution of County-Level Adult Obesity Prevalence",
    x = "Adult obesity prevalence (%)",
    y = "Number of counties"
  ) +
  theme_minimal()

state_summary <- analysis_data %>%
  group_by(State) %>%
  summarise(
    counties = n(),
    mean_obesity = mean(adult_obesity, na.rm = TRUE),
    median_obesity = median(adult_obesity, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(desc(mean_obesity))

knitr::kable(
  head(state_summary, 15),
  digits = 2,
  caption = "States with the highest mean county-level obesity prevalence (top 15)"
) %>%
  kableExtra::kable_styling(full_width = FALSE)

States with the highest mean county-level obesity prevalence (top 15)

State	counties	mean_obesity	median_obesity
Mississippi	83	42.44	42.10
Louisiana	65	41.58	41.85
Alabama	68	41.05	40.40
South Carolina	47	40.23	39.70
Georgia	160	39.94	40.58
Kentucky	121	39.35	39.00
Arkansas	76	39.19	38.38
West Virginia	56	39.15	38.70
Tennessee	96	38.72	39.02
Oklahoma	78	38.56	38.40
Texas	255	38.25	38.30
North Carolina	101	37.96	37.10
Ohio	89	37.93	37.95
Missouri	116	37.88	37.68
Delaware	4	37.56	37.15

library(dplyr)
library(ggplot2)

state_extremes <- bind_rows(
  state_summary %>%
    slice_max(mean_obesity, n = 10) %>%
    mutate(Group = "Highest"),

  state_summary %>%
    slice_min(mean_obesity, n = 10) %>%
    mutate(Group = "Lowest")
) %>%
  arrange(mean_obesity) %>%
  mutate(State = factor(State, levels = State))

ggplot(state_extremes,
       aes(x = State, y = mean_obesity, fill = Group)) +
  geom_col(width = 0.7) +
  coord_flip() +
  labs(
    title = "States with the Highest and Lowest Mean County Obesity",
    x = NULL,
    y = "Mean county obesity prevalence (%)",
    fill = "State group"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "bottom"
  )

6.2 5.2 Correlation Analysis

Pearson correlations are used as an exploratory measure of linear association among numeric variables. Correlation does not establish causation and does not account for other variables in the model.

numeric_data <- analysis_data %>%
  select(where(is.numeric)) %>%
  select(where(~ sum(!is.na(.x)) > 10))

correlation_matrix <- cor(
  numeric_data,
  use = "pairwise.complete.obs"
)

corrplot_available <- requireNamespace("corrplot", quietly = TRUE)

if (corrplot_available) {
  corrplot::corrplot(
    correlation_matrix,
    method = "color",
    type = "upper",
    tl.col = "black",
    tl.cex = 0.75,
    addCoef.col = "black",
    number.cex = 0.5
  )
} else {
  knitr::kable(
    round(correlation_matrix, 2),
    caption = "Pearson correlation matrix"
  )
}

6.3 5.3 Multiple Linear Regression

The baseline model uses adult obesity prevalence as the dependent variable and the available county-level explanatory variables as predictors. The model is specified using the predictors that are successfully identified in the repository and have sufficient non-missing observations.

The general model is:

\[ Obesity_i = \beta_0 + \sum_{k=1}^{K}\beta_k X_{ki} + \epsilon_i \]

where \(Obesity_i\) is the adult obesity prevalence measure for county \(i\), \(X_{ki}\) represents the measured county-level characteristics, and \(\epsilon_i\) is the error term.

candidate_predictors <- c(
  "physical_inactivity",
  "alcohol_use",
  "college_education",
  "unemployment",
  "median_income",
  "food_insecurity",
  "sunlight",
  "environmental_quality",
  "population"
)

# Identify predictors available in the dataset
existing_predictors <- intersect(
  candidate_predictors,
  names(analysis_data)
)

# Check predictor quality
predictor_diagnostics <- purrr::map_dfr(
  existing_predictors,
  function(v) {

    x <- analysis_data[[v]]

    tibble(
      Variable = v,
      Non_Missing = sum(!is.na(x)),
      Distinct_Values = dplyr::n_distinct(x, na.rm = TRUE),
      Data_Type = class(x)[1]
    )
  }
)

knitr::kable(
  predictor_diagnostics,
  caption = "Regression predictor diagnostics"
) %>%
  kableExtra::kable_styling(full_width = FALSE)

Regression predictor diagnostics

Variable	Non_Missing	Distinct_Values	Data_Type
physical_inactivity	3194	756	numeric
college_education	0	0	numeric
unemployment	0	0	numeric
median_income	0	0	numeric
food_insecurity	0	0	numeric
sunlight	3153	3141	numeric
population	3135	3081	numeric

# Select usable numerical predictors
available_predictors <- predictor_diagnostics %>%
  filter(
    Non_Missing >= 100,
    Distinct_Values > 1,
    Data_Type %in% c("numeric", "integer")
  ) %>%
  pull(Variable)

cat("Usable predictors:",
    paste(available_predictors, collapse = ", "))

## Usable predictors: physical_inactivity, sunlight, population

if (length(available_predictors) < 3) {
  stop(
    "Fewer than three usable predictors. ",
    "Review the diagnostic table before continuing."
  )
}

# Create complete-case regression dataset
model_variables <- c(
  "adult_obesity",
  "State",
  "County_key",
  available_predictors
)

model_data <- analysis_data %>%
  select(all_of(model_variables)) %>%
  drop_na() %>%
  mutate(State = factor(State))

if (nrow(model_data) < 100) {
  stop("Regression sample contains fewer than 100 counties.")
}

# Baseline OLS model
baseline_formula <- as.formula(
  paste(
    "adult_obesity ~",
    paste(available_predictors, collapse = " + ")
  )
)

model_ols <- lm(
  baseline_formula,
  data = model_data
)

# State fixed-effects model
state_formula <- as.formula(
  paste(
    "adult_obesity ~",
    paste(available_predictors, collapse = " + "),
    "+ State"
  )
)

model_state <- lm(
  state_formula,
  data = model_data
)

cat("\nRegression sample:",
    nrow(model_data), "counties")

## 
## Regression sample: 3101 counties

cat("\nStates represented:",
    n_distinct(model_data$State))

## 
## States represented: 49

cat("\nBaseline adjusted R-squared:",
    round(summary(model_ols)$adj.r.squared, 4))

## 
## Baseline adjusted R-squared: 0.707

cat("\nState model adjusted R-squared:",
    round(summary(model_state)$adj.r.squared, 4))

## 
## State model adjusted R-squared: 0.7971

6.4 5.4 State Fixed-Effects Model

The geographic model adds state indicators to the baseline regression:

\[ Obesity_{is} = \beta_0 + \sum_{k=1}^{K}\beta_k X_{kis} + \gamma_s + \epsilon_{is} \]

where \(\gamma_s\) represents the state fixed effect. One state is treated as the reference category by the regression software.

The state fixed-effects model tests whether state membership adds explanatory information after adjustment for the measured county-level variables. It does not establish that state location causes differences in obesity prevalence.

6.5 5.5 Model Evaluation

The models are compared using adjusted R-squared and RMSE. A nested-model F-test evaluates whether the state indicators collectively improve the in-sample fit. These measures describe explanatory performance in the analysis sample and should not be confused with out-of-sample predictive validation.

rmse <- function(model) {
  sqrt(mean(residuals(model)^2))
}

model_comparison <- tibble(
  Model = c("Baseline OLS", "OLS with state fixed effects"),
  Observations = c(nobs(model_ols), nobs(model_state)),
  R_Squared = c(summary(model_ols)$r.squared,
                summary(model_state)$r.squared),
  Adjusted_R_Squared = c(summary(model_ols)$adj.r.squared,
                         summary(model_state)$adj.r.squared),
  RMSE = c(rmse(model_ols), rmse(model_state))
)

knitr::kable(
  model_comparison,
  digits = 4,
  caption = "Comparison of baseline and state fixed-effects models"
) %>%
  kableExtra::kable_styling(full_width = FALSE)

Comparison of baseline and state fixed-effects models

Model	Observations	R_Squared	Adjusted_R_Squared	RMSE
Baseline OLS	3101	0.7073	0.7070	2.3016
OLS with state fixed effects	3101	0.8004	0.7971	1.9003

state_model_test <- anova(model_ols, model_state)
knitr::kable(
  broom::tidy(state_model_test),
  digits = 4,
  caption = "Nested-model F-test: additional contribution of state fixed effects"
) %>%
  kableExtra::kable_styling(full_width = FALSE)

Nested-model F-test: additional contribution of state fixed effects

term	df.residual	rss	df	sumsq	statistic	p.value
adult_obesity ~ physical_inactivity + sunlight + population	3097	16427.01	NA	NA	NA	NA
adult_obesity ~ physical_inactivity + sunlight + population + State	3049	11197.89	48	5229.118	29.6625	0

6.6 5.6 Regression Coefficients and Uncertainty

The coefficients are interpreted as adjusted county-level associations. A coefficient represents the estimated difference in adult obesity prevalence associated with a one-unit difference in a predictor, holding the other included variables constant. Units must be considered when interpreting the size of each coefficient.

Heteroscedasticity-consistent standard errors (HC3) are reported as a robustness measure.

robust_results <- lmtest::coeftest(
  model_state,
  vcov. = sandwich::vcovHC(model_state, type = "HC3")
)

coefficient_table <- broom::tidy(robust_results) %>%
  rename(
    Estimate = estimate,
    Std_Error_HC3 = std.error,
    t_value = statistic,
    p_value = p.value
  ) %>%
  mutate(
    CI_Lower = Estimate - qt(0.975, df.residual(model_state)) * Std_Error_HC3,
    CI_Upper = Estimate + qt(0.975, df.residual(model_state)) * Std_Error_HC3
  )

knitr::kable(
  coefficient_table,
  digits = 4,
  caption = "State fixed-effects regression coefficients with HC3 robust standard errors"
) %>%
  kableExtra::kable_styling(full_width = FALSE, font_size = 13)

State fixed-effects regression coefficients with HC3 robust standard errors

term	Estimate	Std_Error_HC3	t_value	p_value	CI_Lower	CI_Upper
(Intercept)	-1.8414	NaN	NaN	NaN	NaN	NaN
physical_inactivity	0.6658	NaN	NaN	NaN	NaN	NaN
sunlight	0.0003	NaN	NaN	NaN	NaN	NaN
population	0.0000	NaN	NaN	NaN	NaN	NaN
StateArizona	2.0828	NaN	NaN	NaN	NaN	NaN
StateArkansas	0.4114	NaN	NaN	NaN	NaN	NaN
StateCalifornia	2.1268	NaN	NaN	NaN	NaN	NaN
StateColorado	-0.3809	NaN	NaN	NaN	NaN	NaN
StateConnecticut	-0.0820	NaN	NaN	NaN	NaN	NaN
StateDelaware	3.0589	NaN	NaN	NaN	NaN	NaN
StateDistrict of Columbia	-1.0196	NaN	NaN	NaN	NaN	NaN
StateFlorida	0.6223	NaN	NaN	NaN	NaN	NaN
StateGeorgia	1.2014	NaN	NaN	NaN	NaN	NaN
StateIdaho	3.0474	NaN	NaN	NaN	NaN	NaN
StateIllinois	2.3759	NaN	NaN	NaN	NaN	NaN
StateIndiana	1.2373	NaN	NaN	NaN	NaN	NaN
StateIowa	2.1367	NaN	NaN	NaN	NaN	NaN
StateKansas	1.8379	NaN	NaN	NaN	NaN	NaN
StateKentucky	-0.2502	NaN	NaN	NaN	NaN	NaN
StateLouisiana	0.6253	NaN	NaN	NaN	NaN	NaN
StateMaine	4.3100	NaN	NaN	NaN	NaN	NaN
StateMaryland	2.9480	NaN	NaN	NaN	NaN	NaN
StateMassachusetts	-0.3727	NaN	NaN	NaN	NaN	NaN
StateMichigan	5.4586	NaN	NaN	NaN	NaN	NaN
StateMinnesota	3.0931	NaN	NaN	NaN	NaN	NaN
StateMississippi	0.2031	NaN	NaN	NaN	NaN	NaN
StateMissouri	1.8512	NaN	NaN	NaN	NaN	NaN
StateMontana	2.4744	NaN	NaN	NaN	NaN	NaN
StateNebraska	2.8708	NaN	NaN	NaN	NaN	NaN
StateNevada	1.7372	NaN	NaN	NaN	NaN	NaN
StateNew Hampshire	3.9369	NaN	NaN	NaN	NaN	NaN
StateNew Jersey	-0.6874	NaN	NaN	NaN	NaN	NaN
StateNew Mexico	0.2666	NaN	NaN	NaN	NaN	NaN
StateNew York	2.4980	NaN	NaN	NaN	NaN	NaN
StateNorth Carolina	0.9091	NaN	NaN	NaN	NaN	NaN
StateNorth Dakota	2.5510	NaN	NaN	NaN	NaN	NaN
StateOhio	3.0466	NaN	NaN	NaN	NaN	NaN
StateOklahoma	1.0174	NaN	NaN	NaN	NaN	NaN
StateOregon	5.5529	NaN	NaN	NaN	NaN	NaN
StatePennsylvania	3.3849	NaN	NaN	NaN	NaN	NaN
StateRhode Island	-1.4941	NaN	NaN	NaN	NaN	NaN
StateSouth Carolina	2.8095	NaN	NaN	NaN	NaN	NaN
StateSouth Dakota	2.4899	NaN	NaN	NaN	NaN	NaN
StateTennessee	-2.7053	NaN	NaN	NaN	NaN	NaN
StateTexas	0.6439	NaN	NaN	NaN	NaN	NaN
StateUtah	1.7079	NaN	NaN	NaN	NaN	NaN
StateVermont	3.0854	NaN	NaN	NaN	NaN	NaN
StateVirginia	1.2315	NaN	NaN	NaN	NaN	NaN
StateWashington	5.9157	NaN	NaN	NaN	NaN	NaN
StateWest Virginia	0.2676	NaN	NaN	NaN	NaN	NaN
StateWisconsin	4.9081	NaN	NaN	NaN	NaN	NaN
StateWyoming	1.9160	NaN	NaN	NaN	NaN	NaN

6.7 5.7 Regression Diagnostics

Regression diagnostics are used to evaluate whether the fitted models show important departures from the assumptions of linear regression. The analysis examines residual patterns, approximate normality, heteroscedasticity, and multicollinearity. Statistical significance should be interpreted cautiously because counties may not be statistically independent.

old_par <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
plot(model_state)

par(old_par)

vif_values <- tryCatch(
  car::vif(model_state),
  error = function(e) NULL
)

if (!is.null(vif_values)) {
  if (is.matrix(vif_values)) {
    vif_output <- as.data.frame(vif_values) %>%
      tibble::rownames_to_column("Variable")
  } else {
    vif_output <- tibble(
      Variable = names(vif_values),
      VIF = as.numeric(vif_values)
    )
  }

  knitr::kable(
    vif_output,
    digits = 3,
    caption = "Variance inflation factor diagnostics"
  ) %>%
    kableExtra::kable_styling(full_width = FALSE)
} else {
  cat("VIF could not be calculated for this model specification.")
}

Variance inflation factor diagnostics

Variable	GVIF	Df	GVIF^(1/(2*Df))
physical_inactivity	2.512	1	1.585
sunlight	12.270	1	3.503
population	1.096	1	1.047
State	29.607	48	1.036

6.8 5.8 Residual Geographic Variation

A residual is the difference between observed and model-predicted obesity prevalence:

\[ Residual_i = Observed_i - Predicted_i \]

Positive residuals indicate that the observed county prevalence is higher than the model predicts. Negative residuals indicate that observed prevalence is lower than predicted. Residuals are examined to identify counties with larger deviations from the model.

model_data <- model_data %>%
  mutate(
    fitted_obesity = fitted(model_state),
    residual = residuals(model_state),
    absolute_residual = abs(residual)
  )

residual_summary <- model_data %>%
  summarise(
    Mean_residual = mean(residual),
    SD_residual = sd(residual),
    Mean_absolute_residual = mean(absolute_residual),
    Minimum_residual = min(residual),
    Maximum_residual = max(residual)
  )

knitr::kable(
  residual_summary,
  digits = 3,
  caption = "Summary of residuals from the state fixed-effects model"
) %>%
  kableExtra::kable_styling(full_width = FALSE)

Summary of residuals from the state fixed-effects model

Mean_residual	SD_residual	Mean_absolute_residual	Minimum_residual	Maximum_residual
0	1.901	1.448	-10.788	7.618

ggplot(model_data, aes(x = fitted_obesity, y = residual)) +
  geom_point(alpha = 0.45) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(
    title = "Residuals Versus Fitted Values",
    x = "Predicted adult obesity prevalence (%)",
    y = "Residual (percentage points)"
  ) +
  theme_minimal()

largest_residuals <- model_data %>%
  mutate(County = analysis_data$County_key[match(
    paste(State, County_key), paste(analysis_data$State, analysis_data$County_key)
  )]) %>%
  arrange(desc(absolute_residual)) %>%
  select(State, County, adult_obesity, fitted_obesity, residual) %>%
  slice_head(n = 15)

knitr::kable(
  largest_residuals,
  digits = 2,
  caption = "Counties with the largest absolute model residuals"
) %>%
  kableExtra::kable_styling(full_width = FALSE)

Counties with the largest absolute model residuals

State	County	adult_obesity	fitted_obesity	residual
Virginia	Virginia|fallschurchcity	18.00	28.79	-10.79
Pennsylvania	Pennsylvania|montour	27.10	37.47	-10.37
California	California|sanfrancisco	19.00	29.17	-10.17
New York	New York|richmond	28.95	38.38	-9.43
Virginia	Virginia|fairfaxcity	22.70	31.66	-8.96
New York	New York|queens	29.10	37.39	-8.29
New York	New York|nassau	28.70	36.53	-7.83
Virginia	Virginia|portsmouthcity	45.35	37.73	7.62
Alabama	Alabama|marion	37.70	45.28	-7.58
North Carolina	North Carolina|orange	23.45	30.89	-7.44
Maryland	Maryland|montgomery	24.40	31.65	-7.25
New Mexico	New Mexico|mckinley	41.70	34.63	7.07
New York	New York|westchester	25.50	32.44	-6.94
New York	New York|kings	32.05	38.80	-6.75
North Carolina	North Carolina|ashe	32.30	38.83	-6.53

7 6. Results

7.1 6.1 County-Level Variation in Adult Obesity

The descriptive analysis reports the distribution of the county-level adult obesity outcome. The final interpretation should describe the range, central tendency, and variability shown in the generated summary and histogram.

The analysis includes 3194 counties with an available outcome measure. The mean county-level adult obesity prevalence is 36.94%, and the standard deviation is 4.28 percentage points.

7.2 6.2 Regression Results

The baseline OLS model has an adjusted R-squared of 0.707. The state fixed-effects model has an adjusted R-squared of 0.797. The corresponding RMSE values are 2.302 and 1.9 percentage points, respectively.

The nested-model F-test evaluates whether adding the state indicators improves model fit. Its result should be interpreted together with the change in adjusted R-squared and RMSE, rather than relying only on the p-value.

7.3 6.3 Geographic Differences After Adjustment

The state fixed-effects model evaluates whether state indicators collectively explain additional variation after accounting for the measured county characteristics. A statistically significant nested-model test would indicate that state indicators add explanatory information in this sample. It would not demonstrate that every state differs from every other state, nor would it identify the specific local mechanisms responsible for the differences.

The residual analysis identifies counties where observed prevalence differs most from model predictions. These residual differences may reflect unmeasured county characteristics, measurement error, differences in source years, model misspecification, or spatially patterned factors.

7.4 6.4 Model Diagnostics

The diagnostic plots and VIF results should be reviewed before final interpretation. Potential issues include nonlinearity, heteroscedasticity, influential observations, and multicollinearity. HC3 robust standard errors help address heteroscedasticity in coefficient uncertainty estimates, but they do not correct spatial dependence between neighboring counties.

8 7. Discussion

The purpose of this study is to determine whether geographic differences in adult obesity prevalence remain after controlling for measured county-level characteristics. The comparison between the baseline OLS model and the state fixed-effects model provides the primary evidence for evaluating the research question.

Complete this discussion after reviewing the knitted results. Specifically, explain:

1. Whether the county-level obesity outcome shows substantial variation.
2. Which measured predictors have the largest estimated associations, including their direction, units, and confidence intervals.
3. Whether the state fixed-effects model improves adjusted R-squared and RMSE.
4. Whether the nested-model F-test supports the conclusion that state indicators add explanatory information.
5. What the residual diagnostics suggest about unexplained county-level differences.

The findings should be compared with the literature reviewed in Section 2. Avoid describing an association as a causal effect. A statistically significant coefficient does not necessarily indicate a large or practically important association.

9 8. Limitations

Several limitations should be considered when interpreting the findings.

First, the study uses aggregated county-level data. County-level associations cannot automatically be generalized to individuals; doing so would risk an ecological fallacy.

Second, the outcome is calculated as the arithmetic mean of male and female IHME 2009 obesity estimates. It is not a population-weighted overall prevalence estimate. This definition should be kept explicit throughout the paper.

Third, the source datasets may represent different years. Differences in the timing of obesity, physical activity, socioeconomic, food-related, and environmental measurements may affect the estimated relationships.

Fourth, the analysis relies on county-name and state-name standardization for files without a consistently available FIPS identifier. Although state-plus-county matching reduces ambiguity, a verified five-digit FIPS crosswalk would be preferable for the final data integration.

Fifth, the selected variables do not capture every factor associated with obesity. Potentially relevant omitted characteristics include healthcare access, transportation infrastructure, demographic composition, housing, local policy, and neighborhood-level resources.

Sixth, state fixed effects capture average differences between states but do not directly test spatial autocorrelation among neighboring counties. A separate spatial analysis using validated county geometries and a defensible spatial-weights matrix would be needed to assess neighboring-county dependence.

Finally, OLS estimates may be sensitive to influential observations, multicollinearity, heteroscedasticity, and measurement error. The diagnostic results should therefore be considered alongside the coefficient estimates.

10 9. Conclusion

This study examines county-level geographic differences in adult obesity prevalence after adjustment for measured behavioral, socioeconomic, food-related, and environmental characteristics. It uses a sequential modeling strategy that compares a baseline OLS regression with a model containing state fixed effects.

The final conclusion should be completed after the data and model outputs have been reviewed. It should state how much variation the baseline model explains, whether the state fixed-effects model adds explanatory value, which measured associations are most important, and what geographic differences remain unexplained.

The findings should be interpreted as county-level observational associations. They may help identify geographic patterns that warrant additional public-health investigation, but they do not establish causal effects or explain all sources of obesity differences across communities.

11 10. Reproducibility and Final Checks

The analysis is designed to be reproducible from the public project repository. Before submission, confirm that the repository files have not changed, that each source column was mapped correctly, and that the final complete-case sample is appropriate for the research question.

final_checks <- tibble(
  Check = c(
    "Outcome variable present",
    "At least 3 predictors included",
    "Regression sample has at least 100 counties",
    "More than one state represented",
    "Baseline and state models use same observations",
    "Model coefficients are finite"
  ),
  Passed = c(
    "adult_obesity" %in% names(model_data),
    length(available_predictors) >= 3,
    nrow(model_data) >= 100,
    n_distinct(model_data$State) > 1,
    nobs(model_ols) == nobs(model_state),
    all(is.finite(coef(model_state)))
  )
)

knitr::kable(
  final_checks,
  caption = "Reproducibility and model validation checks"
) %>%
  kableExtra::kable_styling(full_width = FALSE)

Reproducibility and model validation checks

Check	Passed
Outcome variable present	TRUE
At least 3 predictors included	TRUE
Regression sample has at least 100 counties	TRUE
More than one state represented	TRUE
Baseline and state models use same observations	TRUE
Model coefficients are finite	TRUE

12 References

Ahern, M., Brown, C., & Dukas, S. (2011). A national study of the association between food environments and county-level health outcomes. The Journal of Rural Health. [Verify volume, issue, pages, and DOI before submission.]

An, R., Li, X., & Jiang, N. (2017). Geographical variations in the environmental determinants of physical inactivity among U.S. adults. [Verify complete publication details before submission.]

Centers for Disease Control and Prevention. (n.d.). North America Land Data Assimilation System (NLDAS) daily sunlight (kJ/m²), 1979–2011. https://wonder.cdc.gov/controller/datarequest/D80

Feeding America. (n.d.). Map the Meal Gap: County-level food insecurity data. https://www.feedingamerica.org/research/map-the-meal-gap/by-county

Hammarström, A., Janlert, U., & Winefield, H. (2014). Length of unemployment and health-related outcomes. [Verify complete publication details before submission.]

Holick, M. F., & Wacker, M. (2013). Sunlight and vitamin D. [Verify publisher and edition or article details before submission.]

Institute for Health Metrics and Evaluation. (n.d.). United States alcohol use prevalence by county, 2002–2012. Global Health Data Exchange. http://ghdx.healthdata.org/record/ihme-data/united-states-alcohol-use-prevalence-county-2002-2012

Institute for Health Metrics and Evaluation. (n.d.). United States physical activity and obesity prevalence by county. Global Health Data Exchange. https://ghdx.healthdata.org/

U.S. Census Bureau. (2019). Small Area Income and Poverty Estimates (SAIPE): State and county estimates. https://www.census.gov/data/datasets/2019/demo/saipe/2019-state-and-county.html

U.S. Department of Agriculture, Economic Research Service. (n.d.). County-level data sets. https://www.ers.usda.gov/data-products/county-level-data-sets/download-data/

U.S. Environmental Protection Agency. (n.d.). Environmental Quality Index (EQI), 2006–2010. https://edg.epa.gov/EPADataCommons/public/ORD/CPHEA/EQI_2006_2010/

---

Submission note: Complete the discussion and conclusion using the actual output from the final knit. Verify every literature reference’s full bibliographic details and add in-text citations where the cited studies are discussed. The final paper should not be submitted with unresolved source mappings, incomplete references, or unreviewed model diagnostics.