Praktikum dan Soal Latihan

A. Praktikum

1. Diagnostik lengkap OLS pada data mtcars heteroskedastisitas (Breusch–Pagan), multikolinearitas (VIF), pengaruh (Cook’s D), normalitas (QQ-plot).

Ordinary Least Squares (OLS) adalah metode estimasi parameter yang paling umum digunakan dalam analisis regresi linier. Metode ini didasarkan pada prinsip bahwa garis regresi terbaik adalah garis yang memiliki jumlah kuadrat selisih antara nilai observasi dan nilai prediksi yang paling kecil. Berikut OLS pada data mtcars:

Simulasi Data

data(mtcars)
model_ols <- lm(mpg ~ wt + hp + disp, data = mtcars)
summary(model_ols)
## 
## Call:
## lm(formula = mpg ~ wt + hp + disp, data = mtcars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.891 -1.640 -0.172  1.061  5.861 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.105505   2.110815  17.579  < 2e-16 ***
## wt          -3.800891   1.066191  -3.565  0.00133 ** 
## hp          -0.031157   0.011436  -2.724  0.01097 *  
## disp        -0.000937   0.010350  -0.091  0.92851    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.639 on 28 degrees of freedom
## Multiple R-squared:  0.8268, Adjusted R-squared:  0.8083 
## F-statistic: 44.57 on 3 and 28 DF,  p-value: 8.65e-11

1) Heteroskedastisitas (Breusch–Pagan)

#UjiBreusch-Pagan untuk heteroskedastisitas
 lmtest::bptest(model_ols)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_ols
## BP = 0.9459, df = 3, p-value = 0.8143

2) Multikolinearitas (VIF)

 # Variance Inflation Factor (VIF)
 car::vif(model_ols)
##       wt       hp     disp 
## 4.844618 2.736633 7.324517

3) Pengaruh (Cook’s D)

cooks_d<- cooks.distance(model_ols)

which(cooks_d > 4/length(cooks_d))
## Chrysler Imperial    Toyota Corolla     Maserati Bora 
##                17                20                31
#Plot
plot(cooks_d, main="Cook's Distance", ylab="D", xlab="Observasi")
abline(h=4/length(cooks_d),lty=2)

4) Normalitas (QQ-plot)

#Normalitas Residual
residual<- rstandard(model_ols)
shapiro.test(residual) 
## 
##  Shapiro-Wilk normality test
## 
## data:  residual
## W = 0.92255, p-value = 0.02439
#QQ-plot
qqnorm(residual, col="darkgreen",main="Normal QQ Plot") 
qqline(residual, col="purple")

Pada hasil diatas Diagnostok menggunakan OLS uji dengan menggunakan data mtcars diperoleh bahwa:

  1. Pada uji heteroskedastisitas, jika p-value < 0.05, maka terdapat heteroskedastisitas (varian residual tidak konstan) dan jika p-value > 0.05, maka asumsi homoskedastisitas terpenuhi. Berdasarkan hasil uji diperoleh nilai p-value (0.8143) > 0.05 sehingga dapat disimpulkan bahwa tidak terdapat heteroskedastisitas atau tidak ada bukti kuat bahwa varian residual berubah seiring prediktor.

  2. Pada uji multikolinearitas, apabila VIF > 10 maka indikasi kuat multikolinearitas, jika VIF antara 5–10 maka perlu perhatian dan jika VIF < 5 maka aman. Berdasarkan hasil uji yang diperoleh bahwa nilai VIF untuk wt(4.844618) < 5 dan hp(2.736633) < 5 korelasi/aman serta disp(7.324517) < 10 nilai korelasi cukup tinggi, perlu perhatian tetapi belum masuk kategori berbahaya. Semua variabel (wt, hp, disp) masih dalam batas aman. Artinya, hasil regresi bisa dipercaya dan interpretasi koefisiennya masih stabil.

  3. Pada uji pengaruh (Cook’s D) apabila Cook’s D > 1 atau jauh lebih tinggi dari yang lain maka berpengaruh besar terhadap model. Berdasarkan plot yang diperoleh model cukup stabil karena sebagian besar observasi tidak berpengaruh besar.Namun, ada beberapa observasi ke- (17, 20, dan 31) yang perlu diteliti lebih lanjut apakah benar merupakan outlier atau data valid.

  4. Pada uji normalitas (QQ-plot) apabila terdapat titik-titik menyimpang dari garis diagonal maka mengindikasikan non-normalitas residual. Berdasarkan hasil uji Shapiro-Wilk Karena p-value < 0.05 artinya residual tidak sepenuhnya berdistribusi normal.Titik-titik residual mengikuti garis diagonal di bagian tengah (distribusi residual mendekati normal pada nilai rata-rata). Namun, di bagian ekor atas (right tail) terlihat penyimpangan yang cukup jelas mengindikasikan adanya heavy tail atau beberapa outlier pada residual.

