Tugas Individu Praktikum PSD

Data

Menggunakan dataset yang terdapat di R, yaitu dataset mtcars. Dataset mtcars merupakan kumpulan data bawaan di R yang berisi pengukuran pada 11 atribut berbeda untuk 32 mobil berbeda. Peubah yang digunakan:

Peubah

y = mpg (Miles/(US) gallon)

x1 = cyl (Number of cylinders)

x2 = disp (Displacement (cu.in.))

x3 = hp (Gross horsepower)

x4 = wt (Weight (1000 lbs))

data("mtcars")
y <- mtcars$mpg
x1 <- mtcars$cyl
x2 <- mtcars$disp
x3 <- mtcars$hp
x4 <- mtcars$wt

Ekplorasi

# Sebaran peubah y (mpg)
hist(y, col = "pink4")

boxplot(y, col = "pink4")

Matriks Korelasi

library(tidyverse)

## ── 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.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 = "pink4") +
  theme_minimal() +
  labs(title = "Hubungan Antar Peubah",
       subtitle = "Peubah respon : Production.",
       x = NULL, y = NULL)

Dari eksplorasi menggunakan plot di atas, terlihat bahwa peubah x1, x2, x3, dan x4 memiliki nilai korelasi yang 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 
## -4.0562 -1.4636 -0.4281  1.2854  5.8269 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 40.82854    2.75747  14.807 1.76e-14 ***
## x1          -1.29332    0.65588  -1.972 0.058947 .  
## x2           0.01160    0.01173   0.989 0.331386    
## x3          -0.02054    0.01215  -1.691 0.102379    
## x4          -3.85390    1.01547  -3.795 0.000759 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.513 on 27 degrees of freedom
## Multiple R-squared:  0.8486, Adjusted R-squared:  0.8262 
## F-statistic: 37.84 on 4 and 27 DF,  p-value: 1.061e-10

Didapatkan model regresi klasik dengan R-squared sebesar 84.86%. 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 x4 (wt).

\[Y=40.82854-1.29332X_1+0.01160X_2-0.02054X_3-3.85390X_4\]

Multikolinearitas

car::vif(modelreg)

##        x1        x2        x3        x4 
##  6.737707 10.373286  3.405983  4.848016

Metode VIF digunakan untuk mendeteksi multikolinearitas yang terjadi ketika VIF > 10. Diperoleh bahwa terdapat multikolinearitas pada peubah x2 (disp) yang memiliki VIF = 10.373286 > 10. Sehingga akan dicoba membuat model tanpa peubah x2 dan model ditangani menggunakan metode penyusutan, yaitu Regresi Ridge dan Lasso.

Model tanpa x2

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

## 
## Call:
## lm(formula = y ~ x1 + x3 + x4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9290 -1.5598 -0.5311  1.1850  5.8986 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 38.75179    1.78686  21.687  < 2e-16 ***
## x1          -0.94162    0.55092  -1.709 0.098480 .  
## x3          -0.01804    0.01188  -1.519 0.140015    
## x4          -3.16697    0.74058  -4.276 0.000199 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.512 on 28 degrees of freedom
## Multiple R-squared:  0.8431, Adjusted R-squared:  0.8263 
## F-statistic: 50.17 on 3 and 28 DF,  p-value: 2.184e-11

car::vif(modelreg2)

##       x1       x3       x4 
## 4.757456 3.258481 2.580486

Dugaan Model Regresi Klasik: \[Y=38.75179-0.94162X_1-0.01804X_3-3.16697X_4\]

Diperoleh bahwa VIF > 10 sehingga tidak terdapat multikolinearitas pada model regresi klasik tanpa peubah x2 (disp).

Pengujian Asumsi

1. Uji Multikolinearitas

car::vif(modelreg2)

##       x1       x3       x4 
## 4.757456 3.258481 2.580486

Nilai VIF pada kedua peubah penjelas < 10. Sehingga tidak terjadi multikolinearitas pada kedua peubah penjelas tersebut.

2. Sisaan menyebar normal

plot(modelreg2, 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(modelreg2), mu = 0,)

## 
##  One Sample t-test
## 
## data:  resid(modelreg2)
## t = 2.9189e-16, df = 31, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.8605801  0.8605801
## sample estimates:
##    mean of x 
## 1.231654e-16

