Ridge Regression

# Ridge Regression Example 1: House Prices with Multicollinearity
# Compare OLS (unstable) vs Ridge (stable) predictions

library(glmnet)      # Ridge regression
library(MASS)        # For ridge (alternative)
library(car)         # VIF calculation
library(ggplot2)     # Visualization
library(dplyr)       # Data manipulation
library(caret)       # Cross-validation

set.seed(2025)

# === STEP 1: Simulate House Price Data with Multicollinearity ===
n <- 200
p <- 15

# Create correlated predictor matrix
# Simulate from multivariate normal with high correlations
library(MASS)
cor_matrix <- matrix(0.7, p, p)  # High correlation (0.7) between all predictors
diag(cor_matrix) <- 1

X <- mvrnorm(n, mu=rep(50, p), Sigma=cor_matrix * 100)
colnames(X) <- paste0("X", 1:p)

# True coefficients (sparse: only first 5 matter)
true_beta <- c(10, 8, -6, 5, 4, rep(0, 10))  

# Generate outcome with noise
y <- X %*% true_beta + rnorm(n, 0, sd=20)

data <- data.frame(y = as.vector(y), X)

cat("=== Data Simulation Complete ===")
cat("\nn =", n, ", p =", p)
cat("\nTrue non-zero coefficients: 5 out of", p)
cat("\nCorrelation between predictors: 0.7 (high multicollinearity)\n")

# === STEP 2: Check Multicollinearity (VIF) ===
model_ols_full <- lm(y ~ ., data=data)
vif_values <- vif(model_ols_full)

cat("\n=== Variance Inflation Factors(VIF) ===")
print(round(vif_values, 2))
cat("\nVIF > 10 indicates severe multicollinearity")
cat("\nMean VIF:", round(mean(vif_values), 2))
cat("\nMax VIF:", round(max(vif_values), 2), "\n")

if (max(vif_values) > 10) {
  cat("\n*** SEVERE MULTICOLLINEARITY DETECTED ***")
  cat("\nOLS will be unstable; Ridge regression recommended\n")
}

# === STEP 3: Train-Test Split ===
set.seed(123)
train_idx <- sample(1:n, size=0.7*n)
train_data <- data[train_idx, ]
test_data <- data[-train_idx, ]

cat("\n=== Train-Test Split ===")
cat("\nTraining set:", nrow(train_data), "observations")
cat("\nTest set:", nrow(test_data), "observations\n")

# Prepare matrices for glmnet
X_train <- as.matrix(train_data[, -1])
y_train <- train_data$y
X_test <- as.matrix(test_data[, -1])
y_test <- test_data$y

# === STEP 4: OLS Regression (Baseline) ===
model_ols <- lm(y ~ ., data=train_data)
summary(model_ols)

# OLS predictions on test set
pred_ols_test <- predict(model_ols, newdata=test_data)
mse_ols_test <- mean((y_test - pred_ols_test)^2)
rmse_ols_test <- sqrt(mse_ols_test)
r2_ols_test <- 1 - mse_ols_test / var(y_test)

cat("\n=== OLS Performance(Test Set) ===")
cat("\nMSE:", round(mse_ols_test, 2))
cat("\nRMSE:", round(rmse_ols_test, 2))
cat("\nR²:", round(r2_ols_test, 3), "\n")

# OLS coefficient instability
cat("\n=== OLS Coefficients(Unstable due to collinearity) ===")
print(round(coef(model_ols)[-1], 3))  # Exclude intercept

# === STEP 5: Ridge Regression with CV for λ Selection ===

# Fit ridge with cross-validation to select λ
# alpha=0 means ridge (alpha=1 is lasso, alpha in (0,1) is elastic net)
cv_ridge <- cv.glmnet(X_train, y_train, 
                      alpha=0,           # Ridge (L2 penalty)
                      nfolds=10,         # 10-fold CV
                      standardize=TRUE)  # Auto-standardize predictors

cat("\n=== Ridge Regression: Cross-Validation for λ ===")

# Plot CV curve
plot(cv_ridge, main="Ridge Regression: CV Error vs Log(λ)")
abline(v=log(cv_ridge$lambda.min), col="red", lty=2)
abline(v=log(cv_ridge$lambda.1se), col="blue", lty=2)
legend("topright", 
       legend=c("λ min(lowest CV error)", "λ 1SE (simplest within 1 SE)"),
       col=c("red", "blue"), lty=2)

cat("\nλ that minimizes CV error:", round(cv_ridge$lambda.min, 4))
cat("\nλ within 1 SE(more regularization):", round(cv_ridge$lambda.1se, 4))
cat("\nMin CV MSE:", round(min(cv_ridge$cvm), 2), "\n")

