PENALIZED REGRESSION: BRIDGE, RIDGE, dan LASSO (Life Expectancy)
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
data <- read_excel(
"D:/SEMESETER 6/KAPITA SELEKTA/life_expectancy_bridge.xlsx"
)
head(data)## # A tibble: 6 × 8
## life_expectancy gdp schooling income_composition_of_resou…¹ adult_mortality
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 65.5 41.9 10.3 0.511 22
## 2 63.1 1158. 8.5 0.513 25
## 3 67.5 1625. 11.2 0.565 19
## 4 75 9644. 13.1 0.787 123
## 5 69.7 919. 11.3 0.625 161
## 6 74 419. 12.3 0.692 146
## # ℹ abbreviated name: ¹income_composition_of_resources
## # ℹ 3 more variables: infant_deaths <dbl>, bmi <dbl>,
## # percentage_expenditure <dbl>3. DEFINISI VARIABEL
X <- as.matrix(data[, c(
"gdp", "schooling" , "income_composition_of_resources",
"adult_mortality", "infant_deaths", "bmi", "percentage_expenditure"
)])
Y <- data$life_expectancy
n <- length(Y)
# Standardisasi
X_std <- scale(X)
X_design <- cbind(Intercept = 1, X_std)4. MODEL OLS
model_ols <- lm(life_expectancy ~ gdp + schooling + income_composition_of_resources +
adult_mortality + infant_deaths + bmi + percentage_expenditure,
data = data)
yhat_ols <- predict(model_ols)
coef_ols <- coef(model_ols)
max(vif(model_ols))## [1] 5.662507vif(model_ols)## gdp schooling
## 5.662507 3.169110
## income_composition_of_resources adult_mortality
## 3.013710 1.342304
## infant_deaths bmi
## 1.070506 1.588180
## percentage_expenditure
## 5.2729705. 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 22.5773 4.7516 0.7583 0.7576 14363.163 14415.245 5.6625
## 2 Bridge γ=0.5 91.6675 22.5780 4.7516 0.7583 0.7576 7524.795 7571.091 5.6625
## 3 Bridge γ=1 95.8338 22.5783 4.7517 0.7583 0.7576 7524.823 7571.119 5.6625
## 4 Bridge γ=1.5 87.5012 22.5808 4.7519 0.7583 0.7576 7525.098 7571.394 5.6625
## 5 Bridge γ=2 58.3375 22.5883 4.7527 0.7582 0.7575 7525.899 7572.194 5.6625
## 6 Ridge 0.7293 22.6711 4.7614 0.7573 0.7566 7534.703 7580.999 5.6625
## 7 LASSO 0.0363 22.5805 4.7519 0.7583 0.7576 7525.066 7571.362 5.6625MENAMPILKAN 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
## (Intercept) Intercept 54.8727
## gdp gdp 0.0000
## schooling schooling 0.9036
## income_composition_of_resources income_composition_of_resources 9.4620
## adult_mortality adult_mortality -0.0296
## infant_deaths infant_deaths -0.0022
## bmi bmi 0.0557
## percentage_expenditure percentage_expenditure 0.0001
## Bridge_0.5 Bridge_1 Bridge_1.5 Bridge_2 Ridge
## (Intercept) 69.5318 69.5318 69.5318 69.5318 69.5318
## gdp 0.4380 0.4235 0.4275 0.4432 0.4828
## schooling 3.0304 3.0208 2.9813 2.9298 2.7804
## income_composition_of_resources 2.0128 2.0117 2.0218 2.0466 2.0797
## adult_mortality -3.7331 -3.7229 -3.6924 -3.6511 -3.4914
## infant_deaths -0.2588 -0.2607 -0.2703 -0.2792 -0.2837
## bmi 1.1026 1.1009 1.1133 1.1382 1.1878
## percentage_expenditure 0.2574 0.2764 0.2935 0.3012 0.3230
## LASSO
## (Intercept) 69.5318
## gdp 0.4265
## schooling 3.0232
## income_composition_of_resources 2.0079
## adult_mortality -3.7135
## infant_deaths -0.2475
## bmi 1.0947
## percentage_expenditure 0.2630VIF 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: 7 × 8
## Variabel OLS `Bridge γ=0.5` `Bridge γ=1` `Bridge γ=1.5` `Bridge γ=2` Ridge
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 gdp 5.66 5.66 5.66 5.66 5.66 5.66
## 2 schooling 3.17 3.17 3.17 3.17 3.17 3.17
## 3 income_co… 3.01 3.01 3.01 3.01 3.01 3.01
## 4 adult_mor… 1.34 1.34 1.34 1.34 1.34 1.34
## 5 infant_de… 1.07 1.07 1.07 1.07 1.07 1.07
## 6 bmi 1.59 1.59 1.59 1.59 1.59 1.59
## 7 percentag… 5.27 5.27 5.27 5.27 5.27 5.27
## # ℹ 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))