Rats-PreAnalysis.Rmd

---
title:  36-315 Final Project Data Pre-analysis 
author:  Arthur Jakobsson, Alex Cheng, Liz Chu, Kevin Ren 
date:  November 18, 2022 
output:
  html_document:
    toc: yes
    toc_float: yes
    code_folding: show
  pdf_document:
    toc: yes
urlcolor: blue
---

# Rat Sightings!

#Loading in our Libraries
Before we can do anything, we're going to load in all the libaries that we will use for the rest of the notebook. In the case that you find a new library you have to install, please do it in this large cell here rather than in assorted cells below.

```{r} 
# library imports
library(tigris)
library(dplyr)
library(leaflet)
library(tidyverse)
library(sp)
library(ggmap)
library(maptools)
library(broom)
library(httr)
library(rgdal)
library(gridExtra)
library(stringr)
library(ggseas)
library(geosphere)
library(stringr)
library(hydroTSM)
library(vcd)


library(magrittr)
library(knitr)
# distm(c(lon1, lat1), c(lon2, lat2), fun = distHaversine)
```


#Loading In Data

Our data will come from the 'Raw' folder, and will often come in the form of csv files. We will also do some map cleaning data in this cell so that it won't have to be done later. This cell should also onlybe loaded in once, and if there is new data to be loaded in please do so in this cell.
```{r, include=FALSE}
rats <- read.csv(file = 'Raw/Rat_Sightings.csv')
# # https://data.cityofnewyork.us/widgets/g642-4e55?mobile_redirect=true
health_inspection <- read.csv(file = 'Raw/DOHMH_New_York_City_Restaurant_Inspection_Results.csv')
# # https://data.cityofnewyork.us/Health/DOHMH-New-York-City-Restaurant-Inspection-Results/rs6k-p7g6
drinking_water_inspection <- read.csv(file = 'Raw/Self-Reported_Drinking_Water_Tank_Inspection_Results.csv')
# # https://data.cityofnewyork.us/Health/Self-Reported-Drinking-Water-Tank-Inspection-Resul/gjm4-k24g/data 

mayors <- read.csv(file = 'Raw/mayors.csv')

events <- read.csv(file = 'Raw/nycevents.csv')

# knicks <- read.csv(file = 'Raw/knicks.csv')
# # @ALEX
# 
# tax2019 <- read.csv(file = 'Raw/nyc2019taxreturns.csv')
# # @ALEX
# 
# pop2022 <- read.csv(file = 'Raw/nyc2022population.csv')
# # @ALEX
# 
# barLoc <- read.csv(file = 'Raw/nycbarlocations.csv')
# # @ALEX
# 
# dailyCovid <- read.csv(file = 'Raw/nycdailycovid.csv')
# # @ALEX
# 
# filmLoc <- read.csv(file = 'Raw/nycfilmlocations.csv')
# # @ALEX
# 
# zipcodenames <- read.csv(file = 'Raw/zipcodenames.csv')

subway_entrances <- read.csv(file = 'Raw/SubwayEntrances.csv')
subway_entrances <- subway_entrances %>% add_column(longitude = NA)
subway_entrances <- subway_entrances %>% add_column(latitude = NA)
lats = 1:nrow(subway_entrances)
longs = 1:nrow(subway_entrances)
for (i in (1:nrow(subway_entrances))) {
  longs[i] = as.numeric(sub(".", "", scan(text = subway_entrances$the_geom[i], what = "")[2]))
  lats [i] = as.numeric(str_sub(scan(text = subway_entrances$the_geom[i], what = "")[3],1,-2))
}
subway_entrances$longitude = longs
subway_entrances$latitude = lats

# 
# #NYC map lines
r <- GET('http://data.beta.nyc//dataset/0ff93d2d-90ba-457c-9f7e-39e47bf2ac5f/resource/35dd04fb-81b3-479b-a074-a27a37888ce7/download/d085e2f8d0b54d4590b1e7d1f35594c1pediacitiesnycneighborhoods.geojson')
nyc_neighborhoods <- readOGR(content(r,'text'), 'OGRGeoJSON', verbose = F)
# # https://rpubs.com/jhofman/nycmaps
# 
# 
# # clean up rat dataset (remove empty/useless columns)
# rats = subset(rats, select = -c(Agency, Agency.Name,Complaint.Type,  Descriptor ,  Vehicle.Type ,  Taxi.Company.Borough ,  Taxi.Pick.Up.Location , Bridge.Highway.Name ,  Bridge.Highway.Direction ,  Road.Ramp ,  Bridge.Highway.Segment, Facility.Type, Park.Facility.Name) )
# # clean up health inspection dataset (remove empty/useless columns)
health_inspection = subset(health_inspection, select = -c(Location.Point, Zip.Codes, Community.Districts, Borough.Boundaries, City.Council.Districts, Police.Precincts))
```




