1 Introduction

1.1 Our Purpose

Data crime.csv merupakan data sosio-demografi dari 47 negara bagian US di tahun 1960. Secara umum, data ini berisi tentang tingkat kejahatan atau kriminalitas di negara bagian US pada tahun 1960, yang juga didukung dengan data-data lain yang diduga mendukung tingkat kriminalitas tersebut, seperti :

  • gdp : gross domestic product -> produk domestik bruto, dimana erat kaitannya dengan kesehatan ekonomi suatu wilayah.
  • inequality : income inequality -> kesenjangan pendapatan.
  • state_pop : state population -> populasi penduduk di negara bagian US. dan variabel-variabel lainnya.

Pada data ini, terdapat juga variabel crime_rate yang menggambarkan ukuran tingkat kejahatan di Amerika Serikat pada tahun 1960 untuk berbagai negara bagian. Disinilah saya, sebagai seorang data analis, akan melihat bagaimana kondisi sosio-ekonomi (yang tercermin pada variabel-variabel sosio-demografi tersebut) dapat mencerminkan tingkat kejahatan di suatu negara bagian.

1.2 Preparation

Kita panggil setiap library yang dibutuhkan dalam proses pembuatan model regresi linear hingga evaluasi dan hasil interpretasi model nantinya.

# EDA
library(dplyr)
library(GGally)

# evaluasi model
library(lmtest)
library(car)
library(MLmetrics)

2 Data Preparation

2.1 Read Dataset

Kita mulai dengan membaca dataset crime.csv dan kita simpan dalam variabel crime.

crime <- read.csv("crime.csv")
head(crime)

Kemudian kita cek setiap tipe data dengan fungsi str.

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 ...

Berikut ini adalah penjelasan terkait variabel-variabel tersebut :

  • percent_m: persentase pria berusia 14-24
  • is_south: apakah lokasi ini berada di negara bagian selatan? 1 untuk Ya, 0 untuk Tidak
  • mean_education: rata-rata usia pendidikan
  • police_exp60: pengeluaran (anggaran) kepolisian di tahun 1960
  • police_exp59: pengeluaran (anggaran) kepolisian di tahun 1959
  • labour_participation: tingkat partisipasi penduduk usia kerja
  • m_per1000f: jumlah pria per 1.000 wanita
  • state_pop: populasi penduduk di negara bagian
  • nonwhites_per1000: jumlah penduduk tidak berkulit putih per 1.000 orang
  • unemploy_m24: tingkat pengangguran pria perkotaan berusia 14-24 tahun
  • unemploy_m39: tingkat pengangguran pria perkotaan berusia 35-39 tahun
  • gdp: produk domestik bruto per kepala keluarga
  • inequality: kesenjangan pendapatan
  • prob_prison: kemungkinan (probabilitas) hukuman penjara / orang dipenjara
  • time_prison: rata-rata orang menjalani hukuman penjara
  • crime_rate: tingkat kejahatan dalam kategori yang tidak ditentukan

Berdasarkan penjelasan dari setiap variabel tersebut, terdapat 1 variabel yang belum sesuai, yaitu is_south yang seharusnya bertipe factor.

2.2 Data Cleansing

Ubah tipe data variabel is_south dari integer menjadi factor.

crime <- crime %>% mutate(is_south = as.factor(is_south))
str(crime)
#> 'data.frame':    47 obs. of  16 variables:
#>  $ percent_m           : int  151 143 142 136 141 121 127 131 157 140 ...
#>  $ is_south            : Factor w/ 2 levels "0","1": 2 1 2 1 1 1 2 2 2 1 ...
#>  $ 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 ...

3 Exploratory Data Analysis

Disini kita mulai memeriksa korelasi antara setiap variabel. Khususnya kita akan melihat variabel apa saja yang akan memiliki korelasi yang kuat dengan variabel target kita nantinya yaitu variabel crime_rate.

ggcorr(crime, hjust = 1, layout.exp = 3, label = TRUE)

Berdasarkan grafik eksplorasi antar variabel di atas, dengan target variabel crime_rate didapatkan insight sebagai berikut :

  1. Variabel police_exp60 dan police_exp59 memiliki korelasi kuat secara positif karena nilainya 0.7
  2. Variabel gdp memiliki korelasi sedang secara positif dengan nilai 0.4
  3. Variabel prob_prison memiliki korelasi sedang secara negatif dengan nilai -0.4
  4. Variabel-variabel lainnya bisa disimpulkan memiliki korelasi yang lemah (dengan nilai korelasi di bawah 0.4) terhadap variabel crime_rate
  • Dengan demikian, semakin besar nilai police_exp60, police_exp59 dan gdp, akan semakin besar pula nilai crime_rate
  • Namun, semakin besar nilai prob_prison, semakin kecil nilai crime_rate

4 Modeling

Selanjutnya untuk melakukan proses analisis lebih lanjut dengan target variabel crime_rate, kita membuat model regresi linear dengan 3 metode berikut ini.

