1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
 
def diagnose_data_for_regularization(X, y, feature_names=None):
    """
    Comprehensive data diagnostics to guide regularization choice.
    
    Parameters:
    -----------
    X : array of shape (n, p), feature matrix
    y : array of shape (n,), target vector
    feature_names : optional list of feature names
    
    Returns:
    --------
    dict : Diagnostic results with recommendations
    """
    n, p = X.shape
    
    print("=" * 70)
    print("DATA DIAGNOSTICS FOR REGULARIZATION SELECTION")
    print("=" * 70)
    
    # 1. Dimensionality Analysis
    print("\n1. DIMENSIONALITY ANALYSIS")
    print("-" * 40)
    print(f"   Samples (n): {n}")
    print(f"   Features (p): {p}")
    print(f"   Ratio n/p: {n/p:.2f}")
    
    if n/p > 10:
        dim_regime = "low-dimensional"
        dim_advice = "All methods viable; choose based on sparsity belief"
    elif n/p > 1:
        dim_regime = "moderate"
        dim_advice = "Elastic Net recommended for stability"
    else:
        dim_regime = "high-dimensional (p > n)"
        dim_advice = "Elastic Net or careful Lasso; Ridge if no sparsity needed"
    
    print(f"   Regime: {dim_regime}")
    print(f"   Advice: {dim_advice}")
    
    # 2. Correlation Structure
    print("\n2. CORRELATION STRUCTURE")
    print("-" * 40)
    
    # Compute correlation matrix
    corr_matrix = np.corrcoef(X.T)
    
    # Get upper triangle (excluding diagonal)
    upper_tri = corr_matrix[np.triu_indices_from(corr_matrix, k=1)]
    
    # Statistics
    max_corr = np.max(np.abs(upper_tri))
    mean_abs_corr = np.mean(np.abs(upper_tri))
    high_corr_pairs = np.sum(np.abs(upper_tri) > 0.7)
    very_high_corr_pairs = np.sum(np.abs(upper_tri) > 0.9)
    
    print(f"   Max |correlation|: {max_corr:.3f}")
    print(f"   Mean |correlation|: {mean_abs_corr:.3f}")
    print(f"   Pairs with |r| > 0.7: {high_corr_pairs} "
          f"({100*high_corr_pairs/len(upper_tri):.1f}%)")
    print(f"   Pairs with |r| > 0.9: {very_high_corr_pairs}")
    
    if very_high_corr_pairs > 0:
        corr_advice = "Strong grouping needed → Elastic Net (α ≤ 0.5)"
    elif high_corr_pairs > p/2:
        corr_advice = "Moderate grouping needed → Elastic Net (α ≈ 0.5-0.7)"
    elif mean_abs_corr > 0.3:
        corr_advice = "Some correlation → Elastic Net or Lasso"
    else:
        corr_advice = "Low correlation → Lasso appropriate"
    
    print(f"   Advice: {corr_advice}")
    
    # 3. Condition Number (Multicollinearity)
    print("\n3. MULTICOLLINEARITY (Condition Number)")
    print("-" * 40)
    
    XtX = X.T @ X
    eigenvalues = np.linalg.eigvalsh(XtX)
    condition_number = np.sqrt(eigenvalues.max() / (eigenvalues.min() + 1e-10))
    
    print(f"   Condition number: {condition_number:.1f}")
    
    if condition_number > 1000:
        cond_advice = "Severe multicollinearity → Ridge or Elastic Net (low α)"
    elif condition_number > 30:
        cond_advice = "Moderate multicollinearity → Elastic Net recommended"
    else:
        cond_advice = "Well-conditioned → All methods viable"
    
    print(f"   Advice: {cond_advice}")
    
    # 4. Response-Feature Relationships
    print("\n4. RESPONSE-FEATURE RELATIONSHIPS")
    print("-" * 40)
    
    # Marginal correlations with response
    marginal_corrs = np.array([np.corrcoef(X[:, j], y)[0, 1] for j in range(p)])
    
    # Count seemingly relevant features
    sig_threshold = 2 / np.sqrt(n)  # Approximate significance threshold
    n_significant = np.sum(np.abs(marginal_corrs) > sig_threshold)
    
    print(f"   Features with |r(X_j, y)| > {sig_threshold:.3f}: {n_significant}/{p}")
    print(f"   Max |marginal correlation|: {np.max(np.abs(marginal_corrs)):.3f}")
    
    if n_significant < 0.1 * p:
        sparsity_advice = "Likely sparse model → Lasso or Elastic Net (α ≥ 0.7)"
    elif n_significant < 0.3 * p:
        sparsity_advice = "Moderately sparse → Elastic Net (α ≈ 0.5)"
    else:
        sparsity_advice = "Many relevant features → Ridge or Elastic Net (α ≤ 0.3)"
    
    print(f"   Sparsity inference: {sparsity_advice}")
    
    # 5. Final Recommendation
    print("\n" + "=" * 70)
    print("FINAL RECOMMENDATION")
    print("=" * 70)
    
    # Simple scoring system
    score_lasso = 0
    score_ridge = 0
    score_enet = 0
    
    if n/p > 10:
        score_lasso += 1
    else:
        score_enet += 1
    
    if very_high_corr_pairs > 0:
        score_enet += 2
        score_ridge += 1
    elif high_corr_pairs > p/2:
        score_enet += 1
    else:
        score_lasso += 1
    
    if condition_number > 100:
        score_ridge += 2
        score_enet += 1
    elif condition_number > 30:
        score_enet += 1
    
    if n_significant < 0.1 * p:
        score_lasso += 1
        score_enet += 0.5
    else:
        score_ridge += 1
        score_enet += 0.5
    
    if score_enet >= max(score_lasso, score_ridge):
        recommended = "Elastic Net"
        alpha_rec = 0.5 - 0.1 * (very_high_corr_pairs > 0) - 0.1 * (condition_number > 100)
        alpha_rec = max(0.2, min(0.8, alpha_rec))
        print(f"   Primary: Elastic Net (suggested α ≈ {alpha_rec:.1f})")
    elif score_lasso > score_ridge:
        recommended = "Lasso"
        print(f"   Primary: Lasso")
    else:
        recommended = "Ridge"
        print(f"   Primary: Ridge Regression")
    
    print(f"\n   Scores: Lasso={score_lasso:.1f}, Ridge={score_ridge:.1f}, "
          f"Elastic Net={score_enet:.1f}")
    print("\n   Always validate with cross-validation on your specific data!")
    
    return {
        'n': n, 'p': p, 'ratio': n/p,
        'max_corr': max_corr, 'mean_corr': mean_abs_corr,
        'condition_number': condition_number,
        'n_significant': n_significant,
        'recommended': recommended
    }
 
