Formal Report

Accessing the Data

library(tidyverse)
library(RSocrata)
library(lubridate)

if(!requireNamespace("devtools")) install.packages("devtools")
devtools::install_github("dkahle/ggmap", ref = "tidyup")
ggmap(get_googlemap())
register_google(key = "insert your google api key here")

## Install the required package with:
## install.packages("RSocrata) if not already installed

nyc_crime = read.socrata(
  "https://data.cityofnewyork.us/resource/9s4h-37hy.json",
  app_token = "xxxx",
  email     = "xxxx",
  password  = "xxxx" )

The raw, acquired dataset has 6,036,805 observations and 35 variables. Broadly speaking, the variables contain information on the exact date, time and location of crime, description of crime, demographic information of the victim and suspect, and police department information.

Subsetting Data

We were interested in 2014 to 2017 data pertaining to sex-related, weapons-related and drug-related felony. Therefore, we extracted this information from the data.

Below is a list of the penal codes in the dataset that we used to select the felonies of interest. Following that is the code we used to subset our data.

#This code was used to retrieve and assess the observations with missing date of crime 
#occurrence

crime_missing_cp = nyc_crime %>% 
  summarize_all(funs(sum(is.na(.))))

#Cleaning and subsetting the data

nyc_felony_crimes = nyc_crime %>% 
  janitor::clean_names() %>% 
  mutate(year = year(cmplnt_fr_dt))%>%
  mutate_if(is.character, tolower) %>% 
  filter(year %in% 2014:2017) %>% 
  filter(law_cat_cd == "felony") %>% 
  select(- station_name, - transit_district, - hadevelopt, - patrol_boro, - housing_psa, 
         - juris_desc)

saveRDS(nyc_felony_crimes, file = "./datasets/nyc_felony_crimes.rds")

#Selecting crimes of interest
sex_drug_weapons = nyc_felony_crimes %>% 
  filter(pd_cd %in% c(178, 694, 697, 176, 180, 153, 157, 177, 168, 159, 166, 164, 179, 155,
                      586, 696, # Sex related felony crimes
                      
                      ## Drug related felony crimes
                      500, 501, 502, 503, 505, 507, 510, 512, 514, 515, 519, 520, 521, 523,
                      524, 529, 530, 531, 532, 568, 570,
                      
                      ### Weapon related felony crimes
                      781, 792, 793, 796)) %>% 
  
  mutate(longitude = as.numeric(longitude),
         latitude = as.numeric(latitude)) %>% 
  
  select(cmplnt_num, boro_nm, cmplnt_fr_dt, cmplnt_to_dt, cmplnt_fr_tm, latitude, longitude,
         ky_cd, ofns_desc, pd_cd, pd_desc, vic_race, vic_sex, vic_age_group, year, 
         prem_typ_desc) %>% 
  
  mutate(
             boro_nm = if_else(boro_nm == "staten island", "staten_island", boro_nm),
             
         crime_group = if_else(pd_cd %in% c(178, 694, 697, 176, 180, 153, 157, 177, 168, 159,
                                            166, 164, 179, 155, 586, 696), "Sex-Related", 
                               
                       if_else(pd_cd %in% c(500, 501, 502, 503, 505, 507, 510, 512, 514, 515,
                                            519, 520, 521, 523, 524, 529, 530, 531, 532, 568, 
                                            570), "Drug-Related", 
                               
                       if_else(pd_cd %in% c(781, 792, 793, 796), "Weapon-Related", pd_cd)))
         )

saveRDS(sex_drug_weapons, file = "./datasets/sex_drug_weapons.rds")

Reverse Geocoding

The dataset comes without zip code or neighborhood information; there are only coordinates of the exact location where the crime happened. This makes aggregation beyond clustering of the geographic points impossible. Therefore, we used Google Maps’ API to reverse geocode and get the exact address of the crime occurrence. We limited ourselves to 2017 data only to reduce the burden when using Google’s API. We then subsequently took out the zip codes for mapping. The full code on how we reverse geocoded the longitude and latitude information is below.

#### Felony - Drug Offenses
#### Sampling only year 2017 and selecting all felony drug rape reports

### Drug offenses in 2017
drug_2017 = sex_drug_weapons %>% 
            filter(crime_group == "Drug-Related" & 
            year == 2017)