4.1 Simple Linear Regression

Disini, kita membuat model regresi linear dengan menggunakan 1 variabel prediktor. Kita memilih variabel dengan korelasi paling kuat. Namun kita melihat terdapat 2 variabel dengan korelasi yang sama kuat (dengan nilai 0.7) yaitu police_exp60 dan police_exp59. Jika kita melihat pada bagian EDA pada saat memahami masing-masing variabel, kedua variabel tersebut mengacu pada nilai yang sama yaitu police expenditure atau pengeluaran / anggaran kepolisian yang berada di tahun yang berbeda yaitu tahun 1960 dan 1959. Oleh karena itu, kita cukup memilih salah satu dari kedua variabel tersebut untuk kita gunakan dalam simple linear regression model ini. Kita coba menggunakan variabel prediktor police_exp60.

model_crime_simple <- lm(crime_rate ~ police_exp60, data = crime)

summary(model_crime_simple)
#> 
#> Call:
#> lm(formula = crime_rate ~ police_exp60, data = crime)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -586.91 -155.63   32.52  139.58  568.84 
#> 
#> Coefficients:
#>              Estimate Std. Error t value     Pr(>|t|)    
#> (Intercept)   144.464    126.693   1.140         0.26    
#> police_exp60    8.948      1.409   6.353 0.0000000934 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 283.9 on 45 degrees of freedom
#> Multiple R-squared:  0.4728, Adjusted R-squared:  0.4611 
#> F-statistic: 40.36 on 1 and 45 DF,  p-value: 0.00000009338

Interpretasi model model_crime_simple :

  • Intercept : jika tidak ada nilai police_exp60, maka nilai crime_rate adalah 144.464
  • Slope : setiap kenaikan 1 police_exp60, menambahkan nilai crime_rate sebesar 8.948 (koefisien positif)
  • P-value : police_exp60 memiliki nilai p-value < 0.05 (alpha), dan memiliki signifikansi codes (tanda bintang), hal ini menunjukkan bahwa police_exp60 merupakan prediktor yang signifikan (berpengaruh)
  • R-squared : nilai Multiple R-squared (1 prediktor) 0.4728 sehingga dikatakan model ini belum cukup baik (model yang baik nilai R-squared mendekati 1)

4.2 Multiple Linear Regression

4.2.1 Model dengan beberapa variabel prediktor

Disini, kita mencoba menggabungkan beberapa variabel yang memiliki korelasi yang kuat dan sedang hasil dari EDA yang sudah kita lakukan sebelumnya. Jadi kita akan menggunakan police_exp60 (police_exp59 dapat diwakilkan dengan variabel ini), gdp dan prob_prison sebagai variabel prediktor.

model_crime_multiple <- lm(crime_rate ~ police_exp60 + gdp + prob_prison, data = crime)
summary(model_crime_multiple)
#> 
#> Call:
#> lm(formula = crime_rate ~ police_exp60 + gdp + prob_prison, data = crime)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -706.64  -85.35   14.62  137.45  506.92 
#> 
#> Coefficients:
#>                Estimate Std. Error t value   Pr(>|t|)    
#> (Intercept)    887.4747   341.7691   2.597     0.0128 *  
#> police_exp60    11.3728     2.2076   5.152 0.00000615 ***
#> gdp             -1.4740     0.7202  -2.047     0.0469 *  
#> prob_prison  -3709.4279  2139.6191  -1.734     0.0901 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 273.7 on 43 degrees of freedom
#> Multiple R-squared:  0.5318, Adjusted R-squared:  0.4991 
#> F-statistic: 16.28 on 3 and 43 DF,  p-value: 0.0000003246

Interpretasi model model_crime_multiple :

  • Intercept : jika tidak ada nilai police_exp60, gdp, dan prob_prison, maka nilai crime_rate adalah 887.4747
  • Slope dan p-value (signifikansi) :
    • variabel police_exp60 merupakan slope positif, dimana secara statistik signifikan (dengan p-value lebih kecil 0.05)
    • variabel gdp dan prob_prison merupakan slope negatif, dimana secara statistik tidak signifikan (dengan p-value lebih besar 0.05)
  • R-squared : nilai Adjusted R-squared (lebih dari 1 prediktor) 0.4991 sehingga dikatakan model ini belum cukup baik

4.2.2 Model dengan semua variabel sebagai prediktor

Disini, kita mencoba menggunakan semua variabel sebagai prediktor untuk memprediksi nilai crime_rate.

