Perbandingan Model Regresi Klasik, Regresi Ridge, dan Regresi Lasso
Pingkan Febbe Fiorela Kereh - G1401211087
2023-10-06
### Packages
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.3 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(glmnet)## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Loaded glmnet 4.1-8library(lmtest)## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(car)## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## somelibrary(corrplot)## corrplot 0.92 loadedlibrary(leaps)Input Data
Data yang digunakan adalah data Air Quality Index di New Delhi, yang bersumber dari kaggle.com. Data ini terdiri dari 72 observasi dan 12 peubah di dalamnya.
aqi <- read.csv("D:/COLLEGE/SMT 5/Pengantar Sains Data/aqi.csv")
glimpse(aqi)## Rows: 72
## Columns: 12
## $ X <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, …
## $ AQI <dbl> 30.2, 28.2, 26.6, 25.0, 26.0, 26.0, 26.0, 24.0, 23.0, …
## $ CO <dbl> 198.6027, 197.6013, 198.6027, 201.9405, 205.2784, 208.…
## $ datetime <chr> "2022-10-21:18", "2022-10-21:19", "2022-10-21:20", "20…
## $ no2 <dbl> 0.04685717, 0.04645554, 0.04685717, 0.04819594, 0.0488…
## $ o3 <dbl> 55.78995, 54.93164, 54.64554, 55.07469, 55.78995, 56.5…
## $ pm10 <dbl> 10.486722, 10.719325, 11.155578, 11.116206, 10.405250,…
## $ pm25 <dbl> 5.637410, 4.618169, 3.520902, 2.225919, 1.979471, 1.79…
## $ so2 <dbl> 0.3874302, 0.4097819, 0.4023313, 0.3762543, 0.3390014,…
## $ timestamp_local <chr> "2022-10-21T23:00:00", "2022-10-22T00:00:00", "2022-10…
## $ timestamp_utc <chr> "2022-10-21T18:00:00", "2022-10-21T19:00:00", "2022-10…
## $ ts <int> 1666375200, 1666378800, 1666382400, 1666386000, 166638…Cek apakah ada missing value
sum(is.na(aqi))## [1] 0Filtering Data
Dilakukan filtering data untuk memilih beberapa peubah yang relevan dengan analisis yang akan dilakukan.
y <- aqi$AQI
x1 <- aqi$CO
x2 <- aqi$no2
x3 <- aqi$o3
x4 <- aqi$pm10
x5 <- aqi$pm25
x6 <- aqi$so2
data <- cbind(y,x1,x2,x3,x4,x5,x6)
df <- as.data.frame(data)
head(df)## y x1 x2 x3 x4 x5 x6
## 1 30.2 198.6027 0.04685717 55.78995 10.486722 5.637410 0.3874302
## 2 28.2 197.6013 0.04645554 54.93164 10.719325 4.618169 0.4097819
## 3 26.6 198.6027 0.04685717 54.64554 11.155578 3.520902 0.4023313
## 4 25.0 201.9405 0.04819594 55.07469 11.116206 2.225919 0.3762543
## 5 26.0 205.2784 0.04886533 55.78995 10.405250 1.979471 0.3390014
## 6 26.0 208.6163 0.05020411 56.50520 9.402518 1.792677 0.3129244Eksplorasi Data
Matriks Korelasi
df.cor = cor(df)
corrplot.mixed(df.cor, lower = "ellipse", upper = "number",
tl.pos = "lt",diag = "l",tl.col = "black")
Berdasarkan matriks korelasi yang dihasilkan, peubah x3 yaitu O3
memiliki korelasi yang positif dan paling tinggi terhadap y (AQI). Dapat
dilihat pula peubah x1,x2,x3,dan x4 memiliki hubungan yang positif
terhadap y dan pada peubah x5 dan x6 memiliki hubungan yang cenderung
negatif terhadap y.
Regresi Klasik
Model Regresi Awal
model.rlb <- lm(y ~., data = df)
summary(model.rlb)##
## Call:
## lm(formula = y ~ ., data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73339 -0.25673 0.03901 0.24063 0.52176
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.46525 3.49796 -0.419 0.67668
## x1 0.01142 0.01674 0.682 0.49769
## x2 16.32158 5.28647 3.087 0.00297 **
## x3 0.45725 0.01196 38.241 < 2e-16 ***
## x4 -0.09697 0.02868 -3.381 0.00123 **
## x5 1.36334 0.08151 16.726 < 2e-16 ***
## x6 -9.45886 1.68392 -5.617 4.36e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3027 on 65 degrees of freedom
## Multiple R-squared: 0.991, Adjusted R-squared: 0.9902
## F-statistic: 1194 on 6 and 65 DF, p-value: < 2.2e-16Berdasarkan hasil perhitungan di atas, didapatkan model regresi sebagai berikut : \[yduga=-1.46525+0.01142x1+16.32158x2+0.45725x3-0.09697x4+1.36334x5-9.45886x6+e\]. Berdasarkan ringkasan model dapat diketahui bahwa hasil uji F memiliki p-value < \(\alpha\) (5%). Artinya, minimal terdapat satu variabel yang berpengaruh nyata terhadap model.
Diperoleh pula hasil uji t-parsial yaitu didapat beberapa parameter regresi yang signifikan dengan memiliki p-value < \(\alpha\) (5%) yaitu koefisien peubah x2,x3,x4,x5, dan x6 pada taraf nyata 5%.
Selanjutnya dapat dilihat juga nilai R^2=99.1. Artinya, sebesar 99.1% keragaman y dapat dijelaskan oleh peubah-peubah prediktornya.
Pengujian Asumsi
Prosedur Eksploratif
#Sisaan dan fitted value
sisaan<- residuals(model.rlb)
fitValue<- predict(model.rlb)
#Normal Q-Q Plot
qqnorm(sisaan);qqline(sisaan, col = "blue", lwd = 2)
#Plot Residual vs Fitted Value
plot(fitValue, sisaan, xlab = "Residual", ylab = "Fitted Values", main = "Residual vs Fitted Values", pch = 20);abline(h = 0, lwd = 1,col = "coral")
#Plot Residual vs Order
plot(sisaan, type = 'o', xlab = "Residual", ylab = "Order",
main = "Residual vs Order", pch = 20);abline(h = 0, lwd = 1, col = "purple")
#ACF dan PACF
par(mfrow = c(1,2))
acf(sisaan)
pacf(sisaan)
1)Berdasarkan Normal Q-Q Plot yang dihasilkan, terlihat bahwa sebaran titik-titik data cenderung menjauhi garis yang berarti mengindikasikan sisaan tidak menyebar normal. Namun hal ini akan diperiksa kembali melalui uji formal.
2)Berdasarkan Plot Residual vs Fitted Values yang dihasilkan, terlihat bahwa sisaan tidak menyebar secara acak dan memiliki lebar pita yang tidak sama yang berarti ragam sisaan tidak homogen. Namun hal ini akan diperiksa kembali melalui uji formal.
3)Berdasarkan Plot Residual vs Order yang dihasilkan, terlihat bahwa sisaan tidak menyebar secara acak melainkan memiliki pola naik turun (fluktuatif) yang berarti sisaan tidak saling bebas atau terindikasi adanya autokorelasi.
4)Berdasarkan plot ACF dan PACF yang dihasilkan, terlihat bahwa selain pada lag 0, tidak terdapat garis vertikal yang melewati batas signifikan. Oleh karena itu dapat disimpulkan bahwa tidak terdapat autokorelasi. Hal ini berbanding terbalik dengan hasil plot residual vs order dan hal ini akan diperiksa kembali melalui uji formal.
Uji Formal
#Nilai harapan sisaan = 0
t.test(model.rlb$residuals,
mu = 0,
conf.level = 0.95)##
## One Sample t-test
##
## data: model.rlb$residuals
## t = -4.8306e-17, df = 71, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.06806181 0.06806181
## sample estimates:
## mean of x
## -1.648891e-18\(H_0\) : Nilai harapan sisaan sama dengan nol
\(H_1\) : Nilai harapan sisaan tidak sama dengan nol
Berdasarkan hasil perhitungan di atas, diperoleh p-value = 1 lebih besar dari 0.05 sehingga tak tolak H0 yang artinya nilai harapan sisaan sama dengan nol.
#Uji Normalitas
ks.test(model.rlb$residuals, "pnorm", mean=mean(model.rlb$residuals), sd=sd(model.rlb$residuals))##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: model.rlb$residuals
## D = 0.084463, p-value = 0.652
## alternative hypothesis: two-sided\(H_0\) : Sisaan menyebar normal
\(H_1\) : Sisaan tidak menyebar normal
Berdasarkan hasil perhitungan di atas, diperoleh p-value = 0.652 lebih besar dari 0.05 sehingga tak tolak H0 yang artinya sisaan menyebar normal.
#Uji Kehomogenan Ragam
bptest(model.rlb)##
## studentized Breusch-Pagan test
##
## data: model.rlb
## BP = 11.684, df = 6, p-value = 0.0694\(H_0\) : Ragam sisaan homogen
\(H_1\) : Ragam sisaan tidak homogen
Berdasarkan hasil perhitungan di atas, diperoleh p-value = 0.0694 lebih besar dari 0.05 sehingga tak tolak H0 yang artinya ragam sisaan homogen.
#Uji Autokorelasi
dwtest(model.rlb)##
## Durbin-Watson test
##
## data: model.rlb
## DW = 1.9256, p-value = 0.1541
## alternative hypothesis: true autocorrelation is greater than 0\(H_0\) : Tidak terdapat autokorelasi
\(H_1\) : Terdapat autokorelasi
Berdasarkan hasil perhitungan di atas, diperoleh p-value = 0.1541 lebih besar dari 0.05 sehingga tak tolak H0 yang artinya tidak terdapat autokorelasi.
#Uji Multikolinieritas
vif(model.rlb)## x1 x2 x3 x4 x5 x6
## 4.429015 2.591475 4.582356 1.495781 2.218728 4.382814Berdasarkan hasil di atas, nilai vif < 10 pada masing-masing peubah prediktor sehingga tidak mengindikasikan adanya multikolinearitas.
Variable Selection
Best Subset Selection
library(olsrr)##
## Attaching package: 'olsrr'## The following object is masked from 'package:datasets':
##
## riversmodel.bs <- ols_step_best_subset(model.rlb)
model.bs## Best Subsets Regression
## --------------------------------
## Model Index Predictors
## --------------------------------
## 1 x3
## 2 x3 x5
## 3 x3 x5 x6
## 4 x3 x4 x5 x6
## 5 x2 x3 x4 x5 x6
## 6 x1 x2 x3 x4 x5 x6
## --------------------------------
##
## Subsets Regression Summary
## -----------------------------------------------------------------------------------------------------------------------------------
## Adj. Pred
## Model R-Square R-Square R-Square C(p) AIC SBIC SBC MSEP FPE HSP APC
## -----------------------------------------------------------------------------------------------------------------------------------
## 1 0.9479 0.9472 0.9462 308.7050 157.3962 -51.4749 164.2262 35.5058 0.5068 0.0071 0.0551
## 2 0.9838 0.9833 0.9782 51.3849 75.4439 -131.5020 84.5506 11.2267 0.1624 0.0023 0.0176
## 3 0.9884 0.9878 0.9858 20.2210 53.5390 -151.9911 64.9223 8.1751 0.1198 0.0017 0.0130
## 4 0.9894 0.9888 0.9872 14.5514 48.6642 -156.2645 62.3242 7.5432 0.1120 0.0016 0.0122
## 5 0.9909 0.9903 0.9892 5.4651 39.4001 -163.7491 55.3368 6.5500 0.0985 0.0014 0.0107
## 6 0.9910 0.9902 0.989 7.0000 40.8868 -161.9558 59.1002 6.6051 0.1005 0.0014 0.0109
## -----------------------------------------------------------------------------------------------------------------------------------
## AIC: Akaike Information Criteria
## SBIC: Sawa's Bayesian Information Criteria
## SBC: Schwarz Bayesian Criteria
## MSEP: Estimated error of prediction, assuming multivariate normality
## FPE: Final Prediction Error
## HSP: Hocking's Sp
## APC: Amemiya Prediction CriteriaBerdasarkan hasil variable selection menggunakan metode best subset, model ke-5 menghasilkan nilai adjusted R-square paling tinggi dan nilai AIC terkecil. Sehingga model ke-5 merupakan model terbaik yang terdiri dari 5 peubah prediktor yaitu x2, x3, x4, x5 dan x6.
Forward Stepwise Selection
model.fwd <- ols_step_forward_p(model.rlb)
model.fwd##
## Selection Summary
## --------------------------------------------------------------------------
## Variable Adj.
## Step Entered R-Square R-Square C(p) AIC RMSE
## --------------------------------------------------------------------------
## 1 x3 0.9479 0.9472 308.7050 157.3962 0.7022
## 2 x5 0.9838 0.9833 51.3849 75.4439 0.3948
## 3 x6 0.9884 0.9878 20.2210 53.5390 0.3369
## 4 x4 0.9894 0.9888 14.5514 48.6642 0.3236
## 5 x2 0.9909 0.9903 5.4651 39.4001 0.3015
## --------------------------------------------------------------------------Berdasarkan hasil metode forward di atas, peubah yang pertama kali dimasukkan ke dalam model yaitu X3 dilanjutkan dengan X5,x6,x4 dan X2. Penambahan X5 setelah X3 dans seterusnya menghasilkan perubahan kriteria yang besar (adj R-square yang menjadi semakin tinggi) dan model menjadi lebih baik.
Backward Stepwise Selection
model.bwd <- ols_step_backward_p(model.rlb)
model.bwd##
##
## Elimination Summary
## -----------------------------------------------------------------------
## Variable Adj.
## Step Removed R-Square R-Square C(p) AIC RMSE
## -----------------------------------------------------------------------
## 1 x1 0.9909 0.9903 5.4651 39.4001 0.3015
## -----------------------------------------------------------------------Berdasarkan hasil metode backward di atas, hanya terdapat satu peubah penjelas yang dikeluarkan dari model penuh, yaitu peubah X1. Model akhir yang diperoleh dari metode ini sama dengan model akhir (step 5) dari forward selection.
Model Hasil Variable Selection
Dari ketiga metode variable selection sebelumnya, menghasilkan model yang sama yaitu dengan peubah x2,x3,x4,x5,dan x6. Sehingga diperoleh model baru sebagai berikut :
new.model <- lm(y ~ x2+x3+x4+x5+x6, data = df)
summary(new.model)##
## Call:
## lm(formula = y ~ x2 + x3 + x4 + x5 + x6, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73766 -0.26339 0.03505 0.23931 0.55032
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.81836 1.00697 0.813 0.419315
## x2 17.14110 5.12717 3.343 0.001368 **
## x3 0.46048 0.01094 42.092 < 2e-16 ***
## x4 -0.10275 0.02729 -3.765 0.000357 ***
## x5 1.36303 0.08118 16.791 < 2e-16 ***
## x6 -9.90519 1.54522 -6.410 1.79e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3015 on 66 degrees of freedom
## Multiple R-squared: 0.9909, Adjusted R-squared: 0.9903
## F-statistic: 1445 on 5 and 66 DF, p-value: < 2.2e-16Berdasarkan hasil perhitungan di atas, didapatkan model regresi sebagai berikut : \[yduga=0.81836+17.14110x2+0.46048 x3-0.10275x4+1.36303x5-9.90519x6+e\]. Berdasarkan ringkasan model dapat diketahui bahwa hasil uji F memiliki p-value < \(\alpha\) (5%). Artinya, minimal terdapat satu variabel yang berpengaruh nyata terhadap model.
Diperoleh pula hasil uji t-parsial yaitu didapat seluruh parameter regresi kecuali intercept memiliki p-value < \(\alpha\) (5%) sehingga signifikan pada taraf nyata 5%.
Selanjutnya dapat dilihat juga nilai R^2=99.09%. Artinya, sebesar 99.09% keragaman y dapat dijelaskan oleh peubah-peubah prediktornya.
Regresi Ridge
x <- cbind(x1,x2,x3,x4,x5,x6)
y <- df$yCV
cv.r<-cv.glmnet(x,y,alpha=0);plot(cv.r)
Best Model Ridge
best.lr <-cv.r$lambda.min
bestridge <- glmnet(x,y,alpha=0,lambda=best.lr);coef(bestridge)## 7 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) -9.25162966
## x1 0.08367469
## x2 34.55907448
## x3 0.33780048
## x4 -0.05536372
## x5 1.07011621
## x6 -10.52153546Model Regresi Ridge :
\[y_duga=-9.25162966+0.08367469x134.55907448x2+0.33780048x3-0.05536372x4+1.07011621x5-10.52153546x6+e\] Interpretasi : - \(Intercept\) : nilai dugaan rataan y (AQI) ketika seluruh peubah x bernilai 0 sebesar -9.25162966. - \(x1 (CO)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar 0.08367469 pada setiap kenaikan satu satuan variabel x1 (CO). - \(x2 (no2)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar 34.55907448 pada setiap kenaikan satu satuan variabel x2 (no2). - \(x3 (o3)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar 0.33780048 pada setiap kenaikan satu satuan variabel x3 (o3). - \(x4 (pm10)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar -0.05536372 pada setiap kenaikan satu satuan variabel x4 (pm10). - \(x5 (pm25)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar 1.07011621 pada setiap kenaikan satu satuan variabel x5 (pm25). - \(x6 (so2)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar -10.52153546 pada setiap kenaikan satu satuan variabel x6 (so2).
Fungsi R-square
rsq<-function(bestmodel,bestlambda,x,y){
#y duga
y.duga <- predict(bestmodel, s = bestlambda, newx = x)
#JKG dan JKT
jkt <- sum((y - mean(y))^2)
jkg <- sum((y.duga- y)^2)
#find R-Squared
rsq <- 1 - jkg/jkt
return(rsq)
}R-square Ridge
rsq(bestridge,best.lr,x,y)## [1] 0.9744175Diperoleh nilai R-Squre model ridge sebesar 97.44%
Regresi Lasso
CV
cv.l<-cv.glmnet(x,y,alpha=1);plot(cv.l)
Best Model Lasso
best.ll<-cv.l$lambda.min
bestlasso<-glmnet(x,y,alpha=1,lambda=best.ll);coef(bestlasso)## 7 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) -1.79850521
## x1 0.01216225
## x2 15.29552875
## x3 0.45776935
## x4 -0.08576635
## x5 1.32314730
## x6 -8.94032751Model Regresi Lasso :
\[y_duga=-1.79850521+0.01216225x1+15.29552875x2+0.45776935x3-0.08576635x4+1.323147301x5-8.94032751x6+e\] Interpretasi : - \(Intercept\) : nilai dugaan rataan y (AQI) ketika seluruh peubah x bernilai 0 sebesar -1.79850521. - \(x1 (CO)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar 0.01216225 pada setiap kenaikan satu satuan variabel x1 (CO). - \(x2 (no2)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar 15.29552875 pada setiap kenaikan satu satuan variabel x2 (no2). - \(x3 (o3)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar 0.45776935 pada setiap kenaikan satu satuan variabel x3 (o3). - \(x4 (pm10)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar -0.08576635 pada setiap kenaikan satu satuan variabel x4 (pm10). - \(x5 (pm25)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar 1.323147301 pada setiap kenaikan satu satuan variabel x5 (pm25). - \(x6 (so2)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar -8.94032751 pada setiap kenaikan satu satuan variabel x6 (so2).
R-square Lasso
rsq(bestlasso,best.ll,x,y)## [1] 0.9909643Diperoleh nilai R-Squre model ridge sebesar 99.09%
Membandingkan R-Square
#### R-Square model regresi awal
summary(model.rlb)$r.squared## [1] 0.9910094#### R-Square model hasil variable selection
summary(new.model)$r.squared## [1] 0.990945#### R-Square model regresi ridge
rsq(bestridge,best.lr,x,y)## [1] 0.9744175#### R-Square model regresi lasso
rsq(bestlasso,best.ll,x,y)## [1] 0.9909643Berdasarkan nilai R-Square dari masing-masing model, diketahui model regresi awal memiliki nilai R-Square paling tinggi, sehingga model awal (model regresi klasik awal) merupakan model regresi terbaik.
Interpretasi Model Terbaik
\[yduga=-1.46525+0.01142x1+16.32158x2+0.45725x3-0.09697x4+1.36334x5-9.45886x6+e\].
Interpretasi : - \(Intercept\) : nilai dugaan rataan y (AQI) ketika seluruh peubah x bernilai 0 sebesar -1.46525. - \(x1 (CO)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar 0.01142 pada setiap kenaikan satu satuan variabel x1 (CO). - \(x2 (no2)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar 16.32158 pada setiap kenaikan satu satuan variabel x2 (no2). - \(x3 (o3)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar 0.45725 pada setiap kenaikan satu satuan variabel x3 (o3). - \(x4 (pm10)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar -0.09697 pada setiap kenaikan satu satuan variabel x4 (pm10). - \(x5 (pm25)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar 1.36334 pada setiap kenaikan satu satuan variabel x5 (pm25). - \(x6 (so2)\) : nilai dugaan perubahan dugaan rataan y (AQI) sebesar -9.45886 pada setiap kenaikan satu satuan variabel x6 (so2).