### Sampling to reduce the dataset to avoid getting rejected by Google API
drug_2017_1 = drug_2017 %>% 
              sample_n(2290)

drug_2017_zip_1 = drug_2017_1 %>% 
  mutate(zip = map2(.x = longitude, .y = latitude, ~ revgeocode(c(lon = .x, lat = .y))))

saveRDS(drug_2017_zip_1, file = "datasets/drug_2017_zip_1.rds")

#### Getting the data not sample. Second half of the data
drug_2017_2 = anti_join(drug_2017, drug_2017_zip_1, by = "cmplnt_num")

#### Iterating to get the data
drug_2017_zip_2 = drug_2017_2 %>% 
  mutate(zip = map2(.x = longitude, .y = latitude, ~ revgeocode(c(lon = .x, lat = .y))))

saveRDS(drug_2017_zip_2, file = "datasets/drug_2017_zip_2.rds")

### Joining to have the complete dataset

complete_drug_zip_2017 = bind_rows(drug_2017_zip_1, drug_2017_zip_2) 

saveRDS(complete_drug_zip_2017, file = "datasets/complete_drug_zip_2017.rds")

#### Felony - Sex Offenses

#The same process will be used for the sex-related crimes. There are 1993 observations; 
#therefore, there's no need to split the data before reverse geocoding.

sex_2017 = sex_drug_weapons %>% 
           filter(crime_group == "Sex-Related" & 
           year == 2017)

complete_sex_zip_2017 = sex_2017 %>% 
  mutate(zip = map2(.x = longitude, .y = latitude, ~ revgeocode(c(lon = .x, lat = .y))))

saveRDS(complete_sex_zip_2017, file = "datasets/complete_sex_zip_2017.rds")

#### Felony - Weapon Offenses
weapon_2017 = sex_drug_weapons %>% 
              filter(crime_group == "Weapon-Related" & 
              year == 2017)

complete_weapon_zip_2017 = weapon_2017 %>% 
  mutate(zip = map2(.x = longitude, .y = latitude, ~ revgeocode(c(lon = .x, lat = .y))))

saveRDS(complete_weapon_zip_2017, file = "datasets/complete_weapon_zip_2017.rds")

#### 2017 Dataset With Zip Codes
#Dataset with the 3 crimes together with their zipcodes
all_sex_drug_weapon_zip_2017 = bind_rows(complete_drug_zip_2017, complete_sex_zip_2017, 
                                         complete_weapon_zip_2017) 

saveRDS(all_sex_drug_weapon_zip_2017, file = "datasets/all_sex_drug_weapon_zip_2017.rds")

Usage

After cleaning, subsetting and reverse geocoding, we saved the resulting dataset in rds format. You can download the resulting datasets from this Google Drive. The code on how to read the datasets can also be found in our repository file called acquiringdata.Rmd

### nyc_felony_crimes dataset
readRDS(file = "./datasets/nyc_felony_crimes.rds")

### sex_drug_weapons dataset
readRDS(file = "./datasets/sex_drug_weapons.rds")

### Drug offenses for 2017 with zipcodes
readRDS(complete_drug_zip_2017, file = "datasets/complete_drug_zip_2017.rds")

### Sex offenses for 2017 with zipcodes
readRDS(complete_sex_zip_2017, file = "datasets/complete_sex_zip_2017.rds")

### Weapon offenses for 2017 with zipcodes
readRDS(complete_weapon_zip_2017, file = "datasets/complete_weapon_zip_2017.rds")

### Weapons, sex, drug for 2017 with zipcodes
readRDS(file = "datasets/all_sex_drug_weapon_zip_2017.rds")

Data Dictionary

The list below provides all 17 variables in the datasets with their brief descriptions:

cmplnt_num: randomly generated persistent ID for each complaint

boro_nm: the name of the borough in which the incident occurred

cmplnt_fr_dt: exact start date of occurrence for the reported incident

cmplnt_to_dt: exact end date of occurrent for the reported incident

cmplnt_fr_tm: exact time of occurrence for the reported incident

latitude: midblock latitude coordinate for Global Coordinate System, WGS 1984, decimal degrees (EPSG 4326)

longitude: midblock longitude coordinate for Global Coordinate System, WGS 1984, decimal degrees (EPSG 4326)

ky_cd: three-digit offense classification code

