Do rat complaints differ systematically between NYC and Chicago

Introduction

In the United States, many cities use centralized non-emergency reporting systems like 311, through which residents can report issues including rodent activity. These systems are often viewed as indicators of urban conditions, the complaint volumes are shaped by socioeconomic factors, trust in public institutions, and ease of reporting. Neighborhoods with similar environmental conditions may have different levels of reported complaints, depending on institutional design and resident engagement. Understanding how reporting systems produce observed outcomes is essential for interpreting administrative data and for evaluating the effectiveness of urban policy interventions. This paper examines reported rat activity in two of the largest and densest U.S. cities: New York City and Chicago. Both cities face well established rodent challenges, and maintain extensive 311 reporting infrastructures.The enforcement mechanisms, and public usage of these systems differ in meaningful ways. New York City’s 311 platform emphasizes resident initiated reporting and public transparency, whereas Chicago’s system emphasizes on service requests tied to scheduled baiting and inspection activity. These institutional differences raise an important empirical question: do observed differences in rat complaints reflect true differences in rat prevalence, or differences in how residents and governments engage with reporting systems? To address this question, the study constructs a monthly panel dataset of rat related 311 complaints at the ZIP Code Tabulation Area (ZCTA) level from 2010 through 2025. The analysis combines complaint data with weather conditions, local unemployment rates, and neighborhood characteristics drawn from the American Census. To account for systematic heterogeneity across neighborhoods, ZCTAs are grouped into socioeconomically similar clusters using principal component analysis and k-means clustering. This approach allows comparisons between New York City and Chicago within broadly comparable neighborhood types, reducing confounding from structural differences in income, population density, and racial composition.

Data

The primary outcome variable in this study is the monthly count of rat-related service requests recorded through 311 systems in New York City and Chicago. For both cities, individual service requests are timestamped and geolocated using latitude and longitude coordinates. Records with missing or invalid spatial information are removed. Each observation is assigned to a ZIP Code Tabulation Area (ZCTA) using spatial joins with U.S. Census Bureau ZCTA shapefiles. Service requests are then aggregated to the ZCTA–month level, producing a balanced structure at the neighborhood scale. The resulting dependent variable, rat_count, measures the total number of rat-related 311 requests reported within each ZCTA in a given month. Although both datasets capture rodent-related complaints, it is important to note that they reflect slightly different institutional processes. New York City’s data emphasize resident sightings, while Chicago’s data emphasize service requests tied to baiting and inspection activity. These differences are central to the interpretation of cross-city comparisons in reported outcomes. Weather conditions are used to account for biological and behavioral activity of rats. Daily temperature data are obtained from the National Oceanic and Atmospheric Administration for a central monitoring station in each city. For New York City, observations come from station Central Park, and for Chicago from station Midway Airport. Daily maximum and minimum temperatures are converted Celsius to Fahrenheit and aggregated to the monthly level. For each city and month, average maximum temperature and average minimum temperature are computed. Unemployment rates are included to account for potential correlations between economic stress, time spent in residential environments, and the likelihood to submit 311 complaints. Because unemployment rates vary only at the city month level, they capture broad macroeconomic conditions. Local economic conditions are measured using monthly unemployment rates obtained from the Federal Reserve Economic Data database.

#importing 
library(DT)
library(dplyr)
library(ggplot2)
library(forecast)
library(tseries)
library(lubridate)
library(readr)
library(rnoaa)
library(zoo)
library(tidyr)
library(tidygeocoder)
library(purrr)  
library(sf)
#temp
temp_data_nyc <- meteo_pull_monitors(
  monitors = "USW00094728",
  var = c("TMAX","TMIN"),
  date_min = "2010-01-01",
  date_max = "2025-10-01"
)
temp_data_chi <- meteo_pull_monitors(
  monitors = "USW00094849",
  var = c("TMAX","TMIN"),
  date_min = "2010-01-01",
  date_max = "2025-10-01"
)
#rat
rat_data_nyc <- read_csv("C:/Users/David/Downloads/Rat_Sightings_20251101.csv")
rat_data_chi <-read_csv("C:/Users/David/Downloads/311_Service_Requests_-_Rodent_Baiting_-_Historical_20251124.csv")

unemp_data_nyc <- read_csv("C:/Users/David/Downloads/NEWY636URN.csv")
unemp_data_chi <- read_csv("C:/Users/David/Downloads/CHIC917URN.csv")
names(unemp_data_nyc)[names(unemp_data_nyc) == "NEWY636URN"] <- "nyc_unemp"
names(unemp_data_chi)[names(unemp_data_chi) == "CHIC917URN"] <- "chi_unemp"

#import cencus data for group
library(tidycensus)
library(dplyr)
library(purrr)
library(tigris)
library(sf)  

options(tigris_use_cache = TRUE)



chi_counties <- c("031", "043")  
nyc_counties <- c("005", "047", "061", "081", "085")  



zcta_sf <- zctas(year = 2023, class = "sf")

chi_counties_sf <- counties(state = "IL", cb = TRUE, class = "sf") %>%
  filter(COUNTYFP %in% chi_counties)

nyc_counties_sf <- counties(state = "NY", cb = TRUE, class = "sf") %>%
  filter(COUNTYFP %in% nyc_counties)



chicago_zctas <- zcta_sf[st_intersects(zcta_sf, chi_counties_sf, sparse = FALSE) %>% rowSums() > 0, ]$ZCTA5CE20
nyc_zctas <- zcta_sf[st_intersects(zcta_sf, nyc_counties_sf, sparse = FALSE) %>% rowSums() > 0, ]$ZCTA5CE20



acs_vars <- c(
  pop = "B01001_001",
  income = "B19013_001",
  pov_total = "B17001_001",
  pov_below = "B17001_002",
  white = "B02001_002",
  black = "B02001_003",
  asian = "B02001_005",
  hispanic = "B03003_003"
)

 years <- 2011:2023  

acs_panel <- map_df(years, function(y) {
  message("Downloading year: ", y)
  
  get_acs(
    geography = "zcta",
    variables = acs_vars,
    year = y,
    output = "wide"
  ) %>%
    filter(GEOID %in% c(chicago_zctas, nyc_zctas)) %>%
    mutate(
      year = y,
      city = case_when(
        GEOID %in% chicago_zctas ~ "Chicago",
        GEOID %in% nyc_zctas ~ "New York"
      ),
      poverty_rate = 100 * pov_belowE / pov_totalE
    )
})
acs_panel_clean <- na.omit(acs_panel)

chicago_geoids <- unique(acs_panel_clean$GEOID[acs_panel_clean$city == "Chicago"])
nyc_geoids <- unique(acs_panel_clean$GEOID[acs_panel_clean$city == "New York"])
zcta_file <- st_read("C:/Users/David/Downloads/tl_2023_us_zcta520/tl_2023_us_zcta520.shp")

Reading layer `tl_2023_us_zcta520' from data source 
  `C:\Users\David\Downloads\tl_2023_us_zcta520\tl_2023_us_zcta520.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 33791 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -176.6967 ymin: -14.37378 xmax: 145.8305 ymax: 71.34132
Geodetic CRS:  NAD83

chicago_id <- intersect(
  chicago_geoids,
  zcta_file$ZCTA5CE20
)

length(chicago_id)

[1] 228

head(chicago_id)

[1] "46394" "60803" "60187" "60164" "60173" "60162"

zcta_chicago <- zcta_file[zcta_file$ZCTA5CE20 %in% chicago_id, ]

setdiff(chicago_geoids, zcta_file$ZCTA5CE20)

character(0)

nyc_id <- intersect(
  nyc_geoids,
  zcta_file$ZCTA5CE20
)
length(nyc_id)

[1] 207

