---
title: "Fused Ridge Regression on High-Dimensional (Wide) Data"
author: "ANTsR Team"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Fused Ridge Regression on High-Dimensional (Wide) Data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)
library(ANTsR)
library(glmnet)
library(ggplot2)
library(reshape2)
```

## Introduction

In neuroimaging, we often encounter "wide" datasets where the number of predictors $M$ (voxels, vertices, or brain signatures) exceeds the number of observations $N$. Furthermore, the response variable is often a continuous clinical score that we wish to model across various severity thresholds.

This vignette demonstrates how to use `fusedRidge` to handle high-dimensional data by:
1.  Smoothing coefficients across thresholds using a Graph Laplacian (Fusion).
2.  Employing a two-stage "Sparse Refit" to identify a robust subset of features.
3.  Avoiding "Data Leakage" using **Nested Cross-Validation**.

## The Problem: Data Leakage in Two-Stage Models

A common pitfall in two-stage modeling (Select -> Refit) is selecting features using the *entire* dataset and then reporting the cross-validation performance of the second stage. This is **Data Leakage**. Because the features were chosen with knowledge of the test samples, the resulting AUC will be optimistically biased.

**The Best Solution:** Perform the entire two-stage procedure (including feature selection) inside an outer cross-validation loop.

## Realistic Synthetic Experiment

### 1. Data Generation

We simulate $N=100$ subjects and $M=500$ features. Only 20 features have a true signal. The signal varies smoothly across 5 binarization thresholds.

```{r data_gen}
set.seed(42)
N <- 100
M <- 500
J <- 5 # Number of thresholds

# Generate features with some correlation
X <- matrix(rnorm(N * M), N, M)
# Add some structure: features 1-50 are correlated
X[, 1:50] <- X[, 1:50] + 0.5 * rnorm(N)

# Define true coefficients: Sparse and Smooth across thresholds
# Features 10-29 are active
beta_true <- matrix(0, M, J)
for (j in 1:J) {
  # Effect grows and then shrinks across thresholds
  effect_size <- sin(seq(0.5, 2.5, length.out = J))[j]
  beta_true[10:29, j] <- effect_size * seq(1, 0.5, length.out = 20)
}

# Generate latent continuous outcome
y_cont <- X %*% beta_true[, floor(J/2)] + rnorm(N, sd = 0.5)
thresholds <- quantile(y_cont, probs = seq(0.2, 0.8, length.out = J))

# Helper to calculate AUC
get_auc <- function(prob, actual) {
  if (length(unique(actual)) < 2) return(NA)
  p_obj <- pROC::roc(actual, prob, quiet = TRUE)
  return(as.numeric(pROC::auc(p_obj)))
}
```

### 2. Nested Cross-Validation (The "No Leakage" Approach)

We compare three models:
1.  **Standard Ridge:** No fusion ($\lambda_2 = 0$), no sparsity.
2.  **Fused Ridge (Dense):** Fusion ($\lambda_2 = 0.5$), no sparsity.
3.  **Fused Ridge (Sparse):** Fusion + Two-stage `topK=20` refit.

```{r nested_cv}
n_folds <- 5
folds <- sample(rep(1:n_folds, length.out = N))
results_list <- list()