ofns_desc: description of offense corresponding with key code

pd_cd: three-digit internal classification code (more granular than Key Code)

pd_desc: description of internal classification corresponding with PD code
vic_race: victim’s race description

vic_sex: victim’s sex description (D=Business/Organization, E=PSNY/People of the State of New York, F=Female, M=Male)

year: year the incident occurred

prem_typ_desc: specific description of premises; grocery store, residence, street, etc.

crime_group: sex-related felony offenses, drug-related felony offenses, weapon-related felony offenses

zip: reverse geocoded address of the incident

Victim Demographics

For our first exploratory analysis, we analyzed demographic information for victims of sex-related, weapon-related and drug-related felonies.

felony_crime = readRDS(file = "./datasets/sex_drug_weapons.rds")

felony_graph =
  felony_crime %>%
  filter(!is.na(vic_age_group), vic_age_group != "unknown")

Note that in the table above, “People of the State of NY” indicates that the information was not provided by the city, presumably to protect the victims’ identity and privacy. Also of note is that young women are disproportionately affected by sex-relate felonies.

The code used to build the table above is shown below.

felony_crime = readRDS(file = "./datasets/sex_drug_weapons.rds")
overall_felony = readRDS(file = "./datasets/nyc_felony_crimes.rds")
complete_data = inner_join(overall_felony, felony_crime, by = "cmplnt_num") %>%
  janitor::clean_names(dat = .) 

complete_data %>%
  group_by(crime_group) %>%
  summarise(number = n()) %>%
  arrange(desc(number)) %>%
  knitr::kable()

complete_data %>%
  filter(crime_group == "Weapon-Related") %>%
  filter(!is.na(vic_race_y), vic_race_y != "unknown") %>% 
  group_by(vic_race_y) %>%
  summarise(number = n()) %>% arrange(desc(number)) %>% 
  knitr::kable() 

complete_data %>%
  filter(crime_group == "Weapon-Related") %>%
  filter(!is.na(vic_age_group_y), vic_age_group_y != "unknown") %>% 
  group_by(vic_age_group_y) %>%
  summarise(number = n()) %>% arrange(desc(number)) %>% 
  knitr::kable() 

complete_data %>%
  filter(crime_group == "Weapon-Related") %>%
  filter(!is.na(vic_sex_y), vic_sex_y != "unknown") %>% 
  group_by(vic_sex_y) %>%
  summarise(number = n()) %>% arrange(desc(number)) %>% 
  knitr::kable() 

complete_data %>%
  filter(crime_group == "Drug-Related") %>%
  filter(!is.na(vic_race_y), vic_race_y != "unknown") %>% 
  group_by(vic_race_y) %>%
  summarise(number = n()) %>% arrange(desc(number)) %>% 
  knitr::kable()

complete_data %>%
  filter(crime_group == "Drug-Related") %>%
  filter(!is.na(vic_age_group_y), vic_age_group_y != "unknown") %>% 
  group_by(vic_age_group_y) %>%
  summarise(number = n()) %>% arrange(desc(number)) %>% 
  knitr::kable() 

complete_data %>%
  filter(crime_group == "Drug-Related") %>%
  filter(!is.na(vic_sex_y), vic_sex_y != "unknown") %>% 
  group_by(vic_sex_y) %>%
  summarise(number = n()) %>% arrange(desc(number)) %>% 
  knitr::kable() 

complete_data %>%
  filter(crime_group == "Sex-Related") %>%
  filter(!is.na(vic_race_y), vic_race_y != "unknown") %>% 
  group_by(vic_race_y) %>%
  summarise(number = n()) %>% arrange(desc(number)) %>% 
  knitr::kable() 

complete_data %>%
  filter(crime_group == "Sex-Related") %>%
  filter(!is.na(vic_age_group_y), vic_age_group_y != "unknown") %>% 
  group_by(vic_age_group_y) %>%
  summarise(number = n()) %>% arrange(desc(number)) %>% 
  knitr::kable()

complete_data %>%
  filter(crime_group == "Sex-Related") %>%
  filter(!is.na(vic_sex_y), vic_sex_y != "unknown") %>% 
  group_by(vic_sex_y) %>%
  summarise(number = n()) %>% arrange(desc(number)) %>% 
  knitr::kable() 

f = female, m = male, e = people of new york (non-binary gender), d = business/organization

