Mini Project

Identify your research question by picking an outcome of interest and some covariates with data available at the county/state level.

In 2019, to what extent do state-level poverty rates and low educational attainment (adults aged 25+ without a high school diploma) predict state-level variations in Hepatitis A rates?

Data Source: CDC NCHHSTP AtlasPlus

library(tigris)
## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## βœ” forcats   1.0.1     βœ” readr     2.2.0
## βœ” ggplot2   4.0.3     βœ” stringr   1.6.0
## βœ” lubridate 1.9.5     βœ” tibble    3.3.0
## βœ” purrr     1.2.2     βœ” tidyr     1.3.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(here)
## here() starts at /Users/kymoore/Downloads/SISMID_2026_spatial_statistics_module-main-3
library(tmap)

#poverty <- read.csv("/Users/kymoore/Downloads/Mini Project/MapData-3.csv")
#hepa <- read.csv("/Users/kymoore/Downloads/Mini Project/MapData-2.csv")
#education <- read.csv("/Users/kymoore/Downloads/Mini Project/MapData-4.csv")
poverty <- read.csv("/Users/kymoore/Downloads/MapData-3.csv", skip=9)
hepa <- read.csv("/Users/kymoore/Downloads/MapData-2.csv", skip=9)
education <- read.csv("/Users/kymoore/Downloads/MapData-4.csv", skip=9)



#here::i_am("SISMID.Gentrification Mapping Lab.2026.Rmd")

#Data Wrangling

library(dplyr)
library(stringr)

#Tasks
#1. Identify the variables needed from the downloaded datasets.
#2. Check issues related to missing data.
  #check documentation
  #if its missing becuase its such a low number you can set it (ex: 5 set it to 2.5) otherwise set it to NA
#3. Examine summary statistics of the variables and create some exploratory plots (boxplot, histogram, scatter plot).

# 1. Prepare Dataset A (Select and Rename)
# Syntax: new_name = old_name
poverty_sub <- poverty %>%
  select(
    FIPS, 
    poverty_p = Percent,
    poverty_num = Numerator, 
    pov_den = Denominator
  )

# 2. Prepare Dataset B (Select and Rename)
education_sub <- education %>%
  select(
    FIPS, 
    education_p = Percent,
    education_num = Numerator, 
    education_den = Denominator
  )

# 3. Merge everything into the Main Dataset
hepa_final <- hepa %>%
  left_join(poverty_sub, by = "FIPS") %>%
  left_join(education_sub, by = "FIPS")

hepa_final
#remove punctuation
hepa_final_1 <- hepa_final %>% 
  # 1. Clean and convert numeric columns
  mutate(across(
    c(Cases, Rate.per.100000, poverty_p, poverty_num, education_p, education_num), 
    ~ {
      x <- na_if(.x, "Data not available")
      x <- str_remove_all(x, "[%, ]")
      as.numeric(x)
    }
  )) %>% 
  # 2. Convert fips to numeric and keep only states + DC (< 60)
  filter(as.numeric(FIPS) < 60)

# Replaces with NA
#hepa_final_1 <- hepa_final %>% 
#  mutate(across(c(Cases, Rate.per.100000, poverty_p, poverty_num, education_p, education_num), ~as.numeric(na_if(.x, "Data not available"))))


hist(hepa_final_1$Cases)

hist(hepa_final_1$Rate.per.100000)

#right skewed
summary(hepa_final_1$Cases)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0    16.0   128.0   369.5   375.0  3392.0
boxplot(hepa_final_1$Cases)

#outliers


plot(x = hepa_final_1$education_p, y = hepa_final_1$Rate.per.100000,
     main = "Low Educational Attainment vs. HepA Rate",
     xlab = "% of Population 25 years and older w/o HS diploma",
     ylab = "HepA Rate per 100,000",
     col = "black",          # Change point color
     pch = 19)              # Change point shape to solid circle

# Optional: Add a linear regression trend line
#abline(lm(education_p ~ Rate.per.100000, data = hepa_final_1), col = "red", lwd = 2)
# Corrected formula (predicts Rate from education):
abline(lm(Rate.per.100000 ~ education_p, data = hepa_final_1), col = "red", lwd = 2)

