Perbandingan Performa Model : Regresi Klasik vs Ridge vs Lasso

Pada kesempatan ini akan dilakukan perbandingan performa model dari beberapa jenis model regresi yakni regresi klasik, regresi ridge dan regresi lasso. Dasar perbandingan ini ingin mengetahui manakah yang lebih baik dalam menjabarkan permasalahan variable selection. Pembandingan dilakukan dengan melihat nilai R-square pada setiap modelnya.

Data

Data yang digunakan merupakan data mengenai indeks baca dan faktor-faktor yang diduga berpengaruh kepadanya. Data ini merupakan data cross-section dengan amatan Provinsi di Indonesia.

Input data

library(readr)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

data<- rio::import("https://raw.githubusercontent.com/fax17/PSD/main/Regresi%20Ridge%2C%20Lasso%20dan%20Klasik/datap7.csv")
data

##    Indeks_Alibaca   IPM Universitas_Negeri Universitas_Swasta Univ
## 1           34,37  71,9                  7                 89   96
## 2           35,73 71,74                  3                211  214
## 3           38,57 72,39                  5                 87   92
## 4           38,71    73                  2                 81   83
## 5           37,32 71,26                  1                 40   41
## 6           36,06 70,02                  2                 99  101
## 7           37,41 71,21                  2                 15   17
## 8           30,59 69,57                  3                 71   74
## 9           41,97  71,3                  2                 14   16
## 10          54,76 75,48                  2                 32   34
## 11          58,16 80,76                  4                272  276
## 12          39,47 72,03                 12                376  388
## 13           33,3 71,73                  9                241  250
## 14           56,2 79,99                  5                104  109
## 15          33,19  71,5                 17                324  341
## 16          40,81 72,44                  2                111  113
## 17          44,58 75,38                  4                 56   60
## 18          33,64 68,14                  1                 54   55
## 19          29,83 65,23                  4                 61   65
## 20          28,63 67,65                  4                 46   50
## 21          33,86 70,91                  1                 23   24
## 22             37 70,72                  3                 46   49
## 23          46,01 76,61                  5                 49   54
## 24          42,86 71,15                  2                  8   10
## 25           40,2 72,99                  4                 47   51
## 26          31,55  69,5                  1                 36   37
## 27          38,82 71,66                  5                187  192
## 28          34,37  71,2                  2                 51   53
## 29          34,99 68,49                  1                 11   12
## 30          32,92 65,73                  1                 19   20
## 31          33,52 69,45                  3                 30   33
## 32          31,33  68,7                  1                 18   19
## 33          28,25  64,7                  2                 20   22
## 34           19,9 60,84                  3                 53   56
##    Jumlah_Perpustakaan
## 1                 5442
## 2                11629
## 3                 5732
## 4                 1900
## 5                 4276
## 6                 3995
## 7                 1317
## 8                 5412
## 9                 1307
## 10                 774
## 11                2103
## 12               22122
## 13               22028
## 14                3365
## 15               17450
## 16                7808
## 17                2163
## 18                5766
## 19                4813
## 20                3279
## 21                1623
## 22                3843
## 23                3025
## 24                 236
## 25                3000
## 26                2032
## 27                9037
## 28                2286
## 29                1224
## 30                1614
## 31                2411
## 32                 700
## 33                 406
## 34                 322

Data preprocessing

str(data)

## 'data.frame':    34 obs. of  6 variables:
##  $ Indeks_Alibaca     : chr  "34,37" "35,73" "38,57" "38,71" ...
##  $ IPM                : chr  "71,9" "71,74" "72,39" "73" ...
##  $ Universitas_Negeri : int  7 3 5 2 1 2 2 3 2 2 ...
##  $ Universitas_Swasta : int  89 211 87 81 40 99 15 71 14 32 ...
##  $ Univ               : int  96 214 92 83 41 101 17 74 16 34 ...
##  $ Jumlah_Perpustakaan: int  5442 11629 5732 1900 4276 3995 1317 5412 1307 774 ...

#ubah data yang strukturnya masih tidak tepat
data <- data %>%
  mutate(Indeks_Alibaca = as.numeric(gsub(",", ".", Indeks_Alibaca)))
data <- data %>%
  mutate(IPM = as.numeric(gsub(",", ".", IPM)))

note: disini data peubah univ tidak digunakan karena hanya berupa penjumlahan univ negeri dan swastanya.

Pembentukan matriks untuk pemodelan

