Park, Proximity, and Perception: Rat Sightings and 311 Reporting Behavior in New York City

I. Introduction

Rats have been a persistent and highly visible problem in New York City for centuries. The infamous “Pizza Rat” offered viral evidence of what New Yorkers encounter regularly. While the moment was humorous to the outside, for city residents, these rodents are more than just a nuisance. They pose public health risks, contaminate food, and damage property. Rats are commonly found in alleys, parks, subway stations, and, at times, inside apartment buildings.

In 2003, NYC 311 was launched to divert non-emergency calls from 911. Since then this hotline number has evolved into a comprehensive reporting system, allowing residents to submit service requests via phone, the 311 website, mobile apps, text messages, and even social media. These reports cover a wide range of issues, including noise complaints, infrastructure concerns, and even rodent sightings.

Rodent-related calls provide a rich source of information on the spatial and temporal patterns of reported rat activity across New York City. By analyzing these calls alongside urban and socio-economic features, we can explore the factors associated with rat presence. At the same time, it is important to consider whether 311 rodents calls primarily reflect underlying rat activity or instead are shaped by human reporting behavior.

Our broader analysis examines several features, including urban infrastructure, socio-economic factors, environmental conditions, language and demographics, and public space and land use.

This section of the analysis focuses specifically on public parks and addresses the question:

Do parks lead to more rat sightings due to increased visibility, or fewer sightings due to differences in land use and human interaction patterns?

II. Data Acquisition and Preparation

This analysis draws on three data sources from NYC Open Data:

1. Parks Properties
2. 311 Rodent Complaints
3. Neighborhood Tabulation Area (NTA) Boundaries

The Parks Properties dataset provides information on properties managed partially or entirely by NYC Parks. It includes a variety of managed spaces beyond traditional parks and green areas, such as playgrounds, plazas, and recreational facilities. While the data captures all public park locations across the city, including acreage, park name, and neighborhood identifiers, not all entries represent natural green space, which is an important distinction for this analysis.

#Downloading the NYC Parks Data
#Create directory folder if it doesn't exist
if(!dir.exists(file.path("data", "final"))){
  dir.create(file.path("data", "final"), showWarnings = FALSE, recursive = TRUE)
}


# Define paths and Parks URL
PARKS_FILE <- file.path("data","final","parks.geojson")
PARKS_URL <- "https://data.cityofnewyork.us/resource/enfh-gkve.geojson"
                      

if(!file.exists(PARKS_FILE)) {
  download.file(PARKS_URL, destfile = PARKS_FILE, mode="wb")
  message("Downloaded PARKS FILE.")
} else {
  message("PARKS FILE already exists; skipping download.")
}

parks_data <- st_read(PARKS_FILE, quiet=TRUE)
head(parks_data)

The 311 Rodent Complaints dataset includes records of rodent-related service requests submitted since 2003. Each report includes the date the service request was created,its location, and the type of complaint.

#Downloading the full 311 rats dataset from NYC Open Data
#Looping through pages of the API and saving locally
if (!dir.exists(file.path("data", "final"))){
  dir.create(file.path("data", "final"), showWarnings = FALSE, recursive = TRUE)
}

#Define API endpoint
base_url <-"https://data.cityofnewyork.us/resource/cvf2-zn8s.geojson"
limit <- 50000 #number of rows per request
offset <- 0   #start from the beginning
page <-1 #page counter
all_data<-list()

repeat {
  file_path <- file.path("data/final", paste0("311rodents_", page, ".geojson"))
  
  if (!file.exists(file_path)) {
    req <- request(base_url) |>
      req_url_query(`$limit` = limit, `$offset` = offset)
    
    resp <- req |> req_perform()
    resp_check_status(resp)
    
    writeBin(resp_body_raw(resp), file_path)
    message("Downloaded page ", page)
  } else {
    message("File already exists: page ", page, "; skipping download")
  }

  # Read page into sf
  page_data <- st_read(file_path, quiet = TRUE)
  all_data[[page]] <- page_data

  # If fewer rows than limit, stop
  if (nrow(page_data) < limit) break

  page <- page + 1
  offset <- offset + limit
}

rats_full <- bind_rows(all_data)

To enable neighborhood-level analysis, we also obtained Neighborhood Tabulation Area (NTA) boundaries from NYC Open Data. These were used to assign each park and 311 rodent complaint to an NTA, allowing aggregation and comparison of park characteristics and rodent sightings across neighborhoods.