```{r}
# maps setup

#nice map (manhattan, brooklyn, queens, a bit of bronx)
left =-74.03 
bottom = 40.64
right = -73.87
top = 40.85
nyc_coords <- c(left, bottom, right, top)

#full map (all boroughs)
leftF = -74.2
bottomF = 40.55
rightF = -73.87
topF = 40.85
nyc_coordsF <- c(leftF, bottomF, rightF, topF)

#just dowtown manhattan
leftM = -74.03
bottomM = 40.69
rightM = -73.94
topM = 40.81
nyc_coordsM <- c(leftM, bottomM, rightM, topM)

nyc_map <- get_stamenmap(nyc_coords, maptype = "terrain", zoom = 11)
nyc_mapF <- get_stamenmap(nyc_coordsF, maptype = "terrain", zoom = 11)
nyc_mapM <- get_stamenmap(nyc_coordsM, maptype = "terrain", zoom = 11)

nyc_neighborhoods_df <- tidy(nyc_neighborhoods) # https://rpubs.com/jhofman/nycmaps

# @here I tried to make this map larger but I couldn't figure out how
nyc_mapPolygon <- ggmap(nyc_map) +
  geom_polygon(data=nyc_neighborhoods_df, aes(x=long, y=lat, group=group), color="blue", fill=NA)

nyc_mapFPolygon <- ggmap(nyc_mapF) +
  geom_polygon(data=nyc_neighborhoods_df, aes(x=long, y=lat, group=group), color="blue", fill=NA)

nyc_mapMPolygon <- ggmap(nyc_mapM) + 
  geom_polygon(data=nyc_neighborhoods_df, aes(x=long, y=lat, group=group), color="blue", fill=NA)

# nyc_mapPolygon
# nyc_mapMPolygon
# nyc_mapFPolygon
```

```{r}
# importing rats

ratSubset <- subset(rats, Longitude<right & Latitude<top & Longitude > left & Latitude >bottom)
ratSubsetF <- subset(rats, Longitude<rightF & Latitude<topF & Longitude > leftF & Latitude >bottomF)
ratSubsetM <- subset(rats, Longitude<rightM & Latitude<topM & Longitude > leftM & Latitude >bottomM)

health_inspection[, c("SCORE")] <- sapply(health_inspection[, c("SCORE")], as.integer)
health_inspection["SCORE"][is.na(health_inspection["SCORE"])] <- 0

healthSubset <- subset(health_inspection, Longitude<right & Latitude<top & Longitude > left & Latitude >bottom & SCORE>50)
healthSubsetF <- subset(health_inspection, Longitude<rightF & Latitude<topF & Longitude > leftF & Latitude >bottomF & SCORE>50)
healthSubsetM <- subset(health_inspection, Longitude<rightM & Latitude<topM & Longitude > leftM & Latitude >bottomM & SCORE>50)

rat_map <- ggmap(nyc_map) +
  geom_point(data=ratSubset, aes(x=Longitude, y = Latitude), alpha=0.2, size =0.01, color = "coral3")

rat_mapF <- ggmap(nyc_mapF) +
  geom_point(data=ratSubsetF, aes(x=Longitude, y = Latitude), alpha=0.2, size =0.01, color = "coral3")

rat_mapM <- ggmap(nyc_mapM) +
  geom_point(data=ratSubsetM, aes(x=Longitude, y = Latitude), alpha=0.2, size =0.01, color = "coral3")

# rat_map
# rat_mapF
# rat_mapM


ratHealthScore_map <- ggmap(nyc_map) +
  geom_point(data=ratSubset, aes(x=Longitude, y = Latitude), alpha=0.2, size =0.01, color = "chocolate3")+
  geom_point(data=healthSubset, aes(x=Longitude, y = Latitude, color = SCORE), size = 0.5, alpha=0.2) +
  scale_color_distiller(palette = "PiYG")

ratHealthScore_mapF <- ggmap(nyc_mapF) +
  geom_point(data=ratSubsetF, aes(x=Longitude, y = Latitude), alpha=0.2, size =0.01, color = "chocolate3")+
  geom_point(data=healthSubsetF, aes(x=Longitude, y = Latitude, color = SCORE), size = 0.5, alpha=0.2) +
  scale_color_distiller(palette = "PiYG")

ratHealthScore_mapM <- ggmap(nyc_mapM) +
  geom_point(data=ratSubsetM, aes(x=Longitude, y = Latitude), alpha=0.2, size =0.1, color = "chocolate3") +
  geom_point(data=healthSubsetM, aes(x=Longitude, y = Latitude, color = SCORE), size = 0.5, alpha=0.2) + 
  scale_color_distiller(palette = "PiYG") + 
  geom_polygon(data=nyc_neighborhoods_df, aes(x=long, y=lat, group=group), color="blue", fill="white", alpha=0.3)

```