\(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(modelreg2)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelreg2
## BP = 2.9351, df = 3, p-value = 0.4017

\(p-value=0.4017 > 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)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

dwtest(modelreg2)

## 
##  Durbin-Watson test
## 
## data:  modelreg2
## DW = 1.6441, p-value = 0.1002
## alternative hypothesis: true autocorrelation is greater than 0

#ACF dan PACF identifikasi autokorelasi
sisaan = modelreg2$residuals
par(mfrow = c(1,2))
acf(sisaan)
pacf(sisaan)

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

Seleksi Peubah

Metode Backward

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

## Start:  AIC=62.66
## y ~ x1 + x3 + x4
## 
##        Df Sum of Sq    RSS    AIC
## <none>              176.62 62.665
## - x3    1    14.551 191.17 63.198
## - x1    1    18.427 195.05 63.840
## - x4    1   115.354 291.98 76.750

summary(backward)

## 
## Call:
## lm(formula = y ~ x1 + x3 + x4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9290 -1.5598 -0.5311  1.1850  5.8986 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 38.75179    1.78686  21.687  < 2e-16 ***
## x1          -0.94162    0.55092  -1.709 0.098480 .  
## x3          -0.01804    0.01188  -1.519 0.140015    
## x4          -3.16697    0.74058  -4.276 0.000199 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.512 on 28 degrees of freedom
## Multiple R-squared:  0.8431, Adjusted R-squared:  0.8263 
## F-statistic: 50.17 on 3 and 28 DF,  p-value: 2.184e-11

Metode Forward

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

## Start:  AIC=115.94
## y ~ 1
## 
##        Df Sum of Sq     RSS     AIC
## + x4    1    847.73  278.32  73.217
## + x1    1    817.71  308.33  76.494
## + x3    1    678.37  447.67  88.427
## <none>              1126.05 115.943
## 
## Step:  AIC=73.22
## y ~ x4
## 
##        Df Sum of Sq    RSS    AIC
## + x1    1    87.150 191.17 63.198
## + x3    1    83.274 195.05 63.840
## <none>              278.32 73.217
## 
## Step:  AIC=63.2
## y ~ x4 + x1
## 
##        Df Sum of Sq    RSS    AIC
## + x3    1    14.551 176.62 62.665
## <none>              191.17 63.198
## 
## Step:  AIC=62.66
## y ~ x4 + x1 + x3

summary(forward)

## 
## Call:
## lm(formula = y ~ x4 + x1 + x3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9290 -1.5598 -0.5311  1.1850  5.8986 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 38.75179    1.78686  21.687  < 2e-16 ***
## x4          -3.16697    0.74058  -4.276 0.000199 ***
## x1          -0.94162    0.55092  -1.709 0.098480 .  
## x3          -0.01804    0.01188  -1.519 0.140015    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.512 on 28 degrees of freedom
## Multiple R-squared:  0.8431, Adjusted R-squared:  0.8263 
## F-statistic: 50.17 on 3 and 28 DF,  p-value: 2.184e-11

Merode Stepwise

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

## Start:  AIC=115.94
## y ~ 1
## 
##        Df Sum of Sq     RSS     AIC
## + x4    1    847.73  278.32  73.217
## + x1    1    817.71  308.33  76.494
## + x3    1    678.37  447.67  88.427
## <none>              1126.05 115.943
## 
## Step:  AIC=73.22
## y ~ x4
## 
##        Df Sum of Sq     RSS     AIC
## + x1    1     87.15  191.17  63.198
## + x3    1     83.27  195.05  63.840
## <none>               278.32  73.217
## - x4    1    847.73 1126.05 115.943
## 
## Step:  AIC=63.2
## y ~ x4 + x1
## 
##        Df Sum of Sq    RSS    AIC
## + x3    1    14.551 176.62 62.665
## <none>              191.17 63.198
## - x1    1    87.150 278.32 73.217
## - x4    1   117.162 308.33 76.494
## 
## Step:  AIC=62.66
## y ~ x4 + x1 + x3
## 
##        Df Sum of Sq    RSS    AIC
## <none>              176.62 62.665
## - x3    1    14.551 191.17 63.198
## - x1    1    18.427 195.05 63.840
## - x4    1   115.354 291.98 76.750

summary(stepwise)

## 
## Call:
## lm(formula = y ~ x4 + x1 + x3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9290 -1.5598 -0.5311  1.1850  5.8986 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 38.75179    1.78686  21.687  < 2e-16 ***
## x4          -3.16697    0.74058  -4.276 0.000199 ***
## x1          -0.94162    0.55092  -1.709 0.098480 .  
## x3          -0.01804    0.01188  -1.519 0.140015    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.512 on 28 degrees of freedom
## Multiple R-squared:  0.8431, Adjusted R-squared:  0.8263 
## F-statistic: 50.17 on 3 and 28 DF,  p-value: 2.184e-11

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

Regresi Ridge (glmnet)

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"    "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 <- mtcars$mpg

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) 37.703388902
## x1          -0.921014161
## x2          -0.001918732
## x3          -0.018202005
## x4          -2.735655859