Length of Reported Incident by Borough and Felony

For our second exploratory analysis, We examined whether the length of time between when the crime started and ended differed by borough and felony type. Examining the length of time between when the crime started and ended can serve as a proxy indicator of the severity of the crime. Longer times may mean the crime is more severe, harder to resolve, more violent, and may require more resources to deal with. Furthermore, differences in the length of reported felonies may have implications for law enforcement officials, policymakers, and urban residents.

time_data = felony_crime %>%
  mutate(crime_group = forcats::fct_relevel(crime_group, "Drug-Related"),
         boro_nm = forcats::fct_relevel(boro_nm, "manhattan")) %>% 
  janitor::clean_names() %>% 
  mutate(time_diff2 = (as.numeric(cmplnt_to_dt - cmplnt_fr_dt, units = "days", na.rm = TRUE))) %>% 
  select(time_diff2, boro_nm, crime_group) %>% 
  mutate(time_diff2 = if_else(is.na(time_diff2), 0, time_diff2)) 
# if value is missing it means the crime doesn't have an end date. It happened at a point in time. 

time_data %>% rename(`Crime group` = crime_group) %>% 
group_by(`Crime group`) %>%
  summarise(Count = n(), 
            `Average time in days` = mean(time_diff2),
             `Standard deviation` = sd(time_diff2)) %>% knitr::kable()

time_data %>% filter(!is.na(boro_nm)) %>% 
  rename(Borough = boro_nm) %>% 
  group_by(Borough) %>%
  summarise(Count = n(), `Average time in days` = mean(time_diff2), 
            `Standard deviation` = sd(time_diff2)) %>% 
  knitr::kable(digits = 2) 