model_crime_all <- lm(crime_rate ~ ., data = crime)
summary(model_crime_all)
#> 
#> Call:
#> lm(formula = crime_rate ~ ., data = crime)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -395.74  -98.09   -6.69  112.99  512.67 
#> 
#> Coefficients:
#>                        Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)          -5984.2876  1628.3184  -3.675 0.000893 ***
#> percent_m                8.7830     4.1714   2.106 0.043443 *  
#> is_south1               -3.8035   148.7551  -0.026 0.979765    
#> mean_education          18.8324     6.2088   3.033 0.004861 ** 
#> police_exp60            19.2804    10.6110   1.817 0.078892 .  
#> police_exp59           -10.9422    11.7478  -0.931 0.358830    
#> labour_participation    -0.6638     1.4697  -0.452 0.654654    
#> m_per1000f               1.7407     2.0354   0.855 0.398995    
#> state_pop               -0.7330     1.2896  -0.568 0.573845    
#> nonwhites_per1000        0.4204     0.6481   0.649 0.521279    
#> unemploy_m24            -5.8271     4.2103  -1.384 0.176238    
#> unemploy_m39            16.7800     8.2336   2.038 0.050161 .  
#> gdp                      0.9617     1.0367   0.928 0.360754    
#> inequality               7.0672     2.2717   3.111 0.003983 ** 
#> prob_prison          -4855.2658  2272.3746  -2.137 0.040627 *  
#> time_prison             -3.4790     7.1653  -0.486 0.630708    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 209.1 on 31 degrees of freedom
#> Multiple R-squared:  0.8031, Adjusted R-squared:  0.7078 
#> F-statistic: 8.429 on 15 and 31 DF,  p-value: 0.0000003539

Interpretasi model model_crime_all : - Signifikansi : - terdapat 4 variabel yang secara statistik signifikan (p-value < 0.05) yaitu : percent_m, mean_education, inequality, dan prob_prison - variabel police_exp60 dan police_exp59 memiliki nilai p-value > 0.05, yang menandakan variabel ini tidak signifikan untuk model ini - R-squared : nilai Adjusted R-squared 0.7078 sehingga dikatakan model ini sudah baik (mendekati 1)

Berdasarkan hasil EDA, variabel police_exp60 dan police_exp59 seharusnya merupakan variabel yang korelasinya paling tinggi terhadap crime_rate. Namun, pada model ini yang menggunakan seluruh variabel, justru kedua variabel prediktor tersebut nilainya tidak signifikan.

Oleh karena itu, kita coba membuat model dengan menggunakan semua variabel, namun menghilangkan salah satu variabel tersebut di atas, dalam hal ini kita coba menghilangkan variabel police_exp59.

crime_clean <- crime %>% select(-c('police_exp59'))
model_crime_clean_all <- lm(crime_rate ~ ., data = crime_clean)
summary(model_crime_clean_all)
#> 
#> Call:
#> lm(formula = crime_rate ~ ., data = crime_clean)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -442.55 -116.46    8.86  118.26  473.49 
#> 
#> Coefficients:
#>                        Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)          -6379.0457  1568.9382  -4.066 0.000291 ***
#> percent_m                8.9861     4.1571   2.162 0.038232 *  
#> is_south1                5.6688   148.0997   0.038 0.969705    
#> mean_education          17.7304     6.0824   2.915 0.006445 ** 
#> police_exp60             9.6526     2.3921   4.035 0.000317 ***
#> labour_participation    -0.2801     1.4079  -0.199 0.843538    
#> m_per1000f               1.8218     2.0293   0.898 0.376026    
#> state_pop               -0.7836     1.2857  -0.609 0.546523    
#> nonwhites_per1000        0.2446     0.6187   0.395 0.695239    
#> unemploy_m24            -5.4163     4.1784  -1.296 0.204164    
#> unemploy_m39            16.9362     8.2148   2.062 0.047441 *  
#> gdp                      0.9072     1.0329   0.878 0.386292    
#> inequality               7.2708     2.2564   3.222 0.002921 ** 
#> prob_prison          -4285.1992  2183.8679  -1.962 0.058484 .  
#> time_prison             -1.1276     6.6919  -0.168 0.867251    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 208.6 on 32 degrees of freedom
#> Multiple R-squared:  0.7976, Adjusted R-squared:  0.709 
#> F-statistic: 9.006 on 14 and 32 DF,  p-value: 0.0000001673

Interpretasi model model_crime_clean_all : - Signifikansi : - terdapat 5 variabel yang secara statistik signifikan (p-value < 0.05) yaitu : percent_m, mean_education, police_exp60, unemploy_m39 dan inequality. - variabel prediktor police_exp60 merupakan variabel yang paling signifikan (nilai p-value terkecil). Hal ini sesuai dengan nilai korelasi yang tinggi untuk variabel ini pada saat EDA. - R-squared : nilai Adjusted R-squared 0.709 sehingga dikatakan model ini sudah baik (mendekati 1).

Jika melihat nilai R-squared pada model model_crime_all dan model_crime_clean_all tidak jauh berbeda, sehingga kedua model tersebut sama baiknya. Yang membedakan adalah variabel prediktor yang membentuk kedua model tersebut.

4.3 Step-wise Regression for Feature Selection

