Packages

As school of thought dictate, we should first install the packages to be used in this research by using the following syntax “install.packages(”The Package name“)” one after the other.

• Gist : Uncomment the codes below for a faster installation of all the packages used

# install.packages(c("flexdashboard", "DT", "pandocfilters", "googleVis", "R.rsp", "rmarkdown", "knitr", "plyr", "plotrix", "tidyr","dplyr", "ggplot2","plotly", "MASS","reshape2", "Amelia","mice","stringi","ROCR", "XML", "httr", "AUC", "magrittr", "rvest", "rgdal", "purrr", "DBI",
# "leaflet", "tigris", "RPostgreSQL" , "data.table", "tibble", "cowplot", "Hmisc", "sp", "maptools", "PerformanceAnalytics" , "ggcorrplot", "ggalt",
# "ggmap", "ggthemes", "readr", "webshot", "caTools", "caret", "rpart", "agricolae", "sqldf"))

suppressMessages(library("rmarkdown"))
#setwd("C:/Users/HP/Desktop")
library(flexdashboard)
library(DT)
library(prettydoc)
suppressMessages(library(R.rsp))
suppressMessages(library(pandocfilters))
suppressPackageStartupMessages(library(googleVis))
op <- options(gvis.plot.tag='chart')

Cont’d…

suppressMessages(library(knitr))
suppressMessages(library(plyr))
suppressMessages(library(plotrix))
suppressMessages(library(tidyr))
suppressMessages(library(dplyr))
suppressMessages(library(ggplot2))
suppressMessages(library(plotly))
suppressMessages(library(MASS))
suppressMessages(library(reshape2))
suppressMessages(library(Amelia))
suppressMessages(library(mice))
suppressMessages(library(stringi))
suppressMessages(library(ROCR))
suppressMessages(library(XML))
suppressMessages(library(httr))
suppressMessages(library(AUC))
suppressMessages(library(magrittr))
suppressMessages(library(rvest))
suppressMessages(library(rgdal))
suppressMessages(library(purrr))
suppressMessages(library(DBI))
suppressMessages(library(leaflet))
suppressMessages(library(tigris))
suppressMessages(library(RPostgreSQL))
suppressMessages(library(data.table))
suppressMessages(library(tibble))
suppressMessages(library(cowplot))
suppressMessages(library(Hmisc))
suppressMessages(library(sp))
suppressMessages(library(maptools))
suppressMessages(library(PerformanceAnalytics))
suppressMessages(library(ggcorrplot))
suppressMessages(library(ggalt))
suppressMessages(library(ggmap))
suppressMessages(library(ggthemes))
suppressMessages(library(readr))
suppressMessages(library(webshot))
suppressMessages(library(caTools))
suppressMessages(library(randomForest))
suppressMessages(library(caret))
suppressMessages(library(rpart))
suppressMessages(library(agricolae))
library(htmltools)
library(RColorBrewer);

PostgreSQL Connection.

It will be very difficult to work or combine all the needed dataset from all the above dataset. Therefore to make it a little bit easier, we will implore the help hand of SQL (Structured Query Language). A relational database like ORACLE SQL, MySQL or PostgreSQL would be nice in doing it better. I outed to use PostreSQL to access and extract the require data.

Connect your PostreSQL and R using the code below.

drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv, dbname = "capstone",
                 host = "localhost", port = 5432,
                 user = "postgres", password = "oracle")

Web Scrapping

The divorce came as a table from New York state health department website as follows.Here is the glimpse after cleaning: As mentioned above, this research project did used (heavily) relational database like PostgreSQL, MySQL, OracleSQL and SQLIte (Structured Query Language) for most of the query. Most (if not all) of the dataset used in the analysis are not readily available. I therefore sorted dataset from various government agencies. Some of the datase was scrapped for websites e.g the divorce dataset for year 2015 was scrapped from: https://www.health.ny.gov/statistics/vital_statistics/2015/table50.htm while others were downloaded and the required rows and columns were extracted.

div_2015 <- ("https://www.health.ny.gov/statistics/vital_statistics/2015/table50.htm") %>% 
  map(possibly(~html_session(.) %>% 
                 read_html %>% 
                 html_table(fill = TRUE) %>% 
                 .[[1]], 
               NULL))


div_2014 <- ("https://www.health.ny.gov/statistics/vital_statistics/2014/table50.htm") %>% 
  map(possibly(~html_session(.) %>% 
                 read_html %>% 
                 html_table(fill = TRUE) %>% 
                 .[[1]], 
               NULL))


div_2013 <- ("https://www.health.ny.gov/statistics/vital_statistics/2013/table50.htm") %>% 
  map(possibly(~html_session(.) %>% 
                 read_html %>% 
                 html_table(fill = TRUE) %>% 
                 .[[1]], 
               NULL))



div_2012 <- ("https://www.health.ny.gov/statistics/vital_statistics/2012/table50.htm") %>% 
  map(possibly(~html_session(.) %>% 
                 read_html %>% 
                 html_table(fill = TRUE) %>% 
                 .[[1]], 
               NULL))



div_2011 <- ("https://www.health.ny.gov/statistics/vital_statistics/2011/table50.htm") %>% 
  map(possibly(~html_session(.) %>% 
                 read_html %>% 
                 html_table(fill = TRUE) %>% 
                 .[[1]], 
               NULL))

div_2010 <- ("https://www.health.ny.gov/statistics/vital_statistics/2010/table50.htm") %>% 
  map(possibly(~html_session(.) %>% 
                 read_html %>% 
                 html_table(fill = TRUE) %>% 
                 .[[1]], 
               NULL))

# After scrapping the table, we realized that they are in list format. We will therefore convert it to accepatable format; first as a dataframe, subsetting unwanted rows and adding year column to the dataset as it was not provided in the scrapped table.

div_2010 <- as.data.frame(div_2010)
div_2011 <- as.data.frame(div_2011)
div_2012 <- as.data.frame(div_2012)
div_2013 <- as.data.frame(div_2013)
div_2014 <- as.data.frame(div_2014)
div_2015 <- as.data.frame(div_2015)


div_2010 <- div_2010[-c(1, 2,3,4,5, 11, 12), ]
div_2011 <- div_2011[-c(1, 2,3,4,5, 11, 12), ]
div_2012 <- div_2012[-c(1, 2,3,4,5, 11, 12), ]
div_2013 <- div_2013[-c(1, 2,3,4,5, 11, 12), ]
div_2014 <- div_2014[-c(1, 2,3,4,5, 11, 12), ]
div_2015 <- div_2015[-c(1, 2,3,4,5, 11, 12), ]

div_2010[["year"]] <- rep(2010, nrow(div_2010))
div_2011[["year"]] <- rep(2011, nrow(div_2011))
div_2012[["year"]] <- rep(2012, nrow(div_2012))
div_2013[["year"]] <- rep(2013, nrow(div_2013))
div_2014[["year"]] <- rep(2014, nrow(div_2014))
div_2015[["year"]] <- rep(2015, nrow(div_2015))
#The column names were not in conformitting with the standard R Programming column naming, we will therefore rename the columns.

Merge the divorce dataset from 2010 to 2015.

divorce <- Reduce(function(...) merge(..., all=T), list(div_2010, div_2011,
                            div_2012, div_2013, div_2014, div_2015))
colnames(divorce) <- c("county", "total", "cruelty", "abandonment", "imprisonment",
                       "adultery", "afterlegalseparation", "afterseparationbyagreement", "notstated", "year")
kable(head(divorce[, 1:8]))

Above is the scrapped and cleaned New York state divorce dataset from the year 2010 to 2015. Also, other dataset like the number Arrest, Mental Health, Police, Unemployment, Poverty, Crimes, Longtitude and Latitude was included in the analysis.

Arrest<- "https://docs.google.com/spreadsheets/d/e/2PACX-1vTQbIoMp7kHRauK90knVVWxL0ZQgSyJo5w53ufpeoJhv3NEKZ7e52Z9fU5gI2yj4Q50kCt2S9-WIe8t/pub?gid=0&single=true&output=csv"

LongLat <- "https://docs.google.com/spreadsheets/d/e/2PACX-1vRdUPece17flbIDQQ-prWnhh9ZOntvmHz9Saao-qC6EvV5VODzQ6TUfacf4oeB_w6EnDq0xmEhjw7X5/pub?gid=0&single=true&output=csv"

