Penentuan Model Terbaik Menggunakan Variable Selection, Ridge Regression dan Lasso Regression
Nazuwa Aulia
2023-09-29
Library yang digunakan
library(readxl)
library(stats)
library(base)
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, unionlibrary(corrplot)## corrplot 0.92 loadedlibrary(ggplot2)
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(randtests)
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(glmnet)## Loading required package: Matrix## Loaded glmnet 4.1-8library(lmridge)##
## Attaching package: 'lmridge'## The following object is masked from 'package:car':
##
## vifImport data yang digunakan
Data yang digunakan bersumber dari kaggle dengan judul Student Performance. Terdapat total 10000 data namun pada analisis ini hanya 1000 data yang digunakan. Adapun peubah bebas yang digunakan sebanyak 4.
Peubah Bebas (X)
Hours Studied (X1), jumlah total jam yang dihabiskan siswa untuk belajar. Previous Score (X2), nilai yang diperoleh siswa dalam tes sebelumnya. Sleep Hours (X3), jumlah rata-rata jam tidur yang dimiliki siswa per hari. Sample Question Papers Practiced (X4), jumlah soal latihan yang dikerjakan oleh siswa.
Peubah Respon (Y)
Performance Index, ukuran kinerja keseluruhan setiap siswa. Indeks kinerja menggambarkan kinerja akademik siswa dan telah dibulatkan ke bilangan bulat terdekat. Indeks ini berkisar dari 10 hingga 100, dengan nilai yang lebih tinggi menunjukkan kinerja yang lebih baik.
datapsd <- read_excel("C:/Users/Nazuwa Aulia/OneDrive/Documents/College/Semester 5/Pengantar Sains Data/student performance.xlsx", sheet = "Sheet2")
datapsd## # A tibble: 1,000 × 5
## y x1 x2 x3 x4
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 91 7 99 9 1
## 2 65 4 82 4 2
## 3 45 8 51 7 2
## 4 36 5 52 5 2
## 5 66 7 75 8 5
## 6 61 3 78 9 6
## 7 63 7 73 5 6
## 8 42 8 45 4 6
## 9 61 5 77 8 2
## 10 69 4 89 4 0
## # ℹ 990 more rowsstr(datapsd)## tibble [1,000 × 5] (S3: tbl_df/tbl/data.frame)
## $ y : num [1:1000] 91 65 45 36 66 61 63 42 61 69 ...
## $ x1: num [1:1000] 7 4 8 5 7 3 7 8 5 4 ...
## $ x2: num [1:1000] 99 82 51 52 75 78 73 45 77 89 ...
## $ x3: num [1:1000] 9 4 7 5 8 9 5 4 8 4 ...
## $ x4: num [1:1000] 1 2 2 2 5 6 6 6 2 0 ...Eksplorasi Data
summary(datapsd)## y x1 x2 x3
## Min. : 10.00 Min. :1.000 Min. :40.00 Min. :4.000
## 1st Qu.: 40.75 1st Qu.:3.000 1st Qu.:55.00 1st Qu.:5.000
## Median : 56.00 Median :5.000 Median :70.00 Median :7.000
## Mean : 55.41 Mean :4.898 Mean :69.94 Mean :6.516
## 3rd Qu.: 71.00 3rd Qu.:7.000 3rd Qu.:85.25 3rd Qu.:8.000
## Max. :100.00 Max. :9.000 Max. :99.00 Max. :9.000
## x4
## Min. :0.000
## 1st Qu.:2.000
## Median :5.000
## Mean :4.484
## 3rd Qu.:7.000
## Max. :9.000Korelasi Peubah
Pair plot
ggpairs(datapsd)
Plot pair diatas menunjukkan menunjukkan hubungan antar peubah. Pada tabel nilai korelasi terdapat tanda bintang (*) yang menandakan adanya hubungan linier yang kuat atau signifikan dengan peubah respons, yaitu X1 dan X2. Meskipun nilai korelasi cenderung lemah, tidak ada peubah penjelas yang dieliminasi pada model awal karena akan dilakukan variable selection untuk melihat model terbaik.
Matriks korelasi untuk setiap peubah numerik
corr_matx <- datapsd %>% select_if(is.numeric) %>% cor() %>% round(3)
corr_matx## y x1 x2 x3 x4
## y 1.000 0.395 0.919 0.061 0.052
## x1 0.395 1.000 0.019 -0.009 0.001
## x2 0.919 0.019 1.000 0.027 0.026
## x3 0.061 -0.009 0.027 1.000 -0.030
## x4 0.052 0.001 0.026 -0.030 1.000corrplot(corr_matx,
method = "color",
type = "lower",
tl.cex = 0.5,
tl.col = "black",
addCoef.col = "#2F2F2F",
addCoefasPercent = FALSE,
number.cex = 0.5,
diag=F)
Model Regresi Linear Berganda
model.reg <- lm(y ~., data = datapsd)
summary(model.reg)##
## Call:
## lm(formula = y ~ ., data = datapsd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.4160 -1.3150 -0.0403 1.2852 7.4600
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -33.816104 0.391403 -86.397 <2e-16 ***
## x1 2.829452 0.024672 114.684 <2e-16 ***
## x2 1.021915 0.003707 275.635 <2e-16 ***
## x3 0.461903 0.037863 12.199 <2e-16 ***
## x4 0.197141 0.022865 8.622 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.03 on 995 degrees of freedom
## Multiple R-squared: 0.9892, Adjusted R-squared: 0.9891
## F-statistic: 2.274e+04 on 4 and 995 DF, p-value: < 2.2e-16Model regresi linear berganda yang diperoleh adalah :
\[ \hat{y}=-33.816104+2.829452X_{1}+1.021915X_{2}+0.461903X_{i}+0.197141 X_{4} \] Selain itu dapat dilihat nilai Multiple R-squared yang diperoleh sebesar 0.9892
Pendeteksian Multikolinearitas
car::vif(model.reg)## x1 x2 x3 x4
## 1.000450 1.001840 1.001763 1.001602Multikolinieritas terdeteksi apabila nilai VIF lebih dari 10. Pada output diatas dapat dilihat bahwa hasil semua nilai VIF untuk setiap peubah kurang dari 10. Artinya, tidak ada multikolinieritas pada setiap peubah penjelas.
Uji Asumsi
Harapan Sisaan = 0
plot(model.reg,1, pch=20)
t.test(model.reg$residuals,mu = 0,conf.level = 0.95)##
## One Sample t-test
##
## data: model.reg$residuals
## t = 1.9322e-15, df = 999, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.1257231 0.1257231
## sample estimates:
## mean of x
## 1.237894e-16Nilai p-value = 1 > 0.05 sehingga tak tolak H0 artinya nilai harapan sisaan = 0. Secara eksploratif juga dapat dilihat pada plot dimana sebagian besar sisaan berada di sekitar 0. Sehingga dapat disimpulkan bahwa terdapat cukup bukti untuk menyatakan nilai harapan sisaan = 0 dan asumsi terpenuhi
Kehomogenan Sisaan
plot(model.reg,1, pch=20)
bptest(model.reg)##
## studentized Breusch-Pagan test
##
## data: model.reg
## BP = 3.3844, df = 4, p-value = 0.4957Nilai p-value = 0.4957 > 0.05 sehingga tak tolak H0 artinya ragam sisaan homogen. Secara eksploratif juga dapat dilihat pada plot dimana lebar pita sama untuk setiap nilai dugaan dan tidak membentuk seperti corong. Sehingga dapat disimpulkan bahwa terdapat cukup bukti untuk menyatakan ragam sisaan homogen dan asumsi terpenuhi
Kesalingbebasan Sisaan
plot(x = 1:dim(datapsd)[1],
y = model.reg$residuals,
type = 'b',
ylab = "Residuals",
xlab = "Observation", pch=20)
dwtest(model.reg)##
## Durbin-Watson test
##
## data: model.reg
## DW = 1.9771, p-value = 0.3592
## alternative hypothesis: true autocorrelation is greater than 0runs.test(model.reg$residuals)##
## Runs Test
##
## data: model.reg$residuals
## statistic = -0.063277, runs = 500, n1 = 500, n2 = 500, n = 1000,
## p-value = 0.9495
## alternative hypothesis: nonrandomnessNilai p-value = 0.3592 > 0.05 sehingga tak tolak H0 artinya sisaan saling bebas dan tidak terdapat autokorelasi. Secara eksploratif juga dapat dilihat pada plot dimana sisaan menyebar tak berpola, sehingga asumsi sisaan saling bebas terpenuhi. Sehingga dapat disimpulkan bahwa terdapat cukup bukti untuk menyatakan sisaan saling bebas dan asumsi terpenuhi
Normalitas Sisaan
plot(model.reg,2,pch=20)
ks.test(model.reg$residuals, "pnorm", mean=mean(model.reg$residuals), sd=sd(model.reg$residuals))## Warning in ks.test.default(model.reg$residuals, "pnorm", mean =
## mean(model.reg$residuals), : ties should not be present for the
## Kolmogorov-Smirnov test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: model.reg$residuals
## D = 0.021432, p-value = 0.7479
## alternative hypothesis: two-sidedNilai p-value = 0.7479 > 0.05 sehingga tak tolak H0 artinya sisaan menyebar normal. Secara eksploratif juga dapat dilihat pada plot dimana banyak amatan telah mendekati garis qq-plot normal, sehingga asumsi sisaan menyebar normal terpenuhi. Sehingga dapat disimpulkan bahwa terdapat cukup bukti untuk menyatakan sisaan menyebar normal dan asumsi terpenuhi
Variable Selection
Forward Selection
forward_model <- stepAIC(model.reg, direction = "forward")## Start: AIC=1421.13
## y ~ x1 + x2 + x3 + x4Backward Selection
backward_model <- step(model.reg, direction = "backward")## Start: AIC=1421.13
## y ~ x1 + x2 + x3 + x4
##
## Df Sum of Sq RSS AIC
## <none> 4101 1421.1
## - x4 1 306 4407 1491.2
## - x3 1 613 4714 1558.5
## - x1 1 54204 58304 4073.7
## - x2 1 313107 317207 5767.6Stepwise Selection
stepwise_model <- stepAIC(model.reg, direction = "both")## Start: AIC=1421.13
## y ~ x1 + x2 + x3 + x4
##
## Df Sum of Sq RSS AIC
## <none> 4101 1421.1
## - x4 1 306 4407 1491.2
## - x3 1 613 4714 1558.5
## - x1 1 54204 58304 4073.7
## - x2 1 313107 317207 5767.6Ridge Regression (glmnet)
x<-data.matrix(datapsd[, c('x1', 'x2', 'x3', 'x4')])
y<-datapsd$yCV
cv.r<-cv.glmnet(x,y,alpha=0);plot(cv.r)
Best Model
best.lr<-cv.r$lambda.min
bestridge<-glmnet(x,y,alpha=0,lambda=best.lr);coef(bestridge)## 5 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) -26.5793011
## x1 2.6010253
## x2 0.9365830
## x3 0.4416825
## x4 0.1930866R-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.9823767Lasso Regression
CV
cv.l<-cv.glmnet(x,y,alpha=1);plot(cv.l)
Best Model
best.ll<-cv.l$lambda.min
bestlasso<-glmnet(x,y,alpha=1,lambda=best.ll);coef(bestlasso)## 5 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) -33.0685544
## x1 2.8038169
## x2 1.0183117
## x3 0.4216845
## x4 0.1730727R-Squared Lasso
rsq(bestlasso,best.ll,x,y)## [1] 0.9891326Ridge Regression (lmridge)
lmr<-lmridge(y~.,datapsd,scaling="centered");plot(lmr);vif(lmr)
## x1 x2 x3 x4
## k=0 0.00015 0 0.00035 0.00013summary(lmr)##
## Call:
## lmridge.default(formula = y ~ ., data = datapsd, scaling = "centered")
##
##
## Coefficients: for Ridge parameter K= 0
## Estimate Estimate (Sc) StdErr (Sc) t-value (Sc) Pr(>|t|)
## Intercept -33.8161 -33.8161 0.7297 -46.3444 < 2.2e-16 ***
## x1 2.8294 2.8295 0.0247 114.7416 < 2.2e-16 ***
## x2 1.0219 1.0219 0.0037 275.7735 < 2.2e-16 ***
## x3 0.4619 0.4619 0.0378 12.2054 < 2.2e-16 ***
## x4 0.1971 0.1971 0.0229 8.6265 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Ridge Summary
## R2 adj-R2 DF ridge F AIC BIC
## 0.98920 0.98910 3.99999 22763.38879 1419.13130 8346.51755
## Ridge minimum MSE= 0.002576261 at K= 0
## P-value for F-test ( 3.99999 , 996 ) = 0
## -------------------------------------------------------------------Kesimpulan
Berdasarkan hasil analisis yang dilakukan melalui variable selection, ridge regression dan lasso regression model terbaik yang dihasilkan merupakan model awal yang memiliki 4 peubah bebas dengan nilai R-Squared nya sebesar 98.92%