Berdasarkan hasil analisis dari kedua metode model tersebut di atas, sementara ini dapat disimpulkan bahwa untuk dapat meningkatkan ukuran kebaikan model (Goodness of Fit), yang ditandai dengan nilai R-squared, kita dapat menambahkan lebih banyak variabel prediktor. Hal ini ditunjukkan dengan model yang nilai R-squared tinggi adalah model dengan prediktor yang paling banyak.

Salah satu teknik lain yang dapat kita gunakan dalam memilih variabel prediktor dalam membentuk model adalah dengan algoritma step-wise regression. Kita akan mencoba menggunakan metode ini dengan parameter direction = backward() dengan menggunakan data crime_clean. Kita menggunakan data crime_clean agar kita cukup menggunakan police_exp60 dan tidak menggunakan police_exp59, karena kedua nilai variabel tersebut secara konsep nilainya sama dan berpengaruh signifikan terhadap variabel crime_rate.

model_crime_clean_step <- step(object = model_crime_clean_all,
                               direction = "backward")
#> Start:  AIC=513.95
#> crime_rate ~ percent_m + is_south + mean_education + police_exp60 + 
#>     labour_participation + m_per1000f + state_pop + nonwhites_per1000 + 
#>     unemploy_m24 + unemploy_m39 + gdp + inequality + prob_prison + 
#>     time_prison
#> 
#>                        Df Sum of Sq     RSS    AIC
#> - is_south              1        64 1392929 511.95
#> - time_prison           1      1236 1394101 511.99
#> - labour_participation  1      1723 1394588 512.00
#> - nonwhites_per1000     1      6802 1399667 512.18
#> - state_pop             1     16168 1409033 512.49
#> - gdp                   1     33582 1426447 513.07
#> - m_per1000f            1     35080 1427945 513.12
#> <none>                              1392865 513.95
#> - unemploy_m24          1     73136 1466001 514.35
#> - prob_prison           1    167590 1560455 517.29
#> - unemploy_m39          1    185009 1577874 517.81
#> - percent_m             1    203389 1596254 518.35
#> - mean_education        1    369864 1762729 523.01
#> - inequality            1    451937 1844802 525.15
#> - police_exp60          1    708738 2101603 531.28
#> 
#> Step:  AIC=511.95
#> crime_rate ~ percent_m + mean_education + police_exp60 + labour_participation + 
#>     m_per1000f + state_pop + nonwhites_per1000 + unemploy_m24 + 
#>     unemploy_m39 + gdp + inequality + prob_prison + time_prison
#> 
#>                        Df Sum of Sq     RSS    AIC
#> - time_prison           1      1319 1394247 509.99
#> - labour_participation  1      2646 1395574 510.04
#> - nonwhites_per1000     1      8949 1401878 510.25
#> - state_pop             1     16166 1409095 510.49
#> - gdp                   1     36125 1429054 511.15
#> - m_per1000f            1     36467 1429396 511.16
#> <none>                              1392929 511.95
#> - unemploy_m24          1     86999 1479928 512.80
#> - prob_prison           1    171381 1564310 515.40
#> - unemploy_m39          1    196372 1589301 516.15
#> - percent_m             1    206121 1599050 516.43
#> - mean_education        1    371159 1764088 521.05
#> - inequality            1    534611 1927540 525.22
#> - police_exp60          1    728570 2121499 529.72
#> 
#> Step:  AIC=509.99
#> crime_rate ~ percent_m + mean_education + police_exp60 + labour_participation + 
#>     m_per1000f + state_pop + nonwhites_per1000 + unemploy_m24 + 
#>     unemploy_m39 + gdp + inequality + prob_prison
#> 
#>                        Df Sum of Sq     RSS    AIC
#> - labour_participation  1      3019 1397266 508.09
#> - nonwhites_per1000     1      7996 1402243 508.26
#> - state_pop             1     19634 1413881 508.65
#> - gdp                   1     35276 1429524 509.17
#> - m_per1000f            1     40680 1434928 509.34
#> <none>                              1394247 509.99
#> - unemploy_m24          1     85946 1480194 510.80
#> - unemploy_m39          1    195095 1589343 514.15
#> - percent_m             1    206909 1601157 514.50
#> - prob_prison           1    223309 1617557 514.98
#> - mean_education        1    381593 1775840 519.36
#> - inequality            1    537046 1931294 523.31
#> - police_exp60          1    764536 2158784 528.54
#> 
#> Step:  AIC=508.09
#> crime_rate ~ percent_m + mean_education + police_exp60 + m_per1000f + 
#>     state_pop + nonwhites_per1000 + unemploy_m24 + unemploy_m39 + 
#>     gdp + inequality + prob_prison
#> 
#>                     Df Sum of Sq     RSS    AIC
#> - nonwhites_per1000  1      6963 1404229 506.33
#> - state_pop          1     23381 1420648 506.87
#> - gdp                1     34787 1432053 507.25
#> - m_per1000f         1     41289 1438555 507.46
#> <none>                           1397266 508.09
#> - unemploy_m24       1     84385 1481652 508.85
#> - unemploy_m39       1    197849 1595115 512.32
#> - prob_prison        1    221812 1619078 513.02
#> - percent_m          1    226201 1623468 513.15
#> - mean_education     1    395884 1793150 517.82
#> - inequality         1    534370 1931637 521.32
#> - police_exp60       1    834362 2231628 528.10
#> 
#> Step:  AIC=506.33
#> crime_rate ~ percent_m + mean_education + police_exp60 + m_per1000f + 
#>     state_pop + unemploy_m24 + unemploy_m39 + gdp + inequality + 
#>     prob_prison
#> 
#>                  Df Sum of Sq     RSS    AIC
#> - state_pop       1     22345 1426575 505.07
#> - gdp             1     32142 1436371 505.39
#> - m_per1000f      1     36808 1441037 505.54
#> <none>                        1404229 506.33
#> - unemploy_m24    1     86373 1490602 507.13
#> - unemploy_m39    1    205814 1610043 510.76
#> - prob_prison     1    218607 1622836 511.13
#> - percent_m       1    307001 1711230 513.62
#> - mean_education  1    389502 1793731 515.83
#> - inequality      1    608627 2012856 521.25
#> - police_exp60    1   1050202 2454432 530.57
#> 
#> Step:  AIC=505.07
#> crime_rate ~ percent_m + mean_education + police_exp60 + m_per1000f + 
#>     unemploy_m24 + unemploy_m39 + gdp + inequality + prob_prison
#> 
#>                  Df Sum of Sq     RSS    AIC
#> - gdp             1     26493 1453068 503.93
#> <none>                        1426575 505.07
#> - m_per1000f      1     84491 1511065 505.77
#> - unemploy_m24    1     99463 1526037 506.24
#> - prob_prison     1    198571 1625145 509.20
#> - unemploy_m39    1    208880 1635455 509.49
#> - percent_m       1    320926 1747501 512.61
#> - mean_education  1    386773 1813348 514.35
#> - inequality      1    594779 2021354 519.45
#> - police_exp60    1   1127277 2553852 530.44
#> 
#> Step:  AIC=503.93
#> crime_rate ~ percent_m + mean_education + police_exp60 + m_per1000f + 
#>     unemploy_m24 + unemploy_m39 + inequality + prob_prison
#> 
#>                  Df Sum of Sq     RSS    AIC
#> <none>                        1453068 503.93
#> - m_per1000f      1    103159 1556227 505.16
#> - unemploy_m24    1    127044 1580112 505.87
#> - prob_prison     1    247978 1701046 509.34
#> - unemploy_m39    1    255443 1708511 509.55
#> - percent_m       1    296790 1749858 510.67
#> - mean_education  1    445788 1898855 514.51
#> - inequality      1    738244 2191312 521.24
#> - police_exp60    1   1672038 3125105 537.93
summary(model_crime_clean_step)
#> 
#> Call:
#> lm(formula = crime_rate ~ percent_m + mean_education + police_exp60 + 
#>     m_per1000f + unemploy_m24 + unemploy_m39 + inequality + prob_prison, 
#>     data = crime_clean)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -444.70 -111.07    3.03  122.15  483.30 
#> 
#> Coefficients:
#>                 Estimate Std. Error t value     Pr(>|t|)    
#> (Intercept)    -6426.101   1194.611  -5.379 0.0000040395 ***
#> percent_m          9.332      3.350   2.786      0.00828 ** 
#> mean_education    18.012      5.275   3.414      0.00153 ** 
#> police_exp60      10.265      1.552   6.613 0.0000000826 ***
#> m_per1000f         2.234      1.360   1.642      0.10874    
#> unemploy_m24      -6.087      3.339  -1.823      0.07622 .  
#> unemploy_m39      18.735      7.248   2.585      0.01371 *  
#> inequality         6.133      1.396   4.394 0.0000863344 ***
#> prob_prison    -3796.032   1490.646  -2.547      0.01505 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 195.5 on 38 degrees of freedom
#> Multiple R-squared:  0.7888, Adjusted R-squared:  0.7444 
#> F-statistic: 17.74 on 8 and 38 DF,  p-value: 0.0000000001159