Y <- data.matrix(data[, 1])
x1 <- data.matrix(data[, 2]) 
x2 <- data.matrix(data[, 3])
x3 <- data.matrix(data[, 4])
x4 <- data.matrix(data[, 6])

X <- cbind(x1,x2,x3,x4)

Dalam pembentukan matriks disini ditetapkan kode peubah sebagai berikut

Y : indeks baca (poin)

x1 : IPM (poin)

x2 : Jumlah Universitas Negeri

x3 : Jumlah Universitas Swasta

x4 : Jumlah Perpustakaan

Model Regresi Klasik

Model awal

Model_klasik <- lm(Y~x1+x2+x3+x4)
summary(Model_klasik)

## 
## Call:
## lm(formula = Y ~ x1 + x2 + x3 + x4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0088 -2.0779 -0.6022  1.6091  8.4031 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -9.118e+01  1.092e+01  -8.351 3.33e-09 ***
## x1           1.826e+00  1.558e-01  11.724 1.59e-12 ***
## x2          -2.367e-01  2.834e-01  -0.835    0.410    
## x3           1.298e-02  1.304e-02   0.996    0.328    
## x4          -3.090e-04  2.065e-04  -1.496    0.145    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.102 on 29 degrees of freedom
## Multiple R-squared:  0.8644, Adjusted R-squared:  0.8457 
## F-statistic: 46.21 on 4 and 29 DF,  p-value: 3.544e-12

Multikolinearitas

# Menghitung VIF
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

vif_results <- vif(Model_klasik)
print(vif_results)

##       x1       x2       x3       x4 
## 1.274154 3.099745 5.077081 4.604487

Beberapa sumber mengatakan bahwa nilai VIF> 5 memiliki potensi untuk terjadi multikolinearitas. Sehingga akan dilakukan pemodelan regresi dengan ridge dan lasso sebagai pembanding dari model regresi klasik yang telah dibentuk ini yang akan dilakukan seleksi peubahnya jika Vif>5 tersebut.

Model tanpa x3

modelnew <- lm(Y ~ x1+x2+x4)
car::vif(modelnew)

##       x1       x2       x4 
## 1.055117 2.748074 2.658884

Terlihat bahwa sudah tidak terdapat multikolinearitas lagi pada peubah.

Pengujian Asumsi

Uji Asumsi Normalitas

ks.test(modelnew$residuals, "pnorm", mean=mean(modelnew$residuals), sd=sd(modelnew$residuals))

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  modelnew$residuals
## D = 0.14806, p-value = 0.4058
## alternative hypothesis: two-sided

