Chicago Kidnapping Crime: Go Low or High?
Aathirah Valmay
2023-08-18
OBJECTIVES
Tujuan dari latihan ini adalah membuat model time series dari data Chicago Crime dataset (https://data.cityofchicago.org/Public-Safety/Crimes-2001-to-present/ijzp-q8t2). Kita akan melakukan beberapa tahap dengan memahami data terlebih dahulu, kemudian akan kita buat pola berdasarkan waktu time series dengan beberapa persyaratan yang harus dipatuhi. Kemudian akan kita analisis pola data time seriesnya apakah mengandung tren dan seasonal, yang mana akan mempengaruhi model apa yang akan dipilih untuk pemodelan time series kali ini. Harapan dari latihan kali ini adalah mendapatkan hasil error dari model sekecil mungkin. Berikut library yang akan digunakan untuk latihan kali ini.
library(dplyr) # data wrangling
library(lubridate) # date manipulation
library(padr) # complete data frame
library(zoo) # Missing value imputation
library(forecast) # time series library
library(TTR) # for Simple moving average function
library(MLmetrics) # calculate error
library(tseries) # adf.test
library(fpp) # data for forecasting: principles and practice
library(TSstudio) # mempercantik visualisasi timeseries
library(ggplot2)
library(tidyr)1. Read Data & Data Understanading
Kita akan membaca data csv dan menyimpannya menjadi variabel
crimes
crimes <- read.csv("Crimes_2001_to_Present.csv")head(crimes)Deskripsi kolom
Berikut adalah keterangan dari tiap kolom yang ada pada data
crimes
ID : Unique identifier for the record
Case Number : The Chicago Police Department RD Number
(Records Division Number), which is unique to the incident
Date : Date when the incident occurred. this is sometimes a
best estimate Block : The partially redacted address where
the incident occurred, placing it on the same block as the actual
address IUCR : The Illinois Uniform Crime Reporting code.
This is directly linked to the Primary Type and Description
Primary Type : The primary description of the IUCR code
Description : The secondary description of the IUCR code, a
subcategory of the primary description Location Description
: Description of the location where the incident occurred
Arrest : Indicates whether an arrest was made
Domestic : Indicates whether the incident was
domestic-related as defined by the Illinois Domestic Violence Act
Beat : Indicates the beat where the incident occurred. A
beat is the smallest police geographic area – each beat has a dedicated
police beat car. Three to five beats make up a police sector, and three
sectors make up a police district. The Chicago Police Department has 22
police districts District : Indicates the police district
where the incident occurred Ward : The ward (City Council
district) where the incident occurred Community Area :
Indicates the community area where the incident occurred. Chicago has 77
community areas FBI Code : Indicates the crime
classification as outlined in the FBI’s National Incident-Based
Reporting System (NIBRS) X Coordinate : The x coordinate of
the location where the incident occurred in State Plane Illinois East
NAD 1983 projection. This location is shifted from the actual location
for partial redaction but falls on the same block
Y Coordinate : The y coordinate of the location where the
incident occurred in State Plane Illinois East NAD 1983 projection. This
location is shifted from the actual location for partial redaction but
falls on the same block Year : Year the incident occurred
Updated On : Date and time the record was last updated
Latitude : The latitude of the location where the incident
occurred. This location is shifted from the actual location for partial
redaction but falls on the same block Longitude : The
longitude of the location where the incident occurred. This location is
shifted from the actual location for partial redaction but falls on the
same block Location : The location where the incident
occurred in a format that allows for creation of maps and other
geographic operations on this data portal. This location is shifted from
the actual location for partial redaction but falls on the same
block
Mengecek Uniqueness crimes dari kolom Primary Type untuk
keperluan memilih tipe crime yang akan kita analisis
table(crimes$Primary.Type)#>
#> ARSON
#> 5 10810
#> ASSAULT BATTERY
#> 403830 1106867
#> BURGLARY CONCEALED CARRY LICENSE VIOLATION
#> 301058 1118
#> CRIM SEXUAL ASSAULT CRIMINAL DAMAGE
#> 20283 694847
#> CRIMINAL SEXUAL ASSAULT CRIMINAL TRESPASS
#> 7082 167306
#> DECEPTIVE PRACTICE DOMESTIC VIOLENCE
#> 275658 1
#> GAMBLING HOMICIDE
#> 10643 12463
#> HUMAN TRAFFICKING INTERFERENCE WITH PUBLIC OFFICER
#> 97 12279
#> INTIMIDATION KIDNAPPING
#> 3835 5796
#> LIQUOR LAW VIOLATION MOTOR VEHICLE THEFT
#> 11785 296606
#> NARCOTICS NON-CRIMINAL
#> 542685 155
#> NON-CRIMINAL (SUBJECT SPECIFIED) NON - CRIMINAL
#> 6 23
#> OBSCENITY OFFENSE INVOLVING CHILDREN
#> 663 43278
#> OTHER NARCOTIC VIOLATION OTHER OFFENSE
#> 112 378078
#> PROSTITUTION PUBLIC INDECENCY
#> 57235 129
#> PUBLIC PEACE VIOLATION RITUALISM
#> 34840 24
#> ROBBERY SEX OFFENSE
#> 221818 25308
#> STALKING THEFT
#> 3972 1248456
#> WEAPONS VIOLATION
#> 88583Data Aggregation
Mari kita memilih satu jenis crime saja, kali ini kita akan memilih
jenis crime KIDNAPPING
crime_kd <- crimes %>%
filter(Primary.Type == "KIDNAPPING")head(crime_kd)Karena kita telah memilih jenis crime KIDNAPPING, maka
kita memerlukan data agregasi berdasarkan waktu. Sebelumnya kita
perhatikan dulu frekuensi tiap tahun dari jenis crime kidnapping
berikut
table(crime_kd$Year)#>
#> 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016
#> 933 829 705 482 389 338 328 358 121 2 3 2 3 2 156 202
#> 2017 2018 2019 2020 2021 2022 2023
#> 190 172 173 120 87 117 84Berdasarkan data di atas, tidak menunjukkan data yang signifikan ditunjukkan dalam frekuensi harian. Maka dari itu, kita akan mengagregasi data tersebut berdasarkan bulan. Pertama-tama kita extract terlebih dahulu bulan dan tahun dari tanggal kejadian
crime_kd <- crime_kd %>%
mutate(Date = mdy_hms(Date))
crime_kdcrime_kd$Month <- format(as.Date(crime_kd$Date), "%Y-%m")crime_kd <- crime_kd %>%
mutate(Month = ym(Month))
crime_kdcrime_agg <- crime_kd %>%
group_by(Month) %>%
summarise(kidnap = n())
crime_agg2. Data Wrangling
Kita perlu memenuhi persyaratan pada file untuk melakukan time series sebagai berikut
1. Data time series wajib urut
Kita akan mengurutkan terlebih dahulu data crime_agg hasil agregasi
data crime_kd
crime_agg <- crime_agg %>% arrange(Month)Mengecek apakah terdapat missing values
crime_agg$Year %>% is.na() %>% sum()#> [1] 02. Tidak boleh ada waktu yang terlewat
Kita akan menggunakan seq.Date untuk mengecek apakah ada
waktu yang terlewat berdasarkan bulan, maka kita akan menggunakan
by = "month"
all(seq.Date(from = min(crime_agg$Month), to = max(crime_agg$Month), by = "month")==crime_agg$Month)#> [1] FALSEHasil di atas menunjukkan terdapat waktu yang terlewat, maka dari itu
akan kita pad
crime_agg <- crime_agg %>% pad()
crime_aggKita cek ulang apakah padnya berhasil
all(seq.Date(from = min(crime_agg$Month), to = max(crime_agg$Month), by = "month")==crime_agg$Month)#> [1] TRUEHasil sudah menunjukkan TRUE, maka data sudah terurut dan tidak ada yang terlewat
3. No Missing Value
Kita akan mengecek data yang hilang menggunakan
anyNA
anyNA(crime_agg)#> [1] TRUESetelah padding, tentunya pasti ada data yang hilang. Karena memang tidak ada crime pada waktu-waktu tersebut, maka akan kita isi menggunakan nilai nol.
crime_agg <- crime_agg %>%
mutate(kidnap = na.fill(kidnap,fill=0))
crime_agganyNA(crime_agg)#> [1] FALSESetelah kita cek ulang, seluruh data telah terisi dan kita siap untuk
memproses data crime_agg menjadi data time series
3. Data Pre-processing
Untuk memodelkan time series di R, kita perlu mengubah objek menjadi bentuk time series atau ts()
Kita perlu menyiapkan parameter data, yaitu data
crime_agg, start, yaitu awal mula waktu yaitu tahun 2001
bulan Januari, dan frequency, yaitu jumlah frekuensi dari
pola yang telah kita buat, karena data kita adalah bulanan, sementara
kita akan memprediksi secara tahunan, maka frekuensi yang kita buat
adalah 12
# your code here
crime_ts <- ts(data = crime_agg$kidnap,
start = c(2001,1),
frequency = 12
)
crime_ts#> Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
#> 2001 75 57 72 49 56 119 107 88 91 59 74 86
#> 2002 72 57 73 73 105 75 55 78 67 63 59 52
#> 2003 47 37 58 57 84 68 57 56 71 48 56 66
#> 2004 53 53 50 45 38 33 35 23 31 39 35 47
#> 2005 26 23 33 24 41 34 35 41 28 30 35 39
#> 2006 32 30 33 31 36 34 30 18 22 31 21 20
#> 2007 31 21 35 30 28 27 24 23 30 21 32 26
#> 2008 30 33 28 33 35 29 34 20 30 31 24 31
#> 2009 25 26 18 20 21 8 0 2 0 0 0 1
#> 2010 0 1 0 0 0 0 1 0 0 0 0 0
#> 2011 0 1 0 1 0 0 0 1 0 0 0 0
#> 2012 0 0 0 0 0 0 1 0 1 0 0 0
#> 2013 0 0 0 0 0 1 1 0 0 0 0 1
#> 2014 1 0 0 1 0 0 0 0 0 0 0 0
#> 2015 0 1 13 14 14 8 15 21 17 16 18 19
#> 2016 23 10 15 19 14 26 15 22 17 15 15 11
#> 2017 8 13 18 15 19 16 23 13 22 11 16 16
#> 2018 16 14 14 12 17 18 14 15 12 7 16 17
#> 2019 8 14 19 13 21 13 21 10 13 4 19 18
#> 2020 14 11 8 3 11 15 8 11 14 12 7 6
#> 2021 3 6 10 8 9 8 5 9 7 6 9 7
#> 2022 7 5 7 5 11 10 15 10 14 11 12 10
#> 2023 12 5 6 5 22 13 19 2Mengecek kelas dari crime_ts
class(crime_ts)#> [1] "ts"4. Exploratory Data Analysis
Visualisasi
Kita akan membuat plot terhadap data time series yang telah kita buat
crime_ts %>% autoplot()
Insight dari grafik di atas adalah: - Terdapat penurunan drastis dari rentang tahun 2001 sampai 2010 - Terdapat nilai kriminal kidnapping yang sangat rendah dari rentang tahun 2010 sampai 2015 - Terdapat jumlah kriminal kidnapping yang cenderung tetap dari tahun 2015 sampai dengan 2020
Kita akan melakukan subsetting dari tahun yang paling terakhir, yaitu 2015 sampai dengan 2023.
window(crime_ts, start = c(2015, 1), end = c(2023,8)) %>%
autoplot()
berdasarkan data di atas, kita bisa lebih mudah untuk membuat model prediksi dikarenakan tidak ada data yang janggal atau seasonal pada urutan seluruh data
Subset data crime_ts:
crime_ts <- window(crime_ts, start = c(2015, 1), end = c(2023,8))Decomposition
Kita akan melakukan dekomposisi pada data untuk melihat apakah data tersebut memiliki tren atau pun seasonal dengan plot yang lebih jelas
crime_decom <- decompose(crime_ts)
crime_decom %>% autoplot()
Terlihat pada data di atas, row data menunjukkan plot
asli data. Dari data tersebut kita dapat mengetahui bahwa plot tersebut
bersifat additive. Kemudian untuk row trend masih belum
menunjukkan tren yang baik. pada row seasonal menunjukkan
adanya pola berulang dari data tersebut. Kemudian row
remainder adalah pola data yang tak tertangkap oleh data
trend dan seasonal.
Maka dari itu, kita akan subset ulang menjadi lebih singkat lagi yaitu dari tahun 2020 sampai dengan tahun 2023
crime_ts <- window(crime_ts, start = c(2020, 1), end = c(2023,8))Berikut kita lihat lagi dekomposisi dari hasil subset terbaru
crime_decom <- decompose(crime_ts)
crime_decom %>% autoplot()
Kita akan menggunakan hasil subset terbaru di atas, dikarenakan row
trend menunjukkan hasil yang lebih baik dari sebelumnya,
yaitu menunjukkan adanya tren naik dan turun. Sekaligus terdapat data
seasonal yang menunjukkan pola berulang pada data
5. Cross Validation
Pada tahap ini kita akan membagi data menjadi data train dan data tes. Kita akan memilih data test adalah data kriminal kidnapping dari tahun 2023 yang berjumlah 8 bulan saja. Maka dari itu kita akan mengambil data train adalah data selain tahun 2023.
#membuat data test
crime_test <- tail(crime_ts,8)
#membuat data train
crime_train <- head(crime_ts,-8)crime_test#> Jan Feb Mar Apr May Jun Jul Aug
#> 2023 12 5 6 5 22 13 19 2crime_train#> Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
#> 2020 14 11 8 3 11 15 8 11 14 12 7 6
#> 2021 3 6 10 8 9 8 5 9 7 6 9 7
#> 2022 7 5 7 5 11 10 15 10 14 11 12 106. Build Model
Kita akan membuat model Exponential Smoothing dan Mole HoltWinters dengan pertimbangan model memiliki tren dan seasonal.
Exponential Smoothing
- Membuat Model
Kita akan membuat model menggunakan fungsi ets() dengan
parameter AZZ yang tandanya plot timeseries yang kita
miliki additive dan memiliki tren dan seasonal
crime_ets <- ets(y = crime_train,
model = "AZZ")tuning model dengan alpha = 0.1 sebagai pembanding batas minimum
crime_ets01 <- ets(y = crime_train,
model = "AZZ",
alpha = 0.1)tuning model dengan alpha = 0.9 sebagai pembanding batas maksimal
crime_ets09 <- ets(y = crime_train,
model = "AZZ",
alpha = 0.9)Melihat Hasil Model SES
Kita dapat melihat hasil model dengan memanggil kolom fitted
# Hasil Model SES
crime_ets$fitted#> Jan Feb Mar Apr May Jun Jul
#> 2020 10.655774 11.649670 11.456589 10.429300 8.221332 9.047145 10.816318
#> 2021 8.959820 7.188577 6.835335 7.775865 7.842477 8.186490 8.131066
#> 2022 7.450342 7.316502 6.628043 6.738588 6.221884 7.641928 8.342741
#> Aug Sep Oct Nov Dec
#> 2020 9.979315 10.282660 11.387444 11.569494 10.211452
#> 2021 7.200521 7.735322 7.516786 7.066001 7.640781
#> 2022 10.321261 10.225783 11.347470 11.244203 11.468824- Evaluasi hasil forecast
Kita akan membandingkan ketiga hasil model ets di atas yang telah kita buat, dan akan kita bandingkan MAPE (Mean Average Percentage Error)
accuracy(crime_ets$fitted,crime_ts)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set 0.03519162 3.121312 2.431561 -14.88481 35.35886 0.1617885 0.6850659accuracy(crime_ets01$fitted,crime_ts)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set 0.07167957 3.192858 2.58273 -16.28343 37.48661 0.3211142 0.6911783accuracy(crime_ets09$fitted,crime_ts)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set -0.107251 3.366726 2.849386 -11.61175 37.12904 -0.2460957 0.9511432Berdasarkan hasil evaluasi forecast di atas MAPE yang paling kecil adalah model tanpa tuning, yaitu mencapai angka 35.3%.
Model HoltWinters
Model HoltWinters adalah model yang akan kita buat untuk data crime. Model ini dipilih karena kemampuan tuning dan bisa digunakan pada model yang memiliki tren dan seasonal
- Membuat Model
Berikut adalah model dasar holtwinters yang akan kita buat
crime_hw <- HoltWinters(x = crime_train,
seasonal = "additive")- Melakukan forecast
Forecast dibuat berdasarkan model holtwinters yang telah kita buat
yaitu crime_hw
crime_forecast <- forecast(crime_hw, h = 12)- Evaluasi hasil forecast
Evaluasi hasil forecast adalah data mean dari forecast dari model holtwinters
# Please type your code
crime_forecast$mean#> Jan Feb Mar Apr May Jun Jul
#> 2023 7.967139 8.385932 11.997989 10.863914 14.879777 13.885013 14.749688
#> Aug Sep Oct Nov Dec
#> 2023 13.788546 16.520511 14.645371 14.833524 12.916871Selanjutnya kita akan menghitung error dari model holtwinters kita
accuracy(crime_forecast$mean,crime_test)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set -1.56475 6.190776 5.415599 -99.10708 121.1926 -0.2438704 0.6163527Berdasarkan hasil MAPE, menunjukkan nilai yang sangat besar yaitu 121.1%. Hal ini sangat tidak diinginkan, maka dari itu kita perlu menambahkan tuning khususnya karena data kita memiliki pola tren dan seasonal
crime_hw_iterate <- HoltWinters(x = crime_train,
seasonal = "additive",
alpha = 0.001,
beta = 0.001,
gamma = 0.0001)
crime_forecast_iterate <- forecast(crime_hw_iterate, h = 12)
accuracy(crime_forecast_iterate$mean,crime_test)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set 7.218427 10.40642 7.920906 35.77622 70.90017 -0.1721666 1.007171Setelah melakukan beberapa iterasi, ditemukan MAPE senilai 70.9%. Error sudah menurun dibandingkan model sebelumnya, akan tetapi masih terlampau sangat besar dibandingkan model exponential smoothing yang telah kita buat.
SUMMARY
Berdasarkan data dan model yang telah dibuat, dapat disimpulkan bahwa
perlu ada analisis lebih mendalam terkait model crime dengan tipe
kidnapping. Terutama pada tahap subsetting, sangat menentukan bagaimana
kita memilih pola tersebut memiliki pola tren dan seasonal. Berdasarkan
data crime_ts yang telah disubset menjadi lebih singkat
dari data tahun 2020-2023, dapat dilihat bahwa data tersebut memiliki
tren naik dan turun, yang mana data tren turun tersebut dijadikan data
train. Hal ini dapat menimbulkan hasil yang berbeda apabila kita memilih
data tahun 2015-2020, dikarenakan pada rentang waktu tersebut
kemungkinan besar tidak memiliki tren, namun bisa jadi memiliki pola
seasonal.
Berdasarkan dua model di atas, exponential smoothing dan model
HoltWinters, model yang lebih baik adalah model exponential smoothing
tanpa tuning. Hal ini mungkin terjadi karena ada kemungkinan kesalahan
perkiraan tren pada model crime_ts di rentang tahun
2020-2023. Kemungkinan lebih baik adalah apabila model dianggap tidak
memiliki tren dan seasonal sama sekali namun pernyataan ini perlu
dibuktikan dengan metode model lain, contohnya model Simple Moving
Average.