plot(x = hepa_final_1$poverty_p, y = hepa_final_1$Rate.per.100000,
     main = "Poverty vs. HepA Rate",
     xlab = "% of Households living below the federal poverty level",
     ylab = "HepA Rate per 100,000",
     col = "black",          # Change point color
     pch = 19)              # Change point shape to solid circle

# Optional: Add a linear regression trend line
#abline(lm(poverty_p ~ Rate.per.100000, data = hepa_final_1), col = "red", lwd = 2)
abline(lm(Rate.per.100000 ~ poverty_p, data = hepa_final_1), col = "red", lwd = 2)

### Calculating the standardized differences
hepa_final_1 <- hepa_final_1 %>%
  mutate(
    # Standardized Hepatitis A Rate (Z-score)
    rate_z = (Rate.per.100000 - mean(Rate.per.100000, na.rm = TRUE)) / sd(Rate.per.100000, na.rm = TRUE),
    
     # Standardized Education
    edu_z = (education_p - mean(education_p, na.rm = TRUE)) / sd(education_p, na.rm = TRUE),
    
    # Standardized Poverty
    poverty_z = (poverty_p - mean(poverty_p, na.rm = TRUE)) / sd(poverty_p, na.rm = TRUE),
    
    # If you are taking the difference between two standardized proportions (e.g., Poverty vs Education):
    poverty_edu_diff_z = poverty_z - edu_z
  )

hist(hepa_final_1$rate_z)

hist(hepa_final_1$edu_z)

hist(hepa_final_1$poverty_z)

hist(hepa_final_1$poverty_edu_diff_z)

# 1. Set up a grid of 1 row and 3 columns
par(mfrow = c(2, 2))

# 2. Draw each histogram
hist(hepa_final_1$rate_z, 
     main = "Standardized Hepatitis A Rate", 
     xlab = "Z-score", 
     col = "lightblue")

hist(hepa_final_1$poverty_z, 
     main = "Standardized Poverty", 
     xlab = "Z-score", 
     col = "salmon")

hist(hepa_final_1$edu_z, 
     main = "Standardized Education", 
     xlab = "Z-score", 
     col = "lightgreen")

hist(hepa_final_1$poverty_edu_diff_z, 
     main = "Standardized Covariate Differences", 
     xlab = "Z-score", 
     col = "purple")

# 3. Reset plotting grid back to default (1x1)
par(mfrow = c(1, 1))

#Creating Geospatial Datasets and Mapping

library (tigris)
options(tigris_use_cache = TRUE)


#option cb = TRUE to download the coarser version for speed
#we choose the year 2020
#dat_sf <- counties(state = c("MD", "PA", "NY", "NJ"), cb = TRUE, year = 2020)
#all_states_sf <- states(cb = TRUE, year = 2020)

#plot(all_states_sf)
#plot(all_states_sf["STATEFP"])

# 1. Plot all attributes for US states + DC only
#all_states_sf %>% 
 # filter(as.numeric(STATEFP) < 60) %>% 
  #plot()

# 2. Plot ONLY the STATEFP attribute map
#all_states_sf %>% 
 # filter(as.numeric(STATEFP) < 60) %>% 
#  select(STATEFP) %>% 
 # plot()


# Pulls US states with AK & HI automatically resized and shifted
#us_shifted <- states(cb = TRUE, resolution = "20m") %>% 
 # shift_geometry() # Shifts AK and HI under the Lower 48!

# 1. Load spatial geometry for US states and shift AK/HI
us_shifted <- states(cb = TRUE, resolution = "20m", year = 2020) %>% 
  filter(as.numeric(STATEFP) < 60) %>%  # Keep only states + DC
  shift_geometry()

# Plot using ggplot2
ggplot(us_shifted) +
  geom_sf(fill = "grey80", color = "white", size = 0.2) +
  theme_void() + # Removes empty margins completely 
  theme(legend.position = "none") # Drops the legend box

# 1. Standardize FIPS to a 2-digit character vector with leading zeros
hepa_final_1 <- hepa_final_1 %>%
  mutate(FIPS = str_pad(as.character(FIPS), width = 2, side = "left", pad = "0"))