Interpretasi model model_crime_clean_step :

  • Model step-wise ini berhenti dengan nilai AIC = 503.93 dengan menggunakan 8 variabel prediktor, yaitu :
    • m_per1000f, unemploy_m24, prob_prison, unemploy_m39, percent_m, mean_education, inequality, police_exp60.
  • Dari 8 variabel prediktor tersebut, saat digunakan menjadi prediktor untuk model regresi linear, hanya 6 variabel yang siginifikan.
  • Nilai R-squared untuk model yang dibuat menggunakan metode ini 0.7444.

5 Model Evaluation

Berdasarkan proses pembuatan model di atas, kita setidaknya memiliki 4 model regresi linear yang telah kita buat, yaitu :

  1. model_crime_simple : model linear dengan 1 variabel prediktor, yaitu menggunakan variabel police_exp60
  2. model_crime_multiple : model linear dengan beberapa variabel prediktor (police_exp60, gdp, dan prob_prison)
  3. model_crime_clean_all : model linear dengan seluruh variabel prediktor (selain variabel police_exp59)
  4. model_crime_clean_step : model linear yang dibangun dengan metode step-wise regression

5.1 R-Squared

Dengan kita melihat hasil dari R-squared (yang merepresentasikan Goodness of Fit dari suatu model), dari setiap model tersebut, model dengan banyak prediktor akan menghasilkan R-squared yang lebih baik. Berikut hasil R-squared dari tiap model tersebut.