Mental <- "https://docs.google.com/spreadsheets/d/e/2PACX-1vT87ZxCxB9WAXFNmOWmxAV0aBlpyiE22Jqu_sgTtXpgpgHCzmskdboIiDtcJVgramrtVoTyqT2lzFeT/pub?gid=0&single=true&output=csv"


CrimesNY <- "https://docs.google.com/spreadsheets/d/e/2PACX-1vQ8KAFelT2l0AXFAabbbZSaeoqmQJRt-11kt7Jm_LKaHGbHa_ndUQidLj8Op_d19thmsyxoe9p5Dl3o/pub?gid=0&single=true&output=csv"

 LawEnf <- "https://docs.google.com/spreadsheets/d/e/2PACX-1vQauvaPEbqzKX1Vr9ojJuwqjYZFYTyMgclAOLxeUg9ithYryjbQz0foYmAe52jk3A7xBUl33vpDfksj/pub?gid=0&single=true&output=csv"

 
UnemploymentData <- "https://docs.google.com/spreadsheets/d/e/2PACX-1vSUjXqUWfUlOIQ_hYncSP8-BlOEVnzNQ5YDg49HD63IZDoZuw-2h8DiXt6hBuHTKeGAa9F2l4c23fwN/pub?gid=0&single=true&output=csv"


PovData <- "https://docs.google.com/spreadsheets/d/e/2PACX-1vSpp6mizldwEOSIXYPIuLH4Y8Ho3WvfkdAY8mTCcNnOwFExK7xo4q7bv6qUW-v2Z49rZIHJxo-3iG1X/pub?gid=0&single=true&output=csv"

Read datasets into R

# Note, the dataset contains column(s) of data we needed. We will load all of them and rename them to meet up with the standard.

nypd <- read.csv("C:/Users/HP/Documents/GitHub/Capstone_Project_2/NYPD_Complaint_Data_Historic.csv",
                 sep = ",", header=TRUE, stringsAsFactors = FALSE)


arrest <- read.csv(url(Arrest),header = TRUE, sep = ",", col.names= c("county", "year", "total",
                                          "felonytotal", "drugfelony","violentfelony", "dwifelony",
                                          "otherfelony", "misdemeanortotal", "drugmisd",
                                         "dwimisd", "propertymisd", "othermisd"))



mental_health <- read.csv(url(Mental), header = TRUE, sep = ",", col.names= c("rowcreateddatetime",
                          "year", "omhregioncode", "omhregionlabel", "county", "agegroup", "ratecodegroup", "countycode",
                          "countOfrecipientbycounty", "total", "paidclaimtotal"))



police <- read.csv(url(LawEnf), header = TRUE, sep = ",", stringsAsFactors = FALSE,
                   col.names=  c("county", "year", "policetotal"))



unemployment <- read.csv(UnemploymentData, header = TRUE, sep = ",", stringsAsFactors = TRUE,
                     col.names= c('year',   'month',    'county', 'total', 'benefitamountsdollars'))


lat_long <- read.csv(url(LongLat), header = TRUE, sep = ",")


crimes_ny <- read.csv(url(CrimesNY),header = TRUE, sep = ",")
crimes_ny <- crimes_ny[, 1:7]
colnames(crimes_ny)= c('county','year','population','indexcount', 'violentcount','propertycount','firearmcount')



poverty <- read_csv(url(PovData), locale = locale(encoding = "UTF-8"))
pov_sub <- poverty[, c(1, 2, 3, 4, 6, 15, 24)]
colnames(pov_sub) <- c("year","state", "countyid", "county", "allageproverty", "ueighteenpoverty", "fivetoseventenpov")

NYPD HISTORICAL COMPLAINTS

The New York City Department of Police, NYPD historical dataset that was downloaded from the city open data source was first analyzed as the stratum of the entire state as the NYC dataset contains all the require variables for further scrutiny. Below is the glimpse of the dataset

• Download The NYPD Historical Complaint Here!

• You’re required to download the NYPD Historical Complaint data and change path below to where you downloaded it. I can not include it here, its about 1.4GB in size!

• NB: The NYPD dataset is very large, i will therefore divide the table into chuck so are to fit the screen as (option=width) does not seem to work!

kable(head(nypd[, 1:8]))

print(paste0(".....Last columns"))

## [1] ".....Last columns"

kable(head(nypd[, 17:24]))

The arrest table…

kable(head(arrest[, 1:8]))

The mental health table

DT::datatable(mental_health, options = list(
  pageLength = 10
))

• NYPD Historical Complaint Data

kable(head(police));

Change the columns to months after dcast. Here is Unemployment table after cleaning..

# The unemployment dataset will be reshaped and have their columns renamed appropriately.

unemployment2 <- unemployment[, -5]
unemployment_dcast <- dcast(unemployment2, year + county ~ month, value.var="total")

names(unemployment_dcast)[names(unemployment_dcast)== 1] <- "jan"
names(unemployment_dcast)[names(unemployment_dcast)== 2] <- "feb"
names(unemployment_dcast)[names(unemployment_dcast)== 3] <- "mar"
names(unemployment_dcast)[names(unemployment_dcast)== 4] <- "apr"
names(unemployment_dcast)[names(unemployment_dcast)== 5] <- "may"
names(unemployment_dcast)[names(unemployment_dcast)== 6] <- "jun"
names(unemployment_dcast)[names(unemployment_dcast)== 7] <- "jul"
names(unemployment_dcast)[names(unemployment_dcast)== 8] <- "aug"
names(unemployment_dcast)[names(unemployment_dcast)== 9] <- "sep"
names(unemployment_dcast)[names(unemployment_dcast)== 10] <- "oct"
names(unemployment_dcast)[names(unemployment_dcast)== 11] <- "nov"
names(unemployment_dcast)[names(unemployment_dcast)== 12] <- "dec"

unemployment_dcast <- tbl_df(unemployment_dcast)



DT::datatable(unemployment_dcast, options = list(
  pageLength = 10
))

New York State crimes count.

DT::datatable(crimes_ny, options = list(
  pageLength = 10
))

Crimes County By County Table

Load Data Into PostreSQL

Let start by checking if any of the dataset do exist in the database, and if it exist, we will remove them and write (or save) our datasets in the PostgreSQL Automatically.

dbRemoveTable(con, "arrest")

## [1] TRUE

dbWriteTable(con, "arrest", arrest, row.names=FALSE)

## [1] TRUE

dbRemoveTable(con, "mental_health")

## [1] TRUE

dbWriteTable(con, "mental_health", mental_health, row.names=FALSE)

## [1] TRUE

dbRemoveTable(con, "police")

## [1] TRUE

dbWriteTable(con, "police", police, row.names=FALSE)

## [1] TRUE

dbRemoveTable(con, "unemployment_dcast")

## [1] TRUE

dbWriteTable(con, "unemployment_dcast", unemployment_dcast, row.names=FALSE)

## [1] TRUE

dbRemoveTable(con, "divorce")

## [1] TRUE

dbWriteTable(con, "divorce", divorce, row.names=FALSE)

## [1] TRUE

dbRemoveTable(con, "crimes_ny")

## [1] TRUE

dbWriteTable(con, "crimes_ny", crimes_ny, row.names=FALSE)

## [1] TRUE

dbRemoveTable(con, "pov_sub")

## [1] TRUE

dbWriteTable(con, "pov_sub", pov_sub, row.names=FALSE)

## [1] TRUE

If your connection is right, you should see that the above table are now in your database. But some of the data type may be wrong, use the code below to alter the it to desirable data type.

