Crime Rate in US
Justinus Dedy Handyka Simanjuntak
3/26/2022
Introduction
kriminalitas adalah perilaku atau tindakan yang melanggar norma hukum legal atau formal. Contohnya perampokan, pembunuhan, penipuan dan lainnya. Pada kali ini saya akan melakukan prediksi tingkat kejahatan crime rate menggunakan Regresi Linear dan akan membandingkan model mana yang terbaik. Data yang akan digunakan adalah data laporan kriminalitas di beberapa negara bagian di United States pada tahun 1960.
Libarary
library(tidyverse)
library(dplyr)
library(GGally)
library(lmtest)
library(car)
library(MLmetrics)Data Preparation
crime <- read.csv("crime.csv") %>% select(-X)
names(crime) <- c("percent_m", "is_south", "mean_education", "police_exp60", "police_exp59", "labour_participation", "m_per1000f", "state_pop", "nonwhites_per1000", "unemploy_m24", "unemploy_m39", "gdp", "inequality", "prob_prison", "time_prison", "crime_rate")
head(crime)Exploratory Data Analysis
str(crime)## 'data.frame': 47 obs. of 16 variables:
## $ percent_m : int 151 143 142 136 141 121 127 131 157 140 ...
## $ is_south : int 1 0 1 0 0 0 1 1 1 0 ...
## $ mean_education : int 91 113 89 121 121 110 111 109 90 118 ...
## $ police_exp60 : int 58 103 45 149 109 118 82 115 65 71 ...
## $ police_exp59 : int 56 95 44 141 101 115 79 109 62 68 ...
## $ labour_participation: int 510 583 533 577 591 547 519 542 553 632 ...
## $ m_per1000f : int 950 1012 969 994 985 964 982 969 955 1029 ...
## $ state_pop : int 33 13 18 157 18 25 4 50 39 7 ...
## $ nonwhites_per1000 : int 301 102 219 80 30 44 139 179 286 15 ...
## $ unemploy_m24 : int 108 96 94 102 91 84 97 79 81 100 ...
## $ unemploy_m39 : int 41 36 33 39 20 29 38 35 28 24 ...
## $ gdp : int 394 557 318 673 578 689 620 472 421 526 ...
## $ inequality : int 261 194 250 167 174 126 168 206 239 174 ...
## $ prob_prison : num 0.0846 0.0296 0.0834 0.0158 0.0414 ...
## $ time_prison : num 26.2 25.3 24.3 29.9 21.3 ...
## $ crime_rate : int 791 1635 578 1969 1234 682 963 1555 856 705 ...Deskripsi kolom:
percent_m: percentage of males aged 14-24is_south: whether it is in a Southern state. 1 for Yes, 0 for No.mean_education: mean years of schoolingpolice_exp60: police expenditure in 1960police_exp59: police expenditure in 1959labour_participation: labour force participation ratem_per1000f: number of males per 1000 femalesstate_pop: state populationnonwhites_per1000: number of non-whites resident per 1000 peopleunemploy_m24: unemployment rate of urban males aged 14-24unemploy_m39: unemployment rate of urban males aged 35-39gdp: gross domestic product per headinequality: income inequalityprob_prison: probability of imprisonmenttime_prison: average time served in prisonscrime_rate: crime rate in an unspecified category
Dapat dilihat pada variabel is_south tipe datanya masih belum sesuai sehingga kita harus mengubahnya ke tipe data yang lebih sesuai. Kemudian kita dapat melakukan pengecekan terhadap missing value
# ubah tipe data
crime <- crime %>%
mutate(is_south=as.factor(is_south))
#cek missing value
anyNA(crime)## [1] FALSECek korelasi variable menggunakan fungsi ggcorr()
ggcorr(crime, label = T, hjust = 1, layout.exp = 3)
Pada grafik diatas variabel yang berkorelasi kuat dengan crime rate adalah police_exp59 dan police_exp60
Model Regresi Linear
Sebelum masuk kedalam pembuatan model kita splitting data menjadi data train dan data test
RNGkind(sample.kind = "Rounding")
set.seed(123)
intrain <- sample(nrow(crime), nrow(crime)*0.8)
c_train <- crime[intrain, ]
c_test <- crime[-intrain, ]Selanjutnya kita akan membuat model linear sederhana dengan salah satu variabel yang memiliki korelasi tertinggi yaitu police_exp59
model_crime <- lm(crime_rate ~ police_exp59, c_train)
summary(model_crime)##
## Call:
## lm(formula = crime_rate ~ police_exp59, data = c_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -592.18 -150.07 17.86 147.27 588.68
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 194.588 155.075 1.255 0.218
## police_exp59 8.966 1.775 5.050 1.38e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 316.9 on 35 degrees of freedom
## Multiple R-squared: 0.4215, Adjusted R-squared: 0.405
## F-statistic: 25.5 on 1 and 35 DF, p-value: 1.38e-05kita mendapatkan nilai R-squared sebesar 0.4215. Model dapat menjelaskan 42.15% dari data, sedangkan sisanya sebanyak 57.85% dijelaskan oleh variabel atau faktor-faktor lain yang tidak kita masukkan ke dalam model. Idealnya kita harus mendapat nilai R-squared yang lebih baik lagi dengan cara menambah kan variabel - variabel lain.
Selanjutnya akan dicoba pemilihan variabel prediktor secara automatis menggunakan step-wise regression dengan metode backward elimination
model_full <- lm(crime_rate ~ ., c_train)
model_step <- step(model_full, direction = "backward")## Start: AIC=414.49
## crime_rate ~ percent_m + is_south + mean_education + police_exp60 +
## police_exp59 + labour_participation + m_per1000f + state_pop +
## nonwhites_per1000 + unemploy_m24 + unemploy_m39 + gdp + inequality +
## prob_prison + time_prison
##
## Df Sum of Sq RSS AIC
## - labour_participation 1 653 1143026 412.52
## - time_prison 1 1921 1144294 412.56
## - police_exp59 1 2875 1145248 412.59
## - state_pop 1 3010 1145383 412.59
## - is_south 1 8757 1151130 412.78
## - nonwhites_per1000 1 13116 1155489 412.92
## - gdp 1 13142 1155515 412.92
## - police_exp60 1 28035 1170408 413.39
## - m_per1000f 1 45287 1187660 413.93
## <none> 1142373 414.49
## - prob_prison 1 70404 1212777 414.71
## - percent_m 1 97608 1239981 415.53
## - unemploy_m24 1 103796 1246169 415.71
## - mean_education 1 200844 1343218 418.49
## - unemploy_m39 1 228152 1370525 419.23
## - inequality 1 244314 1386687 419.67
##
## Step: AIC=412.52
## crime_rate ~ percent_m + is_south + mean_education + police_exp60 +
## police_exp59 + m_per1000f + state_pop + nonwhites_per1000 +
## unemploy_m24 + unemploy_m39 + gdp + inequality + prob_prison +
## time_prison
##
## Df Sum of Sq RSS AIC
## - police_exp59 1 2351 1145376 410.59
## - time_prison 1 2678 1145704 410.60
## - state_pop 1 4131 1147157 410.65
## - is_south 1 8388 1151414 410.79
## - gdp 1 12638 1155664 410.92
## - nonwhites_per1000 1 13220 1156246 410.94
## - police_exp60 1 27439 1170465 411.39
## - m_per1000f 1 52705 1195731 412.18
## <none> 1143026 412.52
## - prob_prison 1 71104 1214130 412.75
## - percent_m 1 103075 1246101 413.71
## - unemploy_m24 1 121084 1264110 414.24
## - mean_education 1 200303 1343329 416.49
## - unemploy_m39 1 227521 1370547 417.23
## - inequality 1 246070 1389096 417.73
##
## Step: AIC=410.59
## crime_rate ~ percent_m + is_south + mean_education + police_exp60 +
## m_per1000f + state_pop + nonwhites_per1000 + unemploy_m24 +
## unemploy_m39 + gdp + inequality + prob_prison + time_prison
##
## Df Sum of Sq RSS AIC
## - state_pop 1 3598 1148975 408.71
## - time_prison 1 6058 1151435 408.79
## - is_south 1 8688 1154064 408.87
## - nonwhites_per1000 1 11027 1156404 408.95
## - gdp 1 11041 1156418 408.95
## <none> 1145376 410.59
## - prob_prison 1 68965 1214341 410.76
## - m_per1000f 1 88286 1233663 411.34
## - percent_m 1 102660 1248036 411.77
## - unemploy_m24 1 122681 1268058 412.36
## - mean_education 1 203773 1349149 414.65
## - unemploy_m39 1 233268 1378645 415.45
## - inequality 1 245177 1390554 415.77
## - police_exp60 1 358986 1504363 418.68
##
## Step: AIC=408.71
## crime_rate ~ percent_m + is_south + mean_education + police_exp60 +
## m_per1000f + nonwhites_per1000 + unemploy_m24 + unemploy_m39 +
## gdp + inequality + prob_prison + time_prison
##
## Df Sum of Sq RSS AIC
## - time_prison 1 3983 1152958 406.84
## - gdp 1 9092 1158067 407.00
## - is_south 1 9896 1158871 407.03
## - nonwhites_per1000 1 12118 1161093 407.10
## <none> 1148975 408.71
## - prob_prison 1 68307 1217282 408.84
## - percent_m 1 105640 1254614 409.96
## - unemploy_m24 1 142101 1291075 411.02
## - m_per1000f 1 152198 1301173 411.31
## - mean_education 1 206104 1355079 412.81
## - unemploy_m39 1 246722 1395697 413.91
## - inequality 1 274184 1423159 414.63
## - police_exp60 1 402260 1551235 417.81
##
## Step: AIC=406.84
## crime_rate ~ percent_m + is_south + mean_education + police_exp60 +
## m_per1000f + nonwhites_per1000 + unemploy_m24 + unemploy_m39 +
## gdp + inequality + prob_prison
##
## Df Sum of Sq RSS AIC
## - is_south 1 10049 1163007 405.16
## - gdp 1 10978 1163937 405.19
## - nonwhites_per1000 1 17994 1170952 405.41
## <none> 1152958 406.84
## - percent_m 1 111548 1264506 408.25
## - prob_prison 1 142080 1295039 409.14
## - unemploy_m24 1 144421 1297379 409.20
## - m_per1000f 1 148215 1301173 409.31
## - mean_education 1 202307 1355265 410.82
## - unemploy_m39 1 251866 1404825 412.15
## - inequality 1 271850 1424809 412.67
## - police_exp60 1 398984 1551942 415.83
##
## Step: AIC=405.16
## crime_rate ~ percent_m + mean_education + police_exp60 + m_per1000f +
## nonwhites_per1000 + unemploy_m24 + unemploy_m39 + gdp + inequality +
## prob_prison
##
## Df Sum of Sq RSS AIC
## - gdp 1 11334 1174342 403.52
## - nonwhites_per1000 1 13901 1176909 403.60
## <none> 1163007 405.16
## - percent_m 1 103127 1266134 406.30
## - unemploy_m24 1 135954 1298961 407.25
## - m_per1000f 1 148204 1311211 407.60
## - prob_prison 1 213993 1377000 409.41
## - mean_education 1 231928 1394935 409.89
## - unemploy_m39 1 242747 1405754 410.17
## - inequality 1 262973 1425981 410.70
## - police_exp60 1 389276 1552284 413.84
##
## Step: AIC=403.52
## crime_rate ~ percent_m + mean_education + police_exp60 + m_per1000f +
## nonwhites_per1000 + unemploy_m24 + unemploy_m39 + inequality +
## prob_prison
##
## Df Sum of Sq RSS AIC
## - nonwhites_per1000 1 9466 1183808 401.81
## <none> 1174342 403.52
## - percent_m 1 95312 1269653 404.40
## - unemploy_m24 1 143492 1317833 405.78
## - m_per1000f 1 177963 1352305 406.74
## - mean_education 1 227959 1402300 408.08
## - prob_prison 1 242508 1416850 408.46
## - unemploy_m39 1 249410 1423751 408.64
## - inequality 1 286285 1460626 409.59
## - police_exp60 1 616530 1790872 417.13
##
## Step: AIC=401.81
## crime_rate ~ percent_m + mean_education + police_exp60 + m_per1000f +
## unemploy_m24 + unemploy_m39 + inequality + prob_prison
##
## Df Sum of Sq RSS AIC
## <none> 1183808 401.81
## - unemploy_m24 1 134893 1318700 403.81
## - percent_m 1 149388 1333195 404.21
## - m_per1000f 1 169039 1352847 404.75
## - mean_education 1 220608 1404416 406.14
## - prob_prison 1 237275 1421082 406.57
## - unemploy_m39 1 239975 1423783 406.64
## - inequality 1 461989 1645796 412.00
## - police_exp60 1 1025300 2209107 422.90summary(model_full)##
## Call:
## lm(formula = crime_rate ~ ., data = c_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -356.59 -131.61 -6.69 127.93 415.48
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7334.2000 2327.6113 -3.151 0.00482 **
## percent_m 7.1664 5.3500 1.340 0.19472
## is_south1 -79.9456 199.2543 -0.401 0.69231
## mean_education 18.0178 9.3770 1.921 0.06835 .
## police_exp60 12.4150 17.2940 0.718 0.48074
## police_exp59 -4.3954 19.1200 -0.230 0.82041
## labour_participation -0.2101 1.9179 -0.110 0.91381
## m_per1000f 3.3102 3.6280 0.912 0.37191
## state_pop -0.4162 1.7697 -0.235 0.81632
## nonwhites_per1000 0.4357 0.8874 0.491 0.62850
## unemploy_m24 -7.8713 5.6984 -1.381 0.18170
## unemploy_m39 23.0303 11.2456 2.048 0.05328 .
## gdp 0.6426 1.3074 0.492 0.62816
## inequality 6.4449 3.0411 2.119 0.04616 *
## prob_prison -3304.2217 2904.4492 -1.138 0.26808
## time_prison 1.8181 9.6743 0.188 0.85273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 233.2 on 21 degrees of freedom
## Multiple R-squared: 0.812, Adjusted R-squared: 0.6777
## F-statistic: 6.046 on 15 and 21 DF, p-value: 0.0001131summary(model_step)##
## Call:
## lm(formula = crime_rate ~ percent_m + mean_education + police_exp60 +
## m_per1000f + unemploy_m24 + unemploy_m39 + inequality + prob_prison,
## data = c_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -394.74 -125.65 27.95 128.24 465.39
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7289.928 1383.644 -5.269 1.33e-05 ***
## percent_m 7.636 4.062 1.880 0.0706 .
## mean_education 16.603 7.269 2.284 0.0301 *
## police_exp60 9.219 1.872 4.925 3.41e-05 ***
## m_per1000f 3.712 1.857 2.000 0.0553 .
## unemploy_m24 -7.495 4.196 -1.786 0.0849 .
## unemploy_m39 22.563 9.471 2.382 0.0242 *
## inequality 5.625 1.702 3.306 0.0026 **
## prob_prison -3837.725 1619.981 -2.369 0.0250 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 205.6 on 28 degrees of freedom
## Multiple R-squared: 0.8052, Adjusted R-squared: 0.7495
## F-statistic: 14.46 on 8 and 28 DF, p-value: 4.236e-08Metode step-wise regression menghasilkan formula optimum dengan nilai AIC yang terendah, dimana semakin rendah nilai AIC tersebut, maka nilai observasi yang tidak tertangkap semakin kecil. Dari model_step kita mendapatkan nilai adjusted R-squared sebesar 0.7495 lebih baik dibandingkan model_full yaitu 0.6777.
Predict
pred_crime1 <- predict(model_full, c_test)
pred_crime2 <- predict(model_step, c_test)
data.frame(MAPE1 = MAPE(pred_crime1, c_train$crime_rate),
MAPE2 = MAPE(pred_crime2, c_train$crime_rate)
)## Warning in y_true - y_pred: longer object length is not a multiple of shorter
## object length
## Warning in y_true - y_pred: longer object length is not a multiple of shorter
## object lengthBerdasarkan nilai MAPE kedual model tersebut, model step memiliki MAPE yang lebih kecil yaitu 0.4531014 dibanding model full.
Evaluasi Model
Normalitas
shapiro.test(model_step$residuals)##
## Shapiro-Wilk normality test
##
## data: model_step$residuals
## W = 0.98774, p-value = 0.9498Dapat kita lihat bahwa p-value = 0.9498 yang berarti > 0.05 sehingga H0 diterima yang berarti residual atau error terdistribusi secara normal.
Homocesdasticity
Kita gunakan uji statistik Breusch-Pagan untuk melihat apakah error bersifat konstan atau tidak Penentuan :
p-value > 0.05: Tidak menolak H0 atau error bersifat konstan (Homocesdasticity)
p-value < 0.05: Menolak H0 atau error bersifat tidak konstan (Heteroscesdasticity
bptest(model_step)##
## studentized Breusch-Pagan test
##
## data: model_step
## BP = 11.85, df = 8, p-value = 0.158Dari uji diatas menunjukkan p-value > 0.05 sehingga H0 diterima dan dapat diartikan residual tidak membentuk pola tertentu sehingga kita bisa menyimpulkan bahwa model sudah memiliki error yang bersifat konstan.
Multicollinearity
library(car)
vif(model_step)## percent_m mean_education police_exp60 m_per1000f unemploy_m24
## 2.298626 5.609506 2.990798 2.276976 4.760707
## unemploy_m39 inequality prob_prison
## 5.522660 4.416847 1.392422Tidak ada nilai sama dengan atau lebih dari 10 sehingga tidak ditemukan Multicollinearity antar variabel (antar variabel prediktor saling independen).
Conclusion
Dari kedua model yang kita dapat, model step adalah model terbaik dengan nilai adjusted R-squared sebesar 0.7495 dan nilai MAPE sebesar 0.4531014. Dan juga melalui evaluasi model dengan menguji(Normalitas, Homocesdasticity, Multicollinearity) dapat dilihat bahwa model ini memiliki kriteria yang baik