Berdasarkan hasil uji heteroskedastisitas, multikolinearitas, pengaruh (Cook’s D) dan normalitas QQ Plot, dapat diberikan beberapa saran sebagai berikut:

  1. Jika ujinya signifikan pada uji heteroskedastisitas, pertimbangkan melakukan transformasi variabel atau gunakan Weighted Least Squares (WLS).

  2. jika nilai VIF yang diperoleh tinggi, pertimbangkan menghapus atau menggabungkan variabel, atau gunakan PCA / regularisasi (Ridge/Lasso).

  3. Setelah meninjau apakah observasi tersebut terdapat outlier, ada data yang salah dalam menginput atau memang datanya valid maka dapat mempertimbangkan penggunaan robust regression.

  4. Pada Uji Normalitas (QQ-plot) jika penyimpangan signifikan, pertimbangkan maka bisa dipertimbangkan transformasi data (log, square root, dll.) atau menggunakan metode robust regression.

2. Robust SE (HC3) vs WLS pada data simulasi heteroskedastik; bandingkan standar error dan p-value.

set.seed(123)

n <- 100
x <- runif(n, 0, 10)
error_sd<-0.6+0.2*x
e <- rnorm(n, mean = 0, sd = error_sd)  # varians error meningkat seiring x
y <- 5 + 2 * x + e

data_sim <- data.frame(y, x)

library(sandwich)
library(broom)
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, union
# OLS model
fit_ols <- lm(y ~ x, data = data_sim)

# WLS model dengan bobot invers varians (1 / x^2)
fit_wls <- lm(y ~ x, data = data_sim, weights = 1 / (x + 1)^2)  # transformasi agar bobot tidak ekstrem

#uji heteoskedastisitas
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(fit_ols)
## 
##  studentized Breusch-Pagan test
## 
## data:  fit_ols
## BP = 13.24, df = 1, p-value = 0.000274
# Robust SE (HC3) dari OLS
robust_se <- sqrt(diag(vcovHC(fit_ols, type = "HC3")))

# Hitung p-value robust secara manual
robust_p <- 2 * pt(-abs(coef(fit_ols) / robust_se), df = df.residual(fit_ols))

# Tidy output
tidy_ols <- tidy(fit_ols) %>%
  mutate(Robust_SE = robust_se,
         Robust_p = robust_p)

tidy_wls <- tidy(fit_wls) %>%
  rename(WLS_SE = std.error, WLS_p = p.value)

# Gabungkan
comp_stats <- left_join(
  tidy_ols %>% select(term, Robust_SE, Robust_p),
  tidy_wls %>% select(term, WLS_SE, WLS_p),
  by = "term"
)
print(comp_stats)
## # A tibble: 2 × 5
##   term        Robust_SE Robust_p WLS_SE    WLS_p
##   <chr>           <dbl>    <dbl>  <dbl>    <dbl>
## 1 (Intercept)    0.251  6.21e-36 0.132  2.89e-60
## 2 x              0.0629 3.55e-53 0.0486 2.89e-63

Berdasarkan hasil yang diperoleh nilai Standar Error Robust SE lebih besar karena mengakomodasi heteroskedastisitas sedangkan pada nilai Standar error WLS jauh lebih kecil karena bobot memperkecil pengaruh observasi dengan varians besar.Kedua p-value menunjukkan signifikansi yang sangat tinggi, tetapi p-value dari WLS lebih kecil secara signifikan dari pada p-value Robust SE.

3. Spline vs polinomial: modelkan mpg ~ wt dan bandingkan AIC/BIC + 5-fold CV.

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.5
## ✔ ggplot2   3.5.2     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.3.0
## ✔ purrr     1.1.0     ✔ tidyr     1.3.1
## ── 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
library(splines)
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
# Data
data(mtcars)

# Model polinomial orde 2
fit_poly2 <- lm(mpg ~ poly(wt, 2, raw = TRUE), data = mtcars)