```{r}
#HEAT MAP VERSION

#SUBWAY STATIONS
ggmap(nyc_map) + 
  geom_density_2d_filled(data=ratSubset, aes(x = Longitude, y = Latitude, fill = after_stat(level)), alpha = 0.4) + 
  geom_point(data=subway_entrances, aes(x=longitude, y = latitude), color="red",  size = 0.6, alpha=0.2) + 
  scale_color_distiller(palette = "PiYG")

#Restaurants (subset, score worse than 50)
ggmap(nyc_map) +
  geom_density_2d_filled(data=ratSubset, aes(x = Longitude, y = Latitude, fill = after_stat(level)), alpha = 0.4) + 
  geom_point(data=healthSubset, aes(x=Longitude, y = Latitude, color = SCORE), size = 0.3, alpha=0.1) + 
  scale_color_distiller(palette = "YlOrRd")
```

```{r}

# ratSubsetM <- subset(rats, Longitude<rightM & Latitude<topM & Longitude > leftM & Latitude >bottomM)

ratSubway_mapF <- ggmap(nyc_mapF) +
  geom_point(data=ratSubsetF, aes(x=Longitude, y = Latitude), alpha=0.05, size =0.1, color = "chocolate3") +
  geom_point(data=subway_entrances, aes(x=longitude, y = latitude), color="red",  size = 1, alpha=0.5) + 
  scale_color_distiller(palette = "PiYG") # + 
  # geom_polygon(data=nyc_neighborhoods_df, aes(x=long, y=lat, group=group), color="blue", fill=NA, alpha=0.3)

# ratHealthScore_map
# ratHealthScore_mapF
# ratHealthScore_mapM
ratSubway_mapF
```



```{r}
# system.time({
# minSubwayDist = function(longitude, latitude)
# {
#   min =  as.integer(.Machine$integer.max)
#   for (i in (1:nrow(subway_entrances)))
#   {
#     distance = distm(c(longitude, latitude), c(longs[i], lats[i]), fun = distHaversine)[1]
#     if(distance<min)
#     {
#       min=distance
#     }
#   }
#   return(min)
# }
# 
# ratSubsetM <- ratSubsetM %>% drop_na(Longitude)
# ratSubsetM <- ratSubsetM %>% drop_na(Latitude)
# # ratsub = rats[1:100,]
# ratSubsetM <- ratSubsetM %>% add_column(minSubwayDistance = NA)
# 
# mindists = 1:nrow(ratSubsetM)
# for(i in (1:nrow(ratSubsetM)))
# {
#   mindists[i] = minSubwayDist(ratSubsetM$Longitude[i], ratSubsetM$Latitude[i])
# }
# ratSubsetM$minSubwayDistance = mindists
# })
```

```{r}
#I want to now do a correlation between "region" and (rat count/area of region). nyc_neighborhoods_df is stored in the same way as states are in HW09.

#should also do proportion of "bad restaurants" to "good restaurants" per region
```