Dugaan Model Regresi Ridge:

\[Y=37.703388902-0.921014161X_1-0.001918732X_2-0.018202005X_3-2.735655859X_4\]

Interpretasi:

Intersep: nilai rataan peubah respon (mpg) ketika seluruh peubah penjelas bernilai 0 sebesar 37.703388902
X1 (cyl): -0.921014161, artinya jika cyl (Number of cylinders) bertambah satu satuan, maka mpg turun dengan rata-rata sebesar 0.921014161
X2 (dips): -0.001918732, artinya jika disp (Displacement (cu.in.)) bertambah satu satuan, maka mpg turun dengan rata-rata sebesar 0.001918732
X3 (hp): -0.018202005, artinya jika hp (Gross horsepower) bertambah satu satuan, maka mpg turun dengan rata-rata sebesar 0.018202005
X4 (wt): -2.735655859, artinya jika wt (Weight (1000 lbs)) bertambah satu satuan, maka mpg turun dengan rata-rata sebesar 2.735655859

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

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) 38.35869594
## x1          -0.93185789
## x2           .         
## x3          -0.01715355
## x4          -3.10388967

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.

R-Square Lasso

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

## [1] 0.8427049

Regresi Ridge (lmridge)

lmr <- lmridge(mpg~cyl+disp+hp+wt,data=mtcars,scaling="centered")
plot(lmr)

summary(lmr)

## 
## Call:
## lmridge.default(formula = mpg ~ cyl + disp + hp + wt, data = mtcars, 
##     scaling = "centered")
## 
## 
## Coefficients: for Ridge parameter K= 0 
##           Estimate Estimate (Sc) StdErr (Sc) t-value (Sc) Pr(>|t|)    
## Intercept  40.8285       40.8285      6.1179       6.6736   <2e-16 ***
## cyl        -1.2933       -1.2933      0.6441      -2.0081   0.0547 .  
## disp        0.0116        0.0116      0.0115       1.0073   0.3227    
## hp         -0.0205       -0.0205      0.0119      -1.7219   0.0965 .  
## wt         -3.8539       -3.8539      0.9972      -3.8648   0.0006 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Ridge Summary
##        R2    adj-R2  DF ridge         F       AIC       BIC 
##   0.84860   0.83240   4.00000  39.24578  61.52554 178.29203 
## Ridge minimum MSE= 1.409444 at K= 0 
## P-value for F-test ( 4 , 28.00001 ) = 4.268799e-11 
## -------------------------------------------------------------------

Membandingkan R-Square

Regresi Klasik

summary(modelreg2)$r.squared

## [1] 0.84315

Regresi Ridge

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

## [1] 0.8395815

Regresi Lasso

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

## [1] 0.8427049

Membandingkan Residual Standard Error

Regresi Klasik

summary(modelreg2)$sigma

## [1] 2.511548

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.53995779147227"

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.515109350162"

Model Regresi Terbaik

Model	R-Square	Residual Standard Error
Regresi Klasik	\(84.31 \%\)	\(2.511548\)
Regresi Ridge	\(83.95\%\)	\(2.53995779147227\)
Regresi Lasso	\(84.27\%\)	\(2.515109350162\)

Model Regresi Klasik memiliki R-Square paling besar dan RSE paling kecil. Jadi, model terbaik adalah Model Regresi Klasik.

\[Y=38.75179-0.94162X_1-0.01804X_3-3.16697X_4\]

Interpretasi Model Terbaik

Intersep: nilai rataan peubah respon (mpg) ketika seluruh peubah penjelas bernilai 0 sebesar 38.75179
X1 (cyl): -0.94162, artinya jika cyl (Number of cylinders) bertambah satu satuan, maka mpg turun dengan rata-rata sebesar 0.94162
X3 (hp): -0.01804 , artinya jika hp (Gross horsepower) bertambah satu satuan, maka mpg turun dengan rata-rata sebesar 0.01804
X4 (wt): -3.16697, artinya jika wt (Weight (1000 lbs)) bertambah satu satuan, maka mpg turun dengan rata-rata sebesa 3.16697