# menggunakan Multiple R-squared karena hanya 1 variabel prediktor 
summary(model_crime_simple)$r.squared
#> [1] 0.4727999
# menggunakan Adjusted R-squared karena menggunakan > 1 variabel prediktor
summary(model_crime_multiple)$adj.r.squared
#> [1] 0.4991476
# menggunakan Adjusted R-squared karena menggunakan > 1 variabel prediktor
summary(model_crime_clean_all)$adj.r.squared
#> [1] 0.7090155
# menggunakan Adjusted R-squared karena menggunakan > 1 variabel prediktor
summary(model_crime_clean_step)$adj.r.squared
#> [1] 0.7443692

Berdasarkan evaluasi berdasarkan R-squared, model dengan metode step-wise adalah metode yang terbaik dibandingkan dengan seluruh model yang telah dibangun, dengan nilai 0.7444.

Dengan demikian untuk proses evaluasi model selanjutnya, kita cukup membandingkan 2 model dengan R-squared tertinggi yaitu model_crime_clean_all dan model_crime_clean_step.

5.2 Metrics Error

Untuk evaluasi performa model, kita akan menghitung nilai error dari hasil prediksi atas kedua model yang akan kita evaluasi. Tujuannya adalah kita akan melihat mana hasil prediksi model yang menghasilkan error paling kecil.

Sebelumnya kita memiliki kumpulan data baru yang mencatat variabel sosio-ekonomi yang sama dari pengamatan (observasi) yang berbeda. Data tersebut terdapat pada crime_test.csv, kemudian kita akan baca data tersebut dan simpan kedalam objek bernama crime_test.

crime_test <- read.csv("crime_test.csv")
head(crime_test)

Kemudian kita hitung nilai prediction dari kedua model tersebut dan kita simpan dalam setiap object.

predict_model_step <- predict(model_crime_clean_step, newdata = crime_test)
head(predict_model_step)
#>         1         2         3         4         5         6 
#>  745.0201 1012.3317  949.8039 1190.7017  772.4885  686.1097
crime_test <- crime_test %>% mutate(is_south = as.factor(is_south))
predict_model_all <- predict(model_crime_clean_all, newdata = crime_test)
head(predict_model_all)
#>         1         2         3         4         5         6 
#>  692.8391 1008.3053  922.3480 1200.5146  781.0148  699.6570

Kemudian kita hitung nilai error dari kedua prediksi model tersebut dan kita bandingkan hasilnya nanti.

5.2.1 Mean Absolute Error (MAE)

Mean Absolute Error (MAE) menunjukkan rata-rata dari nilai absolut error. MAE bisa diinterpretasikan sebagai seberapa besar penyimpangan hasil prediksi terhadap nilai aktualnya. Kita hitung masing-masing nilai MAE dari kedua model yang akan kita bandingkan.

  1. model_crime_clean_step
MAE(predict_model_step, crime_test$crime_rate)
#> [1] 180.7295
  1. model_crime_clean_all
MAE(predict_model_all, crime_test$crime_rate)
#> [1] 166.445

Berdasarkan nilai MAE tersebut di atas, model_crime_clean_all memiliki nilai error yang lebih kecil dibandingkan model model_crime_clean_step, sehingga dapat disimpulkan berdasarkan evaluasi error ini model_crime_clean_all lebih baik.

Namun untuk mengetahui apakah nilai MAE yang diperoleh cukup besar/tidak maka perlu dibandingkan dengan range dari target variabelnya.

range(crime_test$crime_rate)
#> [1]  373 1674

Karena nilai MAE relatif cukup kecil dibandingkan range data, maka model ini memiliki error yang cukup kecil.

5.2.2 Mean Squared Error (MSE)

Mean Squared Error (MSE) menunjukkan selisih kuadrat dari hasil prediksi dan nilai aktual kemudian dipangkatkan. Kita hitung masing-masing nilai MSE dari kedua model yang akan kita bandingkan.

  1. model_crime_clean_step
MSE(predict_model_step, crime_test$crime_rate)
#> [1] 46447.42
  1. model_crime_clean_all