# 2. Join your clean spatial geometry with your non-spatial dataset
# Note: Ensure state FIPS codes match (e.g., us_shifted$STATEFP == hepa_final_1$FIPS)
hepa_sf <- us_shifted %>% 
  left_join(hepa_final_1, by = c("STATEFP" = "FIPS"))

Now to map the outcome, using the mapping package tmap’s routines tm_shape and tm_fill. tmap builds the maps as objects making it easier to customize parts of the map (legends, credits, and the layout) and to build multi-map figures.

hepa_map <- tm_shape(hepa_sf) + 
  tm_fill('Rate.per.100000', 
style='quantile', 
palette='BuPu', 
title='HepA Rates \n Per 100,000 people') + 
  tm_borders(alpha=0.7) + 
  tm_credits('Quantile (Equal-Frequency) Class Intervals', 
             position=c('RIGHT', 'BOTTOM')) + 
  tm_layout(main.title="HepA Rates in the US",
            inner.margins = c(0.1, 0.1, 0.05, 0.05), 
            main.title.size=1.2, legend.title.size=0.5,  
            legend.text.size=0.5)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_fill()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## β„Ή Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
hepa_map
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "BuPu" is named
## "brewer.bu_pu"
## Multiple palettes called "bu_pu" found: "brewer.bu_pu", "matplotlib.bu_pu". The first one, "brewer.bu_pu", is returned.

Next, map standardized log poverty.

poverty_map <- tm_shape(hepa_sf) + 
  tm_fill('poverty_z', 
style='quantile', 
palette='BuPu', 
title='Standardized Log \n Poverty (%)') + # "\n" moves text to the next line 
  tm_borders(alpha=0.2) + 
  tm_credits('Quantile (Equal-Frequency) Class Intervals', 
             position=c('RIGHT', 'BOTTOM')) + 
  tm_layout(main.title="% of Households living below the federal poverty level in the US",
            inner.margins = c(0.1, 0.1, 0.05, 0.05), 
            main.title.size=1.2, legend.title.size=0.5,  
            legend.text.size=0.5)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_fill()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## β„Ή Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
poverty_map
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "BuPu" is named
## "brewer.bu_pu"
## Multiple palettes called "bu_pu" found: "brewer.bu_pu", "matplotlib.bu_pu". The first one, "brewer.bu_pu", is returned.
## 
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## β„Ή Set the tmap option `component.autoscale = FALSE` to disable rescaling.

Plot standardized log educational attainment.

edu_map <- tm_shape(hepa_sf) + 
  tm_fill('edu_z', 
style='quantile', 
palette='BuPu', 
title='Standardized Log \n Low Educational Attainment') + #reverse coding?? 
  tm_borders(alpha=0.2) + 
  tm_credits('Quantile (Equal-Frequency) Class Intervals', 
             position=c('RIGHT', 'BOTTOM')) + 
  tm_layout(main.title="% of Population 25 years and older w/o HS diploma",
            inner.margins = c(0.1, 0.1, 0.05, 0.05), 
            main.title.size=1.2, legend.title.size=0.5,  
            legend.text.size=0.5)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_fill()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## β„Ή Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
edu_map
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "BuPu" is named
## "brewer.bu_pu"
## Multiple palettes called "bu_pu" found: "brewer.bu_pu", "matplotlib.bu_pu". The first one, "brewer.bu_pu", is returned.

###Spatial Exploratory Analysis

Defining Spatial Proximity