# Model spline (3 knots)
fit_spline <- lm(mpg ~ bs(wt, df = 4), data = mtcars)

# Perbandingan AIC/BIC
model_comp <- tibble(
  Model = c("Polynomial_2", "Spline_df4"),
  AIC   = c(AIC(fit_poly2), AIC(fit_spline)),
  BIC   = c(BIC(fit_poly2), BIC(fit_spline))
)
print(model_comp)
## # A tibble: 2 × 3
##   Model          AIC   BIC
##   <chr>        <dbl> <dbl>
## 1 Polynomial_2  158.  164.
## 2 Spline_df4    162.  171.
# 5-fold CV RMSE 
set.seed(123)
train_control <- trainControl(method = "cv", number = 5)

cv_poly2 <- train(mpg ~ poly(wt, 2, raw = TRUE),
                  data = mtcars, method = "lm",
                  trControl = train_control)

cv_spline <- train(mpg ~ bs(wt, df = 4),
                   data = mtcars, method = "lm",
                   trControl = train_control)

cv_results <- tibble(
  Model = c("Polynomial_2", "Spline_df4"),
  RMSE  = c(mean(cv_poly2$resample$RMSE),
            mean(cv_spline$resample$RMSE))
)
print(cv_results)
## # A tibble: 2 × 2
##   Model         RMSE
##   <chr>        <dbl>
## 1 Polynomial_2  2.91
## 2 Spline_df4    3.08

a) AIC

AIC memberikan keseimbangan antara kecocokan model dengan kompleksitasnya. Model dengan nilai AIC yang lebih rendah dianggap lebih baik. Dalam kasus ini, model polinomial memiliki AIC yang lebih rendah (158 vs 162), menunjukkan bahwa model ini lebih baik dalam menyeimbangkan kecocokan data dengan jumlah parameternya dibandingkan dengan model spline.

b) BIC

Sama seperti AIC, BIC juga menghukum model yang lebih kompleks, tetapi dengan penalti yang lebih besar untuk setiap parameter tambahan. Model dengan nilai BIC yang lebih rendah juga dianggap lebih baik. Di sini, model polinomial memiliki BIC yang lebih rendah (164 vs 171), yang mengkonfirmasi temuan dari AIC bahwa model ini adalah pilihan yang lebih baik dari kedua model yang diuji.

c) RMSE

RMSE adalah ukuran standar untuk mengevaluasi seberapa baik model memprediksi nilai target. Ini mengukur rata-rata selisih kuadrat antara nilai prediksi model dan nilai sebenarnya. Nilai RMSE yang lebih rendah menunjukkan bahwa model memiliki performa prediksi yang lebih baik.Dari hasil ini, model polinomial (derajat 2) memiliki RMSE yang lebih rendah (2.91) dibandingkan dengan model spline (3.08). Ini berarti, secara rata-rata, kesalahan prediksi model polinomial lebih kecil daripada model spline. Dengan kata lain, model polinomial adalah pilihan yang lebih akurat untuk memprediksi nilai dalam kasus ini

4. AR(1) GLS: gunakan kode GLS di atas, ganti x dengan sinyal musiman (mis. sin(2pit/12)), interpretasikan̂𝜌 ; bandingkan SE OLS vs GLS. **

# Package yang diperlukan
library(nlme)
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
# 1. Buat data simulasi musiman
set.seed(123)
T <- 120
t <- 1:T
x <- sin(2*pi*t/12)               # sinyal musiman (periode 12)
u <- arima.sim(list(ar = 0.6), n = T) # error AR(1) dengan rho=0.6
y <- 2 + 1.5*x + u                # model: y = 2 + 1.5*sin + error

# 2. Estimasi OLS
ols_fit <- lm(y ~ x)

# 3. Estimasi GLS dengan korelasi AR(1)
gls_fit <- gls(y ~ x,
               correlation = corAR1(form = ~ t),
               method = "ML")

# 4. Bandingkan Standard Error
se_compare <- data.frame(
  Term     = names(coef(ols_fit)),
  OLS_SE   = coef(summary(ols_fit))[, "Std. Error"],
  GLS_SE   = coef(summary(gls_fit))[, "Std.Error"]
)
se_compare
##                    Term    OLS_SE    GLS_SE
## (Intercept) (Intercept) 0.1009513 0.1961441
## x                     x 0.1427668 0.2017475

