PENDAHULUAN

Kita akan belajar menggunakan model regresi linear, menggunakan dataset crime.csv. Kita ingin mengetahui hubungan antar variabel, khususnya variabel crime_rate dengan variabel lainnya.

PERSIAPAN DATA

Gunakan library dibawah ini :

library(tidyverse)
library(caret)
library(plotly)
library(data.table)
library(GGally)
library(car)
library(scales)
library(lmtest)

Load dataset crime.csv

crime <- read.csv("data_input/crime.csv")
crime

Cek apakah terdapat missing value

colSums(is.na(crime))
##    X    M   So   Ed  Po1  Po2   LF  M.F  Pop   NW   U1   U2  GDP Ineq Prob Time 
##    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0 
##    y 
##    0

Sesuaikan nama kolom dari dataset crime.csv

crime <- read.csv("data_input/crime.csv")
colnames(crime) <- c("no", "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")
crime

Cek tipe data apakah sudah sesuai atau belum

str(crime)
## 'data.frame':    47 obs. of  17 variables:
##  $ no                  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ 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 ...

EXPLORATORY DATA ANALYSIS

Cek kolerasi antar variabel

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

Dari grafik diatas, dapat dilihat banyak variabel yang memiliki korelasi dengan variabel crime_rate.

MODELING

Train-Test Split

Sebelum kita membuat model, kita perlu membagi data menjadi dataset train dan dataset test. Kita akan menggunakan dataset train untuk melatih model regresi linier. Dataset test akan digunakan sebagai pembanding dan melihat apakah model menjadi overfit dan tidak dapat memprediksi data baru yang belum terlihat selama fase training. Kita akan menggunakan 70% dari data sebagai data train dan sisanya sebagai data test.

set.seed(123)
samplesize <- round(0.7 * nrow(crime), 0)
index <- sample(seq_len(nrow(crime)), size = samplesize)

data_train <- crime[index, ]
data_test <- crime[-index, ]

Regresi Linear

set.seed(123)
crime_lm <- lm(crime_rate ~ ., data = data_train)

summary(crime_lm)
## 
## Call:
## lm(formula = crime_rate ~ ., data = data_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -252.38  -66.57  -31.65   76.80  399.58 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -9727.4733  2265.9745  -4.293 0.000559 ***
## no                      -3.6008     3.0504  -1.180 0.255084    
## percent_m               10.5988     4.4936   2.359 0.031395 *  
## is_south               -19.1557   221.9756  -0.086 0.932301    
## mean_education          23.8519     7.8869   3.024 0.008061 ** 
## police_exp60            21.2341    11.7678   1.804 0.090020 .  
## police_exp59           -16.1373    13.4265  -1.202 0.246891    
## labour_participation    -4.5126     2.1510  -2.098 0.052150 .  
## m_per1000f               7.7615     2.8602   2.714 0.015337 *  
## state_pop                1.0530     1.7923   0.587 0.565078    
## nonwhites_per1000        1.2015     0.9941   1.209 0.244335    
## unemploy_m24           -11.3288     5.1038  -2.220 0.041241 *  
## unemploy_m39            14.5241     9.7595   1.488 0.156144    
## gdp                      0.7919     1.2543   0.631 0.536720    
## inequality               4.3596     2.8243   1.544 0.142234    
## prob_prison            -40.5129  3869.7654  -0.010 0.991776    
## time_prison             12.3693     9.5434   1.296 0.213325    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 197.6 on 16 degrees of freedom
## Multiple R-squared:  0.8908, Adjusted R-squared:  0.7816 
## F-statistic: 8.158 on 16 and 16 DF,  p-value: 6.385e-05

Dari summary(crime_lm) menunjukkan banyak informasi. Tapi untuk saat ini, lebih baik fokus pada Pr(>|t|). Kolom ini menunjukkan tingkat signifikansi variabel terhadap model. Jika nilainya < 0.05, maka kita dapat dengan aman mengasumsikan bahwa variabel tersebut berpengaruh signifikan terhadap model (artinya koefisien yang diestimasi tidak berbeda dengan 0), begitu pula sebaliknya.

Dengan demikian, kita dapat membuat model yang lebih sederhana dengan cara menghapus variabel yang memiliki nilai Pr(>|t|) > 0.05, karena variabel tersebut tidak berpengaruh signifikan terhadap model kita.

crime2 <- crime %>%
  select(percent_m, mean_education, m_per1000f, unemploy_m24, crime_rate)

data_train2 <- crime2[index, ]
data_test2 <- crime2[-index, ]

set.seed(123)
crime_lm2 <- lm(crime_rate ~ ., data = data_train2)

summary(crime_lm2)
## 
## Call:
## lm(formula = crime_rate ~ ., data = data_train2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -576.85 -314.11  -37.03  157.29  926.35 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -2537.864   2619.161  -0.969    0.341
## percent_m          2.687      7.467   0.360    0.722
## mean_education    11.502      8.489   1.355    0.186
## m_per1000f         2.213      3.036   0.729    0.472
## unemploy_m24      -3.685      5.015  -0.735    0.469
## 
## Residual standard error: 416.6 on 28 degrees of freedom
## Multiple R-squared:  0.1503, Adjusted R-squared:  0.02889 
## F-statistic: 1.238 on 4 and 28 DF,  p-value: 0.3175

Kita telah menghapus variabel tidak signifikan. Untuk melihat apakah tindakan ini mempengaruhi model kita, kita dapat memeriksa nilai Adjusted R-Squared dari dua model sebelumnya.

Model pertama dengan variabel lengkap Adjusted R-squared sebesar 0.7816 (78.16%), artinya model dapat menjelaskan 78.16% varians pada variabel target.

