psd

Data

Menggunakan dataset yang terdapat di R, yaitu dataset mtcars. Dataset mtcars merupakan kumpulan data bawaan di R yang berisi pengukuran pada 7 atribut berbeda untuk200 baris data yang digunakan:

Peubah

y = Fartility (Kesuburan)

x1 = Examination (telah menjalani pemeriksaan militer)

x2 = Agriculture (lahan pertanian)

x3 = Catholic(agama Catholic)

x4 = Infant.Mortality (kematian bayi per 1000 kelahiran.))

data("swiss")
y <- swiss$Fertility
x1 <- swiss$Examination
x2 <- swiss$Agriculture
x3 <- swiss$Catholic
x4 <- swiss$Infant.Mortality

Ekplorasi

# Sebaran peubah y (Fertility)
hist(y, col = "navy")

boxplot(y, col = "navy")

Matriks Korelasi

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.2.3

## Warning: package 'ggplot2' was built under R version 4.2.3

## Warning: package 'tibble' was built under R version 4.2.3

## Warning: package 'tidyr' was built under R version 4.2.3

## Warning: package 'readr' was built under R version 4.2.3

## Warning: package 'dplyr' was built under R version 4.2.3

## Warning: package 'forcats' was built under R version 4.2.3

## Warning: package 'lubridate' was built under R version 4.2.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; 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 = "red") +
  theme_minimal() +
  labs(title = "Hubungan Antar Peubah",
       subtitle = "Peubah respon : kesuburan",
       x = NULL, y = NULL)

Dari eksplorasi menggunakan plot di atas, terlihat bahwa peubah x1, x2, x3, dan x4 memiliki nilai korelasi yang cukup tinggi terhadap y. Peubah-peubah tersebut yang akan digunakan sebagai peubah penjelas dalam tahapan analisis berikutnya.Peubah responnya adalah y

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 
## -23.9194  -3.5530  -0.6489   6.5956  14.1767 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 59.60267   13.04246   4.570 4.25e-05 ***
## x1          -0.96805    0.25284  -3.829 0.000423 ***
## x2          -0.04759    0.08032  -0.593 0.556688    
## x3           0.02611    0.03843   0.679 0.500551    
## x4           1.39597    0.46259   3.018 0.004315 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.82 on 42 degrees of freedom
## Multiple R-squared:  0.5448, Adjusted R-squared:  0.5014 
## F-statistic: 12.57 on 4 and 42 DF,  p-value: 8.272e-07

Didapatkan model regresi klasik dengan R-squared sebesar 54.48%. Model ini memiliki p-value < 0.05 artinya cukup bukti untuk menyatakan bahwa minimal ada satu peubah penjelas yang berpengaruh signifikan terhadap peubah respon. Berdasarkan hasil uji t, peubah penjelas yang berpengaruh yaitu peubah x1 (Examination) danx4 (catholic).

\[Y=59.60267-0.96805X1-0.04759X2+0.02611x3+1.39597X4\]

Interpretasi:

• Intersep: nilai rataan peubah respon (Fertility) ketika seluruh peubah penjelas bernilai 0 sebesar 40.82854

• X1 (Examination): -0.96805, artinya jika Examination bertambah satu satuan, maka Fertility turun dengan rata-rata sebesar 0.96805

• X2 (Agriculture): -0.04759, artinya jika Agriculture bertambah satu satuan, maka Fertility turun dengan rata-rata sebesar 0.04759

• X3 (catholic): 0.02611, artinya jika catholic bertambah satu satuan, maka Fertility naik dengan rata-rata sebesar 0.02611

• X4 (Infant.Mortality): 1.39597, artinya jika Infant.Mortality bertambah satu satuan, maka Fertility naik dengan rata-rata sebesa 1.39597

Multikolinearitas

car::vif(modelreg)

##       x1       x2       x3       x4 
## 2.405776 1.967676 1.518401 1.073400

Metode VIF digunakan untuk mendeteksi multikolinearitas yang terjadi ketika VIF > 10. Diperoleh bahwa tidak terdapat multikolinearitas pada semua peubah.