# Example usage with simulated data
np.random.seed(42)
 
# Simulate data with correlated features
n, p = 200, 50
# Create correlation blocks
X = np.random.randn(n, p)
for i in range(0, p, 5):
    base = np.random.randn(n)
    for j in range(min(5, p - i)):
        X[:, i + j] = 0.8 * base + 0.2 * X[:, i + j]
 
# Sparse true model
beta_true = np.zeros(p)
beta_true[[0, 5, 10, 15, 20]] = [2, -1.5, 1, -0.8, 0.5]
y = X @ beta_true + 0.3 * np.random.randn(n)
 
# Run diagnostics
results = diagnose_data_for_regularization(X, y)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160

import numpy as np
from sklearn.linear_model import LassoCV, ElasticNetCV
from sklearn.model_selection import cross_val_predict
from sklearn.metrics import mean_squared_error
 
def check_lasso_sufficiency(X, y, n_bootstrap=50):
    """
    Check if pure Lasso is sufficient or if Elastic Net is needed.
    
    Tests for:
    1. Lasso selection stability across bootstraps
    2. CV performance comparison
    3. Coefficient variance
    
    Returns recommendation based on diagnostics.
    """
    np.random.seed(42)
    n, p = X.shape
    
    print("=" * 60)
    print("LASSO SUFFICIENCY CHECK")
    print("=" * 60)
    
    # 1. Selection Stability Analysis
    print("\n1. Selection Stability (Bootstrap Analysis)")
    print("-" * 40)
    
    feature_selections = np.zeros((n_bootstrap, p))
    
    for b in range(n_bootstrap):
        # Bootstrap sample
        idx = np.random.choice(n, size=n, replace=True)
        X_boot, y_boot = X[idx], y[idx]
        
        lasso = LassoCV(cv=5, random_state=b)
        lasso.fit(X_boot, y_boot)
        
        feature_selections[b] = (np.abs(lasso.coef_) > 1e-6).astype(int)
    
    # Selection frequency for each feature
    selection_freq = feature_selections.mean(axis=0)
    
    # Count features that are sometimes selected
    ever_selected = np.sum(selection_freq > 0)
    always_selected = np.sum(selection_freq > 0.95)
    inconsistent = np.sum((selection_freq > 0.1) & (selection_freq < 0.9))
    
    print(f"   Features ever selected: {ever_selected}")
    print(f"   Features always selected (>95%): {always_selected}")
    print(f"   Inconsistently selected (10-90%): {inconsistent}")
    
    stability_ok = inconsistent < 0.1 * ever_selected
    print(f"   Stability assessment: {'GOOD' if stability_ok else 'POOR - consider Elastic Net'}")
    
    # 2. CV Performance Comparison
    print("\n2. Cross-Validation Performance")
    print("-" * 40)
    
    # Lasso
    lasso = LassoCV(cv=5, random_state=42)
    lasso.fit(X, y)
    y_pred_lasso = cross_val_predict(lasso, X, y, cv=5)
    mse_lasso = mean_squared_error(y, y_pred_lasso)
    
    # Elastic Net with α = 0.5
    enet = ElasticNetCV(l1_ratio=0.5, cv=5, random_state=42)
    enet.fit(X, y)
    y_pred_enet = cross_val_predict(enet, X, y, cv=5)
    mse_enet = mean_squared_error(y, y_pred_enet)
    
    print(f"   Lasso CV MSE: {mse_lasso:.4f}")
    print(f"   Elastic Net (α=0.5) CV MSE: {mse_enet:.4f}")
    
    improvement = (mse_lasso - mse_enet) / mse_lasso * 100
    print(f"   Elastic Net improvement: {improvement:.1f}%")
    
    performance_ok = improvement < 5  # Less than 5% improvement from EN
    print(f"   Performance assessment: {'Lasso OK' if performance_ok else 'Elastic Net better'}")
    
    # 3. Coefficient Variance (Stability)
    print("\n3. Coefficient Variance")
    print("-" * 40)
    
    lasso_coefs = []
    enet_coefs = []
    
    for b in range(n_bootstrap):
        idx = np.random.choice(n, size=n, replace=True)
        X_boot, y_boot = X[idx], y[idx]
        
        lasso = LassoCV(cv=3, random_state=b)
        lasso.fit(X_boot, y_boot)
        lasso_coefs.append(lasso.coef_)
        
        enet = ElasticNetCV(l1_ratio=0.5, cv=3, random_state=b)
        enet.fit(X_boot, y_boot)
        enet_coefs.append(enet.coef_)
    
    lasso_coefs = np.array(lasso_coefs)
    enet_coefs = np.array(enet_coefs)
    
    lasso_mean_std = np.mean(np.std(lasso_coefs, axis=0))
    enet_mean_std = np.mean(np.std(enet_coefs, axis=0))
    
    print(f"   Lasso mean coefficient std: {lasso_mean_std:.4f}")
    print(f"   Elastic Net mean coefficient std: {enet_mean_std:.4f}")
    
    variance_ok = lasso_mean_std < 1.5 * enet_mean_std
    print(f"   Variance assessment: {'Lasso OK' if variance_ok else 'Elastic Net more stable'}")
    
    # Final Recommendation
    print("\n" + "=" * 60)
    print("RECOMMENDATION")
    print("=" * 60)
    
    issues = []
    if not stability_ok:
        issues.append("selection instability")
    if not performance_ok:
        issues.append("performance gap")
    if not variance_ok:
        issues.append("coefficient variance")
    
    if len(issues) == 0:
        print("   ✓ Pure Lasso is sufficient for this data")
        return "lasso"
    else:
        print(f"   ✗ Issues detected: {', '.join(issues)}")
        print("   → Recommend Elastic Net")
        return "elastic_net"
 