# Downloads and prepares NTA polygons 
nta_path <- "data/final/nta_2020.geojson" 
 
# Reads the saved NTA GeoJSON as an sf object (polygons). 
# st_transform(4326) ensures everything is in WGS84 lat/lon (EPSG:4326),  
# which matches the coordinates used in the rat 311 data. 
nta <- st_read(nta_path, quiet = TRUE) |> 
  st_transform(4326) 

All datasets were downloaded directly from their respective NYC Open Data API endpoints and saved locally as GeoJSON files within a structured data directory. This approach ensures reproducibility, minimizes repeated API calls, and allows verification that the data was successfully and completely downloaded before analysis.

The 311 dataset was initially explored using a smaller subset of observations pulled via the NYC Open Data API. The full dataset was subsequently downloaded in pages of 50,000 rows per request to create a complete local copy in a responsible and efficient manner.

Once downloaded, all datasets were parsed, cleaned and processed. Park acreage values were converted to numeric format, duplicate rows were removed, and both the parks and 311 complaint datasets were prepared for spatial joins with NTA boundaries. Some park geometries were initially invalid or consisted of multipolygons features. These issues were corrected using standard spatial cleaning procedures to ensure reliable spatial joins and aggregation.

During exploratory mapping, certain well-known large parks, such as Central Park, did not appear explicitly in the Parks Properties dataset used for this analysis. While the dataset contains a broad list of parks managed by NYC Parks, representation may vary due to data reporting or classification practices. Aside from these naming and representation discrepancies, all available records in the dataset were successfully processed and included in the analysis.

Together, these steps enabled the integration of multiple data sources for a robust neighborhood-level analysis.

#Attach NTA info to full 311 rat dataset

files <- list.files("data/final", pattern = "311rodents_.*\\.geojson$", full.names = TRUE)
all_pages <- lapply(files, st_read, quiet = TRUE)
rats_full <- bind_rows(all_pages)

#Convert to sf points for spatial operations
rats_sf <- rats_full |>
  mutate(
    lon = as.numeric(longitude),
    lat = as.numeric(latitude))|>
  drop_na(lon, lat) |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

#Read NTA polygons
nta <- st_read("data/final/nta_2020.geojson", quiet = TRUE) |>
  st_transform(4326)

#Spatial join: add NTA info to 311 rat reports
rats_with_nta <- st_join(rats_sf, nta[, c("NTA2020", "NTAName", "BoroName")])

#Drop geometry if needed to create a DF or tibble
rats_with_nta_df <- rats_with_nta %>% st_drop_geometry()

#Quick check
glimpse(rats_with_nta_df)

#Attach NTA info to parks dataset
parks_sf <- st_read("data/final/parks.geojson", quiet = TRUE) |>
  st_transform(4326)|>
   st_make_valid()

#Make sure NTA polygons are ready
nta <- st_read("data/final/nta_2020.geojson", quiet = TRUE) %>%
  st_transform(4326)
  
#Spatial join to add NTA info to parks (for mapping)
parks_with_nta_sf <- st_join(
  parks_sf,
  nta |> select(NTA2020, NTAName, BoroName),
  left = TRUE
)

#Clean parks_with_nta_df
parks_with_nta_df <- parks_with_nta_sf |>
  st_drop_geometry() |> # drop geometry first
  select(NTAName, BoroName, signname, acres) |>  # keep only relevant columns
  mutate(acres_numeric = as.numeric(gsub(",", "", acres))) |>  # convert acres to numeric
  distinct(NTAName, BoroName, signname, acres_numeric, .keep_all = TRUE)  # remove duplicate rows

#Check the result
glimpse(parks_with_nta_df)

III. Exploratory Analysis

Spotlight on NYC Parks: Distribution, Size, and Categories

New York has over 1,700 parks, playgrounds, and recreation facilities managed by NYC Parks, covering more than 30,000 acres ¹. The map below shows the spatial distribution of the parks included in the dataset, which represents a subset of approximately 1,000 of these parks.

leaflet() |>
  addProviderTiles("CartoDB.Positron") %>%  
  addPolygons(
    data = parks,
    fillColor = "forestgreen",
    weight = 1,
    opacity = 1,
    fillOpacity = 0.5,
    color = "darkgreen",
    group = "Parks",
    label = ~signname  # show the park name on hover
  ) |>
  addLayersControl(
    overlayGroups = c("Parks"),
    options = layersControlOptions(collapsed = FALSE)
  )