plot(modelreg, 2);

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

\(H_0 : E[\varepsilon]=0\) \(H_1 : E[\varepsilon]\ne0\)

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

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

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

4. Ragam sisaan homogen

\(H_0: Var[\varepsilon]=\sigma^2I\) \(H_1: Var[\varepsilon]\ne \sigma^2I\)

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

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

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

5. Deteksi autokorelasi (tidak terjadi autokorelasi)

\(H_0 : Cov[\varepsilon_i,\varepsilon_j]=0\) (tidak terjadi autokorelasi pada sisaan) \(H_0 : Cov[\varepsilon_i,\varepsilon_j]\neq0\) (terjadi autokorelasi pada sisaan)

# Uji Durbin Watson
library(lmtest)

## Warning: package 'lmtest' was built under R version 4.2.3

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

dwtest(modelreg)

## 
##  Durbin-Watson test
## 
## data:  modelreg
## DW = 1.0238, p-value = 3.723e-05
## alternative hypothesis: true autocorrelation is greater than 0

#ACF dan PACF identifikasi autokorelasi

sisaan = modelreg$residuals
par(mfrow = c(1,2))
acf(sisaan)
pacf(sisaan)

\(p-value= 3.723e-05 < 0.1\) (tolak \(H_0\)), artinya terjadi autkorelasi pada sisaan pada taraf 5%.

Seleksi Peubah

Metode Backward

backward <- step(modelreg, direction="backward", scope=formula(lm(y ~ x1+x2+x3+x4)), trace=1)

## Start:  AIC=209.36
## y ~ x1 + x2 + x3 + x4
## 
##        Df Sum of Sq    RSS    AIC
## - x2    1     27.31 3294.9 207.75
## - x3    1     35.92 3303.5 207.87
## <none>              3267.6 209.36
## - x4    1    708.49 3976.1 216.58
## - x1    1   1140.44 4408.0 221.43
## 
## Step:  AIC=207.75
## y ~ x1 + x3 + x4
## 
##        Df Sum of Sq    RSS    AIC
## - x3    1     33.49 3328.4 206.22
## <none>              3294.9 207.75
## - x4    1    794.74 4089.7 215.91
## - x1    1   1507.71 4802.6 223.46
## 
## Step:  AIC=206.22
## y ~ x1 + x4
## 
##        Df Sum of Sq    RSS    AIC
## <none>              3328.4 206.22
## - x4    1    855.16 4183.6 214.97
## - x1    1   2604.04 5932.4 231.39

summary(backward)

## 
## Call:
## lm(formula = y ~ x1 + x4)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.3107  -4.5149  -0.3636   6.5116  14.6796 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  56.0810     9.6026   5.840 5.79e-07 ***
## x1           -0.9493     0.1618  -5.867 5.29e-07 ***
## x4            1.4900     0.4432   3.362  0.00161 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.697 on 44 degrees of freedom
## Multiple R-squared:  0.5363, Adjusted R-squared:  0.5152 
## F-statistic: 25.44 on 2 and 44 DF,  p-value: 4.541e-08

Metode Forward

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

## Start:  AIC=238.35
## y ~ 1
## 
##        Df Sum of Sq    RSS    AIC
## + x1    1   2994.39 4183.6 214.97
## + x3    1   1543.29 5634.7 228.97
## + x4    1   1245.51 5932.4 231.39
## + x2    1    894.84 6283.1 234.09
## <none>              7178.0 238.34
## 
## Step:  AIC=214.97
## y ~ x1
## 
##        Df Sum of Sq    RSS    AIC
## + x4    1    855.16 3328.4 206.22
## <none>              4183.6 214.97
## + x2    1    110.83 4072.7 215.71
## + x3    1     93.91 4089.7 215.91
## 
## Step:  AIC=206.22
## y ~ x1 + x4
## 
##        Df Sum of Sq    RSS    AIC
## <none>              3328.4 206.22
## + x3    1    33.489 3294.9 207.75
## + x2    1    24.881 3303.5 207.87

