CP Mini 3
Amritha Karthikeyan
2025-10-02
Cleaning up data
Import the hospital POI data (link) and make sure it is tidy and ready for analysis.
library(sf)## Warning: package 'sf' was built under R version 4.3.3## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUElibrary(dplyr)## Warning: package 'dplyr' was built under R version 4.3.1##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)## Warning: package 'tidyr' was built under R version 4.3.1library(stringr)## Warning: package 'stringr' was built under R version 4.3.1library(tigris)## Warning: package 'tigris' was built under R version 4.3.3## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.library(tidycensus)## Warning: package 'tidycensus' was built under R version 4.3.3library(tmap)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.3Cleaning up data
Import the hospital POI data (link) and make sure it is tidy and ready for analysis.
Variables chosen
Median Household Income (ACS: B19013)
I chose this as one of my variables because median household income can often determine the kind of resources a family has access to, like the ability to afford health insurance and transit to hospitals.
Insurance Coverage B27010_017
I chose this as another one of my variables because its an indicator of whether or not people have direct access to care. The number of hospitals in uninsured areas versus insured areas can also highlight the inequities in health access.
#Download ACS data
# get counties
ga_counties <- counties(state = "GA", cb = TRUE, year = 2023) %>%
st_transform(st_crs(hosp))## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%hosp_with_county <- st_join(hosp, ga_counties["NAME"])
counties_list <- unique(hosp_with_county$NAME)
#download acs data
vars <- c(
median_income = "B19013_001",
uninsured = "B27010_017",
uninsured_base = "B27010_001",
total_pop = "B01001_001"
)
acs <- get_acs(
geography = "tract",
state = "GA",
county = counties_list,
variables = vars,
year = 2023,
survey = "acs5",
geometry = TRUE
) %>%
select(GEOID, variable, estimate, geometry) %>%
tidyr::spread(variable, estimate)## Getting data from the 2019-2023 5-year ACS## Warning: • You have not set a Census API key. Users without a key are limited to 500
## queries per day and may experience performance limitations.
## ℹ For best results, get a Census API key at
## http://api.census.gov/data/key_signup.html and then supply the key to the
## `census_api_key()` function to use it throughout your tidycensus session.
## This warning is displayed once per session.## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.## | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 20% | |============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%#Tidy up data
#checking for acs na data
acs_na_check <- acs %>%
st_set_geometry(NULL) %>%
summarize(
na_income = sum(is.na(median_income)),
na_uninsured = sum(is.na(uninsured)),
na_unins_base = sum(is.na(uninsured_base)),
na_total_pop = sum(is.na(total_pop))
)
#acs_na_check
#join acs and poi data
#same projected CRS
projected_crs <- 5070 # NAD83 / Albers Equal Area
acs_proj <- st_transform(acs, projected_crs)
hosp_proj <- st_transform(hosp, projected_crs)
hosp_tract <- st_join(hosp_proj, acs_proj %>% select(GEOID))
hosp_counts <- hosp_tract %>%
st_set_geometry(NULL) %>%
count(GEOID, name = "hosp_count")
acs_joined <- acs_proj %>%
left_join(hosp_counts, by = "GEOID") %>%
mutate(hosp_count = ifelse(is.na(hosp_count), 0, hosp_count))
acs_joined <- acs_joined %>%
mutate(
hosp_per_100k = ifelse(total_pop > 0, hosp_count / total_pop * 1e5, NA_real_),
hosp_density = hosp_count / (as.numeric(st_area(geometry)) / 1e6) # per km²
)
#acs_joined#Define your operational measure of healthcare access How many hospitals are within X miles of each tract?: This is the most suited measure because it determines how accessible hospitals are for people in these census tracts. Now I can compare the amount of hospitals in wealthier neighborhoods to amount of hospitals in lower-income neighborhoods.
#compute
tract_pts <- st_centroid(acs_proj)## Warning: st_centroid assumes attributes are constant over geometriesbuffers <- st_buffer(tract_pts, dist = 5*1609.34)
hosp_within_buffer <- st_join(hosp_proj, buffers, join = st_within)
hosp_counts_buffer <- hosp_within_buffer %>%
st_set_geometry(NULL) %>%
count(GEOID, name = "hosp_within_5mi")
acs_access <- tract_pts %>%
left_join(hosp_counts_buffer, by = "GEOID") %>%
mutate(hosp_within_5mi = ifelse(is.na(hosp_within_5mi), 0, hosp_within_5mi))
#acs_access# Hospitals and Median Household Income with hospital count
tm_shape(acs_proj) +
tm_fill("median_income",
style = "fixed",
breaks = c(0,
quantile(acs_proj$median_income, 0.2, na.rm = TRUE),
quantile(acs_proj$median_income, 0.8, na.rm = TRUE),
max(acs_proj$median_income, na.rm = TRUE)),
palette = c("pink", "white", "lightblue"),
title = "Low vs High Household Income") +
tm_shape(hosp_proj) +
tm_dots(size = 0.3, col = "black", legend.show = TRUE) +
tm_add_legend(type = "symbol",
labels = "Hospital",
col = "black",
size = 0.95,
title = "Hospitals") +
tm_title("Hospital Count and Median Household Income") ## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_fill()`: instead of `style = "fixed"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'breaks', '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>)'
## [tm_dots()] Argument `legend.show` unknown.
## [v3->v4] `tm_add_legend()`: use `type = "symbols"` instead of `type =
## "symbol"`.
## [v3->v4] `tm_add_legend()`: use `fill` instead of `col` for the fill color of
## symbols.
## [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.
tmap_mode("view")## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`Looking at this plot there seems to be more hospitals in the middle income to high income households compared to the low income.
acs_access$income_group <- cut(acs_access$median_income,
breaks = quantile(acs_access$median_income, probs = seq(0,1,0.25), na.rm=TRUE),
include.lowest = TRUE,
labels = c("Low", "Lower-Mid", "Upper-Mid", "High"))
ggplot(
acs_access %>% filter(!is.na(income_group), !is.na(hosp_within_5mi)),
aes(x = income_group, y = hosp_within_5mi, fill = income_group)
) +
geom_boxplot() +
labs(
title = "Hospital Access by Income Quartile",
x = "Income Group",
y = "Hospitals within 5 miles"
) +
theme_minimal()
It seems like the low, lower-mid, and upper-mid all have similar medians when it comes to hospital axis while the high income households have the higher median which means they are in closer proximity to hospitals.
Let’s look at the amount of insured people in these areas.
acs_proj <- acs_proj %>%
mutate(insured_rate = 1 - (uninsured / uninsured_base)) %>%
mutate(insured_group2 = ifelse(insured_rate < median(insured_rate, na.rm = TRUE),
"Low coverage", "High coverage"))
tmap_mode("view")## ℹ tmap modes "plot" - "view"tm_shape(acs_proj) +
tm_fill(
col = "insured_group2",
palette = c("white", "lightblue"),
title = "Insurance Coverage",
border.col = NA, # turn borders off
lwd = 0
) +
tm_shape(hosp_proj) +
tm_dots(col = "black", size = 0.30) +
tm_title("Hospitals and Insurance Coverage (Low vs High)")##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'Registered S3 method overwritten by 'jsonify':
## method from
## print.json jsonliteacs_access <- acs_access %>%
mutate(insured_rate = 1 - (uninsured / uninsured_base))
acs_access <- acs_access %>%
mutate(insured_group = ntile(insured_rate, 4)) %>%
mutate(insured_group = factor(insured_group,
labels = c("Low", "Lower-Mid", "Upper-Mid", "High")))
library(ggplot2)
ggplot(
acs_access %>% filter(!is.na(insured_rate), !is.na(hosp_within_5mi)),
aes(x = insured_group, y = hosp_within_5mi, fill = insured_group)
) +
geom_boxplot() +
scale_fill_brewer(palette = "Blues") +
labs(
title = "Hospital Access by Insurance Coverage Quartile",
x = "Insurance Coverage Group",
y = "Hospitals within 5 miles"
) +
theme_minimal()
Similar to the incom chart, the median for hospital access by insurance
coverage for the High insurance group is far higher than the rest of the
categories. # Regression
acs_model <- acs_proj %>%
st_set_geometry(NULL) %>%
mutate(insured_rate = 1 - (uninsured / uninsured_base)) %>%
left_join(
acs_access %>%
st_set_geometry(NULL) %>%
select(GEOID, hosp_within_5mi),
by = "GEOID"
)
model2 <- glm(
hosp_within_5mi ~ median_income + insured_rate,
data = acs_model,
family = "poisson"
)
summary(model2)##
## Call:
## glm(formula = hosp_within_5mi ~ median_income + insured_rate,
## family = "poisson", data = acs_model)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.962e-03 5.382e-01 0.006 0.99561
## median_income 1.102e-06 2.828e-07 3.896 9.79e-05 ***
## insured_rate 1.529e+00 5.539e-01 2.761 0.00576 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4337.4 on 1253 degrees of freedom
## Residual deviance: 4307.1 on 1251 degrees of freedom
## (16 observations deleted due to missingness)
## AIC: 8059.6
##
## Number of Fisher Scoring iterations: 5These results show that for every $10,000 increase in income, the expected hospital access increases by about 1.1%.
Tracts with higher insurance coverage have 4.6x more hospitals within 5 miles than those with lower rates of insurance coverage.
This regression model shows that both wealthier and better-insured areas have more hospitals nearby than those that are poorer and worse-insured areas. Hospital access is indeed inequitable across the Metro Atlanta area.