{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":29653,"databundleVersionId":2420395}],"dockerImageVersionId":31328,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# 🧠 Research-Grade Radiogenomics Framework\n## MRI → Radiomics → MGMT Methylation Prediction (Upgraded Pipeline)\n\n> **Dataset:** RSNA-MICCAI Brain Tumor Radiogenomic Classification  \n> **Task:** Predict MGMT promoter methylation status from multi-modal MRI  \n> **Framework:** Anti-leakage, reproducible, thesis-ready radiogenomics pipeline  \n> **Author:** [Your Name] | **Date:** [Date]\n\n---\n\n## 🔬 What is Radiogenomics?\nRadiogenomics studies the **statistical associations between quantitative imaging features (radiomics)** and **genomic/molecular characteristics** of tumors. Instead of invasive biopsies, we extract imaging biomarkers and learn to predict molecular status non-invasively.\n\n**Clinical relevance:**  \n- MGMT (O6-methylguanine-DNA methyltransferase) promoter methylation is a key prognostic biomarker in glioblastoma (GBM)  \n- MGMT-methylated tumors respond better to temozolomide chemotherapy (Stupp et al., 2005)  \n- Standard detection requires invasive biopsy — non-invasive imaging prediction is highly valuable\n\n**The Previous Notebook** : https://www.kaggle.com/code/hamzaharmanhusni/radiogenomics-research?scriptVersionId=307379674\n\n## 📋 Pipeline Overview (CRISP-DM Adapted)\n```\n1. DATA UNDERSTANDING  → Load labels, audit class balance, inspect DICOMs\n2. PATIENT-LEVEL SPLIT → Train/Val split BEFORE any preprocessing (anti-leakage)\n3. PREPROCESSING       → Correct DICOM loading (RescaleSlope, InstanceNumber sort)\n                          Percentile normalization, consistent resize\n4. ROI APPROXIMATION   → Multi-slice center crop + modality-aware intensity mask\n5. FEATURE EXTRACTION  → First-order + GLCM texture + shape features (per modality)\n6. FEATURE SELECTION   → Correlation filter → ANOVA/MI → LASSO → RF importance\n                          [All selection fitted ONLY on training data]\n7. MODELING            → Stratified K-Fold CV | LR + RF | No leakage\n8. EVALUATION          → ROC-AUC, F1, bootstrapped CI\n9. RADIOGENOMIC ANALYSIS → Statistical associations + biological interpretation\n```\n\n---\n## ⚠️ Audit Summary — Issues Fixed vs Original Notebook\n\n| # | Issue | Severity | Status |\n|---|-------|----------|--------|\n| 1 | DICOM loaded without RescaleSlope/Intercept correction | 🔴 Critical | ✅ Fixed |\n| 2 | Slice ordering by filename not InstanceNumber | 🔴 Critical | ✅ Fixed |\n| 3 | Feature selection & scaling fitted BEFORE train/val split | 🔴 Critical (Leakage) | ✅ Fixed |\n| 4 | RF importance computed on full data before CV | 🔴 Critical (Leakage) | ✅ Fixed |\n| 5 | Confusion matrix on training data, not CV folds | 🔴 Critical | ✅ Fixed |\n| 6 | LASSO used for classification (regression loss) | 🔴 Critical | ✅ Fixed → LogisticRegression L1 |\n| 7 | Per-slice z-score normalization (slice-level bias) | 🟠 High | ✅ Fixed → percentile norm |\n| 8 | No patient-level stratified split before preprocessing | 🔴 Critical | ✅ Fixed |\n| 9 | GLCM with levels=256 on float images → very slow / wrong | 🟠 High | ✅ Fixed → 64 levels |\n| 10 | Scaler re-fit after consensus feature selection (double leak) | 🔴 Critical | ✅ Fixed |\n| 11 | Missing modality NaN fill with global median (leakage) | 🟡 Medium | ✅ Fixed → train-median only |\n| 12 | No reproducibility guard for numpy global state | 🟡 Medium | ✅ Fixed |\n| 13 | `BASE_PATH` / `TRAIN_PATH` / `PALETTE` / `MODALITIES` referenced before definition | 🔴 Critical | ✅ Fixed |\n| 14 | No bootstrap confidence intervals on final metrics | 🟡 Medium | ✅ Added |\n| 15 | `statsmodels` imported mid-cell without install check | 🟡 Medium | ✅ Fixed |","metadata":{}},{"cell_type":"markdown","source":"---\n## A. 🛠️ Setup & Environment","metadata":{}},{"cell_type":"code","source":"!ls /kaggle/input/competitions","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T11:34:01.525648Z","iopub.execute_input":"2026-03-29T11:34:01.526112Z","iopub.status.idle":"2026-03-29T11:34:01.667484Z","shell.execute_reply.started":"2026-03-29T11:34:01.526060Z","shell.execute_reply":"2026-03-29T11:34:01.666130Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION A: SETUP\n# Install any missing packages first (Kaggle usually has these)\n# ─────────────────────────────────────────────────────────────\nimport subprocess, sys\n\ndef install(pkg):\n    subprocess.check_call([sys.executable, '-m', 'pip', 'install', pkg, '-q'])\n\nfor pkg in ['pydicom', 'statsmodels']:\n    install(pkg)\n\n# ── Core Libraries ──────────────────────────────────────────\nimport os\nimport warnings\nimport random\nimport hashlib\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec\nimport seaborn as sns\nimport cv2\nwarnings.filterwarnings('ignore')\n\n# ── Medical Imaging ──────────────────────────────────────────\nimport pydicom\nfrom pydicom.errors import InvalidDicomError\n\n# ── Image Processing ─────────────────────────────────────────\nfrom skimage import measure\nfrom skimage.feature import graycomatrix, graycoprops\nfrom scipy import stats\n\n# ── Machine Learning ─────────────────────────────────────────\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.dummy import DummyClassifier\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.feature_selection import SelectKBest, f_classif, mutual_info_classif\nfrom sklearn.model_selection import StratifiedKFold, cross_validate, StratifiedShuffleSplit\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.metrics import (\n    roc_auc_score, accuracy_score,\n    confusion_matrix, ConfusionMatrixDisplay,\n    roc_curve, classification_report, f1_score\n)\nfrom statsmodels.stats.multitest import multipletests\nfrom collections import Counter\n\n# ── Reproducibility ──────────────────────────────────────────\n# CRITICAL: Set ALL random seeds for full reproducibility\nSEED = 42\nnp.random.seed(SEED)\nrandom.seed(SEED)\nos.environ['PYTHONHASHSEED'] = str(SEED)\n\n# ── Dataset Paths ─────────────────────────────────────────────\n# FIX: BASE_PATH and TRAIN_PATH must be defined BEFORE any use\n# Adjust this to your Kaggle input path:\nBASE_PATH  = '/kaggle/input/competitions/rsna-miccai-brain-tumor-radiogenomic-classification'\nTRAIN_PATH = os.path.join(BASE_PATH, 'train')\n\n# ── Experiment Constants ──────────────────────────────────────\n# FIX: MODALITIES and PALETTE defined here, not scattered across cells\nMODALITIES   = ['FLAIR', 'T1w', 'T1wCE', 'T2w']\nTARGET_SIZE  = (128, 128)\nN_SLICES     = 5      # Number of central slices to aggregate per patient\nN_PATIENTS   = 50     # Increase to None for full dataset run\nCV_FOLDS     = 5      # Stratified K-Fold splits\nVAL_FRACTION = 0.20   # Held-out validation set fraction\nTOP_K        = 20     # Top-K features per selection method\n\n# ── Styling ──────────────────────────────────────────────────\nPALETTE = ['#7b68ee', '#ff6b9d', '#ffd93d', '#6bcb77', '#4d96ff', '#ff6b6b']\nplt.rcParams.update({\n    'figure.facecolor': '#0f0f1a',\n    'axes.facecolor':   '#1a1a2e',\n    'axes.edgecolor':   '#444466',\n    'axes.labelcolor':  '#e0e0ff',\n    'xtick.color':      '#b0b0cc',\n    'ytick.color':      '#b0b0cc',\n    'text.color':       '#e0e0ff',\n    'grid.color':       '#2a2a4a',\n    'grid.linestyle':   '--',\n    'grid.alpha':        0.5,\n    'font.family':      'DejaVu Sans',\n    'font.size':         11,\n    'figure.dpi':        100,\n})\n\nprint('✅ Environment setup complete')\nprint(f'   SEED={SEED} | TARGET_SIZE={TARGET_SIZE} | N_SLICES={N_SLICES}')\nprint(f'   MODALITIES: {MODALITIES}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T11:34:27.670553Z","iopub.execute_input":"2026-03-29T11:34:27.671047Z","iopub.status.idle":"2026-03-29T11:34:41.172001Z","shell.execute_reply.started":"2026-03-29T11:34:27.670972Z","shell.execute_reply":"2026-03-29T11:34:41.170775Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## B. 📊 Data Understanding\nLoad clinical labels, audit class distribution, verify data integrity.","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION B.1: LOAD LABELS & DATA AUDIT\n# ─────────────────────────────────────────────────────────────\n\ntrain_df = pd.read_csv(os.path.join(BASE_PATH, 'train_labels.csv'))\nprint('=== Training Labels ===')\nprint(train_df.head(10))\nprint(f'\\nShape: {train_df.shape}')\n\n# ── Missing Value Audit ───────────────────────────────────────\nprint('\\n=== Missing Value Audit ===')\nprint(train_df.isna().sum())\n\n# Drop rows with missing MGMT label\nbefore = len(train_df)\ntrain_df = train_df.dropna(subset=['MGMT_value']).reset_index(drop=True)\ntrain_df['MGMT_value'] = train_df['MGMT_value'].astype(int)\nprint(f'\\nDropped {before - len(train_df)} rows with missing MGMT label')\nprint(f'Final dataset: {len(train_df)} patients')\n\n# ── Class Distribution ────────────────────────────────────────\nprint('\\n=== Class Distribution ===')\ncounts = train_df['MGMT_value'].value_counts().sort_index()\nprint(counts)\nmethylated_pct = counts[1] / counts.sum() * 100\nprint(f'\\n📊 Methylated: {methylated_pct:.1f}% | Unmethylated: {100-methylated_pct:.1f}%')\nprint('💡 Note: ~balanced; use stratified splits and report F1 as primary metric')\n\n# ── Verify DICOM folder existence for all patients ────────────\nprint('\\n=== DICOM Folder Audit ===')\nmissing_folders = []\nfor pid in train_df['BraTS21ID']:\n    patient_dir = os.path.join(TRAIN_PATH, str(pid).zfill(5))\n    if not os.path.isdir(patient_dir):\n        missing_folders.append(pid)\nprint(f'Patients with missing DICOM folder: {len(missing_folders)}')\nif missing_folders:\n    print(f'IDs: {missing_folders[:5]} ...')\n    # Remove patients with no data\n    train_df = train_df[~train_df['BraTS21ID'].isin(missing_folders)].reset_index(drop=True)\n    print(f'Remaining: {len(train_df)} patients')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T11:35:50.755772Z","iopub.execute_input":"2026-03-29T11:35:50.756727Z","iopub.status.idle":"2026-03-29T11:35:56.895710Z","shell.execute_reply.started":"2026-03-29T11:35:50.756684Z","shell.execute_reply":"2026-03-29T11:35:56.894834Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION B.2: CLASS DISTRIBUTION VISUALIZATION\n# ─────────────────────────────────────────────────────────────\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\nfig.suptitle('MGMT Methylation Status — Class Distribution',\n             fontsize=14, fontweight='bold', color='#e0e0ff')\n\ncounts = train_df['MGMT_value'].value_counts().sort_index()\n\n# Bar chart\nbars = axes[0].bar(\n    ['Unmethylated (0)', 'Methylated (1)'],\n    counts.values,\n    color=PALETTE[:2], edgecolor='white', linewidth=0.5\n)\nfor bar, val in zip(bars, counts.values):\n    axes[0].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 2,\n                 str(val), ha='center', va='bottom', fontweight='bold', color='white')\naxes[0].set_title('Absolute Count', color='#e0e0ff')\naxes[0].set_ylabel('Number of Patients')\naxes[0].grid(axis='y', alpha=0.3)\n\n# Pie chart\naxes[1].pie(\n    counts.values,\n    labels=['Unmethylated\\n(MGMT−)', 'Methylated\\n(MGMT+)'],\n    colors=PALETTE[:2],\n    autopct='%1.1f%%',\n    startangle=90,\n    textprops={'color': 'white', 'fontsize': 12}\n)\naxes[1].set_title('Proportion', color='#e0e0ff')\n\nplt.tight_layout()\nplt.savefig('class_distribution.png', dpi=150, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T11:37:05.559927Z","iopub.execute_input":"2026-03-29T11:37:05.560822Z","iopub.status.idle":"2026-03-29T11:37:06.226905Z","shell.execute_reply.started":"2026-03-29T11:37:05.560781Z","shell.execute_reply":"2026-03-29T11:37:06.225783Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## C. ✂️ Patient-Level Train/Validation Split\n\n**⚠️ CRITICAL ANTI-LEAKAGE STEP — This MUST happen before ANY preprocessing, feature extraction, or normalization.**\n\n### Why patient-level split?\n- If the same patient's data appears in both train and validation, the model may learn patient-specific artifacts rather than generalizable features.\n- Normalization statistics (mean, std, percentiles) MUST be computed only on training patients; applying them to validation is correct and necessary but they must NOT be recomputed using validation data.\n\n### Why before preprocessing?\n- Any normalization, imputation, or scaling fitted on the full dataset constitutes data leakage — the model indirectly \"sees\" validation labels during training.","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION C: PATIENT-LEVEL STRATIFIED SPLIT\n# FIX: Split FIRST — this is the most critical anti-leakage step.\n# ─────────────────────────────────────────────────────────────\n\n# Use StratifiedShuffleSplit to maintain class balance in both sets\nsplitter = StratifiedShuffleSplit(\n    n_splits=1, test_size=VAL_FRACTION, random_state=SEED\n)\n\nX_idx = np.arange(len(train_df))\ny_all = train_df['MGMT_value'].values\n\ntrain_idx, val_idx = next(splitter.split(X_idx, y_all))\n\ntrain_meta = train_df.iloc[train_idx].reset_index(drop=True)\nval_meta   = train_df.iloc[val_idx].reset_index(drop=True)\n\nprint('=== Patient-Level Split ===')\nprint(f'Total patients:      {len(train_df)}')\nprint(f'Training patients:   {len(train_meta)}  '\n      f'({train_meta[\"MGMT_value\"].mean()*100:.1f}% methylated)')\nprint(f'Validation patients: {len(val_meta)}  '\n      f'({val_meta[\"MGMT_value\"].mean()*100:.1f}% methylated)')\nprint('\\n✅ No patient appears in both splits — confirmed.')\noverlap = set(train_meta['BraTS21ID']) & set(val_meta['BraTS21ID'])\nassert len(overlap) == 0, f'LEAKAGE DETECTED: {overlap}'\n\n# Apply N_PATIENTS subset AFTER split, within each group (stratified)\nif N_PATIENTS is not None:\n    n_train = int(N_PATIENTS * (1 - VAL_FRACTION))\n    n_val   = int(N_PATIENTS * VAL_FRACTION)\n    train_meta = (\n        train_meta.groupby('MGMT_value', group_keys=False)\n        .apply(lambda x: x.sample(min(n_train // 2, len(x)), random_state=SEED))\n        .reset_index(drop=True)\n    )\n    val_meta = (\n        val_meta.groupby('MGMT_value', group_keys=False)\n        .apply(lambda x: x.sample(min(n_val // 2, len(x)), random_state=SEED))\n        .reset_index(drop=True)\n    )\n    print(f'\\n📊 Subset mode: {len(train_meta)} train | {len(val_meta)} val patients')\n\ny_train_meta = train_meta['MGMT_value'].values\ny_val_meta   = val_meta['MGMT_value'].values","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T11:37:49.696012Z","iopub.execute_input":"2026-03-29T11:37:49.697190Z","iopub.status.idle":"2026-03-29T11:37:49.725051Z","shell.execute_reply.started":"2026-03-29T11:37:49.697143Z","shell.execute_reply":"2026-03-29T11:37:49.723659Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## D. 🔧 Correct DICOM Loading & Preprocessing\n\n### Bugs fixed in original:\n1. **RescaleSlope/Intercept not applied** — raw `pixel_array` is in \"stored\" values; the true Hounsfield Units (or MRI signal) require `pixel_value * RescaleSlope + RescaleIntercept`. Skipping this step means different scanners produce incomparable intensities.\n2. **Slice ordering by filename** — DICOM filenames are arbitrary and do NOT guarantee anatomical order. The correct attribute is `InstanceNumber` (or `SliceLocation`). Wrong order leads to incorrect spatial aggregation.\n3. **Per-slice z-score normalization** — computing z-score on each slice independently destroys inter-slice intensity relationships and introduces slice-level statistical bias. Percentile-based normalization (e.g., clamp to [1st, 99th] percentile per volume) is more robust and radiomics-valid.","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION D: CORRECT DICOM LOADING\n# ─────────────────────────────────────────────────────────────\n\ndef load_dicom_sorted_volume(folder: str) -> np.ndarray | None:\n    \"\"\"\n    Load all DICOM slices from a folder and return a spatially-ordered 3D volume.\n\n    FIX 1: Apply RescaleSlope and RescaleIntercept correctly.\n             pixel_value_true = stored_value * RescaleSlope + RescaleIntercept\n             This is mandatory for cross-scanner comparability.\n\n    FIX 2: Sort slices by InstanceNumber (DICOM tag), NOT by filename.\n             Filenames are arbitrary; InstanceNumber is the anatomical z-position.\n\n    Returns:\n        3D numpy array of shape (n_slices, H, W) in float32, or None if loading fails.\n    \"\"\"\n    if not os.path.exists(folder):\n        return None\n\n    dcm_paths = [os.path.join(folder, f)\n                 for f in os.listdir(folder) if f.endswith('.dcm')]\n    if not dcm_paths:\n        return None\n\n    slices = []\n    for path in dcm_paths:\n        try:\n            ds = pydicom.dcmread(path)\n            # FIX 2: Extract InstanceNumber for sorting; fall back to filename\n            instance_num = int(getattr(ds, 'InstanceNumber', 0))\n            # FIX 1: Apply RescaleSlope and RescaleIntercept\n            slope     = float(getattr(ds, 'RescaleSlope',     1.0))\n            intercept = float(getattr(ds, 'RescaleIntercept', 0.0))\n            img = ds.pixel_array.astype(np.float32) * slope + intercept\n            slices.append((instance_num, img))\n        except (InvalidDicomError, Exception):\n            continue\n\n    if not slices:\n        return None\n\n    # FIX 2: Sort by InstanceNumber (anatomical order)\n    slices.sort(key=lambda x: x[0])\n    volume = np.stack([s[1] for s in slices], axis=0)  # (N_slices, H, W)\n    return volume\n\n\ndef normalize_volume_percentile(volume: np.ndarray,\n                                p_low: float = 1.0,\n                                p_high: float = 99.0) -> np.ndarray:\n    \"\"\"\n    FIX 3: Percentile-based volume-level normalization.\n\n    WHY: Per-slice z-score normalization (original approach) produces a different\n    statistical distribution per slice, destroying the inter-slice intensity\n    relationships that encode tumor heterogeneity. Volume-level percentile clipping\n    preserves relative intensities and is robust to outlier voxels.\n\n    This normalizes to [0, 1] using the 1st and 99th percentile of the whole volume,\n    which is standard in MRI radiomics (Zwanenburg et al., 2020).\n    \"\"\"\n    lo = np.percentile(volume, p_low)\n    hi = np.percentile(volume, p_high)\n    if hi == lo:\n        return np.zeros_like(volume)\n    vol_norm = np.clip(volume, lo, hi)\n    vol_norm = (vol_norm - lo) / (hi - lo)\n    return vol_norm.astype(np.float32)\n\n\ndef resize_slice(img: np.ndarray, target_size: tuple = TARGET_SIZE) -> np.ndarray:\n    \"\"\"Resize a 2D slice to target_size using area interpolation (anti-aliasing).\"\"\"\n    return cv2.resize(img, target_size, interpolation=cv2.INTER_AREA).astype(np.float32)\n\n\ndef get_central_slices(volume: np.ndarray, n: int = N_SLICES) -> list:\n    \"\"\"\n    Select n central slices from a 3D volume.\n\n    Rationale: GBM tumors are predominantly found in the central regions\n    of the brain volume; central slices maximize the probability of capturing\n    the tumor mass (Bakas et al., 2017).\n    \"\"\"\n    total = volume.shape[0]\n    mid   = total // 2\n    half  = n // 2\n    start = max(0, mid - half)\n    end   = min(total, start + n)\n    return [volume[i] for i in range(start, end)]\n\n\ndef load_patient_volume(\n    patient_id: str | int,\n    modalities: list = MODALITIES,\n    n_slices:   int  = N_SLICES\n) -> dict:\n    \"\"\"\n    Full pipeline: load → sort by InstanceNumber → percentile normalize → resize → aggregate.\n\n    Returns:\n        dict: modality → averaged preprocessed slice (H×W float32), or None if failed.\n    \"\"\"\n    result = {}\n    for mod in modalities:\n        folder = os.path.join(TRAIN_PATH, str(patient_id).zfill(5), mod)\n        vol = load_dicom_sorted_volume(folder)\n        if vol is None:\n            result[mod] = None\n            continue\n        # Volume-level normalization (FIX 3)\n        vol_norm = normalize_volume_percentile(vol)\n        # Central slices (correct spatial ordering via InstanceNumber)\n        central  = get_central_slices(vol_norm, n=n_slices)\n        # Resize each slice, then average (pseudo-3D representation)\n        resized  = [resize_slice(s) for s in central]\n        result[mod] = np.mean(resized, axis=0) if resized else None\n    return result\n\n\nprint('✅ Corrected DICOM loading functions defined')\nprint('   Fixes applied: RescaleSlope/Intercept | InstanceNumber sort | Volume percentile norm')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T11:38:10.380843Z","iopub.execute_input":"2026-03-29T11:38:10.381356Z","iopub.status.idle":"2026-03-29T11:38:10.402362Z","shell.execute_reply.started":"2026-03-29T11:38:10.381320Z","shell.execute_reply":"2026-03-29T11:38:10.401137Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION D.2: PREPROCESSING VISUALIZATION\n# ─────────────────────────────────────────────────────────────\n\n# Use a training patient for visualization (no leakage risk here)\ntest_pid = train_meta['BraTS21ID'].iloc[0]\n\n# Load raw (uncorrected — for comparison only)\nraw_folder = os.path.join(TRAIN_PATH, str(test_pid).zfill(5), 'FLAIR')\nraw_vol    = None\nif os.path.exists(raw_folder):\n    dcm_files_raw = sorted([f for f in os.listdir(raw_folder) if f.endswith('.dcm')])\n    if dcm_files_raw:\n        mid_raw = len(dcm_files_raw) // 2\n        ds_raw  = pydicom.dcmread(os.path.join(raw_folder, dcm_files_raw[mid_raw]))\n        raw_vol = ds_raw.pixel_array.astype(np.float32)\n\n# Load corrected\ncorrected_volume = load_patient_volume(test_pid)\nproc_img = corrected_volume.get('FLAIR')\n\nif raw_vol is not None and proc_img is not None:\n    fig, axes = plt.subplots(1, 3, figsize=(16, 4))\n    fig.suptitle(f'Preprocessing Comparison — Patient {test_pid} (FLAIR)',\n                 fontsize=13, fontweight='bold', color='#e0e0ff')\n\n    # Normalize raw for display only\n    raw_disp = (raw_vol - raw_vol.min()) / (raw_vol.max() - raw_vol.min() + 1e-8)\n    axes[0].imshow(raw_disp, cmap='magma')\n    axes[0].set_title('1. Raw DICOM (no correction)', color='#e0e0ff')\n    axes[0].axis('off')\n\n    axes[1].imshow(proc_img, cmap='magma')\n    axes[1].set_title('2. Corrected: RescaleSlope+Intercept\\n+ InstanceNumber sort + Percentile norm',\n                      color='#e0e0ff', fontsize=9)\n    axes[1].axis('off')\n\n    axes[2].hist(raw_disp.flatten(), bins=60, alpha=0.6, color=PALETTE[0],\n                 label='Raw (display norm)', density=True)\n    axes[2].hist(proc_img.flatten(), bins=60, alpha=0.6, color=PALETTE[1],\n                 label='Corrected (percentile norm)', density=True)\n    axes[2].set_title('Intensity Histogram Comparison', color='#e0e0ff')\n    axes[2].legend(fontsize=8)\n    axes[2].set_xlabel('Normalized Intensity')\n    axes[2].grid(True, alpha=0.3)\n\n    plt.tight_layout()\n    plt.savefig('preprocessing.png', dpi=150, bbox_inches='tight')\n    plt.show()\n\nprint('✅ Preprocessing visualization complete')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T11:38:28.675662Z","iopub.execute_input":"2026-03-29T11:38:28.677146Z","iopub.status.idle":"2026-03-29T11:38:50.875660Z","shell.execute_reply.started":"2026-03-29T11:38:28.677084Z","shell.execute_reply":"2026-03-29T11:38:50.874108Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## E. 🎯 Tumor Region Approximation (Pseudo-ROI)\n\n**Real radiomics requires expert tumor segmentation** (e.g., BraTS annotations). Our pragmatic approach uses two complementary heuristics:\n\n1. **Center crop** — 60% central region where GBM tumors predominantly occur  \n2. **Modality-aware intensity mask** — FLAIR/T2w hyperintense regions (edema), T1wCE enhancing regions\n\n**Trade-off:** This approach may include normal brain tissue and cannot distinguish tumor core from edema. For production research, deep learning segmenters (nnU-Net, SynthSeg) should be used.","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION E: TUMOR REGION APPROXIMATION\n# ─────────────────────────────────────────────────────────────\n\ndef get_center_crop(img: np.ndarray, crop_frac: float = 0.6) -> np.ndarray:\n    \"\"\"\n    Extract central crop of the image.\n    Brain tumors in GBM patients predominantly occur in central/peri-ventricular regions.\n    \"\"\"\n    h, w  = img.shape\n    ch, cw = int(h * crop_frac), int(w * crop_frac)\n    top  = (h - ch) // 2\n    left = (w - cw) // 2\n    return img[top:top+ch, left:left+cw]\n\n\ndef get_intensity_mask(img: np.ndarray, percentile: float = 75.0) -> np.ndarray:\n    \"\"\"\n    Binary mask: pixels brighter than the given percentile.\n    In FLAIR, tumor/edema appear hyperintense relative to normal brain.\n    Returns mask of same spatial dimensions as input.\n    \"\"\"\n    threshold = np.percentile(img, percentile)\n    return (img > threshold).astype(np.float32)\n\n\ndef apply_roi(\n    img:           np.ndarray,\n    modality:      str   = 'FLAIR',\n    use_intensity: bool  = True,\n    crop_frac:     float = 0.6\n) -> np.ndarray:\n    \"\"\"\n    Combined ROI: center crop followed by modality-aware intensity mask.\n\n    Modality-specific thresholds:\n    - FLAIR/T2w : 75th percentile (hyperintense edema/tumor)\n    - T1wCE     : 70th percentile (enhancing tumor, less edema)\n    - T1w       : center crop only (homogeneous anatomy, less tumor contrast)\n    \"\"\"\n    cropped = get_center_crop(img, crop_frac=crop_frac)\n    if use_intensity and modality in ['FLAIR', 'T2w', 'T1wCE']:\n        pct = 70 if modality == 'T1wCE' else 75\n        mask = get_intensity_mask(cropped, pct)\n        return cropped * mask\n    return cropped\n\n\n# ── Visualize ROI approximation (training patient only) ───────\nif proc_img is not None:\n    center_crop    = get_center_crop(proc_img)\n    intensity_mask = get_intensity_mask(proc_img)\n    roi_img        = apply_roi(proc_img, modality='FLAIR')\n\n    fig, axes = plt.subplots(1, 4, figsize=(18, 4))\n    fig.suptitle('Tumor Region Approximation Strategies (FLAIR)',\n                 fontsize=13, fontweight='bold', color='#e0e0ff')\n\n    titles = ['Original Preprocessed', 'Center Crop (60%)',\n              'Intensity Mask (>75th pct)', 'Combined ROI']\n    imgs   = [proc_img, center_crop, intensity_mask, roi_img]\n    cmaps  = ['magma', 'magma', 'binary_r', 'magma']\n\n    for ax, title, im, cmap in zip(axes, titles, imgs, cmaps):\n        ax.imshow(im, cmap=cmap)\n        ax.set_title(title, color='#e0e0ff', fontsize=10)\n        ax.axis('off')\n\n    plt.tight_layout()\n    plt.savefig('roi_approximation.png', dpi=150, bbox_inches='tight')\n    plt.show()\n\nprint('✅ ROI approximation functions defined')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T11:40:05.425520Z","iopub.execute_input":"2026-03-29T11:40:05.426203Z","iopub.status.idle":"2026-03-29T11:40:06.127912Z","shell.execute_reply.started":"2026-03-29T11:40:05.426164Z","shell.execute_reply":"2026-03-29T11:40:06.126300Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## F. 🔬 Radiomics Feature Extraction\n\nWe extract three categories of quantitative imaging features:\n\n| Category | Features | Biological Relevance |\n|---|---|---|\n| **Intensity** | Mean, SD, skewness, kurtosis, energy, entropy, P25/P75, IQR | Tumor cellularity, necrosis, heterogeneity |\n| **Texture (GLCM)** | Contrast, correlation, homogeneity, dissimilarity, energy | Tissue heterogeneity, architectural complexity |\n| **Shape/Morphology** | Area fraction, solidity, eccentricity, perimeter, compactness | Tumor invasiveness, geometry |\n\n**Fix:** GLCM now uses 64 levels (not 256) — dramatically faster and radiomics-standard.","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION F: RADIOMICS FEATURE EXTRACTION\n# ─────────────────────────────────────────────────────────────\n\ndef extract_intensity_features(roi: np.ndarray, prefix: str = '') -> dict:\n    \"\"\"\n    First-order intensity statistics on ROI foreground pixels.\n\n    Excludes background (zero/masked pixels) to focus on tumor signal.\n    These features quantify the distribution of signal intensities within\n    the tumor region — a direct proxy for tissue heterogeneity.\n    \"\"\"\n    flat = roi.flatten()\n    flat = flat[flat > 0]  # Exclude background\n    if len(flat) < 10:\n        return {f'{prefix}_{k}': np.nan for k in\n                ['mean','std','skewness','kurtosis','energy','entropy','p25','p75','iqr']}\n    p25, p75 = np.percentile(flat, [25, 75])\n    hist, _  = np.histogram(flat, bins=64, density=True)\n    hist     = hist[hist > 0]\n    return {\n        f'{prefix}_mean':     float(np.mean(flat)),\n        f'{prefix}_std':      float(np.std(flat)),\n        f'{prefix}_skewness': float(stats.skew(flat)),\n        f'{prefix}_kurtosis': float(stats.kurtosis(flat)),\n        f'{prefix}_energy':   float(np.sum(flat ** 2)),\n        f'{prefix}_entropy':  float(-np.sum(hist * np.log2(hist + 1e-10))),\n        f'{prefix}_p25':      float(p25),\n        f'{prefix}_p75':      float(p75),\n        f'{prefix}_iqr':      float(p75 - p25),\n    }\n\n\ndef extract_glcm_features(roi: np.ndarray, prefix: str = '') -> dict:\n    \"\"\"\n    Gray-Level Co-occurrence Matrix (GLCM) texture features.\n\n    FIX: Use 64 levels instead of 256. GLCM with 256 levels on float images\n    is astronomically slow and statistically sparse. 64 levels is the\n    radiomics community standard (Zwanenburg et al., 2020, IBSI guidelines).\n\n    Computed at 4 orientations (0°, 45°, 90°, 135°) and averaged for\n    rotation invariance. Two distances (1, 3 pixels) capture multi-scale texture.\n    \"\"\"\n    # Quantize to 64 levels for GLCM (standard radiomics practice)\n    img_q = (roi * 63).astype(np.uint8)\n    if img_q.max() == 0:\n        return {f'{prefix}_glcm_{k}': np.nan for k in\n                ['contrast','correlation','homogeneity','dissimilarity','energy']}\n\n    glcm = graycomatrix(\n        img_q, distances=[1, 3],\n        angles=[0, np.pi/4, np.pi/2, 3*np.pi/4],\n        levels=64, symmetric=True, normed=True\n    )\n    props = ['contrast', 'correlation', 'homogeneity', 'dissimilarity', 'energy']\n    return {\n        f'{prefix}_glcm_{prop}': float(graycoprops(glcm, prop).mean())\n        for prop in props\n    }\n\n\ndef extract_shape_features(roi: np.ndarray, prefix: str = '') -> dict:\n    \"\"\"\n    Morphological features from the ROI mask.\n\n    These proxy the geometric properties of the hyperintense tumor region.\n    Note: True shape features require exact segmentation masks.\n    These are useful as approximate shape biomarkers.\n    \"\"\"\n    binary = (roi > roi.mean()).astype(np.uint8)\n    area_frac = float(binary.sum() / binary.size)\n    regions = measure.regionprops(binary)\n    if not regions:\n        return {\n            f'{prefix}_area_frac':    area_frac,\n            f'{prefix}_solidity':     np.nan,\n            f'{prefix}_eccentricity': np.nan,\n            f'{prefix}_perimeter':    np.nan,\n            f'{prefix}_compactness':  np.nan,\n        }\n    r = max(regions, key=lambda x: x.area)\n    perimeter   = r.perimeter if r.perimeter > 0 else 1.0\n    compactness = (4 * np.pi * r.area) / (perimeter ** 2)\n    return {\n        f'{prefix}_area_frac':    area_frac,\n        f'{prefix}_solidity':     float(r.solidity),\n        f'{prefix}_eccentricity': float(r.eccentricity),\n        f'{prefix}_perimeter':    float(perimeter),\n        f'{prefix}_compactness':  float(compactness),\n    }\n\n\ndef extract_all_features(patient_id: str | int) -> dict:\n    \"\"\"\n    Master feature extraction for one patient.\n    Runs the full pipeline: load → preprocess → ROI → extract features.\n    \"\"\"\n    volume      = load_patient_volume(patient_id)\n    all_features = {'BraTS21ID': patient_id}\n\n    for mod in MODALITIES:\n        img = volume.get(mod)\n        if img is None:\n            # Fill missing modality with NaN placeholders\n            dummy = {f'{mod}_mean': np.nan, f'{mod}_std': np.nan,\n                     f'{mod}_skewness': np.nan, f'{mod}_kurtosis': np.nan,\n                     f'{mod}_energy': np.nan, f'{mod}_entropy': np.nan,\n                     f'{mod}_p25': np.nan, f'{mod}_p75': np.nan, f'{mod}_iqr': np.nan,\n                     f'{mod}_glcm_contrast': np.nan, f'{mod}_glcm_correlation': np.nan,\n                     f'{mod}_glcm_homogeneity': np.nan, f'{mod}_glcm_dissimilarity': np.nan,\n                     f'{mod}_glcm_energy': np.nan,\n                     f'{mod}_area_frac': np.nan, f'{mod}_solidity': np.nan,\n                     f'{mod}_eccentricity': np.nan, f'{mod}_perimeter': np.nan,\n                     f'{mod}_compactness': np.nan}\n            all_features.update(dummy)\n            continue\n        roi = apply_roi(img, modality=mod)\n        all_features.update(extract_intensity_features(roi, prefix=mod))\n        all_features.update(extract_glcm_features(roi,      prefix=mod))\n        all_features.update(extract_shape_features(roi,     prefix=mod))\n\n    return all_features\n\n\nprint(f'🔬 Feature extraction functions defined')\nprint(f'   Features per modality: 9 (intensity) + 5 (GLCM) + 5 (shape) = 19')\nprint(f'   Total features (4 modalities): 76')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T11:40:32.627910Z","iopub.execute_input":"2026-03-29T11:40:32.628965Z","iopub.status.idle":"2026-03-29T11:40:32.652233Z","shell.execute_reply.started":"2026-03-29T11:40:32.628924Z","shell.execute_reply":"2026-03-29T11:40:32.650990Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION F.2: EXTRACT FEATURES — TRAIN AND VALIDATION SEPARATELY\n# CRITICAL: Extract from train and val subsets independently.\n# No statistics from validation data must influence training preprocessing.\n# ─────────────────────────────────────────────────────────────\n\ndef run_extraction(meta_df: pd.DataFrame, label: str = 'train') -> pd.DataFrame:\n    \"\"\"Extract features for all patients in meta_df and return a feature DataFrame.\"\"\"\n    records = []\n    failed  = []\n    n = len(meta_df)\n    print(f'⏳ Extracting features for {n} {label} patients...')\n\n    for i, row in meta_df.iterrows():\n        pid = row['BraTS21ID']\n        try:\n            feats = extract_all_features(pid)\n            feats['MGMT_value'] = row['MGMT_value']\n            records.append(feats)\n        except Exception as e:\n            failed.append((pid, str(e)))\n\n        if (len(records) + len(failed)) % 10 == 0:\n            done = len(records) + len(failed)\n            print(f'  → {done}/{n} | OK:{len(records)} | Failed:{len(failed)}')\n\n    df = pd.DataFrame(records)\n    print(f'\\n✅ {label.upper()} extraction complete: {df.shape} | Failed: {len(failed)}')\n    if failed:\n        print(f'   Failed IDs: {[f[0] for f in failed[:5]]}')\n    return df\n\n\nfeat_train_df = run_extraction(train_meta, label='train')\nfeat_val_df   = run_extraction(val_meta,   label='val')\n\nprint(f'\\n📋 Train feature matrix: {feat_train_df.shape}')\nprint(f'📋 Val   feature matrix: {feat_val_df.shape}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T11:41:13.932786Z","iopub.execute_input":"2026-03-29T11:41:13.934251Z","iopub.status.idle":"2026-03-29T11:52:23.895338Z","shell.execute_reply.started":"2026-03-29T11:41:13.934204Z","shell.execute_reply":"2026-03-29T11:52:23.893338Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## G. 🔍 Anti-Leakage Feature Selection\n\n**Critical rule:** ALL selection steps (correlation filter, ANOVA, mutual information, LASSO, RF importance) are fitted ONLY on training data. The resulting feature set is then applied to validation data.\n\nAdditionally:\n- **FIX:** `StandardScaler` is fitted on training data only, then applied to both train and val.\n- **FIX:** `LassoCV` is replaced with `LogisticRegression(penalty='l1')` — LASSO is for regression, not classification.\n- **FIX:** Median imputation for NaN uses training-data medians only.","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION G: ANTI-LEAKAGE FEATURE SELECTION PIPELINE\n# ─────────────────────────────────────────────────────────────\n\nMETA_COLS = ['BraTS21ID', 'MGMT_value']\nfeat_cols = [c for c in feat_train_df.columns if c not in META_COLS]\n\nX_train_raw = feat_train_df[feat_cols].copy()\ny_train     = feat_train_df['MGMT_value'].values\nX_val_raw   = feat_val_df[feat_cols].copy()\ny_val       = feat_val_df['MGMT_value'].values\n\nprint(f'Feature matrix — Train: {X_train_raw.shape} | Val: {X_val_raw.shape}')\n\n# ── NaN Audit ────────────────────────────────────────────────\nnan_counts = X_train_raw.isna().sum().sort_values(ascending=False)\nprint(f'\\nMissing values (top 10):')\nprint(nan_counts.head(10))\n\n# FIX: Compute medians ONLY from training data, apply to both\ntrain_medians = X_train_raw.median()\nX_train_raw   = X_train_raw.fillna(train_medians)\nX_val_raw     = X_val_raw.fillna(train_medians)  # Use TRAIN medians for val\n\nprint(f'\\n✅ After median imputation (train medians applied to both sets):')\nprint(f'   Train NaN count: {X_train_raw.isna().sum().sum()}')\nprint(f'   Val   NaN count: {X_val_raw.isna().sum().sum()}')\n\n# ── Step 1: Correlation Filtering ────────────────────────────\nprint('\\n━━ Step 1: Correlation Filtering ━━')\ncorr_matrix = X_train_raw.corr().abs()\nupper_tri   = corr_matrix.where(\n    np.triu(np.ones(corr_matrix.shape), k=1).astype(bool)\n)\nto_drop   = [col for col in upper_tri.columns if any(upper_tri[col] > 0.90)]\nX_train_d = X_train_raw.drop(columns=to_drop)\nX_val_d   = X_val_raw.drop(columns=to_drop)\nfeat_names_d = X_train_d.columns.tolist()\nprint(f'  Original: {X_train_raw.shape[1]} → Dropped: {len(to_drop)} → Remaining: {len(feat_names_d)}')\n\n# ── Step 2: Standardize (FIT on train ONLY) ──────────────────\nscaler        = StandardScaler()\nX_train_sc    = scaler.fit_transform(X_train_d)   # fit + transform on train\nX_val_sc      = scaler.transform(X_val_d)          # transform only on val\n\n# ── Step 3: ANOVA F-test (train data only) ───────────────────\nprint('\\n━━ Step 2: ANOVA F-test ━━')\nanova_sel = SelectKBest(score_func=f_classif, k='all')\nanova_sel.fit(X_train_sc, y_train)\nanova_df = pd.DataFrame({\n    'feature': feat_names_d,\n    'f_score': anova_sel.scores_,\n    'p_value': anova_sel.pvalues_\n}).sort_values('f_score', ascending=False)\nsig_anova = anova_df[anova_df['p_value'] < 0.05]\nprint(f'  Significant (p<0.05): {len(sig_anova)}')\nprint(anova_df.head(10).to_string(index=False))\n\n# ── Step 4: Mutual Information (train only) ──────────────────\nprint('\\n━━ Step 3: Mutual Information ━━')\nmi_scores = mutual_info_classif(X_train_sc, y_train, random_state=SEED)\nmi_df = pd.DataFrame({\n    'feature':  feat_names_d,\n    'mi_score': mi_scores\n}).sort_values('mi_score', ascending=False)\nprint(mi_df.head(10).to_string(index=False))\n\n# ── Step 5: L1-Logistic Regression (FIX: not LASSO regression) ──\nprint('\\n━━ Step 4: L1 Logistic Regression (sparse feature selection) ━━')\n# FIX: LassoCV is a regression model. For classification, use LogisticRegression(penalty='l1')\nl1_lr = LogisticRegression(penalty='l1', solver='liblinear', C=0.1,\n                            max_iter=2000, random_state=SEED)\nl1_lr.fit(X_train_sc, y_train)\nl1_coefs = pd.Series(np.abs(l1_lr.coef_[0]), index=feat_names_d)\nl1_selected = l1_coefs[l1_coefs > 0].sort_values(ascending=False)\nprint(f'  Selected by L1 (coef≠0): {len(l1_selected)}')\nprint(l1_selected.head(15))\n\n# ── Step 6: Random Forest Importance (train only) ────────────\nprint('\\n━━ Step 5: Random Forest Feature Importance ━━')\nrf_fs = RandomForestClassifier(n_estimators=200, random_state=SEED, n_jobs=-1, max_depth=4)\nrf_fs.fit(X_train_sc, y_train)\nrf_importances = pd.Series(rf_fs.feature_importances_, index=feat_names_d)\nrf_top20 = rf_importances.sort_values(ascending=False).head(20)\nprint(rf_top20)\n\n# ── Consensus Feature Set ─────────────────────────────────────\nanova_top  = set(anova_df.head(TOP_K)['feature'])\nmi_top     = set(mi_df.head(TOP_K)['feature'])\nrf_top_set = set(rf_top20.index)\nl1_set     = set(l1_selected.index)\n\nvote_counter = Counter()\nfor feat_set in [anova_top, mi_top, rf_top_set, l1_set]:\n    for f in feat_set:\n        vote_counter[f] += 1\n\nconsensus_features = [f for f, v in vote_counter.items() if v >= 2]\nif not consensus_features:\n    print('⚠️  No consensus features (threshold 2). Relaxing to top-10 ANOVA.')\n    consensus_features = anova_df.head(10)['feature'].tolist()\n\nprint(f'\\n📌 Consensus features (≥2 methods): {len(consensus_features)}')\nprint(consensus_features)\n\n# FIX: Scaler already fitted above — just index to selected columns\n# We need to re-extract columns from the ALREADY scaled data, not re-fit scaler\nconsensus_idx = [feat_names_d.index(f) for f in consensus_features]\nX_train_sel   = X_train_sc[:, consensus_idx]\nX_val_sel     = X_val_sc[:, consensus_idx]\n\nprint(f'\\n✅ Feature selection complete (NO leakage from validation data)')\nprint(f'   Train: {X_train_sel.shape} | Val: {X_val_sel.shape}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T12:19:46.479628Z","iopub.execute_input":"2026-03-29T12:19:46.484713Z","iopub.status.idle":"2026-03-29T12:19:47.447392Z","shell.execute_reply.started":"2026-03-29T12:19:46.484543Z","shell.execute_reply":"2026-03-29T12:19:47.445608Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Feature Selection Summary Plot ────────────────────────────\nfig, axes = plt.subplots(1, 3, figsize=(20, 6))\nfig.suptitle('Feature Selection — Multi-Method Ranking (Training Data Only)',\n             fontsize=14, fontweight='bold', color='#e0e0ff')\n\ntop15_anova = anova_df.head(15)\naxes[0].barh(top15_anova['feature'][::-1], top15_anova['f_score'][::-1], color=PALETTE[0])\naxes[0].set_title('ANOVA F-score (Top 15)', color='#e0e0ff')\naxes[0].set_xlabel('F-score')\naxes[0].grid(axis='x', alpha=0.3)\n\ntop15_mi = mi_df.head(15)\naxes[1].barh(top15_mi['feature'][::-1], top15_mi['mi_score'][::-1], color=PALETTE[1])\naxes[1].set_title('Mutual Information (Top 15)', color='#e0e0ff')\naxes[1].set_xlabel('MI Score')\naxes[1].grid(axis='x', alpha=0.3)\n\nrf_top15 = rf_importances.sort_values(ascending=False).head(15)\naxes[2].barh(rf_top15.index[::-1], rf_top15.values[::-1], color=PALETTE[2])\naxes[2].set_title('RF Importance (Top 15)', color='#e0e0ff')\naxes[2].set_xlabel('Importance')\naxes[2].grid(axis='x', alpha=0.3)\n\nfor ax in axes:\n    ax.tick_params(labelsize=8)\n\nplt.tight_layout()\nplt.savefig('feature_selection.png', dpi=150, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T12:19:55.874367Z","iopub.execute_input":"2026-03-29T12:19:55.875276Z","iopub.status.idle":"2026-03-29T12:19:57.728507Z","shell.execute_reply.started":"2026-03-29T12:19:55.875220Z","shell.execute_reply":"2026-03-29T12:19:57.727162Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## H. 🧬 Radiogenomic Association Analysis\n\nThis section is the scientific core of radiogenomics: **statistically linking imaging features to MGMT methylation status**.\n\nWe use:\n- **Mann-Whitney U test** (non-parametric, robust to non-normality)\n- **Rank-biserial correlation** as effect size\n- **Benjamini-Hochberg FDR correction** to control false discovery rate (critical for high-dimensional testing)\n\nNote: Association analysis is conducted on training data only.","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION H.1: STATISTICAL ASSOCIATION ANALYSIS\n# Conducted on TRAINING data only (no leakage)\n# ─────────────────────────────────────────────────────────────\n\nprint('═══ RADIOGENOMIC ASSOCIATION ANALYSIS (Training Data) ═══')\n\nX_train_df = pd.DataFrame(X_train_raw, columns=feat_names_d) if hasattr(X_train_raw, 'columns') \\\n    else pd.DataFrame(scaler.inverse_transform(X_train_sc), columns=feat_names_d)\n\n# Use the raw (unscaled) decorrelated features for interpretable mean comparison\nX_train_df = X_train_d.copy()\nX_train_df['MGMT_value'] = y_train\n\nresults = []\nfor feat in feat_names_d:\n    vals   = X_train_df[feat].values\n    group0 = vals[y_train == 0]\n    group1 = vals[y_train == 1]\n    group0 = group0[~np.isnan(group0)]\n    group1 = group1[~np.isnan(group1)]\n    if len(group0) < 3 or len(group1) < 3:\n        continue\n    stat, pval = stats.mannwhitneyu(group0, group1, alternative='two-sided')\n    n1, n2     = len(group0), len(group1)\n    effect     = 1 - (2 * stat) / (n1 * n2)  # rank-biserial correlation\n    results.append({\n        'feature':     feat,\n        'u_statistic': stat,\n        'p_value':     pval,\n        'effect_size': abs(effect),\n        'mean_unmeth': float(group0.mean()),\n        'mean_meth':   float(group1.mean()),\n        'delta_mean':  float(group1.mean() - group0.mean()),\n    })\n\nassoc_df = pd.DataFrame(results).sort_values('p_value')\n\n# FDR correction (Benjamini-Hochberg)\n_, fdr_pvals, _, _ = multipletests(assoc_df['p_value'].values, method='fdr_bh')\nassoc_df['fdr_pval']    = fdr_pvals\nassoc_df['significant'] = assoc_df['fdr_pval'] < 0.05\n\nprint('Top 15 features associated with MGMT methylation:')\ndisp_cols = ['feature','p_value','fdr_pval','effect_size','mean_unmeth','mean_meth','significant']\nprint(assoc_df[disp_cols].head(15).to_string(index=False, float_format='{:.4f}'.format))\nprint(f'\\n🔬 Significant after FDR (q<0.05): {assoc_df[\"significant\"].sum()} features')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T12:20:37.992925Z","iopub.execute_input":"2026-03-29T12:20:37.993468Z","iopub.status.idle":"2026-03-29T12:20:38.083430Z","shell.execute_reply.started":"2026-03-29T12:20:37.993431Z","shell.execute_reply":"2026-03-29T12:20:38.082390Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Volcano Plot: Radiogenomic Association ────────────────────\n\nfig, axes = plt.subplots(1, 2, figsize=(18, 6))\nfig.suptitle('Radiogenomic Association: Imaging Features vs MGMT Methylation (Training Data)',\n             fontsize=13, fontweight='bold', color='#e0e0ff')\n\n# Volcano Plot\nax = axes[0]\nlog_p    = -np.log10(assoc_df['p_value'].clip(1e-10))\neffect   = assoc_df['effect_size']\ncol_pts  = ['#ff6b6b' if sig else '#444466' for sig in assoc_df['significant']]\n\nax.scatter(effect, log_p, c=col_pts, alpha=0.8, s=40, edgecolors='none')\nax.axhline(-np.log10(0.05), color='#ffa726', linestyle='--', alpha=0.7, label='p=0.05')\nax.set_xlabel('Effect Size (rank-biserial |r|)')\nax.set_ylabel('-log10(p-value)')\nax.set_title('Volcano Plot — Feature–MGMT Associations', color='#e0e0ff')\nax.legend()\nax.grid(True, alpha=0.3)\n\nfor _, row in assoc_df.head(5).iterrows():\n    xi = row['effect_size']\n    yi = -np.log10(max(row['p_value'], 1e-10))\n    ax.annotate(row['feature'].replace('_','\\n'), (xi, yi), fontsize=7, color='#ffd700',\n                xytext=(xi+0.02, yi+0.1), arrowprops=dict(arrowstyle='->', color='#888'))\n\n# Box plots: top 5 features by p-value\nax2 = axes[1]\ntop5 = assoc_df.head(5)['feature'].tolist()\ndata_meth   = [X_train_d.loc[y_train == 1, f].dropna().values for f in top5]\ndata_unmeth = [X_train_d.loc[y_train == 0, f].dropna().values for f in top5]\n\npos0 = np.arange(len(top5)) * 2\npos1 = pos0 + 0.7\n\nbp0 = ax2.boxplot(data_unmeth, positions=pos0, widths=0.6, patch_artist=True,\n                  boxprops=dict(facecolor='#7b68ee', alpha=0.7),\n                  medianprops=dict(color='white', linewidth=2),\n                  whiskerprops=dict(color='#7b68ee'), capprops=dict(color='#7b68ee'),\n                  flierprops=dict(marker='o', color='#7b68ee', alpha=0.3, markersize=3))\nbp1 = ax2.boxplot(data_meth, positions=pos1, widths=0.6, patch_artist=True,\n                  boxprops=dict(facecolor='#ff6b6b', alpha=0.7),\n                  medianprops=dict(color='white', linewidth=2),\n                  whiskerprops=dict(color='#ff6b6b'), capprops=dict(color='#ff6b6b'),\n                  flierprops=dict(marker='o', color='#ff6b6b', alpha=0.3, markersize=3))\n\nax2.set_xticks(pos0 + 0.35)\nax2.set_xticklabels([f.replace('_','\\n') for f in top5], fontsize=8)\nax2.set_title('Top 5 Features by Group', color='#e0e0ff')\nax2.legend([bp0['boxes'][0], bp1['boxes'][0]], ['Unmethylated', 'Methylated'],\n           loc='upper right', fontsize=9)\nax2.grid(axis='y', alpha=0.3)\n\nplt.tight_layout()\nplt.savefig('radiogenomic_association.png', dpi=150, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T12:20:45.515728Z","iopub.execute_input":"2026-03-29T12:20:45.516712Z","iopub.status.idle":"2026-03-29T12:20:47.081174Z","shell.execute_reply.started":"2026-03-29T12:20:45.516669Z","shell.execute_reply":"2026-03-29T12:20:47.080051Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION H.2: MULTI-MODALITY FUSION ANALYSIS\n# Evaluate each modality in isolation vs combined — using CV on training data\n# ─────────────────────────────────────────────────────────────\n\nprint('═══ MULTI-MODALITY FUSION ANALYSIS ═══')\n\ncv = StratifiedKFold(n_splits=CV_FOLDS, shuffle=True, random_state=SEED)\n\ndef evaluate_modality_subset(modality_list, X_df, y, cv, label=''):\n    mod_feats = [c for c in X_df.columns if any(c.startswith(m+'_') for m in modality_list)]\n    if not mod_feats:\n        return {'label': label, 'auc': np.nan, 'auc_std': np.nan}\n    Xm = X_df[mod_feats].fillna(X_df[mod_feats].median()).values\n    pipe = Pipeline([\n        ('scaler', StandardScaler()),\n        ('clf', LogisticRegression(max_iter=1000, random_state=SEED, C=0.1))\n    ])\n    sc = cross_validate(pipe, Xm, y, cv=cv, scoring=['roc_auc','accuracy'], n_jobs=-1)\n    return {\n        'label':   label,\n        'auc':     sc['test_roc_auc'].mean(),\n        'auc_std': sc['test_roc_auc'].std(),\n        'acc':     sc['test_accuracy'].mean(),\n        'n_feats': len(mod_feats)\n    }\n\nfusion_results = []\nfor mod in MODALITIES:\n    res = evaluate_modality_subset([mod], X_train_d, y_train, cv, f'{mod} only')\n    fusion_results.append(res)\n    print(f'  {mod:8s}: AUC={res[\"auc\"]:.3f}±{res[\"auc_std\"]:.3f}')\n\nres_all = evaluate_modality_subset(MODALITIES, X_train_d, y_train, cv, 'All Modalities')\nfusion_results.append(res_all)\nprint(f'  {\"All fused\":8s}: AUC={res_all[\"auc\"]:.3f}±{res_all[\"auc_std\"]:.3f}')\n\nfusion_df = pd.DataFrame(fusion_results)\n\nfig, ax = plt.subplots(figsize=(12, 5))\nx_pos = np.arange(len(fusion_df))\nbar_colors = PALETTE[:4] + ['#ffd700']\nbars = ax.bar(x_pos, fusion_df['auc'], color=bar_colors[:len(fusion_df)], alpha=0.85,\n              yerr=fusion_df['auc_std'], capsize=5, edgecolor='white', linewidth=0.5)\nax.set_xticks(x_pos)\nax.set_xticklabels(fusion_df['label'], fontsize=11)\nax.set_ylabel('ROC-AUC (CV mean ± std)')\nax.set_title('Multi-Modality Fusion Analysis', color='#e0e0ff', fontweight='bold')\nax.set_ylim(0.3, 1.0)\nax.axhline(0.5, color='#888', linestyle='--', alpha=0.6, label='Random baseline')\nfor bar, val, std in zip(bars, fusion_df['auc'], fusion_df['auc_std']):\n    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + std + 0.01,\n            f'{val:.3f}', ha='center', fontsize=10, fontweight='bold', color='white')\nax.legend()\nax.grid(axis='y', alpha=0.3)\nplt.tight_layout()\nplt.savefig('modality_fusion.png', dpi=150, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T12:20:55.590325Z","iopub.execute_input":"2026-03-29T12:20:55.591342Z","iopub.status.idle":"2026-03-29T12:21:00.123921Z","shell.execute_reply.started":"2026-03-29T12:20:55.591300Z","shell.execute_reply":"2026-03-29T12:21:00.122659Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## I. 🤖 Predictive Modeling with Stratified K-Fold CV\n\n**Fixes applied:**\n- All models evaluated within CV folds on training data (no leakage from validation set during CV)\n- Final held-out validation set used ONCE at the end for unbiased performance estimate\n- `Pipeline` used so scaler is re-fitted inside each fold (true anti-leakage CV)","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION I: MODELING — STRATIFIED K-FOLD CV\n# FIX: Use Pipeline so scaler is fit inside each fold (not before CV)\n# FIX: Evaluate on training CV folds; hold-out val used ONCE at the end\n# ─────────────────────────────────────────────────────────────\n\ncv = StratifiedKFold(n_splits=CV_FOLDS, shuffle=True, random_state=SEED)\n\nmodels = {\n    'Dummy (baseline)':   DummyClassifier(strategy='stratified', random_state=SEED),\n    'Logistic Regression': LogisticRegression(max_iter=2000, C=0.1, random_state=SEED),\n    'Random Forest':       RandomForestClassifier(n_estimators=300, max_depth=5,\n                                                  random_state=SEED, n_jobs=-1),\n}\n\n# Feature configs: all decorrelated vs consensus selected\n# NOTE: X_train_sc and X_train_sel are already scaled from the train-only scaler.\n# For proper CV, we SHOULD use raw features inside a pipeline. However,\n# to allow fair comparison of feature sets without double-fitting,\n# we wrap in a Pipeline with an identity step and use pre-scaled data.\n# ALTERNATIVELY (more correct for CV): use Pipeline(scaler + clf) on X_train_d.\n\nfeature_configs = {\n    'All Features':      (X_train_d, feat_names_d),      # raw train data\n    'Selected Features': (X_train_d[consensus_features], consensus_features),\n}\n\nresults_table = []\nprint('═══ MODEL EVALUATION (Stratified K-Fold CV on Training Data) ═══')\nprint(f'CV: {CV_FOLDS}-fold stratified | All scalers fit inside each fold\\n')\n\nfor feat_name, (X_feat_raw, feat_list) in feature_configs.items():\n    print(f'\\n── Features: {feat_name} ({X_feat_raw.shape[1]} features) ──')\n    for model_name, clf in models.items():\n        # Pipeline: scaler fit inside each fold (true anti-leakage CV)\n        pipe = Pipeline([\n            ('scaler', StandardScaler()),\n            ('clf', clf)\n        ])\n        cv_res = cross_validate(\n            pipe, X_feat_raw.values, y_train, cv=cv,\n            scoring=['roc_auc', 'accuracy', 'f1'], n_jobs=-1\n        )\n        auc     = cv_res['test_roc_auc'].mean()\n        auc_std = cv_res['test_roc_auc'].std()\n        acc     = cv_res['test_accuracy'].mean()\n        f1      = cv_res['test_f1'].mean()\n        results_table.append({\n            'Features': feat_name, 'Model': model_name,\n            'AUC': auc, 'AUC_std': auc_std, 'Accuracy': acc, 'F1': f1\n        })\n        print(f'  {model_name:30s}: AUC={auc:.3f}±{auc_std:.3f}  Acc={acc:.3f}  F1={f1:.3f}')\n\nresults_df = pd.DataFrame(results_table)\nprint('\\n✅ CV Modeling complete')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T12:25:50.594052Z","iopub.execute_input":"2026-03-29T12:25:50.594836Z","iopub.status.idle":"2026-03-29T12:25:55.235598Z","shell.execute_reply.started":"2026-03-29T12:25:50.594795Z","shell.execute_reply":"2026-03-29T12:25:55.234168Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Performance Comparison Plot ───────────────────────────────\n\nfig, axes = plt.subplots(1, 2, figsize=(16, 6))\nfig.suptitle('Model Comparison: Before vs After Feature Selection (CV on Training Data)',\n             fontsize=14, fontweight='bold', color='#e0e0ff')\n\nfor ax_i, metric in enumerate(['AUC', 'F1']):\n    ax = axes[ax_i]\n    pivot = results_df.pivot(index='Model', columns='Features', values=metric)\n    x   = np.arange(len(pivot.index))\n    w   = 0.35\n    cols = pivot.columns.tolist()\n    for j, col in enumerate(cols):\n        bars = ax.bar(x + (j-0.5)*w, pivot[col], w, label=col,\n                     color=PALETTE[j], alpha=0.85, edgecolor='white', linewidth=0.5)\n        if metric == 'AUC':\n            stds = results_df[results_df['Features']==col].set_index('Model')['AUC_std']\n            ax.errorbar(x + (j-0.5)*w, pivot[col],\n                        yerr=[stds.get(m, 0) for m in pivot.index],\n                        fmt='none', ecolor='white', capsize=4, alpha=0.7)\n    ax.set_xticks(x)\n    ax.set_xticklabels(pivot.index, rotation=15, ha='right', fontsize=9)\n    ax.set_ylabel(metric)\n    ax.set_title(metric, color='#e0e0ff')\n    ax.set_ylim(0.3, 1.05)\n    ax.axhline(0.5, color='#888', linestyle='--', alpha=0.4)\n    ax.legend(fontsize=8)\n    ax.grid(axis='y', alpha=0.3)\n\nplt.tight_layout()\nplt.savefig('model_comparison.png', dpi=150, bbox_inches='tight')\nplt.show()\n\nprint('\\n📋 Summary Table:')\nprint(results_df.to_string(index=False, float_format='{:.3f}'.format))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T12:26:16.007788Z","iopub.execute_input":"2026-03-29T12:26:16.008186Z","iopub.status.idle":"2026-03-29T12:26:16.851336Z","shell.execute_reply.started":"2026-03-29T12:26:16.008146Z","shell.execute_reply":"2026-03-29T12:26:16.850218Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## J. 📈 Final Evaluation on Held-Out Validation Set\n\n**This is the unbiased estimate of generalization performance.**  \nThe validation set has NOT been used in any prior step (no normalization fitting, no feature selection, no model selection).\n\nWe also compute **bootstrapped 95% confidence intervals** on ROC-AUC — required for publication-quality reporting.","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION J: FINAL EVALUATION ON HELD-OUT VALIDATION SET\n# FIX: Validation set used only HERE, not in any previous step.\n# FIX: Bootstrapped confidence intervals added.\n# FIX: Confusion matrix from val predictions (not full data).\n# ─────────────────────────────────────────────────────────────\n\n# Best model selected from CV results (lowest leakage risk: use CV AUC)\nbest_row  = results_df.loc[results_df['AUC'].idxmax()]\nprint(f'Best CV configuration:')\nprint(f'  Model:    {best_row[\"Model\"]}')\nprint(f'  Features: {best_row[\"Features\"]}')\nprint(f'  CV AUC:   {best_row[\"AUC\"]:.3f} ± {best_row[\"AUC_std\"]:.3f}')\n\n# Determine final feature set\nfinal_feat_cols = consensus_features if best_row['Features'] == 'Selected Features' else feat_names_d\nX_train_final   = X_train_d[final_feat_cols].values\nX_val_final     = X_val_d[final_feat_cols].values\n\n# Build and train final model (on all training data)\nfinal_pipe = Pipeline([\n    ('scaler', StandardScaler()),\n    ('clf', RandomForestClassifier(n_estimators=300, max_depth=5, random_state=SEED, n_jobs=-1))\n])\nfinal_pipe.fit(X_train_final, y_train)\n\n# Evaluate on held-out validation set\ny_val_prob = final_pipe.predict_proba(X_val_final)[:, 1]\ny_val_pred = final_pipe.predict(X_val_final)\n\nval_auc = roc_auc_score(y_val, y_val_prob)\nval_f1  = f1_score(y_val, y_val_pred)\nval_acc = accuracy_score(y_val, y_val_pred)\n\nprint(f'\\n=== VALIDATION SET RESULTS ===')\nprint(f'  ROC-AUC:  {val_auc:.3f}')\nprint(f'  F1-Score: {val_f1:.3f}')\nprint(f'  Accuracy: {val_acc:.3f}')\n\n# ── Bootstrapped 95% CI on AUC ────────────────────────────────\n# Required for publication-quality reporting\nN_BOOTSTRAP = 1000\nboot_aucs   = []\nrng = np.random.RandomState(SEED)\nn_val = len(y_val)\nfor _ in range(N_BOOTSTRAP):\n    idx = rng.randint(0, n_val, size=n_val)\n    if len(np.unique(y_val[idx])) < 2:\n        continue\n    boot_aucs.append(roc_auc_score(y_val[idx], y_val_prob[idx]))\n\nci_low  = np.percentile(boot_aucs, 2.5)\nci_high = np.percentile(boot_aucs, 97.5)\nprint(f'\\n  Bootstrapped 95% CI on AUC: [{ci_low:.3f}, {ci_high:.3f}]')\nprint(f'  (n_bootstrap={N_BOOTSTRAP})')\nprint('\\n📋 Classification Report (Validation Set):')\nprint(classification_report(y_val, y_val_pred, target_names=['MGMT− (0)', 'MGMT+ (1)']))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T12:26:31.338537Z","iopub.execute_input":"2026-03-29T12:26:31.339747Z","iopub.status.idle":"2026-03-29T12:26:33.588990Z","shell.execute_reply.started":"2026-03-29T12:26:31.339708Z","shell.execute_reply":"2026-03-29T12:26:33.587963Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Evaluation Plots (Validation Set) ────────────────────────\n\nfig, axes = plt.subplots(1, 3, figsize=(20, 6))\nfig.suptitle('Model Evaluation — Random Forest (Held-Out Validation Set)',\n             fontsize=14, fontweight='bold', color='#e0e0ff')\n\n# ROC Curve (val set)\nfpr, tpr, _ = roc_curve(y_val, y_val_prob)\naxes[0].plot(fpr, tpr, color=PALETTE[0], linewidth=2.5,\n             label=f'Val AUC={val_auc:.3f} [{ci_low:.3f}, {ci_high:.3f}]')\naxes[0].plot([0,1],[0,1], 'r--', alpha=0.5, label='Random')\naxes[0].fill_between(fpr, tpr, alpha=0.1, color=PALETTE[0])\naxes[0].set_xlabel('False Positive Rate')\naxes[0].set_ylabel('True Positive Rate')\naxes[0].set_title('ROC Curve (Validation Set)\\nwith 95% Bootstrap CI', color='#e0e0ff')\naxes[0].legend(fontsize=9)\naxes[0].grid(True, alpha=0.3)\n\n# Confusion Matrix (val set — FIX: not full data)\ncm = confusion_matrix(y_val, y_val_pred)\ncm_disp = ConfusionMatrixDisplay(cm, display_labels=['MGMT−', 'MGMT+'])\ncm_disp.plot(ax=axes[1], colorbar=False, cmap='Blues')\naxes[1].set_title('Confusion Matrix (Validation Set)', color='#e0e0ff')\naxes[1].set_facecolor('#1a1a2e')\n\n# Bootstrap AUC distribution\naxes[2].hist(boot_aucs, bins=40, color=PALETTE[2], alpha=0.85, edgecolor='white')\naxes[2].axvline(val_auc, color='white', linewidth=2, linestyle='-', label=f'Observed AUC={val_auc:.3f}')\naxes[2].axvline(ci_low,  color='#ff6b6b', linewidth=1.5, linestyle='--', label=f'95% CI lower={ci_low:.3f}')\naxes[2].axvline(ci_high, color='#ff6b6b', linewidth=1.5, linestyle='--', label=f'95% CI upper={ci_high:.3f}')\naxes[2].set_xlabel('Bootstrap AUC')\naxes[2].set_ylabel('Frequency')\naxes[2].set_title('Bootstrap AUC Distribution (n=1000)', color='#e0e0ff')\naxes[2].legend(fontsize=8)\naxes[2].grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.savefig('evaluation.png', dpi=150, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T12:26:41.666131Z","iopub.execute_input":"2026-03-29T12:26:41.666450Z","iopub.status.idle":"2026-03-29T12:26:42.973736Z","shell.execute_reply.started":"2026-03-29T12:26:41.666423Z","shell.execute_reply":"2026-03-29T12:26:42.972507Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## K. 🔍 Radiogenomic Interpretation Dashboard","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION K: INTERPRETATION & BIOLOGICAL RELEVANCE\n# ─────────────────────────────────────────────────────────────\n\n# Feature importance from the final model\nfinal_rf        = final_pipe.named_steps['clf']\nrf_imp_final    = pd.Series(final_rf.feature_importances_, index=final_feat_cols)\nrf_imp_sorted   = rf_imp_final.sort_values(ascending=True)\n\nfig = plt.figure(figsize=(18, 12))\ngs  = gridspec.GridSpec(2, 2, figure=fig, hspace=0.4, wspace=0.35)\nfig.suptitle('Radiogenomics Interpretation Dashboard',\n             fontsize=15, fontweight='bold', color='#e0e0ff', y=1.01)\n\n# Panel 1: Feature Importance (Final Model)\nax1 = fig.add_subplot(gs[0, 0])\ncolors_bar = []\nfor f in rf_imp_sorted.index:\n    if f.startswith('FLAIR'):  colors_bar.append(PALETTE[0])\n    elif f.startswith('T1wCE'):colors_bar.append(PALETTE[2])\n    elif f.startswith('T1w'):  colors_bar.append(PALETTE[1])\n    else:                      colors_bar.append(PALETTE[3])\nax1.barh(range(len(rf_imp_sorted)), rf_imp_sorted.values,\n         color=colors_bar, alpha=0.85, edgecolor='none')\nax1.set_yticks(range(len(rf_imp_sorted)))\nax1.set_yticklabels([f.replace('_glcm_','\\nGLCM_').replace('_','\\n')\n                     for f in rf_imp_sorted.index], fontsize=7)\nax1.set_xlabel('Feature Importance')\nax1.set_title('RF Feature Importance (Final Model)', color='#e0e0ff', fontsize=10)\nlegend_patches = [plt.Rectangle((0,0),1,1,fc=PALETTE[i],label=m)\n                  for i,m in enumerate(MODALITIES)]\nax1.legend(handles=legend_patches, loc='lower right', fontsize=7)\nax1.grid(axis='x', alpha=0.3)\n\n# Panel 2: Correlation heatmap (selected features)\nax2 = fig.add_subplot(gs[0, 1])\nfeat_subset = final_feat_cols[:min(12, len(final_feat_cols))]\ncorr_sub    = X_train_d[feat_subset].corr()\nsns.heatmap(corr_sub, ax=ax2, cmap='coolwarm', center=0,\n            annot=True, fmt='.2f', annot_kws={'size': 7},\n            linewidths=0.5, cbar_kws={'shrink': 0.8})\nax2.set_title('Feature Correlation (Training Data)', color='#e0e0ff', fontsize=10)\nax2.tick_params(labelsize=7)\n\n# Panel 3: Effect size vs Feature Importance scatter\nax3 = fig.add_subplot(gs[1, 0])\nshared_feats = list(set(rf_imp_final.index) & set(assoc_df['feature']))\nif shared_feats:\n    assoc_sub = assoc_df.set_index('feature').loc[shared_feats]\n    rf_sub    = rf_imp_final.loc[shared_feats]\n    sc = ax3.scatter(rf_sub.values, assoc_sub['effect_size'].values,\n                     c=assoc_sub['p_value'].apply(lambda p: -np.log10(p+1e-10)),\n                     cmap='hot', s=60, alpha=0.85, edgecolors='none')\n    plt.colorbar(sc, ax=ax3, label='-log10(p-value)')\n    ax3.set_xlabel('RF Importance')\n    ax3.set_ylabel('Effect Size (rank-biserial |r|)')\n    ax3.set_title('Predictive Power vs Statistical Effect',\n                  color='#e0e0ff', fontsize=10)\n    for feat in assoc_sub.nsmallest(3, 'p_value').index:\n        ax3.annotate(feat.replace('_','\\n'), (rf_imp_final[feat], assoc_sub.loc[feat,'effect_size']),\n                     fontsize=7, color='#ffd700')\n    ax3.grid(True, alpha=0.3)\n\n# Panel 4: Bootstrap AUC distribution (from validation)\nax4 = fig.add_subplot(gs[1, 1])\nax4.hist(boot_aucs, bins=40, color=PALETTE[4], alpha=0.85, edgecolor='white')\nax4.axvline(val_auc,  color='white',   linewidth=2, label=f'AUC={val_auc:.3f}')\nax4.axvline(ci_low,   color='#ff6b6b', linewidth=1.5, linestyle='--', label=f'95% CI')\nax4.axvline(ci_high,  color='#ff6b6b', linewidth=1.5, linestyle='--')\nax4.axvline(0.5,      color='#888',    linewidth=1, linestyle=':',  label='Random')\nax4.set_xlabel('Bootstrap AUC')\nax4.set_title('Validation AUC with 95% Bootstrap CI', color='#e0e0ff', fontsize=10)\nax4.legend(fontsize=8)\nax4.grid(True, alpha=0.3)\n\nplt.savefig('interpretation_dashboard.png', dpi=150, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T12:26:47.115237Z","iopub.execute_input":"2026-03-29T12:26:47.115575Z","iopub.status.idle":"2026-03-29T12:26:50.572081Z","shell.execute_reply.started":"2026-03-29T12:26:47.115547Z","shell.execute_reply":"2026-03-29T12:26:50.570284Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION K.2: BIOLOGICAL INTERPRETATION SUMMARY\n# ─────────────────────────────────────────────────────────────\n\nBIOLOGICAL_ANNOTATIONS = {\n    'mean':               'Average signal intensity → tissue cellularity, edema volume',\n    'entropy':            'Information heterogeneity → intratumoral heterogeneity (ITH), key MGMT correlate',\n    'skewness':           'Signal asymmetry → necrotic vs viable tumor fraction',\n    'kurtosis':           'Tail behavior → outlier enhancement zones, micro-necrosis',\n    'std':                'Signal variance → heterogeneous microenvironment',\n    'glcm_contrast':      'Textural contrast → architectural irregularity, cellular density variation',\n    'glcm_correlation':   'Structural regularity → tissue organization (lower in aggressive GBM)',\n    'glcm_homogeneity':   'Textural smoothness → uniform tissue regions',\n    'glcm_energy':        'Textural uniformity → homogeneous tumor sub-regions',\n    'glcm_dissimilarity': 'Local heterogeneity → boundary irregularity',\n    'area_frac':          'Proportion of hyperintense region → tumor/edema burden',\n    'solidity':           'Convexity → tumor margin regularity (irregular in MGMT-unmethylated)',\n    'eccentricity':       'Shape elongation → infiltrative growth pattern',\n    'compactness':        'Shape circularity → degree of tumor boundary regularity',\n    'iqr':                'Interquartile intensity range → robust heterogeneity measure',\n}\n\nprint('═══════════════════════════════════════════════════════════')\nprint('  RADIOGENOMICS BIOLOGICAL INTERPRETATION')\nprint('═══════════════════════════════════════════════════════════\\n')\n\nprint('Top features and their biological relevance:\\n')\nfor feat in assoc_df.head(10)['feature']:\n    for key, bio in BIOLOGICAL_ANNOTATIONS.items():\n        if feat.endswith(key):\n            mod  = feat.split('_')[0]\n            pval = assoc_df.loc[assoc_df['feature']==feat, 'p_value'].values[0]\n            eff  = assoc_df.loc[assoc_df['feature']==feat, 'effect_size'].values[0]\n            print(f'🔬 [{mod}] {feat}')\n            print(f'   Biological: {bio}')\n            print(f'   p={pval:.4f} | effect_size={eff:.3f}\\n')\n            break\n\nprint(\"\"\"\n📖 LITERATURE CONTEXT:\n───────────────────────────────────────────────────────────\n1. ENTROPY / HETEROGENEITY:\n   MGMT-methylated GBMs exhibit greater FLAIR signal heterogeneity,\n   reflecting mixtures of viable tumor, necrosis, and edema.\n   Entropy features are among the most discriminative radiomics\n   biomarkers for MGMT prediction (Kickingereder et al., 2016;\n   Lohmann et al., 2021).\n\n2. TEXTURE (GLCM CONTRAST, DISSIMILARITY):\n   MGMT-unmethylated tumors show more homogeneous enhancement (T1wCE)\n   due to uniform BBB disruption. Higher GLCM contrast in MGMT+\n   cases reflects heterogeneous vascular architecture.\n\n3. SHAPE (ECCENTRICITY, SOLIDITY):\n   MGMT-unmethylated GBMs are often more infiltrative, leading to\n   irregular tumor shapes (lower solidity, higher eccentricity).\n   MGMT-methylated tumors may exhibit more circumscribed morphology.\n\n4. FLAIR/T2w SIGNAL (MEAN, IQR):\n   MGMT methylation influences tumor water content and cellularity,\n   reflected in T2/FLAIR signal levels. Higher FLAIR mean often\n   correlates with greater peritumoral edema extent.\n\n5. MULTI-MODALITY ADVANTAGE:\n   Each MRI sequence captures complementary biological processes.\n   FLAIR: edema; T1wCE: BBB disruption; T2w: infiltration.\n   Multi-modal fusion consistently outperforms single-modality\n   approaches in MGMT prediction tasks.\n\"\"\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T12:26:57.562379Z","iopub.execute_input":"2026-03-29T12:26:57.563827Z","iopub.status.idle":"2026-03-29T12:26:57.585148Z","shell.execute_reply.started":"2026-03-29T12:26:57.563783Z","shell.execute_reply":"2026-03-29T12:26:57.584313Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## L. 🏆 Final Summary & Thesis Recommendations","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION L: FINAL SUMMARY\n# ─────────────────────────────────────────────────────────────\n\nprint('\\n🏆 FINAL PERFORMANCE SUMMARY')\nprint('═' * 72)\nprint('Cross-Validation Results (Training Data):')\nsummary = results_df.copy()\nsummary['AUC (mean±std)'] = summary.apply(\n    lambda r: f\"{r['AUC']:.3f} ± {r['AUC_std']:.3f}\", axis=1\n)\nprint(summary[['Features','Model','AUC (mean±std)','Accuracy','F1']]\n      .to_string(index=False, float_format='{:.3f}'.format))\n\nprint(f'\\n─── Held-Out Validation Set (UNBIASED ESTIMATE) ───')\nprint(f'  ROC-AUC:        {val_auc:.3f}')\nprint(f'  95% Bootstrap CI: [{ci_low:.3f}, {ci_high:.3f}]')\nprint(f'  F1-Score:       {val_f1:.3f}')\nprint(f'  Accuracy:       {val_acc:.3f}')\n\nprint(\"\"\"\n═══════════════════════════════════════════════════════════\n  SUGGESTIONS FOR FURTHER RESEARCH (THESIS-LEVEL)\n═══════════════════════════════════════════════════════════\n\n1. PROPER TUMOR SEGMENTATION\n   Replace intensity-based ROI with BraTS ground truth segmentation\n   masks or nnU-Net automated segmentation. This is the single largest\n   improvement for feature specificity.\n\n2. FULL IBSI-COMPLIANT FEATURE EXTRACTION\n   Use PyRadiomics (van Griethuysen et al., 2017) for standardized\n   IBSI-compliant features including GLRLM, GLSZM, NGTDM features.\n   This enables comparability with published radiomics studies.\n\n3. DEEP LEARNING FEATURES\n   Extract CNN features (transfer learning from 3D ResNet/DenseNet\n   pretrained on BraTS) as complementary deep radiomics features.\n   Combine with handcrafted features via late fusion.\n\n4. TRUE 3D VOLUMETRIC ANALYSIS\n   Replace 2D slice aggregation with full 3D feature extraction\n   (3D GLCM, volumetric shape descriptors). Requires more compute\n   but captures axial heterogeneity.\n\n5. MRI HARMONIZATION\n   Apply ComBat or WhiteStripe harmonization to remove scanner/site\n   batch effects — critical for multi-center radiomics.\n\n6. MULTIVARIATE RADIOGENOMIC MODELS\n   Integrate imaging features with clinical covariates (age, IDH status,\n   tumor grade) in a multivariate Cox regression or ensemble model.\n\n7. EXTERNAL VALIDATION\n   Validate on an independent cohort (e.g., TCGA-GBM, EGD dataset)\n   to assess model generalizability — required for publication.\n\n8. EXPLAINABILITY (XAI)\n   Apply SHAP values to explain per-patient predictions and identify\n   which imaging regions drive MGMT prediction.\n\"\"\")\n\nprint('\\n📁 Output files generated:')\nfor fname in ['class_distribution.png','preprocessing.png','roi_approximation.png',\n              'feature_selection.png','radiogenomic_association.png','modality_fusion.png',\n              'model_comparison.png','evaluation.png','interpretation_dashboard.png']:\n    print(f'   • {fname}')\n\nprint('\\n✅ Research-grade radiogenomics framework complete!')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T12:29:28.734470Z","iopub.execute_input":"2026-03-29T12:29:28.734950Z","iopub.status.idle":"2026-03-29T12:29:28.752825Z","shell.execute_reply.started":"2026-03-29T12:29:28.734905Z","shell.execute_reply":"2026-03-29T12:29:28.751511Z"}},"outputs":[],"execution_count":null}]}