Menentukan Model Regresi Terbaik

Dinda Khamila Nurfatimah

2023-09-28

Data dan Eksplorasi Data

Input Data

library(rio)
data <- import("https://raw.githubusercontent.com/DindaKhamila/PSD/main/Data%20Tugas/StudentPerformance.csv")
y <- data$'Performance Index'
x1 <- data$'Hours Studied'
x2 <- data$'Previous Scores'
x3 <- data$'Sleep Hours'
x4 <- data$'Sample Question Papers Practiced'

head(data)

##   Hours Studied Previous Scores Sleep Hours Sample Question Papers Practiced
## 1             7              99           9                                1
## 2             4              82           4                                2
## 3             8              51           7                                2
## 4             5              52           5                                2
## 5             7              75           8                                5
## 6             3              78           9                                6
##   Performance Index
## 1                91
## 2                65
## 3                45
## 4                36
## 5                66
## 6                61

Ekplorasi Data

Histogram

# Sebaran peubah y (Performance Index)
hist(y, col = "skyblue")

Boxplot

# Sebaran peubah y (Performance Index) 
boxplot(y, col = "skyblue")

Matriks Korelasi

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 errors

df <- data.frame(y,x1,x2,x3,x4)
df %>%
  as_tibble() %>%
  cor() %>%
  ggcorrplot::ggcorrplot(type = "upper", lab = TRUE, lab_size = 2, colors = "skyblue") +
  theme_minimal() +
  labs(title = "Hubungan Antar Peubah",
       subtitle = "Peubah respon : Performance Index",
       x = NULL, y = NULL)

Regresi Klasik

Pemodelan

modelreg <- lm(y ~ x1 + x2 + x3 + x4)
summary(modelreg)

## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.3299 -1.3831 -0.0062  1.3701  8.4864 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -33.763726   0.126841 -266.19   <2e-16 ***
## x1            2.853429   0.007962  358.40   <2e-16 ***
## x2            1.018584   0.001189  857.02   <2e-16 ***
## x3            0.476333   0.012153   39.19   <2e-16 ***
## x4            0.195198   0.007189   27.15   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.061 on 9995 degrees of freedom
## Multiple R-squared:  0.9885, Adjusted R-squared:  0.9885 
## F-statistic: 2.147e+05 on 4 and 9995 DF,  p-value: < 2.2e-16

Didapatkan model regresi klasik dengan R-squared sebesar 98.85 %, artinya sebesar 98.85% keragaman nilai performance index dapat dijelaskan oleh semua peubah penjelas. Hasil tersebut menunjukkan hasil yang sangat baik. Semua peubah penjelas memiliki p-value < 0.05 artinya cukup bukti untuk menyatakan bahwa peubah penjelas berpengaruh signifikan terhadap peubah respon.

Model Regresi Klasik yang dihasilkan: \[Y=-33.763726+2.853429X1+1.018584X2+0.476333X3+0.195198X4\]

Interpretasi:

Intersep : nilai rataan peubah respon (Performance Index) ketika seluruh peubah penjelas bernilai 0 sebesar -33.763726
X1 (Hours Studied) : 2.853429, artinya jika Hours Studied bertambah satu satuan, maka Performance Index naik dengan rata-rata sebesar 2.853429
X2 (Previous Scores) : 1.018584, artinya jika Previous Scores bertambah satu satuan, maka Performance Index naik dengan rata-rata sebesar 1.018584
X3 (Sleep Hours) : 0.476333, artinya jika Sleep Hours bertambah satu satuan, maka Performance Index naik dengan rata-rata sebesar 0.476333
X4 (Sample Question Papers Practiced) : 0.195198, artinya jika Sample Question Papers Practiced bertambah satu satuan, maka Performance Index naik dengan rata-rata sebesar 0.195198

Pengujian Asumsi

#sisaan dan fitted value
sisaan<- residuals(modelreg)
fitValue<- predict(modelreg)