for (f in 1:n_folds) {
  train_idx <- which(folds != f)
  test_idx <- which(folds == f)
  
  X_train <- X[train_idx, ]
  y_train <- y_cont[train_idx]
  X_test <- X[test_idx, ]
  y_test_cont <- y_cont[test_idx]
  
  # binarize test labels
  y_test_mat <- matrix(0, length(test_idx), J)
  for (j in 1:J) y_test_mat[,j] <- ifelse(y_test_cont >= thresholds[j], 1, 0)
  
  # 1. Standard Ridge (Independent per threshold)
  fit_ridge <- fusedRidge(X_train, y_train, thresholds, lambda2 = 0)
  pred_ridge <- predict(fit_ridge, X_test, type = "response")
  
  # 2. Fused Ridge (Dense)
  fit_fused <- fusedRidge(X_train, y_train, thresholds, lambda1 = 0.5, lambda2 = 0.5)
  pred_fused <- predict(fit_fused, X_test, type = "response")
  
# 3. Fused Ridge (Sparse/Two-Stage) - Compare different topK values
  for (tk in c(5, 10, 25, 50)) {
    fit_sparse <- fusedRidge(X_train, y_train, thresholds, lambda1 = 0.5, lambda2 = 0.5, topK = tk)
    pred_sparse <- predict(fit_sparse, X_test, type = "response")
    
    auc_tk <- sapply(1:J, function(j) get_auc(pred_sparse[,j], y_test_mat[,j]))
    results_list[[length(results_list) + 1]] <- data.frame(
      Threshold = 1:J,
      Method = paste0("Fused-Sparse (topK=", tk, ")"),
      AUC = auc_tk,
      Fold = f
    )
  }
  
# 4. Independent Lasso (Sparsity without Fusion)
  pred_lasso_mat <- matrix(0, length(test_idx), J)
  for (j in 1:J) {
    y_train_j <- ifelse(y_train >= thresholds[j], 1, 0)
    # Check if we have both classes
    if (length(unique(y_train_j)) > 1) {
      cv_lasso <- glmnet::cv.glmnet(X_train, y_train_j, family = "binomial", alpha = 1)
      pred_lasso_mat[, j] <- as.vector(predict(cv_lasso, X_test, s = "lambda.min", type = "response"))
    } else {
      pred_lasso_mat[, j] <- mean(y_train_j)
    }
  }
  
  # Calculate AUCs for all methods
  results_list[[length(results_list) + 1]] <- data.frame(
    Threshold = 1:J,
    Method = "Ridge (Independent)",
    AUC = sapply(1:J, function(j) get_auc(pred_ridge[,j], y_test_mat[,j])),
    Fold = f
  )
  results_list[[length(results_list) + 1]] <- data.frame(
    Threshold = 1:J,
    Method = "Lasso (Independent)",
    AUC = sapply(1:J, function(j) get_auc(pred_lasso_mat[,j], y_test_mat[,j])),
    Fold = f
  )
  results_list[[length(results_list) + 1]] <- data.frame(
    Threshold = 1:J,
    Method = "Fused-Dense",
    AUC = sapply(1:J, function(j) get_auc(pred_fused[,j], y_test_mat[,j])),
    Fold = f
  )
}

results_df <- do.call(rbind, results_list)
summary_df <- aggregate(AUC ~ Threshold + Method, data = results_df, FUN = mean)
```

## Results and Visualization

### 1. Predictive Performance (AUC)

We plot the mean AUC across thresholds. This allows us to see how **Fusion** (sharing information across thresholds) and **Sparsity** (feature selection) interact. Note how **Lasso (Independent)** may select features, but because it treats each threshold in isolation, it often suffers from higher variance and lower AUC compared to the **Fused** models.

```{r plot_auc}
# Order methods for better plotting
summary_df$Method <- factor(summary_df$Method, levels = c(
  "Ridge (Independent)", "Lasso (Independent)", "Fused-Dense", 
  "Fused-Sparse (topK=5)", "Fused-Sparse (topK=10)", 
  "Fused-Sparse (topK=25)", "Fused-Sparse (topK=50)"
))

