Upload retrain_with_pca.py with huggingface_hub

retrain_with_pca.py

@@ -0,0 +1,279 @@
+#!/usr/bin/env python3
+"""
+Retrain all EEG + Braak models with PCA dimensionality reduction
+================================================================
+Compares: No PCA vs PCA(95%) vs PCA(99%) vs PCA+SelectKBest
+Also retrains EEG without MMSE to find pure EEG biomarkers
+
+Author: Satyawan Singh — Infonova Solutions
+"""
+import os, json, pickle, warnings
+import numpy as np
+import pandas as pd
+warnings.filterwarnings('ignore')
+
+if not hasattr(np, 'trapz'):
+    np.trapz = np.trapezoid
+
+from sklearn.decomposition import PCA
+from sklearn.preprocessing import StandardScaler, LabelEncoder
+from sklearn.model_selection import StratifiedKFold, cross_val_predict, cross_val_score
+from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier, VotingClassifier
+from sklearn.svm import SVC
+from sklearn.feature_selection import SelectKBest, f_classif
+from sklearn.metrics import accuracy_score, roc_auc_score, classification_report
+from sklearn.pipeline import Pipeline
+
+OUTPUT_DIR = '/Users/satyawansingh/Documents/alzheimer-research-complete/models/pca_analysis'
+os.makedirs(OUTPUT_DIR, exist_ok=True)
+
+results = {}
+
+# ══════════════════════════════════════════════════════════════
+# 1. DEEP EEG MODEL (515 features, 88 subjects)
+# ══════════════════════════════════════════════════════════════
+print("=" * 70)
+print("  MODEL 1: DEEP EEG — AD vs FTD vs Control (88 subjects)")
+print("=" * 70)
+
+df_eeg = pd.read_csv('/Users/satyawansingh/Documents/alzheimer-research-complete/models/eeg_deep_analysis/deep_features.csv')
+labels = np.load('/Users/satyawansingh/Documents/alzheimer-research-complete/models/eeg_deep_analysis/labels.npy', allow_pickle=True)
+
+feature_cols = [c for c in df_eeg.columns if c not in ['subject', 'group']]
+X_raw = df_eeg[feature_cols].fillna(0).values
+y = LabelEncoder().fit_transform(labels)  # AD=0, C=1, F=2
+
+# Binary: AD vs non-AD
+y_binary = (labels == 'A').astype(int)
+
+print(f"Features: {X_raw.shape[1]}, Samples: {X_raw.shape[0]}")
+print(f"Classes: AD={sum(labels=='A')}, Control={sum(labels=='C')}, FTD={sum(labels=='F')}")
+
+cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
+
+# --- With MMSE vs Without MMSE ---
+mmse_col_idx = feature_cols.index('mmse') if 'mmse' in feature_cols else None
+age_col_idx = feature_cols.index('age') if 'age' in feature_cols else None
+
+eeg_only_cols = [i for i, c in enumerate(feature_cols) if c not in ['mmse', 'age', 'sex']]
+X_eeg_only = X_raw[:, eeg_only_cols]
+
+print(f"\nWith MMSE+demographics: {X_raw.shape[1]} features")
+print(f"Pure EEG only: {X_eeg_only.shape[1]} features")
+
+for data_name, X_data in [("All features", X_raw), ("Pure EEG only", X_eeg_only)]:
+    print(f"\n{'─' * 60}")
+    print(f"  {data_name}")
+    print(f"{'─' * 60}")
+    
+    for task_name, y_task in [("3-class", y), ("AD vs non-AD", y_binary)]:
+        print(f"\n  Task: {task_name}")
+        
+        configs = {
+            'No reduction (SelectKBest 120)': ('skb', 120),
+            'PCA 95% variance': ('pca', 0.95),
+            'PCA 99% variance': ('pca', 0.99),
+            'PCA 30 components': ('pca', 30),
+            'PCA 15 components': ('pca', 15),
+            'PCA 95% + SelectKBest 20': ('pca_skb', (0.95, 20)),
+        }
+        
+        for config_name, (method, param) in configs.items():
+            scaler = StandardScaler()
+            X_scaled = scaler.fit_transform(X_data)
+            
+            if method == 'skb':
+                k = min(param, X_data.shape[1])
+                selector = SelectKBest(f_classif, k=k)
+                X_reduced = selector.fit_transform(X_scaled, y_task)
+            elif method == 'pca':
+                if isinstance(param, float):
+                    pca = PCA(n_components=param, random_state=42)
+                else:
+                    pca = PCA(n_components=min(param, X_data.shape[1], X_data.shape[0]-1), random_state=42)
+                X_reduced = pca.fit_transform(X_scaled)
+            elif method == 'pca_skb':
+                pca_param, skb_k = param
+                pca = PCA(n_components=pca_param, random_state=42)
+                X_pca = pca.fit_transform(X_scaled)
+                skb_k = min(skb_k, X_pca.shape[1])
+                selector = SelectKBest(f_classif, k=skb_k)
+                X_reduced = selector.fit_transform(X_pca, y_task)
+            
+            n_components = X_reduced.shape[1]
+            
+            # Train GradientBoosting
+            gb = GradientBoostingClassifier(n_estimators=200, max_depth=3, 
+                                             learning_rate=0.05, random_state=42)
+            
+            try:
+                y_pred = cross_val_predict(gb, X_reduced, y_task, cv=cv)
+                acc = accuracy_score(y_task, y_pred)
+                
+                if len(np.unique(y_task)) == 2:
+                    y_proba = cross_val_predict(gb, X_reduced, y_task, cv=cv, method='predict_proba')
+                    auc = roc_auc_score(y_task, y_proba[:, 1])
+                    print(f"    {config_name:40s} → {n_components:3d} dims  Acc={acc:.1%}  AUC={auc:.3f}")
+                    key = f"deep_eeg_{data_name}_{task_name}_{config_name}"
+                    results[key] = {'acc': acc, 'auc': auc, 'dims': n_components}
+                else:
+                    y_proba = cross_val_predict(gb, X_reduced, y_task, cv=cv, method='predict_proba')
+                    auc = roc_auc_score(y_task, y_proba, multi_class='ovr')
+                    print(f"    {config_name:40s} → {n_components:3d} dims  Acc={acc:.1%}  AUC={auc:.3f}")
+                    key = f"deep_eeg_{data_name}_{task_name}_{config_name}"
+                    results[key] = {'acc': acc, 'auc': auc, 'dims': n_components}
+            except Exception as e:
+                print(f"    {config_name:40s} → FAILED: {e}")
+
+
+# ══════════════════════════════════════════════════════════════
+# 2. BASIC EEG MODEL (851 features, 88 subjects) 
+# ══════════════════════════════════════════════════════════════
+print(f"\n\n{'=' * 70}")
+print("  MODEL 2: BASIC EEG CLASSIFIER (851 features, 88 subjects)")
+print("=" * 70)
+
+df_basic = pd.read_csv('/Users/satyawansingh/Documents/alzheimer-research-complete/models/eeg_ad_classifier/eeg_features.csv')
+y_basic = np.load('/Users/satyawansingh/Documents/alzheimer-research-complete/models/eeg_ad_classifier/eeg_labels.npy', allow_pickle=True)
+
+X_basic = df_basic.fillna(0).values
+le = LabelEncoder()
+y_basic_enc = le.fit_transform(y_basic)
+y_basic_binary = (y_basic == 'A').astype(int)
+
+print(f"Features: {X_basic.shape[1]}, Samples: {X_basic.shape[0]}")
+
+for task_name, y_task in [("3-class", y_basic_enc), ("AD vs non-AD", y_basic_binary)]:
+    print(f"\n  Task: {task_name}")
+    
+    configs = {
+        'No reduction (SelectKBest 100)': ('skb', 100),
+        'PCA 95% variance': ('pca', 0.95),
+        'PCA 99% variance': ('pca', 0.99),
+        'PCA 30 components': ('pca', 30),
+        'PCA 15 components': ('pca', 15),
+        'PCA 10 components': ('pca', 10),
+    }
+    
+    for config_name, (method, param) in configs.items():
+        scaler = StandardScaler()
+        X_scaled = scaler.fit_transform(X_basic)
+        
+        if method == 'skb':
+            k = min(param, X_basic.shape[1])
+            selector = SelectKBest(f_classif, k=k)
+            X_reduced = selector.fit_transform(X_scaled, y_task)
+        elif method == 'pca':
+            if isinstance(param, float):
+                pca = PCA(n_components=param, random_state=42)
+            else:
+                pca = PCA(n_components=min(param, X_basic.shape[1], X_basic.shape[0]-1), random_state=42)
+            X_reduced = pca.fit_transform(X_scaled)
+        
+        n_components = X_reduced.shape[1]
+        
+        gb = GradientBoostingClassifier(n_estimators=200, max_depth=3,
+                                         learning_rate=0.05, random_state=42)
+        try:
+            y_pred = cross_val_predict(gb, X_reduced, y_task, cv=cv)
+            acc = accuracy_score(y_task, y_pred)
+            
+            if len(np.unique(y_task)) == 2:
+                y_proba = cross_val_predict(gb, X_reduced, y_task, cv=cv, method='predict_proba')
+                auc = roc_auc_score(y_task, y_proba[:, 1])
+                print(f"    {config_name:40s} → {n_components:3d} dims  Acc={acc:.1%}  AUC={auc:.3f}")
+            else:
+                y_proba = cross_val_predict(gb, X_reduced, y_task, cv=cv, method='predict_proba')
+                auc = roc_auc_score(y_task, y_proba, multi_class='ovr')
+                print(f"    {config_name:40s} → {n_components:3d} dims  Acc={acc:.1%}  AUC={auc:.3f}")
+            
+            results[f"basic_eeg_{task_name}_{config_name}"] = {'acc': acc, 'auc': auc, 'dims': n_components}
+        except Exception as e:
+            print(f"    {config_name:40s} → FAILED: {e}")
+
+
+# ══════════════════════════════════════════════════════════════
+# 3. BRAAK PREDICTOR (141 features, 174 subjects)
+# ══════════════════════════════════════════════════════════════
+print(f"\n\n{'=' * 70}")
+print("  MODEL 3: BRAAK STAGE PREDICTOR (141 features, 174 donors)")
+print("=" * 70)
+
+df_braak = pd.read_csv('/Users/satyawansingh/Documents/alzheimer-research-complete/models/braak_predictor/donor_cell_features.csv')
+df_braak = df_braak.dropna(subset=['braak_numeric'])
+
+feature_cols_braak = [c for c in df_braak.columns if c not in 
+                      ['donor_label', 'braak_numeric', 'total_cells']]
+X_braak = df_braak[feature_cols_braak].fillna(0).values
+y_braak = df_braak['braak_numeric'].astype(int).values
+
+# 3-group: Low (0-I), Mid (II-III), High (IV-VI)
+y_braak_grouped = np.where(y_braak <= 1, 0, np.where(y_braak <= 3, 1, 2))
+
+print(f"Features: {X_braak.shape[1]}, Samples: {X_braak.shape[0]}")
+print(f"Groups: Low={sum(y_braak_grouped==0)}, Mid={sum(y_braak_grouped==1)}, High={sum(y_braak_grouped==2)}")
+
+cv_braak = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
+
+for config_name, (method, param) in {
+    'No reduction (SelectKBest 50)': ('skb', 50),
+    'PCA 95% variance': ('pca', 0.95),
+    'PCA 30 components': ('pca', 30),
+    'PCA 15 components': ('pca', 15),
+    'PCA 10 components': ('pca', 10),
+}.items():
+    scaler = StandardScaler()
+    X_scaled = scaler.fit_transform(X_braak)
+    
+    if method == 'skb':
+        k = min(param, X_braak.shape[1])
+        selector = SelectKBest(f_classif, k=k)
+        X_reduced = selector.fit_transform(X_scaled, y_braak_grouped)
+    elif method == 'pca':
+        if isinstance(param, float):
+            pca = PCA(n_components=param, random_state=42)
+        else:
+            pca = PCA(n_components=min(param, X_braak.shape[1], X_braak.shape[0]-1), random_state=42)
+        X_reduced = pca.fit_transform(X_scaled)
+    
+    n_components = X_reduced.shape[1]
+    
+    gb = GradientBoostingClassifier(n_estimators=200, max_depth=3,
+                                     learning_rate=0.05, random_state=42)
+    try:
+        y_pred = cross_val_predict(gb, X_reduced, y_braak_grouped, cv=cv_braak)
+        acc = accuracy_score(y_braak_grouped, y_pred)
+        y_proba = cross_val_predict(gb, X_reduced, y_braak_grouped, cv=cv_braak, method='predict_proba')
+        auc = roc_auc_score(y_braak_grouped, y_proba, multi_class='ovr')
+        print(f"  {config_name:40s} → {n_components:3d} dims  Acc={acc:.1%}  AUC={auc:.3f}")
+        results[f"braak_{config_name}"] = {'acc': acc, 'auc': auc, 'dims': n_components}
+    except Exception as e:
+        print(f"  {config_name:40s} → FAILED: {e}")
+
+
+# ══════════════════════════════════════════════════════════════
+# 4. PCA VARIANCE ANALYSIS
+# ══════════════════════════════════════════════════════════════
+print(f"\n\n{'=' * 70}")
+print("  PCA VARIANCE ANALYSIS — How many components needed?")
+print("=" * 70)
+
+for name, X_data in [("Deep EEG (515)", X_raw), ("Basic EEG (851)", X_basic), ("Braak (141)", X_braak)]:
+    X_s = StandardScaler().fit_transform(X_data)
+    pca_full = PCA(random_state=42)
+    pca_full.fit(X_s)
+    cumvar = np.cumsum(pca_full.explained_variance_ratio_)
+    
+    for threshold in [0.80, 0.90, 0.95, 0.99]:
+        n = np.argmax(cumvar >= threshold) + 1
+        print(f"  {name:20s}  {threshold:.0%} variance → {n:3d} components (out of {X_data.shape[1]})")
+    print()
+
+
+# Save results
+with open(os.path.join(OUTPUT_DIR, 'pca_comparison_results.json'), 'w') as f:
+    json.dump(results, f, indent=2)
+
+print(f"\n{'=' * 70}")
+print("  COMPLETE — Results saved to {OUTPUT_DIR}")
+print("=" * 70)