head(nyc_id)

[1] "10305" "10314" "10451" "10463" "10475" "10453"

zcta_nyc <- zcta_file[zcta_file$ZCTA5CE20 %in% nyc_id, ]

setdiff(nyc_geoids, zcta_file$ZCTA5CE20)

character(0)

#aggregate by zip and monthly
rat_data_nyc <- rat_data_nyc|>
 select(
  -`Unique Key`,
  -`Closed Date`,
  -Agency,
  -`Agency Name`,
  -`Complaint Type`,
  -Descriptor,
  -`Location Type`,
  -`Incident Address`,
  -`Street Name`,
  -`Cross Street 1`,
  -`Cross Street 2`,
  -`Intersection Street 1`,
  -`Intersection Street 2`,
  -`Address Type`,
  -City,
  -Landmark,
  -`Facility Type`,
  -Status,
  -`Due Date`,
  -`Resolution Action Updated Date`,
  -`Community Board`,
  -Borough,
  -`X Coordinate (State Plane)`,
  -`Y Coordinate (State Plane)`,
  -`Park Facility Name`,
  -`Park Borough`,
  -`Vehicle Type`,
  -`Taxi Company Borough`,
  -`Taxi Pick Up Location`,
  -`Bridge Highway Name`,
  -`Bridge Highway Direction`,
  -`Road Ramp`,
  -`Bridge Highway Segment`,
)
rat_data_nyc$Date <- as.Date(
  rat_data_nyc$`Created Date`,
  format = "%Y %b %d %I:%M:%S %p"
)
rat_data_nyc <- rat_data_nyc|>
  select(
    -`Created Date`
  )
rat_data_nyc$Month <- format(rat_data_nyc$Date, "%Y-%m")

library(sf)

rat_data_nyc <- rat_data_nyc |>
  drop_na()

rats_sf_nyc <- st_as_sf(
  rat_data_nyc,
  coords = c("Longitude", "Latitude"),
  crs = 4326,   # WGS84 (standard for lat/lon)
  remove = FALSE
)

rats_sf_nyc <- st_transform(rats_sf_nyc, st_crs(zcta_nyc))

rat_data_nyc <- st_join(
  rats_sf_nyc,
  zcta_nyc[, "ZCTA5CE20"],  # only keep the ZCTA ID column
  join = st_within
)
colnames(rat_data_nyc)[colnames(rat_data_nyc) == "ZCTA5CE20"] <- "ZCTA"

rat_data_nyc <- rat_data_nyc |>
  drop_na()
rat_data_nyc$Month <- format(rat_data_nyc$Date, "%Y-%m")

rat_data_nyc <- rat_data_nyc |>
  group_by(Month, ZCTA) |>
  summarise(rat_count = n()) |>
  arrange(Month)

any(is.na(rat_data_nyc))

[1] FALSE

#aggregate nyc temp 
temp_data_nyc$date <- as.Date(temp_data_nyc$date)
temp_data_nyc <- temp_data_nyc |>
  mutate(
    tmax_c = tmax / 10,
    tmin_c = tmin / 10,
    tmax_f = tmax_c * 9/5 + 32,
    tmin_f = tmin_c * 9/5 + 32
  )
temp_data_nyc$date <- as.Date(
  temp_data_nyc$date,
  format = "%Y %b %d %I:%M:%S %p"
)
temp_data_nyc$Month <- format(temp_data_nyc$date, "%Y-%m")

temp_data_nyc <- temp_data_nyc |>
  group_by(Month) |>
  summarise(
    avg_tmax = mean(tmax_f, na.rm = TRUE),
    avg_tmin = mean(tmin_f, na.rm = TRUE),
    .groups = "drop"
  )

any(is.na(temp_data_nyc))

[1] FALSE

unemp_data_nyc$observation_date <- as.Date(unemp_data_nyc$observation_date)

unemp_data_nyc$observation_date <- as.Date(
  unemp_data_nyc$observation_date,
  format = "%Y %b %d %I:%M:%S %p"
)
unemp_data_nyc$Month <- format(unemp_data_nyc$observation_date, "%Y-%m")

unemp_data_nyc <- unemp_data_nyc |>
  group_by(Month) |>
  summarise(
    average_unemp=mean(nyc_unemp, na.rm = TRUE),
    .groups = "drop"
  )

unemp_data_nyc <- unemp_data_nyc |>
  filter(Month > 100)
unemp_data_nyc <- unemp_data_nyc %>%
  filter(Month >= "2010-01" & Month <= "2025-10")
any(is.na(unemp_data_nyc))

[1] FALSE

#Joinin all NYC data
NYC_data <- rat_data_nyc|>
  left_join(temp_data_nyc, by = "Month") |>
  left_join(unemp_data_nyc, by = "Month")
any(is.na(NYC_data))

[1] TRUE

#creating variables (Cityvar)
NYC_data$City <- 1

#creating month dummy variables
NYC_data$MonthDate <- as.Date(paste0(NYC_data$Month, "-01"))

NYC_data$MonthNum <- as.numeric(format(NYC_data$MonthDate, "%m"))

NYC_data$MonthFactor <- factor(NYC_data$MonthNum, levels = 1:12, labels = month.abb)

month_dummies <- model.matrix(~ MonthFactor - 1, data = NYC_data)

NYC_data <- cbind(NYC_data, month_dummies)

NYC_data <- NYC_data %>%
  filter( Month <= "2025-08")


any(is.na(NYC_data))

[1] FALSE

NYC_data$ZCTA <- as.character(NYC_data$ZCTA)

rat_data_chi <- rat_data_chi %>%
    filter(!is.na(`Creation Date`))
rat_data_chi <- rat_data_chi |>
  select(
    -Status,
    -`Completion Date`,
    -`Service Request Number`,
    -`Type of Service Request`,
    -`Number of Premises Baited`,
    -`Number of Premises with Garbage`,
    -`Number of Premises with Rats`,
    -`Current Activity`,
    -`Most Recent Action`,
    -`Street Address`,
    -`X Coordinate`,
    -`Y Coordinate`,
    -Ward,
    -`Police District`,
    -`Community Area`,
  )

rat_data_chi$Date <- as.Date(
  rat_data_chi$`Creation Date`,
  format = "%m/%d/%Y"
)
rat_data_chi <- rat_data_chi|>
 select(
    -`Creation Date`
 )

rat_data_chi$Month <- format(rat_data_chi$Date, "%Y-%m")

library(sf)

rat_data_chi <- rat_data_chi |>
  drop_na()

chi_sf_nyc <- st_as_sf(
  rat_data_chi,
  coords = c("Longitude", "Latitude"),
  crs = 4326,   # WGS84 (standard for lat/lon)
  remove = FALSE
)

chi_sf_nyc <- st_transform(chi_sf_nyc, st_crs(zcta_chicago))

rat_data_chi <- st_join(
  chi_sf_nyc,
  zcta_chicago[, "ZCTA5CE20"],  # only keep the ZCTA ID column
  join = st_within
)

colnames(rat_data_chi)[colnames(rat_data_chi) == "ZCTA5CE20"] <- "ZCTA"


rat_data_chi <- rat_data_chi |>
  drop_na()

rat_data_chi <- rat_data_chi |>
  group_by(Month, ZCTA) |>
  summarise(rat_count = n()) |>
  arrange(Month)

any(is.na(rat_data_chi))

[1] FALSE

#fixing chi temp
temp_data_chi$date <- as.Date(temp_data_chi$date)
temp_data_chi <- temp_data_chi |>
  mutate(
    tmax_c = tmax / 10,
    tmin_c = tmin / 10,
    tmax_f = tmax_c * 9/5 + 32,
    tmin_f = tmin_c * 9/5 + 32
  )