# Example usage
np.random.seed(42)
 
# Scenario 1: Lasso-friendly data (independent features)
print("\n" + "=" * 60)
print("TEST CASE: Independent Features (Lasso-friendly)")
print("=" * 60)
X1 = np.random.randn(300, 30)
beta1 = np.zeros(30)
beta1[:5] = [2, -1.5, 1, -0.8, 0.5]
y1 = X1 @ beta1 + 0.3 * np.random.randn(300)
 
result1 = check_lasso_sufficiency(X1, y1)
 
# Scenario 2: Correlated features (Elastic Net needed)
print("\n\n" + "=" * 60)
print("TEST CASE: Correlated Features (Elastic Net needed)")
print("=" * 60)
X2 = np.random.randn(300, 30)
for i in range(0, 30, 6):
    base = np.random.randn(300)
    for j in range(min(6, 30 - i)):
        X2[:, i + j] = 0.9 * base + 0.1 * X2[:, i + j]
 
beta2 = np.zeros(30)
beta2[[0, 6, 12, 18, 24]] = [2, -1.5, 1, -0.8, 0.5]
y2 = X2 @ beta2 + 0.3 * np.random.randn(300)
 
result2 = check_lasso_sufficiency(X2, y2)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158

import numpy as np
from sklearn.linear_model import RidgeCV, LassoCV, ElasticNetCV
from sklearn.model_selection import cross_val_score, KFold
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import matplotlib.pyplot as plt
 