dbGetQuery(con, statement = paste("ALTER TABLE crimes_ny
    ALTER COLUMN population TYPE INT USING population::integer"));

# dbGetQuery(con, statement = paste("ALTER TABLE inmates_dcast
#           ALTER COLUMN values TYPE VARCHAR USING values::VARCHAR"))
# 

Data Cleaning (SQL)

After writing (or saving) them into the database, we can now try to subset the columns and the row we needed from the tables.

police_sub <- dbGetQuery(con, statement = paste(
                       "SELECT year, county, policetotal as police",
                       "FROM police",
                       "WHERE year IN (2010, 2011,2012,2013,2014,2015)",
                      "ORDER BY year"));


arrest_sub <- dbGetQuery(con, statement=paste(
                           "SELECT *",
                           "FROM arrest",
                           "WHERE year IN (2010, 2011,2012,2013,2014,2015)",
                           "ORDER BY county ASC"));



divorce$total <- (gsub(",", "", divorce$total));
divorce_sub <-  dbGetQuery(con, statement=paste(
                           "SELECT total as divorce, county, year",
                           "FROM divorce",
                           "ORDER BY county ASC"));


emp <-  dbGetQuery(con, statement=paste(
                           "SELECT *",
                           "FROM unemployment_dcast",
                           "WHERE year IN (2010, 2011,2012,2013,2014,2015)",
                           "ORDER BY county ASC"));


crimes_ny_sub <-  dbGetQuery(con, statement=paste(
                           "SELECT county, year, population ",
                           "FROM crimes_ny",
                           "WHERE year IN (2010, 2011,2012,2013,2014,2015)",
                           "ORDER BY county ASC"));

Saving them in tbl for faster accessibility.

# For easy and better table ouput
police_sub <- tbl_df(police_sub)
arrest_sub <- tbl_df(arrest_sub)
divorce_sub <- tbl_df(divorce_sub)
crimes_ny_sub <- tbl_df(crimes_ny_sub)
#inmates_sub <- tbl_df(inmates_sub)

dbGetQuery(con, statement = paste("ALTER TABLE mental_health ALTER COLUMN total TYPE INT"));

mental_health_sub <- as.data.frame(dbGetQuery(con, statement = paste(
                       "SELECT year, county, agegroup, total as mentalhealth",
                       "FROM mental_health",
                       "WHERE year BETWEEN 2010 AND 2015",
                       "ORDER BY year, agegroup")));

Finally, after all the hard work done by PostreSQl, we now have our sought after dataset, hoooray!

Lets disconnect PostreSQL and start using R and SQLDF.

#dbListTables(con);

dbDisconnect(con)

## [1] TRUE

detach("package:RPostgreSQL", unload=TRUE)
library(sqldf);

DATA CLEANING (R and SQL)

Upon further scrutiny, we discovered that the poverty table still have some unwanted character in them, we therefore removed them using gsub from grep, a R based.

• Mental Health and Divorce.

mental_sub <- dcast(mental_health_sub, year +county ~ agegroup, value.var="mentalhealth", fun.aggregate = sum)

mental_sub <- mental_sub[, c(1:4)]

divorce_sub$divorce <- (gsub(",", "", divorce_sub$divorce));
divorce_sub1 <- sapply(divorce_sub$divorce, as.numeric)
divorce_sub2 <- divorce_sub[, c(2:3)]
div_extraction <- cbind.data.frame(divorce_sub2, divorce_sub1)
colnames(div_extraction) <- c("county", "year", "divorce")

• Poverty

povselect <- sqldf::sqldf("SELECT * FROM pov_sub 
                          WHERE year BETWEEN 2010 AND 2015 AND state = 36 AND county !='New York'
                          ORDER BY year ASC")

povselect$county <- (gsub("\\(NY)", "", povselect$county));
povselect$county <- (gsub("\\County", "", povselect$county));
povselect$allageproverty <- (gsub(",", "", povselect$allageproverty));
povselect$ueighteenpoverty <- (gsub(",", "", povselect$ueighteenpoverty));
povselect$fivetoseventenpov <- (gsub(",", "", povselect$fivetoseventenpov))
col_pov <- c('allageproverty', 'ueighteenpoverty', 'fivetoseventenpov')
povselect[col_pov] <- sapply(povselect[col_pov], as.numeric)
povertydata <- povselect[, -c(2:3)]

• Employment

emp$jan <- (gsub(",", "", emp$jan))
emp$feb <- (gsub(",", "", emp$feb))
emp$mar <- (gsub(",", "", emp$mar))
emp$apr <- (gsub(",", "", emp$apr))
emp$may <- (gsub(",", "", emp$may))
emp$jun <- (gsub(",", "", emp$jun))
emp$jul <- (gsub(",", "", emp$jul))
emp$aug <- (gsub(",", "", emp$aug))
emp$sep <- (gsub(",", "", emp$sep))
emp$oct <- (gsub(",", "", emp$oct))
emp$nov <- (gsub(",", "", emp$nov))
emp$dec <- (gsub(",", "", emp$dec))



emp1 <- sapply(emp[, c(3:14)], as.data.frame.numeric)
emp2 <- as.data.frame(emp1)


unemp_sub <- cbind.data.frame(emp2,emp[, c(1,2)], stringsAsFactors=FALSE)
colnames(unemp_sub) <- c("jan", "feb", "mar", "apr", "may", "jun", "jul", "aug",
                         "sept", "oct", "nov", "dec", "year", "county")

df_unemp_sub <-sqldf::sqldf("SELECT county, year, (jan + feb + mar + apr + may +
                           jun + jul + aug + sept + oct + nov + dec) as unemployment
                           FROM unemp_sub")

• Converting Character Columns To Date.

nypd[, 2] <- as.Date.character(nypd[, 2], format = "%m/%d/%Y")
nypd[, 4] <- as.Date.character(nypd[, 4], format = "%m/%d/%Y")
nypd[, 6] <- as.Date.character(nypd[, 6], format = "%m/%d/%Y")

Data cleaning of NYPD’s latitude and Longtitude from coordinate to a usable R format.

nypd$Lat_Lon <- (gsub(")","", nypd$Lat_Lon))
nypd$Lat_Lon <- (gsub(",",":", nypd$Lat_Lon))
nypd$Lat_Lon <- (gsub("[[alphanum:]]","", nypd$Lat_Lon))

nypd_new <- subset(nypd, CMPLNT_FR_DT >= '2010-1-1' & CMPLNT_FR_DT <= '2015-12-31', 
                   select=c(CMPLNT_NUM , CMPLNT_FR_DT, CMPLNT_FR_TM, RPT_DT,  OFNS_DESC,
                            CRM_ATPT_CPTD_CD, LAW_CAT_CD,BORO_NM,  Latitude, Longitude,
                            HADEVELOPT, ADDR_PCT_CD, PREM_TYP_DESC))

colnames(nypd_new) <- c("ComplainNum", "IncidentDate", "IncidentTime",
                        "ReportedDate", "OffenseDesc", "CrimeFailSuccces", "CrimeLevel",
                        "Borough", "lat", "Long", "NYCHA", "PolicePrecinct", "BuildingType" )

nypd_new<- as.data.frame(nypd_new)

• Converting the NYPD data to days of week, months and even years to easy accessibility.

nypd_new$DateTime = paste(nypd_new$IncidentDate, nypd_new$IncidentTime, sep="")

nypd_new$dayofweek <- strftime(nypd_new$IncidentDate, '%A') # Days of the week

nypd_new$months <- strftime(nypd_new$IncidentDate, '%B') # Converting To Months Column

nypd_new$year <- strftime(nypd_new$IncidentDate, '%Y') # Converting To Year Column

kable(head(nypd_new[, 10:17]))

Incident Report

• Writing Incident Report With SQL Query

 incidentreport <- sqldf::sqldf("SELECT
                            
                             'At ' || IncidentTime || ',' || dayofweek ||', '|| months || ' ' || year
                            ||', a complaint was filled with NYPD police precinct ' || PolicePrecinct ||
                           ' about ' || OffenseDesc
                         || ' that was committed . The complainant maintained that the incident occurred at a '||
                          BuildingType || ' in ' || Borough || ' with co-ordinate ' || lat|| ' and ' ||Long||
                         ' . The Police described the incident as ' || CrimeFailSuccces || ' and classified it as ' ||
                          CrimeLevel ||' with complain number: ' || ComplainNum
                           as Incident_Report from nypd_new
                                ")
 
 kable(head(incidentreport))

Inferential Statistics (NYC)

• Top Ten NYPD Precinct & Neigbhorhoods with High Crime Incidents From 2010 to 2015.

precinct_count <- plyr::count(nypd_new, vars = "PolicePrecinct")
precinct <- sqldf::sqldf("SELECT PolicePrecinct, freq as NumberOfCrimesHandled
                         FROM precinct_count ORDER BY freq DESC LIMIT 10")
precinct_low <- sqldf::sqldf("SELECT PolicePrecinct,  freq as NumberOfCrimesHandled
                             FROM precinct_count ORDER BY freq ASC LIMIT 10")
lat= c(40.67124, 40.84156, 40.82885, 40.81687 , 40.75792, 40.86400 , 40.68564, 40.64174, 40.84865, 40.65247)
long= c(-73.88831, -73.86918, -73.92980, -73.91769,-73.99111 , -73.89342, -73.91079, -74.07568, -73.91168, -73.94426)

lat2 <- c(0,40.77472,40.58751,40.67687,40.58230,40.73599,40.75520,40.50930,40.81077,40.72350)
long2 <- c(0,-73.97468,-73.81171,-74.00806,-74.16894,-73.74580,-73.89195,-74.24090,-73.88334,-73.94862)

precinct <- cbind.data.frame(precinct, lat, long)
precinct_low <- cbind.data.frame(precinct, lat2, long2)

kable(head(precinct, 10))

MAPS

• Map showing all New York Counties, least and most crime area in New York City.

• click on the pop-up to see the county name.

neigbourhood_json <- GET('http://data.beta.nyc//dataset/0ff93d2d-90ba-457c-9f7e-39e47bf2ac5f/resource/35dd04fb-81b3-479b-a074-a27a37888ce7/download/d085e2f8d0b54d4590b1e7d1f35594c1pediacitiesnycneighborhoods.geojson')

neigbourhood_json2 <- GET('http://data.beta.nyc//dataset/0ff93d2d-90ba-457c-9f7e-39e47bf2ac5f/resource/35dd04fb-81b3-479b-a074-a27a37888ce7/download/d085e2f8d0b54d4590b1e7d1f35594c1pediacitiesnycneighborhoods.geojson')

• Precinct juridiction with low crime
• Precinct juridiction with high crime

neighborhoods_low <- readOGR(content(neigbourhood_json2,'text'), 'OGRGeoJSON', verbose = F)

precinct_low2 <- precinct_low
coordinates(precinct_low2) <- ~long2 + lat2
proj4string(precinct_low2) <- proj4string(neighborhoods_low)
combining_low <- over(precinct_low2, neighborhoods_low)
precint_low3 <- cbind(precinct_low, combining_low)
precint_low3 <- precint_low3[, -c(3,4, 8)]

map1 <- leaflet(neighborhoods_low) %>%
  addTiles() %>% 
  addPolygons(popup = ~neighborhoods_low, color = "#03F", fill = TRUE,
              fillColor = "brown", opacity = 0.5, fillOpacity = 0.2) %>% 
  addMarkers(~long2, ~lat2, popup = ~neighborhoods_low, data = precint_low3) %>%
  addProviderTiles("Stamen.Toner", options = providerTileOptions()) %>%
  setView(-73.92, 40.75, zoom = 9.5)




## Map 2

# Neigbourhood with high crime rate

 neighborhoods <- readOGR(content(neigbourhood_json,'text'), 'OGRGeoJSON', verbose = F)

precinct2 <- precinct
coordinates(precinct2) <- ~long + lat
proj4string(precinct2) <- proj4string(neighborhoods)
combining <- over(precinct2, neighborhoods)
precint3 <- cbind(precinct, combining)
precint3 <- precint3[, -8]

map2 <- leaflet(neighborhoods) %>%
  addTiles() %>% 
  addPolygons(popup = ~neighborhood, color = "#03F", fill = TRUE,
              fillColor = "brown", opacity = 0.5, fillOpacity = 0.2) %>% 
  addMarkers(~long, ~lat, popup = ~neighborhood, data = precint3) %>%
  addProviderTiles("HERE.normalday", options = providerTileOptions()) %>%
  setView(-73.92, 40.75, zoom = 10.5)


# New York state counties

#df_maps <- df_incident[1:62, ]  

counties_map <- leaflet(lat_long) %>%
  addTiles() %>%
  addMarkers(lng = lat_long$long, lat = lat_long$lat, popup = lat_long$county, data =lat_long) %>%
  addProviderTiles("Esri.WorldImagery", options = providerTileOptions()) %>%
  setView(-73.82, 40.75, zoom = 5.9)


counties_map

map1

map2

# Table (nyc high and low)
kable(precint_low3[, 1:6])

kable(precint3)

From analysis and map above, we can deduce that we have 5 out 10 locations with the highest number of crimes complaint in Bronx, 3 in Brooklyn, 1 each in Manhattan and Staten Island, but none in Queens! Does this have anything to do with poverty level?

• What Day Has Most Crime?

daysofweek <- plyr::count(nypd_new, vars = "dayofweek")
offenses <- plyr::count(nypd_new, vars = "OffenseDesc")
daysofweek[, -1] <- sapply(daysofweek[, -1], as.numeric)
daysofweek$percentages <- round(daysofweek$freq/sum(daysofweek$freq)*100)
 
days <- paste(daysofweek$dayofweek, daysofweek$percentages) # add percents to labels
days <- paste(days,"%",sep="") 
 
pie(daysofweek$freq,labels = days, col=rainbow(length(days)),
    main="Pie Chart of Crime Complaints In Days Of The Week") 

print(paste0("Most crimes are carried out on FRIDAYS with the total being: ", max(daysofweek$freq)))

## [1] "Most crimes are carried out on FRIDAYS with the total being: 462851"

print(paste0("While the least crimes were carried out on SUNDAYS with the total being: ", min(daysofweek$freq)))

## [1] "While the least crimes were carried out on SUNDAYS with the total being: 387899"

• What Month(s) Has Most Crimes?

theme_set(theme_bw())
incident_months <- plyr::count(nypd_new, vars = "months")
incident_months[, 2] <- sapply(incident_months[, 2],as.numeric)
incident_months2 <- sqldf::sqldf("SELECT months as MonthOfIncident, freq as IncidentFrequency
                                FROM incident_months ORDER BY freq DESC")
kable(head(incident_months2))

ggplot(incident_months2, aes(x=MonthOfIncident, y=IncidentFrequency)) + 
  geom_point(size=5) + 
  geom_segment(aes(x=MonthOfIncident, 
                   xend=MonthOfIncident, 
                   y=0, 
                   yend=IncidentFrequency)) + 
  labs(title="NYPD Crime Complaints Per Month ", 
       subtitle="FROM 2010 To 2015") + 
  theme(axis.text.x = element_text(angle=90, vjust=0.9)) + coord_flip()

The Most crimes was carried out in July follow by August and May in Third place , while February has the least numbers of crimes complaints.

• What Time Of The Day Has Most Crime?

incident_time <- plyr::count(nypd_new, vars = "IncidentTime")
incident_time <- sqldf::sqldf("SELECT IncidentTime as TimeOfIncident, freq as IncidentFrequency
                              FROM incident_time ORDER BY freq DESC")
kable(head(incident_time))

Most Crimes are carried out around 12 Noon, 3PM, 6PM and 4PM, while the least crimes are carried out in the early morning toward the morning.

• Misdemeanor Vs Felony

Crime_cat <- plyr::count(nypd_new, vars = "CrimeLevel")

sprintf(" Misdemeanor crimes are '%i' times greater than Felony crimes in New York City", round(Crime_cat[2, 2]/Crime_cat[1, 2]))

## [1] " Misdemeanor crimes are '2' times greater than Felony crimes in New York City"

• Location (Area/Place) Where Most Crime Occured?

Building_cat <- plyr::count(nypd_new, vars = "BuildingType")
Building_cat <- sqldf::sqldf("SELECT BuildingType as TypeOfBuilding, freq as IncidentFrequency
                             FROM Building_cat ORDER BY freq DESC LIMIT 20")

Building_cat[, -1] <- sapply(Building_cat[, -1], as.numeric)


Building_cat$percentages <- round(Building_cat$IncidentFrequency/sum(Building_cat$IncidentFrequency)*100)
days <- paste(Building_cat$TypeOfBuilding, Building_cat$percentages) # add percents to labels
days <- paste(days,"%",sep="") # ad % to labels
pie(Building_cat$IncidentFrequency,labels = days, col= c("purple", "violetred1", "green3","cornsilk", "cyan", "white"),
   main="Top 20 Building/Area Where Crimes Were Carried Out", border = FALSE, init.angle = 315) 

kable(head(Building_cat))

Crimes were most carried out on the street, followed by an apartmental houses, where residential houses came third!

At this time, we will be combining the required datasets into one and write it into PostgreSQL for further scrutiny.

df <- sqldf::sqldf("SELECT c.county, c.population, c.year,  a.felonytotal,a.misdemeanortotal,
a.total as totalcrimes,p.policetotal as lawenforcementtotal ,a.othermisd, d.divorce, m.child as mentalchild,
m.adult as mentaladult, (m.child+m.adult) as mentaltotal, u.unemployment ,a.violentfelony, a.drugfelony,
a.dwifelony, a.otherfelony, a.drugmisd, a.dwimisd, a.propertymisd FROM crimes_ny_sub c
LEFT JOIN arrest_sub a ON a.county=c.county AND a.year=c.year
LEFT JOIN div_extraction d ON d.year=c.year AND a.county=d.county
LEFT JOIN df_unemp_sub u ON u.year=a.year AND u.county=d.county
LEFT JOIN mental_sub m ON m.year=a.year AND m.county=a.county
LEFT JOIN police p ON p.year=a.year AND p.county=d.county
ORDER BY c.year, c.county ASC")

df <- sqldf::sqldf("SELECT * from df WHERE county NOT IN ('Out-of-State', 'St. Lawrence') ")
df <- cbind.data.frame(df, povertydata[, 3:5])

Inferential Statistics (NY)

A summary statistics of the datasets

summary(df)

##     county            population           year       felonytotal   
##  Length:372         Min.   :   4678   Min.   :2010   Min.   :   11  
##  Class :character   1st Qu.:  50431   1st Qu.:2011   1st Qu.:  326  
##  Mode  :character   Median :  89674   Median :2012   Median :  529  
##                     Mean   : 316077   Mean   :2012   Mean   : 2529  
##                     3rd Qu.: 233454   3rd Qu.:2014   3rd Qu.: 1512  
##                     Max.   :2640252   Max.   :2015   Max.   :27551  
##                                                      NA's   :6      
##  misdemeanortotal  totalcrimes    lawenforcementtotal   othermisd      
##  Min.   :   52    Min.   :   66   Min.   :   0.0      Min.   :   13.0  
##  1st Qu.:  890    1st Qu.: 1214   1st Qu.: 131.0      1st Qu.:  329.8  
##  Median : 1382    Median : 1934   Median : 224.0      Median :  510.5  
##  Mean   : 6448    Mean   : 8977   Mean   : 550.8      Mean   : 2167.5  
##  3rd Qu.: 4055    3rd Qu.: 5520   3rd Qu.: 471.0      3rd Qu.: 1512.2  
##  Max.   :68603    Max.   :95647   Max.   :4751.0      Max.   :21854.0  
##  NA's   :6        NA's   :6       NA's   :44          NA's   :6        
##     divorce         mentalchild      mentaladult      mentaltotal     
##  Min.   :    5.0   Min.   :   365   Min.   :  1082   Min.   :   2162  
##  1st Qu.:  158.8   1st Qu.:  2215   1st Qu.:  6798   1st Qu.:   9067  
##  Median :  255.5   Median :  6826   Median : 14301   Median :  21757  
##  Mean   :  931.6   Mean   : 30786   Mean   : 71983   Mean   : 102768  
##  3rd Qu.:  667.5   3rd Qu.: 23274   3rd Qu.: 44093   3rd Qu.:  67605  
##  Max.   :15666.0   Max.   :290846   Max.   :915912   Max.   :1159220  
##  NA's   :8         NA's   :12       NA's   :12       NA's   :12       
##   unemployment    violentfelony      drugfelony       dwifelony    
##  Min.   :  1000   Min.   :   1.0   Min.   :   0.0   Min.   :  1.0  
##  1st Qu.:  9275   1st Qu.:  57.0   1st Qu.:  43.0   1st Qu.: 34.0  
##  Median : 15450   Median : 103.5   Median :  85.5   Median : 56.0  
##  Mean   : 48497   Mean   : 698.0   Mean   : 483.0   Mean   : 96.7  
##  3rd Qu.: 37825   3rd Qu.: 370.0   3rd Qu.: 223.0   3rd Qu.:108.0  
##  Max.   :459700   Max.   :9507.0   Max.   :7679.0   Max.   :615.0  
##  NA's   :8        NA's   :6        NA's   :6        NA's   :6      
##   otherfelony         drugmisd           dwimisd        propertymisd  
##  Min.   :    7.0   Min.   :    3.00   Min.   :  19.0   Min.   :    9  
##  1st Qu.:  167.2   1st Qu.:   65.25   1st Qu.: 186.2   1st Qu.:  234  
##  Median :  279.5   Median :  122.50   Median : 306.5   Median :  480  
##  Mean   : 1251.5   Mean   : 1471.95   Mean   : 620.5   Mean   : 2188  
##  3rd Qu.:  824.2   3rd Qu.:  464.00   3rd Qu.: 713.5   3rd Qu.: 1212  
##  Max.   :13590.0   Max.   :27784.00   Max.   :3943.0   Max.   :30197  
##  NA's   :6         NA's   :6          NA's   :6        NA's   :6      
##  allageproverty   ueighteenpoverty fivetoseventenpov
##  Min.   :   490   Min.   :   123   Min.   :    80   
##  1st Qu.:  7311   1st Qu.:  2140   1st Qu.:  1423   
##  Median : 12260   Median :  3510   Median :  2343   
##  Mean   : 48550   Mean   : 15236   Mean   : 10440   
##  3rd Qu.: 32620   3rd Qu.:  9326   3rd Qu.:  6185   
##  Max.   :615319   Max.   :206071   Max.   :138276   
## 

We saw that some the columns in the summary are not in the right format and have some of NAs and are not in the right/accetable datatype for a perfect analysis. We will therefore convert them back to numerica datatype (a more acceptable datatype for analysis)

df[, -1] <- sapply(df[, -1], as.numeric) #Convert some df column to numeric

A glimpse of the dataset with the NAS

missmap(df, legend = TRUE, col = c("wheat","darkred", col=c('yellow', 'darkgreen')),
        main ="Plot Showing The Missing Values Per Observation",
        y.cex = 1, x.cex = 0.8,csvar = NULL, tsvar =NULL, rank.order = TRUE)

The above is a graphical representation of all the missing values.The darkgreen area shows the no missing values while the other color indicates the present of missing values in each columns.

A mulitivariate Imputation by Chained Equations (mice package) is said to be the best approach to tackle end-on the menace missing observation may cause. It is well know that if proper care is not taken in dealing with missing observation, the analysis would likely be biased. To avoid this, the mice package in R will use a fully conditional specification were each incomplete variable is imputed by a separate model at which the algoritm will impute mixes of continous, binary, ordered or unordered data based on the data type.

• Filling The Missing Values

df_missing_val <- mice(df, m=1, method='cart', seed = 500, printFlag=FALSE, where = is.na(df))

## Warning: Number of logged events: 87

df_complete <- complete(df_missing_val, action='long', include =FALSE)

If possible, we will show map of the location, and as googlevis requirement, the longtitude and latitude needed for the map must be in “lat:long” format

• Cleaning and adding map longitute and latititudes.

df_complete <- merge(df_complete , lat_long, all=T)
df_complete$Lat_Long = paste(df_complete$lat, df_complete$long, sep=":")
df_cleaned <- df_complete[, -c(2,3)]

df_pca <- df_complete[-c(235, 302),  c(1, 4,6, 7, 8, 9, 11, 12,13, 15, 23, 24, 25)]

DT::datatable(df_cleaned, options = list(
  pageLength = 10
))

Checking out for outliers (an observation point that is abnormal distant from other observations) in an analysis one of the best way to avoid variablity in the measurement or experimantal error.

outliers_check <- function(check) {
  check < quantile(check, 0.25, na.rm = TRUE) - 1.5 * IQR(check, na.rm = TRUE) |
    check > quantile(check, 0.75, na.rm = TRUE) + 1.5 * IQR(check, na.rm = TRUE)
}


outlier_obs <- df_cleaned %>%
  group_by(year, county) %>%
  mutate(Flagged = is.na(felonytotal) | outliers_check(felonytotal) | felonytotal <= 0) %>%
  ungroup()

#outliers <- subset(outlier_obs, Flagged == TRUE | year==2015, select=c(year, county, population, lawenforcementtotal, felonytotal, divorce))

A boxplot showing observation with oultiers

theme_set(theme_bw())


ct <- ggplot(df_cleaned, aes(county, totalcrimes))
ct + geom_boxplot() + 
  geom_dotplot(binaxis='y', 
               stackdir='center', 
               dotsize = .15, 
               fill="red") +
  theme(axis.text.x = element_text(angle=75, vjust=0.6)) + 
  labs(title="Box plot + Dot plot of", 
       subtitle="County Vs Total Crimes",
       x="New York State Counties",
       y="Total Crimes (Felony & Misdemeanor) ")

df_corr <- df_cleaned[ , -c(1,3, 24, 25, 26)]


corr_coef <- cor(df_corr, method = "kendall", use = "complete")

ggcorrplot(corr_coef, hc.order = TRUE, 
           type = "lower", 
           lab = TRUE, 
           lab_size = 3, 
           method="circle", 
           colors = c("tomato2", "white", "springgreen3"), 
           title="Correlogram of variables", 
           ggtheme=theme_bw)

The data set would be used in two different way.

df_incident <- sqldf::sqldf("SELECT 

county as County, year as Year, population as Populations, 
(((felonytotal + misdemeanortotal)/population)*100000) as TotalCrimes,
((mentalchild/population)*10000) as ChildMental, ((mentaladult/population)*10000) as AdultMental,
((mentaltotal/population)*10000) as TotalMental,
((lawenforcementtotal/population)*100000) as lawenforcementtotal,
((divorce/population)*100000) as DivorcePerPopPerK, ((drugfelony/population)*100000) as DrugFelonyPerPopPerK, 
((violentfelony/population)*100000) as ViolentTotal, ((dwifelony/population)*100000) as DWIFelony,
((otherfelony/population)*100000) as OtherFelony, ((felonytotal/population)*100000) as TotalFelony, 
((propertymisd/population)*100000) as PropertyMisdemeanor, 
((drugmisd/population)*100000) as DrugMisdemeanor, ((dwimisd/population)*100000) as DwiMisdemeanor,
((othermisd/population)*100000) as OtherMisdemeanor,
((misdemeanortotal)/population)*100000
as TotalMisdemeanor, 
((unemployment/population)*100000) as UnemploymentTotal,
((allageproverty/population)*100000) as AllPersonInPoverty ,
((ueighteenpoverty/population)*100000) as Under18InPoverty,
((fivetoseventenpov/population)*100000) as FiveTo17InPoverty,lat, long, Lat_long
FROM df_cleaned
WHERE county NOT IN ('Out-of-State', 'St. Lawrence') 
ORDER BY Year, county ASC

")
df_incident <- as.data.frame(df_incident)
kable(head(df_incident[, 1:9]))

kable(head(df_incident[, 10:18]))

kable(head(df_incident[, 19:26]))

Charts

theme_set(theme_bw())

poverty_inc <- sqldf::sqldf("SELECT 
                        Year, AllPersonInPoverty, County
                        FROM df_incident
                        ORDER BY AllPersonInPoverty DESC LIMIT 30
                        ")

poverty_plot <- ggplot(poverty_inc, aes(x=AllPersonInPoverty, y=reorder(County, +AllPersonInPoverty), fill=as.factor(Year))) +
geom_point(colour="blue", size=5, alpha=.8) +
scale_fill_brewer(palette="Blues", breaks=rev(levels(poverty_inc$Year))) +
  labs(title="Chart of Top 30 Poverty  In New York State Per County
       \n (Incident Rate Per 100K Persons From 2010 To 2015))")




## Divorce Chart

div_inc <- sqldf::sqldf("SELECT 
                        Year, DivorcePerPopPerK, County
                        FROM df_incident
                        ORDER BY DivorcePerPopPerK DESC LIMIT 30
                        ")

divorce_plot <- ggplot(div_inc, aes(x=DivorcePerPopPerK, y=reorder(County, +DivorcePerPopPerK), fill=as.factor(Year))) +
geom_point(colour="pink", size=5, alpha=.8) +
scale_fill_brewer(palette="Blues", breaks=rev(levels(div_inc$Year))) +
  labs(title="Chart of Top 30 Divorce  In New York State Per County
       \n (Incident Rate Per 100K Persons From 2010 To 2015))")



## Unemployment Chart



unemp_inc <- sqldf::sqldf("SELECT 
                        Year, UnemploymentTotal, County
                        FROM df_incident
                        ORDER BY UnemploymentTotal DESC LIMIT 30
                        ")


unemployment_plot <- ggplot(unemp_inc, aes(x=UnemploymentTotal, y=reorder(County, +UnemploymentTotal), fill=as.factor(Year))) +
geom_point(colour="brown", size=5, alpha=.8) +
scale_fill_brewer(palette="Blues", breaks=rev(levels(unemp_inc$Year))) +
  labs(title="Chart of Top 30 Unemployment In New York State Per County
       \n (Incident Rate Per 100K Persons From 2010 To 2015)")



## Mental Health Chart

mentalhealth_inc <- sqldf::sqldf("SELECT 
                        Year, TotalMental, County
                        FROM df_incident
                        ORDER BY TotalMental DESC LIMIT 30
                        ")

mental_plot <- ggplot(mentalhealth_inc, aes(x=TotalMental, y=reorder(County, +TotalMental), fill=as.factor(Year))) +
geom_point(colour="purple", size=5, alpha=.8, shape = 21) +
scale_fill_brewer(palette="Blues", breaks=rev(levels(mentalhealth_inc$Year))) +
  labs(title="Chart of Top 30 Mental Health In New York State Per County
       \n (Incident Rate Per 100K Persons From 2010 To 2015)")


poverty_plot #poverty plot
divorce_plot #divorce plot
unemployment_plot #unemployment plot
mental_plot #mental health plot

ggplot(df_incident, aes(TotalFelony)) + geom_density(aes(color=as.factor(Year), fill=as.factor(Year)),alpha=0.1) +
  labs(x="Total Felony", y="Density", title="Felony Total Per 100k Population Density") +
  theme(legend.position = "top")+ theme_economist()

misd_dens_plt <- ggplot(df_incident, aes(TotalMisdemeanor)) + geom_density(aes(color=as.factor(Year),
                        fill=as.factor(Year)),alpha=0.1) + labs(x="Misdemeanor Total", y="Density",
                        title="Misdemeanor Total Per 100k Population Density") + theme(legend.position = "top")+
                       theme_economist()

misd_dens_plt

felony <- dplyr::select(df_incident, TotalFelony, Year, UnemploymentTotal) %>%
    mutate(Year = as.factor(Year)) %>%
ggplot(aes(x = Year, y = TotalFelony)) + 
    geom_jitter(aes(col = Year), alpha = 0.8) +
    labs(title = "x") +
    theme_light()

misd <- dplyr::select(df_incident, TotalMisdemeanor, Year, UnemploymentTotal) %>%
    mutate(Year = as.factor(Year)) %>%
ggplot(aes(x = Year, y = TotalMisdemeanor)) + 
    geom_jitter(aes(col = Year), alpha = 0.8) +
    labs(title = "x") +
    theme_light()

p <- plot_grid(felony, misd) 
title <- ggdraw() + draw_label("A Plot of Felony and Misdemeanor Incident Per Population 100k", fontface='bold')
plot_grid(title, p, ncol=1, rel_heights = c(0.1, 2.2))

df_incident$size <- df_incident$TotalCrimes
colors <- c('#4AC6B7', '#1972A4', '#965F8A', '#FF7070', '#C61951')

comparison <- plot_ly(df_incident[1:1000, ], x = ~TotalMisdemeanor, y = ~TotalFelony, z = ~lawenforcementtotal,
                      color = ~Year, size = ~size, colors = colors,
             marker = list(symbol = 'circle', sizemode = 'diameter'), sizes = c(1, 50),
             text = ~paste('County:', County, '<br>Population Total:', Populations,
                           '<br>Total Unemployment Pop in 100K:', UnemploymentTotal)) %>%
  layout(title = 'Comparison Among Police, Felony and Misdemeanor Crimes In New York State',
         scene = list(xaxis = list(title = 'Misdemeanor Pop At Risk Per 1K',
                      gridcolor = 'rgb(255, 255, 255)',
                      type = 'log',
                      zerolinewidth = 1,
                      ticklen = 5,
                      gridwidth = 2),
               yaxis = list(title = ' Felony Crime Pop At Risk Per 100K',
                      gridcolor = 'rgb(255, 255, 255)',
                      zerolinewidth = 1,
                      ticklen = 5,
                      gridwith = 2),
               zaxis = list(title = 'Police Pop Per 100K',
                            gridcolor = 'rgb(255, 255, 255)',
                            type = 'log',
                            zerolinewidth = 1,
                            ticklen = 5,
                            gridwith = 2)),
         paper_bgcolor = 'rgb(243, 243, 243)',
         plot_bgcolor = 'rgb(243, 243, 243)')

comparison

qqnorm(df_cleaned$totalcrimes)
qqline(df$totalcrimes, col = 2)

counties <- (c("Albany", "Allegany", "Bronx", "Broome", "Cattaraugus", "Cayuga",
         "Chautauqua", "Chemung", "Chenango", "Clinton", "Columbia", "Cortland",
         "Delaware", "Dutchess", "Erie", "Essex", "Franklin","Fulton", "Genesee", 
         "Greene", "Hamilton", "Herkimer", "Jefferson", "Kings", "Lewis","Livingston",
         " Madison", "Monroe", "Montgomery", "Nassau", "New York", "Niagara", "Oneida",
         " Onondaga", " Ontario", "Orange", "Orleans", "Oswego", "Otsego",  "Putnam", "Queens",
         "Rensselaer", "Richmond", "Rockland","Saratoga", "Schenectady", "Schoharie", "Schuyler",
         "Seneca", "St Lawrence",  "Steuben", "Suffolk",  "Sullivan", "Tioga", "Tompkins", "Ulster",
         "Warren", "Washington", "Wayne", "Westchester", "Wyoming", "Yates"))


lab <- c(1, 2, 3, 4, 5,  6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
         26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49,
         52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64)
dummy <- cbind(counties, lab)
kable(head(dummy))

df_cleaned <- df_cleaned[-c(235, 302), -c(4, 5, 8, 12, 14:20,24:26)]


df_cleaned$county <- factor(df_cleaned$county,
                    levels=counties,
                    labels= lab)

df_cleaned$year <- factor(df_cleaned$year,
                     levels = c("2010", "2011", "2012", "2013", "2014", "2015"),
                     labels = c(1, 2, 3, 4, 5, 6));

#df_cleaned$totalcrimes <- as.factor(df_cleaned$totalcrimes)


set.seed(123)

splits <- sample.split(df_cleaned$county, SplitRatio = 0.8)
train_set <- subset(df_cleaned, splits==TRUE)
test_set <- subset(df_cleaned, splits==FALSE)

#train_set[, 2:19] <- scale(train_set[, 2:19])
#test_set[, 2:19] <- scale(test_set[, 2:19])

Regression

regressor <-  glm(totalcrimes~., family=gaussian(link = "identity"),data=train_set)
stepwise <- step(regressor, direction = "both")

## Start:  AIC=5002.1
## totalcrimes ~ county + population + year + lawenforcementtotal + 
##     divorce + mentalchild + mentaladult + unemployment + allageproverty + 
##     ueighteenpoverty + fivetoseventenpov
## 
##                       Df   Deviance    AIC
## - lawenforcementtotal  1  241942427 5000.3
## - divorce              1  242655053 5001.1
## <none>                    241814047 5002.1
## - allageproverty       1  245092297 5004.1
## - ueighteenpoverty     1  248825615 5008.5
## - year                 5  257390382 5010.5
## - mentalchild          1  253561895 5014.1
## - unemployment         1  253903448 5014.5
## - fivetoseventenpov    1  254127113 5014.8
## - mentaladult          1  256384468 5017.4
## - population           1  320388196 5083.1
## - county              58 1098802081 5332.7
## 
## Step:  AIC=5000.26
## totalcrimes ~ county + population + year + divorce + mentalchild + 
##     mentaladult + unemployment + allageproverty + ueighteenpoverty + 
##     fivetoseventenpov
## 
##                       Df   Deviance    AIC
## - divorce              1  242705306 4999.2
## <none>                    241942427 5000.3
## + lawenforcementtotal  1  241814047 5002.1
## - allageproverty       1  245277509 5002.3
## - ueighteenpoverty     1  248945811 5006.7
## - year                 5  258021041 5009.2
## - mentalchild          1  253617214 5012.2
## - unemployment         1  254031095 5012.6
## - fivetoseventenpov    1  254342989 5013.0
## - mentaladult          1  256596366 5015.6
## - population           1  320413512 5081.1
## - county              58 1107065647 5332.9
## 
## Step:  AIC=4999.19
## totalcrimes ~ county + population + year + mentalchild + mentaladult + 
##     unemployment + allageproverty + ueighteenpoverty + fivetoseventenpov
## 
##                       Df   Deviance    AIC
## <none>                    242705306 4999.2
## + divorce              1  241942427 5000.3
## + lawenforcementtotal  1  242655053 5001.1
## - allageproverty       1  248206467 5003.8
## - ueighteenpoverty     1  250067204 5006.0
## - year                 5  259171855 5008.6
## - unemployment         1  254031855 5010.6
## - mentalchild          1  254634736 5011.3
## - fivetoseventenpov    1  254805543 5011.5
## - mentaladult          1  256868278 5013.9
## - population           1  327927630 5086.0
## - county              58 2175493407 5530.2

Analysis of Variance (ANOVA)

Analysis of Variance (ANOVA) above is used to compare the means of the variance for statistical significance.

anova(stepwise, test="F")

# 
fit <- aov(totalcrimes~., data=train_set)
summary(fit)

##                      Df    Sum Sq   Mean Sq  F value   Pr(>F)    
## county               58 1.080e+11 1.863e+09 1710.044  < 2e-16 ***
## population            1 4.378e+08 4.378e+08  401.888  < 2e-16 ***
## year                  5 2.327e+07 4.654e+06    4.272 0.000993 ***
## lawenforcementtotal   1 2.568e+06 2.568e+06    2.358 0.126070    
## divorce               1 2.233e+06 2.233e+06    2.050 0.153654    
## mentalchild           1 1.845e+07 1.845e+07   16.940 5.44e-05 ***
## mentaladult           1 3.540e+07 3.540e+07   32.500 3.77e-08 ***
## unemployment          1 1.471e+07 1.471e+07   13.507 0.000298 ***
## allageproverty        1 6.444e+06 6.444e+06    5.916 0.015793 *  
## ueighteenpoverty      1 1.628e+05 1.628e+05    0.149 0.699407    
## fivetoseventenpov     1 1.231e+07 1.231e+07   11.304 0.000911 ***
## Residuals           222 2.418e+08 1.089e+06                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# HSD.test(fit, 'unemployment ')

Residual Plot

residual <- data.frame(Fitted = fitted(stepwise), Residuals = resid(stepwise), Treatment = train_set$totalcrimes)

plot_res <- ggplot(residual, aes(Fitted, Residuals, colour = Treatment)) + geom_point()

ggplotly(plot_res)

## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`

Confident Interval

confint(stepwise)

## Waiting for profiling to be done...

##                           2.5 %        97.5 %
## (Intercept)        3.514137e+04  5.257861e+04
## county2           -4.509359e+04 -3.026503e+04
## county3            1.759665e+05  2.561847e+05
## county4           -1.847632e+04 -1.218690e+04
## county5           -4.003689e+04 -2.698067e+04
## county6           -3.933880e+04 -2.645562e+04
## county7           -2.996834e+04 -2.008484e+04
## county8           -3.761863e+04 -2.524088e+04
## county9           -4.475257e+04 -3.004656e+04
## county10          -3.902694e+04 -2.619976e+04
## county11          -4.312707e+04 -2.893036e+04
## county12          -4.447702e+04 -2.975810e+04
## county13          -4.543456e+04 -3.042556e+04
## county14          -6.699662e+03 -3.683733e+03
## county15           6.884225e+04  1.050451e+05
## county16          -4.706535e+04 -3.154038e+04
## county17          -4.421168e+04 -2.960197e+04
## county18          -4.386743e+04 -2.948005e+04
## county19          -4.344640e+04 -2.916084e+04
## county20          -4.478454e+04 -2.993710e+04
## county21          -5.172073e+04 -3.468023e+04
## county22          -4.270715e+04 -2.869432e+04
## county23          -3.289370e+04 -2.202968e+04
## county24           2.831810e+05  4.268649e+05
## county25          -4.877178e+04 -3.274365e+04
## county26          -4.276384e+04 -2.872511e+04
## county28           4.858833e+04  7.553175e+04
## county29          -4.467252e+04 -3.006745e+04
## county30           9.342428e+04  1.472886e+05
## county31           1.852829e+05  2.661927e+05
## county32          -1.728764e+04 -1.161247e+04
## county33          -1.381400e+04 -9.091988e+03
## county36           4.654261e+03  1.032043e+04
## county37          -4.618999e+04 -3.096345e+04
## county38          -3.264768e+04 -2.201878e+04
## county39          -4.273305e+04 -2.860299e+04
## county40          -3.802003e+04 -2.554821e+04
## county41           2.211987e+05  3.367049e+05
## county42          -2.767832e+04 -1.873803e+04
## county43           1.653561e+04  2.676194e+04
## county44          -6.870157e+03 -2.116914e+03
## county45          -1.868452e+04 -1.235097e+04
## county46          -2.706671e+04 -1.781711e+04
## county47          -4.634840e+04 -3.126450e+04
## county48          -4.977473e+04 -3.336613e+04
## county49          -4.742532e+04 -3.170056e+04
## county52          -3.630333e+04 -2.444980e+04
## county53          -3.631639e+04 -2.459065e+04
## county54           1.126232e+05  1.757783e+05
## county55          -3.973076e+04 -2.648663e+04
## county56          -4.529073e+04 -3.044249e+04
## county57          -3.595419e+04 -2.424635e+04
## county58          -2.260695e+04 -1.485317e+04
## county59          -4.202072e+04 -2.810953e+04
## county60          -4.285657e+04 -2.874644e+04
## county61          -3.851339e+04 -2.601193e+04
## county62           6.553923e+04  1.023905e+05
## county63          -4.659207e+04 -3.131903e+04
## county64          -4.876006e+04 -3.267136e+04
## population        -1.316622e-01 -8.400117e-02
## year2             -2.078425e+01  8.704594e+02
## year3             -2.181291e+01  8.911599e+02
## year4              1.632421e+02  1.072794e+03
## year5              1.180490e+02  1.037611e+03
## year6             -5.025215e+02  4.932371e+02
## mentalchild       -6.769890e-02 -1.742043e-02
## mentaladult        8.618134e-03  2.902459e-02
## unemployment       1.103547e-02  4.501041e-02
## allageproverty    -1.598438e-01 -1.112728e-02
## ueighteenpoverty   8.897740e-02  6.283368e-01
## fivetoseventenpov -7.998222e-01 -2.084640e-01

Prediction

pred <- predict(regressor, test_set, type="response")
actual_pred <- cbind.data.frame(test_set$totalcrimes, pred)
actual_pred <- actual_pred[1:20, ]
kable(actual_pred)

Check for Accuracy

MAPE The mean absolute percentage error (MAPE), also known as mean absolute percentage deviation (MAPD), is a measure of prediction accuracy of a forecasting method in statistics. Represented mathematically as:

\(M=\quad \frac { 100 }{ n } \sum _{ t=1 }^{ n }{ \left| \frac { { A }_{ t }\quad -\quad { F }_{ t } }{ { A }_{ t } } \right| }\) ,

Where,

\({ A }_{ t }\quad =\quad Actual\quad Value\\ \quad \quad \quad \quad \quad \quad { F }_{ t }\quad =\quad Predicted\quad Value.\)

min_max_accuracy <- mean(apply(actual_pred, 1, min) / apply(actual_pred, 1, max))*100 
print(paste('The Percentage Accuracy Is: ', min_max_accuracy))

## [1] "The Percentage Accuracy Is:  87.4774508797211"

mape <- mean(abs((actual_pred$pred - actual_pred$`test_set$totalcrimes`))/actual_pred$pred)*100 
print(paste('The MAPE Is: ', mape ))

## [1] "The MAPE Is:  16.5027632643134"

Principal Component Analysis

To obtain the best model, the best way is to adopt a dimensionality reduction. The PCA, like factor analysis, is a mathematical procedure that transforms a number of correlated variables into a smaller number of uncorrelated variables.

prin_comp <- prcomp(train_set[, -c(1,3)], scale. = T)
biplot(prin_comp, scale = 0)

std_dev <- prin_comp$sdev
pr_var <- std_dev^2
pca_var <- pr_var/sum(pr_var)

plot(pca_var, xlab = "Principal Component",
             ylab = "Proportion of Variance Explained",
             type = "b")

plot(cumsum(pca_var), xlab = "Principal Component",
              ylab = "Cumulative Proportion of Variance Explained",
              type = "b")

pca2 <- df_pca[, -1]
pca1 <- prcomp(pca2, scale. = TRUE)

cir <- function(center = c(0, 0), npoints = 100) {
    r = 1
    t = seq(0, 2 * pi, length = npoints)
    x = center[1] + r * cos(t)
    y = center[1] + r * sin(t)
    return(data.frame(x = x, y = y))
}
corcir = cir(c(0, 0), npoints = 100)

correlations = as.data.frame(cor(pca2, pca1$x))

arrows = data.frame(x1 = c(0, 0, 0, 0), y1 = c(0, 0, 0, 0), x2 = correlations$PC1, 
    y2 = correlations$PC2)


ggplot() + geom_path(data = corcir, aes(x = x, y = y), colour = "black") + 
    geom_segment(data = arrows, aes(x = x1, y = y1, xend = x2, yend = y2), colour = "blue") + 
    geom_text(data = correlations, aes(x = PC1, y = PC2, label = rownames(correlations))) + 
    geom_hline(yintercept = 0, colour = "brown") + geom_vline(xintercept = 0, 
    colour = "brown") + xlim(-1.1, 1.1) + ylim(-1.1, 1.1) + labs(x = "pc1 aixs", 
    y = "pc2 axis") + ggtitle("Circle of correlations")

PCA Visual

#train set
train_pca <- data.frame(TotalCrimes = train_set$totalcrimes, prin_comp$x)
train_pca <- train_pca[,1:8]
pca_model <- rpart(TotalCrimes ~ .,data = train_pca, method = "anova")

#test set
test_pca <- predict(prin_comp, newdata = test_set)
test_pca <- as.data.frame(test_pca)
test_pca <- test_pca[,1:8]
#prediction
prediction <- predict(pca_model, test_pca)

INTERPRETATIONS:

• Generalized Linear Model: The analysis revealed that variables like: years, mentalchild, mentaladult, unemployment, allageproverty, ueighteenpoverty, fivetoseventenpov and population have P-values less than 0.05, which means they are statistically significant.

• Residual Plot: Unpattern nature of the residuals conform to the accuracy of model.

• Analysis of Variance: Supports the GLM model but struck out the variable “ueighteenpoverty”.

• Confidence Interval: The CI test indicates that at 95% the confidence interval would contain the true population mean as none of the interval is zero.

• Principal Component Analysis: After the dimensionality reduction, PCA also confirms the best model.

DECISION: Since the calculated P-values are less than the tabulated values (\(??\) = 0.05), we therefore reject the alternative hypothesis, H_o.

CONCLUSION: We therefore conclude that there is a positive relationship between crimes and factors like Economic and Psychological as confirmed by variables like: Unemployment, Mental Health, Poverty, County (or location) and Year. It also means; as the factors increases, so also is crime rate and vice-versa. We can therefore confirm that not only increase and decrease in population affect the crime rate, but other factors like economic and psychological also contribute significantly to crime rate in a geographical area. While security and social factors may not contribute significantly to crime rate, but a further analysis may prove otherwise.

Bibliography

The Origins of attachment theory- John Bowlby and Mary Ainsworth

http://www.sakkyndig.com/psykologi/artikler/Bowlbyatt.pdf

ftp://statgen.ncsu.edu/pub/thorne/molevoclass/AtchleyOct19.pdf

https://data.cityofnewyork.us/Public-Safety/NYPD-Complaint-Data-Current-YTD/5uac-w243

http://data.beta.nyc/dataset/pediacities-nyc-neighborhoods/resource/35dd04fb-81b3-479b-a074-a27a37888ce7

http://www.kevjohnson.org/making-maps-in-r/

http://leaflet-extras.github.io/leaflet-providers/preview/

http://nyscommunityaction.org/wp-content/uploads/2018/04/Comparitive-Look-@-County-Profiles.pdf

https://stackoverflow.com/questions/19449307/how-to-solve-error-in-as-posixct-numericx2l-origin-must-be-supp

http://r-statistics.co/Top50-Ggplot2-Visualizations-MasterList-R-Code.html

https://data.ny.gov/en/browse?q=nypd

http://www1.nyc.gov/site/nypd/stats/crime-statistics/crime-statistics-landing.page

https://data.ny.gov/en/browse?q=New%20York%20State%20Population%20Data%3A%20Beginning%202003&anonymous=true&sortBy=relevance

https://www.health.ny.gov/statistics/vital_statistics/2006/table50.htm

http://www.postgresqltutorial.com/postgresql-alter-table/

https://www.statmethods.net/graphs/scatterplot.html

http://nyscommunityaction.org/wp-content/uploads/2018/04/Comparitive-Look-@-County-Profiles.pdf