In the NYC Parks dataset, parks are represented as spatial features rather than one record per park name. That means a single park may consist of multiple polygons. The table below is a borough-level summary of NYC parks, distinguishing between total park features and unique park names.

kable(
  borough_summary,
  col.names = c("Borough", "Park Features", "Unique Parks", "Total Acres", "Mean Park Size (acres)"),
  align = "c",
  caption = "Summary of NYC Parks by Borough"
)

Summary of NYC Parks by Borough

Borough	Park Features	Unique Parks	Total Acres	Mean Park Size (acres)
Bronx	141	124	4014.9	28.47
Staten Island	85	77	2130.9	25.07
Queens	266	243	1374.2	5.17
Brooklyn	363	341	853.3	2.35
Manhattan	145	140	261.1	1.80

New York City parks vary widely in type and size. To better understand how park space may influence rat sightings, the focus is on park types by total acreage. The dataset classifies parks into 19 different categories. Although only a few flagship parks exist, they occupy a large portion of the city’s parkland. In contrast, smaller types like playgrounds, lots, and strips are numerous, but they occupy relatively less acreage. This suggests they may have a smaller ecological footprint. In fact, the top five park categories (Flagship Parks, Community Parks, Nature Areas, Neighborhood Parks, and Recreational Field/Courts) account for 7,682 acres and represent over half of the total acreage in the dataset.

#Plotting park type and total acreage
parks_by_type <- parks |>
  mutate(acres_numeric = as.numeric(gsub(",", "", acres))) |>
  st_drop_geometry() |>
  group_by(typecategory) |>
  summarise(
    total_acres = sum(acres_numeric, na.rm = TRUE),
    .groups = "drop"
  ) |>
  arrange(desc(total_acres))

ggplot(parks_by_type, aes(x = reorder(typecategory, -total_acres), y = total_acres)) +
  geom_col(fill = "forestgreen", alpha = 0.5) +
  scale_y_continuous(labels = comma) +
  labs(
    x = "Park Type",
    y = "Total Acres",
    title = "Total Park Acreage by Type in NYC"
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.title = element_text(face = "bold"),
    panel.grid.major.x = element_blank(), 
    panel.grid.minor = element_blank()
  )

# Join NTA summary back to spatial polygons for mapping
nta_sf <- nta |>
  left_join(
    nta_summary,
    by = c("NTAName", "BoroName")
  )

Rat Activity Across NYC Neighborhoods

The choropleth map visualizes the total 311 rat complaints across New York City at the NTA level. Each NTA polygon is shaded according to the number of reported incidents, with darker shades indicating higher volumes of complaints. While parks are distributed throughout the city,NTAs with the highest volume of rat complaints tend to be densely built neighborhoods rather than areas dominated by large park spaces.

#Creating rat complaints by NTA with Parks Overlay
ggplot() +
  # Choropleth: NTAs colored by total rat complaints
  geom_sf(data = nta_sf, aes(fill = total_calls), color = NA) +
  
  # Smooth blue gradient with lighter ends
  scale_fill_gradientn(
  colors = c("#E0F4FA", "#C0E4F8", "#A0D4F6", "#80C4F4", "#60B4F2", "#4094E0", "#1E3F8B"),  # light → dark
  trans = "log",
  na.value = "grey90",
  labels = function(x) round(x)
)+
  # Overlay parks as contrasting red polygons
  geom_sf(data = parks_sf, fill = "green", color = "darkgreen", alpha = 0.7) +
  
  
  labs(
    title = "Rat Complaints Across NYC Neighborhoods and Parks",
    fill = "311 Rat Reports"
  ) +
   theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 10), 
    legend.title = element_text(face = "bold", size = 10),
    legend.text = element_text(size = 10)
  )

To explore whether these spatial patterns are related to park coverage, the relationship between total park acreage and 311 rat complaints at the NTA level was analysed.

The scatterplot below shows the relationship between total park acreage and 311 rat calls per NTA. A log scale was applied to the x-axis to account for the large variability in total park acreage across neighborhoods. Most NTAs fall roughly between 0.5 and 500 acres, with a few extreme cases extending beyond 500 acres.While the points are concentrated among low-to-medium park acreage NTAs, neighborhoods with very large park acreage do not show remarkably higher rat calls. Most points are concentrated among low-to-medium park acreage NTAs, indicating the trend is influenced by a few large parks.