summary(forward)

## 
## Call:
## lm(formula = y ~ x1 + x4)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.3107  -4.5149  -0.3636   6.5116  14.6796 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  56.0810     9.6026   5.840 5.79e-07 ***
## x1           -0.9493     0.1618  -5.867 5.29e-07 ***
## x4            1.4900     0.4432   3.362  0.00161 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.697 on 44 degrees of freedom
## Multiple R-squared:  0.5363, Adjusted R-squared:  0.5152 
## F-statistic: 25.44 on 2 and 44 DF,  p-value: 4.541e-08

Merode Stepwise

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

## Start:  AIC=238.35
## y ~ 1
## 
##        Df Sum of Sq    RSS    AIC
## + x1    1   2994.39 4183.6 214.97
## + x3    1   1543.29 5634.7 228.97
## + x4    1   1245.51 5932.4 231.39
## + x2    1    894.84 6283.1 234.09
## <none>              7178.0 238.34
## 
## Step:  AIC=214.97
## y ~ x1
## 
##        Df Sum of Sq    RSS    AIC
## + x4    1    855.16 3328.4 206.22
## <none>              4183.6 214.97
## + x2    1    110.83 4072.7 215.71
## + x3    1     93.91 4089.7 215.91
## - x1    1   2994.39 7178.0 238.34
## 
## Step:  AIC=206.22
## y ~ x1 + x4
## 
##        Df Sum of Sq    RSS    AIC
## <none>              3328.4 206.22
## + x3    1     33.49 3294.9 207.75
## + x2    1     24.88 3303.5 207.87
## - x4    1    855.16 4183.6 214.97
## - x1    1   2604.04 5932.4 231.39

summary(stepwise)

## 
## Call:
## lm(formula = y ~ x1 + x4)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.3107  -4.5149  -0.3636   6.5116  14.6796 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  56.0810     9.6026   5.840 5.79e-07 ***
## x1           -0.9493     0.1618  -5.867 5.29e-07 ***
## x4            1.4900     0.4432   3.362  0.00161 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.697 on 44 degrees of freedom
## Multiple R-squared:  0.5363, Adjusted R-squared:  0.5152 
## F-statistic: 25.44 on 2 and 44 DF,  p-value: 4.541e-08

Seleksi peubah menggunakan metode backward, forward, dan stepwise hasilnya sama, yaitu menggunakan peubah x1, x2, x3, dan x4. \[R-square= 53,63\%\]

Regresi Ridge (glmnet)

Packages

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

## Warning: package 'glmnet' was built under R version 4.2.3

## 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"    "lmtest"    "zoo"       "lubridate" "forcats"  
##  [7] "stringr"   "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"   
## [13] "ggplot2"   "tidyverse" "stats"     "graphics"  "grDevices" "utils"    
## [19] "datasets"  "methods"   "base"

Peubah

x <- cbind(x1,x2,x3,x4)
y <- swiss$Fertility

CV

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) 57.09085071
## x1          -0.84677851
## x2          -0.02192213
## x3           0.03234276
## x4           1.34357923

Dugaan Model Regresi Ridge:

\[Y=57.09085071-0.84677851X1 -0.02192213X2+0.03234276X3+1.34357923X4\]

Interpretasi:

• Intersep: tidak dapat diinterpretasikan nilai tersebut tidak bermakna

• X1 (Examination):-0.84677851, artinya jika Examination bertambah satu satuan, maka Fertility turun dengan rata-rata sebesar -0.84677851

• X2 (Agriculture): -0.02192213, artinya jika Agriculture bertambah satu satuan, maka Fertility turun dengan rata-rata sebesar -0.02192213

• X3 (catholic):0.03234276, artinya jika catholic bertambah satu satuan, maka Fertility naik dengan rata-rata sebesar 0.03234276

• X4 (Infant.Mortality): 1.34357923, artinya jika Infant.Mortality bertambah satu satuan, maka Fertility naik dengan rata-rata sebesar 1.34357923

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.5418479