Sementara itu, model kita yang lebih sederhana Adjusted R-squared sebesar 0.02889 atau 2.889%, artinya model hanya dapat menjelaskan 2.889% varians pada variabel target. Ini menunjukkan bahwa sangat signifikan menghapus variabel yang tidak kita butuhkan.

EVALUASI

Model Performance

Performance model pertama (dengan variabel lengkap)

lm_pred <- predict(crime_lm, newdata = data_test %>% select(-crime_rate))

#RMSE of train dataset
RMSE(pred = crime_lm$fitted.values, obs = data_train$crime_rate)
## [1] 137.5803
# RMSE of test dataset
RMSE(pred = lm_pred, obs = data_test$crime_rate)
## [1] 385.2774

Performance model kedua (dengan variabel yang lebih sederhana)

lm_pred2 <- predict(crime_lm2, newdata = data_test2 %>% select(-crime_rate))

# RMSE of train dataset
RMSE(pred = crime_lm2$fitted.values, obs = data_train2$crime_rate)
## [1] 383.7879
# RMSE of test dataset
RMSE(pred = lm_pred2, obs = data_test2$crime_rate)
## [1] 305.5912

Model kedua lebih baik dari yang pertama dalam memprediksi data test, sedangkan model pertama lebih baik dalam memprediksi data train, terlihat dari RMSE data train model pertama jauh lebih kecil dibandingkan RMSE data train model yang kedua.

Asumsi

Regresi linier merupakan model parametrik, artinya untuk membuat persamaan model, model tersebut mengikuti beberapa asumsi klasik. Regresi linier yang tidak mengikuti asumsi dapat menyesatkan, atau hanya memiliki penduga yang bias. Untuk bagian ini, kita hanya memeriksa model kedua (model dengan variabel yang dihapus).

1. Linearity

resact <- data.frame(residual = crime_lm2$residuals, fitted = crime_lm2$fitted.values)

resact %>% ggplot(aes(fitted, residual)) + 
  geom_point() + 
  geom_smooth() + 
  geom_hline(aes(yintercept = 0)) + 
    theme(panel.grid = element_blank(), panel.background = element_blank())
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Pola diatas menunjukkan bahwa model kita mungkin tidak cukup linier.

2. Normality Test

shapiro.test(crime_lm2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  crime_lm2$residuals
## W = 0.93145, p-value = 0.03848

Hipotesis nol adalah residual mengikuti distribusi normal. Dengan p-value > 0.05, kita dapat menyimpulkan bahwa hipotesis nol kita terima, dan residual kita mengikuti distribusi normal.

3. Autocorrelation

Autocorrelation dapat dideteksi dengan menggunakan uji durbin watson test, dengan hipotesis nol yaitu tidak terdapat Autocorrelation.

durbinWatsonTest(crime_lm2)
##  lag Autocorrelation D-W Statistic p-value
##    1      -0.2932646      2.520759   0.106
##  Alternative hypothesis: rho != 0

Dengan p-value > 0.05 maka dapat disimpulkan bahwa tidak terdapat Autocorrelation.

4. Heterocedasticity

bptest(crime_lm2)
## 
##  studentized Breusch-Pagan test
## 
## data:  crime_lm2
## BP = 7.7235, df = 4, p-value = 0.1023
resact %>% 
  ggplot(aes(fitted, residual)) + 
  geom_point() + 
  theme_light() + 
  geom_hline(aes(yintercept = 0))

Dengan p-value > 0.05, kita dapat menyimpulkan bahwa tidak terdapat Heterocedasticity dalam model kita.

5. Multicollinearity

vif(crime_lm2)
##      percent_m mean_education     m_per1000f   unemploy_m24 
##       1.449995       1.680526       1.462630       1.313499

Dari data yang ditunjukkan diatas dapat disimpulkan bahwa tidak ada Multikolinearitas (karena tidak > 10).

Kita telah melihat bahwa model kita telah memenuhi beberapa asumsi, termasuk Linearitas, Heterosdastisitas, Autokorelasi, dan Multikolinearitas. Maka dari itu kita tidak perlu melakukan Perbaikan Model.

summary(crime_lm2)
## 
## Call:
## lm(formula = crime_rate ~ ., data = data_train2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -576.85 -314.11  -37.03  157.29  926.35 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -2537.864   2619.161  -0.969    0.341
## percent_m          2.687      7.467   0.360    0.722
## mean_education    11.502      8.489   1.355    0.186
## m_per1000f         2.213      3.036   0.729    0.472
## unemploy_m24      -3.685      5.015  -0.735    0.469
## 
## Residual standard error: 416.6 on 28 degrees of freedom
## Multiple R-squared:  0.1503, Adjusted R-squared:  0.02889 
## F-statistic: 1.238 on 4 and 28 DF,  p-value: 0.3175

KESIMPULAN

Variabel yang berguna untuk mendeskripsikan varians crime_rate adalah state_pop, dan time_prison. Model terakhir kita telah memenuhi asumsi. Meskipun Adjusted R-squared modelnya hanya 0.02889 yang artinya variabel hanya dapat menjelaskan 2.889% varian dari target variabel.

Keakuratan model dalam memprediksi crime_rate diukur dengan RMSE, dengan data train memiliki RMSE 383.7879 dan data test memiliki RMSE 305.5912, menunjukkan bahwa model kami mungkin overfit pada dataset traning.

Kita telah mempelajari bagaimana membangun model regresi linier dan apa yang perlu diperhatikan saat membangun model. Untuk mendapatkan model yang bagus, data crime.csv mungkin perlu ditambahkan lagi, sehingga kita dapat menghasilkan model regresi linear yang outvit.