Abstract
This analysis examines the spatial distribution of PM2.5 air
pollution exposure around JFK Airport and its relationship with
demographic characteristics, particularly minority populations. Using
buffer-based proximity assessment and spatial regression modeling, I
investigated whether environmental justice concerns exist in the study
area. Additionally, I conducted a Modifiable Areal Unit Problem (MAUP)
analysis to test the robustness of findings across multiple geographic
scales.
Key Findings: - High-minority tracts experience 0.24
µg/m³ higher PM2.5 exposure (p < 0.001) - Cohen’s d effect size: 0.74
(medium-large) - Strong spatial autocorrelation detected (Moran’s I =
0.926) - Spatial Error Model (SEM) provides best fit - MAUP sensitivity:
High (CV = 204.9%)
Package Loading and
Global Variables
library(sf)
library(spdep)
library(spatialreg)
library(tidycensus)
library(tigris)
library(tidyverse)
library(tmap)
library(biscale)
library(effsize)
library(coin)
library(ggplot2)
# Set data directory
data_dir <- "/Users/lily/Fall2025/data/"
# Define Coordinate Reference Systems
crs_wgs84 <- 4326
crs_projected_epsg <- 2263 # NAD83 / New York Long Island (ft)
nyc_counties <- c("New York", "Bronx", "Kings", "Queens", "Richmond")
Study Area
Definition
# JFK Airport coordinates (WGS84)
jfk_coords <- data.frame(name = "JFK Airport", lon = -73.7781, lat = 40.6413)
# Convert to spatial point and project
jfk_point <- st_as_sf(jfk_coords, coords = c("lon", "lat"), crs = crs_wgs84) %>%
st_transform(crs = crs_projected_epsg)
# 10-mile buffer around JFK Airport
buffer_distance <- 10 * 5280
study_area <- st_buffer(jfk_point, dist = buffer_distance)
cat("Study area defined:", buffer_distance/5280, "mile radius\n")
## Study area defined: 10 mile radius
Census Data
Acquisition
Census Tract
Demographics
# Get Census Tract Demographics and Geography
census_vars <- c(TotalPop = "B01003_001", MedIncome = "B19013_001")
nyc_tracts <- tidycensus::get_acs(
geography = "tract",
variables = census_vars,
state = "NY",
county = nyc_counties,
year = 2023,
survey = "acs5",
geometry = TRUE,
output = "wide"
) %>%
dplyr::select(GEOID, NAME, TotalPop = TotalPopE, MedIncome = MedIncomeE) %>%
sf::st_transform(crs = crs_projected_epsg) %>%
st_intersection(study_area)
## | | | 0% | |= | 1% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 12% | |========= | 13% | |========== | 14% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================ | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%
cat("Acquired", nrow(nyc_tracts), "census tracts\n")
## Acquired 1085 census tracts
Race and Ethnicity
Data
# Get Race/Ethnicity Data
minority_vars <- c(TotalPop_Race = "B03002_001", WhiteNH_Pop = "B03002_003")
nyc_race_data <- tidycensus::get_acs(
geography = "tract",
variables = minority_vars,
state = "NY",
county = nyc_counties,
year = 2023,
survey = "acs5",
geometry = FALSE,
output = "wide"
) %>%
dplyr::select(GEOID, TotalPop_Race = TotalPop_RaceE, WhiteNH_Pop = WhiteNH_PopE)
cat("Acquired race/ethnicity data\n")
## Acquired race/ethnicity data
Buffer-Based Exposure
Assignment
# Calculate distance from each tract centroid to JFK
nyc_tracts_centroids <- st_centroid(nyc_tracts)
nyc_tracts$Dist_to_JFK <- st_distance(nyc_tracts_centroids, jfk_point) %>% as.numeric()
# Assign exposure based on proximity zones
nyc_tracts <- nyc_tracts %>%
mutate(
PM25_Exposure_Category = case_when(
Dist_to_JFK <= 5280 * 1 ~ "High (0-1 mile)",
Dist_to_JFK <= 5280 * 3 ~ "Moderate (1-3 miles)",
Dist_to_JFK <= 5280 * 5 ~ "Background (3-5 miles)",
TRUE ~ "Reference (>5 miles)"
),
PM25_Estimate = case_when(
Dist_to_JFK <= 5280 * 1 ~ 8.5,
Dist_to_JFK <= 5280 * 3 ~ 7.8,
Dist_to_JFK <= 5280 * 5 ~ 7.2,
TRUE ~ 6.5
)
)
cat("PM2.5 exposure zones assigned\n")
## PM2.5 exposure zones assigned
Data Integration
nyc_sf_buffer <- nyc_tracts %>%
left_join(nyc_race_data, by = "GEOID") %>%
mutate(
Minority_Pop = TotalPop - WhiteNH_Pop,
Minority_Pct = (Minority_Pop / TotalPop) * 100
) %>%
filter(TotalPop > 0, !is.na(Minority_Pct), !is.na(PM25_Estimate)) %>%
select(GEOID, NAME, TotalPop, MedIncome, PM25_Estimate, PM25_Exposure_Category,
Minority_Pop, Minority_Pct, Dist_to_JFK, geometry)
cat("Final dataset:", nrow(nyc_sf_buffer), "census tracts,",
format(sum(nyc_sf_buffer$TotalPop), big.mark = ","), "total population\n")
## Final dataset: 1043 census tracts, 3,546,252 total population
Descriptive Statistics
by Exposure Zone
exposure_zone_summary <- nyc_sf_buffer %>%
st_drop_geometry() %>%
group_by(PM25_Exposure_Category) %>%
summarise(
N_Tracts = n(),
Total_Population = sum(TotalPop, na.rm = TRUE),
Mean_PM25 = mean(PM25_Estimate),
Mean_Minority_Pct = mean(Minority_Pct, na.rm = TRUE),
Median_Income = median(MedIncome, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(desc(Mean_PM25))
knitr::kable(exposure_zone_summary,
caption = "Summary Statistics by PM2.5 Exposure Zone",
digits = 2,
format.args = list(big.mark = ","))
Summary Statistics by PM2.5 Exposure Zone
| Moderate (1-3 miles) |
46 |
174,129 |
7.8 |
95.19 |
96,013.5 |
| Background (3-5 miles) |
166 |
518,923 |
7.2 |
91.85 |
84,590.0 |
| Reference (>5 miles) |
831 |
2,853,200 |
6.5 |
73.91 |
81,875.0 |
Interpretation: The table shows that areas closer to
JFK Airport have higher PM2.5 concentrations and higher minority
percentages. The moderate exposure zone (1-3 miles) has the highest
minority percentage (95.2%) and elevated PM2.5 levels (7.8 µg/m³).
Statistical Disparity
Testing
Define Minority
Groups
minority_threshold <- median(nyc_sf_buffer$Minority_Pct, na.rm = TRUE)
cat("Minority Percentage Median Threshold:", round(minority_threshold, 2), "%\n")
## Minority Percentage Median Threshold: 86.88 %
nyc_sf_buffer <- nyc_sf_buffer %>%
mutate(
Minority_Group = factor(
ifelse(Minority_Pct >= minority_threshold, "High Minority", "Low Minority"),
levels = c("Low Minority", "High Minority")
)
)
Descriptive
Statistics by Group
descriptive_stats_buffer <- nyc_sf_buffer %>%
st_drop_geometry() %>%
group_by(Minority_Group) %>%
summarise(
N_Tracts = n(),
Mean_PM25 = mean(PM25_Estimate, na.rm = TRUE),
SD_PM25 = sd(PM25_Estimate, na.rm = TRUE),
Median_PM25 = median(PM25_Estimate, na.rm = TRUE),
.groups = "drop"
)
knitr::kable(descriptive_stats_buffer,
caption = "PM2.5 Statistics by Minority Status",
digits = 3)
PM2.5 Statistics by Minority Status
| Low Minority |
521 |
6.546 |
0.187 |
6.5 |
| High Minority |
522 |
6.791 |
0.428 |
6.5 |
Welch’s T-Test
t_test_buffer <- t.test(PM25_Estimate ~ Minority_Group,
data = nyc_sf_buffer,
var.equal = FALSE)
print(t_test_buffer)
##
## Welch Two Sample t-test
##
## data: PM25_Estimate by Minority_Group
## t = -11.941, df = 712.53, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Low Minority and group High Minority is not equal to 0
## 95 percent confidence interval:
## -0.2845308 -0.2041802
## sample estimates:
## mean in group Low Minority mean in group High Minority
## 6.546449 6.790805
Result: High-minority tracts have significantly
higher PM2.5 exposure (t = -11.94, p < 0.001).
Cohen’s d Effect
Size
cohens_d <- cohen.d(PM25_Estimate ~ Minority_Group, data = nyc_sf_buffer)
print(cohens_d)
##
## Cohen's d
##
## d estimate: -0.7390159 (medium)
## 95 percent confidence interval:
## lower upper
## -0.8646135 -0.6134183
Result: Cohen’s d = -0.74 (medium-large effect),
indicating a substantial difference in PM2.5 exposure between low and
high minority areas.
Spatial Autocorrelation
Analysis
# Create spatial weights matrix
tracts_nb_buffer <- poly2nb(nyc_sf_buffer, queen = TRUE)
# Remove isolated tracts
isolated_tracts <- which(card(tracts_nb_buffer) == 0)
if (length(isolated_tracts) > 0) {
cat("Removing", length(isolated_tracts), "isolated tract(s)\n")
nyc_sf_buffer <- nyc_sf_buffer[-isolated_tracts, ]
tracts_nb_buffer <- poly2nb(nyc_sf_buffer, queen = TRUE)
}
## Removing 1 isolated tract(s)
tracts_listw_buffer <- nb2listw(tracts_nb_buffer, style = "W", zero.policy = TRUE)
# Global Moran's I test
moran_test_buffer <- moran.test(nyc_sf_buffer$PM25_Estimate,
tracts_listw_buffer,
alternative = "two.sided")
print(moran_test_buffer)
##
## Moran I test under randomisation
##
## data: nyc_sf_buffer$PM25_Estimate
## weights: tracts_listw_buffer
##
## Moran I statistic standard deviate = 49.727, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.9259343190 -0.0009606148 0.0003474406
Interpretation: Moran’s I = 0.926 (p < 0.001)
indicates very strong positive spatial autocorrelation. PM2.5 exposure
is highly clustered spatially, justifying the use of spatial regression
models.
Spatial Regression
Modeling
Ordinary Least
Squares (OLS) Baseline
formula_buffer <- PM25_Estimate ~ Minority_Pct + MedIncome
ols_model_buffer <- lm(formula_buffer, data = nyc_sf_buffer)
summary(ols_model_buffer)
##
## Call:
## lm(formula = formula_buffer, data = nyc_sf_buffer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52486 -0.21500 -0.12699 0.08446 1.24401
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.991e+00 5.566e-02 107.63 < 2e-16 ***
## Minority_Pct 5.631e-03 4.536e-04 12.41 < 2e-16 ***
## MedIncome 2.812e-06 3.874e-07 7.26 7.68e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3272 on 1023 degrees of freedom
## (16 observations deleted due to missingness)
## Multiple R-squared: 0.1417, Adjusted R-squared: 0.14
## F-statistic: 84.42 on 2 and 1023 DF, p-value: < 2.2e-16
Lagrange Multiplier
Tests
lm_tests_buffer <- lm.RStests(ols_model_buffer, tracts_listw_buffer,
test = "all", zero.policy = TRUE)
print(lm_tests_buffer)
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = formula_buffer, data = nyc_sf_buffer)
## test weights: tracts_listw_buffer
##
## RSerr = 2097.4, df = 1, p-value < 2.2e-16
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = formula_buffer, data = nyc_sf_buffer)
## test weights: tracts_listw_buffer
##
## RSlag = 1100, df = 1, p-value < 2.2e-16
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = formula_buffer, data = nyc_sf_buffer)
## test weights: tracts_listw_buffer
##
## adjRSerr = 1001.8, df = 1, p-value < 2.2e-16
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = formula_buffer, data = nyc_sf_buffer)
## test weights: tracts_listw_buffer
##
## adjRSlag = 4.3399, df = 1, p-value = 0.03723
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = formula_buffer, data = nyc_sf_buffer)
## test weights: tracts_listw_buffer
##
## SARMA = 2101.7, df = 2, p-value < 2.2e-16
Interpretation: All LM tests are highly significant
(p < 0.001), confirming the need for spatial regression models.
Spatial Error Model
(SEM)
sem_model_buffer <- errorsarlm(formula_buffer,
data = nyc_sf_buffer,
listw = tracts_listw_buffer,
zero.policy = TRUE)
summary(sem_model_buffer)
##
## Call:
## errorsarlm(formula = formula_buffer, data = nyc_sf_buffer, listw = tracts_listw_buffer,
## zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4268248 -0.0093148 -0.0045955 0.0053407 0.4714222
##
## Type: error
## Regions with no neighbours included:
## 266
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.5765e+00 5.7735e-02 113.9078 < 2.2e-16
## Minority_Pct 7.0153e-04 2.6842e-04 2.6136 0.008961
## MedIncome -9.4460e-08 1.2820e-07 -0.7368 0.461232
##
## Lambda: 0.95843, LR test value: 2335.2, p-value: < 2.22e-16
## Asymptotic standard error: 0.0075307
## z-value: 127.27, p-value: < 2.22e-16
## Wald statistic: 16198, p-value: < 2.22e-16
##
## Log likelihood: 859.3885 for error model
## ML residual variance (sigma squared): 0.0080171, (sigma: 0.089538)
## Number of observations: 1026
## Number of parameters estimated: 5
## AIC: -1708.8, (AIC for lm: 624.44)
Key Findings: - Minority_Pct coefficient: 0.000702
(p = 0.009) - statistically significant - Lambda (spatial error
parameter): 0.958 (p < 0.001) - very strong spatial dependence - AIC:
-1708.8 (much better than OLS: 624.4)
Spatial Lag Model
(SAR)
sar_model_buffer <- lagsarlm(formula_buffer,
data = nyc_sf_buffer,
listw = tracts_listw_buffer,
zero.policy = TRUE)
summary(sar_model_buffer)
##
## Call:
## lagsarlm(formula = formula_buffer, data = nyc_sf_buffer, listw = tracts_listw_buffer,
## zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3988860 -0.0740661 -0.0482718 0.0044738 4.2917750
##
## Type: lag
## Regions with no neighbours included:
## 266
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.0610e+00 1.3447e-01 15.3270 < 2.2e-16
## Minority_Pct 1.5072e-03 2.8141e-04 5.3558 8.519e-08
## MedIncome 1.0613e-06 2.3135e-07 4.5877 4.482e-06
##
## Rho: 0.6604, LR test value: 960.81, p-value: < 2.22e-16
## Asymptotic standard error: 0.020761
## z-value: 31.81, p-value: < 2.22e-16
## Wald statistic: 1011.9, p-value: < 2.22e-16
##
## Log likelihood: 172.1858 for lag model
## ML residual variance (sigma squared): 0.037824, (sigma: 0.19448)
## Number of observations: 1026
## Number of parameters estimated: 5
## AIC: -334.37, (AIC for lm: 624.44)
## LM test for residual autocorrelation
## test value: 427.51, p-value: < 2.22e-16
Model
Comparison
model_comparison_buffer <- data.frame(
Model = c("OLS", "SEM", "SAR"),
AIC = c(AIC(ols_model_buffer), AIC(sem_model_buffer), AIC(sar_model_buffer))
)
knitr::kable(model_comparison_buffer,
caption = "Model Comparison by AIC (Lower is Better)",
digits = 2)
Model Comparison by AIC (Lower is Better)
| OLS |
624.44 |
| SEM |
-1708.78 |
| SAR |
-334.37 |
best_model <- model_comparison_buffer$Model[which.min(model_comparison_buffer$AIC)]
Best Model: SEM with AIC = -1708.78
Visualization
Exposure Zone
Map
tmap_mode("plot")
tm_shape(nyc_sf_buffer) +
tm_fill(
col = "PM25_Exposure_Category",
palette = c("High (0-1 mile)" = "#8B0000",
"Moderate (1-3 miles)" = "#FF6347",
"Background (3-5 miles)" = "#FFA500",
"Reference (>5 miles)" = "#FFFFE0"),
title = "JFK Proximity Zone"
) +
tm_borders(col = "gray40", lwd = 0.3) +
tm_layout(
title = "PM2.5 Exposure Zones Around JFK Airport",
legend.position = c("right", "bottom"),
frame = FALSE
) +
tm_compass(type = "8star", position = c("left", "top"), size = 2) +
tm_scale_bar(position = c("left", "bottom"))
PM2.5 Continuous
Map
tm_shape(nyc_sf_buffer) +
tm_fill(
col = "PM25_Estimate",
style = "cont",
palette = "YlOrRd",
title = "PM2.5 (µg/m³)"
) +
tm_borders(col = "gray40", lwd = 0.3) +
tm_layout(
title = "Estimated PM2.5 Exposure by Census Tract",
legend.position = c("right", "bottom"),
frame = FALSE
)
Bivariate Map
bivariate_data_buffer <- bi_class(
nyc_sf_buffer,
x = PM25_Estimate,
y = Minority_Pct,
style = "quantile",
dim = 3
)
bivariate_map_buffer <- ggplot(bivariate_data_buffer) +
geom_sf(aes(fill = bi_class), color = "white", size = 0.1, show.legend = FALSE) +
bi_scale_fill(pal = "DkViolet", dim = 3) +
labs(title = "PM2.5 Exposure and Minority Population Co-Location") +
bi_theme() +
theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"))
legend_buffer <- bi_legend(
pal = "DkViolet",
dim = 3,
xlab = "Higher PM2.5 →",
ylab = "Higher Minority % →",
size = 10
)
print(bivariate_map_buffer)
print(legend_buffer)
Box Plot
Comparison
ggplot(nyc_sf_buffer,
aes(x = Minority_Group, y = PM25_Estimate, fill = Minority_Group)) +
geom_boxplot(alpha = 0.7, outlier.alpha = 0.5) +
geom_jitter(width = 0.2, alpha = 0.15, size = 0.5) +
scale_fill_manual(values = c("Low Minority" = "#4682B4", "High Minority" = "#DC143C")) +
labs(
title = "PM2.5 Exposure by Minority Status",
x = "Demographic Group",
y = "Estimated PM2.5 (µg/m³)",
caption = paste0("Welch's t = ", round(t_test_buffer$statistic, 2),
", p < 0.001\nCohen's d = ", round(cohens_d$estimate, 2),
" (medium-large effect)")
) +
theme_minimal() +
theme(
legend.position = "none",
plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
plot.caption = element_text(hjust = 0.5, size = 10)
)
LISA Cluster
Map
# LISA analysis
local_moran <- localmoran(nyc_sf_buffer$PM25_Estimate, tracts_listw_buffer)
# Classify into categories
nyc_sf_buffer$lisa_category <- case_when(
local_moran[,5] > 0.05 ~ "Not Significant",
local_moran[,1] > 0 & nyc_sf_buffer$PM25_Estimate > mean(nyc_sf_buffer$PM25_Estimate) ~ "High-High",
local_moran[,1] > 0 & nyc_sf_buffer$PM25_Estimate < mean(nyc_sf_buffer$PM25_Estimate) ~ "Low-Low",
local_moran[,1] < 0 & nyc_sf_buffer$PM25_Estimate > mean(nyc_sf_buffer$PM25_Estimate) ~ "High-Low",
TRUE ~ "Low-High"
) %>% factor(levels = c("High-High", "Low-Low", "High-Low", "Low-High", "Not Significant"))
tm_shape(nyc_sf_buffer) +
tm_fill("lisa_category",
palette = c("High-High" = "red", "Low-Low" = "blue",
"High-Low" = "pink", "Low-High" = "lightblue",
"Not Significant" = "grey"),
title = "LISA PM2.5 Clusters") +
tm_borders(col = "gray40", lwd = 0.3) +
tm_layout(title = "Local Spatial Clusters (LISA)", frame = FALSE)
Interpretation: The LISA map reveals distinct
spatial clusters of high PM2.5 exposure (red areas) concentrated near
JFK Airport, confirming local spatial autocorrelation patterns.
MAUP Sensitivity
Analysis
Block Group Level
Analysis
# Get block group demographics
bg_tracts <- tidycensus::get_acs(
geography = "block group",
variables = census_vars,
state = "NY",
county = nyc_counties,
year = 2023,
survey = "acs5",
geometry = TRUE,
output = "wide"
) %>%
dplyr::select(GEOID, NAME, TotalPop = TotalPopE, MedIncome = MedIncomeE) %>%
sf::st_transform(crs = crs_projected_epsg) %>%
st_intersection(study_area)
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 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% | |======================================================================| 99% | |======================================================================| 100%
# Get race data
bg_race_data <- tidycensus::get_acs(
geography = "block group",
variables = minority_vars,
state = "NY",
county = nyc_counties,
year = 2023,
survey = "acs5",
geometry = FALSE,
output = "wide"
) %>%
dplyr::select(GEOID, TotalPop_Race = TotalPop_RaceE, WhiteNH_Pop = WhiteNH_PopE)
# Calculate distances and assign exposure
bg_centroids <- st_centroid(bg_tracts)
bg_tracts$Dist_to_JFK <- st_distance(bg_centroids, jfk_point) %>% as.numeric()
bg_tracts <- bg_tracts %>%
mutate(
PM25_Exposure_Category = case_when(
Dist_to_JFK <= 5280 * 1 ~ "High (0-1 mile)",
Dist_to_JFK <= 5280 * 3 ~ "Moderate (1-3 miles)",
Dist_to_JFK <= 5280 * 5 ~ "Background (3-5 miles)",
TRUE ~ "Reference (>5 miles)"
),
PM25_Estimate = case_when(
Dist_to_JFK <= 5280 * 1 ~ 8.5,
Dist_to_JFK <= 5280 * 3 ~ 7.8,
Dist_to_JFK <= 5280 * 5 ~ 7.2,
TRUE ~ 6.5
)
)
# Integrate demographics
bg_final <- bg_tracts %>%
left_join(bg_race_data, by = "GEOID") %>%
mutate(
Minority_Pop = TotalPop - WhiteNH_Pop,
Minority_Pct = (Minority_Pop / TotalPop) * 100
) %>%
filter(TotalPop > 0, !is.na(Minority_Pct), !is.na(PM25_Estimate)) %>%
select(GEOID, NAME, TotalPop, MedIncome, PM25_Estimate,
Minority_Pct, Dist_to_JFK, geometry)
cat("Block group dataset:", nrow(bg_final), "units\n")
## Block group dataset: 2600 units
ZCTA Level
Analysis
# Get ZCTA demographics
zcta_tracts <- tidycensus::get_acs(
geography = "zcta",
variables = census_vars,
year = 2020,
survey = "acs5",
geometry = TRUE,
output = "wide"
) %>%
dplyr::select(GEOID, NAME, TotalPop = TotalPopE, MedIncome = MedIncomeE) %>%
sf::st_transform(crs = crs_projected_epsg) %>%
st_intersection(study_area)
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 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% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 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%
# Get race data
zcta_race_data <- tidycensus::get_acs(
geography = "zcta",
variables = minority_vars,
year = 2020,
survey = "acs5",
geometry = FALSE,
output = "wide"
) %>%
dplyr::select(GEOID, TotalPop_Race = TotalPop_RaceE, WhiteNH_Pop = WhiteNH_PopE)
# Calculate distances and assign exposure
zcta_centroids <- st_centroid(zcta_tracts)
zcta_tracts$Dist_to_JFK <- st_distance(zcta_centroids, jfk_point) %>% as.numeric()
zcta_tracts <- zcta_tracts %>%
mutate(
PM25_Exposure_Category = case_when(
Dist_to_JFK <= 5280 * 1 ~ "High (0-1 mile)",
Dist_to_JFK <= 5280 * 3 ~ "Moderate (1-3 miles)",
Dist_to_JFK <= 5280 * 5 ~ "Background (3-5 miles)",
TRUE ~ "Reference (>5 miles)"
),
PM25_Estimate = case_when(
Dist_to_JFK <= 5280 * 1 ~ 8.5,
Dist_to_JFK <= 5280 * 3 ~ 7.8,
Dist_to_JFK <= 5280 * 5 ~ 7.2,
TRUE ~ 6.5
)
)
# Integrate demographics
zcta_final <- zcta_tracts %>%
left_join(zcta_race_data, by = "GEOID") %>%
mutate(
Minority_Pop = TotalPop - WhiteNH_Pop,
Minority_Pct = (Minority_Pop / TotalPop) * 100
) %>%
filter(TotalPop > 0, !is.na(Minority_Pct), !is.na(PM25_Estimate)) %>%
select(GEOID, NAME, TotalPop, MedIncome, PM25_Estimate,
Minority_Pct, Dist_to_JFK, geometry)
cat("ZCTA dataset:", nrow(zcta_final), "units\n")
## ZCTA dataset: 114 units
Spatial Regression
at Multiple Scales
# Block Group Analysis
bg_nb <- poly2nb(bg_final, queen = TRUE)
isolated_bg <- which(card(bg_nb) == 0)
if (length(isolated_bg) > 0) {
bg_final <- bg_final[-isolated_bg, ]
bg_nb <- poly2nb(bg_final, queen = TRUE)
}
bg_listw <- nb2listw(bg_nb, style = "W", zero.policy = TRUE)
moran_bg <- moran.test(bg_final$PM25_Estimate, bg_listw, alternative = "two.sided")
ols_bg <- lm(PM25_Estimate ~ Minority_Pct + MedIncome, data = bg_final)
sem_bg <- errorsarlm(PM25_Estimate ~ Minority_Pct + MedIncome,
data = bg_final, listw = bg_listw, zero.policy = TRUE)
# ZCTA Analysis
zcta_nb <- poly2nb(zcta_final, queen = TRUE)
isolated_zcta <- which(card(zcta_nb) == 0)
if (length(isolated_zcta) > 0) {
zcta_final <- zcta_final[-isolated_zcta, ]
zcta_nb <- poly2nb(zcta_final, queen = TRUE)
}
zcta_listw <- nb2listw(zcta_nb, style = "W", zero.policy = TRUE)
moran_zcta <- moran.test(zcta_final$PM25_Estimate, zcta_listw, alternative = "two.sided")
ols_zcta <- lm(PM25_Estimate ~ Minority_Pct + MedIncome, data = zcta_final)
sem_zcta <- errorsarlm(PM25_Estimate ~ Minority_Pct + MedIncome,
data = zcta_final, listw = zcta_listw, zero.policy = TRUE)
Multi-Scale
Comparison Table
maup_comparison <- data.frame(
Scale = c("Census Tract", "Block Group", "ZCTA"),
N_Units = c(nrow(nyc_sf_buffer), nrow(bg_final), nrow(zcta_final)),
Morans_I = c(
round(moran_test_buffer$estimate[1], 3),
round(moran_bg$estimate[1], 3),
round(moran_zcta$estimate[1], 3)
),
OLS_Coef = c(
round(coef(ols_model_buffer)["Minority_Pct"], 6),
round(coef(ols_bg)["Minority_Pct"], 6),
round(coef(ols_zcta)["Minority_Pct"], 6)
),
SEM_Coef = c(
round(coef(sem_model_buffer)["Minority_Pct"], 6),
round(coef(sem_bg)["Minority_Pct"], 6),
round(coef(sem_zcta)["Minority_Pct"], 6)
),
SEM_AIC = c(
round(AIC(sem_model_buffer), 1),
round(AIC(sem_bg), 1),
round(AIC(sem_zcta), 1)
)
)
knitr::kable(maup_comparison,
caption = "Multi-Scale MAUP Sensitivity Analysis",
digits = c(0, 0, 3, 6, 6, 1))
Multi-Scale MAUP Sensitivity Analysis
| Census Tract |
1042 |
0.926 |
0.005631 |
0.000702 |
-1708.8 |
| Block Group |
2599 |
0.957 |
0.004397 |
0.000295 |
-5124.1 |
| ZCTA |
114 |
0.721 |
0.003037 |
-0.000277 |
-8.1 |
# Calculate coefficient of variation
sem_coefs <- maup_comparison$SEM_Coef
coefficient_cv <- (sd(sem_coefs) / mean(sem_coefs)) * 100
Coefficient of Variation: 204.9%
MAUP Conclusion: Findings sensitive to spatial scale
choice (High MAUP sensitivity)
Coefficient
Comparison Plot
coef_plot_data <- data.frame(
Scale = factor(c("Census Tract", "Block Group", "ZCTA"),
levels = c("Block Group", "Census Tract", "ZCTA")),
Coefficient = c(
coef(sem_model_buffer)["Minority_Pct"],
coef(sem_bg)["Minority_Pct"],
coef(sem_zcta)["Minority_Pct"]
),
SE = c(
summary(sem_model_buffer)$Coef["Minority_Pct", "Std. Error"],
summary(sem_bg)$Coef["Minority_Pct", "Std. Error"],
summary(sem_zcta)$Coef["Minority_Pct", "Std. Error"]
)
) %>%
mutate(
Lower_CI = Coefficient - 1.96 * SE,
Upper_CI = Coefficient + 1.96 * SE
)
ggplot(coef_plot_data, aes(x = Scale, y = Coefficient, color = Scale)) +
geom_point(size = 4) +
geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI), width = 0.2) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black", alpha = 0.6) +
scale_color_manual(values = c("Block Group" = "#DC143C",
"Census Tract" = "#FF6347",
"ZCTA" = "#4682B4")) +
labs(
title = "Impact of Spatial Scale on Environmental Justice Finding (SEM)",
subtitle = "Coefficient for Minority Percentage vs. PM2.5 Exposure (95% CI)",
x = "Geographic Scale (MAUP Analysis)",
y = expression(paste(beta, " Coefficient (PM2.5 per % Minority)"))
) +
coord_flip() +
theme_minimal() +
theme(
legend.position = "none",
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5)
)
Interpretation: The plot shows that the minority
percentage coefficient varies across scales, with census tracts and
block groups showing significant positive relationships, while the ZCTA
level shows no significant effect (confidence interval crosses
zero).
Model Diagnostics
Exploratory
Plots
# Histogram by group
p1 <- ggplot(nyc_sf_buffer, aes(x = PM25_Estimate, fill = Minority_Group)) +
geom_histogram(alpha = 0.6, position = "identity", bins = 30) +
labs(title = "PM2.5 Distribution by Minority Status", fill = "Group") +
theme_minimal()
# Density plot
p2 <- ggplot(nyc_sf_buffer, aes(x = PM25_Estimate, color = Minority_Group)) +
geom_density(linewidth = 1.2) +
labs(title = "PM2.5 Density by Minority Status", color = "Group") +
theme_minimal()
# Scatterplot
p3 <- ggplot(nyc_sf_buffer, aes(x = Minority_Pct, y = PM25_Estimate)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm") +
labs(title = "Bivariate Relationship",
x = "Minority Percentage",
y = "PM2.5 Estimate") +
theme_minimal()
print(p1)
print(p2)
print(p3)
Permutation
Test
perm_test <- coin::independence_test(PM25_Estimate ~ Minority_Group,
data = nyc_sf_buffer,
distribution = coin::approximate(nresample = 10000))
print(perm_test)
##
## Approximative General Independence Test
##
## data: PM25_Estimate by
## Minority_Group (Low Minority, High Minority)
## Z = -11.257, p-value < 1e-04
## alternative hypothesis: two.sided
Result: The permutation test confirms the
significant difference between groups (Z = -11.26, p < 0.001),
validating our parametric t-test results.
SEM Residual
Diagnostics
# Extract SEM residuals
sem_residuals <- residuals(sem_model_buffer)
# Align spatial data
nyc_sf_model <- nyc_sf_buffer[names(sem_residuals), ]
# Rebuild neighbors and weights
nb_model <- spdep::poly2nb(nyc_sf_model, queen = TRUE)
listw_model <- spdep::nb2listw(nb_model, style = "W", zero.policy = TRUE)
# Moran's I on residuals
moran_residuals <- spdep::moran.test(sem_residuals, listw_model, zero.policy = TRUE)
print(moran_residuals)
##
## Moran I test under randomisation
##
## data: sem_residuals
## weights: listw_model
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = -2.8303, p-value = 0.9977
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.0543201735 -0.0009765625 0.0003552321
Interpretation: Moran’s I on residuals = -0.0543 (p
= 0.9977). The non-significant result and negative Moran’s I indicates
that the SEM successfully removed spatial autocorrelation from the
residuals.
Q-Q Plot
qqnorm(sem_residuals, main = "Q-Q Plot of SEM Residuals")
qqline(sem_residuals, col = "red", lwd = 2)
Interpretation: The Q-Q plot shows some deviation
from normality in the tails, but the overall pattern is reasonable for
spatial data. The slight deviations suggest some heavy-tailed
distribution, which is common in spatial environmental data.
Summary Table
summary_table <- data.frame(
Metric = c(
"Study Design",
"Total Census Tracts",
"Total Population",
"Median Minority %",
"",
"Mean PM2.5 (Low Minority)",
"Mean PM2.5 (High Minority)",
"Absolute Difference",
"Relative Difference",
"",
"T-statistic",
"P-value",
"Cohen's d Effect Size",
"",
"Moran's I",
"Spatial Clustering",
"Best Model (AIC)",
"",
"Minority % in High Zone (0-3 mi)",
"Minority % in Reference (>5 mi)"
),
Value = c(
"Buffer-based proximity",
as.character(nrow(nyc_sf_buffer)),
format(sum(nyc_sf_buffer$TotalPop), big.mark = ","),
paste0(round(minority_threshold, 1), "%"),
"",
paste0(round(descriptive_stats_buffer$Mean_PM25[1], 2), " µg/m³"),
paste0(round(descriptive_stats_buffer$Mean_PM25[2], 2), " µg/m³"),
paste0("+", round(descriptive_stats_buffer$Mean_PM25[2] - descriptive_stats_buffer$Mean_PM25[1], 2), " µg/m³"),
paste0("+", round((descriptive_stats_buffer$Mean_PM25[2] / descriptive_stats_buffer$Mean_PM25[1] - 1) * 100, 1), "%"),
"",
as.character(round(t_test_buffer$statistic, 2)),
"< 0.001",
paste0(round(abs(cohens_d$estimate), 2), " (medium-large)"),
"",
as.character(round(moran_test_buffer$estimate[1], 3)),
"Yes (p < 0.001)",
best_model,
"",
paste0(round(mean(exposure_zone_summary$Mean_Minority_Pct[1:2]), 1), "%"),
paste0(round(exposure_zone_summary$Mean_Minority_Pct[3], 1), "%")
)
)
knitr::kable(summary_table,
caption = "Comprehensive Analysis Summary",
col.names = c("Metric", "Value"))
Comprehensive Analysis Summary
| Study Design |
Buffer-based proximity |
| Total Census Tracts |
1042 |
| Total Population |
3,543,929 |
| Median Minority % |
86.9% |
|
|
| Mean PM2.5 (Low Minority) |
6.55 µg/m³ |
| Mean PM2.5 (High Minority) |
6.79 µg/m³ |
| Absolute Difference |
+0.24 µg/m³ |
| Relative Difference |
+3.7% |
|
|
| T-statistic |
-11.94 |
| P-value |
< 0.001 |
| Cohen’s d Effect Size |
0.74 (medium-large) |
|
|
| Moran’s I |
0.926 |
| Spatial Clustering |
Yes (p < 0.001) |
| Best Model (AIC) |
SEM |
|
|
| Minority % in High Zone (0-3 mi) |
93.5% |
| Minority % in Reference (>5 mi) |
73.9% |
Conclusions
Main Findings
Environmental Justice Disparity Confirmed:
High-minority census tracts experience significantly higher PM2.5
exposure (0.24 µg/m³ difference, p < 0.001) with a medium-large
effect size (Cohen’s d = 0.74).
Strong Spatial Autocorrelation: Moran’s I =
0.926 indicates highly clustered PM2.5 exposure patterns, necessitating
spatial regression approaches.
Model Selection: Spatial Error Model (SEM)
provides the best fit (AIC = -1708.8), substantially outperforming OLS
(AIC = 624.4) and SAR (AIC = -334.4).
MAUP Sensitivity: High coefficient of variation
(204.9%) across spatial scales indicates that findings are sensitive to
the choice of geographic unit, with the relationship weakening at
coarser (ZCTA) scales.
Limitations
- PM2.5 estimates based on proximity buffers rather than actual
monitoring data
- Cross-sectional design limits causal inference
- MAUP sensitivity suggests results should be interpreted
cautiously
- Residual normality assumptions partially violated
Implications
The analysis provides evidence of environmental justice concerns near
JFK Airport, with minority communities experiencing elevated PM2.5
exposure. However, the high MAUP sensitivity indicates that policy
interventions should consider multiple spatial scales and verify
findings with actual air quality monitoring data.