# === STEP 6: Extract Ridge Coefficients at Optimal λ ===

# Use lambda.min for best prediction
ridge_coef_min <- coef(cv_ridge, s="lambda.min")
cat("\n=== Ridge Coefficients(λ = lambda.min) ===")
print(round(as.vector(ridge_coef_min)[-1], 3))  # Exclude intercept

# Use lambda.1se for more regularization (simpler model)
ridge_coef_1se <- coef(cv_ridge, s="lambda.1se")
cat("\n=== Ridge Coefficients(λ = lambda.1se, more shrinkage) ===")
print(round(as.vector(ridge_coef_1se)[-1], 3))

# === STEP 7: Ridge Predictions on Test Set ===

pred_ridge_min <- predict(cv_ridge, newx=X_test, s="lambda.min")
mse_ridge_min <- mean((y_test - pred_ridge_min)^2)
rmse_ridge_min <- sqrt(mse_ridge_min)
r2_ridge_min <- 1 - mse_ridge_min / var(y_test)

cat("\n=== Ridge Performance(Test Set, λ = lambda.min) ===")
cat("\nMSE:", round(mse_ridge_min, 2))
cat("\nRMSE:", round(rmse_ridge_min, 2))
cat("\nR²:", round(r2_ridge_min, 3))

pred_ridge_1se <- predict(cv_ridge, newx=X_test, s="lambda.1se")
mse_ridge_1se <- mean((y_test - pred_ridge_1se)^2)

cat("\n=== Ridge Performance(Test Set, λ = lambda.1se) ===")
cat("\nMSE:", round(mse_ridge_1se, 2), "\n")

# === STEP 8: Compare OLS vs Ridge ===

cat("\n=== Model Comparison(Test Set MSE) ===")
cat("\nOLS MSE:", round(mse_ols_test, 2))
cat("\nRidge MSE(λ.min):", round(mse_ridge_min, 2))
cat("\nRidge MSE(λ.1se):", round(mse_ridge_1se, 2))
cat("\n\nImprovement(OLS → Ridge):", 
    round((mse_ols_test - mse_ridge_min) / mse_ols_test * 100, 1), "%\n")

if (mse_ridge_min < mse_ols_test) {
  cat("\n*** Ridge OUTPERFORMS OLS(lower test MSE) ***")
  cat("\nMulticollinearity successfully handled by regularization\n")
}

# === STEP 9: Coefficient Path Plot (Ridge Trace) ===

# Fit ridge across sequence of λ values
ridge_path <- glmnet(X_train, y_train, alpha=0, standardize=TRUE)

plot(ridge_path, xvar="lambda", label=TRUE,
     main="Ridge Coefficient Paths",
     xlab="Log(λ)", ylab="Standardized Coefficients")
abline(v=log(cv_ridge$lambda.min), col="red", lty=2)
abline(h=0, col="gray", lty=2)
legend("topright", legend="λ.min", col="red", lty=2)

cat("\n=== Coefficient Shrinkage ===")
cat("\nAs λ increases(moving right), all coefficients shrink toward zero")
cat("\nBut NEVER reach exact zero(unlike lasso)")
cat("\nCoefficients most affected by multicollinearity shrink fastest\n")

# === STEP 10: Effective Degrees of Freedom ===

# df(λ) = trace(H_λ) where H_λ is the hat matrix
# Measures effective number of parameters
edf_min <- cv_ridge$glmnet.fit$df[which(cv_ridge$lambda == cv_ridge$lambda.min)]

cat("\n=== Model Complexity ===")
cat("\nOLS degrees of freedom:", p, "(all predictors)")
cat("\nRidge effective df(λ.min):", round(edf_min, 2))
cat("\nReduction:", round(p - edf_min, 2), "\n")
cat("\nRidge effectively uses fewer parameters due to shrinkage\n")

# === STEP 11: Visualizations ===

# Predicted vs Observed (Test Set)
results <- data.frame(
  Observed = y_test,
  OLS = pred_ols_test,
  Ridge = as.vector(pred_ridge_min)
)

ggplot(results, aes(x=Observed)) +
  geom_point(aes(y=OLS, color="OLS"), alpha=0.6, size=2) +
  geom_point(aes(y=Ridge, color="Ridge"), alpha=0.6, size=2) +
  geom_abline(slope=1, intercept=0, linetype="dashed", color="black") +
  scale_color_manual(values=c("OLS"="red", "Ridge"="blue")) +
  labs(title="Predicted vs Observed House Prices(Test Set)",
       x="Observed Price", y="Predicted Price", color="Model") +
  theme_classic() +
  theme(legend.position=c(0.15, 0.85))