def comprehensive_method_comparison(X, y, cv_folds=5):
    """
    Rigorous comparison of regularization methods using CV.
    
    Compares:
    - Ridge
    - Lasso
    - Elastic Net at various α values
    
    Reports:
    - CV performance (R² and MSE)
    - Number of selected features
    - Stability measures
    """
    np.random.seed(42)
    n, p = X.shape
    
    print("=" * 70)
    print("COMPREHENSIVE REGULARIZATION METHOD COMPARISON")
    print("=" * 70)
    print(f"Data: n={n}, p={p}")
    print(f"Cross-validation: {cv_folds}-fold")
    print()
    
    # Standardize features
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    
    # Methods to compare
    methods = {
        'Ridge': RidgeCV(cv=cv_folds),
        'Lasso': LassoCV(cv=cv_folds, random_state=42, max_iter=10000),
        'Elastic Net (α=0.3)': ElasticNetCV(l1_ratio=0.3, cv=cv_folds, 
                                             random_state=42, max_iter=10000),
        'Elastic Net (α=0.5)': ElasticNetCV(l1_ratio=0.5, cv=cv_folds,
                                             random_state=42, max_iter=10000),
        'Elastic Net (α=0.7)': ElasticNetCV(l1_ratio=0.7, cv=cv_folds,
                                             random_state=42, max_iter=10000),
        'Elastic Net (α=0.9)': ElasticNetCV(l1_ratio=0.9, cv=cv_folds,
                                             random_state=42, max_iter=10000),
    }
    
    results = {}
    
    print(f"{'Method':<25} {'CV R²':>10} {'CV MSE':>12} {'Non-zero':>10}")
    print("-" * 60)
    
    for name, model in methods.items():
        # Cross-validate
        scores_r2 = cross_val_score(model, X_scaled, y, cv=cv_folds, 
                                    scoring='r2')
        scores_mse = -cross_val_score(model, X_scaled, y, cv=cv_folds,
                                       scoring='neg_mean_squared_error')
        
        # Fit on full data for feature count
        model.fit(X_scaled, y)
        n_nonzero = np.sum(np.abs(model.coef_) > 1e-6)
        
        results[name] = {
            'r2_mean': scores_r2.mean(),
            'r2_std': scores_r2.std(),
            'mse_mean': scores_mse.mean(),
            'mse_std': scores_mse.std(),
            'n_nonzero': n_nonzero
        }
        
        print(f"{name:<25} {scores_r2.mean():>10.3f} {scores_mse.mean():>12.4f} "
              f"{n_nonzero:>10}")
    
    # Find best method
    best_method = max(results.keys(), key=lambda k: results[k]['r2_mean'])
    
    print()
    print(f"Best method by CV R²: {best_method}")
    
    # Stability analysis
    print("\n" + "-" * 70)
    print("STABILITY ANALYSIS (Bootstrap)")
    print("-" * 70)
    
    n_bootstrap = 30
    kf = KFold(n_splits=cv_folds, shuffle=True, random_state=42)
    
    stability_results = {}
    
    for name, model_class in [('Lasso', LassoCV), 
                               ('Elastic Net (α=0.5)', ElasticNetCV)]:
        coefs_list = []
        
        for b in range(n_bootstrap):
            idx = np.random.choice(n, size=n, replace=True)
            X_boot = X_scaled[idx]
            y_boot = y[idx]
            
            if name == 'Lasso':
                model = model_class(cv=3, random_state=b, max_iter=10000)
            else:
                model = model_class(l1_ratio=0.5, cv=3, random_state=b, max_iter=10000)
            
            model.fit(X_boot, y_boot)
            coefs_list.append(model.coef_)
        
        coefs_array = np.array(coefs_list)
        mean_coef_std = np.mean(np.std(coefs_array, axis=0))
        selection_freq = np.mean(np.abs(coefs_array) > 1e-6, axis=0)
        n_stable = np.sum((selection_freq > 0.9) | (selection_freq < 0.1))
        
        stability_results[name] = {
            'mean_coef_std': mean_coef_std,
            'n_stable_features': n_stable,
            'selection_freq': selection_freq
        }
        
        print(f"{name}:")
        print(f"   Mean coefficient std: {mean_coef_std:.4f}")
        print(f"   Stably selected/excluded features: {n_stable}/{p}")
    
    return results, stability_results
 
# Generate test data with structured correlations
np.random.seed(42)
n, p = 300, 100
 
# Create block correlation structure
X = np.random.randn(n, p)
for block_start in range(0, p, 10):
    block_end = min(block_start + 10, p)
    base = np.random.randn(n)
    for j in range(block_start, block_end):
        X[:, j] = 0.7 * base + 0.3 * X[:, j]
 
# Sparse true model (one feature per block for first 5 blocks)
beta_true = np.zeros(p)
for b in range(5):
    beta_true[b * 10] = 2 - 0.3 * b
 
y = X @ beta_true + 0.5 * np.random.randn(n)
 
# Run comparison
results, stability = comprehensive_method_comparison(X, y)
 
# Visualization
print("\n" + "=" * 70)
print("SUMMARY RECOMMENDATION")
print("=" * 70)
print()
print("For this data with block correlation structure and sparse signal:")
print("- Ridge achieves good prediction but includes all features")
print("- Lasso is sparse but may be unstable across bootstraps")
print("- Elastic Net (α ≈ 0.5-0.7) often provides best balance")
print("- Final choice should weight prediction vs. interpretability needs")