PENALIZED REGRESSION: BRIDGE, RIDGE, dan LASSO (Body Fat)
Kelompok 5 KAPITA SELEKTA UNIB
2026-02-24
Nama anggota kelompok 4:
1. M. Syarlan
(F1F023005)
2. A’inun Fauziyah (F1F023017)
3. Zahra Hana Andrea
(F1F023027)
4. Chesa Nabilla Azzahra (F1F023040)
1. PERSIAPAN
library(readxl)## Warning: package 'readxl' was built under R version 4.3.3library(car)## Warning: package 'car' was built under R version 4.3.3## Loading required package: carDatalibrary(lmtest)## Warning: package 'lmtest' was built under R version 4.3.3## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(caret)## Warning: package 'caret' was built under R version 4.3.3## Loading required package: ggplot2## Loading required package: latticelibrary(glmnet)## Warning: package 'glmnet' was built under R version 4.3.3## Loading required package: Matrix## Loaded glmnet 4.1-8library(tidyr)##
## Attaching package: 'tidyr'## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpacklibrary(dplyr)## Warning: package 'dplyr' was built under R version 4.3.3##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union2. LOAD DATA
BodyFat_Extended <- read_excel(
"D:/SEMESETER 6/KAPITA SELEKTA/BodyFat - Extended.xlsx"
)
head(BodyFat_Extended)## # A tibble: 6 × 7
## BodyFat Age Weight Height Neck Abdomen Hip
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 12.3 23 70.0 1.72 36.2 85.2 94.5
## 2 6.1 22 78.6 1.84 38.5 83 98.7
## 3 25.3 22 69.8 1.68 34 87.9 99.2
## 4 10.4 26 83.8 1.84 37.4 86.4 101.
## 5 28.7 24 83.6 1.81 34.4 100 102.
## 6 20.9 24 95.4 1.9 39 94.4 108.3. DEFINISI VARIABEL
X <- as.matrix(BodyFat_Extended[, c(
"Age", "Weight", "Height", "Neck", "Abdomen", "Hip"
)])
Y <- BodyFat_Extended$BodyFat
n <- length(Y)
# Standardisasi
X_std <- scale(X)
X_design <- cbind(Intercept = 1, X_std)4. MODEL OLS
model_ols <- lm(
BodyFat ~ Age + Weight + Height + Neck + Abdomen + Hip,
data = BodyFat_Extended
)
yhat_ols <- predict(model_ols)
coef_ols <- coef(model_ols)
max(vif(model_ols))## [1] 27.42996vif(model_ols)## Age Weight Height Neck Abdomen Hip
## 3.145333 27.429958 2.041993 9.018668 15.306452 5.7576875. OBJECTIVE FUNCTION BRIDGE
bridge_objective <- function(beta, X, y, lambda, gamma){
residual <- y - X %*% beta
SSE <- sum(residual^2)
penalty <- lambda * sum(abs(beta[-1])^gamma)
return(SSE + penalty)
}6. CROSS VALIDATION FUNCTION
bridge_cv <- function(X, y, lambda, gamma, K = 5){
folds <- createFolds(y, k = K, list = TRUE)
cv_error <- rep(0, K)
for(i in 1:K){
test_idx <- folds[[i]]
train_idx <- setdiff(1:length(y), test_idx)
fit <- optim(
par = rep(0, ncol(X)),
fn = bridge_objective,
X = X[train_idx, ],
y = y[train_idx],
lambda = lambda,
gamma = gamma,
method = "BFGS"
)
beta_hat <- fit$par
y_pred <- X[test_idx, ] %*% beta_hat
cv_error[i] <- mean((y[test_idx] - y_pred)^2)
}
return(mean(cv_error))
}7. FUNGSI HITUNG VIF PASCA PENALIZED
hitung_vif_post <- function(beta, X_std, y){
vars <- colnames(X_std)[abs(beta[-1]) > 1e-6]
if(length(vars) < 2){
return(NA)
}
df <- as.data.frame(X_std)
colnames(df) <- colnames(X_std)
df$Y <- y
model_subset <- lm(
as.formula(paste("Y ~", paste(vars, collapse="+"))),
data=df
)
return(max(vif(model_subset)))
}8. ESTIMASI SEMUA MODEL
gamma_grid <- c(0.5, 1, 1.5, 2)
results <- data.frame()
coef_list <- list()
pred_list <- list()
############################
# OLS
############################
SSE_ols <- sum((Y - yhat_ols)^2)
df_ols <- length(coef_ols)
results <- rbind(results,
data.frame(
Model="OLS",
Lambda=0,
MSE=mean((Y-yhat_ols)^2),
RMSE=sqrt(mean((Y-yhat_ols)^2)),
R2=summary(model_ols)$r.squared,
R2_adj=summary(model_ols)$adj.r.squared,
AIC=AIC(model_ols),
BIC=BIC(model_ols),
VIF=max(vif(model_ols))
))
coef_list[["OLS"]] <- coef_ols
pred_list[["OLS"]] <- yhat_ols
############################
# BRIDGE
############################
for(g in gamma_grid){
cat("Proses Bridge gamma =", g, "\n")
lambda_grid <- seq(0.01,100,length=25)
cv_vals <- sapply(lambda_grid,
function(lam) bridge_cv(X_design,Y,lam,g))
best_lambda <- lambda_grid[which.min(cv_vals)]
fit <- optim(
par = rep(0,ncol(X_design)),
fn = bridge_objective,
X = X_design,
y = Y,
lambda = best_lambda,
gamma = g,
method="BFGS"
)
beta_hat <- fit$par
yhat <- X_design %*% beta_hat
SSE <- sum((Y-yhat)^2)
df_eff <- sum(abs(beta_hat[-1])>1e-6) + 1
AIC_pen <- n*log(SSE/n) + 2*df_eff
BIC_pen <- n*log(SSE/n) + log(n)*df_eff
r2 <- 1 - SSE/sum((Y-mean(Y))^2)
r2_adj <- 1 - ((1-r2)*(n-1)/(n-df_eff))
vif_post <- hitung_vif_post(beta_hat, X_std, Y)
results <- rbind(results,
data.frame(
Model=paste0("Bridge γ=",g),
Lambda=round(best_lambda,4),
MSE=mean((Y-yhat)^2),
RMSE=sqrt(mean((Y-yhat)^2)),
R2=r2,
R2_adj=r2_adj,
AIC=AIC_pen,
BIC=BIC_pen,
VIF=vif_post
))
coef_list[[paste0("Bridge_",g)]] <- beta_hat
pred_list[[paste0("Bridge γ=",g)]] <- yhat
}## Proses Bridge gamma = 0.5
## Proses Bridge gamma = 1
## Proses Bridge gamma = 1.5
## Proses Bridge gamma = 2############################
# RIDGE
############################
cv_ridge <- cv.glmnet(X_std,Y,alpha=0)
lambda_ridge <- cv_ridge$lambda.min
ridge_pred <- predict(cv_ridge,X_std,s="lambda.min")
coef_ridge <- as.vector(coef(cv_ridge,s="lambda.min"))
SSE_ridge <- sum((Y-ridge_pred)^2)
df_ridge <- sum(abs(coef_ridge[-1])>1e-6)+1
AIC_ridge <- n*log(SSE_ridge/n)+2*df_ridge
BIC_ridge <- n*log(SSE_ridge/n)+log(n)*df_ridge
r2_ridge <- 1-SSE_ridge/sum((Y-mean(Y))^2)
r2adj_ridge <- 1-((1-r2_ridge)*(n-1)/(n-df_ridge))
results <- rbind(results,
data.frame(
Model="Ridge",
Lambda=round(lambda_ridge,4),
MSE=mean((Y-ridge_pred)^2),
RMSE=sqrt(mean((Y-ridge_pred)^2)),
R2=r2_ridge,
R2_adj=r2adj_ridge,
AIC=AIC_ridge,
BIC=BIC_ridge,
VIF=hitung_vif_post(coef_ridge,X_std,Y)
))
coef_list[["Ridge"]] <- coef_ridge
pred_list[["Ridge"]] <- ridge_pred
############################
# LASSO
############################
cv_lasso <- cv.glmnet(X_std,Y,alpha=1)
lambda_lasso <- cv_lasso$lambda.min
lasso_pred <- predict(cv_lasso,X_std,s="lambda.min")
coef_lasso <- as.vector(coef(cv_lasso,s="lambda.min"))
SSE_lasso <- sum((Y-lasso_pred)^2)
df_lasso <- sum(abs(coef_lasso[-1])>1e-6)+1
AIC_lasso <- n*log(SSE_lasso/n)+2*df_lasso
BIC_lasso <- n*log(SSE_lasso/n)+log(n)*df_lasso
r2_lasso <- 1-SSE_lasso/sum((Y-mean(Y))^2)
r2adj_lasso <- 1-((1-r2_lasso)*(n-1)/(n-df_lasso))
results <- rbind(results,
data.frame(
Model="LASSO",
Lambda=round(lambda_lasso,4),
MSE=mean((Y-lasso_pred)^2),
RMSE=sqrt(mean((Y-lasso_pred)^2)),
R2=r2_lasso,
R2_adj=r2adj_lasso,
AIC=AIC_lasso,
BIC=BIC_lasso,
VIF=hitung_vif_post(coef_lasso,X_std,Y)
))
coef_list[["LASSO"]] <- coef_lasso
pred_list[["LASSO"]] <- lasso_pred9. TABEL FINAL
results[sapply(results,is.numeric)] <-
round(results[sapply(results,is.numeric)],4)
print(results)## Model Lambda MSE RMSE R2 R2_adj AIC BIC VIF
## 1 OLS 0.0000 24.9775 4.9978 0.5552 0.5490 2656.352 2688.973 27.4300
## 2 Bridge γ=0.5 8.3425 24.9783 4.9978 0.5552 0.5490 1417.052 1445.595 27.4300
## 3 Bridge γ=1 100.0000 25.5709 5.0568 0.5446 0.5383 1427.275 1455.818 27.4300
## 4 Bridge γ=1.5 16.6750 25.1684 5.0168 0.5518 0.5455 1420.356 1448.900 27.4300
## 5 Bridge γ=2 25.0075 26.1238 5.1111 0.5348 0.5283 1436.601 1465.144 27.4300
## 6 Ridge 0.4410 26.1511 5.1138 0.5343 0.5278 1437.056 1465.600 27.4300
## 7 LASSO 0.0886 25.4602 5.0458 0.5466 0.5424 1421.382 1441.771 7.4658MENAMPILKAN NILAI KOEFISIEN
var_names <- c("Intercept", colnames(X_std))
coef_table <- do.call(cbind, coef_list)
coef_table <- as.data.frame(coef_table)
coef_table$Variable <- var_names
coef_table <- coef_table[, c(ncol(coef_table), 1:(ncol(coef_table)-1))]
coef_table[,-1] <- round(coef_table[,-1], 4)
print(coef_table)## Variable OLS Bridge_0.5 Bridge_1 Bridge_1.5 Bridge_2 Ridge
## (Intercept) Intercept -5.1434 20.2498 20.2499 20.2498 20.2498 20.2498
## Age Age -0.0177 -0.2372 0.0001 -0.0745 0.0238 0.0206
## Weight Weight -0.1892 -2.8177 -0.0004 -1.2717 -0.0654 -0.0541
## Height Height -17.0617 -1.6891 -2.0581 -1.9260 -2.1481 -2.1506
## Neck Neck -1.1927 -4.5940 -4.1943 -4.3477 -3.2205 -3.1916
## Abdomen Abdomen 0.4955 7.1639 4.6369 5.7214 3.7199 3.6893
## Hip Hip 0.7124 4.8949 4.0779 4.4782 4.1160 4.1079
## LASSO
## (Intercept) 20.2498
## Age 0.0000
## Weight 0.0000
## Height -2.0605
## Neck -4.4671
## Abdomen 4.9196
## Hip 4.0633VIF DETAIL PER VARIABEL
hitung_vif_per_variabel <- function(beta, X_std, y, nama_model){
vars <- colnames(X_std)[abs(beta[-1]) > 1e-6]
if(length(vars) < 2){
cat("\nModel", nama_model,
"hanya memiliki <=1 variabel aktif. VIF tidak dapat dihitung.\n")
return(NULL)
}
df <- as.data.frame(X_std)
colnames(df) <- colnames(X_std)
df$Y <- y
model_subset <- lm(
as.formula(paste("Y ~", paste(vars, collapse="+"))),
data=df
)
vif_vals <- vif(model_subset)
hasil <- data.frame(
Model = nama_model,
Variabel = names(vif_vals),
VIF = round(as.numeric(vif_vals),4)
)
return(hasil)
}
############################################################
# MEMBUAT DATA VIF DETAIL SEMUA MODEL
############################################################
vif_all_detail <- data.frame()
# OLS
vif_all_detail <- rbind(
vif_all_detail,
data.frame(
Model = "OLS",
Variabel = names(vif(model_ols)),
VIF = as.numeric(vif(model_ols))
)
)
# Bridge
for(g in gamma_grid){
vif_temp <- hitung_vif_per_variabel(
coef_list[[paste0("Bridge_",g)]],
X_std,
Y,
paste0("Bridge γ=",g)
)
if(!is.null(vif_temp)){
vif_all_detail <- rbind(vif_all_detail, vif_temp)
}
}
# Ridge
vif_ridge <- hitung_vif_per_variabel(
coef_list[["Ridge"]],
X_std,
Y,
"Ridge"
)
if(!is.null(vif_ridge)){
vif_all_detail <- rbind(vif_all_detail, vif_ridge)
}
# LASSO
vif_lasso <- hitung_vif_per_variabel(
coef_list[["LASSO"]],
X_std,
Y,
"LASSO"
)
if(!is.null(vif_lasso)){
vif_all_detail <- rbind(vif_all_detail, vif_lasso)
}
vif_all_detail$VIF <- round(vif_all_detail$VIF,4)
############################################################
# PIVOT TABEL VIF
############################################################
vif_pivot <- vif_all_detail %>%
pivot_wider(
names_from = Model,
values_from = VIF
)
print(vif_pivot)## # A tibble: 6 × 8
## Variabel OLS `Bridge γ=0.5` `Bridge γ=1` `Bridge γ=1.5` `Bridge γ=2` Ridge
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Age 3.15 3.15 3.15 3.15 3.15 3.15
## 2 Weight 27.4 27.4 27.4 27.4 27.4 27.4
## 3 Height 2.04 2.04 2.04 2.04 2.04 2.04
## 4 Neck 9.02 9.02 9.02 9.02 9.02 9.02
## 5 Abdomen 15.3 15.3 15.3 15.3 15.3 15.3
## 6 Hip 5.76 5.76 5.76 5.76 5.76 5.76
## # ℹ 1 more variable: LASSO <dbl>10. PLOT AKTUAL VS PREDIKSI
par(mfrow=c(3,3))
for(name in names(pred_list)){
plot(Y, pred_list[[name]],
main=name,
pch=16)
abline(0,1,col="red",lwd=2)
}
par(mfrow=c(1,1))