temp_data_chi$date <- as.Date(
  temp_data_chi$date,
  format = "%Y %b %d %I:%M:%S %p"
)
temp_data_chi$Month <- format(temp_data_chi$date, "%Y-%m")

temp_data_chi <- temp_data_chi |>
  group_by(Month) |>
  summarise(
    avg_tmax = mean(tmax_f, na.rm = TRUE),
    avg_tmin = mean(tmin_f, na.rm = TRUE),
    .groups = "drop"
  )

any(is.na(temp_data_chi))

[1] FALSE

#fixing unemp chi

unemp_data_chi$observation_date <- as.Date(unemp_data_chi$observation_date)

unemp_data_chi$observation_date <- as.Date(
  unemp_data_chi$observation_date,
  format = "%Y %b %d %I:%M:%S %p"
)
unemp_data_chi$Month <- format(unemp_data_chi$observation_date, "%Y-%m")

unemp_data_chi <- unemp_data_chi |>
  group_by(Month) |>
  summarise(
    average_unemp=mean(chi_unemp, na.rm = TRUE),
    .groups = "drop"
  )

unemp_data_chi <- unemp_data_chi |>
  filter(Month > 100)
unemp_data_chi <- unemp_data_chi %>%
  filter(Month >= "2010-01" & Month <= "2025-10")
any(is.na(unemp_data_chi))

[1] FALSE

CHI_data <- rat_data_chi|>
  left_join(temp_data_chi, by = "Month") |>
  left_join(unemp_data_chi, by = "Month")
CHI_data <- CHI_data %>%
  filter( Month >= "2010-01" & Month <= "2025-08" )

any(is.na(CHI_data))

[1] FALSE

#city var and month dummy
CHI_data$City <- 0

CHI_data$MonthDate <- as.Date(paste0(CHI_data$Month, "-01"))

CHI_data$MonthNum <- as.numeric(format(CHI_data$MonthDate, "%m"))

CHI_data$MonthFactor <- factor(CHI_data$MonthNum, levels = 1:12, labels = month.abb)

month_dummies <- model.matrix(~ MonthFactor - 1, data = CHI_data)

CHI_data <- cbind(CHI_data, month_dummies)

#appending data

combined_data <- rbind(CHI_data, NYC_data)
unique_zcta <-unique(combined_data$ZCTA)
acs_panel_clean$GEOID <- as.character(acs_panel_clean$GEOID)
combined_data$ZCTA <- as.character(combined_data$ZCTA)
df_filtered <- acs_panel_clean[acs_panel_clean$GEOID %in% combined_data$ZCTA, ]

Grouping

To take into account persistent structural differences across neighborhoods, the analysis uses demographic and socioeconomic data from the American Community Survey. ZCTA-level data includes population size, median household income, poverty rates, and racial and ethnic composition. To reduce dimensionality and avoid multicollinearity, principal component analysis is applied to standardized measures of income, poverty, population size, and racial composition. The first three principal components, which capture the majority of variance in neighborhood characteristics, are used as inputs for k-means clustering. ZCTAs are then grouped into five distinct neighborhood types based on their socioeconomic profiles. This grouping strategy allows comparisons between New York City and Chicago within broadly similar neighborhood contexts, rather than relying on direct one to one ZCTA code matching. Each ZCTA is assigned a fixed group label that remains constant over time.
# final dataset
The final dataset is structured as an unbalanced panel indexed by ZCTA and month, spanning January 2010 through August 2025. Observations are kept only when rat complaint data, weather variables, and unemployment rates are jointly observed. Monthly indicator variables are constructed to control for seasonality, with January serving as the omitted reference category. The combined dataset includes observations from both cities, with a binary indicator variable identifying New York City versus Chicago.

#grouping
acs_profile <- df_filtered %>%
  filter(year <= 2017) %>%
  group_by(GEOID, city) %>%
  summarise(
    income = mean(incomeE, na.rm = TRUE),
    poverty = mean(poverty_rate, na.rm = TRUE),
    pct_white = mean(whiteE / popE, na.rm = TRUE),
    pct_black = mean(blackE / popE, na.rm = TRUE),
    pct_asian = mean(asianE / popE, na.rm = TRUE),
    pct_hispanic = mean(hispanicE / popE, na.rm = TRUE),
    log_pop = mean(log(popE), na.rm = TRUE),
    .groups = "drop"
  )

X <- acs_profile %>%
  select(income, poverty, pct_white, pct_black, pct_asian, pct_hispanic, log_pop)
X_scaled <- scale(X)

pca <- prcomp(X_scaled)
summary(pca)

Importance of components:
                          PC1    PC2    PC3    PC4     PC5    PC6     PC7
Standard deviation     1.8507 1.1334 0.9917 0.8373 0.60320 0.4747 0.13005
Proportion of Variance 0.4893 0.1835 0.1405 0.1002 0.05198 0.0322 0.00242
Cumulative Proportion  0.4893 0.6728 0.8133 0.9134 0.96539 0.9976 1.00000

PCs <- pca$x[, 1:3]

set.seed(123)
k <- 5

km <- kmeans(PCs, centers = k, nstart = 50)
acs_profile$geo_group <- km$cluster

aggregate(
  acs_profile[, c("income", "poverty")],
  by = list(acs_profile$geo_group),
  mean
)

  Group.1    income   poverty
1       1  36250.79 30.730686
2       2  61328.96 15.610097
3       3  63384.48 14.993007
4       4  44146.95 24.570097
5       5 108822.59  7.561256

table(acs_profile$geo_group, acs_profile$city)

   
    Chicago New York
  1       9       34
  2       3       37
  3      35       46
  4      15       27
  5      17       37

Interpretation of groups:

The clustering identifies five distinct neighborhood types that differ systematically in income, poverty, and racial composition. These groups represent meaningful socioeconomic and demographic class.
## Group 1

Group 1 is the most economically disadvantaged cluster, with an average household income of approximately $36,000 and a poverty rate exceeding 30 percent. Racial composition is relatively mixed, with no single majority group. The relatively high standard deviations across racial shares indicate moderate heterogeneity within this group.
## Group 2

Group 2 has moderate incomes of about $61,000 and substantially lower poverty rates. This group is notable for its racial diversity, with a substantial Asian population , alongside White and Hispanic residents. The dispersion measures suggest more internal heterogeneity than other groups.
## Group 3

Group 3 exhibits relatively high average income , around $63,000 and low poverty rates . These neighborhoods are predominantly White, with smaller Black, Asian, and Hispanic populations. Compared to Group 2, Group 3 is more racially homogeneous and economically stable.
## Group 4

Group 4 is economically disadvantaged, with average income around $44,000 and poverty rates above 24 percent. These neighborhoods are overwhelmingly Black , with relatively little racial heterogeneity. Low standard deviations across racial shares indicate strong internal homogeneity.
## Group 5 Group 5 is the most affluent cluster, with average income exceeding $108,000 and poverty rates below 8 percent. These neighborhoods are predominantly White, with small Asian and Hispanic minorities and very little internal heterogeneity. This group is economically and demographically distinct from all others.

race_by_group <- acs_profile %>%
  group_by(geo_group) %>%
  summarise(
    white = mean(pct_white, na.rm = TRUE),
    black = mean(pct_black, na.rm = TRUE),
    asian = mean(pct_asian, na.rm = TRUE),
    hispanic = mean(pct_hispanic, na.rm = TRUE),
    n_zcta = n(),
    .groups = "drop"
  )

race_by_group

# A tibble: 5 × 6
  geo_group white  black  asian hispanic n_zcta
      <int> <dbl>  <dbl>  <dbl>    <dbl>  <int>
