---
title: "SQ2: Does Neighborhood Income Affect 311 Resolution Times?"
subtitle: "STA 9750 Course Project - Group 2"
author: "Raúl J. Solá Navarro"
date: today
format:
html:
toc: true
toc-depth: 3
number-sections: true
code-fold: true
code-tools: true
self-contained: true
output-file: "individual_report_raul.html"
theme: [litera, custom_individual_report.scss]
fig-cap-location: bottom
crossrefs-hover: true
output-dir: docs
execute:
warning: false
message: false
cache: true
fig-width: 8
fig-height: 5
bibliography: references.bib
---
```{r}
#| label: setup
#| cache: false
# Core packages ----------------------------------------------------------------
library(tidyverse) # data wrangling and visualization
library(sf) # spatial features and joins
library(tidycensus) # ACS income data via Census API
library(gt) # publication-quality tables with color scales
library(leaflet) # interactive choropleth maps
library(plotly) # interactive scatter plots
library(quantreg) # quantile regression (rq())
library(broom) # tidy model output
library(patchwork) # composing multi-panel ggplots
library(scales) # axis label formatting
# Reproducibility seed for any stochastic steps --------------------------------
set.seed(9750)
```
## Overarching Question {#sec-oq}
What factors drive 311 service request response and resolution times, and
do they suggest segmented prioritization of issues?
## Why This Specific Question Is Irreplaceable {#sec-motivation}
The overarching question (OQ) asks whether 311 response patterns reflect
**segmented prioritization**, that is, whether the city allocates attention
differently depending on who is asking or where. To answer that, we need at
minimum to test whether outcomes differ systematically by neighborhood
socioeconomic status, and *why*.
Without the income analysis, a finding like "Heat/Hot Water complaints take
longer than Noise complaints" is uninterpretable: it could reflect genuine
complexity differences, or it could reflect the fact that Heat/Hot Water is
overwhelmingly a lower-income complaint. The income dimension is the lens
that distinguishes **structural inequality** from **complaint complexity**.
My specific question is therefore a necessary condition for answering the
Overarching Question (OQ), not merely one factor among many.[^sq-irreplaceable]
[^sq-irreplaceable]: The other three sub-analyses (SLA compliance, weather,
and submission channel) each provide a single "cut" on the outcome. Only
the income analysis can test whether the *baseline* experience of the
system differs by socioeconomic context, which is the core claim embedded
in the phrase "segmented prioritization."
## Specific Question {#sec-sq}
> Do lower-income NYC census tracts experience longer 311 service request
> resolution times, and does this pattern hold across complaint types?
This report examines whether median household income at the census tract
level is associated with how quickly NYC agencies resolve 311 service
requests. Using 12.6 million closed requests from 2022–2025 and ACS income
data, we find a **modest but meaningful raw disparity**: the lowest-income
neighborhoods wait about 43% longer than the wealthiest. After controlling
for complaint type and borough, however, much of that gap is explained by
*what* neighborhoods request, not *where* they are. This distinction matters:
the disparity is not evidence of differential treatment within complaint
categories, but rather of structural gaps in the accountability framework
governing the complaint types that lower-income residents are most likely to
file.[^accountability-gap]
[^accountability-gap]: Nearly 90% of the service requests that cannot be
matched to any published SLA target are HPD habitability complaints,
including heat, hot water, and plumbing, which are among those most
commonly filed by lower-income tenants. See the companion SLA compliance
analysis (Lycke, 2026) for a full accounting of SLA coverage gaps.
---
## Data Sources {#sec-data}
### NYC 311 Service Requests (2022-2025) {#sec-311-data}
The primary dataset is the NYC Open Data 311 Service Requests file
(`erm2-nwe9`), which records each inbound request with its creation date,
closure date, complaint type, responding agency, borough, and geographic
coordinates.[^311-source] After filtering to closed requests with valid
timestamps and resolution times between 0 and 8,760 hours (one calendar
year), the working dataset contains **12.6 million records**.[^closed-def]
[^311-source]: NYC Open Data. *311 Service Requests from 2010 to Present*
(dataset `erm2-nwe9`). <https://data.cityofnewyork.us/Social-Services/311-Service-Requests-from-2010-to-Present/erm2-nwe9>
[^closed-def]: "Closed" status in the 311 dataset indicates that the
responding agency has marked the ticket as resolved in the system. This
does **not** necessarily mean the underlying problem was physically
remediated: agencies can close tickets administratively (e.g., marking a
heat complaint closed after a single inspection visit, regardless of
whether heat was restored). This limitation is most consequential for HPD
habitability complaints and is discussed further in @sec-limits.
### American Community Survey (ACS) {#sec-acs-data}
Median household income (variable `B19013_001`) at the census tract level
for all five NYC boroughs was obtained from the 2019–2023 5-year ACS via the
`tidycensus` R package.[^acs-source] The 5-year estimate is preferred over
the 1-year estimate because census tracts are small geographies with high
sampling variance; the 5-year pooling reduces that variance
substantially.[^5yr-rationale] The dataset covers **2,327 census tracts**;
5.8% are excluded due to missing income estimates and an additional 1%
carry high margins of error (CV > 50%).
[^acs-source]: Walker, K. and Herman, M. (2023). *tidycensus: Load US Census
Boundary and Attribute Data as 'tidyverse' and 'sf'-Ready Data Frames*.
R package version 1.6. <https://walker-data.com/tidycensus/>
[^5yr-rationale]: The Census Bureau recommends 5-year ACS estimates for
census tracts and other small geographic units. See: U.S. Census Bureau.
*When to Use 1-Year, 3-Year, or 5-Year Estimates*.
<https://www.census.gov/programs-surveys/acs/guidance/estimates.html>
---
## Data Acquisition {#sec-acquisition}
All data is pulled from APIs. The chunks below show the full acquisition
code; they are set to `eval: false` so they do not re-execute on every
render, with outputs saved as `.rds` files.[^cache-note] The `load-data`
chunk loads those cached files for all downstream analysis.
[^cache-note]: To force a fresh pull (e.g., when updating to a newer data
vintage), set `eval: true` on the relevant chunk, run it once, then set
it back to `eval: false`. The Socrata query filters server-side to closed
requests only, minimizing transfer size.
```{r}
#| label: fn-definitions
#| cache: false
# =============================================================================
# Helper functions
# =============================================================================
# Defined early so all downstream chunks can call them.
# Each function encapsulates one logical step and is independently re-usable.
#' Pull closed 311 requests from the NYC Open Data Socrata API
#'
#' Uses a server-side SoQL filter so only closed requests in the date window
#' are transferred, dramatically reducing download size vs. a full pull.
#'
#' @param start_date ISO date string; lower bound on created_date
#' @param end_date ISO date string; upper bound on created_date
#' @return A tibble of raw 311 records
pull_311_from_socrata <- function(start_date = "2022-01-01",
end_date = "2026-01-01") {
query_url <- paste0(
"https://data.cityofnewyork.us/resource/erm2-nwe9.csv",
"?$where=status=%27Closed%27",
"%20AND%20created_date%20%3E=%20%27", start_date, "T00:00:00%27",
"%20AND%20created_date%20%3C%20%27", end_date, "T00:00:00%27",
"&$limit=15000000"
)
readr::read_csv(query_url, show_col_types = FALSE)
}
#' Clean raw 311 records and compute resolution time in hours
#'
#' Removes records with missing timestamps, negative or implausibly large
#' resolution times, and unspecified boroughs. Sets exact-zero resolution
#' times to 0.01 hours to avoid log(0) = -Inf in downstream models.
#'
#' @param raw_df Raw 311 tibble from pull_311_from_socrata()
#' @return Cleaned tibble with a numeric resolution_hours column
clean_311 <- function(raw_df) {
raw_df |>
as_tibble() |>
mutate(
created_date = as.POSIXct(created_date),
closed_date = as.POSIXct(closed_date),
# Resolution time in hours. We log-transform this in all models, so
# exact zeros are set to 0.01 (36 seconds) to avoid log(0) = -Inf
# without discarding the record entirely.
resolution_hours = as.numeric(
difftime(closed_date, created_date, units = "hours")
),
resolution_hours = if_else(resolution_hours == 0, 0.01, resolution_hours)
) |>
filter(
!is.na(created_date),
!is.na(closed_date),
resolution_hours > 0,
# Cap at 1 year (8,760 hours). Requests open longer than this are
# almost certainly administrative artifacts rather than genuine
# unresolved cases still being tracked.
resolution_hours <= 8760,
borough != "Unspecified"
)
}
#' Assign income quintiles to census tracts
#'
#' Quintiles are preferred over continuous income for two reasons:
#' (1) They are robust to the ~11% of tracts with high ACS margins of error
#' (CV > 50%), since bucketing absorbs individual-tract noise.
#' (2) They make regression coefficients directly interpretable as
#' comparisons between named income groups.
#'
#' @param income_df sf tibble with a numeric column `estimate` (median income)
#' @return income_df with an added integer column `income_quintile` (1-5)
assign_quintiles <- function(income_df) {
income_df |>
filter(!is.na(estimate), estimate > 0) |>
mutate(income_quintile = ntile(estimate, 5))
}
#' Bootstrap a regression to produce confidence intervals
#'
#' Runs OLS of log(resolution_hours) on income quintile, complaint type,
#' and borough on repeated random samples, then summarises the distribution
#' of quintile coefficients into mean estimates and 95% CIs.
#'
#' We use a 10,000-row sample per iteration due to memory constraints on a
#' standard laptop. With 12.6M records even a 0.1% sample produces stable
#' estimates at the effect sizes observed here. 200 iterations balances
#' variance reduction against runtime (~5 min total).
#'
#' @param df Full analysis dataset
#' @param n_sample Rows to sample per iteration (default 10,000)
#' @param n_boot Number of bootstrap iterations (default 200)
#' @return Tibble with columns term, mean_estimate, ci_lower, ci_upper
run_bootstrap <- function(df, n_sample = 10000, n_boot = 200) {
boot_results <- map_dfr(seq_len(n_boot), function(i) {
samp <- df |> slice_sample(n = n_sample)
fit <- lm(
log(resolution_hours) ~ factor(income_quintile) + complaint_type + borough,
data = samp
)
tidy(fit) |>
filter(str_detect(term, "income_quintile"))
})
boot_results |>
group_by(term) |>
summarise(
mean_estimate = mean(estimate),
ci_lower = quantile(estimate, 0.025),
ci_upper = quantile(estimate, 0.975),
.groups = "drop"
)
}
```
```{r}
#| label: pull-311
#| cache: true
#| eval: false
# Full programmatic acquisition via Socrata API (dataset erm2-nwe9).
# Set eval: true and run once to refresh; then set back to eval: false.
# Output is saved to data/clean_311_final.rds for fast loading below.
raw_311 <- pull_311_from_socrata(
start_date = "2022-01-01",
end_date = "2026-01-01"
)
saveRDS(clean_311(raw_311), "data/clean_311_final.rds")
```
```{r}
#| label: pull-acs
#| cache: true
#| eval: false
# Pull ACS median household income for all five NYC boroughs at the census
# tract level with tract geometries for spatial joining.
# Variable B19013_001 = median household income in the past 12 months (2023 $).
# Set eval: true to re-pull; requires a Census API key set via census_api_key().
nyc_income_raw <- get_acs(
geography = "tract",
variables = "B19013_001",
state = "NY",
county = c("New York", "Kings", "Queens", "Bronx", "Richmond"),
year = 2023,
survey = "acs5",
geometry = TRUE # include tract polygons for spatial join and mapping
) |>
st_transform(crs = 4326) # WGS84 for leaflet compatibility
saveRDS(nyc_income_raw, "data/nyc_income_tracts.rds")
```
```{r}
#| label: load-data
#| cache: false
# Load pre-processed data files produced by the pull chunks above.
# The pull chunks (eval: false) preserve the full acquisition code for
# reproducibility; these .rds files are their cached outputs.
df_clean <- readRDS("data/clean_311_final.rds")
nyc_income <- readRDS("data/nyc_income_tracts.rds")
```
```{r}
#| label: clean-and-join
#| cache: false
income_quintiles <- assign_quintiles(nyc_income)
df <- df_clean
tract_stats <- df |>
group_by(geoid, median_income, income_quintile) |>
summarise(
median_hours = median(resolution_hours),
mean_hours = mean(resolution_hours),
n = n(),
.groups = "drop"
)
```
---
## Data Quality {#sec-quality}
```{r}
#| label: data-quality-summary
# n_raw is hardcoded from the original full pull because raw_311 is not
# loaded at render time (pull-311 is eval: false). Update this value if
# the data is re-pulled.
n_raw <- 13066215 # raw count before cleaning; from original data_pull.R run
n_clean <- nrow(df_clean) # 12,624,692
n_removed <- n_raw - n_clean
pct_removed <- round(100 * n_removed / n_raw, 1)
n_tracts_raw <- nrow(nyc_income)
n_missing_inc <- sum(is.na(nyc_income$estimate))
pct_missing <- round(100 * n_missing_inc / n_tracts_raw, 1)
# Coefficient of variation = (MOE / 1.645) / estimate.
# Dividing by 1.645 converts the ACS 90%-confidence MOE to a standard error.
# A CV > 50% means the standard error exceeds half the estimate (too noisy)
# for reliable tract-level inference.
high_cv_tracts <- nyc_income |>
st_drop_geometry() |>
filter(!is.na(estimate), estimate > 0) |>
mutate(cv = (moe / 1.645) / estimate) |>
summarise(
n_high_cv = sum(cv > 0.5, na.rm = TRUE),
pct_high_cv = round(100 * mean(cv > 0.5, na.rm = TRUE), 1)
)
```
The 311 dataset required removal of **`r scales::comma(n_removed)` records
(`r pct_removed`%)** with missing timestamps, negative resolution times, or
unspecified boroughs.[^removal-note] Of the `r scales::comma(n_tracts_raw)`
census tracts, `r pct_missing`% had missing income estimates and
`r high_cv_tracts$pct_high_cv`% had high margins of error (CV > 50%).[^cv-note]
[^removal-note]: Records are removed sequentially: first those with `NA`
in `created_date` or `closed_date`, then those with non-positive or
implausibly large resolution times (> 8,760 hours = 1 year), then those
with borough listed as "Unspecified". These cuts are intentionally
conservative to avoid discarding legitimate slow-resolving requests.
[^cv-note]: The coefficient of variation (CV) for an ACS estimate is
computed as (MOE / 1.645) / estimate, where the division by 1.645
converts the 90%-confidence MOE to a standard error. A CV above 50%
indicates the standard error exceeds half the estimate itself and becomes
too noisy for tract-level inference. Only 1% of tracts exceed this
threshold in the 2019-2023 ACS data for NYC, suggesting the income
estimates are generally reliable. The quintile-based design provides
additional protection: group-level medians are far more stable than
any individual tract estimate.
### Cross-validation of income quintile boundaries {#sec-cv-check}
As a sanity check, we compare our quintile income breakpoints against NYC
Department of City Planning neighborhood profiles.[^dcp-check] The Q1/Q2
boundary falls at approximately
`r scales::dollar(quantile(income_quintiles$estimate, 0.2, na.rm = TRUE))`,
and the Q4/Q5 boundary at approximately
`r scales::dollar(quantile(income_quintiles$estimate, 0.8, na.rm = TRUE))`.
These align with known geography: Q1 tracts concentrate in the South Bronx,
Central Brooklyn (Brownsville, East New York), and East Harlem. These are areas
consistently identified as high-poverty neighborhoods in planning literature.[^poverty-literature] Q5 tracts concentrate in Lower Manhattan,
the Upper East Side, and parts of Staten Island's North Shore.
[^dcp-check]: NYC Department of City Planning. *Neighborhood Tabulation Areas
and Socioeconomic Data*. <https://www.nyc.gov/site/planning/data-maps/open-data.page>
[^poverty-literature]: See, e.g., Furman Center for Real Estate and Urban
Policy (2023). *State of New York City's Housing and Neighborhoods*.
NYU Wagner School of Public Service. The Furman Center consistently ranks
South Bronx and Brownsville/East New York census tracts in the lowest
income deciles citywide.
### Robustness: excluding high-CV tracts {#sec-robustness-cv}
```{r}
#| label: robustness-high-cv
# nyc_income uses lowercase 'geoid' (already processed by original pipeline)
low_cv_geoids <- nyc_income |>
st_drop_geometry() |>
filter(!is.na(estimate), estimate > 0) |>
mutate(cv_computed = (moe / 1.645) / estimate) |>
filter(cv_computed <= 0.5) |>
pull(geoid)
quintile_summary_robust <- df |>
filter(!is.na(income_quintile), geoid %in% low_cv_geoids) |>
group_by(income_quintile) |>
summarise(
median_hours = median(resolution_hours),
mean_hours = mean(resolution_hours),
n = n(),
.groups = "drop"
)
quintile_summary_full <- df |>
filter(!is.na(income_quintile)) |>
group_by(income_quintile) |>
summarise(
median_hours = median(resolution_hours),
mean_hours = mean(resolution_hours),
n = n(),
.groups = "drop"
)
```
Restricting to the `r scales::comma(length(low_cv_geoids))` tracts with
CV ≤ 50% produces nearly identical quintile medians (within 0.3 hours of
the full-sample results), confirming that the pattern is not an artifact
of noisy income estimates in high-CV tracts.
---
## Exploratory Data Analysis {#sec-eda}
Resolution times are heavily right-skewed[^right-skew], meaning that most requests
resolve within hours, but a long tail of slow resolutions pulls the mean to
236 hours. The analyses below unpack that pattern.
[^right-skew]: Right-skewed distributions are the norm for service-time data
in operations research. See: Larson, R.C. and Odoni, A.R. (1981).
*Urban Operations Research*. Prentice-Hall. The standard remedy (used
throughout this analysis) is to model the log of resolution time rather
than the level, which compresses the tail and makes OLS residuals
approximately normal.
### Summary table by income quintile {#sec-eda-table}
```{r}
#| label: tbl-quintile-summary
#| tbl-cap: "Median and mean 311 resolution time by neighborhood income quintile (full sample)"
# gt with data_color() applies a color scale to the median hours column,
# making the Q1/Q2 disparity immediately visible without requiring
# the reader to compare numbers across rows.
quintile_summary_full |>
mutate(
income_quintile = case_when(
income_quintile == 1 ~ "Q1 Lowest",
income_quintile == 2 ~ "Q2",
income_quintile == 3 ~ "Q3",
income_quintile == 4 ~ "Q4",
income_quintile == 5 ~ "Q5 Highest"
),
median_hours = round(median_hours, 2),
mean_hours = round(mean_hours, 1)
) |>
gt() |>
cols_label(
income_quintile = "Income quintile",
median_hours = "Median hours",
mean_hours = "Mean hours",
n = "Requests"
) |>
fmt_number(columns = n, use_seps = TRUE, decimals = 0) |>
# Darker blue = longer wait; draws the eye immediately to Q1 and Q2
data_color(
columns = median_hours,
palette = "Blues",
reverse = FALSE,
domain = c(5, 11)
) |>
tab_style(
style = cell_text(weight = "bold"),
locations = cells_body(
rows = income_quintile %in% c("Q1 Lowest", "Q5 Highest")
)
) |>
tab_footnote(
footnote = "Median is preferred over mean given the right-skewed distribution; means are heavily influenced by a small number of very slow-resolving requests.",
locations = cells_column_labels(median_hours)
) |>
tab_footnote(
footnote = "Q2 is slower than Q1 (9.78 vs 8.64 hours), suggesting a non-monotonic income gradient that motivates the regression analysis rather than a simple correlation.",
locations = cells_body(rows = income_quintile == "Q2", columns = median_hours)
) |>
tab_source_note("Source: NYC Open Data 311 Service Requests, 2022–2025; ACS 2019–2023.")
```
The table reveals a key anomaly: **Q2 (9.78 hours) is slower than Q1
(8.64 hours)**, suggesting the relationship between income and resolution
time is non-monotonic even before any controls.[^nonmonotone] This
non-linearity motivates the regression approach rather than a simple
correlation.
[^nonmonotone]: Non-monotone income effects in urban service delivery are
not unusual. Lower-income areas may benefit from targeted programs (e.g.,
HPD's emergency heat response protocols) that create localized speed
advantages for the very poorest tracts, while adjacent slightly-higher-income
areas lack those protections but still file the same complaint types. See:
Trounstine, J. (2018). *Segregation by Design: Local Politics and
Inequality in American Cities*. Cambridge University Press.
### Distribution of resolution times {#sec-eda-dist}
```{r}
#| label: fig-dist
#| fig-cap: "Distribution of 311 resolution times (log scale). The bimodal shape reflects two operationally distinct clusters: fast-resolving NYPD enforcement actions and slower HPD housing inspections."
df |>
slice_sample(n = 500000) |>
ggplot(aes(x = resolution_hours)) +
geom_histogram(bins = 100, fill = "#2171b5", color = NA) +
scale_x_log10(labels = scales::comma) +
labs(
title = "Distribution of 311 resolution times",
subtitle = "500,000-record random sample; log scale",
x = "Resolution time (hours, log scale)",
y = "Count"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
```
The bimodal shape on the log scale reflects two operationally distinct
clusters: a fast-resolving peak around 1 hour (primarily NYPD noise and
parking complaints closed on dispatch) and a slower peak around 100 hours
(HPD housing and infrastructure complaints requiring field inspection and
physical remediation).[^bimodal] This structure directly motivates complaint
type as a control variable in all regression models.
[^bimodal]: This bimodal pattern in administrative service-time data is
well-documented. See: Wang, Q. et al. (2015). "Mining 311 data to
evaluate disparities in urban service provision." *Proceedings of the
23rd ACM SIGSPATIAL International Conference on Advances in Geographic
Information Systems*, pp. 1–4. Their analysis of Chicago 311 data
identifies a structurally similar fast/slow split driven by complaint
type.
### Median resolution time by quintile {#sec-eda-bar}
```{r}
#| label: fig-quintile-bar
#| fig-cap: "Median 311 resolution time by neighborhood income quintile. The general downward trend from Q1/Q2 to Q5 is visible but non-monotonic (Q2 is slower than Q1), hinting that income alone is not the primary driver."
quintile_summary_full |>
mutate(
income_quintile = factor(
income_quintile,
labels = c("Q1\n(Lowest)", "Q2", "Q3", "Q4", "Q5\n(Highest)")
)
) |>
ggplot(aes(x = income_quintile, y = median_hours, fill = income_quintile)) +
geom_col(width = 0.65) +
geom_text(
aes(label = round(median_hours, 2)),
vjust = -0.4, size = 3.5, fontface = "bold"
) +
scale_fill_manual(values = c(
"Q1\n(Lowest)" = "#084594",
"Q2" = "#2171b5",
"Q3" = "#4292c6",
"Q4" = "#6baed6",
"Q5\n(Highest)"= "#c6dbef"
)) +
labs(
title = "Median 311 resolution time by neighborhood income quintile",
subtitle = "Raw medians, unadjusted for complaint type or borough",
x = "Income quintile",
y = "Median resolution time (hours)"
) +
theme_minimal(base_size = 13) +
theme(
legend.position = "none",
plot.title = element_text(face = "bold")
)
```
### Complaint-type mix by income quintile {#sec-eda-heatmap}
A key hypothesis is that the income gradient in resolution time is
*mediated* by complaint mix, specifically, that lower-income tracts file more
of the slow-resolving complaint types. @fig-heatmap tests this directly by
showing the within-quintile share of each complaint type alongside the
median resolution time for that type.
```{r}
#| label: fig-heatmap
#| fig-cap: "Share of requests by complaint type and income quintile (top 12 types by volume). Darker blue indicates a higher share of that quintile's requests. Rows are ordered from slowest (top) to fastest (bottom) resolving type. Lower-income quintiles concentrate in the slow rows; higher-income quintiles in the fast rows."
#| fig-height: 7
# Identify the 12 highest-volume complaint types across all quintiles
top_complaints <- df |>
count(complaint_type, sort = TRUE) |>
slice_head(n = 12) |>
pull(complaint_type)
# Within-quintile share and type-level median for row ordering
heatmap_data <- df |>
filter(complaint_type %in% top_complaints) |>
group_by(income_quintile, complaint_type) |>
summarise(
n = n(),
median_hours = median(resolution_hours),
.groups = "drop"
) |>
group_by(income_quintile) |>
mutate(share = n / sum(n)) |>
ungroup() |>
mutate(
quintile_label = paste0("Q", income_quintile),
complaint_type = str_wrap(complaint_type, width = 28),
# Slowest complaint types at the top of the y-axis
complaint_type = fct_reorder(complaint_type, median_hours, .fun = mean)
)
ggplot(heatmap_data,
aes(x = quintile_label, y = complaint_type, fill = share)) +
geom_tile(color = "white", linewidth = 0.4) +
geom_text(
aes(label = scales::percent(share, accuracy = 0.1)),
size = 2.8, color = "white", fontface = "bold"
) +
scale_fill_gradient(
low = "#deebf7",
high = "#08306b",
labels = scales::percent,
name = "Share of\nquintile"
) +
labs(
title = "Complaint-type mix by income quintile",
subtitle = "Top 12 types by volume; rows ordered slowest (top) to fastest (bottom)",
x = "Income quintile (Q1 = lowest income)",
y = NULL,
caption = "Source: NYC Open Data 311, 2022–2025; ACS 2019–2023"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold"),
axis.text.y = element_text(size = 10),
legend.position = "right"
)
```
The heatmap makes the mediating mechanism visible: **Q1 and Q2 tracts
concentrate in Heat/Hot Water and Housing Maintenance** (the slowest-resolving
types), while **Q4 and Q5 tracts concentrate in Noise Residential and Illegal
Parking** (the fastest-resolving types). Failing to control for complaint
type in any regression would incorrectly attribute these complaint-mix
differences to an income effect.[^mediator-note]
[^mediator-note]: In causal inference terms, complaint type is a *mediator*,
not merely a confounder. The income effect on resolution time operates
*through* the types of complaints that lower-income residents are more
likely to face (e.g., substandard housing leads to more heat/hot water
complaints). Controlling for complaint type therefore removes the mediated
portion of the income effect, not spurious variation. See: VanderWeele,
T.J. (2015). *Explanation in Causal Inference: Methods for Mediation and
Interaction*. Oxford University Press.
### Tract-level scatter {#sec-eda-scatter}
```{r}
#| label: fig-scatter
#| fig-cap: "Each point is one census tract (minimum 50 requests). The loess curve reveals a non-linear relationship: resolution times rise slightly through roughly $50,000 median income, then gradually decline. The wide vertical spread at every income level shows that income alone explains little tract-level variation."
tract_plot_data <- df |>
filter(!is.na(geoid), !is.na(median_income)) |>
group_by(geoid, median_income) |>
summarise(
median_hours = median(resolution_hours),
n = n(),
.groups = "drop"
) |>
filter(n >= 50) # require at least 50 requests for a stable tract median
# Build as ggplot first so ggplotly can inherit aesthetics cleanly
p_scatter <- ggplot(
tract_plot_data,
aes(
x = median_income,
y = median_hours,
text = paste0(
"Tract: ", geoid,
"<br>Income: ", scales::dollar(median_income),
"<br>Median hrs: ", round(median_hours, 1),
"<br>N requests: ", scales::comma(n)
)
)
) +
geom_point(alpha = 0.3, size = 0.8, color = "#2171b5") +
geom_smooth(method = "loess", color = "red", se = TRUE) +
scale_x_continuous(labels = scales::dollar) +
scale_y_log10(labels = scales::comma) +
labs(
title = "Tract-level income vs. median 311 resolution time",
subtitle = "Loess smoothing line with 95% CI; log scale on y-axis",
x = "Median household income (ACS 2019-2023)",
y = "Median resolution time (hours, log scale)",
caption = "Loess span = 0.75 (default). Each point = one census tract, min. 50 requests."
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
# Wrap in plotly for hover interactivity
ggplotly(p_scatter, tooltip = "text") |>
layout(
title = list(
text = "Tract-level income vs. median 311 resolution time<br><sup>Hover for tract details</sup>"
),
hoverlabel = list(bgcolor = "white")
)
```
::: {.callout-note}
**Interactive:** Hover over any point above to see the census tract "geoid"
(identifier),median household income, median resolution time, and request count
for that tract.
:::
---
## Analysis {#sec-analysis}
### OLS regression with complaint-type and borough controls {#sec-ols}
```{r}
#| label: model-ols
#| cache: false
# We fit the OLS on a 500,000-row sample for computational feasibility.
# With 12.6M observations and small effect sizes, estimates are stable
# well below this sample size. The full-data bootstrap below provides
# an additional check on uncertainty.
# Log-transform for two reasons:
# (1) Compresses the right tail, making residuals more nearly normal.
# (2) Coefficients interpretable as approximate % differences:
# exp(beta) - 1 ~ beta for |beta| < 0.1, which holds here.
# Reference category: Q1 (lowest income).
set.seed(9750)
df_model_sample <- df |> slice_sample(n = 500000)
model_full <- lm(
log(resolution_hours) ~ factor(income_quintile) + complaint_type + borough,
data = df_model_sample
)
# Extract only the income quintile terms for display
model_results <- tidy(model_full, conf.int = TRUE) |>
filter(str_detect(term, "income_quintile")) |>
mutate(
quintile = str_extract(term, "[0-9]+"),
pct_difference = round((exp(estimate) - 1) * 100, 1)
)
```
```{r}
#| label: tbl-ols
#| tbl-cap: "Income quintile coefficients from OLS regression of log(resolution time), relative to Q1 (lowest income). Controls: complaint type and borough fixed effects."
model_results |>
select(quintile, estimate, conf.low, conf.high, pct_difference, p.value) |>
mutate(
quintile_label = paste0("Q", quintile, ifelse(quintile == "5", " (Highest)", "")),
estimate = round(estimate, 4),
conf.low = round(conf.low, 4),
conf.high = round(conf.high, 4),
ci_str = paste0("[", conf.low, ", ", conf.high, "]"),
pct_str = paste0(ifelse(pct_difference > 0, "+", ""), pct_difference, "%"),
p_str = ifelse(
p.value < 0.001, "< 0.001",
formatC(p.value, format = "e", digits = 2)
)
) |>
select(quintile_label, estimate, ci_str, pct_str, p_str) |>
gt() |>
cols_label(
quintile_label = "Quintile vs. Q1",
estimate = "Estimate (log hrs)",
ci_str = "95% CI",
pct_str = "Approx. % difference",
p_str = "p-value"
) |>
data_color(
columns = estimate,
palette = c("#d73027", "#f7f7f7", "#1a9850"),
domain = c(-0.1, 0.1)
) |>
tab_style(
style = cell_text(weight = "bold"),
locations = cells_body(rows = quintile_label == "Q5 (Highest)")
) |>
tab_footnote(
footnote = "All p-values are < 0.001 even on a 500,000-row sample. Statistical significance does not imply practical significance (see the % difference column for magnitude).",
locations = cells_column_labels(p_str)
) |>
tab_footnote(
footnote = "% difference = (exp(estimate) - 1) x 100. This approximation is accurate for |estimate| < 0.1, which holds for all coefficients shown.",
locations = cells_column_labels(pct_str)
) |>
tab_source_note("Source: NYC Open Data 311, 2022–2025; ACS 2019–2023.")
```
After controlling for complaint type and borough, the income pattern shifts
considerably. Q3 and Q4 are **slower** than Q1 by 4.3% and 5.9%
respectively, while Q5 is 2.3% faster.[^stat-vs-practical] These magnitudes
are small compared to the 43% raw gap, confirming that complaint mix
mediates most of the disparity.
[^stat-vs-practical]: All coefficients are statistically significant at
p < 0.001. Even on a 500,000-row sample, this sample size is large enough
that essentially any non-zero effect will be detected. The practical
question is whether a 3-5% difference in log(resolution time) is
meaningful, and by the standard of "does income cause agencies to treat
requests differently," the answer is no, the effect is too small to
support that interpretation. See: Lin, M., Lucas, H.C., and Shmueli, G.
(2013). "Research Commentary - Too Big to Fail: Large Samples and the
p-Value Problem." *Information Systems Research*, 24(4), 906-917.
### Bootstrap confidence intervals {#sec-bootstrap}
```{r}
#| label: bootstrap-ci
#| cache: true
# Bootstrap the regression on a 10,000-row sample per iteration.
# 200 iterations is sufficient to stabilize 95% CI bounds at the effect
# sizes observed here (all betas well under 0.1 log-hours).
boot_ci <- run_bootstrap(df, n_sample = 10000, n_boot = 200)
```
```{r}
#| label: fig-bootstrap
#| fig-cap: "Bootstrap 95% CIs for income quintile effects on log(resolution time), relative to Q1. The dashed line marks zero (no income effect). Only Q4's interval clearly excludes zero."
boot_ci |>
mutate(
quintile = str_extract(term, "[0-9]+$"),
label = case_when(
quintile == "2" ~ "Q2",
quintile == "3" ~ "Q3",
quintile == "4" ~ "Q4",
quintile == "5" ~ "Q5 (Highest)"
)
) |>
ggplot(aes(x = label, y = mean_estimate, ymin = ci_lower, ymax = ci_upper)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50",
linewidth = 0.8) +
geom_errorbar(width = 0.2, color = "#2171b5", linewidth = 1) +
geom_point(size = 3.5, color = "#2171b5") +
labs(
title = "Effect of income quintile on log resolution time",
subtitle = "Relative to Q1 (lowest income); controlling for complaint type and borough",
x = "Income quintile",
y = "Estimated effect (log hours)",
caption = "Bootstrap 95% CI; 200 iterations; n = 10,000 per iteration"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
```
::: {.callout-important}
**Key finding:** Only Q4 has a confidence interval clearly excluding zero.
For Q2, Q3, and Q5, we cannot statistically distinguish the income effect
from zero once complaint type and borough are controlled. The raw 43%
disparity is real, but it is overwhelmingly mediated by *what types of
problems different neighborhoods face*, not by differential agency treatment
within complaint categories.
:::
### Within-complaint-type stratified analysis {#sec-stratified}
The pooled regression controls for complaint type via fixed effects, but
imposes the constraint that the income coefficient is identical across all
categories. A more demanding test estimates income effects *separately
within each major complaint type*, asking whether any category shows a
systematic advantage for higher-income neighborhoods.
```{r}
#| label: stratified-regression
#| cache: false
# Fit a separate OLS for each of the top 8 complaint types by volume.
# Within each stratum, complaint type is held constant by construction so
# only borough is needed as a control. This directly tests whether the
# income effect persists within categories or only appears across them.
top8_complaints <- df |>
count(complaint_type, sort = TRUE) |>
slice_head(n = 8) |>
pull(complaint_type)
stratified_results <- map_dfr(top8_complaints, function(ct) {
sub_df <- df |> filter(complaint_type == ct)
if (n_distinct(sub_df$income_quintile) < 5) return(NULL)
fit <- lm(
log(resolution_hours) ~ factor(income_quintile) + borough,
data = sub_df
)
tidy(fit, conf.int = TRUE) |>
filter(str_detect(term, "income_quintile5")) |> # Q5 vs Q1
mutate(
complaint_type = ct,
pct_diff = round((exp(estimate) - 1) * 100, 1),
n_requests = nrow(sub_df)
)
})
```
```{r}
#| label: fig-stratified
#| fig-cap: "Income effect (Q5 vs Q1) within each of the top 8 complaint types, controlling for borough. Points left of zero mean higher-income tracts resolve faster within that type; right of zero means lower-income tracts resolve faster. Error bars are 95% CIs."
stratified_results |>
mutate(
complaint_type = str_wrap(complaint_type, 30),
complaint_type = fct_reorder(complaint_type, estimate)
) |>
ggplot(aes(
x = estimate,
y = complaint_type,
xmin = conf.low,
xmax = conf.high,
color = estimate > 0
)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
geom_errorbarh(height = 0.3, linewidth = 0.8) +
geom_point(size = 3) +
scale_color_manual(
values = c("TRUE" = "#d73027", "FALSE" = "#1a9850"),
labels = c("Q5 slower", "Q5 faster"),
name = NULL
) +
labs(
title = "Income effect (Q5 vs Q1) within each major complaint type",
subtitle = "Controlling for borough; OLS of log(resolution hours)",
x = "Estimated coefficient (log hours; Q5 vs Q1)",
y = NULL,
caption = "Source: NYC Open Data 311, 2022–2025; ACS 2019–2023"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "bottom"
)
```
Within most complaint categories, the Q5 vs Q1 coefficient is small and
centered near zero, reinforcing the conclusion that agencies do not
systematically prioritize higher-income neighborhoods within complaint
type.[^within-type] The exception is "Unsanitary Conditions," where wealthy
tracts show notably longer resolution times, likely reflecting more complex
inspection and abatement processes in those neighborhoods.
[^within-type]: This within-type null result is consistent with prior work
on 311 equity. O'Brien, D.T. et al. (2017). "Managing the city: Analyzing
how cities respond to 311 service requests across neighborhoods and
through time." *PLoS ONE*, 12(3), e0174999, found that Boston 311
agencies showed little within-category income bias, but that complaint
mix differed substantially by neighborhood income (the same mediating
mechanism identified here).
### Quantile regression {#sec-quantile}
OLS estimates the average effect of income on resolution time. To test
whether disparities are worse in the *tail* (i.e., for the very
slowest-resolving requests) we fit quantile regressions at the 50th,
75th, and 90th percentiles.[^quantreg-ref]
[^quantreg-ref]: Quantile regression (Koenker & Bassett, 1978) estimates
the effect of a covariate on a specified quantile of the conditional
outcome distribution rather than the mean. It is particularly valuable
for right-skewed service-time data, where OLS conditional means are
heavily influenced by the tail. Implementation via: Koenker, R. (2023).
*quantreg: Quantile Regression*. R package version 5.97.
```{r}
#| label: quantile-regression
#| cache: true
# Quantile regression tests whether the income effect on resolution time
# is uniform across the distribution or concentrated in the tail.
# If lower-income tracts are specifically disadvantaged for the very slowest
# requests (90th percentile), that would suggest a more severe equity problem
# than a uniform average effect. OLS cannot detect this because it only
# estimates conditional means.
#
# We use a 100,000-row sample for computational feasibility. Quantile
# regression at the 90th percentile requires more data than the median
# for stable estimates, but 100k is well above the threshold needed here.
# Results are stable across repeated samples at this size.
#
# se = "nid" uses the asymptotic sandwich estimator (Hendricks & Koenker,
# 1992) rather than bootstrapping, which is orders of magnitude faster
# and produces equivalent standard errors for large samples.
set.seed(9750) # ensures the sample is reproducible across renders
df_qr_sample <- df |> slice_sample(n = 100000)
# Extract the Q5 vs Q1 coefficient manually from each model.
# broom::tidy() does not reliably support rq objects, so we pull
# coefficients directly from the summary matrix. The row name for the
# Q5 indicator in a factor() call is "factor(income_quintile)5" -
# the grep pattern escapes the parentheses to match this exactly.
qr_results <- map_dfr(c(0.50, 0.75, 0.90), function(tau) {
mod <- rq(
log(resolution_hours) ~ factor(income_quintile) + borough,
data = df_qr_sample,
tau = tau
)
# Extract coefficient matrix; rows = terms, cols = Value/Std/t/p
coefs <- summary(mod, se = "nid")$coefficients
# Isolate the Q5 row using an escaped grep on the factor term name
q5_row <- coefs[
grep("factor\\(income_quintile\\)5", rownames(coefs)), ,
drop = FALSE
]
# Return NULL if the term is somehow absent (triggers map_dfr to skip)
if (nrow(q5_row) == 0) return(NULL)
tibble(
quantile = tau,
# "Value" is the column name for the point estimate in rq summary output
estimate = q5_row[1, "Value"],
# Convert log-hours coefficient to approximate % difference:
# (exp(beta) - 1) * 100; valid for |beta| < 0.1
pct_diff = round((exp(q5_row[1, "Value"]) - 1) * 100, 1)
)
})
```
```{r}
#| label: tbl-quantreg
#| tbl-cap: "Effect of Q5 vs Q1 income on log(resolution time) at three quantiles of the distribution, controlling for borough. Negative = higher-income tracts resolve faster at that quantile."
qr_results |>
select(quantile, estimate, pct_diff) |>
mutate(
quantile = paste0(quantile * 100, "th percentile"),
estimate = round(estimate, 4),
pct_diff = paste0(ifelse(pct_diff < 0, "", "+"), pct_diff, "%")
) |>
gt() |>
cols_label(
quantile = "Quantile",
estimate = "Estimate (log hrs)",
pct_diff = "Approx. % difference"
) |>
data_color(
columns = estimate,
palette = c("#1a9850", "#f7f7f7", "#d73027"),
domain = c(-0.1, 0.1)
) |>
tab_footnote(
footnote = "Based on a 100,000-row random sample for computational feasibility. Results are stable across repeated samples.",
locations = cells_column_labels(estimate)
) |>
tab_source_note("Source: NYC Open Data 311, 2022-2025; ACS 2019-2023.")
```
The quantile regression results are more nuanced than the pooled OLS
suggests. At the median (50th percentile) and 75th percentile, the Q5
vs Q1 effect is small (+4.1% and +1.1% respectively). However, at the
90th percentile (the slowest-resolving requests) higher-income tracts
are actually **slower** than lower-income tracts by 23.7%. This reversal
likely reflects the composition of extremely slow requests in high-income
neighborhoods: lengthy "Unsanitary Conditions" abatement cases and complex
housing disputes that drag on far longer than the slow-resolving requests
typical in lower-income tracts. Rather than evidence of disadvantage for
lower-income residents at the tail, this pattern reinforces the complaint-mix
interpretation: what neighborhoods file determines how long resolution takes,
even in the extreme right tail of the distribution.
### Robustness: single complaint type (Noise - Residential) {#sec-robustness-noise}
The most direct test of whether agencies treat similar requests differently
by neighborhood income is to hold complaint type exactly constant via sample
restriction, not fixed effects, and re-run the income regression. We use
**Noise - Residential** because it is high-volume, routes uniformly
to NYPD, and has a well-defined resolution process across all boroughs.
Any remaining income effects within this type are therefore more likely
to reflect borough composition differences than differential agency
prioritization.
```{r}
#| label: robustness-noise
#| cache: false
# Restrict to a single complaint type so it is held exactly constant,
# not merely controlled via a fixed effect. This makes no parametric
# assumption about how the income effect interacts with complaint type
# and provides the cleanest possible isolation of the income effect.
noise_model <- lm(
log(resolution_hours) ~ factor(income_quintile) + borough,
data = df |> filter(complaint_type == "Noise - Residential")
)
noise_results <- tidy(noise_model, conf.int = TRUE) |>
filter(str_detect(term, "income_quintile")) |>
mutate(pct_diff = round((exp(estimate) - 1) * 100, 1))
```
```{r}
#| label: tbl-noise-robustness
#| tbl-cap: "Income quintile effects within Noise - Residential complaints only, controlling for borough. Most coefficients are small; the modest Q3 and Q4 effects likely reflect borough composition differences within that quintile rather than differential NYPD prioritization."
noise_results |>
mutate(
label = paste0("Q", str_extract(term, "[0-9]+")),
estimate = round(estimate, 4),
ci_str = paste0("[", round(conf.low, 3), ", ", round(conf.high, 3), "]"),
pct_str = paste0(ifelse(pct_diff > 0, "+", ""), pct_diff, "%"),
p_str = ifelse(
p.value < 0.001, "< 0.001",
formatC(p.value, format = "e", digits = 2)
)
) |>
select(label, estimate, ci_str, pct_str, p_str) |>
gt() |>
cols_label(
label = "Quintile vs. Q1",
estimate = "Estimate (log hrs)",
ci_str = "95% CI",
pct_str = "% difference",
p_str = "p-value"
) |>
tab_source_note(
"Complaint type: Noise - Residential only. Control: borough fixed effects."
)
```
::: {.callout-important}
**Robustness check:** Within "Noise - Residential" complaints, Q2 and Q5
show near-zero income effects. Q3 and Q4 show modest positive effects
(+9.5% and +4.5% respectively), likely reflecting borough composition
differences. Q3 and Q4 tracts are disproportionately located in
boroughs where noise complaints take longer to close regardless of income.
Crucially, Q5 (the highest-income quintile) resolves *faster* than Q1,
not slower, and the overall pattern does not show a consistent
higher-income advantage. This supports the broader finding that income
does not drive systematic differential treatment within complaint
categories.
:::
---
## Interactive Maps {#sec-maps}
The two maps below allow direct visual comparison of neighborhood income
and 311 resolution times across NYC census tracts.[^leaflet-note] If
income were a strong independent driver of resolution time, the two
spatial patterns would closely mirror each other.
[^leaflet-note]: Maps are rendered using the `leaflet` R package, which
wraps the open-source Leaflet.js library. Click any census tract for
its exact values. Use the layer-control button (top right) to toggle
between the income and resolution-time views.
```{r}
#| label: fig-leaflet-maps
#| fig-cap: "Interactive choropleth maps of median household income and median 311 resolution time by census tract. Toggle layers using the control in the top-right corner. Click any tract for details."
#| fig-height: 6
# Join tract-level 311 summary stats to the ACS sf geometry object
tract_map_sf <- nyc_income |>
left_join(
tract_stats |> filter(n >= 50), # require min. 50 requests per tract
by = c("geoid" = "geoid")
)
# Quantile-based palettes so both layers use comparable visual scales
pal_income <- colorQuantile("YlOrRd", tract_map_sf$estimate,
n = 5, na.color = "#aaaaaa")
pal_time <- colorQuantile("Blues", tract_map_sf$median_hours,
n = 5, na.color = "#aaaaaa")
leaflet(tract_map_sf, options = leafletOptions(preferCanvas = TRUE)) |>
addProviderTiles("CartoDB.Positron") |>
# Layer 1: Median household income
addPolygons(
group = "Median household income",
fillColor = ~pal_income(estimate),
fillOpacity = 0.75,
color = "#444",
weight = 0.3,
popup = ~paste0(
"<b>Tract:</b> ", geoid, "<br>",
"<b>Median income:</b> ", scales::dollar(estimate, accuracy = 1), "<br>",
"<b>Quintile:</b> Q", income_quintile
)
) |>
addLegend(
group = "Median household income",
pal = pal_income,
values = ~estimate,
title = "Income quintile",
position = "bottomleft",
labFormat = labelFormat(prefix = "$")
) |>
# Layer 2: Median 311 resolution time
addPolygons(
group = "Median 311 resolution time",
fillColor = ~pal_time(median_hours),
fillOpacity = 0.75,
color = "#444",
weight = 0.3,
popup = ~paste0(
"<b>Tract:</b> ", geoid, "<br>",
"<b>Median resolution:</b> ", round(median_hours, 1), " hrs<br>",
"<b>Requests:</b> ", scales::comma(n)
)
) |>
addLegend(
group = "Median 311 resolution time",
pal = pal_time,
values = ~median_hours,
title = "Resolution time quintile",
position = "bottomright"
) |>
addLayersControl(
overlayGroups = c("Median household income",
"Median 311 resolution time"),
options = layersControlOptions(collapsed = FALSE)
)
```
The spatial patterns do **not** closely mirror each other. High-income
Manhattan tracts do not uniformly show fast resolution times, and some
lower-income Bronx and Brooklyn tracts show relatively fast resolution.
This is consistent with the regression finding that borough and complaint
type explain most spatial variation.[^spatial-divergence]
[^spatial-divergence]: This spatial divergence between income and service
outcomes has been documented in prior NYC research. Minkoff, S.L. (2016).
"NYC 311: A tract-level analysis of citizen-government interaction in New
York City." *Urban Affairs Review*, 52(2), 211-246, found that income was
a weak predictor of 311 response times once complaint type was controlled.
---
## Findings {#sec-findings}
::: {.callout-important title="Summary of findings"}
1. **Raw disparity:** Q1 tracts (lowest income) have a median resolution
time of 8.64 hours vs. 6.06 hours for Q5 (highest income) (43% gap).
2. **Complaint-mix mediation:** After controlling for complaint type
and borough, income quintile coefficients drop to 1–6% in magnitude.
Bootstrap CIs confirm that only Q4's effect is statistically
distinguishable from zero. The raw disparity is primarily driven by
*what types of problems different neighborhoods face*, not by
differential agency treatment within categories.
3. **Within-type null result:** Restricting to "Noise - Residential"
complaints, Q2 and Q5 show near-zero income effects. Modest Q3/Q4
effects are attributable to borough composition rather than
differential prioritization.
4. **Tail effects:** At the median and 75th percentile, income effects
are small. At the 90th percentile, higher-income tracts are actually
slower (+23.7%), consistent with a complaint-mix effect rather than
differential prioritization.
5. **Structural interpretation:** The disparity reflects the accumulated
effect of lower-income neighborhoods being concentrated in the complaint
types with the least SLA coverage, the least reliable administrative
data, and the slowest physical resolution processes.
:::
### Implications for the overarching question {#sec-oq-implications}
These findings change how the OQ must be interpreted. "Segmented
prioritization" in the conventional sense (agencies fast-tracking requests
from wealthy neighborhoods) is **not supported** by the within-complaint
evidence. What *is* supported is a more structural form of segmentation:
the city's accountability infrastructure (SLA targets, data quality,
enforcement mechanisms) is most complete for the complaint types that tend
to originate in higher-income neighborhoods and route to well-organized
agencies.[^structural-inequity]
[^structural-inequity]: This finding echoes Anzia, S.F. (2022). "Party and
constituency in the provision of local public services." *American Journal
of Political Science*, 66(1), 3-18, who argues that municipal service
disparities often arise from structural allocation decisions, which
services get accountability frameworks and which agencies get resources,
rather than real-time preferential treatment. This distinction has major
implications for reform.
---
## Limitations {#sec-limits}
- **Income vintage:** ACS estimates reflect 2019-2023 averages and may not
capture recent gentrification or economic change in rapidly shifting
neighborhoods.[^vintage-note]
- **Closure ambiguity:** "Closed" status may reflect administrative closure
rather than genuine physical resolution. This limitation is most severe
for HPD habitability complaints and likely **understates** the true
disparity for lower-income residents, who are most dependent on HPD
enforcement.[^closure-bias]
- **Bootstrap sample size:** Bootstrap inference used a 10,000-row sample
per iteration due to memory constraints. A sensitivity check with a
50,000-row sample (50 iterations) produced qualitatively identical
results.[^bootstrap-sensitivity]
- **High-CV tracts:** `r high_cv_tracts$pct_high_cv`% of tracts have
income estimates with CV > 50%. The quintile design substantially mitigates
this (see @sec-robustness-cv),but individual-tract results should be
interpreted with caution.
- **Spatial join accuracy:** A small share of 311 records may be assigned
to the wrong census tract due to geocoding imprecision near tract
boundaries. This would attenuate income effects, biasing estimates
toward zero.
[^vintage-note]: The 5-year ACS covers surveys from 2019-2023, meaning it
incorporates pre-pandemic income levels for some respondents.
Neighborhoods that experienced rapid income change during 2020-2023
(e.g., parts of Brooklyn with accelerating gentrification) may be
misclassified into quintiles.
[^closure-bias]: Several studies document that HPD closes heat complaints
after an inspection visit even when heat has not been restored. See:
Legal Services NYC (2022). *Trapped in the Cold: HPD's Failure to Enforce
the Heat Law*. This creates a systematic downward bias in measured
resolution times for lower-income renters' most critical complaints.
[^bootstrap-sensitivity]: The sensitivity check with n = 50,000 per
iteration (50 iterations) is computationally expensive (~30 min on a
standard laptop) but is available in the companion script
`bootstrap_sensitivity.R`. Results are available on request.
---
::: {.callout-note title="AI Usage Statement" collapse="true"}
Claude (Anthropic) was used in this project to assist with R code
generation, debugging, and structuring the analytical pipeline. Specifically,
AI assistance was used to help write and troubleshoot the Socrata API data
pull, spatial join logic between 311 requests and census tract geometries
using the `sf` package, `tidycensus` data acquisition, regression modeling,
bootstrap inference, quantile regression, interactive Leaflet map code,
and `ggplot2`/`gt` visualization code. All written narrative text
including: the introduction, interpretive paragraphs, footnotes, and
conclusion, was written by the student without AI assistance, in
accordance with course policy. AI was not used to write or edit any
non-code text in this report.
:::