Regresi Lasso (glmnet)

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) 57.88784204
## x1          -0.82472514
## x2           .         
## x3           0.01738481
## x4           1.26055196

Tidak ada koefisien yang ditampilkan untuk x2 prediktor karena regresi lasso mengecilkan koefisien hingga nol. Artinya, peubah tersebut dikeluarkan sepenuhnya dari model karena tidak cukup berpengaruh. Dugaan Model Regresi Lasso: Lasso: \[Y= 57.33506810- 0.8247251x1+ 0.01738481x3+ 1.26055196X4\]

R-Square Lasso

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

## [1] 0.5355897

Regresi Ridge (lmridge)

lmr <- lmridge(Fertility~Examination+Agriculture+Catholic+Infant.Mortality,data=swiss,scaling="centered")
plot(lmr)

summary(lmr)

## 
## Call:
## lmridge.default(formula = Fertility ~ Examination + Agriculture + 
##     Catholic + Infant.Mortality, data = swiss, scaling = "centered")
## 
## 
## Coefficients: for Ridge parameter K= 0 
##                  Estimate Estimate (Sc) StdErr (Sc) t-value (Sc) Pr(>|t|)    
## Intercept         59.6027       59.6027     11.0471       5.3953   <2e-16 ***
## Examination       -0.9680       -0.9680      0.2499      -3.8740   0.0004 ***
## Agriculture       -0.0476       -0.0476      0.0794      -0.5995   0.5520    
## Catholic           0.0261        0.0261      0.0380       0.6875   0.4955    
## Infant.Mortality   1.3960        1.3960      0.4572       3.0534   0.0039 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Ridge Summary
##        R2    adj-R2  DF ridge         F       AIC       BIC 
##   0.54480   0.51300   3.99998  12.86457 207.35820 395.71569 
## Ridge minimum MSE= 0.2792002 at K= 0 
## P-value for F-test ( 3.99998 , 43.00004 ) = 5.703745e-07 
## -------------------------------------------------------------------

Membandingkan R-Square

Regresi Klasik

summary(modelreg)$r.squared

## [1] 0.5447723

Regresi Ridge

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

## [1] 0.5418479

Regresi Lasso

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

## [1] 0.5355897

Membandingkan Residual Standard Error

Regresi Klasik

summary(modelreg)$sigma

## [1] 8.820436

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): 8.74522508529952"

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): 8.80475052204174"

Model Regresi Terbaik

Model	R-Square	Residual Standard Error
Regresi Klasik	\(54.47\%\)	\(8.820436\)
Regresi Ridge	\(54.18\%\)	\(8.74522508529952\)
Regresi Lasso	\(53.55\%\)	\(8.80475052204174\)

Model Regresi klasik memiliki R-Square paling besar dan model ridge RSE paling kecil sehingga model ridge tersebut merupakan model terbaik. Hal tersebut karena adanya multikolinearitas dan tipe data peubah yang digunakan merupakan numerik. Berikut model terbaik yang diperoleh:

\[Y=57.09085071-0.84677851X1 -0.02192213X2+0.03234276X3+1.34357923X4\]

Interpretasi Model Terbaik

• Intersep: nilai rataan peubah respon (Fertility) ketika seluruh peubah penjelas bernilai 0 sebesar 57.0908507

• X1 (Examination): -0.96805, artinya jika Examination bertambah satu satuan, maka Fertility turun dengan rata-rata sebesar -0.84677851

• X2 (Agriculture): -0.04759, artinya jika Agriculture bertambah satu satuan, maka Fertility turun dengan rata-rata sebesar -0.02192213

• X3 (catholic): 0.02611, artinya jika catholic bertambah satu satuan, maka Fertility naik dengan rata-rata sebesar 0.03234276

• X4 (Infant.Mortality): 1.39597, artinya jika Infant.Mortality bertambah satu satuan, maka Fertility naik dengan rata-rata sebesa 1.343579