1         1 0.300 0.276  0.0424   0.593      43
2         2 0.413 0.0836 0.352    0.205      40
3         3 0.681 0.111  0.0906   0.247      81
4         4 0.113 0.773  0.0235   0.125      42
5         5 0.779 0.0376 0.135    0.0902     54

race_sd <- acs_profile %>%
  group_by(geo_group) %>%
  summarise(
    sd_white = sd(pct_white, na.rm = TRUE),
    sd_black = sd(pct_black, na.rm = TRUE),
    sd_asian = sd(pct_asian, na.rm = TRUE),
    sd_hispanic = sd(pct_hispanic, na.rm = TRUE)
  )

race_sd

# A tibble: 5 × 5
  geo_group sd_white sd_black sd_asian sd_hispanic
      <int>    <dbl>    <dbl>    <dbl>       <dbl>
1         1   0.129    0.135    0.0401      0.155 
2         2   0.160    0.0888   0.106       0.113 
3         3   0.124    0.109    0.0583      0.134 
4         4   0.0936   0.139    0.0208      0.0992
5         5   0.108    0.0407   0.0879      0.0417

#higher sd means <.05 homogenous
#> .2 mixed

group_summary <- acs_profile %>%
  group_by(geo_group) %>%
  summarise(
    income = mean(income, na.rm = TRUE),
    poverty = mean(poverty, na.rm = TRUE),
    pct_white = mean(pct_white, na.rm = TRUE),
    pct_black = mean(pct_black, na.rm = TRUE),
    pct_asian = mean(pct_asian, na.rm = TRUE),
    pct_hispanic = mean(pct_hispanic, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  )

group_summary

# A tibble: 5 × 8
  geo_group  income poverty pct_white pct_black pct_asian pct_hispanic     n
      <int>   <dbl>   <dbl>     <dbl>     <dbl>     <dbl>        <dbl> <int>
1         1  36251.   30.7      0.300    0.276     0.0424       0.593     43
2         2  61329.   15.6      0.413    0.0836    0.352        0.205     40
3         3  63384.   15.0      0.681    0.111     0.0906       0.247     81
4         4  44147.   24.6      0.113    0.773     0.0235       0.125     42
5         5 108823.    7.56     0.779    0.0376    0.135        0.0902    54

geo_lookup <- acs_profile %>%
  select(GEOID, geo_group)

colnames(geo_lookup)[colnames(geo_lookup) == "GEOID"] <- "ZCTA"

library(plm)   
library(lmtest) 

combined_data <- combined_data|>
  left_join(geo_lookup, by = "ZCTA")

combined_data <-combined_data|>
  drop_na()

group<-split(combined_data, combined_data$geo_group)
group1<- group[["1"]]
any(is.na(group1))

[1] FALSE

group2<- group[["2"]]
any(is.na(group2))

[1] FALSE

group3<- group[["3"]]
any(is.na(group3))

[1] FALSE

group4<- group[["4"]]
any(is.na(group4))

[1] FALSE

group5<- group[["5"]]
any(is.na(group5))

[1] FALSE

summary(group1)

    Month               ZCTA             rat_count         avg_tmax    
 Length:6974        Length:6974        Min.   :  1.00   Min.   :17.11  
 Class :character   Class :character   1st Qu.:  7.00   1st Qu.:48.57  
 Mode  :character   Mode  :character   Median : 13.00   Median :63.72  
                                       Mean   : 21.35   Mean   :62.86  
                                       3rd Qu.: 22.00   3rd Qu.:78.58  
                                       Max.   :596.00   Max.   :90.09  
                                                                       
    avg_tmin      average_unemp         City          MonthDate         
 Min.   :-3.621   Min.   : 3.200   Min.   :0.0000   Min.   :2010-01-01  
 1st Qu.:34.616   1st Qu.: 4.300   1st Qu.:1.0000   1st Qu.:2013-10-01  
 Median :47.612   Median : 5.300   Median :1.0000   Median :2017-05-01  
 Mean   :48.116   Mean   : 6.295   Mean   :0.9012   Mean   :2017-07-21  
 3rd Qu.:63.398   3rd Qu.: 8.500   3rd Qu.:1.0000   3rd Qu.:2021-05-01  
 Max.   :73.342   Max.   :17.600   Max.   :1.0000   Max.   :2025-08-01  
                                                                        
    MonthNum       MonthFactor   MonthFactorJan    MonthFactorFeb   
 Min.   : 1.000   Jun    : 598   Min.   :0.00000   Min.   :0.00000  
 1st Qu.: 3.000   Jul    : 598   1st Qu.:0.00000   1st Qu.:0.00000  
 Median : 6.000   Aug    : 597   Median :0.00000   Median :0.00000  
 Mean   : 6.422   May    : 596   Mean   :0.08446   Mean   :0.08446  
 3rd Qu.: 9.000   Mar    : 592   3rd Qu.:0.00000   3rd Qu.:0.00000  
 Max.   :12.000   Jan    : 589   Max.   :1.00000   Max.   :1.00000  
                  (Other):3404                                      
 MonthFactorMar    MonthFactorApr    MonthFactorMay    MonthFactorJun   
 Min.   :0.00000   Min.   :0.00000   Min.   :0.00000   Min.   :0.00000  
 1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.00000  
 Median :0.00000   Median :0.00000   Median :0.00000   Median :0.00000  
 Mean   :0.08489   Mean   :0.08431   Mean   :0.08546   Mean   :0.08575  
 3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.00000  
 Max.   :1.00000   Max.   :1.00000   Max.   :1.00000   Max.   :1.00000  
                                                                        
 MonthFactorJul    MonthFactorAug   MonthFactorSep    MonthFactorOct   
 Min.   :0.00000   Min.   :0.0000   Min.   :0.00000   Min.   :0.00000  
 1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.00000  
 Median :0.00000   Median :0.0000   Median :0.00000   Median :0.00000  
 Mean   :0.08575   Mean   :0.0856   Mean   :0.08087   Mean   :0.08044  
 3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:0.00000  
 Max.   :1.00000   Max.   :1.0000   Max.   :1.00000   Max.   :1.00000  
                                                                       
 MonthFactorNov    MonthFactorDec      geo_group          geometry   
 Min.   :0.00000   Min.   :0.00000   Min.   :1   MULTIPOINT   :6775  
 1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:1   POINT        : 199  
 Median :0.00000   Median :0.00000   Median :1   epsg:4269    :   0  
 Mean   :0.07872   Mean   :0.07929   Mean   :1   +proj=long...:   0  
 3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:1                       
 Max.   :1.00000   Max.   :1.00000   Max.   :1                       
                                                                     

summary(group3)

    Month               ZCTA             rat_count         avg_tmax    
 Length:10104       Length:10104       Min.   :  1.00   Min.   :17.11  
 Class :character   Class :character   1st Qu.:  3.00   1st Qu.:46.83  
 Mode  :character   Mode  :character   Median :  7.00   Median :63.72  
                                       Mean   : 23.34   Mean   :62.25  
                                       3rd Qu.: 21.00   3rd Qu.:78.37  
                                       Max.   :509.00   Max.   :90.09  
                                                                       
    avg_tmin      average_unemp         City          MonthDate         
 Min.   :-3.621   Min.   : 3.200   Min.   :0.0000   Min.   :2010-01-01  
 1st Qu.:33.697   1st Qu.: 4.400   1st Qu.:1.0000   1st Qu.:2013-08-01  
 Median :46.250   Median : 5.600   Median :1.0000   Median :2016-12-01  
 Mean   :46.613   Mean   : 6.379   Mean   :0.7669   Mean   :2017-03-29  
 3rd Qu.:62.192   3rd Qu.: 8.500   3rd Qu.:1.0000   3rd Qu.:2020-09-01  
 Max.   :73.342   Max.   :17.600   Max.   :1.0000   Max.   :2025-08-01  
                                                                        
    MonthNum       MonthFactor   MonthFactorJan    MonthFactorFeb   
 Min.   : 1.000   Aug    : 901   Min.   :0.00000   Min.   :0.00000  
 1st Qu.: 4.000   Jul    : 899   1st Qu.:0.00000   1st Qu.:0.00000  
 Median : 6.000   Jun    : 885   Median :0.00000   Median :0.00000  
 Mean   : 6.477   May    : 877   Mean   :0.07829   Mean   :0.07937  
 3rd Qu.: 9.000   Apr    : 868   3rd Qu.:0.00000   3rd Qu.:0.00000  
 Max.   :12.000   Mar    : 847   Max.   :1.00000   Max.   :1.00000  
                  (Other):4827                                      
 MonthFactorMar    MonthFactorApr    MonthFactorMay   MonthFactorJun   
 Min.   :0.00000   Min.   :0.00000   Min.   :0.0000   Min.   :0.00000  
 1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.00000  
 Median :0.00000   Median :0.00000   Median :0.0000   Median :0.00000  
 Mean   :0.08383   Mean   :0.08591   Mean   :0.0868   Mean   :0.08759  
 3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.00000  
 Max.   :1.00000   Max.   :1.00000   Max.   :1.0000   Max.   :1.00000  
                                                                       
 MonthFactorJul    MonthFactorAug    MonthFactorSep    MonthFactorOct   
 Min.   :0.00000   Min.   :0.00000   Min.   :0.00000   Min.   :0.00000  
 1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.00000  
 Median :0.00000   Median :0.00000   Median :0.00000   Median :0.00000  
 Mean   :0.08897   Mean   :0.08917   Mean   :0.08343   Mean   :0.08205  
 3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.00000  
 Max.   :1.00000   Max.   :1.00000   Max.   :1.00000   Max.   :1.00000  
                                                                        
 MonthFactorNov    MonthFactorDec     geo_group          geometry   
 Min.   :0.00000   Min.   :0.0000   Min.   :3   MULTIPOINT   :8917  
 1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:3   POINT        :1187  
 Median :0.00000   Median :0.0000   Median :3   epsg:4269    :   0  
 Mean   :0.07779   Mean   :0.0768   Mean   :3   +proj=long...:   0  
 3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:3                       
 Max.   :1.00000   Max.   :1.0000   Max.   :3                       
                                                                    

summary(group5)

    Month               ZCTA             rat_count          avg_tmax    
 Length:5393        Length:5393        Min.   :  1.000   Min.   :17.11  
 Class :character   Class :character   1st Qu.:  2.000   1st Qu.:48.72  
 Mode  :character   Mode  :character   Median :  4.000   Median :65.04  
                                       Mean   :  7.095   Mean   :63.26  
                                       3rd Qu.:  9.000   3rd Qu.:78.93  
                                       Max.   :158.000   Max.   :90.09  
                                                                        
    avg_tmin      average_unemp         City          MonthDate         
 Min.   :-3.621   Min.   : 3.200   Min.   :0.0000   Min.   :2010-01-01  
 1st Qu.:34.616   1st Qu.: 4.400   1st Qu.:1.0000   1st Qu.:2013-11-01  
 Median :48.776   Median : 5.400   Median :1.0000   Median :2017-03-01  
 Mean   :47.794   Mean   : 6.283   Mean   :0.7996   Mean   :2017-06-22  
 3rd Qu.:63.104   3rd Qu.: 8.400   3rd Qu.:1.0000   3rd Qu.:2021-02-01  
 Max.   :73.342   Max.   :17.600   Max.   :1.0000   Max.   :2025-08-01  
                                                                        
    MonthNum       MonthFactor   MonthFactorJan    MonthFactorFeb   
 Min.   : 1.000   Aug    : 511   Min.   :0.00000   Min.   :0.00000  
 1st Qu.: 4.000   Jul    : 497   1st Qu.:0.00000   1st Qu.:0.00000  
 Median : 7.000   Jun    : 492   Median :0.00000   Median :0.00000  
 Mean   : 6.554   Oct    : 484   Mean   :0.07528   Mean   :0.07306  
 3rd Qu.: 9.000   May    : 472   3rd Qu.:0.00000   3rd Qu.:0.00000  
 Max.   :12.000   Sep    : 463   Max.   :1.00000   Max.   :1.00000  
                  (Other):2474                                      
 MonthFactorMar    MonthFactorApr    MonthFactorMay    MonthFactorJun   
 Min.   :0.00000   Min.   :0.00000   Min.   :0.00000   Min.   :0.00000  
 1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.00000  
 Median :0.00000   Median :0.00000   Median :0.00000   Median :0.00000  
 Mean   :0.07881   Mean   :0.08289   Mean   :0.08752   Mean   :0.09123  
 3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.00000  
 Max.   :1.00000   Max.   :1.00000   Max.   :1.00000   Max.   :1.00000  
                                                                        
 MonthFactorJul    MonthFactorAug    MonthFactorSep    MonthFactorOct   
 Min.   :0.00000   Min.   :0.00000   Min.   :0.00000   Min.   :0.00000  
 1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.00000  
 Median :0.00000   Median :0.00000   Median :0.00000   Median :0.00000  
 Mean   :0.09216   Mean   :0.09475   Mean   :0.08585   Mean   :0.08975  
 3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.00000  
 Max.   :1.00000   Max.   :1.00000   Max.   :1.00000   Max.   :1.00000  
                                                                        
 MonthFactorNov    MonthFactorDec      geo_group          geometry   
 Min.   :0.00000   Min.   :0.00000   Min.   :5   MULTIPOINT   :4109  
 1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:5   POINT        :1284  
 Median :0.00000   Median :0.00000   Median :5   epsg:4269    :   0  
 Mean   :0.07806   Mean   :0.07065   Mean   :5   +proj=long...:   0  
 3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:5                       
 Max.   :1.00000   Max.   :1.00000   Max.   :5                       
                                                                     

Empirical Strategy and Model Specification

The objective of this analysis is to assess whether differences in reporting systems contribute to systematic differences in reported rat sightings between New York City and Chicago, conditional on comparable neighborhood characteristics. Because raw comparisons across cities risk confounding institutional differences with underlying socioeconomic and demographic variation, the empirical strategy proceeds in two steps: grouping ZCTA codes into comparable neighborhood types, and then estimating regression models within these matched groups.

panel_data1 <- pdata.frame(group1, index = c("ZCTA", "MonthDate"))
panel_data2 <- pdata.frame(group2, index = c("ZCTA", "MonthDate"))
panel_data3 <- pdata.frame(group3, index = c("ZCTA", "MonthDate"))
panel_data4 <- pdata.frame(group4, index = c("ZCTA", "MonthDate"))
panel_data5 <- pdata.frame(group5, index = c("ZCTA", "MonthDate"))

formula <- rat_count ~ avg_tmax+ avg_tmin + City + average_unemp+
  MonthFactorFeb + MonthFactorMar + MonthFactorApr + MonthFactorMay +
  MonthFactorJun + MonthFactorJul + MonthFactorAug + MonthFactorSep +
  MonthFactorOct + MonthFactorNov + MonthFactorDec

re_model1 <- plm(formula, data = panel_data1, model = "random")
summary(re_model1)

Oneway (individual) effect Random Effect Model 
   (Swamy-Arora's transformation)

Call:
plm(formula = formula, data = panel_data1, model = "random")

Unbalanced Panel: n = 43, T = 1-188, N = 6974

Effects:
                  var std.dev share
idiosyncratic 329.133  18.142 0.822
individual     71.283   8.443 0.178
theta:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.09337 0.84477 0.84517 0.83797 0.84517 0.84517 

Residuals:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-129.60   -6.87   -1.04    0.14    5.42  479.14 

Coefficients:
                 Estimate Std. Error  z-value  Pr(>|z|)    
(Intercept)     25.989595   4.199234   6.1891 6.050e-10 ***
avg_tmax         2.482410   0.131241  18.9148 < 2.2e-16 ***
avg_tmin        -1.567115   0.145347 -10.7819 < 2.2e-16 ***
City           -64.750561   3.717750 -17.4166 < 2.2e-16 ***
average_unemp   -1.274877   0.090625 -14.0677 < 2.2e-16 ***
MonthFactorFeb  -6.193935   1.091036  -5.6771 1.370e-08 ***
MonthFactorMar -13.218440   1.303761 -10.1387 < 2.2e-16 ***
MonthFactorApr -23.618210   1.798755 -13.1303 < 2.2e-16 ***
MonthFactorMay -29.299073   2.369403 -12.3656 < 2.2e-16 ***
MonthFactorJun -33.648822   2.908664 -11.5685 < 2.2e-16 ***
MonthFactorJul -34.153379   3.330873 -10.2536 < 2.2e-16 ***
MonthFactorAug -30.965669   3.173978  -9.7561 < 2.2e-16 ***
MonthFactorSep -28.105409   2.760686 -10.1806 < 2.2e-16 ***
MonthFactorOct -17.187411   2.085653  -8.2408 < 2.2e-16 ***
MonthFactorNov -14.222230   1.465429  -9.7052 < 2.2e-16 ***
MonthFactorDec  -5.897597   1.196574  -4.9287 8.276e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    2947200
Residual Sum of Squares: 2329700
R-Squared:      0.2096
Adj. R-Squared: 0.20789
Chisq: 1770.82 on 15 DF, p-value: < 2.22e-16

re_model2 <- plm(formula, data = panel_data2, model = "random")
summary(re_model2)

Oneway (individual) effect Random Effect Model 
   (Swamy-Arora's transformation)

Call:
plm(formula = formula, data = panel_data2, model = "random")

Unbalanced Panel: n = 40, T = 1-188, N = 5514

Effects:
                 var std.dev share
idiosyncratic 56.565   7.521 0.546
individual    47.054   6.860 0.454
theta:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2612  0.9071  0.9157  0.9097  0.9188  0.9203 

Residuals:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-64.972  -2.744  -0.542   0.070   1.944 129.786 

Coefficients:
                 Estimate Std. Error z-value  Pr(>|z|)    
(Intercept)     32.051792   4.625902  6.9288 4.245e-12 ***
avg_tmax         0.723968   0.069508 10.4156 < 2.2e-16 ***
avg_tmin        -0.395263   0.076434 -5.1713 2.325e-07 ***
City           -45.392471   4.609646 -9.8473 < 2.2e-16 ***
average_unemp   -0.385688   0.041565 -9.2791 < 2.2e-16 ***
MonthFactorFeb  -2.064532   0.540745 -3.8179 0.0001346 ***
MonthFactorMar  -3.753937   0.643158 -5.8367 5.324e-09 ***
MonthFactorApr  -7.270115   0.891171 -8.1579 3.408e-16 ***
MonthFactorMay  -9.128410   1.161796 -7.8572 3.930e-15 ***
MonthFactorJun -10.487175   1.435516 -7.3055 2.762e-13 ***
MonthFactorJul -11.404693   1.650212 -6.9110 4.811e-12 ***
MonthFactorAug -10.556788   1.568262 -6.7315 1.679e-11 ***
MonthFactorSep  -9.272428   1.360035 -6.8178 9.245e-12 ***
MonthFactorOct  -6.368148   1.020717 -6.2389 4.407e-10 ***
MonthFactorNov  -4.793265   0.714603 -6.7076 1.979e-11 ***
MonthFactorDec  -1.818583   0.586242 -3.1021 0.0019215 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    352760
Residual Sum of Squares: 312390
R-Squared:      0.11453
Adj. R-Squared: 0.11211
Chisq: 708.496 on 15 DF, p-value: < 2.22e-16

re_model3 <- plm(formula, data = panel_data3, model = "random")
summary(re_model3)

Oneway (individual) effect Random Effect Model 
   (Swamy-Arora's transformation)

Call:
plm(formula = formula, data = panel_data3, model = "random")

Unbalanced Panel: n = 81, T = 1-188, N = 10104

Effects:
                 var std.dev share
idiosyncratic 451.75   21.25 0.444
individual    566.59   23.80 0.556
theta:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3340  0.9228  0.9330  0.9249  0.9341  0.9350 

Residuals:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-174.78   -8.11   -1.24    0.31    6.64  332.26 

Coefficients:
                 Estimate Std. Error  z-value  Pr(>|z|)    
(Intercept)    -10.116629   4.792704  -2.1108   0.03479 *  
avg_tmax         2.674180   0.109888  24.3355 < 2.2e-16 ***
avg_tmin        -1.618165   0.124455 -13.0020 < 2.2e-16 ***
City           -40.909012   5.594621  -7.3122 2.628e-13 ***
average_unemp   -1.580751   0.089792 -17.6045 < 2.2e-16 ***
MonthFactorFeb  -6.365375   1.090933  -5.8348 5.386e-09 ***
MonthFactorMar -12.314881   1.246568  -9.8790 < 2.2e-16 ***
MonthFactorApr -23.174395   1.666402 -13.9068 < 2.2e-16 ***
MonthFactorMay -31.899873   2.224432 -14.3407 < 2.2e-16 ***
MonthFactorJun -34.983240   2.713585 -12.8919 < 2.2e-16 ***
MonthFactorJul -36.880068   3.093890 -11.9203 < 2.2e-16 ***
MonthFactorAug -31.844389   2.960124 -10.7578 < 2.2e-16 ***
MonthFactorSep -30.503959   2.576840 -11.8377 < 2.2e-16 ***
MonthFactorOct -18.242032   1.966770  -9.2751 < 2.2e-16 ***
MonthFactorNov -13.745140   1.419593  -9.6825 < 2.2e-16 ***
MonthFactorDec  -6.570178   1.187176  -5.5343 3.125e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    5823300
Residual Sum of Squares: 4599200
R-Squared:      0.21039
Adj. R-Squared: 0.20921
Chisq: 2672.85 on 15 DF, p-value: < 2.22e-16

re_model4 <- plm(formula, data = panel_data4, model = "random")
summary(re_model4)

Oneway (individual) effect Random Effect Model 
   (Swamy-Arora's transformation)

Call:
plm(formula = formula, data = panel_data4, model = "random")

Unbalanced Panel: n = 42, T = 52-188, N = 5823

Effects:
                 var std.dev share
idiosyncratic 254.99   15.97 0.605
individual    166.51   12.90 0.395
theta:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.8309  0.8772  0.9082  0.8971  0.9101  0.9101 

Residuals:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-97.564  -6.696  -1.283   0.087   5.016 246.894 

Coefficients:
                 Estimate Std. Error  z-value  Pr(>|z|)    
(Intercept)     10.217518   4.102815   2.4904   0.01276 *  
avg_tmax         1.603357   0.107420  14.9261 < 2.2e-16 ***
avg_tmin        -0.745936   0.121253  -6.1519 7.657e-10 ***
City           -45.173061   4.353607 -10.3760 < 2.2e-16 ***
average_unemp   -0.932429   0.088432 -10.5440 < 2.2e-16 ***
MonthFactorFeb  -4.632376   1.068636  -4.3348 1.459e-05 ***
MonthFactorMar  -6.684170   1.230199  -5.4334 5.529e-08 ***
MonthFactorApr -14.841923   1.646590  -9.0137 < 2.2e-16 ***
MonthFactorMay -20.967027   2.200369  -9.5289 < 2.2e-16 ***
MonthFactorJun -26.562976   2.679166  -9.9146 < 2.2e-16 ***
MonthFactorJul -29.309486   3.054775  -9.5946 < 2.2e-16 ***
MonthFactorAug -23.896960   2.922176  -8.1778 2.891e-16 ***
MonthFactorSep -21.455002   2.543848  -8.4341 < 2.2e-16 ***
MonthFactorOct -11.992444   1.938320  -6.1870 6.131e-10 ***
MonthFactorNov  -8.941400   1.397416  -6.3985 1.569e-10 ***
MonthFactorDec  -5.503209   1.178385  -4.6701 3.010e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    1894400
Residual Sum of Squares: 1492300
R-Squared:      0.21229
Adj. R-Squared: 0.21026
Chisq: 1542.21 on 15 DF, p-value: < 2.22e-16

re_model5 <- plm(formula, data = panel_data5, model = "random")
summary(re_model5)

Oneway (individual) effect Random Effect Model 
   (Swamy-Arora's transformation)

Call:
plm(formula = formula, data = panel_data5, model = "random")

Unbalanced Panel: n = 54, T = 1-188, N = 5393

Effects:
                 var std.dev share
idiosyncratic 44.484   6.670 0.685
individual    20.491   4.527 0.315
theta:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1726  0.8490  0.8857  0.8638  0.8920  0.8932 

Residuals:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-19.754  -3.057  -0.740   0.154   2.052 142.288 

Coefficients:
                Estimate Std. Error z-value  Pr(>|z|)    
(Intercept)    -3.187344   1.575338 -2.0233 0.0430448 *  
avg_tmax        0.417800   0.048628  8.5917 < 2.2e-16 ***
avg_tmin       -0.195785   0.055265 -3.5427 0.0003961 ***
City           -3.688371   1.478417 -2.4948 0.0126025 *  
average_unemp  -0.332520   0.038618 -8.6105 < 2.2e-16 ***
MonthFactorFeb -1.301654   0.482163 -2.6996 0.0069420 ** 
MonthFactorMar -2.021260   0.547994 -3.6885 0.0002256 ***
MonthFactorApr -4.021632   0.737277 -5.4547 4.905e-08 ***
MonthFactorMay -5.014695   0.984145 -5.0955 3.479e-07 ***
MonthFactorJun -5.841519   1.206797 -4.8405 1.295e-06 ***
MonthFactorJul -5.814506   1.378971 -4.2166 2.481e-05 ***
MonthFactorAug -5.607623   1.317342 -4.2568 2.074e-05 ***
MonthFactorSep -4.866720   1.146522 -4.2448 2.188e-05 ***
MonthFactorOct -3.224771   0.866030 -3.7236 0.0001964 ***
MonthFactorNov -2.787656   0.617182 -4.5167 6.280e-06 ***
MonthFactorDec -1.828190   0.521382 -3.5064 0.0004542 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    272630
Residual Sum of Squares: 240820
R-Squared:      0.1172
Adj. R-Squared: 0.11474
Chisq: 726.353 on 15 DF, p-value: < 2.22e-16

Panel Data Structure

he resulting dataset is a ZCTA month panel spanning 2010–2025. The primary outcome variable is the number of reported rat sightings in a given ZCTA code and month. Each observation is associated with a city indicator, neighborhood group membership, weather controls (average monthly minimum and maximum temperatures), local unemployment rates, and month fixed effects to account for seasonality

\[ \{rat_count}_{it}^{(g)} = \beta_0^{(g)} + \beta_1^{(g)} \, \{avg_tmax}_{it} + \beta_2^{(g)} \, \{avg_tmin}_i+ \beta_3^{(g)}\ \{City}_{0,1} + \sum_{m=\text{Feb}}^{\{Dec}} \gamma_m^{(g)} \, \{MonthFactor}_{mt} + \epsilon_{it}^{(g)} \]

- RatCountz,t​ is the number of rat sightings reported in ZCTA code zzz during month ttt.
- _{0,1}​ is an indicator variable equal to 0 if ZCTA code zzz is located in Chicago and 1 if it is located in New York City, which serves as the reference category.

• beta_1^{(g)} , text{avg_tmax}_{it} + _2^{(g)} , text{avg_tmin}_i​ is a vector of time-varying controls, including average monthly maximum and minimum temperatures and the local unemployment rate, which account for environmental and economic conditions that may influence rat activity and reporting behavior.
  

• _m^{(g)} denotes month fixed effects, which capture common seasonal patterns in rat sightings across cities, such as higher reports during warmer months.
  

εz,t_{z,t}εz,t​ is the error term, capturing unobserved factors that vary across ZCTA codes and over time.

To further enhance comparability, the model is estimated separately within each neighborhood group. This approach allows the city coefficient to be interpreted as the difference in reported rat sightings between Chicago and New York City within socioeconomically similar ZCTA, rather than across the full and highly heterogeneous landscape.This group strategy avoids reliance on strong assumptions required by random effects models and mitigates concerns that unobserved ZCTA characteristics such as housing density, sanitation infrastructure, or historical pest prevalence are correlated with city membership.

Results

## 5.1 Temperature Effects

Across all neighborhood groups, temperature variables are statistically significant and display consistent directional effects. Higher average maximum temperature is positively associated with rat sightings, while higher average minimum temperature is negatively associated. However, the magnitude of these effects is moderate. A 10°F increase in average maximum temperature is associated with approximately 4 to 27 additional monthly rat sightings, depending on neighborhood group, while a comparable increase in minimum temperature reduces reported sightings by roughly 2 to 16. These results suggest that temperature influences rat activity and reporting behavior, but weather alone explains only a limited share of overall variation in rat complaints relative to seasonal and institutional factors.

## 5.2 Seasonal Patterns
Month fixed effects are uniformly negative relative to January, which serves as the reference month. This indicates that, conditional on weather and economic conditions, reported rat sightings are highest in January across all neighborhood groups. This pattern likely reflects increased indoor rat activity and heightened reporting during winter months, rather than biological seasonality alone. Because temperature is explicitly controlled for, the month effects capture non-climatic seasonal factors such as reporting behavior, building usage patterns, and enforcement cycles.

## 5.3 Economic Conditions
Local unemployment rates are negatively associated with rat sightings in all specifications. A one-percentage-point increase in unemployment is associated with approximately 0.3 to 1.6 fewer monthly rat reports, depending on neighborhood group. While this relationship is not the primary focus of the analysis, it suggests that economic conditions may influence reporting behavior or service utilization rather than underlying rat prevalence.
## 5.4 City-Level Differences

The estimated coefficient on the city indicator is large, statistically significant, and robust across all model specifications and neighborhood groupings. Because the city variable is coded as 1 for New York City and 0 for Chicago, the negative coefficient implies that, holding constant weather conditions, local economic factors, and seasonal effects, ZCTA codes in New York City consistently report fewer rat sightings than observationally similar ZCTA codes in Chicago.  
<b/>
The magnitude of this difference is substantial. Depending on the neighborhood group, New York City ZCTA codes report between approximately 20 and 50 fewer rat sightings per month relative to Chicago ZCTA. 

The diversity of the city effect across neighborhood groups further supports this interpretation. The estimated New York City reporting gap is largest in lower-income and higher-density neighborhoods, where residents may be more likely to engage with a simple, targeted reporting system. In contrast, the gap is smaller in higher-income ZCTA codes, where residents may rely on private pest control services or alternative reporting channels that are not captured in public 311 data. This pattern is consistent with a mechanism in which institutional design interacts with neighborhood characteristics to shape reporting behavior.

Additionally, the persistence of the city effect after controlling for month fixed effects is particularly important. Month fixed effects absorb broad seasonal patterns in reporting behavior that are common across cities. The fact that a sizable city coefficient remains indicates that the observed differences are not driven by differences in seasonal timing, weather seasonality, or cyclical public awareness campaigns.

Taken together, these findings suggest that the observed disparity in rat sightings between Chicago and New York City is unlikely to reflect differences in environmental conditions. Instead, the results are consistent with the hypothesis that Chicago’s 311 reporting system generates higher reported rat activity through lower reporting barriers and more standardized complaint categorization. While the analysis does not establish causality, it provides strong descriptive evidence that institutional reporting structures play a central role in shaping observed urban pest complaint data.  
<b/>

5.5 Model Fit and Variation

The models explain between 11% and 21% of total variation in rat sightings, depending on neighborhood group. Decomposition of variance components indicates that both within ZCTA and ZCTA variation contribute meaningfully to observed outcomes, supporting the use of panel methods. Importantly, the persistence of the city effect after conditioning on temperature and month fixed effects suggests that observed differences in reported rat sightings between Chicago and New York City are unlikely to be driven by climate or seasonal patterns alone.

library(broom)
library(dplyr)

models <- list(
  Group1 = re_model1,
  Group2 = re_model2,
  Group3 = re_model3,
  Group4 = re_model4,
  Group5 = re_model5
)

coef_df <- bind_rows(
  lapply(names(models), function(n) {
    tidy(models[[n]]) %>% mutate(Group = n)
  })
)

library(ggplot2)

coef_df %>%
  filter(term %in% c("avg_tmax", "avg_tmin", "average_unemp")) %>%
  ggplot(aes(x = term, y = estimate, color = Group)) +
  geom_point(position = position_dodge(width = 0.4), size = 3) +
  geom_errorbar(
    aes(ymin = estimate - 1.96 * std.error,
        ymax = estimate + 1.96 * std.error),
    position = position_dodge(width = 0.4),
    width = 0.2
  ) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(
    title = "Random Effects Estimates Across Groups",
    x = "Covariate",
    y = "Coefficient Estimate"
  ) +
  theme_minimal()

panel_data1$City <- factor(panel_data1$City)
panel_data2$City <- factor(panel_data2$City)
panel_data3$City <- factor(panel_data3$City)
panel_data4$City <- factor(panel_data4$City)
panel_data5$City <- factor(panel_data5$City)





panel_data1$pred <- predict(re_model1)

ggplot(panel_data1, aes(x = pred, y = rat_count, color = City)) +
  geom_point(alpha = 0.4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "black") +
  scale_color_brewer(palette = "Set2") +
  labs(
    title = "Predicted vs Actual Rat Counts (Group 1)",
    x = "Predicted",
    y = "Observed",
    color = "City"
  ) +
  theme_minimal()

panel_data2$pred <- predict(re_model2)

ggplot(panel_data2, aes(x = pred, y = rat_count, color = City)) +
  geom_point(alpha = 0.4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "black") +
  scale_color_brewer(palette = "Set2") +
  labs(
    title = "Predicted vs Actual Rat Counts (Group 2)",
    x = "Predicted",
    y = "Observed",
    color = "City"
  ) +
  theme_minimal()

panel_data3$pred <- predict(re_model3)

ggplot(panel_data3, aes(x = pred, y = rat_count, color = City)) +
  geom_point(alpha = 0.4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "black") +
  scale_color_brewer(palette = "Set2") +
  labs(
    title = "Predicted vs Actual Rat Counts (Group 3)",
    x = "Predicted",
    y = "Observed",
    color = "City"
  ) +
  theme_minimal()

panel_data4$pred <- predict(re_model4)

ggplot(panel_data4, aes(x = pred, y = rat_count, color = City)) +
  geom_point(alpha = 0.4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "black") +
  scale_color_brewer(palette = "Set2") +
  labs(
    title = "Predicted vs Actual Rat Counts (Group 4)",
    x = "Predicted",
    y = "Observed",
    color = "City"
  ) +
  theme_minimal()

panel_data5$pred <- predict(re_model5)

ggplot(panel_data5, aes(x = pred, y = rat_count, color = City)) +
  geom_point(alpha = 0.4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "black") +
  scale_color_brewer(palette = "Set2") +
  labs(
    title = "Predicted vs Actual Rat Counts (Group 5)",
    x = "Predicted",
    y = "Observed",
    color = "City"
  ) +
  theme_minimal()

ggplot(group1, aes(x = MonthDate, y = rat_count, group = ZCTA)) +
  geom_line(alpha = 0.2) +
  stat_summary(aes(group = 1), fun = mean, geom = "line",
               linewidth = 1.2, color = "black") +
  labs(
    title = "Rat Count Trends Over Time (Group 1)",
    y = "Rat Count"
  ) +
  theme_minimal()

ggplot(group2, aes(x = MonthDate, y = rat_count, group = ZCTA)) +
  geom_line(alpha = 0.2) +
  stat_summary(aes(group = 2), fun = mean, geom = "line",
               linewidth = 1.2, color = "black") +
  labs(
    title = "Rat Count Trends Over Time (Group 2)",
    y = "Rat Count"
  ) +
  theme_minimal()

ggplot(group3, aes(x = MonthDate, y = rat_count, group = ZCTA)) +
  geom_line(alpha = 0.2) +
  stat_summary(aes(group = 3), fun = mean, geom = "line",
               linewidth = 1.2, color = "black") +
  labs(
    title = "Rat Count Trends Over Time (Group 3)",
    y = "Rat Count"
  ) +
  theme_minimal()

ggplot(group4, aes(x = MonthDate, y = rat_count, group = ZCTA)) +
  geom_line(alpha = 0.2) +
  stat_summary(aes(group = 4), fun = mean, geom = "line",
               linewidth = 1.2, color = "black") +
  labs(
    title = "Rat Count Trends Over Time (Group 4)",
    y = "Rat Count"
  ) +
  theme_minimal()

ggplot(group5, aes(x = MonthDate, y = rat_count, group = ZCTA)) +
  geom_line(alpha = 0.2) +
  stat_summary(aes(group = 5), fun = mean, geom = "line",
               linewidth = 1.2, color = "black") +
  labs(
    title = "Rat Count Trends Over Time (Group 5)",
    y = "Rat Count"
  ) +
  theme_minimal()

group1$Group <- "Group 1"
group2$Group <- "Group 2"
group3$Group <- "Group 3"
group4$Group <- "Group 4"
group5$Group <- "Group 5"
library(dplyr)

all_groups <- bind_rows(group1, group2, group3, group4, group5)

ggplot(all_groups, aes(x = MonthDate, y = rat_count, group = interaction(Group, ZCTA))) +
  geom_line(alpha = 0.15, aes(color = Group)) +
  stat_summary(
    aes(group = Group, color = Group),
    fun = mean,
    geom = "line",
    linewidth = 1.4
  ) +
  scale_color_brewer(palette = "Set1") +
  labs(
    title = "Rat Count Trends Over Time by Group",
    y = "Rat Count",
    color = "Group"
  ) +
  theme_minimal()

conclusion

This analysis shows clear differences in reported rat sightings between Chicago and New York City, even after accounting for weather, unemployment, and seasonal effects. Chicago ZCTA codes consistently show higher monthly counts, suggesting that institutional and reporting mechanisms, rather than purely environmental factors, shape observed complaint patterns. Neighborhood level heterogeneity further highlights how local demographics interact with reporting systems. These findings indicate the importance of considering institutional context when interpreting urban pest data and suggest that city-level differences in reporting infrastructure can significantly influence the patterns captured in public complaint datasets.