minihw3
2025-10-01
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── 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 errorslibrary(dplyr)
library(sf)## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUElibrary(tmap)
library(leaflet)
library(tidycensus)
library(magrittr)##
## Attaching package: 'magrittr'
##
## The following object is masked from 'package:purrr':
##
## set_names
##
## The following object is masked from 'package:tidyr':
##
## extractlibrary(skimr)
hosp <- st_read('https://raw.githubusercontent.com/ujhwang/urban-analytics-2025/main/Assignment/mini_3/hospital_11counties.geojson')## Reading layer `hospital_11counties' from data source
## `https://raw.githubusercontent.com/ujhwang/urban-analytics-2025/main/Assignment/mini_3/hospital_11counties.geojson'
## using driver `GeoJSON'
## Simple feature collection with 119 features and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -84.73147 ymin: 33.42719 xmax: -83.92052 ymax: 34.24585
## Geodetic CRS: WGS 84hosp <- hosp %>%
select(-c(id, primary_type, types, status))
#It looks quite tidy to me. Unsure about the ID Column, and I'll wait for the types Column
api <- read_lines('senses.txt')
census_api_key(api)## To install your API key for use in future sessions, run this function with `install = TRUE`.acs5_variables <- load_variables(year = 2023, dataset = "acs5", cache = TRUE)
acs5_variables <- acs5_variables %>%
filter(geography == 'tract')
View(acs5_variables)
# Since we're examining equity, I chose race, median income, median age by sex, and type of health insurance (medicaid only). These values I believe are going to be a good representative for seeing if hospitals have some sort of implicit bias within the MSA. Are better rated hospitals located in higher income areas? Are worse rated hospitals rated in medicaid majoirty places? Race correlates with income in the United States. I chose median age by sex to see if they were located where old people are since they need the most healthcare.
county <- get_acs(geography = 'tract', state = 'GA', county = c('Cherokee','Clayton','Cobb','DeKalb','Douglas','Fayette','Forsyth','Fulton','Gwinnett','Henry','Rockdale'), variables = c('B02001_002','B07011_001','B01002_001',"B27010_023", 'B01003_001'), geometry = TRUE)## Getting data from the 2019-2023 5-year ACS
## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |=================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%geo <- county %>%
select(GEOID, geometry)
county <- county %>%
sf::st_drop_geometry()
#pivoting wider to have the variables more visible
county <- county %>%
select(GEOID, NAME, variable, estimate) %>%
pivot_wider(names_from = variable, values_from = estimate) %>%
rename('Race' = B02001_002,
'medinc' = B07011_001,
'sex/age'= B01002_001,
'medicaid' = B27010_023,
'totpop' = B01003_001)
county <- county %>%
mutate(nonwhite = totpop - Race,
whitepop = Race) %>%
select(-Race)
geo <- geo %>%
distinct(GEOID,geometry)
sfcounty <- county %>%
left_join(geo, by = 'GEOID') %>%
sf::st_as_sf()#Checking and setting up for mapping
library(tmap)
library(tigris)## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.install.packages('skimr')## Warning: package 'skimr' is in use and will not be installedlibrary(skimr)
#tigris county lines
counties <- counties(stat = 'GA', year = 2020)## | | | 0% | |======================================================================| 100%## HTTP download failed, trying FTP as fallback...#where are the hospitals are they correct?
map1 <- tmap_mode('plot') +
tm_shape(hosp) +
tm_dots(fill = 'red') +
tm_shape(counties) +
tm_borders(col = 'black', lwd = 2)## ℹ tmap modes "plot" - "view"## ℹ toggle with `tmap::ttm()`
## This message is displayed once per session.print(map1)
#what does my sfcounty df look like?
map2 <- tmap_mode('plot') +
tm_borders () +
tm_shape(sfcounty) +
tm_fill('medinc')## ℹ tmap modes "plot" - "view"print(map2)
#Spatial join time.
st_crs(sfcounty)## Coordinate Reference System:
## User input: NAD83
## wkt:
## GEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]]st_crs(hosp)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]hosp <- st_transform(hosp, st_crs(sfcounty))
tract_hosp <- st_join(sfcounty, hosp, left = TRUE)
#lots of NAs since not one in every county, we remove.
map3 <- tmap_mode('plot') +
tm_borders() +
tm_shape(tract_hosp) +
tm_fill('medinc') +
tm_shape(hosp) +
tm_dots(col = 'red', size = .2, alpha = .8) +
tm_shape(counties) +
tm_borders(col = 'black', lwd = 2)## ℹ tmap modes "plot" - "view"
##
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
##
## [v3->v4] `tm_dots()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').
## [v3->v4] `tm_dots()`: use `fill_alpha` instead of `alpha`.print(map3)
Based on this Pew Research Article, average hospital access in cities is about 4.4 miles for urban areas. Atlanta is definitely an urban area, but I realize that some of these counties are suburban. The nearest hospital in suburban areas is 5.6 miles away. I’m keeping with the urban area 4.4, going to make a buffer of each dot and show how much of Atlanta is more than 4.4 miles away from a hospital. This could show a lack of average care in the US.
#Creating a Buffer
#first need to transform hospital back into projected
hosp_proj <- st_transform(hosp, 26916)
buffer_meters <- 4.4 * 1609.34
hosp_buffers <- st_buffer(hosp_proj, dist = buffer_meters)
# back to CRS
hosp_buffers <- st_transform(hosp_buffers, st_crs(hosp))
#median income shows that theres a lot of overlap of care along midtown through buckhead. it shows some sparsity at the edges.
map4 <- tmap_mode('plot') +
tm_shape(tract_hosp) +
tm_fill('medinc') +
tm_shape(counties) +
tm_borders(col = 'black', lwd = 1) +
tm_shape(hosp) +
tm_dots('red', size = .3, fill_alpha = .5) +
tm_shape(hosp_buffers) +
tm_borders('red', lwd = 1, lty = 'dashed')## ℹ tmap modes "plot" - "view"print(map4)
#i was interested to see the medicaid map and it shows some pockets of no care...
map5 <- tmap_mode('plot') +
tm_shape(tract_hosp) +
tm_fill('medicaid') +
tm_shape(counties) +
tm_borders(col = 'black', lwd = 1) +
tm_shape(hosp) +
tm_dots('red', size = .3, fill_alpha = .5) +
tm_shape(hosp_buffers) +
tm_borders('red', lwd = 1, lty = 'dashed')## ℹ tmap modes "plot" - "view"print(map5)
#turned that into percentages to see how that affected things
tract_hosp <- tract_hosp %>%
mutate('pctmed' = ((medicaid/totpop)*100))
map6 <- tmap_mode('plot') +
tm_shape(tract_hosp) +
tm_fill('pctmed') +
tm_shape(counties) +
tm_borders(col = 'black', lwd = 1) +
tm_shape(hosp_buffers) +
tm_borders('red',lwd = 1, lty = 'dashed')## ℹ tmap modes "plot" - "view"print(map6)
Based on these maps, we see that there are definitely more hospitals in places with higher median income, but that coverage even in places with less median income have access to care. When we take the percentage of those on medicaid, coverage is still within 4.4 miles of a hospital for most census tracts except for one, Tract 103.11. From a visual perspective, it would be hard to say that hospitals are not equitably placed.
#regression -- making a boolean
tract_hosp <- tract_hosp %>%
mutate('avail' = ifelse(is.na(name), 0, 1))
tract_hosp <- tract_hosp %>%
rename('age' = 'sex/age')
m1 <- glm(avail ~ medicaid + rating, data = tract_hosp, family = binomial)## Warning: glm.fit: algorithm did not convergeprint(summary(m1))##
## Call:
## glm(formula = avail ~ medicaid + rating, family = binomial, data = tract_hosp)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.657e+01 1.704e+05 0 1
## medicaid -1.051e-09 5.133e+02 0 1
## rating 2.100e-06 3.937e+04 0 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 0.0000e+00 on 117 degrees of freedom
## Residual deviance: 6.8459e-10 on 115 degrees of freedom
## (1155 observations deleted due to missingness)
## AIC: 6
##
## Number of Fisher Scoring iterations: 25#no significance between medicaid and hospital location
m2 <- glm(avail ~ medicaid + age + medinc, data = tract_hosp, family = binomial)
print(summary(m2))##
## Call:
## glm(formula = avail ~ medicaid + age + medinc, family = binomial,
## data = tract_hosp)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.420e+00 6.351e-01 -2.236 0.0253 *
## medicaid -1.482e-03 1.429e-03 -1.037 0.2998
## age -2.101e-02 1.630e-02 -1.289 0.1975
## medinc 2.449e-07 5.655e-06 0.043 0.9655
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 779.51 on 1262 degrees of freedom
## Residual deviance: 777.13 on 1259 degrees of freedom
## (10 observations deleted due to missingness)
## AIC: 785.13
##
## Number of Fisher Scoring iterations: 5#no significance when we add other factors either ....
m3 <- glm(avail ~ whitepop + nonwhite + age + medinc, data = tract_hosp, family = binomial)
print(summary(m3))##
## Call:
## glm(formula = avail ~ whitepop + nonwhite + age + medinc, family = binomial,
## data = tract_hosp)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.113e+00 7.506e-01 -2.815 0.00487 **
## whitepop 3.120e-04 7.319e-05 4.263 2.02e-05 ***
## nonwhite 1.104e-04 8.046e-05 1.373 0.16985
## age -2.306e-02 1.734e-02 -1.330 0.18341
## medinc -2.932e-06 6.289e-06 -0.466 0.64102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 779.51 on 1262 degrees of freedom
## Residual deviance: 761.14 on 1258 degrees of freedom
## (10 observations deleted due to missingness)
## AIC: 771.14
##
## Number of Fisher Scoring iterations: 5#but once we added if white people are present, then it becomes significant. there are a lot of factors for why this may be.
m4 <- glm(avail ~ rating + whitepop + nonwhite + age + medinc, data = tract_hosp, family = binomial)## Warning: glm.fit: algorithm did not convergeprint(summary(m4))##
## Call:
## glm(formula = avail ~ rating + whitepop + nonwhite + age + medinc,
## family = binomial, data = tract_hosp)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.657e+01 2.724e+05 0 1
## rating -3.341e-09 4.121e+04 0 1
## whitepop -5.269e-13 2.331e+01 0 1
## nonwhite 1.308e-13 2.843e+01 0 1
## age 3.180e-11 5.836e+03 0 1
## medinc -3.796e-14 2.266e+00 0 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 0.0000e+00 on 116 degrees of freedom
## Residual deviance: 6.7879e-10 on 111 degrees of freedom
## (1156 observations deleted due to missingness)
## AIC: 12
##
## Number of Fisher Scoring iterations: 25# adding in rating we see that no values are significant now.
#this shows that hopsitals may be placed equitably since only the white people variable showed significance #plotting
library(ggplot2)
tract_w_hosp <- tract_hosp %>%
filter(avail == 1)
plot1 <- ggplot(tract_hosp, aes(x = totpop, y = medinc)) +
geom_smooth(mapping = aes(x = totpop, y = medinc), method = 'lm') +
geom_point(aes(x = totpop, y = medinc), alpha = .6, color = 'darkblue') +
labs(
title = 'total population and median income in tracts with hospitals',
x = 'total population',
y = 'medinc'
) +
theme_minimal()
print(plot1)## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
# i realize that this actually isn't very helpful
plot2 <- ggplot(tract_w_hosp, aes(x = rating_count, y = rating)) +
geom_point(color = 'darkblue', alpha = .7) +
labs(
title = 'ratings of hospitals',
x = 'rating count',
y = 'rating'
) +
theme_minimal()
print(plot2)
#rating proved to be not very influential
The graphs didn’t prove to be useful since I wanted to look at other
factors. There were other things I could have explored and probably used
the regression lines. Based on the regressions and the maps, it would be
hard to say that hospitals are not equitably placed around the city.
There are some pockets where care is not as easily accessible, but
overall the spread and range of a hospital’s services (4.4 miles)
overlaps across most of the 11 county region in Atlanta. The regressions
showed that age, medicaid population, and median income weren’t
influential. Medicaid was likely not influential since it is an optional
question and therefore under-reported. White population was the only
influential factor with a value of .001 taking everything into account.
However, I added medicaid into the regression too and it made the
variable insignificant. More examination into that relationship should
be considered.