shapiro.test(modelnew$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelnew$residuals
## W = 0.95164, p-value = 0.1375

Berdasarkan Kolmogorov-Smirnov test dan Shapiro-Wilk, residual data menyebar normal dengan p-value > 5%.

Uji Asumsi Homoskedastisitas (Gauss Markov)

lmtest::bptest(modelnew)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelnew
## BP = 0.90202, df = 3, p-value = 0.8249

Karena p-value > 0.05 maka ragam sisaan homogen atau tidak terdapat masalah heteroskedastisitas

randtests::runs.test(modelnew$residuals)

## 
##  Runs Test
## 
## data:  modelnew$residuals
## statistic = 0.69663, runs = 20, n1 = 17, n2 = 17, n = 34, p-value =
## 0.486
## alternative hypothesis: nonrandomness

Karena p-value > 0.05 maka sisaan saling bebas.

Nilai harapan sisaan sama dengan nol (Gauss Markov)

t.test(modelnew$residuals,
       mu = 0,
       conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  modelnew$residuals
## t = -1.0313e-16, df = 33, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.031802  1.031802
## sample estimates:
##     mean of x 
## -5.230319e-17

Karena p-value > 0.05 maka nilai harapan sisaan sama dengan nol

Pemilihan Peubah Penjelas/ Variable Selection

Metode Backward

modelbackward <- step(modelnew, direction="backward", scope=formula(lm(Y ~ x1+x2+x4)), trace=1)

## Start:  AIC=80.71
## Y ~ x1 + x2 + x4
## 
##        Df Sum of Sq     RSS     AIC
## - x2    1      2.71  291.29  79.031
## - x4    1     12.01  300.59  80.099
## <none>               288.58  80.713
## - x1    1   1711.65 2000.23 144.538
## 
## Step:  AIC=79.03
## Y ~ x1 + x4
## 
##        Df Sum of Sq     RSS     AIC
## <none>               291.29  79.031
## - x4    1     59.59  350.88  83.359
## - x1    1   1756.15 2047.45 143.332

summary(modelbackward)

## 
## Call:
## lm(formula = Y ~ x1 + x4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4162 -2.1314 -0.8061  1.7496  8.1357 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -9.473e+01  9.720e+00  -9.746 5.90e-11 ***
## x1           1.875e+00  1.372e-01  13.671 1.14e-14 ***
## x4          -2.410e-04  9.568e-05  -2.518   0.0172 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.065 on 31 degrees of freedom
## Multiple R-squared:  0.8584, Adjusted R-squared:  0.8493 
## F-statistic: 93.99 on 2 and 31 DF,  p-value: 6.914e-14

Berdasarkan metode backward, diperoleh aic 79.031 dan model terbaik adalah model dengan peubah x1 dan x4 saja dengan R-squared: 0.8584

Metode Forward

modelforward <- step(lm(Y ~ 1), direction="forward", scope=formula(modelnew), trace=1)

## Start:  AIC=141.5
## Y ~ 1
## 
##        Df Sum of Sq     RSS     AIC
## + x1    1   1706.81  350.88  83.359
## <none>              2057.69 141.501
## + x4    1     10.24 2047.45 143.332
## + x2    1      2.93 2054.76 143.453
## 
## Step:  AIC=83.36
## Y ~ x1
## 
##        Df Sum of Sq    RSS    AIC
## + x4    1    59.591 291.29 79.031
## + x2    1    50.293 300.59 80.099
## <none>              350.88 83.359
## 
## Step:  AIC=79.03
## Y ~ x1 + x4
## 
##        Df Sum of Sq    RSS    AIC
## <none>              291.29 79.031
## + x2    1    2.7133 288.58 80.713

summary(modelforward)

## 
## Call:
## lm(formula = Y ~ x1 + x4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4162 -2.1314 -0.8061  1.7496  8.1357 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -9.473e+01  9.720e+00  -9.746 5.90e-11 ***
## x1           1.875e+00  1.372e-01  13.671 1.14e-14 ***
## x4          -2.410e-04  9.568e-05  -2.518   0.0172 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.065 on 31 degrees of freedom
## Multiple R-squared:  0.8584, Adjusted R-squared:  0.8493 
## F-statistic: 93.99 on 2 and 31 DF,  p-value: 6.914e-14

Hasil forward selection menunjukkan hasil yang sama dengan backward baik AIC maupun peubah terpilihnya dengan R-squared: 0.8584

####Metode Stepwise

modelstepwise <- step(lm(Y ~ 1), direction="both", scope=formula(modelnew), trace=1)

## Start:  AIC=141.5
## Y ~ 1
## 
##        Df Sum of Sq     RSS     AIC
## + x1    1   1706.81  350.88  83.359
## <none>              2057.69 141.501
## + x4    1     10.24 2047.45 143.332
## + x2    1      2.93 2054.76 143.453
## 
## Step:  AIC=83.36
## Y ~ x1
## 
##        Df Sum of Sq     RSS     AIC
## + x4    1     59.59  291.29  79.031
## + x2    1     50.29  300.59  80.099
## <none>               350.88  83.359
## - x1    1   1706.81 2057.69 141.501
## 
## Step:  AIC=79.03
## Y ~ x1 + x4
## 
##        Df Sum of Sq     RSS     AIC
## <none>               291.29  79.031
## + x2    1      2.71  288.58  80.713
## - x4    1     59.59  350.88  83.359
## - x1    1   1756.15 2047.45 143.332

bestmodel <- summary(modelstepwise)
bestmodel

## 
## Call:
## lm(formula = Y ~ x1 + x4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4162 -2.1314 -0.8061  1.7496  8.1357 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -9.473e+01  9.720e+00  -9.746 5.90e-11 ***
## x1           1.875e+00  1.372e-01  13.671 1.14e-14 ***
## x4          -2.410e-04  9.568e-05  -2.518   0.0172 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.065 on 31 degrees of freedom
## Multiple R-squared:  0.8584, Adjusted R-squared:  0.8493 
## F-statistic: 93.99 on 2 and 31 DF,  p-value: 6.914e-14

Hasil stepwise menunjukkan hasil yang sama dengan metode backward dan forward selection.

R-Squared terbaik

Berdasarkan serangkaian metode regresi klasik dan juga berdasar ketiga model seleksi variabel maka diperoleh nilai terbaik R-square pada model klasik yakni:

R-squared: 0.8584

Ridge Regression

Pemodelan

modelRidge <- glmnet::cv.glmnet(X,Y,alpha=0)
summary(modelRidge)

##            Length Class  Mode     
## lambda     100    -none- numeric  
## cvm        100    -none- numeric  
## cvsd       100    -none- numeric  
## cvup       100    -none- numeric  
## cvlo       100    -none- numeric  
## nzero      100    -none- numeric  
## call         4    -none- call     
## name         1    -none- character
## glmnet.fit  12    elnet  list     
## lambda.min   1    -none- numeric  
## lambda.1se   1    -none- numeric  
## index        2    -none- numeric

print(modelRidge)

## 
## Call:  glmnet::cv.glmnet(x = X, y = Y, alpha = 0) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure    SE Nonzero
## min 0.7085   100   11.65 2.526       4
## 1se 2.3747    87   13.96 3.334       4

Lambda Terbaik

bestlambda <- modelRidge$lambda.min
cat("Lambda Terbaik:", bestlambda, "\n")

## Lambda Terbaik: 0.7085208

koefisien model Ridge

coef(modelRidge, s = bestlambda)

## 5 x 1 sparse Matrix of class "dgCMatrix"
##                        s1
## (Intercept) -8.098394e+01
## V1           1.680118e+00
## V2          -1.877712e-01
## V3           9.912226e-03
## V4          -2.545575e-04

R-squared model Ridge

Yduga <- predict(modelRidge, s = bestlambda, newx = X)
r_squared_ridge <- 1 - sum((Y - Yduga)^2) / sum((Y - mean(Y))^2)
r_squared_ridge

## [1] 0.8577269

Pada regresi Ridge, tidak terdapat peubah yang dihilangkan atau semua peubah dimasukkan dalam model. R-Square yang diperoleh dengan regresi ridge cukup baik yakni:

R−squared: 0.8577269 .

Regresi Lasso

Pemodelan

modelLasso <- glmnet::cv.glmnet(X,Y,alpha=1)
modelLasso

## 
## Call:  glmnet::cv.glmnet(x = X, y = Y, alpha = 1) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure    SE Nonzero
## min 0.0293    60   10.69 2.831       4
## 1se 1.2097    20   13.34 3.238       1

print(modelLasso)

## 
## Call:  glmnet::cv.glmnet(x = X, y = Y, alpha = 1) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure    SE Nonzero
## min 0.0293    60   10.69 2.831       4
## 1se 1.2097    20   13.34 3.238       1

Lambda Terbaik

bestlambdalasso <- modelLasso$lambda.min
cat("Lambda terbaik:", bestlambdalasso, "\n")

## Lambda terbaik: 0.02927617

Koefisien Model Lasso

coef(modelLasso, s = bestlambdalasso)

## 5 x 1 sparse Matrix of class "dgCMatrix"
##                        s1
## (Intercept) -9.149006e+01
## V1           1.830302e+00
## V2          -2.077080e-01
## V3           1.023356e-02
## V4          -2.789346e-04

R-squared model Lasso

Ydugalasso <- predict(modelLasso, s = bestlambdalasso, newx = X)
r_squared_lasso <- 1 - sum((Y - Ydugalasso)^2) / sum((Y - mean(Y))^2)
r_squared_lasso

## [1] 0.8641478

Melalui regresi Lasso terlihat tidak ada juga peubah yang dikeluarkan dari model. R-Square yang diperoleh dengan regresi Lasso terlihat lebih besar yakni :

R−squared:0.8643517

Perbandingan Performa Model

Didasarkan pada perbandingan nilai R-squarenya dapat diperoleh hasil sebagai berikut :

cat("Performa Model Regresi Klasik:", bestmodel$r.squared, "\n")

## Performa Model Regresi Klasik: 0.8584375

cat("Performa Model Regresi Ridge:", r_squared_ridge, "\n")

## Performa Model Regresi Ridge: 0.8577269

cat("Performa Model Regresi Lasso:", r_squared_lasso, "\n")

## Performa Model Regresi Lasso: 0.8641478

Diperoleh model lasso lebih unggul dibanding kedua model lainnya. Pada model lasso ini tidak ada peubah yang dikeluarkan dari model, hal ini dapat terjadi karena kurang dianggap adanya multikolinearitas pada data tiap peubahnya. Hal ini memvalidasi bahwa anggapan disinyalir adanya multikol hanya dugaan saja bukan berarti benar-benar terdeteksi ada multikol. Hasil akhir ini menggambarkan bahwa lasso menjadi model yang paling dapat menjabarkan kondisi indeks alibaca dalam data dengan tidak adanya peubah bebas yang dikeluarkan dalam data.

Sekian Terimakasih