ggplot(nta_summary, aes(x = total_acres, y = total_calls)) +
  geom_point(alpha = 0.6, color = "steelblue") +
  geom_smooth(method = "lm", se = TRUE, color = "darkred") +
  scale_x_continuous(
    trans = "log10",
    breaks = c(0.5, 1, 5, 10, 50, 100, 500, 1000, 2000),
    labels = function(x) {
      ifelse(x < 1, x, scales::comma(x, accuracy = 1))
    }
  ) +
  scale_y_continuous(labels = comma) +
  labs(
    x = "Total Park Acreage (log scale)",
    y = "Total Rat Calls",
    title = "Rat Calls vs. Total Park Acreage per NTA",
    subtitle = "Log scale applied to acreage to account for large variability across neighborhoods"
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size =10),
    plot.subtitle = element_text(hjust = 0.5),
    axis.title = element_text(face = "bold")
  )

Statistical analysis supports this observed pattern. A Spearman rank correlation between log-transformed park acreage and total rat calls is ρ = 0.243 (p < 0.001), indicating a weak but statistically significant association. Linear regression shows that each one-unit increase in log-acres is associated with roughly 316 additional rat calls (p = 0.006), although the model explains only ~2.8% of the variance (R² = 0.028). Together, these results suggest that park acreage has a measurable, but modest effect on reported rat activity, and other neighborhood factors are likely to play a much larger role.

#Correlation and regression analysis for rat calls vs. park acreage
# Log-transform total_acres (add a small constant to avoid log(0))
nta_summary <- nta_summary %>%
  mutate(log_acres = log10(total_acres + 0.1))

# Correlation test (Spearman or Pearson)
cor_test <- cor.test(nta_summary$log_acres, nta_summary$total_calls, method = "spearman")
cor_test

# Linear regression
lm_model <- lm(total_calls ~ log_acres, data = nta_summary)
summary(lm_model)

# Create a data frame with key statistics
# Create the stats table
stats_table <- data.frame(
  Statistic = c(
    "Spearman ρ (log-acres vs. total calls)",
    "Linear regression coefficient (log-acres)",
    "R-squared"
  ),
  Estimate = c(
    0.243,
    316,
    0.028
  ),
  `p-value` = c(
    "< 0.001",
    0.006,
    NA
  ),
  Interpretation = c(
    "Weak positive correlation",
    "Each unit increase in log-acres associated with ~316 additional rat calls",
    "Model explains ~2.8% of variance in total rat calls"
  ),
  stringsAsFactors = FALSE
)

# Display the table centered
kable(
  stats_table,
  caption = "Key Statistical Results for Park Acreage and Rat Calls per NTA",
  align = "c",
  digits = 3
)

Key Statistical Results for Park Acreage and Rat Calls per NTA

Statistic	Estimate	p.value	Interpretation
Spearman ρ (log-acres vs. total calls)	0.243	< 0.001	Weak positive correlation
Linear regression coefficient (log-acres)	316.000	0.006	Each unit increase in log-acres associated with ~316 additional rat calls
R-squared	0.028	NA	Model explains ~2.8% of variance in total rat calls

#Top NTAs with most rat calls
nta_summary |>
  arrange(desc(total_calls)) |>
  slice(1:10) |>
  select(NTAName, BoroName, total_calls, total_parks, total_acres)

# Average calls, parks, and acreage by borough
nta_summary |>
  group_by(BoroName) |>
  summarise(
    avg_calls = mean(total_calls),
    avg_parks = mean(total_parks),
    avg_acres = mean(total_acres)
  )

borough_avg_summary <- nta_summary |>
  group_by(BoroName) |>
  summarise(
    avg_rat_calls = round(mean(total_calls), 0),
    avg_num_parks = round(mean(total_parks), 2),
    avg_park_acres = round(mean(total_acres), 2),
    .groups = "drop"
  )

The table below summarizes average rat calls, park counts, and park acreage by borough. At the borough-level, higher rat call volumes are concentrated in Manhattan and Brooklyn, which have relatively limited park acreage, whereas the Bronx and Staten Island, which have larger average park acreage, report fewer rat incidents. This suggests that park acreage alone does not sufficiently explain rat-reporting patterns across New York City.

knitr::kable(
  borough_avg_summary,
  col.names = c(
    "Borough",
    "Average Rat Calls (311)",
    "Average Number of Parks",
    "Average Park Acreage"
  ),
  align = "c",
  caption = "Average Rat Calls and Park Characteristics by Borough"
)

