1 Introduction

1.1 Tujuan Pelatihan

Tujuan utama dari kursus ini adalah untuk memberikan pengenalan yang komprehensif untuk regression model dengan menggunakan data yang ada.

Regression Model adalah metode untuk menentukan hubungan suatu variable dengan yang lainnya untuk melihat seberapa besar pengaruhnya.

Fungsi dari Analisis Regresi:

  1. Memperbaiki kesalahan.

  2. Memperkirakan masa depan / melakukan prediksi.

  3. Meningkatkan efisiensi.

  4. Memberikan insight baru.

2 Explanatory Data Analysis

2.1 Import Library & Data

Load Library.

library(tidyverse)
library(caret)
library(plotly)
library(data.table)
library(GGally)
library(tidymodels)
library(car)
library(scales)
library(lmtest)
library(carData)
library(dplyr)
library(rsample)

options(scipen = 100, max.print = 1e+06)

Load Dataset.

crime <- read.csv("crime.csv", stringsAsFactors = F) %>% 
  dplyr::select(-X) # take out kolom X

head(crime)

Mengubah Nama variable

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)

Glossary data crime :

  • percent_m: percentage of males aged 14-24
  • is_south: whether it is in a Southern state. 1 for Yes, 0 for No.
  • mean_education: mean years of schooling
  • police_exp60: police expenditure in 1960
  • police_exp59: police expenditure in 1959
  • labour_participation: labour force participation rate
  • m_per1000f: number of males per 1000 females
  • state_pop: state population
  • nonwhites_per1000: number of non-whites resident per 1000 people
  • unemploy_m24: unemployment rate of urban males aged 14-24
  • unemploy_m39: unemployment rate of urban males aged 35-39
  • gdp: gross domestic product per head
  • inequality: income inequality
  • prob_prison: probability of imprisonment
  • time_prison: avg time served in prisons
  • crime_rate: crime rate in an unspecified category

2.2 Data Wrangling