```{r}
# Family/Non-family Rats Visualization: subsetting
rats = subset(rats, rats$Borough != "Unspecified")
rats_fam = subset(rats, grepl("Family", rats$Location.Type, fixed = TRUE))
rats_fam = subset(rats_fam, (Borough == "BRONX" | Borough == "BROOKLYN" | Borough == "MANHATTAN" | Borough == "QUEENS" | Borough == "STATEN ISLAND"))
rats_fam$simple = ifelse(grepl("1-2", rats_fam$Location.Type, fixed = TRUE), "1-2 Family", "3+ Families")

rats_other = subset(rats, !grepl("Family", rats$Location.Type, fixed = TRUE))
rats_other = subset(rats_other, Location.Type != "Other (Explain Below)")
rats_other = subset(rats_other, Location.Type != "Street Area")
rats_other = subset(rats_other, (Borough == "BRONX" | Borough == "BROOKLYN" | Borough == "MANHATTAN" | Borough == "QUEENS" | Borough == "STATEN ISLAND"))
rats_other = rats_other %>% group_by(Location.Type) %>% filter(n() > 300 )
rats_other["Location.Type"][rats_other["Location.Type"] == "Vacant Building" | rats_other["Location.Type"] == "Vacant Lot"] <- "Unoccupied"
rats_other["Location.Type"][rats_other["Location.Type"] == "Government Building" | rats_other["Location.Type"] == "Commercial Building"] <- "Office Building"
```

```{r}
# fam/non-fam visualization
family = ggplot(data = rats_fam, aes(x = Borough)) +
geom_bar(aes(fill = simple)) +
labs(
title = "Rat Sightings Count for Family Dwellings",
x = "Borough",
y = "Count",
fill = "Number of Families"
)
non_family = ggplot(data = rats_other, aes(x = Borough)) + 
geom_bar(aes(fill = Location.Type)) + 
labs(
title = "Rat Sightings Count for Non-Family Dwellings",
subtitle = "Must have at least 300 sightings per location type",
x = "Borough",
y = "Count",
fill = "Location Type"
)
grid.arrange(family, non_family, nrow = 2)
```


```{r}
par(mar = c(5,4,1,10))

mosaicplot(table(rats_other$Borough, rats_other$Location.Type), main = "Mosaic Plot of Non-Family Location Types by Borough", shade=TRUE, las=2)
```


```{r}

# Liz's Graph 1 EDA
# rats is the unedited version of the dataset !!
# rats = read.csv(file = "Rat_Data.csv")
# rats = subset(rats, select = -c(Agency, Agency.Name, Complaint.Type, Descriptor,
#                                 Vehicle.Type, Taxi.Company.Borough, 
#                                 Taxi.Pick.Up.Location , Bridge.Highway.Name, 
#                                 Bridge.Highway.Direction, Road.Ramp, 
#                                 Bridge.Highway.Segment, Facility.Type, 
#                                 Park.Facility.Name))
private = c(
  "1-2 Family Dwelling",
  "1-2 Family Mixed Use Building", 
  "1-2 FamilyDwelling", 
  "1-3 Family Dwelling", 
  "1-3 Family Mixed Use Building", 
  "3+ Family Apartment Building", 
  "3+ Family Apt", 
  "3+ Family Apt.", 
  "3+ Family Apt. Building", 
  "3+ Family Mixed Use Building",
  "3+Family Apt.", 
  "Apartment",
  "Private House",
  "Residence",
  "Residential Building",
  "Residential Property",
  "Single Room Occupancy (SRO)"
)

commercial = c(
  "Cafeteria - Public School",
  "Catering Service", 
  "Commercial Building", 
  "Commercial Property", 
  "Construction Site", 
  "Day Care/Nursery", 
  "Government Building", 
  "Grocery Store",
  "Hospital",
  "Office Building", 
  "Restaurant", 
  "Restaurant/Bar/Deli/Bakery", 
  "Retail Store", 
  "School", 
  "School/Pre-School",
  "Store",
  "Street Fair Vendor", 
  "Summer Camp"
)

public = c(
 "Abandoned Building", 
 "Beach", 
 "Building (Non-Residential)", 
 "Catch Basin/Sewer",
 "Ground",
 "Parking Lot/Garage",
 "Public Garden", 
 "Public Stairs",
 "Street Area", 
 "Vacant Building", 
 "Vacant Lot", 
 "Vacant Lot/Property"
)

other = c(
  "",
  "N/A",
  "Other",
  "Other (Explain Below)"
)

rats$Location.Type[rats$Location.Type %in% private] <- "Private"
rats$Location.Type[rats$Location.Type %in% public] <- "Public"
rats$Location.Type[rats$Location.Type %in% commercial] <- "Commercial"
rats$Location.Type[rats$Location.Type %in% other] <- "Other"
rats$Location.Type <- factor(rats$Location.Type)
rats$Date = as.Date(rats$Created.Date, "%m/%d/%Y")

rats_per_dayL = rats %>% 
  group_by(Date, Location.Type) %>%
  tally() 

names(rats_per_dayL) = c("date", "locationtype", "n_rats")

# ggplot(data=rats_per_day, aes(x=date, y=n_rats, color=locationtype)) + 
#   geom_line(alpha=0.3) + labs(
#     title="Number of Rats Recorded on Each Day", 
#     subtitle="colored by the type of location",
#     x="Date", 
#     y="Number of Rats Recorded"
#   ) + 
#   scale_color_manual("Location Type", 
#                      values = c("Other" = "yellow", 
#                                 "Commercial" = "blue", 
#                                 "Private" = "red", 
#                                 "Public" = "green")) 

mayors$Date.Start <- as.Date(mayors$Date.Start)
mayors$Date.End <- as.Date(mayors$Date.End)
events$Date <- as.Date(events$Date)
events <- events %>%  filter(!row_number() %in% c(1, 2, 5))
```