Berdasarkan hasil diperoleh bahwa Koefisien (Intercept & slope) antara OLS dan GLS relatif sama. SE GLS > SE OLS, menunjukkan bahwa OLS terlalu optimis. SE GLS yang lebih besar menunjukkan bahwa jika kita ingin membuat kesimpulan tentang signifikansi koefisien. Untuk data time series dengan autokorelasi, GLS lebih tepat digunakan karena inference lebih dapat dipercaya.

B. Soal Latihan

1. Mengapa OLS tetap unbiased di bawah heteroskedastisitas, namun tidak efisien? Bagaimana robust SE mengatasi masalah inferensi?

Jawab:

OLS tetap unbiased dan tidak efisien

OLS estimator:

\[ \hat{\beta} = (X'X)^{-1}X'Y = (X'X)^{-1}X'(X\beta + \varepsilon) \\ \quad = \beta + (X'X)^{-1}X'\varepsilon \]

Jika asumsi eksogenitas tetap berlaku:

\[ E[\varepsilon \mid X] = 0, \\ \]

maka:

\[ E[\hat{\beta} \mid X] = \beta \]

Oleh karena itu OLS tetap tak bias (unbiased) meskipun error bersifat heteroskedastik. Karena Heteroskedastisitas hanya mengubah varians error, bukan nilai harapan error. Dalam asumsi klasik (homoskedastisitas, tidak ada autokorelasi), Gauss-Markov theorem menjamin OLS adalah BLUE (Best Linear Unbiased Estimator), yaitu punya varians terkecil di antara semua estimator linear tak bias. Apabila terdapat heteroskedastisitas maka varians error berbeda antar pengamatan:

\[ Var(\varepsilon \mid X) = \sigma_i^2 \]

Sehingga OLS berubah menjadi:

\[ Var(\hat{\beta} \mid X) = (X'X)^{-1} X' \Omega X (X'X)^{-1}, \quad \Omega = diag(\sigma_1^2, \ldots, \sigma_n^2) \]

Meskipun unbiased, estimator OLS menjadi tidak efisien (inefficient) di bawah heteroskedastisitas. Efisiensi mengacu pada seberapa kecil varians dari estimator. Estimator yang paling efisien adalah yang memiliki varians terkecil.Estimator seperti WLS (Weighted Least Squares) atau GLS lebih efisien bila Ω diketahui atau dapat diperkirakan.

Robust SE mengatasi masalah inferensi

Robust Standard Errors (SE) atau juga dikenal sebagai Standard Errors yang dikoreksi untuk Heteroskedastisitas (HCSE), adalah metode untuk mengatasi masalah inferensi yang disebabkan oleh heteroskedastisitas tanpa harus mengubah estimator OLS. Robust SE tidak mengubah koefisien regresi OLS itu sendiri hanya mengoreksi perhitungan varians dari koefisien tersebut. Rumus untuk varians robust SE secara umum adalah:

\[ Var_{robust}(\hat{\beta}) = (X'X)^{-1} \left( \sum_{i=1}^n e_i^2 x_i x_i' \right) (X'X)^{-1} \]

Dengan \(e_i\) adalah error dari observasi ke-\(i\) Dengan menggunakan rumus ini, robust SE memberikan estimasi varians yang konsisten (menjadi lebih akurat seiring dengan bertambahnya jumlah sampel), bahkan ketika varians error tidak konstan. Ini memungkinkan kita untuk melakukan uji hipotesis dan membangun interval kepercayaan yang valid dan dapat diandalkan, sehingga inferensi statistik menjadi benar.Sehingga Robust SE memungkinkan kita untuk tetap menggunakan estimator OLS yang unbiased, sambil memperbaiki masalah varians yang disebabkan oleh heteroskedastisitas, sehingga kita bisa membuat kesimpulan yang tepat dari model regresi.

2. Turunkan bentuk estimator WLS saat varian residual diketahui proporsional terhadap |𝑥|.

Jawab:

Weighted Least Squares (WLS) mengatasi heteroskedastisitas dengan memberikan bobot yang berbeda pada setiap observasi. Bobot ini ditentukan oleh varians residual, sehingga observasi dengan varians kecil (lebih reliabel) diberikan bobot yang lebih besar. Penuruan betuk estimator WLS saat \(\sigma_i^2 = \sigma^2 |x_i|\) sebagai berikut:

Model regresi linear: \[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i \] \[ \mathbb{E}[\varepsilon_i \mid x_i] = 0, \quad \mathrm{Var}(\varepsilon_i \mid x_i) = \sigma^2 |x_i| \] 1) Meminimalkan fungsi berikut:

\[ S(\beta) = \sum_{i=1}^n w_i \varepsilon_i^2 \]

Dengan:

\(w_i\) = bobot

\(\varepsilon_i\) = residual OLS (\(\varepsilon_i = y_i - x_i'\beta\))

Bobot yang optimal untuk WLS adalah invers dari residual, yaitu \(w_i^2 = \sigma_i^2\). Karena \(\sigma_i^2 = \sigma^2 |x_i|\), sehingga diperoleh bobotnya sebagai berikut:

\[ w_i = \frac{1}{\sigma^2 |x_i|} \]

Berdasarkan persamaan diatas karena \(\sigma^2\) merupakan konstanta dan dapat di abaikan dalam proses minimasi karena hanya akan mengalikan seluruh fungsi dengan sebuah konstanta. Oleh karena itu dapat menggunakan bobot sederhana yaitu:

\[ w_i = \frac{1}{|x_i|} \]

Sehingga, minimasi menjadi:

\[ S(\beta) = \sum_{i=1}^n \frac{1}{|x_i|} \left( y_i - x_i' \beta \right)^2 \]

2) Mengubah persamaan kedalam bentuk matriks

Untuk mempermudah penurunan, kita ubah persamaan ke dalam bentuk matriks.

\[ S(\beta) = (y - X\beta)' W (y - X\beta) \] y : Vektor observasi n×1

X : Matriks desain n×k

β : Vektor koefisien k×1

Ω : Matriks kovariansi residual n×n. Dalam kasus WLS, matriks ini adalah matriks diagonal dengan elemen diagonal \(\sigma_i^2\)

Dengan W adalah matriks bobot diagonal n x n, dengan elemen diagonal \(w_i = \frac{1}{|x_i|}\).

\[ W = \begin{pmatrix} \dfrac{1}{|x_1|} & 0 & \cdots & 0 \\ 0 & \dfrac{1}{|x_2|} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \dfrac{1}{|x_n|} \end{pmatrix} \] 3) Meminimalkan fungsi kuadratik

\[ \frac{\partial S}{\partial \beta} = \frac{\partial}{\partial \beta} \big[ (y - X\beta)' W (y - X\beta) \big] = 0 \]

\[ -2X'W(y - X\beta) = 0 \]

\[ -2X'Wy + 2X'WX\beta = 0 \]

\[ X'WX\beta = X'Wy \] \[ \hat{\beta}_{WLS} = (X'WX)^{-1} X'Wy \] \[ \hat{\beta}_{WLS} = \Bigg( X' \begin{pmatrix} \dfrac{1}{|x_1|} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \dfrac{1}{|x_n|} \end{pmatrix} X \Bigg)^{-1} X' \begin{pmatrix} \dfrac{1}{|x_1|} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \dfrac{1}{|x_n|} \end{pmatrix} y \]

Diatas merupakan hasil akhir dari estimator WLS ketika varians residual proporsional pada \({|x_1|}\). setiap observasi \(y_i\) dan \(x_i\) dibagi dengan \(\sqrt{|x_1|}\) sebelum menjalankan regresi OLS

3. Bandingkan Ridge dan Lasso dari sudut pandang bias–variance dan seleksi variabel.

Jawab:

Ridge dan Lasso adalah dua metode regresi yang digunakan untuk mengatasi overfitting dan multikolinearitas. Keduanya menambahkan penalti pada model, tetapi dengan cara yang berbeda, yang memengaruhi trade-off bias-variance dan kemampuan seleksi variabel mereka.

Simulasi data

set.seed(123)

# Data simulasi: banyak prediktor, hanya sebagian relevan
n <- 100
p <- 10
X <- matrix(rnorm(n*p), n, p)
beta_true <- c(3, -2, 0, 0, 1.5, rep(0, 5))
y <- X %*% beta_true + rnorm(n, sd=2)

# Split data: train/test
train <- 1:71
test <- 71:100
X_train <- X[train, ]
y_train <- y[train]
X_test  <- X[test, ]
y_test  <- y[test]

a) Ridge regression

library(glmnet)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## Loaded glmnet 4.1-10
# Ridge (alpha = 0)
ridge_fit <- cv.glmnet(X_train, y_train, alpha=0)
coef(ridge_fit, s="lambda.min")
## 11 x 1 sparse Matrix of class "dgCMatrix"
##              lambda.min
## (Intercept)  0.37541014
## V1           3.13145344
## V2          -1.90510430
## V3          -0.31233744
## V4           0.34527238
## V5           1.58546303
## V6           0.14602918
## V7           0.08855246
## V8           0.06448136
## V9          -0.05729561
## V10          0.08232586
y_pred_ridge <- predict(ridge_fit, X_test, s="lambda.min")
ridge_mse <- mean((y_test - y_pred_ridge)^2)
ridge_mse
## [1] 3.708154

Berdasarkan hasil yang diperoleh dapat dilihat bahwa tidak ada koefisien yang menjadi nol. Ridge hanya mengecilkan nilai koefisien.

b) Lasso regression

# Lasso (alpha = 1)
lasso_fit <- cv.glmnet(X_train, y_train, alpha=1)
coef(lasso_fit, s="lambda.min")
## 11 x 1 sparse Matrix of class "dgCMatrix"
##             lambda.min
## (Intercept)  0.4011800
## V1           3.1189010
## V2          -1.8295334
## V3          -0.1344854
## V4           0.1940877
## V5           1.5263127
## V6           .        
## V7           .        
## V8           .        
## V9           .        
## V10          .
y_pred_lasso <- predict(lasso_fit, X_test, s="lambda.min")
lasso_mse <- mean((y_test - y_pred_lasso)^2)
lasso_mse
## [1] 3.912388

Berdasarkan hasil yang diperoleh bahwa beberapa koefisien menjadi nol. Lasso secara otomatis melakukan seleksi variabel dengan menghilangkan prediktor yang kurang relevan

c) Menampilkan Plot

# Plot koefisien untuk Ridge dan Lasso
plot(ridge_fit, main = "Ridge Regression Coefficients")

plot(lasso_fit, main = "Lasso Regression Coefficients")

Berdasarkan hasil simulasi diatas didapatkan perbandingan Ridge dan Lasso dari sudut pandang bias–variance dan seleksi variabel sebagai berikut:

  • Ridge Mengurangi varians secara signifikan, dengan mengecilkan koefisien secara bersamaan, Ridge membuat model kurang sensitif terhadap fluktuasi data pelatihan tetapi menjaga semua variabel, lebih baik untuk stabilitas prediksi saat semua fitur relevan.

  • Lasso Mengurangi varians secara signifikan dengan mengecilkan koefisien menjadi nol dan melakukan seleksi variabel (untuk membuat model lebih sederhana dan mengurangi varians lebih jauh) dapat membuat bias lebih besar dibanding Ridge. Lasso lebih baik untu kmelakukan seleksi variabel.

4. Hitungan: Diberikan matriks 𝑋 dan vektor 𝑦, hitung ̂ 𝛽 OLS serta 𝐻; identifikasi leverage maksimum.

Jawab:

Diberikan matriks X dan vektor y sebagai berikut:

X = \[\begin{pmatrix}1 & 2 \\1 & 4 \\1 & 6\end{pmatrix}\] , y= \[\begin{pmatrix}5 \\7 \\9\end{pmatrix}\]

1) Menghitung \(\hat{\beta}_{OLS}\) menggunakan rumus:

\[ \hat{\beta} = (X'X)^{-1}X'y \]

# Definisikan matriks X dan vektor y
X <- matrix(c(1, 1, 1, 2, 4, 6), ncol = 2)
y <- as.matrix(c(4, 6, 8))
# 1. (X'X)
XtX <- t(X) %*% X
# 2. Invers dari (X'X)
XtX_inv <- solve(XtX)
# 3. (X'y)
Xty <- t(X) %*% y
# 4. Estimator beta_hat
beta_hat <- XtX_inv %*% Xty
# Menampilkan hasilnya
print(beta_hat)
##      [,1]
## [1,]    2
## [2,]    1

Jadi, Koefisien \(\hat{\beta}_{OLS}\) adalah \(\hat{\beta}_0=1=2\) dan \(\hat{\beta}_1=1\)

2) Menghitung matriks hat (H) dengan menggunakan rumus:

\[ H = X (X'X)^{-1} X' \]

# Hitung matriks H
H <- X %*% XtX_inv %*% t(X)

# Tampilkan matriks H
print(H)
##            [,1]      [,2]       [,3]
## [1,]  0.8333333 0.3333333 -0.1666667
## [2,]  0.3333333 0.3333333  0.3333333
## [3,] -0.1666667 0.3333333  0.8333333

3) mengidentifikasi leverege maksimum

Leverage adalah elemen diagonal dari matriks Hat (H), yang mengukur pengaruh setiap observasi terhadap garis regresi.

# Ambil elemen diagonal dari matriks H
leverage <- diag(H)
# Temukan nilai maksimum leverage
max_leverage <- max(leverage)
# Tampilkan nilai leverage untuk setiap observasi
print(leverage)
## [1] 0.8333333 0.3333333 0.8333333
# Tampilkan nilai leverage maksimum
cat("Nilai leverage maksimum adalah: ", max_leverage, "\n")
## Nilai leverage maksimum adalah:  0.8333333

Sehingga diperoleh nilai leverege maksimumnya adalah 0.83. hasil tersebut menunjukkan bahwa pada obervasi 3 memiliki pengaruh pada model regresi dibandingkan observasi pertama dan kedua.

5. Pilihan Ganda:

a) VIF tinggi menunjukkan?

b) Cook’s D digunakan untuk?

c) Dalam AR(1), kovarians antara \(\varepsilon_t\) dan \(\varepsilon_{t-k}\) adalah?

Jawab:

a) VIF (Variance Inflation Factor) yang tinggi menunjukkan adanya multikolinearitas yang serius. Multikolinearitas terjadi ketika dua atau lebih variabel independen dalam model regresi memiliki korelasi yang tinggi satu sama lain. VIF mengukur seberapa besar varians koefisien regresi terinflasi akibat korelasi ini, yang dapat menyebabkan estimator OLS menjadi tidak stabil dan tidak dapat diandalkan.

b) Cook’s D digunakan untuk mengidentifikasi observasi yang berpengaruh (influential observations). Observasi yang berpengaruh adalah titik data yang, jika dihapus, akan secara signifikan mengubah hasil model regresi (seperti koefisien regresi). Cook’s D mengukur pengaruh gabungan dari leverage dan residual suatu observasi. Nilai Cook’s D yang besar menunjukkan bahwa observasi tersebut mungkin merupakan pencilan yang berpengaruh dan perlu diselidiki lebih lanjut.

c) Model AR(1):

\[ \varepsilon_t = \rho \varepsilon_{t-1} + u_t, \quad u_t \sim i.i.d(0, \sigma_u^2) \]

Dengan Varians dari \(\varepsilon_t\):

\[ \mathrm{Var}(\varepsilon_t) = \frac{\sigma_u^2}{1 - \rho^2} \]

Menghitung kovarians dengan menggunakan sifat stasioneritas dan rekursi:

1) pada \(k =1\): \[ \mathrm{Cov}(\varepsilon_t, \varepsilon_{t-1}) = \mathbb{E}[\varepsilon_t \varepsilon_{t-1}] = \mathbb{E}[(\rho \varepsilon_{t-1} + u_t)\varepsilon_{t-1}] \] Dikarenakan \(u_t\) tidak berkorelasi dengan \(\varepsilon_{t-k}\):

\[ \mathrm{Cov}(\varepsilon_t, \varepsilon_{t-1}) = \rho \, \mathrm{Var}(\varepsilon_{t-1}) = \rho \cdot \frac{\sigma_u^2}{1 - \rho^2}. \]

2) Pada \(k = 2\):

\[ \mathrm{Cov}(\varepsilon_t, \varepsilon_{t-2}) = \rho \cdot \mathrm{Cov}(\varepsilon_{t-1}, \varepsilon_{t-2}) = \rho \cdot \left( \rho \, \frac{\sigma_u^2}{1 - \rho^2} \right) = \rho^2 \, \frac{\sigma_u^2}{1 - \rho^2}. \]

3) Secara umum kovarians untuk \(k \geq 0 :\)

\[ \mathrm{Cov}(\varepsilon_t, \varepsilon_{t-k}) = \rho^k \cdot \sigma_\varepsilon^2 = \rho^k \cdot \frac{\sigma_u^2}{1 - \rho^2} \]

Berdasarkan persamaan diatas diketahui bahwa kovarian menurun secara geometrik dengan \(\rho^k\). Jika \(\rho > 0,\) kovarianas positif dan makin kecil untuk lag lebih besar dan \(\rho < 0\) maka tanda kovariand bergantian positif-negatif.