MSE(predict_model_all, crime_test$crime_rate)
#> [1] 42551.24

5.2.3 Root Mean Squared Error (RMSE)

Root Mean Squared Error (RMSE) menunjukkan bentuk akar kuadrat dari MSE. Kita hitung masing-masing nilai RMSE dari kedua model yang akan kita bandingkan.

  1. model_crime_clean_step
RMSE(predict_model_step, crime_test$crime_rate)
#> [1] 215.5166
  1. model_crime_clean_all
RMSE(predict_model_all, crime_test$crime_rate)
#> [1] 206.2795

Berdasarkan nilai RMSE tersebut di atas, model_crime_clean_all memiliki nilai error yang lebih kecil dibandingkan model model_crime_clean_step, sehingga dapat disimpulkan berdasarkan evaluasi error ini model_crime_clean_all lebih baik.

5.2.4 Mean Absolute Percentage Error (MAPE)

Mean Absolute Percentage Error (MAPE) menunjukkan seberapa besar penyimpangannya dalam bentuk persentase. MAPE memiliki range dalam bentuk persen, semakin kecil nilai MAPE, semakin bagus model yang kita miliki. Kita hitung masing-masing nilai MAPE dari kedua model yang akan kita bandingkan.

  1. model_crime_clean_step
MAPE(predict_model_step, crime_test$crime_rate) * 100
#> [1] 23.39692
  1. model_crime_clean_all
MAPE(predict_model_all, crime_test$crime_rate) * 100
#> [1] 21.31936

Berdasarkan nilai MAPE tersebut di atas, model_crime_clean_all memiliki nilai error yang lebih kecil dibandingkan model model_crime_clean_step, sehingga dapat disimpulkan berdasarkan evaluasi error ini model_crime_clean_all lebih baik.

Kesimpulan atas nilai MAPE atas kedua model tersebut adalah error dari kedua model ini memiliki rata-rata menyimpang sebesar 21-24% dari data aktualnya.

Kesimpulan atas evaluasi model berdasarkan metrics error atas kedua model tersebut, maka model_crime_clean_all lebih baik daripada model_crime_clean_step karena memiliki nilai error yang sedikit lebih kecil.

6 Uji Asumsi

Setelah kita mendapatkan model evaluasi atas model yang telah kita buat, sekarang kita akan melakukan proses uji asumsi untuk memastikan apakah model yang kita buat dianggap sebagai Best Linear Unbiased Estimator (BLUE) model, yaitu model yang dapat memprediksi data baru secara konsisten.

Seperti halnya model evaluasi yang sudah kita lakukan, kita akan melakukan uji asumsi untuk model_crime_clean_all dan model_crime_clean_step.

6.1 Linearity (Linearitas)

Untuk menentukan asumsi linearitas dari model regresi linear dapat dilakukan dengan membuat plot residual vs fitted values. Plot ini merupakan scatter plot dengan sumbu x adalah nilai fitted values (hasil prediksi variabel respon) dan sumbu y adalah nilai residual.

  1. model_crime_clean_step
plot(model_crime_clean_step, which = 1)
abline(h = 200, col = "green")
abline(h = -200, col = "green")

  1. model_crime_clean_all
plot(model_crime_clean_all, which = 1)
abline(h = 200, col = "green")
abline(h = -200, col = "green")

Kesimpulan : kedua model tersebut adalah model linear (lulus uji asumsi linearity karena nilai residual berada di sekitaran nilai 0).

6.2 Normality

Model regresi linear diharapkan menghasilkan error yang terdistribusi secara normal. Dengan begitu, error lebih banyak berkumpul di sekitar angka nol. Untuk itu kita melakukan validasi asumsi normalitas pada kedua model ini menggunakan uji statistik dengan shapiro.test(). Fungsi ini membuat kita memasukkan nilai residual dari model kita.

  1. model_crime_clean_step
shapiro.test(model_crime_clean_step$residuals)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  model_crime_clean_step$residuals
#> W = 0.98511, p-value = 0.8051
  1. model_crime_clean_all
shapiro.test(model_crime_clean_all$residuals)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  model_crime_clean_all$residuals
#> W = 0.98889, p-value = 0.9308

Untuk pengujian normalitas ini, kita menggunakan Shapiro-Wilk hypothesis test:

H0: error berdistribusi normal H1: error TIDAK berdistribusi normal

Kondisi yang diharapkan: H0 - p-value > alpha -> gagal tolak h0 (terima h0) - p-value < alpha -> tolak h0 (terima h1)

Kesimpulan : Kedua model tersebut menghasilkan nilai p-value > 0.05 (lebih besar dari alpha). Hal ini menunjukkan kedua model ini memiliki error yang terdistribusi secara normal.

6.3 Homoscedasticity

Asumsi lain yang perlu kita uji adalah apakah error pada model terdistribusi dengan variasi yang sama/konstan di rentang data yang berbeda atau tidak.