```{r}
ggplot(rats_per_dayL) +
  geom_line(aes(date, n_rats, color=locationtype), alpha=0.3) + labs(
    title="Number of Rats Recorded on Each Day (with Mayors)", 
    subtitle="colored by the type of location",
    x="Date", 
    y="Number of Rats Recorded"
  ) + 
  scale_color_manual("Location Type", 
                     values = c("Other" = "yellow", 
                                "Commercial" = "blue", 
                                "Private" = "orange", 
                                "Public" = "green")) +
  geom_rect( 
    data = mayors,
    aes(xmin = Date.Start, xmax = Date.End, fill = Party),
    ymin = -Inf, ymax = Inf, alpha = 0.1
  ) +
  geom_vline(
    aes(xintercept = as.numeric(Date.Start)),
    data = mayors,
    colour = "grey50", alpha = 0.5
  ) +
  geom_text(
    aes(x = Date.Start+60, y = 110, label = Name),
    data = mayors,
    size = 3, vjust = 0, hjust = 0, nudge_x = 50, angle = 90) +
  geom_segment(data = events, aes(x = Date, y = 40, xend = Date, yend = 55), color = "red") +
  geom_text(data = events, aes(x = Date-50, y = 95, label = Event.Name), angle=90, color = "red")
  scale_fill_manual(values = c("blue", "red"))
```

dts <- strptime(time, "%I:%M:%S %p")
dts <- format(dts,format="%I:%M:%S %p")
actual_dates_vector <- c(actual_dates_vector,dts)

```{r}
dts <- strptime(rats$Time[1], "%I:%M:%S %p")
dts
substr(dts[1], 11,13)
```


```{r}
library(ggplot2)
rats$Time = word(rats$Created.Date, start = 2, end = 3)
rats.time = subset(rats, Time != "12:00:00 AM")
rats.time$Hour = as.numeric(substr(strptime(rats.time$Time, "%I:%M:%S %p"), 11,13))
ggplot(data=rats.time, aes(x=Hour)) +
  geom_bar(aes(fill=Location.Type)) +
  labs(x="Hour of the day", y="Number of rat reportings", title="Number of rat sightings by hour of day")
```

