Patterns of NYC Squirrel Distribution by Nearest Pizza Restaurant
Christopher J Costa
2024-06-03
Background
The New York City Squirrel Census was a citizen science endeavor sponsored by Ever True. It focused on the distribution of the Eastern Grey Squirrel, Sciurus carolinesis, in parks in New York. The data has been made freely available on github and their website. There have been two censuses, the first in 2018 focused on Central Park and cataloged approximately 3400 squirrels. There was a second, smaller census conducted throughout NYC in 2020, identifying 340 squirrels. Attributes for each squirrel were recorded by observers and included: location, fur color, activity, sighting at ground level or above ground. While the squirrel data is interesting on its own, given that NYC wildlife is significantly dependent on human city dwellers, I decided to explore this relationship with a NYC staple in mind: pizza.
The pizza dataset, “Price for a Slice” was developed by Andy Kriebel at the Information Data Lab School. This list compiles NYC pizza restaurants by name, location, price per slice, and, if applicable, Instagram page. This data was made freely available on data.world.
Together, this data presents a unique story of NYC squirrels and their interactions with human city dwellers and how the events of the 2020 COVID-19 pandemic affected this relationship. I sought to examine the distribution of these squirrels among restaurants with a focus on factors that may predict location including: price of a slice, fur color, behavior and spotted location.
Methods
Data
All data used for this analysis has been made publicly available by the owner. See below for sources.
2018 squirrel census data: Ever True, 2018 NYC Squirrel Census. Data is now available on NYC Open Data (https://data.cityofnewyork.us/Environment/2018-Central-Park-Squirrel-Census-Squirrel-Data/vfnx-vebw/about_data).
2020 squirrel census data: Ever True, 2020 NYC Open Data Week (2/28/2020-3/7/2020) (https://www.thesquirrelcensus.com/data)
Price for a slice: Collected by Andy Kriebel and hosted at data.world (https://data.world/makeovermonday/2023w5)
#reading in the 2018 squirrel data
nyc_squirrels <- readr::read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-10-29/nyc_squirrels.csv")
#reading in 2020 squirrel data
squirrel_2020<- read_csv("squirrel-data.csv")
#reading in pizza data
pizza_data<- read_excel("nyc_slice_rawdata.xlsx")Required packages
This work required the use of a number of packages for data management and statistical analysis. Additionally, ggplot2 and ggmap were used for visualization. See below for the full list of required packages.
library(dplyr)
library(plyr)
library(tidyr)
library(readr)
library(readxl)
library(ggplot2)
library(stats)
library(statmod)
library(stringr)
library(ggmap)Data Preparation
After reading in the data, both datasets required modification to make calculations more manageable. Methods were equivalent for both 2018 and 2020 data. The squirrel dataset was simplified to include only squirrel name and location. The pizza dataset required conversion of locations to numeric data.
#simplifying the squirrel dataset to only include necessary information for distance calculations
squirrel_loc<- nyc_squirrels[ , c(1,2,3)]
#simplify pizza data and convert latitude and longitude to numbers
pizza_lat_long<- pizza_data%>%
mutate(location_lat = as.numeric(location_lat, errors= 'coerce'))%>%
mutate(location_lng = as.numeric(location_lng, errors='coerce'))Distance Calculations
Code for this section was inspired in part by code discussed here (https://forum.posit.co/t/finding-the-shortest-distance-in-r/131771). Essentially, the latitude and longitude distances were used along with the Pythagorean theorem to create a list of straight line distances between each squirrel and each pizza restaurant. The minimum function was used to find the shortest distance. These were compiled into a list and associated with the pizza restaurant names. Restaurant names were then joined with the simplified squirrel dataframe.
# Calculate distance between each squirrel and all restaurants
distances <- lapply(1:nrow(squirrel_loc), function(i) {
# Get current squirrel location
squirrel_lat <- squirrel_loc$lat[i]
squirrel_long <- squirrel_loc$long[i]
# Calculate distance to each restaurant using Pythagorean Theorem
restaurant_distances <- sqrt(((pizza_lat_long$location_lat -
squirrel_lat)^2) +
((pizza_lat_long$location_lng
- squirrel_long)^2))
# Find the index of the minimum distance restaurant
min_index <- which.min(restaurant_distances)
# Return the name of the closest restaurant
pizza_lat_long$Name[min_index]
})
# Add a new column with closest pizza restaurant name
squirrel_loc$closest_pizza <- unlist(distances)
# Handle potential errors (e.g., if a squirrel is exactly at a restaurant location)
squirrel_loc$closest_pizza[is.na(squirrel_loc$closest_pizza)] <-
"Tie (equal distance to multiple restaurants)"
#squirrel_loc dataframe now contains nearest pizza restaurant and
#can be used for analysis
#renamed to differentiate from the 2020 data, process is the same
pizza_price_squirrels_2018<- squirrel_loc%>%
rename(c('closest_pizza'='Name'))
##this eliminates any potential doubles which can occur with code
pizza_price_squirrels_2018<- left_join(squirrel_loc, pizza_data, by= "Name", multiple='first')Statistical Analysis
Chi square or Fisher exact test were used for all comparisons of squirrels per restaurant, behavior at restaurant, and fur color at restaurant. Linear regression was used to test for correlations between price for a slice and squirrel number. All analysis conducted in R.
Results
Distribution of Squirrels by Restaurant
Maps showing location of squirrels spotted (yellow dots) and pizza restaurants (red triangles) are displayed below for both 2018 and 2020 data (Figure 1A & 1B, respectively). The 2018 census focused on Central Park while the 2020 census looked at NYC as a whole. Given the number of restaurants with only one squirrel spotted nearby, analysis was limited to restaurants with greater than 5 squirrels identified. For both 2018 (Figure 1C) and 2020 (Figure 1D) data, distribution of squirrels varied by restaurant. Differences in distribution were statistically significant for 2018 (𝜒2 =1431, df=12, p<2.2e-16) and for 2020 (𝜒2 =203.8, df=26, p<2.2e-16).
##conducting distribution analysis
##start by organizing datasets
pizza_price_squirrels_2018<- pizza_price_squirrels_2018%>%
rename(c("unique_squirrel_id"= 'Squirrel_ID'))%>%
full_join(., pizza_data, by="Name", multiple='first')
##create a table to list number of squirrels per restaurant&form dataframe
squirrel_distribution<- tapply(pizza_price_squirrels_2018$Squirrel_ID, pizza_price_squirrels_2018$Name,length)%>%
data.frame(Name= names(.),count=.)%>%
rename(c('count'='squirrel_count'))%>%
arrange(squirrel_count)
##test for differences in proportions of squirrels for each location
test_location<-chisq.test(squirrel_distribution$squirrel_count)###mapping squirrel and pizza restaurant locations
library(ggmap)
API_secret<- 'API_GOES_HERE'
register_google(key=API_secret)
##import NYC map
nyc_map<- get_map(location= 'central park, New York', maptype= 'hybrid',
color= 'color', source= 'google', zoom=14)
##create map and plot points
ggmap(nyc_map)+
geom_point(data=nyc_squirrels, aes(x=long, y=lat), size=0.5, color='yellow')+
geom_point(data=pizza_lat_long, aes(x=location_lng,
y=location_lat), size=1.5, color= "red",
shape='triangle')+
labs(title="Squirrels and Pizza Places-2018 Census Data")
(red triangles).
(red triangles).
Price Analysis
Analysis next focused on the price of a slice of pizza of restaurants identified as near to squirrels in 2018 and 2020. Code is shown for 2018 data and similar code was used for 2020 data. The price distribution of restaurants identified in the 2018 dataset (Figure 2A) found a mean price of $2.91, with $0.55 standard deviation. Prices identified in 2020 (Figure 2B) found a mean price of $2.59 with $1.03 standard deviation.
#descriptives on price with plot
price_summary<-pizza_price_squirrels_2018%>%
summarize(mean= mean(`Price as number`), median= median(`Price as number`),
standard_dev= sd(`Price as number`), IQR= IQR(`Price as number`))
With differences in price patterns of nearby restaurants observed, analysis focused on determining if squirrel distribution was correlated to the price of a slice. This was accomplished by organizing the data by price and tallying squirrels at each price point (code below). Pearson’s linear correlation was then used to test for correlation between price and number of squirrels. For the 2018 data (Figure 3A), there was no correlation between price and squirrel counts (R2 =0.028, p=0.28). For the 2020 data (Figure 3B), there was a significant, strong negative correlation between price and squirrel counts (R2 =0.49, p=0.003).
##create a table to list squirrel numbers as grouped by price
squirrel_price<- tapply(pizza_price_squirrels_2018$Squirrel_ID, pizza_price_squirrels_2018$`Price as number`,
length)%>%
data.frame(Price=names(.), count= .)%>%
mutate(Price =as.numeric(Price, errors='coerce'))%>%
rename(c('count'='squirrel_count'))
squirrel_price_more5<-filter(squirrel_price, squirrel_price$squirrel_count>=5)
Fur Color and Restaurant Distribution
Squirrel primary fur color identified in the censuses was either gray, black or cinnamon. In both datasets, gray squirrels were spotted most often, comprising 81% of 2018 squirrels and 88% of 2020 squirrels. Cinnamon squirrels were the second most common, 13% in 2018 and 7.6% in 2020. Black were the least commonly spotted, 3.4% in 2018 and 3.5% in 2020.
Analyzing the distribution of colors by restaurant showed differing patterns in 2018 and 2020. For the 2018 data (Figure 4A), there were no significant differences in gray (𝜒2= 91, df= 84, p= 0.28), cinnamon (𝜒2= 130, df= 120, p=0,25), or black (𝜒2= 68, df= 56, p= 0.13). However, for 2020 data (Figure 4B) there were significant differences for all fur colors: gray (𝜒2= 91.5, df= 26, p=3.2e-09), cinnamon (𝜒2= 252.3, df= 26, p< 2.2e-16), and black (𝜒2= 153, df= 26, p< 2.2e-16).
##compare coat colors at restaurants
squirrel_behavior_2018<- rename(squirrel_behavior_2018, c('primary_fur_color'='fur_color'))
#new dataframe with counts of color by restaurant with new columns with color data
counts_by_restaurant <- squirrel_behavior_2018 %>%
group_by(Name)%>%
filter(!is.na(fur_color))%>%
dplyr::summarise(sightings=sum(1),
Gray=sum(ifelse(fur_color=="Gray",1,0)),
Cinnamon=sum(ifelse(fur_color=="Cinnamon",1,0)),
Black=sum(ifelse(fur_color=="Black", 1, 0)))%>%
arrange(Name)
Patterns of Behavior by Restaurant location
Census data included behaviors that the squirrels were seen engaging in at time of identification. For the 2018 census, foraging was the most commonly seen activity (47%), followed by eating (25%), running (24%), climbing (21%) and chasing (9.2%). Similarly, foraging was also the most commonly seen activity in 2020 (27%) followed by climbing (25%), eating (20%), running (16%), chasing (9.3%), sitting (2.9%), and digging (1.5%). Notably, there were different behaviors identified in 2018 and 2020 which were excluded from analysis on comparison in behaviors by year discussed below.
Behaviors identified in 2018 did not differ in distribution by restaurant (Figure 5A). This is true for all behaviors: chasing (𝜒2= 117, df= 108, p= 0.26), climbing (𝜒2= 143, df= 132, 0.24), eating (𝜒2= 143, df= 132, p= 0.24), foraging (𝜒2= 156, df= 144, p= 0.23), and running (𝜒2= 143, df= 132, p= 0.24). However, behaviors identified in 2020 (Figure 5B) did differ by location. Notably, all behaviors differed significantly by restaurant: chasing (𝜒2= 39.9, df= 24, p= 0.02), climbing (𝜒2= 122.3, df= 24, p= 3.7e-15), digging (𝜒2= 60, df= 24, p= 6.4e-05), eating (𝜒2= 65.5, df= 24, p=1.0e-05), foraging (𝜒2= 62.6, df= 24, p= 2.7e-05), running (𝜒2= 83.0, df= 24, p= 2.0e-08), and sitting (𝜒2= 60, df=24, p= 6.4e-05).
##work on behaviors of the squirrels by location
##tally behaviors using ones and zeros
squirrel_activities_location2018<- squirrel_behavior_2018%>%
select(Name, fur_color,running,chasing,climbing,eating,foraging)%>%
mutate(running=ifelse(running=="TRUE", 1,0),
chasing=ifelse(chasing=="TRUE",1,0),
climbing=ifelse(climbing=="TRUE",1,0),
eating=ifelse(eating=="TRUE",1,0),
foraging=ifelse(foraging=="TRUE",1,0))
##generate some descriptives as well as the sums of behaviors seen
general_behaviors_2018<-squirrel_activities_location2018%>%
dplyr::summarise(chasing=sum(chasing),climbing=sum(climbing),
eating=sum(eating),foraging=sum(foraging),running=sum(running),
prop_chasing=sum(chasing)/3023,
prop_climbing=sum(climbing)/3023, prop_eating=sum(eating)/3023, prop_foraging=sum(foraging)/3023,
prop_running=sum(running)/3023)
##analyze by nearest restaurant
behavior_counts_2018<-squirrel_activities_location2018%>%
group_by(Name)%>%
dplyr::summarise(running=sum(running),chasing=sum(chasing),
climbing=sum(climbing),
eating=sum(eating),foraging=sum(foraging))
Behaviors by Fur Color
To determine whether squirrel primary fur color played a role in the behavior identified, the distribution of fur colors by behavior were analyzed (code below shown for 2018 data). For the 2018 census data (Figure 6A), distributions of fur color differed significantly for each of the behaviors: chasing (𝜒2= 546.9, df= 3, p< 2.2e-16), climbing (𝜒2= 1115, df= 3, p< 2.2e-16), eating (𝜒2= 1325, df= 3, p< 2.2e-16), foraging (𝜒2= 2586, df= 3, p< 2.2e-16), running (𝜒2= 1285, df= 3, p< 2.2e-16). Similarly, 2020 census data (Figure 6B) also showed significant differences for all behaviors: chasing (𝜒2= 58.2, df= 2, p= 2.3e-13), climbing (𝜒2= 136, df= 2, p< 2.2e-16), digging (𝜒2= 10, df= 2, p= 7.0e-03), eating (𝜒2= 100.4, df= 2, p< 2.2e-16), foraging (𝜒2= 116.6, df= 2, p< 2.2e-16), running (𝜒2= 75.4, df= 2, p< 2.2e-16), sitting (𝜒2= 10.4, df= 2, p= 5.0e-03).
##this analysis uses the same tallied dataframe as above
fur_behavior_2018<-squirrel_activities_location2018%>%
group_by(fur_color)%>%
dplyr::summarise(running=sum(running),chasing=sum(chasing),
climbing=sum(climbing),
eating=sum(eating),foraging=sum(foraging))
Squirrel Location (Ground or Above Ground) by Restaurant and Fur Color
Additional information provided for each of the squirrels identified in the 2018 and 2020 censuses was the location seen, whether at ground level or above (e.g., in a tree, on wires, or on a building). Analysis shifted to determining whether the restaurants and their surrounding environments were influencing sightings and whether the squirrels’ fur color impacted their sighted location. In 2018 census data (Figure 7A) ground level and above ground sightings differed significantly between restaurants for both above ground (𝜒2= 440.3, df= 12, p< 2.2e-16) and ground level (𝜒2= 999.3, df= 12, p< 2.2e-16) sightings. Similar differences were seen in 2020 census data (Figure 7B), and these differences were statistically significant for both above ground (𝜒2= 176, df= 25, p< 2.2e-16) and ground level (𝜒2= 135.4, df= 25, p< 2.2e-16).
When above ground and ground level sightings were analyzed for differences in primary fur color. For 2018 data (Figure 7C), there were significant differences by fur color for both above ground sightings (𝜒2= 964.2, df= 2, p< 2.2e-16) and ground level sightings (𝜒2= 2368, df= 2, p< 2.2e-16). Likewise, for 2020 census data (Figure 7D) significant differences were seen for above ground (𝜒2= 149.2, df= 2, p< 2.2e-16) and ground level sightings (𝜒2= 318.2, df= 2, p< 2.2e-16).
##does ground v. tree position differ by location or color?
#filter out rows lacking locations&tally up by location
squirrel_trees_restuarants_2018<-squirrel_behavior_2018%>%
filter(!is.na(squirrel_behavior_2018$location))%>%
select(c(Name,fur_color,location))%>%
group_by(Name)%>%
dplyr::summarise(total_sightings=sum(1),
Above_ground=sum(ifelse(location=="Above Ground",1,0)),
ground_level=sum(ifelse(location=="Ground Plane",1,0)))%>%
select(Name,Above_ground,ground_level)%>%
arrange(Name)
#filter out rows lacking ground info and tally by color
squirrel_trees_color_2018<-squirrel_behavior_2018%>%
filter(!is.na(fur_color))%>%
group_by(fur_color,location)%>%
tally()%>%
mutate(above_ground=ifelse(location=="Above Ground",n,NA),
ground_level=ifelse(location=="Ground Plane",n,NA))%>%
select(fur_color,above_ground,ground_level)%>%
arrange(fur_color)
Comparison of Behaviors Between 2018 and 2020
Given the COVID-19 pandemic ant its impact on human life in NYC in 2020, whether squirrel behavior differed because of this global pandemic was of interest. This was accomplished by combining datasets from 2018 and 2020 as shown in code segment below, joining on shared behaviors. Aside from the noticeably fewer squirrels identified in 2020 (Figure 8), behavior distribution differed between 2018 and 2020 (𝜒2=23.3, df=4, p=0.0001).
##compare behaviors in 2018 & 2020
behavior_by_year<-dplyr::full_join(general_behaviors_2018, general_behaviors_2020,
by=c('chasing','climbing','eating','foraging','running'))
Discussion
Overall, this work highlights the interconnectedness of the human and nonhuman inhabitants of NYC. This is a delicate and unique ecosystem, informing ways in which the space is co-inhabited. Considering the 2018 squirrel dataset focused in Central Park, there were notable differences restaurants which were closest to the most squirrels. However, this was not reflected in a relationship between slice of a pizza slice and squirrel density. Differences in fur color distributions and behaviors in this dataset are likely reflective of overall distributions of these characteristics among the squirrels. More squirrels were seen at ground level than above ground, due in part to the crowd reported sightings as well as dense tree cover in portions of Central Park. More unique fur color may have aided in sightings among the trees or even on the ground. Presumably, a higher proportion of ground level sightings around restaurants are expected since the same individuals who are reporting the squirrels may be consuming the pizza.
2020 data showed many of the same patterns as the 2018 data did with one notable exception, there was a strong negative correlation between pizza slice price and quantity of squirrels seen. This may be due in part to a wider area of the census and, consequently, a greater diversity of prices. Alternatively, perhaps less expensive pizza is consumed by more individuals leading to more squirrels sighted. Less likely, but still plausible is that this pizza is in some way less desirable to the human consumer and discarded at an increased rate, providing a greater food resource for the squirrels.
When comparing 2018 and 2020 squirrel data, there are differences between the distribution of behaviors. The COVID-19 pandemic and subsequent lockdown had massive effects on human inhabitants of cities worldwide. NYC as the epicenter of the virus in the first waves in the United States had significant changes in its populace and businesses, and presumably, its nonhuman inhabitants as well. It has been noted that rat populations underwent decreases during the pandemic and a redistribution from previously heavily inhabited urban centers to localized sightings near alternative food sources.1 These changes are thought to be due to the commensal relationship between humans and rats. It is possible that a similar relationship underlies the changes in squirrel behaviors seen comparing these datasets.
While an interesting first look at these previously unassociated data sets, there are some limitations including the use of a pizza dataset that was generated in 2023. Conceivably some of the restaurants present in 2018 or 2020 were not included in this dataset and some restaurants on the list may not have existed in the years during which the squirrel data was collected. Squirrel data did have some missing variables which required exclusion of some squirrels from analysis, though numbers were large enough that it should not have an impact on final results.
Many questions arise from this data including whether menu offerings at different restaurants effect squirrel distribution. It may be interesting to explore this question in the future and attempt to determine if these varieties of options drive influx of human and squirrel patrons to these restaurants. With the knowledge that rat populations were significantly altered, it would be worthwhile if the squirrel populations in NYC showed parallel changes.
Ultimately, these data represent an interesting look at NYC, pizza, squirrels, and people. The changes seen during the height of the COVID-19 first wave in NYC underscore that this is a complex relationship that deserves additional attention and exploration.
Citations
Michael H Parsons, Jonathan L Richardson, Yasushi Kiyokawa, Rafal Stryjek, Robert M Corrigan, Michael A Deutsch, Masato Ootaki, Tsutomu Tanikawa, Faith E Parsons, Jason Munshi-South, Rats and the COVID-19 pandemic: considering the influence of social distancing on a global commensal pest, Journal of Urban Ecology, Volume 7, Issue 1, 2021, juab027, https://doi.org/10.1093/jue/juab027
R Core Team (2023). _R: A Language and Environment for Statistical Computing_. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.