Diharapkan error yang dihasilkan oleh model menyebar secara acak atau dengan variasi konstan. Apabila divisualisasikan maka error tidak berpola. Kondisi ini disebut juga sebagai homoscedasticity

Untuk menguji kondisi ini, kita bisa menggunakan fungsi bptest() dari package lmtest dan memasukkan kedua model yang ingin diuji.

  1. model_crime_clean_step
bptest(model_crime_clean_step)
#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  model_crime_clean_step
#> BP = 13.51, df = 8, p-value = 0.09546
  1. model_crime_clean_all
bptest(model_crime_clean_all)
#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  model_crime_clean_all
#> BP = 13.499, df = 14, p-value = 0.4876

Untuk pengujian ini, kita menggunakan Breusch-Pagan hypothesis test:

  • H0: error menyebar konstan atau homoscedasticity
  • H1: error menyebar TIDAK konstan atau heteroscedasticity

Kondisi yang diharapkan: H0

  • p-value > alpha -> gagal tolak h0 (terima h0)
  • p-value < alpha -> tolak h0 (terima h1)

Kesimpulan : Kedua model tersebut menghasilkan nilai p-value > 0.05 (lebih besar dari alpha). Hal ini menunjukkan kedua model ini memiliki residual menyebar secara random atau homoscedasticity.

6.4 Multicollinearity (Multikolinearitas)

Multikolinearitas adalah kondisi adanya korelasi antar prediktor yang kuat.

Hal ini tidak diinginkan karena menandakan prediktor redundan pada model, yang seharusnya dapat dipilih salah satu saja dari variabel yang hubungannya amat kuat tersebut. Harapannya tidak terjadi multikolinearitas.

Dengan menggunakan nilai Variance Inflation Factor (VIF), kita dapat menentukan ada tidaknya multikolinearitas antar variabel prediktor. Nilai VIF yang tinggi menunjukkan korelasi yang tinggi antar variabel prediktor. Anda dapat menggunakan fungsi vif dari package car dan memasukkan kedua model yang ingin diuji.

  1. model_crime_clean_step
vif(model_crime_clean_step)
#>      percent_m mean_education   police_exp60     m_per1000f   unemploy_m24 
#>       2.131963       4.189684       2.560496       1.932367       4.360038 
#>   unemploy_m39     inequality    prob_prison 
#>       4.508106       3.731074       1.381879
  1. model_crime_clean_all
vif(model_crime_clean_all)
#>            percent_m             is_south       mean_education 
#>             2.884547             5.317813             4.893078 
#>         police_exp60 labour_participation           m_per1000f 
#>             5.341037             3.421040             3.779003 
#>            state_pop    nonwhites_per1000         unemploy_m24 
#>             2.532206             4.277326             5.997380 
#>         unemploy_m39                  gdp           inequality 
#>             5.086769            10.496921             8.564507 
#>          prob_prison          time_prison 
#>             2.605654             2.376902

Uji VIF (Variance Inflation Factor) : * nilai VIF > 10: terjadi multikolinearitas pada model * nilai VIF < 10: tidak terjadi multikolinearitas pada model

Jadi, kondisi yang diharapkan adalah semua nilai variabel menghasilkan nilai vif < 10

Kesimpulan : Tidak terdapat Multikolinearitas pada kedua model tersebut karena seluruh variabel prediktor memiliki nilai VIF yang kurang dari 10

7 Conclusion

Model regresi linear ini dibangun untuk menjawab pertanyaan bisnis : nilai-nilai apa saja yang dapat mencerminkan tingkat kejahatan (tingkat kriminalitas) di suatu negara bagian?.

Tingkat kriminalitas pada data ini dicerminkan pada variabel crime_rate.

Berdasarkan hasil proses modeling, kita mendapatkan 2 model yang paling baik dari sisi evaluasi model untuk memprediksi crime_rate yaitu :

  1. model_crime_clean_all : model linear dengan seluruh variabel prediktor (selain variabel police_exp59)
  2. model_crime_clean_step : model linear yang dibangun dengan metode step-wise regression

Berdasarkan hasil evaluasi model, kita mendapatkan hasil evaluasi yang berbeda:

  1. Berdasarkan nilai R-squared, model yang lebih sedikit baik adalah model_crime_clean_step dengan nilai : 0.7444.
  2. Berdasarkan nilai metrics error (MAE, MSE, RMSE dan MAPE), model yang lebih baik adalah model_crime_clean_all dengan nilai error yang sedikit lebih baik.

Berdasarkan hasil uji asumsi, kedua model tersebut sama-sama lolos seluruh uji asumsi

Kesimpulan:

Karena selisih nilai R-squared maupun metrics error tidak berbeda jauh dan kedua model tersebut juga lolos uji asumsi model regresi linear, maka kedua model tersebut sama-sama layak untuk digunakan untuk memprediksi variabel-variabel apa saja yang mempengaruhi nilai crime_rate pada suatu negara bagian.