# Coefficient comparison: OLS vs Ridge
coef_compare <- data.frame(
  Predictor = paste0("X", 1:p),
  OLS = coef(model_ols)[-1],
  Ridge = as.vector(ridge_coef_min)[-1],
  True = true_beta
)

coef_long <- reshape2::melt(coef_compare, id.vars="Predictor")

ggplot(coef_long, aes(x=Predictor, y=value, fill=variable)) +
  geom_bar(stat="identity", position="dodge", alpha=0.7) +
  scale_fill_manual(values=c("OLS"="red", "Ridge"="blue", "True"="green")) +
  labs(title="Coefficient Comparison: OLS vs Ridge vs True",
       x="Predictor", y="Coefficient Value", fill="Model") +
  theme_classic() +
  theme(axis.text.x = element_text(angle=45, hjust=1))

cat("\n=== Visual Insights ===")
cat("\nRidge coefficients(blue) are shrunk toward zero compared to OLS(red)")
cat("\nRidge is closer to true coefficients(green) due to bias-variance tradeoff")
cat("\nOLS coefficients are unstable(high variance) due to multicollinearity\n")

# === STEP 12: Residual Diagnostics (Ridge) ===

residuals_ridge <- y_test - pred_ridge_min

par(mfrow=c(2,2))

# Residuals vs Fitted
plot(pred_ridge_min, residuals_ridge,
     main="Ridge: Residuals vs Fitted(Test Set)",
     xlab="Fitted Values", ylab="Residuals")
abline(h=0, col="red", lty=2)

# Q-Q Plot
qqnorm(residuals_ridge, main="Ridge: Q-Q Plot(Test Set)")
qqline(residuals_ridge, col="red")

# Residuals histogram
hist(residuals_ridge, breaks=15, main="Ridge: Residual Distribution",
     xlab="Residuals", col="lightblue")

# Scale-location
plot(pred_ridge_min, sqrt(abs(residuals_ridge)),
     main="Ridge: Scale-Location(Test Set)",
     xlab="Fitted Values", ylab="√|Residuals|")
abline(h=mean(sqrt(abs(residuals_ridge))), col="red", lty=2)

par(mfrow=c(1,1))

# === STEP 13: Cross-Validation Stability ===

# Repeat CV with different seed to check λ stability
set.seed(456)
cv_ridge2 <- cv.glmnet(X_train, y_train, alpha=0, nfolds=10, standardize=TRUE)

cat("\n=== CV Stability Check ===")
cat("\nλ.min(seed=2025):", round(cv_ridge$lambda.min, 4))
cat("\nλ.min(seed=456):", round(cv_ridge2$lambda.min, 4))
cat("\nDifference:", round(abs(cv_ridge$lambda.min - cv_ridge2$lambda.min), 4))
cat("\nStable λ selection indicates robust model\n")

# === APA-Style Reporting ===

cat("\n" , "="*70, "\n")
cat("=== APA-STYLE REPORT ===")
cat("\n", "="*70, "\n")

cat("\nRidge regression(L2 regularization) was used to predict house prices\n")
cat("from", p, "predictors(N =", n, ", 70/30 train/test split). Predictors\n")
cat("exhibited severe multicollinearity(mean VIF =", round(mean(vif_values), 1), ",\n")
cat("max VIF =", round(max(vif_values), 1), "), causing unstable OLS estimates.\n")
cat("\n")
cat("The ridge penalty parameter λ was selected via 10-fold cross-validation\n")
cat("on the training set(λ =", round(cv_ridge$lambda.min, 4), "). Ridge regression\n")
cat("substantially outperformed OLS on the held-out test set: Ridge MSE =",
    round(mse_ridge_min, 2), ", OLS MSE =", round(mse_ols_test, 2),
    "(improvement:",
    round((mse_ols_test - mse_ridge_min)/mse_ols_test*100, 1), "%).\n")
cat("\n")
cat("Ridge coefficients were shrunk toward zero(regularization) but remained\n")
cat("non-zero for all predictors(L2 penalty does not perform variable selection).\n")
cat("The effective degrees of freedom decreased from", p, "(OLS) to",
    round(edf_min, 1), "(Ridge), indicating reduced model complexity.\n")
cat("\n")
cat("Residual diagnostics on the test set showed approximately normal errors\n")
cat("with no systematic patterns, validating model assumptions. Ridge regression\n")
cat("successfully stabilized predictions in the presence of multicollinearity,\n")
cat("demonstrating the bias-variance tradeoff: accepting small bias(shrinkage)\n")
cat("in exchange for substantial variance reduction and improved prediction.\n")

cat("\n", "="*70, "\n")