Average Rat Calls and Park Characteristics by Borough

Borough	Average Rat Calls (311)	Average Number of Parks	Average Park Acreage
Bronx	1819	3.12	311.14
Brooklyn	2479	5.35	16.67
Manhattan	3137	4.08	7.99
Queens	1028	3.29	23.95
Staten Island	972	3.87	100.82

IV. Rat Sightings and Urban Context Over Time (2018-2025)

#Calls by NTA
rat_calls_by_nta <- rats_2018_2025_final |>
  group_by(NTAName, BoroName) |>
  summarise(total_calls = n(), .groups = "drop")

#Calls by year (baseline framing)
rats_2018_2025_final |>
  mutate(year = lubridate::year(created_date)) |>
  group_by(year) |>
  summarise(total_calls = n())

The plot shows that rat complaints in the top 10 NTAs have generally increased between 2018-2023 but has seen a small decline since. Still, some neighborhoods have experienced significant fluctuations. The NTAs with the highest rat call volumes tend to have fewer total park acreage, which suggests that lower park coverage may coincide with higher rat visibility or prevalence of rats in dense urban areas. Conversely, NTAs with larger or more numerous parks generally report fewer calls, which indicates that park presence alone does not drive sightings. Overall, rat activity is unevenly distributed, highlighting the interaction between park space and human reporting behavior in understanding 311 complaint patterns.

top_ntas <- rats_2018_2025_final %>%
  group_by(NTAName) %>%
  summarise(total = n()) %>%
  slice_max(total, n = 10) %>%
  pull(NTAName)

rat_calls_yearly_top <- rats_2018_2025_final %>%
  mutate(year = lubridate::year(created_date)) %>%
  filter(NTAName %in% top_ntas) %>%
  group_by(NTAName, year) %>%
  summarise(total_calls = n(), .groups = "drop")

ggplot(rat_calls_yearly_top, aes(x = year, y = total_calls, color = NTAName)) +
  geom_line(aes(linetype = NTAName), size = 1.2) +
  geom_point(aes(shape = NTAName), size = 3) +
  scale_color_brewer(palette = "Paired") +
  scale_y_continuous(labels = comma) +
  scale_x_continuous(breaks = 2018:2025) +
  labs(
    title = "Annual Rat Reports Across NYC's 10 Most Affected NTAs (2018-2025)",
    x = "Year",
    y = "Total Rat Calls",
    color = "NTA",
    linetype = "NTA",
    shape = "NTA"
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 10),
    axis.title = element_text(face = "bold"),
    axis.text = element_text(color = "black"),
    legend.position = "right",
    legend.title = element_text(face = "bold"),
    panel.grid.major = element_line(color = "gray90"),
    panel.grid.minor = element_blank()
  )

#Compute distance of each rat report to nearest park
st_crs(rats_2018_2025) <- st_crs(parks_sf)
nearest_idx <- st_nearest_feature(rats_2018_2025, parks_sf)

# Compute the distance to the nearest park
rats_2018_2025 <- rats_2018_2025 |>
  mutate(min_distance_to_park = st_distance(geometry, parks_sf[nearest_idx, ], by_element = TRUE))

# Dropping geometry for tables or regression
rats_2018_2025_df <- rats_2018_2025 |> st_drop_geometry()

# Convert distance to feet
rats_2018_2025 <- rats_2018_2025 %>%
  mutate(min_distance_ft_num = set_units(min_distance_to_park, "ft")|> drop_units())

# Extract year from created_date
rats_2018_2025 <- rats_2018_2025 %>%
  mutate(year = lubridate::year(created_date),
    min_distance_ft_num = set_units(min_distance_to_park, "ft") |> drop_units()
  ) |>
  filter(is.finite(min_distance_ft_num))

This ridge plot visualizes the distribution of rat reports by their distance to the nearest path between 2018 and 2025. Each ridge represents the density of rat complaints at a given distance for that year, with taller, darker regions indicating a higher number of reports. Overall, the plot shows that most rat reports tend to occur relatively close to parks, with density decreasing at greater distances. Comparing across all years, the spatial distribution remains fairly similar. While 2018 has the fewest reports and 2020 saw a dip due to the pandemic, the overall patterns of proximity to parks remained largely unchanged.