library (spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
#Drop missing data
hepa_sf = subset (hepa_sf, !is.na (hepa_sf$Rate.per.100000))

#Because our shapefile here is un-projected, distance operations can be difficult. Alternatively, we can project the data to, for example, the Lambert Conformal Conic North America.
# 1. Project the data to Lambert Conformal Conic North America (ESRI:102004)
# (Recommended for accurate distance and spatial operations across the US)
hepa_sf_proj = st_transform(hepa_sf, crs = "ESRI:102004")
center_proj = st_centroid(hepa_sf_proj)
#Create 1-st order adjacency and distance buffer (100km)
#Plot the projected geometry
#plot (st_geometry(hepa_sf_proj))
nb.1st = poly2nb (hepa_sf_proj)
nb.buffer = dnearneigh(center_proj, d1=0, d2 = 25*1000)
#nb2 = poly2nb (hepa_sf_proj, queen=FALSE) 
#take the union of neighbors
nb = union.nb (nb.1st, nb.buffer) #take the union of neighbors


#A lot of good summaries about the neighbor structure
summary (nb)
## Neighbour list object:
## Number of regions: 51 
## Number of nonzero links: 218 
## Percentage nonzero weights: 8.381392 
## Average number of links: 4.27451 
## 2 regions with no links:
## 50, 51
## 3 disjoint connected subgraphs
## Link number distribution:
## 
##  0  1  2  3  4  5  6  7  8 
##  2  1  5  9 10 10 10  2  2 
## 1 least connected region:
## 37 with 1 link
## 2 most connected regions:
## 6 18 with 8 links

The function nb2mat () then creates a proximity/weight matrix

#Create weight matrix
W= nb2listw(nb, style="W", zero.policy = TRUE)
#W = nb2mat (nb, style = "B", zero.policy = TRUE) #Stopped code - Hawaii (or states separated by water/projected geometry boundaries) do not share a physical border with any other state SO add zero.policy

W
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 51 
## Number of nonzero links: 218 
## Percentage nonzero weights: 8.381392 
## Average number of links: 4.27451 
## 2 regions with no links:
## 50, 51
## 3 disjoint connected subgraphs
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1       S2
## W 49 2401 49 24.57728 202.0046
#center_proj = st_centroid(hepa_sf_proj) ##Extract centroids of polygons
#center = st_transform(center_proj, crs = 4326) #Re-project back to WGS 84
#center = st_coordinates (center)                 

#plot (st_geometry(hepa_sf))
#plot (nb, center, add = TRUE, col = "blue")
#plot (nb2, center, add = TRUE, col = "red", lwd = 1)
#title (main="Blue lines = additional neighbours by queen's case")

## Using kth-nearest distance
#nb.k1 = knn2nb (knearneigh(center_proj,k=1))
#nb.k2 = knn2nb (knearneigh(center_proj,k=2), row.names=row.names(center))

#plot (st_geometry(hepa_sf))
#plot (nb.k2, center, add = TRUE, col = "red", lwd = 1)
#plot (nb.k1, center, add = TRUE, col = "blue", lwd = 1)
#title (main="Blue lines = additional second nearest-neighbour")

## Using buffer distance
#nb.buffer1 = dnearneigh(center_proj, d1=0, d2 = 25*1000)
#nb.buffer2 = dnearneigh(center_proj, d1=0, d2 = 40*1000)

#plot (st_geometry(hepa_sf))
#plot (nb.buffer2, center, add = TRUE, col = "red", lwd = 1)
#plot (nb.buffer1, center, add = TRUE, col = "blue", lwd =2)
#legend ("topright", legend =c("<25 km", "< 40 km"), col= c("blue", "red"), pch = 16)

Monte Carlo-based Moran’s I

moran.mc(hepa_sf$Rate.per.100000, W, nsim = 10000)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  hepa_sf$Rate.per.100000 
## weights: W  
## number of simulations + 1: 10001 
## 
## statistic = 0.19811, observed rank = 9779, p-value = 0.0222
## alternative hypothesis: greater
moran.mc(hepa_sf$poverty_z, W, nsim = 10000) #look at the p value
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  hepa_sf$poverty_z 
## weights: W  
## number of simulations + 1: 10001 
## 
## statistic = 0.50058, observed rank = 10001, p-value = 9.999e-05
## alternative hypothesis: greater
moran.mc(hepa_sf$edu_z, W, nsim = 10000) #look at the p value
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  hepa_sf$edu_z 
## weights: W  
## number of simulations + 1: 10001 
## 
## statistic = 0.52672, observed rank = 10001, p-value = 9.999e-05
## alternative hypothesis: greater

###Spatial Regression Modeling

#For state 𝑠, let π‘Œπ‘  be the number of Hepatitis A cases, 𝑋𝑠 be poverty or low educational attainment, and 𝑃𝑠 be the population size. We consider several Poisson log-linear models of the following form.

#π‘Œπ‘ log(πœ‡π‘ )∼Poisson(πœ‡π‘ )=log(𝑃𝑠)+𝛽0+𝛽1𝑋𝑠+πœƒπ‘ 
#The standard Poisson model assumes no spatial variation (πœƒπ‘ =0 for all 𝑠).

#The independent random intercept model assumes that πœƒπ‘ βˆΌπ‘(0,𝜏2).

#The spatial random intercept model (under the improper CAR model formulation) assumes that πœƒπ‘ =πœ‚π‘ +𝛾𝑠, where πœ‚π‘ βˆΌπ‘(0,𝜏2) and 𝛾𝑠 is a conditional autoregressive model.

library (INLA)
## Warning: package 'INLA' was built under R version 4.5.3
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
# 1. Strip commas and convert population to numeric
hepa_sf <- hepa_sf %>%
  mutate(
    Population = as.numeric(str_remove_all(Population, ",")),
    log_pop = log(Population)
  )

#Some data wrangling first
hepa_sf$areal_ID = 1:nrow (hepa_sf) #Create a county ID from 1 to 226
hepa_sf$log_pop = log (hepa_sf$Population)
hepa_sf$poverty_std = hepa_sf$poverty_p / 0.1 #standardized by the interquartile range (0.1)

#INLA needs a specific summary of adjacency structure
nb2INLA ("adj.txt", nb) #write out a file
G <- inla.read.graph(filename = "adj.txt") #read in to get INLA's graph format

#Fit glm
fit.glm = glm (Cases~poverty_std+offset (log_pop), dat = hepa_sf, family = "poisson")

#Fit exchangeable (random effect model)
#Note: "E" here is the offset and it does not need to be log-transformed
fit.exch = inla (Cases~1+f(areal_ID)+poverty_std, E= hepa_sf$Population, family = "poisson", data = hepa_sf,
                 control.compute = list(dic = TRUE, waic = TRUE, return.marginals.predictor=TRUE))

#Fit improper CAR model (Besag)
fit.iCAR = inla (Cases~1+f(areal_ID,  model = "besag", graph = G)+poverty_std, E = hepa_sf$Population, family = "poisson", data = hepa_sf,control.compute = list(dic = TRUE, waic = TRUE))


library(knitr)

# 1. Extract GLM results
glm_sum <- summary(fit.glm)$coefficients["poverty_std", ]
glm_beta <- glm_sum["Estimate"]
glm_se   <- glm_sum["Std. Error"]
glm_rr   <- exp(glm_beta)
glm_lower <- exp(glm_beta - 1.96 * glm_se)
glm_upper <- exp(glm_beta + 1.96 * glm_se)

# 2. Extract Independent Random Effect (Exchangeable) results
exch_beta  <- fit.exch$summary.fixed["poverty_std", "mean"]
exch_se    <- fit.exch$summary.fixed["poverty_std", "sd"]
exch_rr    <- exp(exch_beta)
exch_lower <- exp(fit.exch$summary.fixed["poverty_std", "0.025quant"])
exch_upper <- exp(fit.exch$summary.fixed["poverty_std", "0.975quant"])

# 3. Extract Spatial (iCAR) results
icar_beta  <- fit.iCAR$summary.fixed["poverty_std", "mean"]
icar_se    <- fit.iCAR$summary.fixed["poverty_std", "sd"]
icar_rr    <- exp(icar_beta)
icar_lower <- exp(fit.iCAR$summary.fixed["poverty_std", "0.025quant"])
icar_upper <- exp(fit.iCAR$summary.fixed["poverty_std", "0.975quant"])

# 4. Combine into a clean Data Frame
results_table <- data.frame(
  Model = c("GLM", "Ind Random Effect", "Spatial"),
  beta = round(c(glm_beta, exch_beta, icar_beta), 2),
  `Std Err` = round(c(glm_se, exch_se, icar_se), 2),
  RR = round(c(glm_rr, exch_rr, icar_rr), 2),
  `95% Interval` = sprintf("(%.2f, %.2f)", 
                           c(glm_lower, exch_lower, icar_lower), 
                           c(glm_upper, exch_upper, icar_upper)),
  check.names = FALSE
)

# Display formatted table in R Console
kable(results_table, align = "lcccc")
Model beta Std Err RR 95% Interval
GLM 0.02 0.00 1.02 (1.02, 1.02)
Ind Random Effect 0.02 0.01 1.02 (1.01, 1.04)
Spatial 0.01 0.01 1.01 (1.00, 1.03)