ggplot(data = time_data, aes(x = time_diff2)) +
  geom_density() +
    labs(
    title = "Density Plot of Length of Reported Felony",
    x = "Length of Reported Felony",
    y = "Density",
    caption = "Data from the crime_data dataset"
    ) +
  scale_x_continuous(breaks = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
                     labels = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10")) +
  xlim(0, 5)

crime_data <- time_data[complete.cases(time_data),] %>% 
    filter(time_diff2 > 0) 

crime_data %>%
  rename('Borough' = boro_nm) %>% 
  group_by(Borough) %>%
  summarise(mean_time = mean(time_diff2)) %>% knitr::kable()

Overall Felony Rates by Borough

For our third exploratory analysis, we chose to examine overall felony rates by borough. Doing so gives us a strong understanding as to how NYPD can prioritize their efforts when stemming crime; it also allows us to understand potential longitudinal trends for each borough. If, for example, felony rates have remained consistently high for a particular borough, laaw enforcement officials have suggestive evidence that the high rates are not a fluke and would be able to mobilize effects in that borough more effectively, readily, and directly.

To study this, we loaded information from the U.S. Census containing information on the total population in each of the boroughs to calculate felony rates. Next, we plotted the overall felony rates over 2014-2017 for each borough by merging the sex_drug_weapons dataset with the data from the U.S. Census.

sex_drug_weapons = readRDS(file = "./datasets/sex_drug_weapons.rds")

url = "https://www.census.gov/quickfacts/fact/table/newyorkcitynewyork,bronxcountybronxboroughnewyork,kingscountybrooklynboroughnewyork,newyorkcountymanhattanboroughnewyork,queenscountyqueensboroughnewyork,richmondcountystatenislandboroughnewyork/PST045217"

nyc_population = read_html(url) %>%  
  html_nodes(css = "table") %>% .[[1]] %>% 
  html_table(header = TRUE) %>% 
  as.tibble() %>% 
  janitor::clean_names()

names(nyc_population)[1:7] = c("estimate_date", "new_york_city", "bronx", "brooklyn", "manhattan", "queens", "staten_island")

nyc_population = nyc_population %>% 
  gather(key = boro_nm, value = population, estimate_date:staten_island) %>% 
  mutate(population = if_else(population == "Population estimates, July 1, 2017, (V2017)", "2017", population),
         population = as.numeric(gsub("," , "", population)))

grouped_df = sex_drug_weapons %>% 
  group_by(boro_nm, year) %>% 
  summarise(number = n())

full = left_join(grouped_df, nyc_population, by = "boro_nm") %>% 
  mutate(crime_rate = (number/population)*100)

x = full %>% 
  filter(!is.na(boro_nm)) %>% 
  ggplot(aes(x = year, y = crime_rate, color = boro_nm)) + 
  geom_point() + geom_line(size = 1) +
  labs(x = "Year",
       y = "Felony Rate",
       legend) + viridis::scale_color_viridis(
      name = "Borough", 
      discrete = TRUE
    ) + theme_classic()

ggplotly(x)

From the graph above, we notice some interesting findings. First, Bronx, by far, has had the highest overall felony rate from 2014-2017. Secondly, the felony rates for Brooklyn and Manhattan are very similar; this applies to Queens and Staten Island as well, which comprise the two boroughs with the lowest overall crime rates from 2014-2017. Third, we can see that felony rates decreased from 2016-2017 for all boroughs except for Brooklyn. The Bronx, in particular, has had a significant decrease in felony rates from 2016 to 2017. Lastly, overall felony rates across the years for all boroughs have tended to remain stable, generally speaking. As with the other exploratory analyses, note that these felonies are strictly sex-related, weapons-related, and drug-related felonies.

To examine this borough-specific trend more closely, we subsetted the graph by type of felony.

### Dataset for 2017 with addresses
crimes_df = readRDS(file = "./datasets/all_sex_drug_weapon_zip_2017.rds")

### Dataset for 4 years without addresses
crimes_14_17 = readRDS(file = "./datasets/sex_drug_weapons.rds")

### Selecting zipcodes from address
crimes_df = crimes_df %>% 
  mutate(zipcode = str_extract(zip, "[0-9][0-9][0-9][0-9][0-9]"))

stratify_df = crimes_14_17 %>% 
  group_by(boro_nm, crime_group, year) %>% 
  summarise(number = n())

population_df = left_join(stratify_df, nyc_population, by = "boro_nm") %>% 
  mutate(crime_rate = (number/population)*100)

population_df %>% 
  filter(!is.na(boro_nm)) %>% 
  ggplot(aes(x = year, y = crime_rate, color = boro_nm)) + 
  geom_point() + geom_line(size = 1) + facet_grid(~crime_group) +
  labs(x = "Year",
       y = "Felony Rate",
       legend) + viridis::scale_color_viridis(
      name = "Borough", 
      discrete = TRUE
    )

The stratification of overall felony rates over the years by type of felony (i.e., sex-related, drug-related, or weapons-related) provides a lot of information amongst the three graphs:

• Broadly speaking, felony rates have decreased across the board from 2016 to 2017 for nearly all boroughs - only Brooklyn has seen a rise in the rate of drug-related felonies from 2015 to 2017.

• Bronx has consistently had the highest crime rate over the last 4 years for all 3 types of felonies. We see a steep reduction, however, in the crime rate for the Bronx between the years of 2016 and 2017 - this change is most notable among the graph showing drug-relaed felonies.

• Drug-related and weapons-related felonies occur significantly more often compared to sex-related felonies. This may not be all that surprising, given the heavy interplay between drugs and weapons.

• Queens and Staten Island have had the lowest felony rates across all years and across all three types of felonies.

Geographical Distribution of Felonies

Our final exploratory analysis examined the geographical distribution of felonies, beyond the borough level. We thought this was critical to investigate because each borough contains such unique and vastly different neighborhoods comprising that borough - solely looking at associated factors at the neighborhood level ignores and masks these different identities within each neighborhood. Examining at this level also gives us more clear-cut information on the absolute numbers of felonies in each neighborhood and provides specific information regarding where law enforcement officials can prioritize efforts.

We first created a visual showing clusters and pockets of the 3 types of felonies in absolute numbers.

year_2014 = filter(crimes_14_17, year == 2014) %>% 
  filter(!is.na(longitude) | !is.na(latitude))

year_2015 = filter(crimes_14_17, year == 2015) %>% 
  filter(!is.na(longitude) | !is.na(latitude))

year_2016 = filter(crimes_14_17, year == 2016) %>% 
  filter(!is.na(longitude) | !is.na(latitude))

year_2017 = filter(crimes_14_17, year == 2017) %>% 
  filter(!is.na(longitude) | !is.na(latitude))

pal <- colorFactor(palette = c("blue", "red", "green"), 
               levels = c("Sex-Related", "Weapon-Related", "Drug-Related"))

leaflet()  %>% 
  addTiles() %>% 
  
#Setting view on NYC. The values are the coordinates of NYC
  setView(lat = 40.7, lng = -74.0, zoom = 11) %>%
  addCircleMarkers(data = year_2014,
                   lng = ~ longitude,
                   lat = ~ latitude, 
                   radius = 1,
                   color = ~pal(crime_group),
                   group = "2014",
                   clusterOptions = markerClusterOptions()) %>% 
  addCircleMarkers(data = year_2015,
                   lng = ~ longitude,
                   lat = ~ latitude, 
                   radius = 1,
                   color = ~pal(crime_group),
                   group = "2015",
                   clusterOptions = markerClusterOptions()) %>% 
  addCircleMarkers(data = year_2016,
                   lng = ~ longitude,
                   lat = ~ latitude, 
                   radius = 1,
                   color = ~pal(crime_group),
                   group = "2016",
                   clusterOptions = markerClusterOptions()) %>% 
  addCircleMarkers(data = year_2017,
                   lng = ~ longitude,
                   lat = ~ latitude, 
                   radius = 1,
                   color = ~pal(crime_group),
                   group = "2017",
                   clusterOptions = markerClusterOptions()) %>%
  addLayersControl(overlayGroups = c("2014", "2015", "2016", "2017")) %>% 
  addLegend(pal = pal, 
            values = c("Sex-Related", "Weapon-Related", "Drug-Related"),
            opacity = .5,
            title = "Crime Group",
            position = "topleft")

We decided to take the visual one step further by incorporating zip codes and creating heat maps to show prevalence of these 3 felony groups across NYC.

group_zip = crimes_df %>% 
  group_by(zipcode) %>% 
  summarize(number = n())

map_json <- geojson_read("./datasets/nyc_zip_map.geojson", what = "sp")

map_w_data = sp::merge(map_json, group_zip, by.x = "postalCode", by.y = "zipcode") 

map_w_data@data <-  
  map_w_data@data %>% 
  mutate(number = as.numeric(number)) %>% 
  mutate(number = if_else(is.na(number), 0, number))

total = tm_shape(map_w_data) +
  # Add title and change palette
  tm_fill(col = "number", 
          title = "Total Crime",
          palette = "Greens") +
  # Add tm_borders()
  tm_borders(col = "grey60", lwd = 0.5) 

drug_group = crimes_df %>% filter(crime_group == "Drug-Related") %>% 
  group_by(zipcode) %>% 
  summarize(number = n())

map_json_drug <- geojson_read("./datasets/nyc_zip_map.geojson", what = "sp")

map_w_drug = sp::merge(map_json_drug, drug_group, by.x = "postalCode", by.y = "zipcode") 

map_w_drug@data <-  
  map_w_drug@data %>% 
  mutate(number = as.numeric(number)) %>% 
  mutate(number = if_else(is.na(number), 0, number))

# Use qtm() to create a choropleth map of gdp
drug = tm_shape(map_w_drug) +
  # Add title and change palette
  tm_fill(col = "number", 
          title = "2017 Drug-Related Crimes",
          palette = "Greens") +
  # Add tm_borders()
  tm_borders(col = "grey60", lwd = 0.5)

weapon_group = crimes_df %>% filter(crime_group == "Weapon-Related") %>% 
  group_by(zipcode) %>% 
  summarize(number = n())

map_json_weapon <- geojson_read("./datasets/nyc_zip_map.geojson", what = "sp")

map_w_weapon = sp::merge(map_json_drug, weapon_group, by.x = "postalCode", by.y = "zipcode") 

map_w_weapon@data <-  
  map_w_weapon@data %>% 
  mutate(number = as.numeric(number)) %>% 
  mutate(number = if_else(is.na(number), 0, number))

# Use qtm() to create a choropleth map of gdp

weapon = tm_shape(map_w_weapon) +
  # Add title and change palette
  tm_fill(col = "number",
          title = "2017 Weapon-Related Crimes",
          palette = "Greens") +
  # Add tm_borders()
  tm_borders(col = "grey60", lwd = 0.5)
  # Add tm_credits()

sex_group = crimes_df %>% filter(crime_group == "Sex-Related") %>% 
  group_by(zipcode) %>% 
  summarize(number = n())

map_json_sex <- geojson_read("./datasets/nyc_zip_map.geojson", what = "sp")

map_w_sex = sp::merge(map_json_sex, sex_group, by.x = "postalCode", by.y = "zipcode") 

map_w_sex@data <-  
  map_w_sex@data %>% 
  mutate(number = as.numeric(number)) %>% 
  mutate(number = if_else(is.na(number), 0, number))

# Use qtm() to create a choropleth map of gdp
sex = tm_shape(map_w_sex) +
  # Add title and change palette
  tm_fill(col = "number",
          title = "2017 Sex-Related Crimes",
          palette = "Greens") +
  # Add tm_borders()
  tm_borders(col = "grey60", lwd = 0.5) 

current.mode <- tmap_mode("plot")
tmap_arrange(total, sex, drug, weapon)

tmap_mode(current.mode)

The heat maps provide an easy way of visually describing the absolute numbers of overall felonies, as well as felonies by each of the 3 types. Overall, we see that these crimes occur mostly in central Brooklyn and Bronx. This finding applies to both sex-related crimes and weapon-related crimes. In looking at the map for just drug-related crimes, we see that the crimes are distributed among 3 areas - central Brooklyn, south Brooklyn, and the Bronx.

Differences in Length of Crime

sex_drug_weapons = readRDS(file = "./datasets/sex_drug_weapons.rds")

test_mean_data = sex_drug_weapons %>%
  mutate(crime_group = forcats::fct_relevel(crime_group, "Drug-Related"),
         boro_nm = forcats::fct_relevel(boro_nm, "manhattan")) %>% 
  janitor::clean_names() %>% 
  mutate(time_diff2 = (as.numeric(cmplnt_to_dt - cmplnt_fr_dt, units = "days", na.rm = TRUE))) %>% 
  select(time_diff2, boro_nm, crime_group) %>% 
  mutate(time_diff2 = if_else(is.na(time_diff2), 0, time_diff2))

 leveneTest(time_diff2 ~ crime_group, data = test_mean_data) %>% 
  broom::tidy() %>% 
  knitr::kable()

kruskal.test(time_diff2 ~ crime_group, data = test_mean_data) %>% 
  broom::tidy() %>% 
  knitr::kable()

dunnTest(time_diff2 ~ crime_group,
              data = test_mean_data,
              #Method of adjustment - Benjamini-Hochberg               
              method = "bh")  

##                      Comparison           Z  P.unadj    P.adj
## 1    Drug-Related - Sex-Related -75.5637359 0.000000 0.000000
## 2 Drug-Related - Weapon-Related  -0.8243442 0.409744 0.409744
## 3  Sex-Related - Weapon-Related  75.8203110 0.000000 0.000000

Association Between Borough and Crime Occurrence

### Making Drug-related Felonies and Manhattan the Reference Groups
sex_drug_weapons = sex_drug_weapons %>%
  mutate(crime_group = forcats::fct_relevel(crime_group, "Drug-Related"),
         boro_nm = forcats::fct_relevel(boro_nm, "manhattan"))

with(sex_drug_weapons, table(boro_nm, crime_group)) %>%
  knitr::kable()

fit = multinom(crime_group ~ boro_nm, data = sex_drug_weapons) %>% 
  broom::tidy() %>% 
  mutate(`Odds Ratio` = exp(estimate))

clean_fit = fit %>% 
            mutate(term = str_replace(term, "\\(", ""),
               term = str_replace(term, "\\)", ""),
               term = str_replace(term, "boro_nm", ""),
               term = toTitleCase(term)) %>% mutate(term = str_replace(term, "_", " ")) %>% 
  rename(`Crime group` = y.level, Term = term, Estimate = estimate, `Standard error` = std.error, `P value` = p.value) 

clean_fit %>% 
  knitr::kable() 

ggplot(clean_fit, aes(x = Term, y = `Odds Ratio`)) +
  geom_segment( aes(x = Term, xend = Term, y = 0, yend = `Odds Ratio`), color = "grey") +
  geom_point(color = "orange", size = 3) + facet_grid(~ `Crime group`) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 80, hjust = 1, size = 8),
    panel.grid.major.x = element_blank(),
    panel.border = element_blank(),
    axis.ticks.x = element_blank()) +
  labs(x = "Borough",
       y = "Odds Ratio")