library(tidyverse)
library(tidyr)
library(leaflet) # interactive maps
library(scales) # axis label formatting
cities500 <- read_csv("500CitiesLocalHealthIndicators.cdc.csv")Healthy Cities GIS Assignment: Health Inequality Across Baltimore’s Census Tracts
Load the libraries and the data
This is the CDC’s 500 Cities: Local Data for Better Health dataset, which reports health indicator estimates for the 500 largest US cities down to the census tract level. It has 810,103 rows and 24 variables.
The GeoLocation variable has (lat, long) format
Split GeoLocation (lat, long) into two columns: lat and long.
latlong <- cities500 |>
mutate(GeoLocation = str_replace_all(GeoLocation, "[()]", "")) |>
separate(GeoLocation, into = c("lat", "long"), sep = ",", convert = TRUE)
head(latlong)# A tibble: 6 × 25
Year StateAbbr StateDesc CityName GeographicLevel DataSource Category
<dbl> <chr> <chr> <chr> <chr> <chr> <chr>
1 2017 CA California Hawthorne Census Tract BRFSS Health Outcom…
2 2017 CA California Hawthorne City BRFSS Unhealthy Beh…
3 2017 CA California Hayward City BRFSS Health Outcom…
4 2017 CA California Hayward City BRFSS Unhealthy Beh…
5 2017 CA California Hemet City BRFSS Prevention
6 2017 CA California Indio Census Tract BRFSS Health Outcom…
# ℹ 18 more variables: UniqueID <chr>, Measure <chr>, Data_Value_Unit <chr>,
# DataValueTypeID <chr>, Data_Value_Type <chr>, Data_Value <dbl>,
# Low_Confidence_Limit <dbl>, High_Confidence_Limit <dbl>,
# Data_Value_Footnote_Symbol <chr>, Data_Value_Footnote <chr>,
# PopulationCount <dbl>, lat <dbl>, long <dbl>, CategoryID <chr>,
# MeasureId <chr>, CityFIPS <dbl>, TractFIPS <dbl>, Short_Question_Text <chr>Filter the dataset (walkthrough example)
Remove the StateDesc that includes the United States, select Prevention as the category, filter for crude prevalence only, and select 2017.
latlong_clean <- latlong |>
filter(StateDesc != "United States") |>
filter(Data_Value_Type == "Crude prevalence") |>
filter(Year == 2017) |>
filter(StateAbbr == "CT") |>
filter(Category == "Unhealthy Behaviors")
head(latlong_clean)# A tibble: 6 × 25
Year StateAbbr StateDesc CityName GeographicLevel DataSource Category
<dbl> <chr> <chr> <chr> <chr> <chr> <chr>
1 2017 CT Connecticut Bridgeport Census Tract BRFSS Unhealthy B…
2 2017 CT Connecticut Danbury City BRFSS Unhealthy B…
3 2017 CT Connecticut Norwalk Census Tract BRFSS Unhealthy B…
4 2017 CT Connecticut Bridgeport Census Tract BRFSS Unhealthy B…
5 2017 CT Connecticut Hartford Census Tract BRFSS Unhealthy B…
6 2017 CT Connecticut Waterbury Census Tract BRFSS Unhealthy B…
# ℹ 18 more variables: UniqueID <chr>, Measure <chr>, Data_Value_Unit <chr>,
# DataValueTypeID <chr>, Data_Value_Type <chr>, Data_Value <dbl>,
# Low_Confidence_Limit <dbl>, High_Confidence_Limit <dbl>,
# Data_Value_Footnote_Symbol <chr>, Data_Value_Footnote <chr>,
# PopulationCount <dbl>, lat <dbl>, long <dbl>, CategoryID <chr>,
# MeasureId <chr>, CityFIPS <dbl>, TractFIPS <dbl>, Short_Question_Text <chr>What variables are included?
names(latlong_clean) [1] "Year" "StateAbbr"
[3] "StateDesc" "CityName"
[5] "GeographicLevel" "DataSource"
[7] "Category" "UniqueID"
[9] "Measure" "Data_Value_Unit"
[11] "DataValueTypeID" "Data_Value_Type"
[13] "Data_Value" "Low_Confidence_Limit"
[15] "High_Confidence_Limit" "Data_Value_Footnote_Symbol"
[17] "Data_Value_Footnote" "PopulationCount"
[19] "lat" "long"
[21] "CategoryID" "MeasureId"
[23] "CityFIPS" "TractFIPS"
[25] "Short_Question_Text" Remove the variables that will not be used
latlong_clean2 <- latlong_clean |>
select(-DataSource, -Data_Value_Unit, -DataValueTypeID, -Low_Confidence_Limit,
-High_Confidence_Limit, -Data_Value_Footnote_Symbol, -Data_Value_Footnote)
head(latlong_clean2)# A tibble: 6 × 18
Year StateAbbr StateDesc CityName GeographicLevel Category UniqueID Measure
<dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 2017 CT Connecticut Bridgep… Census Tract Unhealt… 0908000… Obesit…
2 2017 CT Connecticut Danbury City Unhealt… 918430 Obesit…
3 2017 CT Connecticut Norwalk Census Tract Unhealt… 0955990… Obesit…
4 2017 CT Connecticut Bridgep… Census Tract Unhealt… 0908000… Curren…
5 2017 CT Connecticut Hartford Census Tract Unhealt… 0937000… Obesit…
6 2017 CT Connecticut Waterbu… Census Tract Unhealt… 0980000… Obesit…
# ℹ 10 more variables: Data_Value_Type <chr>, Data_Value <dbl>,
# PopulationCount <dbl>, lat <dbl>, long <dbl>, CategoryID <chr>,
# MeasureId <chr>, CityFIPS <dbl>, TractFIPS <dbl>, Short_Question_Text <chr>My Assignment
Step 1: Filter to my own subset (under 900 observations)
I chose to look inside a single city rather than across many. Baltimore is one of the 500 Cities with tract-level coverage, which lets me ask whether health outcomes vary within one city as much as they do between cities.
My filter keeps four related indicators for Baltimore’s census tracts in 2017: obesity, physical inactivity, diabetes, and current smoking.
# Filter to Baltimore census tracts, 2017, crude prevalence, four chosen measures.
# 200 tracts x 4 measures = 800 observations, which is under the 900 limit.
baltimore <- cities500 |>
filter(StateAbbr == "MD",
CityName == "Baltimore",
GeographicLevel == "Census Tract",
Year == 2017,
Data_Value_Type == "Crude prevalence",
Short_Question_Text %in% c("Obesity", "Physical Inactivity",
"Diabetes", "Current Smoking"))
# Confirm the subset is under 900 observations.
nrow(baltimore)[1] 800# Split GeoLocation into lat and long, strip the commas out of PopulationCount so
# it becomes numeric, then keep only the columns needed.
baltimore_tidy <- baltimore |>
mutate(GeoLocation = str_replace_all(GeoLocation, "[()]", "")) |>
separate(GeoLocation, into = c("lat", "long"), sep = ",", convert = TRUE) |>
mutate(population = as.numeric(str_remove_all(PopulationCount, ","))) |>
select(TractFIPS, lat, long, population, Short_Question_Text, Data_Value)
# Pivot wider so each row is one census tract with all four measures side by side.
# This is what makes it possible to compare measures against each other.
baltimore_wide <- baltimore_tidy |>
pivot_wider(names_from = Short_Question_Text, values_from = Data_Value) |>
rename(obesity = Obesity,
inactivity = `Physical Inactivity`,
diabetes = Diabetes,
smoking = `Current Smoking`) |>
# One tract has suppressed estimates, so drop rows with missing values.
filter(!is.na(obesity), !is.na(inactivity),
!is.na(diabetes), !is.na(smoking), !is.na(lat))
# 199 usable tracts covering roughly 621,000 residents.
nrow(baltimore_wide)[1] 199summary(baltimore_wide$diabetes) Min. 1st Qu. Median Mean 3rd Qu. Max.
2.70 9.55 13.00 13.00 17.05 25.90 Diabetes prevalence runs from 2.7% to 25.9% across tracts in the same city — nearly a tenfold gap.
Step 2: A non-map plot
# Scatterplot of two indicators, with a third (obesity) carried by a colour
# gradient and a fourth (tract population) by point size. Following the Japan
# earthquake tutorial's use of scale_color_gradient for a continuous variable.
ggplot(baltimore_wide, aes(x = inactivity, y = diabetes)) +
geom_point(aes(color = obesity, size = population), alpha = 0.8) +
geom_smooth(method = "lm", se = FALSE, color = "grey30", lty = 2, linewidth = 0.6) +
scale_color_gradient(name = "Obesity (%)", low = "#FFD166", high = "#7A0403") +
scale_size_continuous(name = "Tract population", range = c(1.5, 6),
labels = label_comma()) +
labs(title = "Where Baltimoreans Move Less, They Also Have More Diabetes",
subtitle = "Each point is one census tract; the three indicators rise almost in lockstep",
x = "Adults Physically Inactive (%)",
y = "Adults with Diagnosed Diabetes (%)",
caption = "Source: CDC 500 Cities: Local Data for Better Health, 2017 crude prevalence estimates") +
theme_bw(base_size = 12)
# Quantify how tightly the indicators move together.
cor(baltimore_wide$inactivity, baltimore_wide$diabetes)[1] 0.9231879cor(baltimore_wide$obesity, baltimore_wide$diabetes)[1] 0.9313542Step 3: A first map of the subset
# Baltimore's coordinates, used to centre the map.
baltimore_lat <- 39.2904
baltimore_lon <- -76.6122
# A first leaflet map: one circle per census tract, with radius scaled to
# diabetes prevalence so the harder-hit tracts stand out.
leaflet() |>
setView(lng = baltimore_lon, lat = baltimore_lat, zoom = 12) |>
addProviderTiles("Esri.WorldStreetMap") |>
addCircles(data = baltimore_wide,
lng = ~long, lat = ~lat,
radius = ~diabetes * 12,
color = "#7A0403",
fillOpacity = 0.4)Step 4: Refined map with a mouse-click tooltip
# Build the popup text. <b> bolds a label and <br> forces a line break, so each
# tract's numbers appear on their own line when clicked.
popup_tract <- paste0(
"<b>Census Tract: </b>", baltimore_wide$TractFIPS, "<br>",
"<b>Population: </b>", comma(baltimore_wide$population), "<br>",
"<b>Diabetes: </b>", baltimore_wide$diabetes, "%<br>",
"<b>Obesity: </b>", baltimore_wide$obesity, "%<br>",
"<b>Physically inactive: </b>", baltimore_wide$inactivity, "%<br>",
"<b>Current smoking: </b>", baltimore_wide$smoking, "%"
)
# A continuous colour palette mapped to diabetes prevalence, so colour and not
# just size encodes the value. This also lets me add a legend.
pal <- colorNumeric(palette = c("#FFD166", "#F3722C", "#7A0403"),
domain = baltimore_wide$diabetes)# The refined map: colour encodes diabetes prevalence, radius reinforces it,
# a click popup gives every indicator for that tract, and a legend explains
# the colour scale.
leaflet(baltimore_wide) |>
setView(lng = baltimore_lon, lat = baltimore_lat, zoom = 12) |>
addProviderTiles("Esri.WorldStreetMap") |>
addCircleMarkers(
lng = ~long, lat = ~lat,
radius = ~diabetes * 0.7,
color = "#2B2B2B",
weight = 0.6,
fillColor = ~pal(diabetes),
fillOpacity = 0.8,
popup = popup_tract
) |>
addLegend(position = "bottomright",
pal = pal,
values = ~diabetes,
title = "Diabetes<br>prevalence (%)",
opacity = 0.9)Click any circle to see that tract’s population and all four health indicators.
Step 5: What the plots show
My subset covers 199 Baltimore census tracts, about 621,000 residents, with four 2017 indicators each. The scatterplot in Step 2 plots physical inactivity against diabetes, with obesity as a colour gradient and tract population as point size. The three indicators are almost interchangeable: inactivity correlates with diabetes at 0.92 and obesity with diabetes at 0.93. Tracts do not have one problem in isolation — the light yellow points sit together at the low end and the dark red ones cluster at the high end, so the same neighbourhoods carry all of it at once.
The map shows why that matters, because the pattern is geographic rather than scattered. Diabetes prevalence ranges from 2.7% to 25.9% within a single city, and the darkest circles form a band through East and West Baltimore while the lightest ones sit along the harbour and in the northern neighbourhoods. Clicking a tract makes the contrast concrete: neighbouring tracts a mile apart can differ by fifteen percentage points. What looks like a health statistic is really a map of where advantage is and is not, and a city-level average of about 13% diabetes would hide the entire story.
One honest limitation: these are modelled small-area estimates, not direct measurements. The CDC produces them by applying a statistical model to BRFSS survey responses combined with census demographics, so neighbouring tracts with similar demographics are partly predicted to look alike. The geographic pattern is real, but its sharpness is smoothed by the model that generated it.