```{r}
library(ggplot2)
rats_per_day = rats %>%
  group_by(Date) %>%
  tally()

rats_per_day$Season = time2season(rats_per_day$Date, out.fmt = "seasons")
rats_per_day$Season = ifelse(rats_per_day$Season == "autumm", "autumn", rats_per_day$Season)

ggplot(data=rats_per_day, aes(x=Date, y=n, color = Season)) +
  geom_line(aes(Date, n), alpha=0.3) + labs(
    title="Number of Rats Recorded on Each Day",
    subtitle="colored by the season",
    x="Date",
    y="Number of Rats Recorded"
  ) +
stat_rollapplyr(color = "red", width = 30, align = "left", alpha = 0.5) +
ggtitle("Width = 30")
```
```{r}
dailyCovid <- read.csv(file = 'Raw/nycdailycovid.csv')
dailyCovidtest = mutate(dailyCovid, date_of_interest = as.Date(date_of_interest, format = "%m/%d/%Y"))
names(dailyCovidtest)[names(dailyCovidtest) == "date_of_interest"] <- "Date"
covid_rats = inner_join(x = dailyCovidtest, y = rats_per_day, by = "Date")
ggplot(data = covid_rats, aes(x = n, y = CASE_COUNT)) +
geom_point()
```

```{r}
#this is to remove outliers
Q_covid <- quantile(covid_rats$CASE_COUNT, probs=c(.25, .75), na.rm = FALSE)
iqr_covid <- IQR(covid_rats$CASE_COUNT)
Q_rats <- quantile(covid_rats$n, probs=c(.25, .75), na.rm = FALSE)
iqr_rats <- IQR(covid_rats$n)
covid_rats.iqr = subset(covid_rats,
                        CASE_COUNT < Q_covid[2] + 1.5*iqr_covid &
                        CASE_COUNT > Q_covid[1] - 1.5*iqr_covid &
                        n < Q_rats[2] + 1.5*iqr_rats &
                        n > Q_rats[1] - 1.5*iqr_rats)
#replotting without outliers
ggplot(data = covid_rats.iqr, aes(x = n, y = CASE_COUNT)) +
geom_point()
```
```{r}
raw = lm(CASE_COUNT~n, data = covid_rats)
no_outliers = lm(CASE_COUNT~n, data = covid_rats.iqr)
summary(raw)
summary(no_outliers)
```

```

# DO NOT RUN THIS CHUNK! saving this here just to show how the rats with districts are tagged

findDist <- function(lon, lat) {
  finalres = "COULD NOT IDENTIFY"
  for (x in 1:59) {
    tib = shape[x, 1]
    dist = tib$district
    poly = tib$geometry
    pt = st_point(c(lon, lat))
    res = summary(st_within(pt, poly))
    
    
    if (res[1] == "1") {
      finalres <- dist
    } 
  }
  
  return(finalres)
}

for (i in 1:range(nrow(rats2017))[1]) {
  rats$District[i] = findDist(rats$Longitude[i], rats$Latitude[i])
}
rats = subset(rats, District != "COULD NOT IDENTIFY")
write.csv(rats, "rats_district_tagged.csv", row.names=FALSE)


```

```{r}
# Kevin's Bar chart of rats given borough
borough.counts <- as.data.frame(table(rats$Borough))
names(borough.counts) = c("Borough","Count")
borough.counts <- rownames_to_column(borough.counts)
borough.counts <- borough.counts %>%  filter(!row_number() %in% c(1))
borough.counts
ggplot(data = borough.counts, aes(x=Borough, y=Count)) +
  geom_col(aes(fill=Borough))
```

```{r}
# reads in tax return dataset
tax_returns <- read.csv(file="Raw/nyctax-cleaned.csv")
tax_returns <- subset(tax_returns, Zip.Code != 0)
tax_returns <- subset(tax_returns, Zip.Code != 99999)
```

```{r}
incidents <- rats$Incident.Address
incidents <- removeNumbers(incidents)
head(incidents)
```


```{r}
library(tm)
#first, let's work with the titles
#redefine as Corpus
boroughs = Corpus(VectorSource(rats$Incident.Address))
#change to lowercase
boroughs <- tm_map(boroughs, content_transformer(tolower))
# Remove punctuations
boroughs <- tm_map(boroughs, removePunctuation)
# Eliminate extra white spaces
# boroughs <- tm_map(boroughs, stripWhitespace)
dtm.boroughs <- DocumentTermMatrix(boroughs, control = list(stopwords = TRUE, stemming = TRUE))
boroughs.words <- colnames(dtm.boroughs)
dtm.boroughs.mat = as.matrix(dtm.boroughs)
boroughs.freqs = colSums(dtm.boroughs.mat)
library(wordcloud)
wordcloud(words=boroughs.words, freq=boroughs.freqs, random.order=FALSE, max.words=100, colors=c("red", "blue", "black"))
```