#Diagnostik dengan eksploratif
par(mfrow = c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "skyblue", lwd = 2)
plot(fitValue, sisaan, col = "skyblue", pch = 20, xlab = "Sisaan", ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)
hist(sisaan, col = "skyblue")
plot(seq(1,10000,1), sisaan, col = "skyblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,10000,1), sisaan, col = "pink")
abline(a = 0, b = 0, lwd = 2)

1. Uji Multikolinearitas

car::vif(modelreg)

##       x1       x2       x3       x4 
## 1.000464 1.000254 1.000052 1.000386

Tidak ada nilai VIF yang lebih dari 10, sehingga dapat disimpulkan bahwa tidak terdapat multikolinearitas antar masing-masing peubah penjelas.

2. Sisaan menyebar normal

plot(modelreg, 2);

nortest::ad.test(modelreg$residuals)

## 
##  Anderson-Darling normality test
## 
## data:  modelreg$residuals
## A = 0.26989, p-value = 0.6781

H0: sisaan mengikuti sebaran normal

H1: sisaan tidak mengikuti sebaran normal

Didapatkan nilai p-value (0.6781) > 0.05. Artinya, cukup bukti untuk menyatakan sisaan berdistribusi normal.

Dari QQ-Plot tersebut terlihat bahwa titik-titiknya cenderung mengikuti garis kenormalan. Sehingga dapat disimpulkan bahwa sisaan menyebar normal.

3. Nilai harapan sisaan sama dengan nol

# Uji t
t.test(resid(modelreg), mu = 0,)

## 
##  One Sample t-test
## 
## data:  resid(modelreg)
## t = -1.1215e-15, df = 9999, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.04038966  0.04038966
## sample estimates:
##     mean of x 
## -2.310923e-17

\(p-value=1 > 0.1\) (tak tolak \(H_0\)), artinya nilai harapan sisaan sama dengan nol

4. Ragam sisaan homogen

# Uji Breusch-Pagan
lmtest::bptest(modelreg)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelreg
## BP = 1.5045, df = 4, p-value = 0.8258

H0: sisaan homogen

H1: sisaan tidak homogen

\(p-value=0.8258 > 0.1\) (tolak \(H_0\)), artinya ragam sisaan homogen

5. Deteksi autokorelasi (tidak terjadi autokorelasi)

lmtest::dwtest(modelreg)

## 
##  Durbin-Watson test
## 
## data:  modelreg
## DW = 2.0039, p-value = 0.5764
## alternative hypothesis: true autocorrelation is greater than 0

H0: tidak ada autokorelasi

H1: ada autokorelasi

Berdasarkan hasil DW Test, didapatkan nilai DW=2.0039 dan p-value = 0.5764

\(p-value=0.5764 > 0.1\) (tak tolak \(H_0\)), artinya tidak terjadi autokorelasi pada sisaan pada taraf 5%.

Seleksi Peubah

Metode Stepwise

stepwise <- step(lm(y ~ 1), direction="both", scope=formula(modelreg), trace=1)

## Start:  AIC=59112.28
## y ~ 1
## 
##        Df Sum of Sq     RSS   AIC
## + x2    1   3091353  599501 40939
## + x1    1    515518 3175337 57610
## + x3    1      8541 3682313 59091
## + x4    1      6910 3683945 59096
## <none>              3690855 59112
## 
## Step:  AIC=40939.13
## y ~ x2
## 
##        Df Sum of Sq     RSS   AIC
## + x1    1    547358   52143 16520
## + x3    1      6719  592782 40828
## + x4    1      4797  594704 40861
## <none>               599501 40939
## - x2    1   3091353 3690855 59112
## 
## Step:  AIC=16520.02
## y ~ x2 + x1
## 
##        Df Sum of Sq     RSS   AIC
## + x3    1      6560   45583 15177
## + x4    1      3167   48976 15895
## <none>                52143 16520
## - x1    1    547358  599501 40939
## - x2    1   3123194 3175337 57610
## 
## Step:  AIC=15177.47
## y ~ x2 + x1 + x3
## 
##        Df Sum of Sq     RSS   AIC
## + x4    1      3131   42452 14468
## <none>                45583 15177
## - x3    1      6560   52143 16520
## - x1    1    547199  592782 40828
## - x2    1   3121377 3166960 57585
## 
## Step:  AIC=14467.83
## y ~ x2 + x1 + x3 + x4
## 
##        Df Sum of Sq     RSS   AIC
## <none>                42452 14468
## - x4    1      3131   45583 15177
## - x3    1      6524   48976 15895
## - x1    1    545578  588030 40750
## - x2    1   3119575 3162027 57572

summary(stepwise)

## 
## Call:
## lm(formula = y ~ x2 + x1 + x3 + x4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.3299 -1.3831 -0.0062  1.3701  8.4864 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -33.763726   0.126841 -266.19   <2e-16 ***
## x2            1.018584   0.001189  857.02   <2e-16 ***
## x1            2.853429   0.007962  358.40   <2e-16 ***
## x3            0.476333   0.012153   39.19   <2e-16 ***
## x4            0.195198   0.007189   27.15   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.061 on 9995 degrees of freedom
## Multiple R-squared:  0.9885, Adjusted R-squared:  0.9885 
## F-statistic: 2.147e+05 on 4 and 9995 DF,  p-value: < 2.2e-16

Seleksi peubah menggunakan metode backward, forward, dan stepwise hasilnya sama, yaitu model terbaik dengan peubah x1, x2, x3, dan x4 dengan \[R-square=98.85\%\]

Regresi Ridge

Packages

lapply(c("glmnet","lmridge"),library,character.only=T)[[1]]

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## Loaded glmnet 4.1-8

##  [1] "glmnet"    "Matrix"    "lubridate" "forcats"   "stringr"   "dplyr"    
##  [7] "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"   "tidyverse"
## [13] "rio"       "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [19] "methods"   "base"

Peubah

x <- cbind(x1,x2,x3,x4)
y <- data$'Performance Index'

CV

cv.r <- cv.glmnet(x,y,alpha=0);plot(cv.r)

Pemodelan

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.3120776
## x1            2.6078448
## x2            0.9327929
## x3            0.4416559
## x4            0.1862048

Model Regresi Ridge yang diihasilkan :

\[Y=-26.3120776+2.6078448X1+0.9327929X2+0.4416559X3+0.1862048X4\]

Interpretasi:

Intersep : nilai rataan peubah respon (Performance Index) ketika seluruh peubah penjelas bernilai 0 sebesar -26.3120776
X1 (Hours Studied) : 2.6078448, artinya jika Hours Studied bertambah satu satuan, maka Performance Index naik dengan rata-rata sebesar 2.6078448
X2 (Previous Scores) : 0.9327929, artinya jika Previous Scores bertambah satu satuan, maka Performance Index naik dengan rata-rata sebesar 0.9327929
X3 (Sleep Hours) : 0.4416559, artinya jika Sleep Hours bertambah satu satuan, maka Performance Index naik dengan rata-rata sebesar 0.4416559
X4 (Sample Question Papers Practiced) : 0.1862048, artinya jika Sample Question Papers Practiced bertambah satu satuan, maka Performance Index naik dengan rata-rata sebesar 0.1862048

Model ini memiliki koefisien regresi yang sedikit berbeda dengan model regresi klasik sebelumnya.

R-Square Ridge

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)
}