Cek Tipe Data dan Sesuaikan Tipe Data

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 ...
crime$is_south <- as.factor(crime$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 ...

Dataset terdiri dari 47 obsevation (baris) dan 16 variables (kolom).

Pada kasus ini kita ingin memprediksi ‘inequality’ dari data ‘crime’:

  • Variable target: inequality

  • Variable predictor: percent_m, is_south, mean_education, dan lainnya.

Sebelum melakukan analisa lebih lanjut, eksplorasi data perlu dilakukan terlebih dahulu. Salah satu teknik eksplorasi data awal adalah menggunakan histogram atau boxplot.

hist(crime$inequality)

> Dari grafik histogram diatas, dapat dilihat tidak ada outlier pada variable target kita (inequality).

3 Simple Linear Regression

Model Regresi dengan 1 variable prediktor saja.

3.1 Modeling SLR

Kita akan menggunakan prediktor ‘gdp’ dan target ‘inequality’.

model_slr <- lm(formula = inequality ~ gdp,
                data = crime)
model_slr
#> 
#> Call:
#> lm(formula = inequality ~ gdp, data = crime)
#> 
#> Coefficients:
#> (Intercept)          gdp  
#>    386.0306      -0.3655

Dari model di atas, didapat bahwa intercept = 386.0306 dan koefisien gdp = -0.3655,

sehingga model_slr menghasilkan formula inequality = 386.0306 - 0.3655 * gdp

Dari formula di atas, apa kesimpulan yang bisa kita dapatkan?

  1. Model regresi yang dihasilkan dapat digunakan untuk memprediksi inequality

  2. Ketika gdp-nya 0, maka inequality-nya bernilai 386.0306 (sesuai dengan intercept).

  3. Ketika gdp bertambah satu satuan maka inequality berkurang -0.3655 (sesuai dengan coefficient).

  4. Makin besar nilai gdp, maka makin kecil nilai inequality.

Membuat plot untuk melihat sebaran data berdasarkan model persamaan regresi yang dihasilkan:

plot(crime$gdp, crime$inequality)
abline(model_slr, col = "red")

Salah satu metric untuk mengevaluasi model dari simple linear regression, dilihat multiple r.squarednya yang bisa diperoleh dari summary(model)

summary(model_slr)
#> 
#> Call:
#> lm(formula = inequality ~ gdp, data = crime)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -42.029 -11.754   1.714  12.438  30.006 
#> 
#> Coefficients:
#>              Estimate Std. Error t value            Pr(>|t|)    
#> (Intercept) 386.03058   15.38651   25.09 <0.0000000000000002 ***
#> gdp          -0.36551    0.02881  -12.69 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 18.86 on 45 degrees of freedom
#> Multiple R-squared:  0.7815, Adjusted R-squared:  0.7766 
#> F-statistic: 160.9 on 1 and 45 DF,  p-value: < 0.00000000000000022

Multiple R-squared = 78.15% artinya gdp (variable prediktor) mampu menjelaskan 78% varians/keragaman variable inequality, sisanya dijelaskan oleh variable lain yang tidak dimasukkan ke dalam model.

Multiple R-Squared = Metric yang melihat seberapa baik variable prediktor menangkap variansi dari variable target.

4 Multiple Linear Regression

Sebelum membuat model, kita perlu melakukan split data menjadi data train dan data test.

4.1 Data Correlation

Kita dapat melihat korelasi antar variable:

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

Terdapat beberapa variable yang berkorelasi kuat dengan variable target (inequality) kita, seperti:

  • Korelasi Positive: nonwhites_per1000, percent_m

  • Korelasi Negative: gdp, police_exp59, police_exp60, mean_education

4.2 Splitting Data Train dan Test

Tujuan dari splitting data adalah:

  • Data Train bertujuan untuk membuat model belajar, dan membuat prediksi.

  • Sedangkan Data Test untuk melakukan evaluasi antara model (prediksi) dengan data yang baru.

Dalam proses spliting data, best practice nya adalah 80:20.

RNGkind(sample.kind = "Rounding")
set.seed(998)

index <- sample(nrow(crime),
                nrow(crime)*0.8)

crime_train <- crime[index, ]
crime_test <- crime[-index,]

4.3 Modeling MLR

Model regresi linear dengan lebih dari satu predictor.

Kita akan membuat model regresi untuk memprediksi inequality dengan menggunakan seluruh prediktor yang ada.

set.seed(998)

model_mlr <- lm(formula = inequality ~ . , 
                data = crime_train)
summary(model_mlr)
#> 
#> Call:
#> lm(formula = inequality ~ ., data = crime_train)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -18.940  -5.928   1.800   7.656  16.147 
#> 
#> Coefficients:
#>                       Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)          429.86152  129.49935   3.319 0.003259 ** 
#> percent_m              0.07883    0.39628   0.199 0.844243    
#> is_south1             19.17891    9.37524   2.046 0.053515 .  
#> mean_education        -1.27230    0.45789  -2.779 0.011256 *  
#> police_exp60          -0.14766    0.79948  -0.185 0.855241    
#> police_exp59          -0.02652    0.86727  -0.031 0.975891    
#> labour_participation   0.26948    0.14466   1.863 0.076525 .  
#> m_per1000f            -0.17254    0.16160  -1.068 0.297786    
#> state_pop              0.01824    0.09424   0.194 0.848372    
#> nonwhites_per1000     -0.03046    0.04576  -0.666 0.512874    
#> unemploy_m24           0.23902    0.30116   0.794 0.436263    
#> unemploy_m39           0.09470    0.61172   0.155 0.878446    
#> gdp                   -0.25799    0.06012  -4.291 0.000324 ***
#> prob_prison          -13.53597  171.96277  -0.079 0.938005    
#> time_prison           -0.03239    0.55632  -0.058 0.954116    
#> crime_rate             0.03117    0.01088   2.864 0.009290 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 12.84 on 21 degrees of freedom
#> Multiple R-squared:  0.9385, Adjusted R-squared:  0.8945 
#> F-statistic: 21.36 on 15 and 21 DF,  p-value: 0.000000002112

Ketika ingin menilai suatu model yang memrtimbangkan banyak variable prediktor, gunakan adjusted r squared. Maka nilai dari Adj.R-squared model_mlr adalah 89.45%.

Intrepretasi

Berdasarkan p-value variable yang signifikan berpengaruh terhadap inequality adalah mean_education, gdp, dan crime_rate (yang memiliki tanda bintang).

Sehingga kita hanya akan menggunakan variable yang signifikan saja dan splitting data lagi dalam membuat model.

crime_new <- crime %>% 
  select(mean_education, gdp, crime_rate, inequality)

crime_train2 <- crime_new[index, ]
crime_test2 <- crime_new[-index, ]

set.seed(998)
model_mlr2 <- lm(inequality ~ ., 
                 data = crime_train2)

summary(model_mlr2)
#> 
#> Call:
#> lm(formula = inequality ~ ., data = crime_train2)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -23.5210  -7.5556   0.8991   7.8478  25.7551 
#> 
#> Coefficients:
#>                  Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)    454.124977  20.319590  22.349 < 0.0000000000000002 ***
#> mean_education  -1.128056   0.280638  -4.020             0.000318 ***
#> gdp             -0.316198   0.033920  -9.322      0.0000000000912 ***
#> crime_rate       0.028332   0.005798   4.887      0.0000257518155 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 12.66 on 33 degrees of freedom
#> Multiple R-squared:  0.9061, Adjusted R-squared:  0.8975 
#> F-statistic: 106.1 on 3 and 33 DF,  p-value: < 0.00000000000000022

Kami telah menghapus variabel yang tidak signifikan. Untuk melihat apakah tindakan ini memengaruhi model kita, kita dapat memeriksa nilai Adjusted R-Squared dari dua model tersebut.

Model pertama dengan variabel lengkap nilai Adjusted R-Squared sebesar 0.8945, artinya model tersebut dapat menjelaskan 89.45% varians pada variabel target (inequality).

Sementara itu, model kami yang lebih sederhana Adjusted R-Squared sebesar 89.75%, tidak ada perbedaan besar dengan model sederhana kami. Hal ini menunjukkan bahwa aman untuk menghilangkan variabel yang tidak memiliki nilai koefisien yang signifikan.

5 Evaluation

Tujuan melakukan pemodelan yang bersifat prediktif adalah memprediksi data yang belum dimiliki. Untuk mengetahui model machine learning yang telah dibuat sudah bagus atau tidak dapat melihat error yang dihasilkan. Ketika membuat model, diharapkan error yang dihasilkan sekecil mungkin.

5.1 Model Performance

Menggunakan Root Mean Square Error (RMSE)

Model 1 (model_mlr) performance

crime_pred <- predict(model_mlr,
                      newdata = crime_test %>% select(-inequality))

# RMSE data train model 1
RMSE(pred = model_mlr$fitted.values,
     obs = crime_train$inequality)
#> [1] 9.676469
# RMSE data test model 1
RMSE(pred = crime_pred,
     obs = crime_test$inequality)
#> [1] 21.55809

Model 2 (model_mlr2 / sederhana) performance

crime_pred2 <- predict(model_mlr2,
                      newdata = crime_test2 %>% select(-inequality))

# RMSE data train model 1
RMSE(pred = model_mlr2$fitted.values,
     obs = crime_train2$inequality)
#> [1] 11.95743
# RMSE data test model 1
RMSE(pred = crime_pred2,
     obs = crime_test2$inequality)
#> [1] 21.44447

Model kedua lebih baik dari yang pertama dalam memprediksi dataset test, bahkan hanya dengan margin kecil. Pada kedua model, kesalahan set data train lebih rendah daripada set data test, menunjukkan bahwa model kami mungkin overfit.

5.2 Asumsi

5.2.1 Normality of Residual

Harapannya ketika membuat model linear regression, error/residul yang dihasilkan berdistribusi normal. Artinya error banyak berkumpul disekitar angka 0. Untuk menguji asumsi ini dapat dilakukan:

  1. Visualisasi histogram residual, dengan menggunakan fungsi hist().
hist(model_mlr2$residuals)

> Kesimpulan: residual/error pada hasil prediksi sudah berdistribusi normal, yang mana artinya dari hasil prediksi yang diperoleh model lumayan akurat karena errornya banyak yang mengumpul di titik 0.

  1. Uji statistik menggunakan shapiro.test(). (harapannya pvalue > alpha agar keputusan yang diambil adalah gagal tolak H0)

H0: error berdistribusi normal H1: error tidak berdistribusi normal

harapannya p-value > 0.05

Shapiro-Wilk hypothesis test:

shapiro.test(model_mlr2$residuals)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  model_mlr2$residuals
#> W = 0.97315, p-value = 0.5001

Kesimpulan: p-value nya 0.5001 dengan kata lain H0 diterima / error berdistibusi normal.

5.2.2 Homoscedasticity of Residual

Untuk menguji asumsi ini dapat dilakukan:

  1. Visualisasi dengan scatterplot antara nilai prediksi(fitted values) dengan nilai error
# plot model.fitted dengan residuals
plot(x = model_mlr2$fitted.values, 
     y = model_mlr2$residuals)

abline(h = 0, col = "red", lty = 2)

> Kesimpulan: Variansi error nya menyebar konstan (Homoscedasticity).

  1. Uji statistika dengan Breusch-Pagan dari package lmtest

Breusch-Pagan hypothesis test: (harapannya p-value > alpha agar gagal tolak H0)

bptest(model_mlr2)
#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  model_mlr2
#> BP = 4.8059, df = 3, p-value = 0.1866

Kesimpulan: p-value nya 0.1866 dengan kata lain H0 diterima / variansi error menyebar konstan (Homoscedasticity).

5.2.3 No Multicolinearity

Harapannya pada model linear regression, tidak terjadi multikolinearitas. Multikolinearitas terjadi ketika antar variabel prediktor yang digunakan pada model memiliki hubungan yang kuat. Ada atau tidak multikolinearitas dapat dilihat dari nilai VIF(Variance Inflation Factor):

Ketika nilai VIF lebih dari 10 artinya terjadi multikolinearitas. Harapannya mendapatkan VIF < 10

vif(model_mlr2)
#> mean_education            gdp     crime_rate 
#>       2.296194       2.524987       1.295814

Kesimpulan: Model sudah bagus karena no Multicolinearity.

6 Kesimpulan / Conclusion

Variable yang berguna untuk menggambarkan variansi inequality adalah mean_education, gdp, dan crime_rate.

Model sederhana yang dibuatkan juga sudah memenuhi asumsi (classical assumptions).

Nilai Adjusted R-Squared sebesar 0.8975, artinya model tersebut dapat menjelaskan 89.75% varians pada variabel target (inequality).

Keakuratan model dalam melakukan prediksi inequality diukur dengan RMSE : RMSE train= 11.95743, dan RMSE test= 21.44447 menunjukkan bahwa model tersebut overfit.