Pada contoh ini akan digunakan data yang berasal dari Package:
ISLR2. Data dan contoh penggunaan Ridge Regression ini
diperoleh dari Buku “Introduction to Statistics Learning: with
Application in R” karya Gareth James, Trevor Hastie, Daniela Witten, dan
Robert Tibshirani (https://www.statlearning.com/).
Memanggil packages yang diperlukan
library(ISLR2)## Warning: package 'ISLR2' was built under R version 4.4.3library(glmnet)## Warning: package 'glmnet' was built under R version 4.4.3## Loading required package: Matrix## Loaded glmnet 4.1-8Package: ISLR2 diperlukan untuk memanggil dataset yang
akan digunakan dalam contoh ini. Package: glmnet diperlukan
untuk melakukan analisis dengan penalized likelihood.
Penduga \(\beta\) dalam Ridge Regression diberikan oleh: \[\hat{\mathbf{β}}_{Ridge} =(\mathbf{X}^{T}\mathbf{X}+λ\mathbf{I})^{-1}\mathbf{X}^{T}\mathbf{Y}\]
Dalam hal ini akan digunakan glmnet() dari
Package: glmnet untuk menentukan \(\hat{\mathbf{β}}_{Ridge}\).
Data yang akan digunakan adalah data terkait gaji 322 pemain baseball dengan peubah prediktor sebanyak 19.
data_ridge <- Hitters
str(data_ridge)## 'data.frame': 322 obs. of 20 variables:
## $ AtBat : int 293 315 479 496 321 594 185 298 323 401 ...
## $ Hits : int 66 81 130 141 87 169 37 73 81 92 ...
## $ HmRun : int 1 7 18 20 10 4 1 0 6 17 ...
## $ Runs : int 30 24 66 65 39 74 23 24 26 49 ...
## $ RBI : int 29 38 72 78 42 51 8 24 32 66 ...
## $ Walks : int 14 39 76 37 30 35 21 7 8 65 ...
## $ Years : int 1 14 3 11 2 11 2 3 2 13 ...
## $ CAtBat : int 293 3449 1624 5628 396 4408 214 509 341 5206 ...
## $ CHits : int 66 835 457 1575 101 1133 42 108 86 1332 ...
## $ CHmRun : int 1 69 63 225 12 19 1 0 6 253 ...
## $ CRuns : int 30 321 224 828 48 501 30 41 32 784 ...
## $ CRBI : int 29 414 266 838 46 336 9 37 34 890 ...
## $ CWalks : int 14 375 263 354 33 194 24 12 8 866 ...
## $ League : Factor w/ 2 levels "A","N": 1 2 1 2 2 1 2 1 2 1 ...
## $ Division : Factor w/ 2 levels "E","W": 1 2 2 1 1 2 1 2 2 1 ...
## $ PutOuts : int 446 632 880 200 805 282 76 121 143 0 ...
## $ Assists : int 33 43 82 11 40 421 127 283 290 0 ...
## $ Errors : int 20 10 14 3 4 25 7 9 19 0 ...
## $ Salary : num NA 475 480 500 91.5 750 70 100 75 1100 ...
## $ NewLeague: Factor w/ 2 levels "A","N": 1 2 1 2 2 1 1 1 2 1 ...summary(data_ridge)## AtBat Hits HmRun Runs
## Min. : 16.0 Min. : 1 Min. : 0.00 Min. : 0.00
## 1st Qu.:255.2 1st Qu.: 64 1st Qu.: 4.00 1st Qu.: 30.25
## Median :379.5 Median : 96 Median : 8.00 Median : 48.00
## Mean :380.9 Mean :101 Mean :10.77 Mean : 50.91
## 3rd Qu.:512.0 3rd Qu.:137 3rd Qu.:16.00 3rd Qu.: 69.00
## Max. :687.0 Max. :238 Max. :40.00 Max. :130.00
##
## RBI Walks Years CAtBat
## Min. : 0.00 Min. : 0.00 Min. : 1.000 Min. : 19.0
## 1st Qu.: 28.00 1st Qu.: 22.00 1st Qu.: 4.000 1st Qu.: 816.8
## Median : 44.00 Median : 35.00 Median : 6.000 Median : 1928.0
## Mean : 48.03 Mean : 38.74 Mean : 7.444 Mean : 2648.7
## 3rd Qu.: 64.75 3rd Qu.: 53.00 3rd Qu.:11.000 3rd Qu.: 3924.2
## Max. :121.00 Max. :105.00 Max. :24.000 Max. :14053.0
##
## CHits CHmRun CRuns CRBI
## Min. : 4.0 Min. : 0.00 Min. : 1.0 Min. : 0.00
## 1st Qu.: 209.0 1st Qu.: 14.00 1st Qu.: 100.2 1st Qu.: 88.75
## Median : 508.0 Median : 37.50 Median : 247.0 Median : 220.50
## Mean : 717.6 Mean : 69.49 Mean : 358.8 Mean : 330.12
## 3rd Qu.:1059.2 3rd Qu.: 90.00 3rd Qu.: 526.2 3rd Qu.: 426.25
## Max. :4256.0 Max. :548.00 Max. :2165.0 Max. :1659.00
##
## CWalks League Division PutOuts Assists
## Min. : 0.00 A:175 E:157 Min. : 0.0 Min. : 0.0
## 1st Qu.: 67.25 N:147 W:165 1st Qu.: 109.2 1st Qu.: 7.0
## Median : 170.50 Median : 212.0 Median : 39.5
## Mean : 260.24 Mean : 288.9 Mean :106.9
## 3rd Qu.: 339.25 3rd Qu.: 325.0 3rd Qu.:166.0
## Max. :1566.00 Max. :1378.0 Max. :492.0
##
## Errors Salary NewLeague
## Min. : 0.00 Min. : 67.5 A:176
## 1st Qu.: 3.00 1st Qu.: 190.0 N:146
## Median : 6.00 Median : 425.0
## Mean : 8.04 Mean : 535.9
## 3rd Qu.:11.00 3rd Qu.: 750.0
## Max. :32.00 Max. :2460.0
## NA's :59dim(data_ridge)## [1] 322 20Data yang digunakan memiliki \(n=322\) dengan jumlah peubah bebas adalah
\(p=19\). Dalam hal ini peubah
Salary akan ditetapkan sebagai peubah respon.
Dari hasil str() diketahui terdapat nilai
NA pada data, sehingga perlu dilakukan perlakuan untuk
nilai yang hilang tersebut.
sum(is.na(data_ridge$Salary))## [1] 59data_ridge <- na.omit(data_ridge)
sum(is.na(data_ridge$Salary))## [1] 0dim(data_ridge)## [1] 263 20Dengan mengeluarkan missing value dari data, maka dataset yang digunakan menjadi berdimensi \((262\times20)\).
Setelah missing value dihilangkan selanjutnya akan dibentuk matriks \(\mathbf{X}\) dan \(\mathbf{Y}\).
x <- model.matrix(Salary ~., data = data_ridge)[,-1]
head(x,6)## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun
## -Alan Ashby 315 81 7 24 38 39 14 3449 835 69
## -Alvin Davis 479 130 18 66 72 76 3 1624 457 63
## -Andre Dawson 496 141 20 65 78 37 11 5628 1575 225
## -Andres Galarraga 321 87 10 39 42 30 2 396 101 12
## -Alfredo Griffin 594 169 4 74 51 35 11 4408 1133 19
## -Al Newman 185 37 1 23 8 21 2 214 42 1
## CRuns CRBI CWalks LeagueN DivisionW PutOuts Assists Errors
## -Alan Ashby 321 414 375 1 1 632 43 10
## -Alvin Davis 224 266 263 0 1 880 82 14
## -Andre Dawson 828 838 354 1 0 200 11 3
## -Andres Galarraga 48 46 33 1 0 805 40 4
## -Alfredo Griffin 501 336 194 0 1 282 421 25
## -Al Newman 30 9 24 1 0 76 127 7
## NewLeagueN
## -Alan Ashby 1
## -Alvin Davis 0
## -Andre Dawson 1
## -Andres Galarraga 1
## -Alfredo Griffin 0
## -Al Newman 0y <- data_ridge$Salary
head(y,6)## [1] 475.0 480.0 500.0 91.5 750.0 70.0Penentuan \(\hat{\mathbf{β}}_{Ridge}\) sangat
bergantung pada \(\lambda\) atau
koefisien smoothing. Package: glmnet memungkinkan untuk
menentukan \(\lambda^{optimal}\) dalam
fungsi glmnet(). Lebih lanjut, fungsi ini juga akan
melakukan standarisasi \((x_{i}-\bar{x})\) secara otomatis dengan
menggunakan opsi standardize = TRUE.
grid_lambda <- 10^seq(10, -2, length=100)
model.ridge <- glmnet(x, y, alpha = 0, lambda = grid_lambda, standardize = TRUE)Dalam hal ini model Regresi Ridge dibentuk dengan menggunakan 100 \(\lambda\) di mana nilai \(\lambda\) telah ditetapkan dari \(\lambda=10^{10}\) dan \(\lambda=10^{-2}\). Dengan demikian terdapat 100 model Regresi Ridge dengan 100 \(\lambda\) yang berbeda.
Sebagai contoh, untuk \(\lambda\) ke-50, maka \(\lambda=11.497,57\) dengan \(\hat{\mathbf{β}}_{Ridge}\) sbb:
model.ridge$lambda[50]## [1] 11497.57coef(model.ridge)[,50]## (Intercept) AtBat Hits HmRun Runs
## 407.356050200 0.036957182 0.138180344 0.524629976 0.230701523
## RBI Walks Years CAtBat CHits
## 0.239841459 0.289618741 1.107702929 0.003131815 0.011653637
## CHmRun CRuns CRBI CWalks LeagueN
## 0.087545670 0.023379882 0.024138320 0.025015421 0.085028114
## DivisionW PutOuts Assists Errors NewLeagueN
## -6.215440973 0.016482577 0.002612988 -0.020502690 0.301433531Untuk memperoleh \(\lambda^{optimal}\) maka digunakan:
set.seed(1)
hasil_cv <- cv.glmnet(x,y,alpha=0)
hasil_cv##
## Call: cv.glmnet(x = x, y = y, alpha = 0)
##
## Measure: Mean-Squared Error
##
## Lambda Index Measure SE Nonzero
## min 25.5 100 113880 15595 19
## 1se 1843.3 54 129347 12158 19dalam hal ini \(\lambda^{optimal}=25,5\).
Selanjutnya akan dilakukan estimasi parameter \(\hat{\mathbf{β}}_{Ridge}\) dengan menggunakan \(\lambda^{optimal}=25,5\).
lambda_cv <- hasil_cv$lambda.min
predict(model.ridge, type = "coefficients", s = lambda_cv)[1:20,]## (Intercept) AtBat Hits HmRun Runs
## 81.29339892 -0.68242406 2.78030368 -1.36065813 1.00250210
## RBI Walks Years CAtBat CHits
## 0.70987266 3.37939096 -9.03223765 -0.00209583 0.13509982
## CHmRun CRuns CRBI CWalks LeagueN
## 0.69188442 0.30168573 0.26073898 -0.27939342 53.20308045
## DivisionW PutOuts Assists Errors NewLeagueN
## -122.77518827 0.26392567 0.17030537 -3.68264211 -18.06289111Untuk melakukan pendugaan \(\hat{\mathbf{β}}_{Lasso}\) maka kembali
akan digunakan fungsi glmnet():
model.lasso <- glmnet(x, y, alpha = 1, lambda = grid_lambda, standardize = TRUE)Kembali akan dilihat \(\hat{\mathbf{β}}_{Lasso}\) pada \(\lambda^{optimal}\) untuk Regresi Lasso:
set.seed(1)
hasil_cv <- cv.glmnet(x,y,alpha=1)
hasil_cv##
## Call: cv.glmnet(x = x, y = y, alpha = 1)
##
## Measure: Mean-Squared Error
##
## Lambda Index Measure SE Nonzero
## min 1.84 54 113582 16663 15
## 1se 69.40 15 129171 13470 6Terlihat pada hasil sebelumnya \(\lambda^{optimal}\) untuk Regresi Lasso adalah 1,84. Dan \(\hat{\mathbf{β}}_{Lasso}\) untuk \(\lambda\) ini adalah:
lambda_cv <- hasil_cv$lambda.min
lasso.coef <- predict(model.lasso , type = "coefficients",
s = lambda_cv)[1:20, ]
lasso.coef## (Intercept) AtBat Hits HmRun Runs
## 1.418021e+02 -1.777554e+00 6.176690e+00 2.708846e-01 -7.050538e-02
## RBI Walks Years CAtBat CHits
## 0.000000e+00 5.139758e+00 -1.009351e+01 -9.197719e-03 0.000000e+00
## CHmRun CRuns CRBI CWalks LeagueN
## 5.360893e-01 7.774982e-01 4.115870e-01 -6.255290e-01 3.523070e+01
## DivisionW PutOuts Assists Errors NewLeagueN
## -1.189054e+02 2.784890e-01 2.242065e-01 -2.444054e+00 -1.154487e+00lasso.coef[lasso.coef!=0]## (Intercept) AtBat Hits HmRun Runs
## 1.418021e+02 -1.777554e+00 6.176690e+00 2.708846e-01 -7.050538e-02
## Walks Years CAtBat CHmRun CRuns
## 5.139758e+00 -1.009351e+01 -9.197719e-03 5.360893e-01 7.774982e-01
## CRBI CWalks LeagueN DivisionW PutOuts
## 4.115870e-01 -6.255290e-01 3.523070e+01 -1.189054e+02 2.784890e-01
## Assists Errors NewLeagueN
## 2.242065e-01 -2.444054e+00 -1.154487e+00Terlihat bahwa Regresi Lasso melakukan adjusment pada koefisien parameter model. Beberapa koefisien terseleksi karena \(\hat{{β}}_{Lasso}\) yang dihasilkan adalah 0. Dengan demikian, model Regresi Lasso dapat juga digunakan sebagai metode seleksi peubah.
Pada contoh ini akan coba dilakukan pendugaan parameter \(\beta\) pada Regresi Logistik dengan Penalized Likelihood menggunakan iterasi Newton-Raphson. Data yang digunakan adalah data simulasi.
Pendugaan parameter \(\beta\) dalam Regresi Logistik ketika terdapat regulasi Ridge dapat dilakukan dengan aturan untuk update \(\beta\) sbb: \(\mathbf{\beta}^{baru}=\mathbf{\beta}^{lama}-[-\mathbf{X}^{T}\mathbf{W}\mathbf{X}-\lambda\mathbf{I}]^{-1}[\mathbf{X}^{T}(\mathbf{y}-\mathbf{p})-\lambda\mathbf{\beta}^{lama}]\)
Metode Newton-Rapshon:
Langkah 1. Menentukan \(p(x;\beta)\)
pi_NR <- function(z) {
return(1 / (1 + exp(-z)))
}Langkah 2 Menghitung \([-\mathbf{X}^{T}\mathbf{W}\mathbf{X}-\lambda\mathbf{I}]\)
hessian_ridge <- function(X, y, beta, lambda) {
z <- X %*% beta
p <- pi_NR(z)
D <- diag(as.vector(p * (1 - p)))
hess <- -t(X) %*% D %*% X - lambda * diag(ncol(X)) # Ridge penalty term
return(hess)
}Langkah 3 Menghitung \([\mathbf{X}^{T}(\mathbf{y}-\mathbf{p})-\lambda\mathbf{\beta}^{lama}]\)
gradient_ridge <- function(X, y, beta, lambda) {
z <- X %*% beta
p <- pi_NR(z)
grad <- t(X) %*% (y - p) - lambda * beta # Ridge penalty term
return(grad)
}Langkah 4 Iterasi
newton_raphson_ridge <- function(X, y, lambda, max_iter = 100, tol = 1e-6) {
# Beta (coefficients) awal
beta <- rep(0, ncol(X))
for (i in 1:max_iter) {
grad <- gradient_ridge(X, y, beta, lambda)
hess <- hessian_ridge(X, y, beta, lambda)
# Update beta
beta_new <- beta - solve(hess) %*% grad
# Apakah konvergen
if (sqrt(sum((beta_new - beta)^2)) < tol) {
break
}
beta <- beta_new
}
return(beta)
}Langkah 5 Uji coba data
set.seed(123)
n <- 100
p <- 2
X <- cbind(1, matrix(rnorm(n * p), n, p))
beta_true <- c(0.5, -1, 1)
y <- rbinom(n, 1, pi_NR(X %*% beta_true))
# Menentukan koefisien lambda
lambda <- 0.1
# Iterasi Newton-Raphson dengan penalty ridge
beta_estimated <- newton_raphson_ridge(X, y, lambda)
# Print the estimated coefficients
print(beta_estimated)## [,1]
## [1,] 0.5809084
## [2,] -0.8657742
## [3,] 1.0574137Metode Newton-Rapshon:
Langkah 1. Menentukan \(p(x;\beta)\)
logistic <- function(x) {
return(1 / (1 + exp(-x)))
}Langkah 2 Menghitung \([-\mathbf{X}^{T}\mathbf{W}\mathbf{X}]\)
hessian <- function(X, beta) {
linear_predictor <- X %*% beta
prob <- logistic(linear_predictor)
W <- diag(as.vector(prob * (1 - prob)))
hess <- -t(X) %*% W %*% X
return(hess)
}Langkah 3 Menghitung \([\mathbf{X}^{T}(\mathbf{y}-\mathbf{p})]\)
gradient <- function(X, y, beta) {
linear_predictor <- X %*% beta
prob <- logistic(linear_predictor)
grad <- t(X) %*% (y - prob)
return(grad)
}Langkah 4 Iterasi
newton_raphson_lasso <- function(X, y, lambda, max_iter = 100, tol = 1e-6) {
n <- ncol(X)
beta <- rep(0, n) # Initialize beta
for (iter in 1:max_iter) {
grad <- gradient(X, y, beta)
hess <- hessian(X, beta)
# Lasso penalty (soft-thresholding)
beta_new <- beta - solve(hess) %*% grad
beta_new <- sign(beta_new) * pmax(abs(beta_new) - lambda, 0)
# Check for convergence
if (max(abs(beta_new - beta)) < tol) {
break
}
beta <- beta_new
}
return(beta)
}Langkah 5 Uji coba data
set.seed(123)
n <- 100
p <- 2
X <- cbind(1, matrix(rnorm(n * p), n, p))
beta_true <- c(0.5, -1, 1)
y <- rbinom(n, 1, pi_NR(X %*% beta_true))
lambda <- 0.1
beta_est <- newton_raphson_lasso(X, y, lambda)
print(beta_est)## [,1]
## [1,] 0.4824945
## [2,] -0.7689640
## [3,] 0.9619525Pada contoh ini akan dilakukan penentuan \(\lambda^{optimal}\) dengan menggunakan metode cross-validation untuk kasus pada Ridge Regression.
# Ridge Regression function
ridge_regression <- function(X, y, lambda) {
# Add intercept column to X
X <- cbind(1, X)
# Ridge Regression formula: beta = (X'X + lambda * I)^(-1) * X'y
I <- diag(ncol(X)) # Identity matrix
beta <- solve(t(X) %*% X + lambda * I) %*% t(X) %*% y
return(beta)
}
# Cross-validation function for Ridge Regression
cross_validate_ridge <- function(X, y, lambda, k = 5) {
n <- nrow(X)
fold_size <- floor(n / k)
mse_values <- numeric(k) # To store mean squared error for each fold
# Randomly shuffle the data
set.seed(123)
shuffled_indices <- sample(n)
X <- X[shuffled_indices, ]
y <- y[shuffled_indices]
# Perform k-fold cross-validation
for (i in 1:k) {
# Split data into training and validation sets
start <- (i - 1) * fold_size + 1
end <- min(i * fold_size, n)
validation_indices <- start:end
train_indices <- setdiff(1:n, validation_indices)
X_train <- X[train_indices, ]
y_train <- y[train_indices]
X_val <- X[validation_indices, ]
y_val <- y[validation_indices]
# Train Ridge Regression model
beta <- ridge_regression(X_train, y_train, lambda)
# Make predictions on the validation set
y_pred <- cbind(1, X_val) %*% beta
# Calculate mean squared error (MSE)
mse_values[i] <- mean((y_val - y_pred)^2)
}
# Return the average MSE across all folds
return(mean(mse_values))
}
# Example using mtcars dataset
data(mtcars)
X <- as.matrix(mtcars[, -1]) # Predictors (all columns except the first)
y <- mtcars[, 1] # Response variable (mpg)
# Define a range of lambda values to test
lambda_values <- 10^seq(-2, 2, length = 100)
# Perform cross-validation for each lambda
cv_results <- sapply(lambda_values, function(lambda) {
cross_validate_ridge(X, y, lambda, k = 5)
})
# Find the optimal lambda (the one with the lowest MSE)
optimal_lambda <- lambda_values[which.min(cv_results)]
print(paste("Optimal lambda:", optimal_lambda))## [1] "Optimal lambda: 4.2292428743895"# Fit the final Ridge Regression model using the optimal lambda
final_beta <- ridge_regression(X, y, optimal_lambda)
print("Final coefficients:")## [1] "Final coefficients:"print(final_beta)## [,1]
## 0.261509093
## cyl 0.353717734
## disp -0.005934733
## hp -0.009650392
## drat 1.448268980
## wt -1.775218368
## qsec 0.935844687
## vs -0.074227010
## am 1.524066688
## gear 1.547935542
## carb -0.748453841# Make predictions using the final model
y_pred <- cbind(1, X) %*% final_beta
# Calculate R-squared to evaluate model performance
ss_total <- sum((y - mean(y))^2)
ss_residual <- sum((y - y_pred)^2)
r_squared <- 1 - (ss_residual / ss_total)
print(paste("R-squared:", r_squared))## [1] "R-squared: 0.850010326521097"