rsq(bestridge,best.lr,x,y)

## [1] 0.9814513

Diperoleh nilai R-Square sebesar 0.9814513. Hasil ini berbeda dengan model sebelumnya.

Regresi Lasso

Peubah

x = matrix(c(x1,x2,x3,x4), ncol=4)
y = matrix(y)

CV

cv.l <- cv.glmnet(x,y,alpha=1)
plot(cv.l)

Pemodelan

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.0172257
## V1            2.8280113
## V2            1.0147712
## V3            0.4377285
## V4            0.1727848

Model Regresi Ridge yang diihasilkan :

\[Y=-33.0172257+2.8280113V1+1.0147712V2+0.4377285V3+0.1727848V4\]

Interpretasi:

Intersep : nilai rataan peubah respon (Performance Index) ketika seluruh peubah penjelas bernilai 0 sebesar -33.0172257
V1 (Hours Studied) : 2.8280113, artinya jika Hours Studied bertambah satu satuan, maka Performance Index naik dengan rata-rata sebesar 2.8280113
V2 (Previous Scores) : 1.0147712, artinya jika Previous Scores bertambah satu satuan, maka Performance Index naik dengan rata-rata sebesar 1.0147712
V3 (Sleep Hours) : 0.4377285, artinya jika Sleep Hours bertambah satu satuan, maka Performance Index naik dengan rata-rata sebesar 0.4377285
V4 (Sample Question Papers Practiced) : 0.1727848, artinya jika Sample Question Papers Practiced bertambah satu satuan, maka Performance Index naik dengan rata-rata sebesar 0.1727848

Model ini memiliki koefisien regresi yang sedikit berbeda dengan model regresi klasik sebelumnya.

R-Square Lasso

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) 
}

rsq(bestlasso,best.ll ,x,y)

## [1] 0.9884512

Diperoleh nilai R-Square sebesar 0.9884512. Hasil ini berbeda dengan model sebelumnya.

Perbandingan Model

Regresi Klasik

summary(modelreg)$sigma

## [1] 2.060898

Regresi Ridge

# Prediksi model ridge pada data pelatihan
train_predictionsridge <- predict(bestridge,newx = x)

# Hitung residu (selisih antara prediksi dan nilai sebenarnya)
residualsridge <- y - train_predictionsridge

# Hitung varian residu
dfridge <- length(y) - length(bestridge$beta)
residual_varianceridge <- sum(residualsridge^2) / dfridge

# Hitung RSE
rseridge <- sqrt(residual_varianceridge)

# Tampilkan hasil RSE
print(paste("Residual Standard Error (RSE):",rseridge))

## [1] "Residual Standard Error (RSE): 2.61702186011339"

Regresi Lasso

# Prediksi model Lasso pada data pelatihan
train_predictionsLasso <- predict(bestlasso,newx = x)

# Hitung residu (selisih antara prediksi dan nilai sebenarnya)
residualsLasso <- y - train_predictionsLasso

# Hitung varian residu
dfLasso <- length(y) - length(bestlasso$beta)
residual_varianceLasso <- sum(residualsLasso^2) / dfLasso

# Hitung RSE
rseLasso <- sqrt(residual_varianceLasso)

# Tampilkan hasil RSE
print(paste("Residual Standard Error (RSE):",rseLasso))

## [1] "Residual Standard Error (RSE): 2.06499451318242"

Model Regresi Terbaik

Model	R-Square	Residual Standard Error
Regresi Klasik	\(98.85\%\)	\(2.060898\)
Regresi Ridge	\(98.14\%\)	\(2.61702186011339\)
Regresi Lasso	\(98.84\%\)	\(2.06499451318242\)

Model Regresi Klasik memiliki R-Square paling besar dan RSE paling kecil sehingga model tersebut merupakan model terbaik. Berikut model terbaik yang diperoleh:

\[Y=-33.763726+2.853429X1+1.018584X2+0.476333X3+0.195198X4\]