#Extract and sort unique years to ensure consisten ordering in the ridge plot
years <- sort(unique(rats_2018_2025$year))
custom_blue <- colorRampPalette(c("#6baed6", "#08306b"))(length(years))

#Showing distribution of rat report distances by year
ggplot(rats_2018_2025, aes(x = min_distance_ft_num, y = factor(year), fill = factor(year))) +
  geom_density_ridges(alpha = 0.8, scale = 1.2, color = "white", size = 0.3) +
  geom_vline(
    xintercept = median(rats_2018_2025$min_distance_ft_num, na.rm = TRUE),
    linetype = "dashed",
    color = "red",
    size = 0.5
  ) +
  scale_x_continuous(
    limits = c(0, 2000),
    breaks = seq(0, 2000, 200)
  ) +
  scale_fill_manual(values = custom_blue) +
  labs(
    x = "Distance to Nearest Park (ft)",
    y = "Year",
    title = "Distribution of Rat Reports by Distance to Parks (2018–2025)",
    subtitle = "Ridges show density of rat reports relative to park proximity; dashed line indicates median distance"
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 10),
    plot.subtitle = element_text(hjust = 0.5, size = 9),
    axis.title = element_text(face = "bold", size = 10),
    axis.text.y = element_text(face = "bold", size = 10),
    legend.position = "none"
  )

# Categorize distances
# Categorize distance from nearest park
rats_2018_2025_cat <- rats_2018_2025 |>
  mutate(distance_category = case_when(
    min_distance_ft_num < 100 ~ "< 100 ft",
    min_distance_ft_num < 300 ~ "100–300 ft",
    min_distance_ft_num < 500 ~ "300–500 ft",
    TRUE ~ "> 500 ft"
  ))

While the ridge plot shows that most rat reports occur within approximately 0 and just over 800 feet of parks, which highlights that proximity to parks is a key factor in localized rat activity, the stacked bar chart shows that reports occurring less than 100 feet from parks consistently represents the smallest share citywide. This reflects the difference between local density versus citywide totals. Near each park, rat sightings are concentrated relative to surrounding areas. However, because most of the city’s land area is farther than 100 feet from any park, the majority of reports overall come from the more distant locations. Therefore, proximity matters at the neighborhood level, but citywide, areas farther from parks contribute to the largest number of reports due to their greater spatial extent.

#Analyzing proportion as a stronger answer to whether 311 calls measure real rat activity of human reporting behavior
distance_year_counts <- rats_2018_2025_cat |>
  group_by(year, distance_category) |>
  summarise(n_reports = n(), .groups = "drop")

distance_year_counts <- distance_year_counts |>
  mutate(distance_category = factor(
    distance_category,
    levels = c("< 100 ft", "100–300 ft", "300–500 ft", "> 500 ft")
  ))

# Define custom green palette for distance categories
green_palette <- c(
  "< 100 ft" = "#a8ddb5",  
  "100–300 ft" = "#7bccc4",                 
  "300–500 ft" = "#43a2ca",                  
  "> 500 ft" = "#0868ac"        
)

ggplot(distance_year_counts,
       aes(x = factor(year), y = n_reports, fill = distance_category)) +
  geom_col(position = "fill", width = 0.7, color = "white") +
  scale_y_continuous(labels = percent) +
  scale_fill_manual(values = green_palette) +
  labs(
    x = "Year",
    y = "Share of Reports",
    fill = "Distance to Nearest Park",
    title = "Share of Rat Reports by Distance to Nearest Park (2018–2025)"
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title = element_text(face = "bold", size = 10, hjust = 0.5),
    axis.title = element_text(face = "bold", size = 10),
    axis.text = element_text(face = "bold", size = 10),
    legend.title = element_text(face = "bold", size = 10),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank()
  )

A bootstrap analysis was conducted to compare the number of 311 rat reports by distance from parks. The results show a clear gradient across distance categories. Neighborhoods closest to parks (<100 feet) have an average of 1,013 reports (95% CI: 743-1,259), those 100-300 feet away average 2,884 reports (95% CI: 2,062-3,572), those 300-500 feet away average 4,643 reports (95% CI: 3,244-5,890), and those located more than 500 feet away from a park average 22,236 reports (95% CI: 16,190-27,454).

These results indicate that reported rat activity increases with distance from parks. Neighborhoods farther from parks tend to generate substantially more 311 rat complaints, which suggests that sightings are concentrated in more densely urbanized areas rather than immediately adjacent to park spaces.