ggplot(summary_df, aes(x = Threshold, y = AUC, color = Method, group = Method)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  theme_minimal() +
  scale_color_viridis_d(option = "plasma") +
  labs(title = "Unbiased AUC Comparison (Nested CV)",
       subtitle = "Fusion vs. Independent Selection vs. Ridge",
       y = "Mean Out-of-Fold AUC") +
  ylim(0.5, 1)
```

### 2. Coefficient Trajectories

To avoid the "spaghetti" effect of overlapping lines, we use a heatmap to visualize the coefficients of the selected features across thresholds. This clearly shows both the **sparsity** (which features were selected) and the **smoothness** (how weights change across thresholds).

```{r plot_coefs_heatmap}
# Fit on full data for visualization
final_fit <- fusedRidge(X, y_cont, thresholds, lambda1 = 0.5, lambda2 = 0.5, topK = 25)

# Extract and melt the coefficients for the selected features
selected_idx <- final_fit$selected_features
coefs_selected <- final_fit$coefs_full[selected_idx, ]
colnames(coefs_selected) <- paste0("Threshold ", 1:J)
rownames(coefs_selected) <- paste0("Feature ", selected_idx)

coefs_melted <- melt(coefs_selected)
colnames(coefs_melted) <- c("Feature", "Threshold", "Weight")

# 1. Heatmap Visualization
ggplot(coefs_melted, aes(x = Threshold, y = Feature, fill = Weight)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0) +
  theme_minimal() +
  labs(title = "Coefficient Heatmap (Top 25 Features)",
       subtitle = "Smooth transitions across thresholds visible horizontally",
       x = "", y = "") +
  theme(axis.text.y = element_text(size = 6))
```

For a more detailed look at the individual trajectories of the top 9 most influential features, we use a faceted line plot.

```{r plot_coefs_faceted}
# Identify top 9 features by average absolute weight
avg_weight <- rowMeans(abs(coefs_selected))
top9_idx <- order(avg_weight, decreasing = TRUE)[1:9]
top9_coefs <- coefs_selected[top9_idx, ]

top9_melted <- melt(top9_coefs)
colnames(top9_melted) <- c("Feature", "Threshold", "Weight")

ggplot(top9_melted, aes(x = Threshold, y = Weight, group = Feature)) +
  geom_line(color = "steelblue", size = 1) +
  geom_point(color = "steelblue", size = 2) +
  facet_wrap(~Feature, scales = "free_y") +
  theme_bw() +
  labs(title = "Individual Trajectories (Top 9 Features)",
       subtitle = "Faceted view eliminates line crossing and 'spaghetti' mess",
       x = "Binarization Threshold",
       y = "Coefficient Weight") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
```

## Anatomically-Grounded Example: Predicting Cognitive Decline

In this section, we use the anatomical labels from the **AAL (Automated Anatomical Labeling)** atlas to simulate a more realistic neuroimaging scenario. We assume that a "Cognitive Decline" score is driven by specific brain circuits (e.g., Frontal and Hippocampal regions) and that this signal evolves across the disease spectrum.

### 1. Data Preparation with Anatomical Context

We generate $N=120$ subjects with $M=116$ brain regions.

```{r aal_setup}
data(aal, package = "ANTsR")
region_names <- aal$label_name
M_aal <- length(region_names)

set.seed(88)
N_aal <- 120
X_aal <- matrix(rnorm(N_aal * M_aal), N_aal, M_aal)
colnames(X_aal) <- region_names

# Define signal in specific circuits:
# 1. Hippocampus (Memory)
hippo_idx <- grep("Hippocampus", region_names)
# 2. Frontal Sup (Executive)
frontal_idx <- grep("Frontal_Sup", region_names)
# 3. Temporal Pole
temporal_idx <- grep("Temporal_Pole", region_names)

active_aal <- c(hippo_idx, frontal_idx, temporal_idx)

# Simulate smooth effect across 4 thresholds
J_aal <- 4
beta_aal <- matrix(0, M_aal, J_aal)
for (j in 1:J_aal) {
  # Signal intensifies at higher severity thresholds
  strength <- j * 0.4
  beta_aal[hippo_idx, j] <- strength
  beta_aal[frontal_idx, j] <- strength * 0.8
  beta_aal[temporal_idx, j] <- strength * 0.6
}

# Generate outcome
y_aal_cont <- X_aal %*% beta_aal[, 2] + rnorm(N_aal, sd = 0.3)
thresholds_aal <- quantile(y_aal_cont, probs = c(0.2, 0.4, 0.6, 0.8))
```

### 2. The Power of Sparse Selection

We fit the `fusedRidge` model with `topK = 15`. This forces the model to select a small subset of regions that consistently contribute across the thresholds.

```{r aal_fit}
# Fit sparse model
fit_aal <- fusedRidge(X_aal, y_aal_cont, thresholds_aal, 
                     lambda1 = 0.5, lambda2 = 1.0, topK = 15)

# Which regions were recovered?
recovered_regions <- region_names[fit_aal$selected_features]
print("Top 15 Recovered Regions:")
print(recovered_regions)

# Check if we found our simulated signal
overlap <- intersect(fit_aal$selected_features, active_aal)
cat("\nSuccessfully recovered", length(overlap), "out of", length(active_aal), "true signal regions.\n")
```

### 3. Evaluating Actual Predictions

To see how the sparse model generalizes, we perform a simple train/test split (80/20) and examine the **Actual vs. Predicted** results for the AAL regions.

```{r aal_eval}
# Train/Test Split
train_idx_aal <- sample(1:N_aal, size = 0.8 * N_aal)
test_idx_aal <- setdiff(1:N_aal, train_idx_aal)

# Fit sparse model on training data
fit_aal_train <- fusedRidge(X_aal[train_idx_aal,], y_aal_cont[train_idx_aal], 
                           thresholds_aal, lambda1 = 0.5, lambda2 = 1.0, topK = 15)

# Get actual binarized labels for test set
y_test_mat_aal <- matrix(0, length(test_idx_aal), J_aal)
for (j in 1:J_aal) y_test_mat_aal[,j] <- ifelse(y_aal_cont[test_idx_aal] >= thresholds_aal[j], 1, 0)

# Get predicted probabilities
pred_probs_aal <- predict(fit_aal_train, X_aal[test_idx_aal,], type = "response")

# Calculate AUC per threshold for the test set
test_aucs_aal <- sapply(1:J_aal, function(j) get_auc(pred_probs_aal[,j], y_test_mat_aal[,j]))

# Visualize Prediction Distribution
pred_df_aal <- data.frame(
  Actual = as.factor(as.vector(y_test_mat_aal)),
  Probability = as.vector(pred_probs_aal),
  Threshold = rep(paste0("T", 1:J_aal), each = length(test_idx_aal))
)

ggplot(pred_df_aal, aes(x = Actual, y = Probability, fill = Actual)) +
  geom_boxplot(alpha = 0.7) +
  facet_wrap(~Threshold) +
  theme_minimal() +
  labs(title = "AAL Model: Actual vs. Predicted Probabilities",
       subtitle = paste0("Mean Test AUC: ", round(mean(test_aucs_aal), 3)),
       y = "Predicted Probability of 'High' Decline")
```

### 4. Visualizing Anatomical Patterns
...
We group the recovered coefficients by brain lobe to see if the sparse selection captures the "circuit" structure of the simulated decline.

```{r aal_plot}
# Map regions to lobes for grouping
aal_subset <- aal[fit_aal$selected_features, ]
coefs_aal <- fit_aal$coefs_full[fit_aal$selected_features, ]
colnames(coefs_aal) <- paste0("T", 1:J_aal)
rownames(coefs_aal) <- aal_subset$label_name

# Add lobe information for visualization
aal_subset$Lobe <- "Other"
if("frontal" %in% colnames(aal)) aal_subset$Lobe[aal_subset$frontal == 1] <- "Frontal"
if("temporal" %in% colnames(aal)) aal_subset$Lobe[aal_subset$temporal == 1] <- "Temporal"
if("limbic" %in% colnames(aal)) aal_subset$Lobe[aal_subset$limbic == 1] <- "Limbic"
if("parietal" %in% colnames(aal)) aal_subset$Lobe[aal_subset$parietal == 1] <- "Parietal"

aal_melted <- melt(coefs_aal)
colnames(aal_melted) <- c("Region", "Threshold", "Weight")
# Join lobe info
aal_melted$Lobe <- aal_subset$Lobe[match(aal_melted$Region, aal_subset$label_name)]

ggplot(aal_melted, aes(x = Threshold, y = Region, fill = Weight)) +
  geom_tile() +
  facet_grid(Lobe ~ ., scales = "free_y", space = "free_y") +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
  theme_minimal() +
  labs(title = "Anatomical Signature of Cognitive Decline",
       subtitle = "Sparse Fused Ridge groups features by functional Lobe",
       x = "Disease Severity Threshold",
       y = "AAL Region") +
  theme(strip.text.y = element_text(angle = 0))
```

## Real-World Application: Blood-Brain Barrier (BBB) Permeability

In this final section, we use the **Blood-Brain Barrier** dataset from the `caret` package. This dataset contains 134 molecular descriptors for 208 compounds. The goal is to predict the log concentration ratio of the compound in the brain vs. blood (`logBBB`). 

This is a classic "real-world" problem where we want to identify which molecular properties consistently determine brain permeability across different "penetration" thresholds.

### 1. Data Preparation

```{r bbb_setup}
if (requireNamespace("caret", quietly = TRUE)) {
  data(BloodBrain, package = "caret")
  # bbbDescr contains 134 predictors, logBBB is the response
  X_bbb <- as.matrix(bbbDescr)
  y_bbb <- logBBB
  
  # Remove near-zero variance predictors to keep the model stable
  nzv <- caret::nearZeroVar(X_bbb)
  X_bbb <- X_bbb[, -nzv]
  
  M_bbb <- ncol(X_bbb)
  N_bbb <- nrow(X_bbb)
  
  # Define thresholds for 'High' vs 'Low' permeability
  thresholds_bbb <- quantile(y_bbb, probs = c(0.25, 0.5, 0.75))
  J_bbb <- length(thresholds_bbb)
  
  cat("Dataset loaded:", N_bbb, "compounds,", M_bbb, "molecular descriptors.\n")
} else {
  cat("The 'caret' package is required for this section but is not installed.\n")
}
```

### 2. Identifying Stable Molecular Signatures

We fit a fused ridge model with `topK = 10` to find the 10 most influential molecular descriptors across the permeability spectrum.

```{r bbb_fit}
if (requireNamespace("caret", quietly = TRUE)) {
  # Fit the model
  fit_bbb <- fusedRidge(X_bbb, y_bbb, thresholds_bbb, 
                       lambda1 = 0.1, lambda2 = 0.5, topK = 10)
  
  # Reconstruct the "Permeability Signature"
  coefs_bbb <- fit_bbb$coefs_full[fit_bbb$selected_features, ]
  colnames(coefs_bbb) <- paste0("T", 1:J_bbb)
  
  # Sort by overall influence
  imp_order <- order(rowMeans(abs(coefs_bbb)), decreasing = TRUE)
  coefs_bbb_sorted <- coefs_bbb[imp_order, ]
  
  print("Top Molecular Predictors of BBB Permeability:")
  print(rownames(coefs_bbb_sorted))
}
```

### 3. Evaluating Actual Predictions

We evaluate the BBB model by looking at the separation between predicted probabilities for compounds that actually crossed vs. those that did not, using a hold-out test set.

```{r bbb_eval}
if (requireNamespace("caret", quietly = TRUE)) {
  # Train/Test Split
  train_idx_bbb <- sample(1:N_bbb, size = 0.7 * N_bbb)
  test_idx_bbb <- setdiff(1:N_bbb, train_idx_bbb)
  
  # Fit sparse model
  fit_bbb_train <- fusedRidge(X_bbb[train_idx_bbb,], y_bbb[train_idx_bbb], 
                             thresholds_bbb, lambda1 = 0.1, lambda2 = 0.5, topK = 10)
  
  # Predict on test set
  pred_probs_bbb <- predict(fit_bbb_train, X_bbb[test_idx_bbb,], type = "response")
  
  # Independent Lasso Comparison
  pred_probs_lasso_bbb <- matrix(0, length(test_idx_bbb), J_bbb)
  for (j in 1:J_bbb) {
    y_train_j <- ifelse(y_bbb[train_idx_bbb] >= thresholds_bbb[j], 1, 0)
    cv_lasso <- glmnet::cv.glmnet(X_bbb[train_idx_bbb,], y_train_j, family = "binomial", alpha = 1)
    pred_probs_lasso_bbb[, j] <- as.vector(predict(cv_lasso, X_bbb[test_idx_bbb,], s = "lambda.min", type = "response"))
  }
  
  # Binarize test labels
  y_test_mat_bbb <- matrix(0, length(test_idx_bbb), J_bbb)
  for (j in 1:J_bbb) y_test_mat_bbb[,j] <- ifelse(y_bbb[test_idx_bbb] >= thresholds_bbb[j], 1, 0)
  
  # Summary Results
  test_aucs_bbb <- sapply(1:J_bbb, function(j) get_auc(pred_probs_bbb[,j], y_test_mat_bbb[,j]))
  test_aucs_lasso_bbb <- sapply(1:J_bbb, function(j) get_auc(pred_probs_lasso_bbb[,j], y_test_mat_bbb[,j]))
  
  # Prepare for Plotting
  pred_df_bbb <- data.frame(
    Actual = as.factor(as.vector(y_test_mat_bbb)),
    Probability = as.vector(pred_probs_bbb),
    Threshold = rep(paste0("Q", 1:J_bbb), each = length(test_idx_bbb))
  )
  
  # AUC Comparison Data
  auc_comp_bbb <- data.frame(
    Threshold = rep(paste0("Q", 1:J_bbb), 2),
    Method = rep(c("Fused Ridge", "Independent Lasso"), each = J_bbb),
    AUC = c(test_aucs_bbb, test_aucs_lasso_bbb)
  )
  
  # Density Plot
  p1_bbb <- ggplot(pred_df_bbb, aes(x = Probability, fill = Actual)) +
    geom_density(alpha = 0.5) +
    facet_wrap(~Threshold) +
    theme_minimal() +
    labs(title = "BBB Model: Separation of Predicted Probabilities",
         subtitle = paste0("Fused Ridge Mean Test AUC: ", round(mean(test_aucs_bbb), 3)),
         x = "Predicted Probability of Brain Penetration",
         y = "Density")
  
  # AUC Comparison Plot
  p2_bbb <- ggplot(auc_comp_bbb, aes(x = Threshold, y = AUC, fill = Method)) +
    geom_bar(stat = "identity", position = "dodge") +
    theme_minimal() +
    coord_cartesian(ylim = c(0.5, 1)) +
    labs(title = "BBB Model: AUC Comparison",
         subtitle = "Fused Ridge vs. Independent Lasso across Permeability Quantiles",
         y = "Test AUC")
  
  print(p1_bbb)
  print(p2_bbb)
}
```

### 4. Visualizing the Molecular Signature

The heatmap below reveals which molecular properties (like surface area, hydrogen bond donors, or polarizability) consistently predict crossing the blood-brain barrier.

```{r bbb_plot}
if (requireNamespace("caret", quietly = TRUE)) {
  bbb_melted <- melt(coefs_bbb_sorted)
  colnames(bbb_melted) <- c("Descriptor", "Threshold", "Weight")
  
  ggplot(bbb_melted, aes(x = Threshold, y = Descriptor, fill = Weight)) +
    geom_tile() +
    scale_fill_gradient2(low = "darkgreen", mid = "white", high = "purple") +
    theme_minimal() +
    labs(title = "Molecular Signature of BBB Permeability",
         subtitle = "Fused Ridge identifies descriptors that persist across penetration levels",
         x = "Permeability Quantile Threshold",
         y = "Molecular Descriptor")
}
```

## Conclusion

By using **Nested Cross-Validation**, we have shown that `fusedRidge` provides superior predictive performance on high-dimensional data compared to standard independent ridge regression. The fusion penalty successfully leverages the shared information across thresholds, while the two-stage `topK` approach identifies a robust, interpretable subset of features without suffering from data leakage.