# Remove geometry
distance_no_geom <- distance_year_counts %>% st_set_geometry(NULL)

# Bootstrap function
bootstrap_mean <- function(x, R = 1000) {
  boot_obj <- boot::boot(x, statistic = function(d, i) mean(d[i]), R = R)
  data.frame(
    Mean = mean(boot_obj$t),
    CI_lower = quantile(boot_obj$t, 0.025),
    CI_upper = quantile(boot_obj$t, 0.975)
  )
}

# Apply bootstrap per distance bin
bootstrap_results <- distance_no_geom %>%
  group_by(distance_category) %>%
  summarize(bootstrap_mean(n_reports), .groups = "drop")

bootstrap_results

# Bootstrap results
bootstrap_df <- tibble::tibble(
  distance_category = c("<100 ft", "100–300 ft", "300–500 ft", ">500 ft"),
  mean_reports = c(1013, 2884, 4643, 22236),
  CI_lower = c(743, 2062, 3244, 16190),
  CI_upper = c(1259, 3572, 5890, 27454)
)

# Format numbers with commas
bootstrap_df <- bootstrap_df %>%
  mutate(
    mean_reports = formatC(mean_reports, format = "d", big.mark = ","),
    CI_lower = formatC(CI_lower, format = "d", big.mark = ","),
    CI_upper = formatC(CI_upper, format = "d", big.mark = ",")
  )

# Print table with centered columns in RMarkdown
library(knitr)

knitr::kable(
  bootstrap_df,
  col.names = c("Distance from Park", "Mean Reports", "95% CI Lower", "95% CI Upper"),
  align = "c",
  caption = "Bootstrap Estimates of 311 Rat Reports by Distance from Parks"
)

Bootstrap Estimates of 311 Rat Reports by Distance from Parks

Distance from Park	Mean Reports	95% CI Lower	95% CI Upper
<100 ft	1,013	743	1,259
100–300 ft	2,884	2,062	3,572
300–500 ft	4,643	3,244	5,890
>500 ft	22,236	16,190	27,454

The distribution of rat reports varies by location type and proximity to parks. Residential locations consistently account for the largest number of reports across all distance categories, especially for areas that are more than 500 feet from parks. Commercial and mixed-use buildings and other non-park locations contribute to a smaller, yet notable proportion of reports, whereas open and public spaces have the fewest. This suggests that rat complaints are concentrated in places where people live or work, rather than in parks themselves. This highlights the influence of human presence and reporting behavior on observed rat activity.

#Calculate proportions
location_summary_prop <- location_summary %>%
  group_by(distance_category) %>%
  mutate(prop_reports = n_reports / sum(n_reports)) %>%
  ungroup()

# Plot
ggplot(location_summary_prop, aes(x = distance_category, y = prop_reports, fill = location_group_clean)) +
  geom_col(color = "white") +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_fill_manual(values = custom_colors) +
  labs(
    x = "Distance to Nearest Park",
    y = "Proportion of Rat Reports",
    fill = "Location Type",
    title = "Proportional Rat Reports by Location Type and Distance to Parks (2018–2025)"
  ) +
  theme_minimal(base_size = 10) +
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5),
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.title = element_text(face = "bold"),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank()
  )

V. Conclusion

This analysis suggests that while rat sightings are frequently reported near parks, neighborhoods with greater park acreage or in closer proximity to parks consistently generate fewer 311 rodent complaints than areas that are farther away. The bootstrap results reinforce this pattern, showing a clear increase in reported rat activity as distance from parks increases.

Rather than serving as primary drivers of rat complaints, parks appear to influence where are seen, but not where most reports come from. Parks provide open, visible space that may make occasional rat sightings more noticeable, however they are not the primary drivers of overall reporting volume. Instead, rat complaints are concentrated primarily in residential and mixed-use areas, where housing density and human activity are highest.

In the context of the overarching question, these results highlight that 311 rat complaints are shaped by human presence and reporting behavior, not just by the underlying distribution of rats. Reported patters are shaped mainly by where people live and work, rather by by ecological factors alone.

While parks may increase local visibility of rats, citywide rat complaint patters are primarily driven by urban density and residential and commercial land use. Parks alone do not directly generate higher levels of reported rat activity in New York City.

References

---

Footnotes

1. NYC Parks. (2025). NYC Parks. Retrieved from https://www.nycgovparks.org/about/faq