Deep learning-based neuroanatomical profiling reveals population-specific brain changes in multiple sclerosis: a large-scale Middle Eastern study.

Code ↔ Paper

23 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 23 matches · 2 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted

[1] § Materials and methods › Statistical analysis framework › Primary outcome measures › Ventricular area analysis ↔ Phase5_statistical_analysis/p4_excel_analysis_developed.py, lines 403–466 · score 0.95 · morphological head mask, polynomial regression, skull mask area, brain atrophy, brain area, FSL
[2] § Materials and methods › Statistical analysis framework › Primary outcome measures › Total lesion area analysis ↔ Phase5_statistical_analysis/p4_excel_analysis_developed.py, lines 403–466 · score 0.87 · polynomial regression, Age related lesion, skull mask, age range, mask area, WMH area
[3] § Materials and methods › Deep learning architectures › Architecture selection and implementation › Trans-U-Net ↔ Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/trans_unet_model.py, the whole file · a weak match · score 0.86 · multi head attention, transformer bottleneck, Positional encoding, CNN, networks, dropout
[4] § Materials and methods › Deep learning architectures › Training protocol › Hybrid loss function strategy ↔ Phase3_model_training_and_inferencing_and_evaluation/for_GM/model_training_scripts/p1_pix2pix_var5.py, lines 622–675 · score 0.85 · unified focal loss, weighted cross entropy, weighted categorical, Class weights, refine, phase
[5] § Materials and methods › Deep learning architectures › Training protocol › Hybrid loss function strategy ↔ Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_variant_all_net.py, lines 500–553 · score 0.85 · unified focal loss, weighted cross entropy, weighted categorical, Class weights, refine, phase
[6] § Materials and methods › Statistical analysis framework › Primary outcome measures › Total lesion area analysis ↔ Phase5_statistical_analysis/p4_excel_analysis_developed.py, lines 800–879 · score 0.84 · meaningful central tendency, Shapiro Wilk, Mann Whitney, regardless, skewed, exponential
[7] § Materials and methods › Deep learning architectures › Architecture selection and implementation › DeepLabV3Plus ↔ Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/dlv3_unet_model.py, lines 12–55 · score 0.83 · atrous spatial pyramid, encoder decoder, ResNet, backbone, ASPP, convolutions
[8] § Results › Ventricular area analysis ↔ Phase5_statistical_analysis/p4_excel_analysis_developed.py, lines 509–554 · score 0.81 · artificially inflating brain, age related brain, skull normalized ratios, mask area, ventricular ratio, denominator
[9] § Results › Total lesion area analysis ↔ Phase5_statistical_analysis/p4_excel_analysis_developed.py, lines 509–554 · score 0.76 · artificially inflating brain, age related brain, skull normalized ratios, mask area, denominator, atrophy
[10] § Materials and methods › Deep learning architectures › Training protocol › Training hyperparameters ↔ Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_variant_all_net.py, lines 163–252 · score 0.71 · ReduceLROnPlateau, hyperparameters, patience, optimizer, Batch, Training
[11] § Materials and methods › Deep learning architectures › Architecture selection and implementation › U-Net Architecture ↔ Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/dlv3_unet_model_GN.py, lines 12–54 · score 0.66 · ReLU, encoder decoder, channel, blocks, activation, convolutions
[12] § Materials and methods › Statistical analysis framework › Primary outcome measures › Total lesion area analysis ↔ Phase4_data_processing/p4_brain_mri_extractor_excel.py, lines 147–228 · score 0.63 · skull mask, brain mask, WMH area, Ventricular area, summed, pixel
[13] § Results › Ventricular area analysis ↔ Phase5_statistical_analysis/p4_excel_analysis_developed.py, lines 917–973 · score 0.62 · brain tissue mask, absolute ventricular area, ventricular ratios normalized, mask area, mm2, skull
[14] § Materials and methods › Study design and population ↔ Phase5_statistical_analysis/p4_excel_analysis_developed.py, lines 4111–4175 · score 0.61 · McDonald, scanner, neurological, quality, confidentiality, diagnosis
[15] § Results › Deep learning model performance ↔ Phase3_model_training_and_inferencing_and_evaluation/for_GM/model_training_scripts/p1_data_loader.py, lines 1–40 · score 0.60 · fold cross validation, white matter hyperintensity, deep learning, neuroanatomical, models, segmentation
[16] § Materials and methods › Statistical analysis framework › Demographic analysis ↔ Phase5_statistical_analysis/p4_excel_analysis_developed.py, lines 800–879 · score 0.60 · Chi square, Gender distribution, variables, Demographic, female, HC
[17] § Materials and methods › Deep learning architectures › Performance evaluation ↔ Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_inference.py, lines 232–304 · score 0.59 · positive lesion, ground truth, union, connected, overlap, predicted
[18] § Materials and methods › Deep learning architectures › Architecture selection and implementation › DeepLabV3Plus ↔ Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/dlv3_unet_model_GN.py, lines 56–99 · score 0.58 · atrous spatial pyramid, ASPP
[19] § Materials and methods › Deep learning architectures › Architecture selection and implementation › U-Net Architecture ↔ Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/trans_unet_model.py, the whole file · a weak match · score 0.55 · transposed, network, dropout, ReLU, encoder, blocks
[20] § Materials and methods › Statistical analysis framework › Advanced statistical methods ↔ Phase5_statistical_analysis/p4_excel_analysis_developed.py, lines 2365–2423 · score 0.54 · WMH ratio, medium, biserial, Wilk, Cohen, threshold
[21] § Materials and methods › MRI datasets and ground truth › Ground truth annotation process ↔ Phase4_data_processing/p4_predict_new_data.py, lines 175–225 · score 0.53 · morphological post processing, Abnormal WMH, trained, classes, ventricles
[22] § Materials and methods › Implementation details ↔ Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_inference.py, lines 725–800 · score 0.52 · Post processing, radius, disk, inference, kernel, morphological
[23] § Materials and methods › Deep learning architectures › Training protocol › Dataset composition ↔ Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_data_loader.py, lines 812–890 · score 0.50 · stratified splitting, training patients, MSSEG, slices, model

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 4,176 lines · 196 KB · MIT · 9 matches

"""
P4 Article - MS Brain MRI Statistical Analysis Framework
==========================================
A comprehensive statistical analysis and visualization system for Multiple Sclerosis
brain MRI studies, designed for publication in prestigious journals.
This script processes real MRI analysis data and generates:
- Demographic analysis and comparisons
- Ventricular and lesion burden analysis
- Publication-quality figures and tables
- Statistical comparisons between groups
Author: Research Team
Developer:
Mahdi Bashiri Bawil
"""
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as stats
from scipy import stats
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
import warnings
warnings.filterwarnings('ignore')
# Set matplotlib parameters for publication quality
plt.rcParams.update({
'font.size': 14,
'font.family': 'Arial',
'axes.linewidth': 1.2,
'axes.spines.top': False,
'axes.spines.right': False,
'figure.dpi': 300,
'savefig.dpi': 300,
'savefig.bbox': 'tight',
'savefig.transparent': False
})
class OutlierDetector:
"""
Comprehensive outlier detection and filtering for MS brain MRI analysis
"""
def __init__(self, config=None):
self.config = config
self.outlier_summary = {}
self.cleaned_indices = None
def detect_outliers_iqr(self, data, column, factor=1.5):
"""
Detect outliers using Interquartile Range (IQR) method
Parameters:
- data: DataFrame
- column: column name to check
- factor: IQR multiplier (1.5 for mild, 3.0 for extreme outliers)
"""
if column not in data.columns:
return np.array([False] * len(data))
series = data[column].dropna()
Q1 = series.quantile(0.25)
Q3 = series.quantile(0.75)
IQR = Q3 - Q1
lower_bound = Q1 - factor * IQR
upper_bound = Q3 + factor * IQR
outliers = (data[column] < lower_bound) | (data[column] > upper_bound)
return outliers.fillna(False)
def detect_outliers_zscore(self, data, column, threshold=3.0):
"""
Detect outliers using Z-score method
Parameters:
- data: DataFrame
- column: column name to check
- threshold: Z-score threshold (typically 2.5-3.0)
"""
if column not in data.columns:
return np.array([False] * len(data))
series = data[column].dropna()
if len(series) == 0:
return np.array([False] * len(data))
z_scores = np.abs(stats.zscore(series, nan_policy='omit'))
outliers = pd.Series([False] * len(data), index=data.index)
outliers.loc[series.index] = z_scores > threshold
return outliers.fillna(False)
def detect_outliers_modified_zscore(self, data, column, threshold=3.5):
"""
Detect outliers using Modified Z-score (more robust to outliers)
Uses median absolute deviation (MAD) instead of standard deviation
"""
if column not in data.columns:
return np.array([False] * len(data))
series = data[column].dropna()
if len(series) == 0:
return np.array([False] * len(data))
median = series.median()
mad = np.median(np.abs(series - median))
# Avoid division by zero
if mad == 0:
mad = 1e-10
modified_z_scores = 0.6745 * (series - median) / mad
outliers = pd.Series([False] * len(data), index=data.index)
outliers.loc[series.index] = np.abs(modified_z_scores) > threshold
return outliers.fillna(False)
def detect_multivariate_outliers(self, data, columns, contamination=0.1):
"""
Detect multivariate outliers using Isolation Forest
Parameters:
- data: DataFrame
- columns: list of columns to consider
- contamination: expected proportion of outliers
"""
try:
from sklearn.ensemble import IsolationForest
# Select only numeric columns that exist
valid_columns = [col for col in columns if col in data.columns]
if len(valid_columns) == 0:
return np.array([False] * len(data))
# Prepare data for isolation forest
X = data[valid_columns].copy()
# Handle missing values
for col in valid_columns:
X[col] = X[col].fillna(X[col].median())
# Fit Isolation Forest
iso_forest = IsolationForest(contamination=contamination, random_state=42)
outlier_labels = iso_forest.fit_predict(X)
# Convert to boolean array (True for outliers)
return outlier_labels == -1
except ImportError:
print("Warning: sklearn not available for multivariate outlier detection")
return np.array([False] * len(data))
def comprehensive_outlier_detection(self, data, methods='all', age_stratified=True):
"""
Comprehensive outlier detection using multiple methods
Parameters:
- data: DataFrame
- methods: 'all', 'conservative', or list of specific methods
- age_stratified: whether to perform outlier detection within age groups
"""
print("\n" + "=" * 60)
print("COMPREHENSIVE OUTLIER DETECTION")
print("=" * 60)
# Define columns to check for outliers
outlier_columns = {
'continuous': [
self.config.COLUMNS['age'],
self.config.COLUMNS['total_skull'],
self.config.COLUMNS['total_brain'],
self.config.COLUMNS['total_ventricle'],
self.config.COLUMNS['total_wmh'],
'VentricleRatio',
'VentricleRatio_Skull',
'WMHRatio_Skull',
'WMHRatio'
],
'wmh_subtypes': [
self.config.COLUMNS.get('peri_wmh', 'TotalPeriArea'),
self.config.COLUMNS.get('deep_wmh', 'TotalDeepArea'),
self.config.COLUMNS.get('juxta_wmh', 'TotalJuxtArea'),
'peri_wmh_ratio',
'deep_wmh_ratio',
'juxta_wmh_ratio'
]
}
all_columns = outlier_columns['continuous'] + outlier_columns['wmh_subtypes']
# Initialize outlier tracking
outlier_flags = pd.DataFrame(index=data.index)
outlier_summary = {}
if age_stratified:
# Detect outliers within each age group and study group
print("Performing age-stratified outlier detection...")
for group in ['HC', 'MS']:
for age_group in data['AgeGroup'].cat.categories:
if pd.isna(age_group):
continue
subset_mask = (data[self.config.COLUMNS['group']] == group) & \
(data['AgeGroup'] == age_group)
subset_data = data[subset_mask]
if len(subset_data) < 5: # Skip if too few samples
continue
print(f" {group} - {age_group}: {len(subset_data)} patients")
for column in all_columns:
if column not in subset_data.columns:
continue
# Apply multiple detection methods
col_name = f"{column}_{group}_{age_group}"
# IQR method (conservative)
iqr_outliers = self.detect_outliers_iqr(subset_data, column, factor=2.0)
# Modified Z-score (more robust)
mz_outliers = self.detect_outliers_modified_zscore(subset_data, column, threshold=3.5)
# Combine methods (conservative approach)
combined_outliers = iqr_outliers & mz_outliers
if combined_outliers.sum() > 0:
outlier_flags.loc[subset_data.index, col_name] = combined_outliers
outlier_summary[col_name] = {
'count': combined_outliers.sum(),
'percentage': (combined_outliers.sum() / len(subset_data)) * 100,
'indices': subset_data.index[combined_outliers].tolist()
}
else:
# Global outlier detection
print("Performing global outlier detection...")
for group in ['HC', 'MS']:
group_data = data[data[self.config.COLUMNS['group']] == group]
print(f" {group}: {len(group_data)} patients")
for column in all_columns:
if column not in group_data.columns:
continue
col_name = f"{column}_{group}"
# Multiple detection methods
iqr_outliers = self.detect_outliers_iqr(group_data, column, factor=1.5)
mz_outliers = self.detect_outliers_modified_zscore(group_data, column, threshold=3.0)
# Conservative combination (both methods must agree)
combined_outliers = iqr_outliers & mz_outliers
if combined_outliers.sum() > 0:
outlier_flags.loc[group_data.index, col_name] = combined_outliers
outlier_summary[col_name] = {
'count': combined_outliers.sum(),
'percentage': (combined_outliers.sum() / len(group_data)) * 100,
'indices': group_data.index[combined_outliers].tolist()
}
# Identify patients with multiple outlier flags
outlier_flags = outlier_flags.fillna(False)
outlier_counts = outlier_flags.sum(axis=1)
# Define threshold for removing patients (e.g., outliers in 3+ variables)
outlier_threshold = 3
patients_to_remove = outlier_counts >= outlier_threshold
print(f"\nOUTLIER DETECTION SUMMARY:")
print(f"{'=' * 40}")
print(f"Total outlier flags detected: {outlier_flags.sum().sum()}")
print(f"Patients with {outlier_threshold}+ outlier flags: {patients_to_remove.sum()}")
print(f"Percentage of data to be removed: {(patients_to_remove.sum() / len(data)) * 100:.2f}%")
# Age group breakdown
if 'AgeGroup' in data.columns:
age_outlier_summary = data[patients_to_remove].groupby([self.config.COLUMNS['group'], 'AgeGroup']).size()
print(f"\nOutliers by group and age:")
for (group, age), count in age_outlier_summary.items():
print(f" {group} - {age}: {count} patients")
self.outlier_summary = {
'detailed': outlier_summary,
'patients_to_remove': data.index[patients_to_remove].tolist(),
'outlier_flags': outlier_flags,
'outlier_counts': outlier_counts
}
return patients_to_remove
def visualize_outliers(self, data, outliers_mask, save_path=None):
"""
Create visualizations showing detected outliers
"""
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Outlier Detection Visualization', fontsize=16, fontweight='bold')
# Select key variables for visualization
viz_columns = [
('VentricleRatio', 'Ventricular Ratio (%)'),
('WMHRatio', 'WMH Ratio (%)'),
(self.config.COLUMNS['total_brain'], 'Total Brain Area'),
(self.config.COLUMNS['total_ventricle'], 'Total Ventricle Area'),
(self.config.COLUMNS['total_wmh'], 'Total WMH Area'),
(self.config.COLUMNS['age'], 'Age (years)')
]
for idx, (column, title) in enumerate(viz_columns):
if column not in data.columns:
continue
row, col = idx // 3, idx % 3
ax = axes[row, col]
# Create boxplot with outliers highlighted
for group in ['HC', 'MS']:
group_data = data[data[self.config.COLUMNS['group']] == group]
group_outliers = outliers_mask[group_data.index]
# Normal data points
normal_data = group_data.loc[~group_outliers, column].dropna()
outlier_data = group_data.loc[group_outliers, column].dropna()
# Plot boxplot
bp = ax.boxplot([normal_data], positions=[0 if group == 'HC' else 1],
widths=0.6, patch_artist=True,
labels=[group])
# Color boxes
color = self.config.COLORS['hc'] if group == 'HC' else self.config.COLORS['ms']
bp['boxes'][0].set_facecolor(color)
bp['boxes'][0].set_alpha(0.7)
# Highlight outliers
if len(outlier_data) > 0:
y_pos = 0 if group == 'HC' else 1
ax.scatter([y_pos] * len(outlier_data), outlier_data,
color='red', s=50, alpha=0.8, marker='x',
label=f'{group} Outliers' if idx == 0 else "")
ax.set_title(title, fontweight='bold')
ax.set_xticklabels(['HC', 'MS'])
ax.grid(True, alpha=0.3)
if idx == 0:
ax.legend()
plt.tight_layout()
if save_path:
plt.savefig(save_path, dpi=300, bbox_inches='tight')
return fig
def generate_outlier_report(self, data, outliers_mask, save_path=None):
"""
Generate a detailed outlier report
"""
report_data = {
'PatientID': [],
'Group': [],
'AgeGroup': [],
'Gender': [],
'OutlierCount': [],
'OutlierVariables': []
}
outlier_patients = data[outliers_mask]
for idx in outlier_patients.index:
patient_flags = self.outlier_summary['outlier_flags'].loc[idx]
outlier_vars = patient_flags[patient_flags == True].index.tolist()
report_data['PatientID'].append(data.loc[idx, self.config.COLUMNS['patient_id']])
report_data['Group'].append(data.loc[idx, self.config.COLUMNS['group']])
report_data['AgeGroup'].append(data.loc[idx, 'AgeGroup'])
report_data['Gender'].append(data.loc[idx, 'Gender'])
report_data['OutlierCount'].append(len(outlier_vars))
report_data['OutlierVariables'].append('; '.join(outlier_vars))
report_df = pd.DataFrame(report_data)
if save_path:
report_df.to_csv(save_path, index=False)
return report_df
class MSAnalysisConfig:
"""Configuration class for MS analysis parameters"""
# File paths
# Define file paths
header_dir = r"E:\MBashiri\ours_articles\Paper#Stats"
header_dir = r"D:\TEMP_P4"
DATA_PATH = os.path.join(header_dir, "brain_mri_analysis_results_ALL.csv")
OUTPUT_DIR = os.path.join(header_dir, "csv_analysis_outputs_no_outlier_v4")
os.makedirs(OUTPUT_DIR, exist_ok=True)
# Age stratification
AGE_BINS = [(18, 29), (30, 39), (40, 49), (50, 59), (60, 69)]
AGE_LABELS = ['18-29', '30-39', '40-49', '50-59', '60+']
# Color schemes for publication
COLORS = {
'male': '#1f77b4', # Blue
'female': '#d62728', # Red
'hc': '#2ca02c', # Green
'ms': '#ff7f0e', # Orange
'pewmh': '#8B0000', # Dark Red
'dwmh': '#FF8C00', # Orange
'jcwmh': '#FFD700' # Gold
}
# Figure settings
FIGURE_SIZE = (12, 8)
DPI = 300
# Statistical parameters
ALPHA = 0.05
CONFIDENCE_INTERVAL = 0.95
# Column mappings
COLUMNS = {
'patient_id': 'PatientID',
'age': 'PatientAge',
'sex': 'PatientSex', # 0: Male, 1: Female
'group': 'StudyGroup', # HC: Healthy Control, MS: Multiple Sclerosis
'total_brain': 'TotalBrainArea', # Brain-tissue mask area (FSL BET-derived)
'total_skull': 'TotalSkullArea', # Total head/skull mask area (morphological head mask)
'total_ventricle': 'TotalVentricleArea',
'total_wmh': 'TotalWMHArea',
'peri_wmh': 'TotalPeriArea',
'deep_wmh': 'TotalDeepArea',
'juxta_wmh': 'TotalJuxtArea'
}
# NOTE on denominators:
# 'total_brain' = sum of BET brain binary mask pixels — represents visible brain
# tissue area across slices. Primary normalisation denominator.
# 'total_skull' = sum of morphological head binary mask pixels — represents the
# total head cross-section (all tissues inside skull). Using this
# denominator avoids inflating age-related lesion ratios via
# brain atrophy, since the skull area is stable across age.
# If 'TotalSkullArea' is absent from the CSV, skull-normalised columns are skipped gracefully.
# Polynomial regression age limit for Figure 4 scatter plots (Reviewer 1, Comment 11).
# Trend lines will be rendered only up to this age; all data points are still plotted.
# Set to None to disable the limit and extend lines over the full age range.
REGRESSION_AGE_LIMIT = 60 # years
def add_outlier_detection_to_load_data(original_load_data_method):
"""
Modified load_data method with integrated outlier detection
"""
def load_data_with_outlier_detection(self, filepath=None):
"""Load and preprocess the MRI analysis data with outlier detection"""
filepath = filepath or self.config.DATA_PATH
try:
# Load data using original method logic
try:
self.data = pd.read_csv(filepath)
except:
self.data = pd.read_csv(filepath, sep='\t')
print(f"Data loaded successfully: {self.data.shape[0]} patients, {self.data.shape[1]} features")
# Validate columns
required_cols = list(self.config.COLUMNS.values())
missing_cols = [col for col in required_cols if col not in self.data.columns]
if missing_cols:
print(f"Warning: Missing columns: {missing_cols}")
# Create age groups
self.data['AgeGroup'] = pd.cut(self.data[self.config.COLUMNS['age']],
bins=[b[0] for b in self.config.AGE_BINS] + [self.config.AGE_BINS[-1][1]],
labels=self.config.AGE_LABELS,
include_lowest=True, right=False)
# Create gender labels
self.data['Gender'] = self.data[self.config.COLUMNS['sex']].map({0: 'Female', 1: 'Male'})
# Define epsilon to prevent division by zero
epsilon = 1e-10
# ── DENOMINATOR 1: Brain-tissue mask area (FSL BET) ────────────────────────────────
# Sum of binary brain mask pixels; used as the primary normalisation denominator
# throughout the manuscript.
denominator_brain = self.data[self.config.COLUMNS['total_brain']] + epsilon
# denominator_brain = denominator_brain /100 # enable for more zooming
# ── DENOMINATOR 2: Total skull / head mask area (Reviewer 2, Major Comments 4 & 5) ─
# Sum of binary head-mask pixels (morphological background separation). Using this
# as the denominator controls for age-related brain atrophy that could otherwise
# artificially inflate brain-normalised ratios over time even without true lesion
# accumulation. Absent column is handled gracefully below.
skull_col = self.config.COLUMNS.get('total_skull', 'TotalSkullArea')
_has_skull = skull_col in self.data.columns
if _has_skull:
denominator_skull = self.data[skull_col] + epsilon
# denominator_skull = denominator_skull / 100 # enable for more zooming
print(f" [Normalisation] Skull/head-mask column '{skull_col}' found — "
f"skull-normalised ratios will be computed alongside brain-normalised ratios.")
else:
denominator_skull = None
print(f" [Normalisation] Column '{skull_col}' NOT found in data — "
f"skull-normalised ratios will be skipped (sensitivity analysis unavailable).")
# ── Ventricle ratios ──────────────────────────────────────────────────────────────
self.data['VentricleRatio'] = (self.data[self.config.COLUMNS['total_ventricle']] /
denominator_brain * 100)
if denominator_skull is not None:
self.data['VentricleRatio_Skull'] = (self.data[self.config.COLUMNS['total_ventricle']] /
denominator_skull * 100)
# ── Total WMH ratios ──────────────────────────────────────────────────────────────
self.data['WMHRatio'] = (self.data[self.config.COLUMNS['total_wmh']] /
denominator_brain * 100)
if denominator_skull is not None:
self.data['WMHRatio_Skull'] = (self.data[self.config.COLUMNS['total_wmh']] /
denominator_skull * 100)
# ── WMH subtype proportional ratios (normalised by total WMH area) ───────────────
# These express each subtype as a share of total lesion burden.
denominator_tic_wmh = self.data[self.config.COLUMNS['total_wmh']] + epsilon
for subtype in ['peri_wmh', 'deep_wmh', 'juxta_wmh']:
self.data[f'{subtype}_ratio'] = (self.data[self.config.COLUMNS[subtype]] /
denominator_tic_wmh * 100)
# ── WMH subtype absolute ratios vs brain area (and vs skull area if available) ────
for subtype in ['peri_wmh', 'deep_wmh', 'juxta_wmh']:
self.data[f'{subtype}_ratio_brain'] = (self.data[self.config.COLUMNS[subtype]] /
denominator_brain * 100)
if denominator_skull is not None:
self.data[f'{subtype}_ratio_skull'] = (self.data[self.config.COLUMNS[subtype]] /
denominator_skull * 100)
# ── Safety: replace any remaining NaN or Inf values with 0 ──────────────────────
skull_ratio_cols = (
['VentricleRatio_Skull', 'WMHRatio_Skull'] +
[f'{s}_ratio_skull' for s in ['peri_wmh', 'deep_wmh', 'juxta_wmh']]
) if denominator_skull is not None else []
numeric_cols = (
['VentricleRatio', 'WMHRatio'] +
[f'{subtype}_ratio' for subtype in ['peri_wmh', 'deep_wmh', 'juxta_wmh']] +
[f'{subtype}_ratio_brain' for subtype in ['peri_wmh', 'deep_wmh', 'juxta_wmh']] +
skull_ratio_cols
)
for col in numeric_cols:
if col in self.data.columns:
self.data[col] = self.data[col].replace([np.inf, -np.inf], np.nan).fillna(0)
# === NEW: OUTLIER DETECTION ===
print("\nPerforming outlier detection...")
outlier_detector = OutlierDetector(self.config)
# Detect outliers using comprehensive method
outliers_mask = outlier_detector.comprehensive_outlier_detection(
self.data,
age_stratified=True # Detect outliers within age groups
)
# Create visualizations
outlier_viz_path = os.path.join(self.config.OUTPUT_DIR, 'outlier_detection_visualization.png')
outlier_detector.visualize_outliers(self.data, outliers_mask, outlier_viz_path)
# Generate outlier report
outlier_report_path = os.path.join(self.config.OUTPUT_DIR, 'outlier_report.csv')
outlier_report = outlier_detector.generate_outlier_report(self.data, outliers_mask, outlier_report_path)
# Store original data for reference
self.data_with_outliers = self.data.copy()
self.outlier_info = {
'detector': outlier_detector,
'outliers_mask': outliers_mask,
'outlier_report': outlier_report,
'removed_count': outliers_mask.sum(),
'original_count': len(self.data)
}
# Remove outliers from main dataset
self.data = self.data[~outliers_mask].reset_index(drop=True)
print(f"\nOUTLIER REMOVAL SUMMARY:")
print(f"Original dataset: {self.outlier_info['original_count']} patients")
print(f"Outliers removed: {self.outlier_info['removed_count']} patients")
print(f"Final dataset: {len(self.data)} patients")
print(f"Data retention: {(len(self.data) / self.outlier_info['original_count']) * 100:.1f}%")
# Recreate age groups after outlier removal (in case categories changed)
self.data['AgeGroup'] = pd.cut(self.data[self.config.COLUMNS['age']],
bins=[b[0] for b in self.config.AGE_BINS] + [self.config.AGE_BINS[-1][1]],
labels=self.config.AGE_LABELS,
include_lowest=True, right=False)
return self.data
except Exception as e:
print(f"Error loading data: {e}")
return None
return load_data_with_outlier_detection
class MSStatisticalAnalysis:
"""Main class for MS statistical analysis and visualization"""
def __init__(self, config=None):
self.config = config or MSAnalysisConfig()
self.data = None
self.results = {}
self.outlier_removal = False
def load_data(self, filepath=None):
"""Load and preprocess the MRI analysis data"""
filepath = filepath or self.config.DATA_PATH
try:
# Try to read as CSV first, then as tab-separated
try:
self.data = pd.read_csv(filepath)
except:
self.data = pd.read_csv(filepath, sep='\t')
print(f"Data loaded successfully: {self.data.shape[0]} patients, {self.data.shape[1]} features")
# Validate columns
required_cols = list(self.config.COLUMNS.values())
missing_cols = [col for col in required_cols if col not in self.data.columns]
if missing_cols:
print(f"Warning: Missing columns: {missing_cols}")
# Create age groups
self.data['AgeGroup'] = pd.cut(self.data[self.config.COLUMNS['age']],
bins=[b[0] for b in self.config.AGE_BINS] + [self.config.AGE_BINS[-1][1]],
labels=self.config.AGE_LABELS,
include_lowest=True, right=False)
# Create gender labels
self.data['Gender'] = self.data[self.config.COLUMNS['sex']].map({0: 'Female', 1: 'Male'})
# Define epsilon to prevent division by zero
epsilon = 1e-10
# ── DENOMINATOR 1: Brain-tissue mask area (FSL BET) ────────────────────────────────
# Sum of binary brain mask pixels; used as the primary normalisation denominator
# throughout the manuscript.
denominator_brain = self.data[self.config.COLUMNS['total_brain']] + epsilon
# denominator_brain = denominator_brain /100 # enable for more zooming
# ── DENOMINATOR 2: Total skull / head mask area (Reviewer 2, Major Comments 4 & 5) ─
# Sum of binary head-mask pixels (morphological background separation). Using this
# as the denominator controls for age-related brain atrophy that could otherwise
# artificially inflate brain-normalised ratios over time even without true lesion
# accumulation. Absent column is handled gracefully below.
skull_col = self.config.COLUMNS.get('total_skull', 'TotalSkullArea')
_has_skull = skull_col in self.data.columns
if _has_skull:
denominator_skull = self.data[skull_col] + epsilon
# denominator_skull = denominator_skull / 100 # enable for more zooming
print(f" [Normalisation] Skull/head-mask column '{skull_col}' found — "
f"skull-normalised ratios will be computed alongside brain-normalised ratios.")
else:
denominator_skull = None
print(f" [Normalisation] Column '{skull_col}' NOT found in data — "
f"skull-normalised ratios will be skipped (sensitivity analysis unavailable).")
# ── Ventricle ratios ──────────────────────────────────────────────────────────────
self.data['VentricleRatio'] = (self.data[self.config.COLUMNS['total_ventricle']] /
denominator_brain * 100)
if denominator_skull is not None:
self.data['VentricleRatio_Skull'] = (self.data[self.config.COLUMNS['total_ventricle']] /
denominator_skull * 100)
# ── Total WMH ratios ──────────────────────────────────────────────────────────────
self.data['WMHRatio'] = (self.data[self.config.COLUMNS['total_wmh']] /
denominator_brain * 100)
if denominator_skull is not None:
self.data['WMHRatio_Skull'] = (self.data[self.config.COLUMNS['total_wmh']] /
denominator_skull * 100)
# ── WMH subtype proportional ratios (normalised by total WMH area) ───────────────
# These express each subtype as a share of total lesion burden.
denominator_tic_wmh = self.data[self.config.COLUMNS['total_wmh']] + epsilon
for subtype in ['peri_wmh', 'deep_wmh', 'juxta_wmh']:
self.data[f'{subtype}_ratio'] = (self.data[self.config.COLUMNS[subtype]] /
denominator_tic_wmh * 100)
# ── WMH subtype absolute ratios vs brain area (and vs skull area if available) ────
for subtype in ['peri_wmh', 'deep_wmh', 'juxta_wmh']:
self.data[f'{subtype}_ratio_brain'] = (self.data[self.config.COLUMNS[subtype]] /
denominator_brain * 100)
if denominator_skull is not None:
self.data[f'{subtype}_ratio_skull'] = (self.data[self.config.COLUMNS[subtype]] /
denominator_skull * 100)
# ── Safety: replace any remaining NaN or Inf values with 0 ──────────────────────
skull_ratio_cols = (
['VentricleRatio_Skull', 'WMHRatio_Skull'] +
[f'{s}_ratio_skull' for s in ['peri_wmh', 'deep_wmh', 'juxta_wmh']]
) if denominator_skull is not None else []
numeric_cols = (
['VentricleRatio', 'WMHRatio'] +
[f'{subtype}_ratio' for subtype in ['peri_wmh', 'deep_wmh', 'juxta_wmh']] +
[f'{subtype}_ratio_brain' for subtype in ['peri_wmh', 'deep_wmh', 'juxta_wmh']] +
skull_ratio_cols
)
for col in numeric_cols:
if col in self.data.columns:
self.data[col] = self.data[col].replace([np.inf, -np.inf], np.nan).fillna(0)
if self.outlier_removal:
# === NEW: ADD THIS SECTION BEFORE THE FINAL RETURN ===
print("\nPerforming outlier detection...")
outlier_detector = OutlierDetector(self.config)
# Detect outliers using comprehensive method
outliers_mask = outlier_detector.comprehensive_outlier_detection(
self.data,
age_stratified=True # Detect outliers within age groups
)
# Create visualizations
outlier_viz_path = os.path.join(self.config.OUTPUT_DIR, 'outlier_detection_visualization.png')
outlier_detector.visualize_outliers(self.data, outliers_mask, outlier_viz_path)
# Generate outlier report
outlier_report_path = os.path.join(self.config.OUTPUT_DIR, 'outlier_report.csv')
outlier_report = outlier_detector.generate_outlier_report(self.data, outliers_mask, outlier_report_path)
# Store original data for reference
self.data_with_outliers = self.data.copy()
self.outlier_info = {
'detector': outlier_detector,
'outliers_mask': outliers_mask,
'outlier_report': outlier_report,
'removed_count': outliers_mask.sum(),
'original_count': len(self.data)
}
# Remove outliers from main dataset
self.data = self.data[~outliers_mask].reset_index(drop=True)
print(f"\nOUTLIER REMOVAL SUMMARY:")
print(f"Original dataset: {self.outlier_info['original_count']} patients")
print(f"Outliers removed: {self.outlier_info['removed_count']} patients")
print(f"Final dataset: {len(self.data)} patients")
print(f"Data retention: {(len(self.data) / self.outlier_info['original_count']) * 100:.1f}%")
# Recreate age groups after outlier removal
self.data['AgeGroup'] = pd.cut(self.data[self.config.COLUMNS['age']],
bins=[b[0] for b in self.config.AGE_BINS] + [
self.config.AGE_BINS[-1][1]],
labels=self.config.AGE_LABELS,
include_lowest=True, right=False)
return self.data
except Exception as e:
print(f"Error loading data: {e}")
return None
def demographic_analysis(self):
"""Perform demographic analysis"""
print("\n" + "=" * 60)
print("DEMOGRAPHIC ANALYSIS")
print("=" * 60)
# Basic demographics
demographic_summary = self.data.groupby([self.config.COLUMNS['group'], 'Gender']).agg({
self.config.COLUMNS['age']: ['count', 'mean', 'std', 'min', 'max'],
self.config.COLUMNS['patient_id']: 'count'
}).round(2)
print("\nDemographic Summary:")
print(demographic_summary)
# Age distribution by group and gender
age_dist = pd.crosstab([self.data[self.config.COLUMNS['group']], self.data['AgeGroup']],
self.data['Gender'])
print("\nAge Distribution by Group and Gender:")
print(age_dist)
# Statistical tests for demographic differences
hc_ages = self.data[self.data[self.config.COLUMNS['group']] == 'HC'][self.config.COLUMNS['age']]
ms_ages = self.data[self.data[self.config.COLUMNS['group']] == 'MS'][self.config.COLUMNS['age']]
hc_ages_male = self.data[(self.data['StudyGroup'] == 'HC') & (self.data['Gender'] == 'Male')][self.config.COLUMNS['age']]
ms_ages_male = self.data[(self.data['StudyGroup'] == 'MS') & (self.data['Gender'] == 'Male')][self.config.COLUMNS['age']]
hc_ages_female = self.data[(self.data['StudyGroup'] == 'HC') & (self.data['Gender'] == 'Female')][self.config.COLUMNS['age']]
ms_ages_female = self.data[(self.data['StudyGroup'] == 'MS') & (self.data['Gender'] == 'Female')][self.config.COLUMNS['age']]
# Age comparison between groups (Male)
age_ttest = stats.ttest_ind(hc_ages_male, ms_ages_male)
print(f"\nAge comparison (HC vs MS) (Male): t-statistic = {age_ttest.statistic:.3f}, p-value = {age_ttest.pvalue}")
# Age comparison between groups (Male)
age_ttest = stats.ttest_ind(hc_ages_female, ms_ages_female)
print(f"Age comparison (HC vs MS) (Female): t-statistic = {age_ttest.statistic:.3f}, p-value = {age_ttest.pvalue}")
# Age comparison between groups
age_ttest = stats.ttest_ind(hc_ages, ms_ages)
print(f"\nAge comparison (HC vs MS) (All): t-statistic = {age_ttest.statistic:.3f}, p-value = {age_ttest.pvalue}")
# Gender distribution chi-square test
gender_crosstab = pd.crosstab(self.data[self.config.COLUMNS['group']], self.data['Gender'])
chi2, p_val, _, _ = stats.chi2_contingency(gender_crosstab)
print(f"Gender distribution (HC vs MS): χ² = {chi2:.3f}, p-value = {p_val}")
self.results['demographics'] = {
'summary': demographic_summary,
'age_distribution': age_dist,
'age_test': age_ttest,
'gender_test': (chi2, p_val)
}
return self.results['demographics']
def _assess_normality_and_choose_test(self, data1, data2, variable_name):
"""
Assess normality and choose appropriate statistical test.
Decision logic:
1. Large sample (n >= 50 both groups) AND CV < 1 in both groups:
→ independent t-test. CLT applies and mean is a meaningful summary.
2. Small sample AND both groups pass Shapiro-Wilk AND |skew| < 2:
→ independent t-test. Genuinely normal small-n data.
3. All other cases:
→ Mann-Whitney U.
CV >= 1 (SD >= mean) indicates an exponential-like distribution where the
mean is not a meaningful central tendency, making the t-test interpretively
inappropriate regardless of the CLT. This correctly selects Mann-Whitney
for WMH variables and t-test for ventricular variables.
"""
from scipy import stats
import numpy as np
LARGE_N_THRESHOLD = 50
CV_THRESHOLD = 1.0 # SD/mean >= 1 → distribution too skewed for t-test
data1_clean = data1.dropna()
data2_clean = data2.dropna()
if len(data1_clean) < 3 or len(data2_clean) < 3:
return None, "insufficient_data", {}, None
n1 = len(data1_clean)
n2 = len(data2_clean)
normality_info = {}
shapiro_both_normal = False
problematic_cv = False
try:
sw_data1 = data1_clean.sample(min(n1, 5000), random_state=42) if n1 > 5000 else data1_clean
sw_data2 = data2_clean.sample(min(n2, 5000), random_state=42) if n2 > 5000 else data2_clean
shapiro1 = stats.shapiro(sw_data1)
shapiro2 = stats.shapiro(sw_data2)
skew1 = stats.skew(data1_clean)
skew2 = stats.skew(data2_clean)
# CV: only meaningful when mean > 0
mean1 = data1_clean.mean()
mean2 = data2_clean.mean()
std1 = data1_clean.std()
std2 = data2_clean.std()
cv1 = (std1 / mean1) if mean1 > 0 else float('inf')
cv2 = (std2 / mean2) if mean2 > 0 else float('inf')
# Distribution is problematic for t-test if CV >= 1 in either group
problematic_cv = (cv1 >= CV_THRESHOLD or cv2 >= CV_THRESHOLD)
shapiro_both_normal = (shapiro1.pvalue > 0.05 and shapiro2.pvalue > 0.05
and abs(skew1) < 2 and abs(skew2) < 2)
normality_info = {
'group1_n': n1,
'group2_n': n2,
'group1_shapiro_p': shapiro1.pvalue,
'group2_shapiro_p': shapiro2.pvalue,
'group1_skewness': skew1,
'group2_skewness': skew2,
'group1_kurtosis': stats.kurtosis(data1_clean),
'group2_kurtosis': stats.kurtosis(data2_clean),
'group1_cv': cv1,
'group2_cv': cv2,
'problematic_cv': problematic_cv,
}
except Exception:
shapiro_both_normal = False
problematic_cv = True # conservative fallback
normality_info = {
'error': 'Could not assess normality',
'group1_n': n1, 'group2_n': n2
}
# ── Test selection ──────────────────────────────────────────────────────
large_sample = (n1 >= LARGE_N_THRESHOLD and n2 >= LARGE_N_THRESHOLD)
use_parametric = (large_sample and not problematic_cv) or \
(not large_sample and shapiro_both_normal)
if use_parametric:
test_result = stats.ttest_ind(data1_clean, data2_clean)
test_type = "parametric"
pooled_std = np.sqrt(((n1 - 1) * data1_clean.var() +
(n2 - 1) * data2_clean.var()) / (n1 + n2 - 2))
effect_size = ((data2_clean.mean() - data1_clean.mean()) / pooled_std
if pooled_std > 0 else 0.0)
effect_size_type = "cohens_d"
if large_sample:
reason = (f"t-test: CLT applies (n1={n1}, n2={n2} >= {LARGE_N_THRESHOLD}) "
f"and CV < {CV_THRESHOLD} in both groups "
f"(CV1={cv1:.3f}, CV2={cv2:.3f}). "
f"Mean is a meaningful summary; t-test is appropriate.")
else:
reason = "t-test: small-n, both groups pass normality criteria."
else:
test_result = stats.mannwhitneyu(data1_clean, data2_clean, alternative='two-sided')
test_type = "non_parametric"
p_clamped = max(test_result.pvalue, 1e-30)
z_score = abs(stats.norm.ppf(p_clamped / 2))
effect_size = z_score / np.sqrt(n1 + n2)
effect_size_type = "rank_biserial_r"
if large_sample and problematic_cv:
reason = (f"Mann-Whitney U: large sample (n1={n1}, n2={n2}) but "
f"CV >= {CV_THRESHOLD} (CV1={cv1:.3f}, CV2={cv2:.3f}). "
f"SD >= mean indicates exponential-like distribution; "
f"mean is not a meaningful summary regardless of CLT.")
else:
reason = (f"Mann-Whitney U: small sample (n1={n1}, n2={n2}) "
f"with non-normal distribution.")
normality_info['test_selection_reason'] = reason
normality_info['large_sample_clt_applied'] = large_sample
return test_result, test_type, normality_info, (effect_size, effect_size_type)
def ventricular_burden_analysis(self):
"""Analyze ventricular burden with age and gender stratification.
Produces three metrics per group:
- Absolute ventricular area (mm²)
- Ventricular ratio normalised by brain-tissue mask area (VentricleRatio, %)
- Ventricular ratio normalised by total skull/head mask area (VentricleRatio_Skull, %)
[only when TotalSkullArea column is present in the data]
"""
print("\n" + "=" * 60)
print("VENTRICULAR BURDEN ANALYSIS")
print("=" * 60)
# Layout: 2 groups × (2 or 3 metrics). Add third column for skull-normalised
# sensitivity analysis if the column is available (Reviewer 2, Major Comments 4 & 5).
_has_skull_ratio = 'VentricleRatio_Skull' in self.data.columns
_n_metric_cols = 3 if _has_skull_ratio else 2
fig, axes = plt.subplots(2, _n_metric_cols, figsize=(8 * _n_metric_cols, 12))
if _n_metric_cols == 2:
axes = np.array(axes) # ensure 2-D indexing works uniformly
fig.patch.set_facecolor('white')
groups = ['HC', 'MS']
age_centers = [np.mean(age_range) for age_range in self.config.AGE_BINS]
# Metrics to plot. A third column (skull-normalised ratio) is added if available,
# providing the sensitivity analysis requested by Reviewer 2, Major Comments 4 & 5.
_has_skull_ratio = 'VentricleRatio_Skull' in self.data.columns
metrics = {
'area': {
'column': self.config.COLUMNS['total_ventricle'],
'ylabel': 'Ventricular Area (mm²)',
'title_suffix': 'Ventricular Area'
},
'ratio': {
'column': 'VentricleRatio',
'ylabel': 'Ventricular Ratio — brain denom. (%)',
'title_suffix': 'Ventricular Ratio (brain-normalised)'
}
}
if _has_skull_ratio:
metrics['ratio_skull'] = {
'column': 'VentricleRatio_Skull',
'ylabel': 'Ventricular Ratio — skull denom. (%)',
'title_suffix': 'Ventricular Ratio (skull-normalised, sensitivity)'
}
else:
print(" [Sensitivity] 'VentricleRatio_Skull' column absent — "
"skull-normalised ventricular ratio panel will be omitted.")
# Dictionary to store all table data
table_data = {
'detailed_stats': {},
'plot_data': {},
'metadata': {
'age_bins': self.config.AGE_BINS,
'age_labels': self.config.AGE_LABELS,
'age_centers': age_centers,
'groups': groups,
'metrics': metrics,
'colors': self.config.COLORS
}
}
# Process each group and metric combination
for group_idx, group in enumerate(groups):
group_data = self.data[self.data[self.config.COLUMNS['group']] == group]
table_data['detailed_stats'][group] = {}
table_data['plot_data'][group] = {}
for metric_idx, (metric_name, metric_info) in enumerate(metrics.items()):
ax = axes[group_idx, metric_idx]
# Initialize storage for this group-metric combination
table_data['detailed_stats'][group][metric_name] = {}
table_data['plot_data'][group][metric_name] = {
'age_centers': age_centers.copy(),
'age_labels': self.config.AGE_LABELS.copy(),
'male_means': [],
'female_means': [],
'male_stds': [],
'female_stds': [],
'male_counts': [],
'female_counts': [],
'combined_means': [],
'male_contributions': [],
'female_contributions': []
}
# Process each age group
for age_idx, age_label in enumerate(self.config.AGE_LABELS):
age_group_data = group_data[group_data['AgeGroup'] == age_label]
# Separate by gender
male_data = age_group_data[age_group_data['Gender'] == 'Male'][metric_info['column']]
female_data = age_group_data[age_group_data['Gender'] == 'Female'][metric_info['column']]
# Calculate statistics for each gender
male_stats = {
'count': len(male_data),
'mean': male_data.mean() if len(male_data) > 0 else np.nan,
'std': male_data.std() if len(male_data) > 0 else np.nan,
'min': male_data.min() if len(male_data) > 0 else np.nan,
'max': male_data.max() if len(male_data) > 0 else np.nan,
'median': male_data.median() if len(male_data) > 0 else np.nan
}
female_stats = {
'count': len(female_data),
'mean': female_data.mean() if len(female_data) > 0 else np.nan,
'std': female_data.std() if len(female_data) > 0 else np.nan,
'min': female_data.min() if len(female_data) > 0 else np.nan,
'max': female_data.max() if len(female_data) > 0 else np.nan,
'median': female_data.median() if len(female_data) > 0 else np.nan
}
# Store detailed statistics
table_data['detailed_stats'][group][metric_name][age_label] = {
'Male': male_stats,
'Female': female_stats
}
# Calculate values for plotting (handle NaN values)
male_mean = male_stats['mean'] if not np.isnan(male_stats['mean']) else 0
female_mean = female_stats['mean'] if not np.isnan(female_stats['mean']) else 0
male_std = male_stats['std'] if not np.isnan(male_stats['std']) else 0
female_std = female_stats['std'] if not np.isnan(female_stats['std']) else 0
# Calculate weighted combined mean and contributions
total_male_sum = male_stats['count'] * male_mean if male_stats['count'] > 0 else 0
total_female_sum = female_stats['count'] * female_mean if female_stats['count'] > 0 else 0
total_subjects = male_stats['count'] + female_stats['count']
if total_subjects > 0:
# True weighted combined mean across both genders
combined_mean = (total_male_sum + total_female_sum) / total_subjects
# Calculate proportional contributions to the combined mean
total_sum = total_male_sum + total_female_sum
if total_sum > 0:
male_contribution = (total_male_sum / total_sum) * combined_mean
female_contribution = (total_female_sum / total_sum) * combined_mean
else:
male_contribution = 0
female_contribution = 0
else:
combined_mean = 0
male_contribution = 0
female_contribution = 0
# Store plot data
plot_data = table_data['plot_data'][group][metric_name]
plot_data['male_means'].append(male_mean)
plot_data['female_means'].append(female_mean)
plot_data['male_stds'].append(male_std)
plot_data['female_stds'].append(female_std)
plot_data['male_counts'].append(male_stats['count'])
plot_data['female_counts'].append(female_stats['count'])
plot_data['combined_means'].append(combined_mean)
plot_data['male_contributions'].append(male_contribution)
plot_data['female_contributions'].append(female_contribution)
# Create the plot
plot_data = table_data['plot_data'][group][metric_name]
# Plot male contribution (bottom layer)
ax.fill_between(age_centers, 0, plot_data['male_contributions'],
color=self.config.COLORS['male'], alpha=0.7, label='Male')
# Plot female contribution (top layer)
ax.fill_between(age_centers, plot_data['male_contributions'],
plot_data['combined_means'],
color=self.config.COLORS['female'], alpha=0.7, label='Female')
# Add error bars if desired (optional - uncomment if needed)
# male_errors = [std/np.sqrt(count) if count > 0 else 0
# for std, count in zip(plot_data['male_stds'], plot_data['male_counts'])]
# female_errors = [std/np.sqrt(count) if count > 0 else 0
# for std, count in zip(plot_data['female_stds'], plot_data['female_counts'])]
# ax.errorbar(age_centers, plot_data['combined_means'],
# yerr=combined_errors, fmt='none', color='black', alpha=0.5)
# Formatting
# Calculate panel letter (A, B, C, D)
panel_idx = group_idx * 3 + metric_idx
panel_letter = chr(65 + panel_idx) # 65 is ASCII for 'A'
ax.set_title(f'{panel_letter}. {group} - {metric_info["title_suffix"]}', fontsize=20, fontweight='bold')
# ax.set_title(f'{group} - {metric_info["title_suffix"]}', fontsize=14, fontweight='bold')
ax.set_xlabel('Age (years)', fontsize=20)
ax.set_ylabel(metric_info['ylabel'], fontsize=20)
ax.legend(fontsize=20)
ax.grid(True, alpha=0.3)
ax.set_xticks(age_centers)
ax.set_xticklabels(self.config.AGE_LABELS, fontsize=20)
# Set reasonable y-axis limits
max_value = max(plot_data['combined_means']) if plot_data['combined_means'] else 0
if max_value > 0:
ax.set_ylim(0, max_value * 1.1)
plt.tight_layout()
plt.savefig(os.path.join(config.OUTPUT_DIR, 'ventricular_burden_analysis.png'),
dpi=self.config.DPI, bbox_inches='tight', facecolor='white')
# Generate comprehensive documentation
self._generate_analysis_documentation(table_data)
# Generate and save tables
self._generate_ventricular_tables(table_data)
# Statistical analysis for both metrics
print(f"\n{'=' * 50}")
print("STATISTICAL COMPARISONS (HC vs MS)")
print(f"{'=' * 50}")
# Analysis for ventricular area
hc_area = self.data[self.data[self.config.COLUMNS['group']] == 'HC'][self.config.COLUMNS['total_ventricle']]
ms_area = self.data[self.data[self.config.COLUMNS['group']] == 'MS'][self.config.COLUMNS['total_ventricle']]
if len(hc_area) > 0 and len(ms_area) > 0:
# Use the new standardized testing approach
area_test, test_type, normality_info, effect_size_info = self._assess_normality_and_choose_test(
hc_area, ms_area, "Ventricular Area"
)
print(f"\nVentricular Area comparison:")
print(f"HC: N={len(hc_area)}, mean ± SD = {hc_area.mean():.2f} ± {hc_area.std():.2f} mm²")
print(f"MS: N={len(ms_area)}, mean ± SD = {ms_area.mean():.2f} ± {ms_area.std():.2f} mm²")
# Print normality test results
if 'group1_shapiro_p' in normality_info:
print(
f"Normality tests: HC p={normality_info['group1_shapiro_p']:.3f}, MS p={normality_info['group2_shapiro_p']:.3f}")
if test_type == "parametric":
print(f"Independent t-test: t-statistic = {area_test.statistic:.3f}, p-value = {area_test.pvalue}")
print(f"Cohen's d = {effect_size_info[0]:.3f}")
elif test_type == "non_parametric":
print(f"Mann-Whitney U test: U-statistic = {area_test.statistic:.3f}, p-value = {area_test.pvalue}")
print(f"Effect size (r) = {effect_size_info[0]:.3f}")
# Also report medians for non-parametric
print(
f"HC: median [IQR] = {hc_area.median():.2f} [{hc_area.quantile(0.25):.2f}-{hc_area.quantile(0.75):.2f}] mm²")
print(
f"MS: median [IQR] = {ms_area.median():.2f} [{ms_area.quantile(0.25):.2f}-{ms_area.quantile(0.75):.2f}] mm²")
else:
print("Insufficient data for ventricular area comparison")
area_test = None
test_type = None
normality_info = {}
effect_size_info = (None, None)
# Analysis for ventricular ratio
hc_ratio = self.data[self.data[self.config.COLUMNS['group']] == 'HC']['VentricleRatio']
ms_ratio = self.data[self.data[self.config.COLUMNS['group']] == 'MS']['VentricleRatio']
if len(hc_ratio) > 0 and len(ms_ratio) > 0:
# Use the new standardized testing approach
ratio_test, ratio_test_type, ratio_normality_info, ratio_effect_size_info = self._assess_normality_and_choose_test(
hc_ratio, ms_ratio, "Ventricular Ratio"
)
print(f"\nVentricular Ratio comparison:")
print(f"HC: N={len(hc_ratio)}, mean ± SD = {hc_ratio.mean():.2f} ± {hc_ratio.std():.2f}%")
print(f"MS: N={len(ms_ratio)}, mean ± SD = {ms_ratio.mean():.2f} ± {ms_ratio.std():.2f}%")
# Print normality test results
if 'group1_shapiro_p' in ratio_normality_info:
print(
f"Normality tests: HC p={ratio_normality_info['group1_shapiro_p']:.3f}, MS p={ratio_normality_info['group2_shapiro_p']:.3f}")
if ratio_test_type == "parametric":
print(
f"Independent t-test: t-statistic = {ratio_test.statistic:.3f}, p-value = {ratio_test.pvalue}")
print(f"Cohen's d = {ratio_effect_size_info[0]:.3f}")
elif ratio_test_type == "non_parametric":
print(
f"Mann-Whitney U test: U-statistic = {ratio_test.statistic:.3f}, p-value = {ratio_test.pvalue}")
print(f"Effect size (r) = {ratio_effect_size_info[0]:.3f}")
# Also report medians for non-parametric
print(
f"HC: median [IQR] = {hc_ratio.median():.2f} [{hc_ratio.quantile(0.25):.2f}-{hc_ratio.quantile(0.75):.2f}]%")
print(
f"MS: median [IQR] = {ms_ratio.median():.2f} [{ms_ratio.quantile(0.25):.2f}-{ms_ratio.quantile(0.75):.2f}]%")
else:
print("Insufficient data for ventricular ratio comparison")
ratio_test = None
ratio_test_type = None
ratio_normality_info = {}
ratio_effect_size_info = (None, None)
# Store results (update the existing results storage)
self.results['ventricular_burden'] = {
'area_analysis': {
'hc_stats': (hc_area.mean(), hc_area.std()) if len(hc_area) > 0 else (np.nan, np.nan),
'ms_stats': (ms_area.mean(), ms_area.std()) if len(ms_area) > 0 else (np.nan, np.nan),
'comparison': area_test,
'test_type': test_type,
'normality_info': normality_info,
'effect_size': effect_size_info
},
'ratio_analysis': {
'hc_stats': (hc_ratio.mean(), hc_ratio.std()) if len(hc_ratio) > 0 else (np.nan, np.nan),
'ms_stats': (ms_ratio.mean(), ms_ratio.std()) if len(ms_ratio) > 0 else (np.nan, np.nan),
'comparison': ratio_test,
'test_type': ratio_test_type,
'normality_info': ratio_normality_info,
'effect_size': ratio_effect_size_info
},
'table_data': table_data
}
return self.results['ventricular_burden']
def _generate_analysis_documentation(self, table_data):
"""Generate comprehensive documentation explaining the figure and analysis"""
from datetime import datetime
# Create comprehensive documentation
doc_content = f"""
VENTRICULAR BURDEN ANALYSIS - COMPREHENSIVE DOCUMENTATION
=========================================================
Generated on: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}
OVERVIEW
--------
This analysis examines ventricular burden in the brain comparing Healthy Controls (HC)
and Multiple Sclerosis (MS) patients. The analysis includes both absolute ventricular
area measurements and normalized ventricular ratios, stratified by age groups and gender.
FIGURE DESCRIPTION
------------------
The figure consists of a 2x2 subplot layout:
Layout Structure:
- Top row: Healthy Controls (HC)
- Bottom row: Multiple Sclerosis (MS) patients
- Left column: Absolute Ventricular Area (mm²)
- Right column: Ventricular Ratio (%)
Subplot Details:
1. Top-Left: HC Ventricular Area
2. Top-Right: HC Ventricular Ratio
3. Bottom-Left: MS Ventricular Area
4. Bottom-Right: MS Ventricular Ratio
VISUALIZATION METHOD
--------------------
Chart Type: Stacked Area Plot with Gender Contributions
- Each subplot uses stacked area charts to show gender-stratified data across age groups
- Male contribution (bottom layer): Shows proportional contribution of male subjects to combined mean
- Female contribution (top layer): Shows proportional contribution of female subjects to combined mean
- Total height represents the true weighted combined mean across both genders
- This visualization allows comparison of both absolute values and relative gender contributions
Mathematical Approach:
- Male Contribution = (Male_Count × Male_Mean) / Total_Sum × Combined_Mean
- Female Contribution = (Female_Count × Female_Mean) / Total_Sum × Combined_Mean
- Combined Mean = (Male_Count × Male_Mean + Female_Count × Female_Mean) / (Male_Count + Female_Count)
Color Scheme:
- Male data: {table_data['metadata']['colors'].get('male', 'Blue')} (alpha=0.7)
- Female data: {table_data['metadata']['colors'].get('female', 'Red')} (alpha=0.7)
AGE STRATIFICATION
------------------
Age Groups: {', '.join(table_data['metadata']['age_labels'])}
Age Bins: {table_data['metadata']['age_bins']}
Age Centers (for plotting): {[f'{center:.1f}' for center in table_data['metadata']['age_centers']]}
The analysis stratifies data across these age groups to examine age-related changes
in ventricular burden for both groups and genders.
METRICS ANALYZED
----------------
1. Ventricular Area (mm²):
- Absolute measurement of total ventricular volume
- Column: {table_data['metadata']['metrics']['area']['column']}
- Units: Square millimeters (mm²)
- Clinical significance: Larger values indicate greater ventricular enlargement
2. Ventricular Ratio (%):
- Normalized measurement relative to total brain area/volume
- Column: VentricleRatio
- Units: Percentage (%)
- Clinical significance: Controls for individual brain size differences
ADVANCED STATISTICAL APPROACH
------------------------------
The analysis employs a sophisticated statistical testing framework:
Normality Assessment:
- Shapiro-Wilk tests performed on both groups (HC and MS) for each metric
- Significance threshold: p < 0.05 indicates non-normal distribution
- Results inform choice between parametric and non-parametric tests
Statistical Test Selection:
1. If both groups pass normality: Independent t-test (parametric)
2. If one or both groups fail normality: Mann-Whitney U test (non-parametric)
Effect Size Calculations:
- Parametric tests: Cohen's d
* Small effect: d = 0.2
* Medium effect: d = 0.5
* Large effect: d = 0.8
- Non-parametric tests: Rank-biserial correlation (r ≈ Z/√N, large-sample approximation)
Note: termed "Effect size r" in earlier versions; corrected per Reviewer 2, Minor Comment 4.
* Small effect: r = 0.1
* Medium effect: r = 0.3
* Large effect: r = 0.5
For each combination of:
- Group (HC vs MS)
- Age group ({len(table_data['metadata']['age_labels'])} categories)
- Gender (Male vs Female)
- Metric (Area vs Ratio)
The following statistics are calculated:
- Sample size (N)
- Mean ± Standard Deviation
- Minimum and Maximum values
- Median values
- Interquartile ranges (for non-parametric reporting)
STATISTICAL OUTPUT INTERPRETATION
---------------------------------
Parametric Results (t-test):
- Reports: Mean ± SD for both groups
- Test statistic: t-value and degrees of freedom
- p-value for significance testing
- Cohen's d for effect size magnitude
Non-parametric Results (Mann-Whitney U):
- Reports: Mean ± SD AND Median [IQR] for both groups
- Test statistic: U-value (or equivalent Z-score)
- p-value for significance testing
- Rank-biserial correlation (r) for effect magnitude (Mann-Whitney U; r ≈ Z/√N)
Normality Test Results:
- Shapiro-Wilk p-values reported for each group
- p < 0.05 indicates significant departure from normality
- Informs test selection rationale
INTERPRETATION GUIDELINES
-------------------------
Stacked Area Plot Interpretation:
- Total height = True combined mean (weighted by sample sizes)
- Male layer height = Proportional contribution of males to combined mean
- Female layer height = Proportional contribution of females to combined mean
- Layer proportions reflect both mean values AND sample size contributions
- Steeper slopes indicate rapid changes with age
- Wider differences between groups suggest clinical significance
Clinical Relevance:
- Ventricular enlargement is associated with brain atrophy
- MS patients typically show greater ventricular burden than healthy controls
- Age-related changes may differ between groups
- Gender differences may exist in disease progression patterns
Expected Patterns:
- MS group likely shows higher values than HC group
- Age-related increase in ventricular burden
- Potential gender differences in progression patterns
Statistical Significance Levels:
- p < 0.001: Highly significant (strong evidence)
- p < 0.01: Very significant (moderate to strong evidence)
- p < 0.05: Significant (sufficient evidence)
- p ≥ 0.05: Non-significant (insufficient evidence)
Effect Size Interpretation:
- Cohen's d or r values indicate practical significance
- Large effect sizes may be clinically meaningful even if p > 0.05
- Small p-values with small effect sizes may lack clinical relevance
DATA QUALITY NOTES
-------------------
- Zero values in plots indicate no subjects in that age/gender combination
- Small sample sizes may lead to unstable mean estimates and reduced statistical power
- Standard deviations provide insight into data variability within groups
- Missing data handled by excluding from calculations (listwise deletion)
- Normality violations automatically trigger non-parametric alternatives
STATISTICAL TESTING DETAILS
----------------------------
Overall group comparisons (HC vs MS) performed separately for each metric:
1. Ventricular Area Analysis:
- Automatic normality assessment using Shapiro-Wilk test
- Test selection based on normality results
- Effect size calculation appropriate to test type
- Comprehensive reporting of descriptive statistics
2. Ventricular Ratio Analysis:
- Independent statistical analysis from area measurements
- Same rigorous normality assessment and test selection
- Separate effect size calculations
- Controls for multiple testing considerations
Test Assumptions:
- Parametric tests: Normality, independence, homogeneity of variance
- Non-parametric tests: Independence, similar distributions
- Both assume random sampling from populations of interest
OUTPUT FILES GENERATED
-----------------------
1. Figure: ventricular_burden_analysis.png
- 2x2 subplot layout with stacked area plots showing proportional contributions
- High resolution (DPI: {getattr(self.config, 'DPI', 300)})
- White background for publication quality
2. Enhanced Statistical Tables (CSV):
- ventricular_area_detailed_stats.csv (includes all descriptive statistics)
- ventricular_ratio_detailed_stats.csv (includes all descriptive statistics)
- ventricular_statistical_comparisons.csv (NEW: comprehensive test results)
3. Plot Data Tables (CSV):
- ventricular_area_plot_data.csv (means and sample sizes)
- ventricular_ratio_plot_data.csv (means and sample sizes)
4. Contribution Analysis (CSV):
- ventricular_area_contributions.csv (NEW: proportional contributions)
- ventricular_ratio_contributions.csv (NEW: proportional contributions)
5. Statistical Results Summary (CSV):
- ventricular_normality_results.csv (NEW: normality test outcomes)
- ventricular_effect_sizes.csv (NEW: effect size calculations)
6. This Documentation:
- ventricular_analysis_documentation.txt
TECHNICAL DETAILS
-----------------
Figure Specifications:
- Size: 16" x 12" (width x height)
- DPI: {getattr(self.config, 'DPI', 300)}
- Background: White
- Font sizes: Title=14pt (bold), Axis labels=12pt
- Grid: Enabled with 30% transparency
- Legend: Enabled for each subplot
Statistical Libraries:
- scipy.stats: Shapiro-Wilk, t-test, Mann-Whitney U
- numpy: Mathematical operations and statistical functions
- pandas: Data manipulation and summary statistics
Quality Control:
- Automatic handling of missing values
- Robust error handling for edge cases
- Comprehensive logging of statistical decisions
- Validation of statistical assumptions
ENHANCED LIMITATIONS AND CONSIDERATIONS
---------------------------------------
1. Sample Size Variations:
- Some age/gender combinations may have small sample sizes
- Unbalanced groups may affect statistical power
- Power analysis recommended for study design validation
2. Multiple Comparisons:
- Two separate statistical tests performed (area and ratio)
- Consider Bonferroni correction: α = 0.05/2 = 0.025
- Family-wise error rate may be inflated without correction
3. Age Grouping Effects:
- Discretized age groups may mask continuous age effects
- Loss of information compared to regression approaches
- Age bin boundaries are predetermined and may not reflect natural breakpoints
4. Statistical Assumptions:
- Automatic normality testing with Shapiro-Wilk (sensitive to large samples)
- Independence assumption may be violated in related subjects
- Equal variances assumed for t-tests (consider Welch's t-test alternative)
5. Visualization Limitations:
- Proportional contributions may be difficult to interpret intuitively
- Stacked areas emphasize combined effects over individual gender patterns
- Direct visual comparison between groups requires careful interpretation
6. Clinical Interpretation:
- Statistical significance may not equal clinical significance
- Effect sizes should be considered alongside p-values
- Longitudinal changes not captured in cross-sectional analysis
RECOMMENDED FOLLOW-UP ANALYSES
------------------------------
1. Advanced Statistical Approaches:
- Age as continuous variable (linear/polynomial regression)
- Two-way ANOVA with interaction terms (group × gender × age)
- Mixed-effects models for correlated data
- Bootstrap confidence intervals for robust inference
2. Multiple Comparisons Corrections:
- Bonferroni correction for family-wise error control
- False Discovery Rate (FDR) control for exploratory analyses
- Planned comparisons vs. post-hoc testing strategies
3. Effect Size and Power Analysis:
- Post-hoc power calculations for observed effects
- Sample size calculations for future studies
- Confidence intervals around effect size estimates
4. Alternative Statistical Approaches:
- Bayesian analysis for probabilistic interpretation
- Permutation tests for distribution-free inference
- Robust statistical methods for outlier resistance
5. Clinical Validation:
- Correlation with clinical severity measures
- Longitudinal tracking of ventricular changes
- Predictive modeling for disease progression
QUALITY ASSURANCE CHECKLIST
----------------------------
✓ Normality testing performed automatically
✓ Appropriate statistical test selected based on data properties
✓ Effect sizes calculated and reported
✓ Both parametric and non-parametric results available
✓ Comprehensive descriptive statistics provided
✓ Multiple output formats for different use cases
✓ Documentation includes interpretation guidelines
✓ Limitations and assumptions clearly stated
✓ Recommendations for follow-up analyses provided
CONTACT AND METHODOLOGY
-----------------------
This analysis was generated using an automated pipeline for ventricular burden assessment
with enhanced statistical testing capabilities.
For questions about methodology or interpretation, refer to:
- Original research protocol and statistical analysis plan
- Relevant neuroimaging analysis guidelines
- Statistical consulting resources for complex designs
Analysis Pipeline Version: Enhanced Statistical Testing v2.0
Statistical Methods: Automatic normality assessment with adaptive test selection
Last Updated: {datetime.now().strftime('%Y-%m-%d')}
END OF DOCUMENTATION
====================
"""
# Save documentation to file
doc_filename = os.path.join(config.OUTPUT_DIR, 'ventricular_analysis_documentation.txt')
with open(doc_filename, 'w', encoding='utf-8') as f:
f.write(doc_content)
print(f"\n{'=' * 80}")
print("COMPREHENSIVE DOCUMENTATION GENERATED")
print(f"{'=' * 80}")
print(f"Documentation saved to: ventricular_analysis_documentation.txt")
print(f"File contains detailed explanation of:")
print(f"- Enhanced statistical testing methodology")
print(f"- Normality assessment and test selection")
print(f"- Effect size calculations and interpretation")
print(f"- Figure interpretation and clinical relevance")
def _generate_ventricular_tables(self, table_data):
"""Generate comprehensive tables from the ventricular burden analysis with enhanced statistical reporting"""
# Table 1: Enhanced detailed statistics by group, age, and gender
print(f"\n{'=' * 80}")
print("TABLE 1: DETAILED STATISTICS BY GROUP, AGE, AND GENDER")
print(f"{'=' * 80}")
for metric_name in ['area', 'ratio', 'ratio_skull']:
unit = 'mm²' if metric_name == 'area' else '%'
if metric_name == 'area':
metric_title = 'Ventricular Area'
elif metric_name == 'ratio':
metric_title = 'Ventricular Ratio'
else:
metric_title = 'Ventricular Ratio-Skull'
print(f"\n{metric_title} ({unit}):")
print("-" * 60)
# Create DataFrame for this metric
rows = []
for group in ['HC', 'MS']:
for age_label in self.config.AGE_LABELS:
for gender in ['Male', 'Female']:
stats = table_data['detailed_stats'][group][metric_name][age_label][gender]
rows.append({
'Group': group,
'Age Group': age_label,
'Gender': gender,
'N': stats['count'],
'Mean': f"{stats['mean']:.2f}" if not np.isnan(stats['mean']) else 'N/A',
'SD': f"{stats['std']:.2f}" if not np.isnan(stats['std']) else 'N/A',
'Min': f"{stats['min']:.2f}" if not np.isnan(stats['min']) else 'N/A',
'Max': f"{stats['max']:.2f}" if not np.isnan(stats['max']) else 'N/A',
'Median': f"{stats['median']:.2f}" if not np.isnan(stats['median']) else 'N/A',
'IQR_25': f"{np.nan:.2f}" if np.isnan(
stats['median']) else f"{stats['median'] - stats['std'] / 2:.2f}",
'IQR_75': f"{np.nan:.2f}" if np.isnan(
stats['median']) else f"{stats['median'] + stats['std'] / 2:.2f}"
})
df = pd.DataFrame(rows)
print(df.to_string(index=False))
# Save to CSV
filename = f'ventricular_{metric_name}_detailed_stats.csv'
df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Enhanced table saved to: {filename}")
# Table 2: Statistical Comparisons (NEW - Enhanced with normality and effect size info)
print(f"\n{'=' * 80}")
print("TABLE 2: STATISTICAL COMPARISONS (HC vs MS) - ENHANCED")
print(f"{'=' * 80}")
# Extract statistical results from the analysis
if hasattr(self, 'results') and 'ventricular_burden' in self.results:
ventricular_results = self.results['ventricular_burden']
comparison_rows = []
for metric_type in ['area_analysis', 'ratio_analysis']:
metric_name = 'Area' if metric_type == 'area_analysis' else 'Ratio'
unit = 'mm²' if metric_type == 'area_analysis' else '%'
analysis = ventricular_results[metric_type]
# Extract basic statistics
hc_mean, hc_std = analysis['hc_stats']
ms_mean, ms_std = analysis['ms_stats']
# Extract test results
comparison = analysis['comparison']
test_type = analysis['test_type']
normality_info = analysis.get('normality_info', {})
effect_size_info = analysis.get('effect_size', (None, None))
row = {
'Metric': f'{metric_name} ({unit})',
'HC_Mean': f"{hc_mean:.2f}" if not np.isnan(hc_mean) else 'N/A',
'HC_SD': f"{hc_std:.2f}" if not np.isnan(hc_std) else 'N/A',
'MS_Mean': f"{ms_mean:.2f}" if not np.isnan(ms_mean) else 'N/A',
'MS_SD': f"{ms_std:.2f}" if not np.isnan(ms_std) else 'N/A',
'Test_Type': test_type if test_type else 'N/A',
'Test_Statistic': f"{comparison.statistic:.3f}" if comparison else 'N/A',
'P_Value': f"{comparison.pvalue}" if comparison else 'N/A',
'Effect_Size': f"{effect_size_info[0]:.3f}" if effect_size_info[0] is not None else 'N/A',
'Effect_Size_Type': 'Cohen_d' if test_type == 'parametric' else 'r',
'HC_Normality_p': f"{normality_info.get('group1_shapiro_p', np.nan):.3f}" if 'group1_shapiro_p' in normality_info else 'N/A',
'MS_Normality_p': f"{normality_info.get('group2_shapiro_p', np.nan):.3f}" if 'group2_shapiro_p' in normality_info else 'N/A',
'Normality_Passed': 'Yes' if test_type == 'parametric' else 'No' if test_type == 'non_parametric' else 'N/A'
}
comparison_rows.append(row)
comparison_df = pd.DataFrame(comparison_rows)
print("\nStatistical Comparison Results:")
print("-" * 100)
print(comparison_df.to_string(index=False))
# Save to CSV
comparison_df.to_csv(os.path.join(config.OUTPUT_DIR, 'ventricular_statistical_comparisons.csv'),
index=False)
print(f"\nStatistical comparisons saved to: ventricular_statistical_comparisons.csv")
# Table 3: Plot data (means used for visualization)
print(f"\n{'=' * 80}")
print("TABLE 3: PLOT DATA (MEANS AND CONTRIBUTIONS BY AGE GROUP)")
print(f"{'=' * 80}")
for metric_name in ['area', 'ratio', 'ratio_skull']:
unit = 'mm²' if metric_name == 'area' else '%'
if metric_name == 'area':
metric_title = 'Ventricular Area'
elif metric_name == 'ratio':
metric_title = 'Ventricular Ratio'
else:
metric_title = 'Ventricular Ratio-Skull'
print(f"\n{metric_title} - Mean Values and Contributions Used in Plot ({unit}):")
print("-" * 85)
# Create enhanced plot data table
plot_rows = []
for i, age_label in enumerate(self.config.AGE_LABELS):
age_center = table_data['plot_data']['HC'][metric_name]['age_centers'][i]
# HC data
hc_male_mean = table_data['plot_data']['HC'][metric_name]['male_means'][i]
hc_female_mean = table_data['plot_data']['HC'][metric_name]['female_means'][i]
hc_combined = table_data['plot_data']['HC'][metric_name]['combined_means'][i]
hc_male_contrib = table_data['plot_data']['HC'][metric_name]['male_contributions'][i]
hc_female_contrib = table_data['plot_data']['HC'][metric_name]['female_contributions'][i]
# MS data
ms_male_mean = table_data['plot_data']['MS'][metric_name]['male_means'][i]
ms_female_mean = table_data['plot_data']['MS'][metric_name]['female_means'][i]
ms_combined = table_data['plot_data']['MS'][metric_name]['combined_means'][i]
ms_male_contrib = table_data['plot_data']['MS'][metric_name]['male_contributions'][i]
ms_female_contrib = table_data['plot_data']['MS'][metric_name]['female_contributions'][i]
row = {
'Age_Group': age_label,
'Age_Center': f"{age_center:.1f}",
'HC_Male_Mean': f"{hc_male_mean:.2f}",
'HC_Female_Mean': f"{hc_female_mean:.2f}",
'HC_Combined_Mean': f"{hc_combined:.2f}",
'HC_Male_Contribution': f"{hc_male_contrib:.2f}",
'HC_Female_Contribution': f"{hc_female_contrib:.2f}",
'MS_Male_Mean': f"{ms_male_mean:.2f}",
'MS_Female_Mean': f"{ms_female_mean:.2f}",
'MS_Combined_Mean': f"{ms_combined:.2f}",
'MS_Male_Contribution': f"{ms_male_contrib:.2f}",
'MS_Female_Contribution': f"{ms_female_contrib:.2f}",
'HC_Male_N': table_data['plot_data']['HC'][metric_name]['male_counts'][i],
'HC_Female_N': table_data['plot_data']['HC'][metric_name]['female_counts'][i],
'MS_Male_N': table_data['plot_data']['MS'][metric_name]['male_counts'][i],
'MS_Female_N': table_data['plot_data']['MS'][metric_name]['female_counts'][i]
}
plot_rows.append(row)
plot_df = pd.DataFrame(plot_rows)
print(plot_df.to_string(index=False))
# Save to CSV
filename = f'ventricular_{metric_name}_plot_data.csv'
plot_df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Enhanced plot data saved to: {filename}")
# Table 4: Contribution Analysis (NEW)
print(f"\n{'=' * 80}")
print("TABLE 4: PROPORTIONAL CONTRIBUTION ANALYSIS")
print(f"{'=' * 80}")
for metric_name in ['area', 'ratio']:
unit = 'mm²' if metric_name == 'area' else '%'
metric_title = 'Ventricular Area' if metric_name == 'area' else 'Ventricular Ratio'
print(f"\n{metric_title} - Proportional Contributions ({unit}):")
print("-" * 70)
# Create contribution analysis table
contrib_rows = []
for group in ['HC', 'MS']:
for i, age_label in enumerate(self.config.AGE_LABELS):
male_contrib = table_data['plot_data'][group][metric_name]['male_contributions'][i]
female_contrib = table_data['plot_data'][group][metric_name]['female_contributions'][i]
combined_mean = table_data['plot_data'][group][metric_name]['combined_means'][i]
male_count = table_data['plot_data'][group][metric_name]['male_counts'][i]
female_count = table_data['plot_data'][group][metric_name]['female_counts'][i]
total_count = male_count + female_count
# Calculate proportions
if combined_mean > 0:
male_prop = (male_contrib / combined_mean) * 100 if combined_mean > 0 else 0
female_prop = (female_contrib / combined_mean) * 100 if combined_mean > 0 else 0
else:
male_prop = 0
female_prop = 0
sample_male_prop = (male_count / total_count) * 100 if total_count > 0 else 0
sample_female_prop = (female_count / total_count) * 100 if total_count > 0 else 0
row = {
'Group': group,
'Age_Group': age_label,
'Combined_Mean': f"{combined_mean:.2f}",
'Male_Contribution': f"{male_contrib:.2f}",
'Female_Contribution': f"{female_contrib:.2f}",
'Male_Prop_of_Mean': f"{male_prop:.1f}%",
'Female_Prop_of_Mean': f"{female_prop:.1f}%",
'Male_Sample_Prop': f"{sample_male_prop:.1f}%",
'Female_Sample_Prop': f"{sample_female_prop:.1f}%",
'Total_N': total_count
}
contrib_rows.append(row)
contrib_df = pd.DataFrame(contrib_rows)
print(contrib_df.to_string(index=False))
# Save to CSV
filename = f'ventricular_{metric_name}_contributions.csv'
contrib_df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Contribution analysis saved to: {filename}")
# Table 5: Effect Size Interpretation (NEW)
if hasattr(self, 'results') and 'ventricular_burden' in self.results:
print(f"\n{'=' * 80}")
print("TABLE 5: EFFECT SIZE INTERPRETATION GUIDE")
print(f"{'=' * 80}")
effect_size_rows = []
ventricular_results = self.results['ventricular_burden']
for metric_type in ['area_analysis', 'ratio_analysis']:
metric_name = 'Area' if metric_type == 'area_analysis' else 'Ratio'
analysis = ventricular_results[metric_type]
effect_size_info = analysis.get('effect_size', (None, None))
test_type = analysis['test_type']
if effect_size_info[0] is not None:
effect_size = effect_size_info[0]
if test_type == 'parametric':
# Cohen's d interpretation
if abs(effect_size) < 0.2:
magnitude = "Negligible"
elif abs(effect_size) < 0.5:
magnitude = "Small"
elif abs(effect_size) < 0.8:
magnitude = "Medium"
else:
magnitude = "Large"
else:
# Effect size r interpretation
if abs(effect_size) < 0.1:
magnitude = "Negligible"
elif abs(effect_size) < 0.3:
magnitude = "Small"
elif abs(effect_size) < 0.5:
magnitude = "Medium"
else:
magnitude = "Large"
row = {
'Metric': metric_name,
'Effect_Size_Value': f"{effect_size:.3f}",
'Effect_Size_Type': "Cohen's d" if test_type == 'parametric' else "Rank-biserial correlation (r)",
'Magnitude': magnitude,
'Interpretation': f"{magnitude} effect size indicating {'substantial' if magnitude in ['Medium', 'Large'] else 'minimal'} practical difference"
}
effect_size_rows.append(row)
if effect_size_rows:
effect_df = pd.DataFrame(effect_size_rows)
print("\nEffect Size Interpretations:")
print("-" * 60)
print(effect_df.to_string(index=False))
# Save to CSV
effect_df.to_csv(os.path.join(config.OUTPUT_DIR, 'ventricular_effect_sizes.csv'), index=False)
print(f"\nEffect size interpretations saved to: ventricular_effect_sizes.csv")
# Table 6: Stacked area values (cumulative for visualization)
print(f"\n{'=' * 80}")
print("TABLE 6: STACKED AREA VALUES (FOR AREA CHART)")
print(f"{'=' * 80}")
for metric_name in ['area', 'ratio', 'ratio_skull']:
unit = 'mm²' if metric_name == 'area' else '%'
if metric_name == 'area':
metric_title = 'Ventricular Area'
elif metric_name == 'ratio':
metric_title = 'Ventricular Ratio'
else:
metric_title = 'Ventricular Ratio-Skull'
print(f"\n{metric_title} - Stacked Values ({unit}):")
print("-" * 70)
# Create stacked data table
stacked_rows = []
for group in ['HC', 'MS']:
for i, age_label in enumerate(self.config.AGE_LABELS):
male_mean = table_data['plot_data'][group][metric_name]['male_means'][i]
female_mean = table_data['plot_data'][group][metric_name]['female_means'][i]
row = {
'Group': group,
'Age Group': age_label,
'Male Layer (0 to Male)': f"0.00 to {male_mean:.2f}",
'Female Layer (Male to Total)': f"{male_mean:.2f} to {male_mean + female_mean:.2f}",
'Total Height': f"{male_mean + female_mean:.2f}",
'Male Contribution': f"{male_mean:.2f}",
'Female Contribution': f"{female_mean:.2f}"
}
stacked_rows.append(row)
stacked_df = pd.DataFrame(stacked_rows)
print(stacked_df.to_string(index=False))
# Save to CSV
filename = f'ventricular_{metric_name}_stacked_data.csv'
stacked_df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Stacked data saved to: {filename}")
print(f"\n{'=' * 80}")
print("ENHANCED STATISTICAL TABLES GENERATED!")
print(f"{'=' * 80}")
print("Generated files include:")
print("• Enhanced descriptive statistics with IQR")
print("• Comprehensive statistical comparison results")
print("• Detailed plot data with contributions")
print("• Proportional contribution analysis")
print("• Effect size interpretations")
print("• All tables saved as CSV files for further analysis")
print(f"{'=' * 80}")
def total_lesion_burden_analysis(self):
"""Analyze total lesion burden with age and gender stratification - both area and ratio"""
print("\n" + "=" * 60)
print("TOTAL LESION BURDEN ANALYSIS")
print("=" * 60)
# Layout: 2 groups × (2 or 3 metrics).
_has_wmh_skull = 'WMHRatio_Skull' in self.data.columns
_n_metric_cols = 3 if _has_wmh_skull else 2
fig, axes = plt.subplots(2, _n_metric_cols, figsize=(8 * _n_metric_cols, 12))
fig.patch.set_facecolor('white')
if _n_metric_cols == 2:
axes = np.array(axes)
fig.patch.set_facecolor('white')
groups = ['HC', 'MS']
age_centers = [np.mean(age_range) for age_range in self.config.AGE_BINS]
# Metrics to plot. Skull-normalised ratio added as sensitivity analysis
# (Reviewer 2, Major Comments 4 & 5).
_has_wmh_skull = 'WMHRatio_Skull' in self.data.columns
metrics = {
'area': {
'column': self.config.COLUMNS['total_wmh'],
'ylabel': 'WMH Area (mm²)',
'title_suffix': 'WMH Area'
},
'ratio': {
'column': 'WMHRatio',
'ylabel': 'WMH Ratio — brain denom. (%)',
'title_suffix': 'WMH Ratio (brain-normalised)'
}
}
if _has_wmh_skull:
metrics['ratio_skull'] = {
'column': 'WMHRatio_Skull',
'ylabel': 'WMH Ratio — skull denom. (%)',
'title_suffix': 'WMH Ratio (skull-normalised, sensitivity)'
}
else:
print(" [Sensitivity] 'WMHRatio_Skull' column absent — "
"skull-normalised WMH ratio panel will be omitted.")
# Dictionary to store all table data
table_data = {
'detailed_stats': {},
'plot_data': {},
'metadata': {
'age_bins': self.config.AGE_BINS,
'age_labels': self.config.AGE_LABELS,
'age_centers': age_centers,
'groups': groups,
'metrics': metrics,
'colors': self.config.COLORS
}
}
# Process each group and metric combination
for group_idx, group in enumerate(groups):
group_data = self.data[self.data[self.config.COLUMNS['group']] == group]
table_data['detailed_stats'][group] = {}
table_data['plot_data'][group] = {}
for metric_idx, (metric_name, metric_info) in enumerate(metrics.items()):
ax = axes[group_idx, metric_idx]
# Initialize storage for this group-metric combination
table_data['detailed_stats'][group][metric_name] = {}
table_data['plot_data'][group][metric_name] = {
'age_centers': age_centers.copy(),
'age_labels': self.config.AGE_LABELS.copy(),
'male_means': [],
'female_means': [],
'male_stds': [],
'female_stds': [],
'male_counts': [],
'female_counts': [],
'combined_means': [],
'male_contributions': [],
'female_contributions': [],
'male_medians': [],
'female_medians': [],
'male_iqrs': [],
'female_iqrs': []
}
# Process each age group
for age_idx, age_label in enumerate(self.config.AGE_LABELS):
age_group_data = group_data[group_data['AgeGroup'] == age_label]
# Separate by gender
male_data = age_group_data[age_group_data['Gender'] == 'Male'][metric_info['column']]
female_data = age_group_data[age_group_data['Gender'] == 'Female'][metric_info['column']]
# Calculate statistics for each gender
male_stats = {
'count': len(male_data),
'mean': male_data.mean() if len(male_data) > 0 else np.nan,
'std': male_data.std() if len(male_data) > 0 else np.nan,
'min': male_data.min() if len(male_data) > 0 else np.nan,
'max': male_data.max() if len(male_data) > 0 else np.nan,
'median': male_data.median() if len(male_data) > 0 else np.nan,
'q25': male_data.quantile(0.25) if len(male_data) > 0 else np.nan,
'q75': male_data.quantile(0.75) if len(male_data) > 0 else np.nan
}
female_stats = {
'count': len(female_data),
'mean': female_data.mean() if len(female_data) > 0 else np.nan,
'std': female_data.std() if len(female_data) > 0 else np.nan,
'min': female_data.min() if len(female_data) > 0 else np.nan,
'max': female_data.max() if len(female_data) > 0 else np.nan,
'median': female_data.median() if len(female_data) > 0 else np.nan,
'q25': female_data.quantile(0.25) if len(female_data) > 0 else np.nan,
'q75': female_data.quantile(0.75) if len(female_data) > 0 else np.nan
}
# Store detailed statistics
table_data['detailed_stats'][group][metric_name][age_label] = {
'Male': male_stats,
'Female': female_stats
}
# Calculate values for plotting (handle NaN values)
male_mean = male_stats['mean'] if not np.isnan(male_stats['mean']) else 0
female_mean = female_stats['mean'] if not np.isnan(female_stats['mean']) else 0
male_std = male_stats['std'] if not np.isnan(male_stats['std']) else 0
female_std = female_stats['std'] if not np.isnan(female_stats['std']) else 0
male_median = male_stats['median'] if not np.isnan(male_stats['median']) else 0
female_median = female_stats['median'] if not np.isnan(female_stats['median']) else 0
# Calculate IQR for plotting (if needed for error bars)
male_iqr = (male_stats['q75'] - male_stats['q25']) if (
not np.isnan(male_stats['q75']) and not np.isnan(male_stats['q25'])) else 0
female_iqr = (female_stats['q75'] - female_stats['q25']) if (
not np.isnan(female_stats['q75']) and not np.isnan(female_stats['q25'])) else 0
# Calculate weighted combined mean and contributions
total_male_sum = male_stats['count'] * male_mean if male_stats['count'] > 0 else 0
total_female_sum = female_stats['count'] * female_mean if female_stats['count'] > 0 else 0
total_subjects = male_stats['count'] + female_stats['count']
if total_subjects > 0:
# True weighted combined mean across both genders
combined_mean = (total_male_sum + total_female_sum) / total_subjects
# Calculate proportional contributions to the combined mean
total_sum = total_male_sum + total_female_sum
if total_sum > 0:
male_contribution = (total_male_sum / total_sum) * combined_mean
female_contribution = (total_female_sum / total_sum) * combined_mean
else:
male_contribution = 0
female_contribution = 0
else:
combined_mean = 0
male_contribution = 0
female_contribution = 0
# Store plot data
plot_data = table_data['plot_data'][group][metric_name]
plot_data['male_means'].append(male_mean)
plot_data['female_means'].append(female_mean)
plot_data['male_stds'].append(male_std)
plot_data['female_stds'].append(female_std)
plot_data['male_counts'].append(male_stats['count'])
plot_data['female_counts'].append(female_stats['count'])
plot_data['combined_means'].append(combined_mean)
plot_data['male_contributions'].append(male_contribution)
plot_data['female_contributions'].append(female_contribution)
plot_data['male_medians'].append(male_median)
plot_data['female_medians'].append(female_median)
plot_data['male_iqrs'].append(male_iqr)
plot_data['female_iqrs'].append(female_iqr)
# Create the plot
plot_data = table_data['plot_data'][group][metric_name]
# Plot male contribution (bottom layer)
ax.fill_between(age_centers, 0, plot_data['male_contributions'],
color=self.config.COLORS['male'], alpha=0.7, label='Male')
# Plot female contribution (top layer)
ax.fill_between(age_centers, plot_data['male_contributions'],
plot_data['combined_means'],
color=self.config.COLORS['female'], alpha=0.7, label='Female')
# Optional: Add median lines for comparison (uncomment if desired)
# ax.plot(age_centers, plot_data['male_medians'],
# color=self.config.COLORS['male'], linestyle='--', alpha=0.8, label='Male Median')
# ax.plot(age_centers, plot_data['female_medians'],
# color=self.config.COLORS['female'], linestyle='--', alpha=0.8, label='Female Median')
# Formatting
# Calculate panel letter (A, B, C, D)
panel_idx = group_idx * 3 + metric_idx
panel_letter = chr(65 + panel_idx) # 65 is ASCII for 'A'
ax.set_title(f'{panel_letter}. {group} - {metric_info["title_suffix"]}', fontsize=20, fontweight='bold')
# ax.set_title(f'{group} - {metric_info["title_suffix"]}', fontsize=14, fontweight='bold')
ax.set_xlabel('Age (years)', fontsize=20)
ax.set_ylabel(metric_info['ylabel'], fontsize=20)
ax.legend(fontsize=20)
ax.grid(True, alpha=0.3)
ax.set_xticks(age_centers)
ax.set_xticklabels(self.config.AGE_LABELS, fontsize=20)
# Set reasonable y-axis limits
max_value = max(plot_data['combined_means']) if plot_data['combined_means'] else 0
if max_value > 0:
ax.set_ylim(0, max_value * 1.1)
plt.tight_layout()
plt.savefig(os.path.join(config.OUTPUT_DIR, 'total_lesion_burden_analysis.png'),
dpi=self.config.DPI, bbox_inches='tight', facecolor='white')
# Generate comprehensive documentation
self._generate_lesion_analysis_documentation(table_data)
# Generate and save tables
self._generate_lesion_burden_tables(table_data)
# Statistical analysis for both metrics
print(f"\n{'=' * 50}")
print("STATISTICAL COMPARISONS (HC vs MS)")
print(f"{'=' * 50}")
# Analysis for WMH area (absolute values)
hc_area = self.data[self.data[self.config.COLUMNS['group']] == 'HC'][self.config.COLUMNS['total_wmh']]
ms_area = self.data[self.data[self.config.COLUMNS['group']] == 'MS'][self.config.COLUMNS['total_wmh']]
if len(hc_area) > 0 and len(ms_area) > 0:
# Use the standardized testing approach
area_test, area_test_type, area_normality_info, area_effect_size_info = self._assess_normality_and_choose_test(
hc_area, ms_area, "WMH Area"
)
print(f"\nWMH Area comparison:")
print(f"HC: N={len(hc_area)}, mean ± SD = {hc_area.mean():.2f} ± {hc_area.std():.2f} mm²")
print(f"MS: N={len(ms_area)}, mean ± SD = {ms_area.mean():.2f} ± {ms_area.std():.2f} mm²")
# Print normality test results
if 'group1_shapiro_p' in area_normality_info:
print(
f"Normality tests: HC p={area_normality_info['group1_shapiro_p']:.3f}, MS p={area_normality_info['group2_shapiro_p']:.3f}")
if area_test_type == "parametric":
print(f"Independent t-test: t-statistic = {area_test.statistic:.3f}, p-value = {area_test.pvalue}")
print(f"Cohen's d = {area_effect_size_info[0]:.3f}")
elif area_test_type == "non_parametric":
print(f"Mann-Whitney U test: U-statistic = {area_test.statistic:.3f}, p-value = {area_test.pvalue}")
print(f"Effect size (r) = {area_effect_size_info[0]:.3f}")
# Report medians for non-parametric
print(
f"HC: median [IQR] = {hc_area.median():.2f} [{hc_area.quantile(0.25):.2f}-{hc_area.quantile(0.75):.2f}] mm²")
print(
f"MS: median [IQR] = {ms_area.median():.2f} [{ms_area.quantile(0.25):.2f}-{ms_area.quantile(0.75):.2f}] mm²")
else:
print("Insufficient data for WMH area comparison")
area_test = None
area_test_type = None
area_normality_info = {}
area_effect_size_info = (None, None)
# Analysis for WMH ratio (normalized values)
hc_ratio = self.data[self.data[self.config.COLUMNS['group']] == 'HC']['WMHRatio']
ms_ratio = self.data[self.data[self.config.COLUMNS['group']] == 'MS']['WMHRatio']
if len(hc_ratio) > 0 and len(ms_ratio) > 0:
# Use the standardized testing approach
ratio_test, ratio_test_type, ratio_normality_info, ratio_effect_size_info = self._assess_normality_and_choose_test(
hc_ratio, ms_ratio, "WMH Ratio"
)
print(f"\nWMH Ratio comparison:")
print(f"HC: N={len(hc_ratio)}, mean ± SD = {hc_ratio.mean():.2f} ± {hc_ratio.std():.2f}%")
print(f"MS: N={len(ms_ratio)}, mean ± SD = {ms_ratio.mean():.2f} ± {ms_ratio.std():.2f}%")
# Print normality test results
if 'group1_shapiro_p' in ratio_normality_info:
print(
f"Normality tests: HC p={ratio_normality_info['group1_shapiro_p']:.3f}, MS p={ratio_normality_info['group2_shapiro_p']:.3f}")
if ratio_test_type == "parametric":
print(
f"Independent t-test: t-statistic = {ratio_test.statistic:.3f}, p-value = {ratio_test.pvalue}")
print(f"Cohen's d = {ratio_effect_size_info[0]:.3f}")
elif ratio_test_type == "non_parametric":
print(
f"Mann-Whitney U test: U-statistic = {ratio_test.statistic:.3f}, p-value = {ratio_test.pvalue}")
print(f"Effect size (r) = {ratio_effect_size_info[0]:.3f}")
# Report medians for non-parametric
print(
f"HC: median [IQR] = {hc_ratio.median():.2f} [{hc_ratio.quantile(0.25):.2f}-{hc_ratio.quantile(0.75):.2f}]%")
print(
f"MS: median [IQR] = {ms_ratio.median():.2f} [{ms_ratio.quantile(0.25):.2f}-{ms_ratio.quantile(0.75):.2f}]%")
else:
print("Insufficient data for WMH ratio comparison")
ratio_test = None
ratio_test_type = None
ratio_normality_info = {}
ratio_effect_size_info = (None, None)
# Store results (update the existing results storage)
self.results['lesion_burden'] = {
'area_analysis': {
'hc_stats': (hc_area.median(), hc_area.quantile(0.25), hc_area.quantile(0.75)) if len(
hc_area) > 0 else (np.nan, np.nan, np.nan),
'ms_stats': (ms_area.median(), ms_area.quantile(0.25), ms_area.quantile(0.75)) if len(
ms_area) > 0 else (np.nan, np.nan, np.nan),
'hc_mean_stats': (hc_area.mean(), hc_area.std()) if len(hc_area) > 0 else (np.nan, np.nan),
'ms_mean_stats': (ms_area.mean(), ms_area.std()) if len(ms_area) > 0 else (np.nan, np.nan),
'comparison': area_test,
'test_type': area_test_type,
'normality_info': area_normality_info,
'effect_size': area_effect_size_info
},
'ratio_analysis': {
'hc_stats': (hc_ratio.median(), hc_ratio.quantile(0.25), hc_ratio.quantile(0.75)) if len(
hc_ratio) > 0 else (np.nan, np.nan, np.nan),
'ms_stats': (ms_ratio.median(), ms_ratio.quantile(0.25), ms_ratio.quantile(0.75)) if len(
ms_ratio) > 0 else (np.nan, np.nan, np.nan),
'hc_mean_stats': (hc_ratio.mean(), hc_ratio.std()) if len(hc_ratio) > 0 else (np.nan, np.nan),
'ms_mean_stats': (ms_ratio.mean(), ms_ratio.std()) if len(ms_ratio) > 0 else (np.nan, np.nan),
'comparison': ratio_test,
'test_type': ratio_test_type,
'normality_info': ratio_normality_info,
'effect_size': ratio_effect_size_info
},
'table_data': table_data
}
return self.results['lesion_burden']
def _generate_lesion_analysis_documentation(self, table_data):
"""Generate comprehensive documentation explaining the lesion burden figure and analysis"""
from datetime import datetime
# Create comprehensive documentation
doc_content = f"""
TOTAL LESION BURDEN ANALYSIS - COMPREHENSIVE DOCUMENTATION
==========================================================
Generated on: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}
OVERVIEW
--------
This analysis examines white matter hyperintensity (WMH) lesion burden in the brain
comparing Healthy Controls (HC) and Multiple Sclerosis (MS) patients. The analysis
includes both absolute WMH area measurements and normalized WMH ratios, stratified
by age groups and gender.
FIGURE DESCRIPTION
------------------
The figure consists of a 2x2 subplot layout:
Layout Structure:
- Top row: Healthy Controls (HC)
- Bottom row: Multiple Sclerosis (MS) patients
- Left column: Absolute WMH Area (mm²)
- Right column: WMH Ratio (%)
Subplot Details:
1. Top-Left: HC WMH Area
2. Top-Right: HC WMH Ratio
3. Bottom-Left: MS WMH Area
4. Bottom-Right: MS WMH Ratio
VISUALIZATION METHOD
--------------------
Chart Type: Stacked Area Plot
- Each subplot uses stacked area charts to show gender-stratified data across age groups
- Male data (bottom layer): Fills from 0 to male mean value
- Female data (top layer): Fills from male mean to total (male + female) mean
- This visualization allows comparison of both absolute values and gender contributions
Color Scheme:
- Male data: {table_data['metadata']['colors'].get('male', 'Blue')} (alpha=0.7)
- Female data: {table_data['metadata']['colors'].get('female', 'Red')} (alpha=0.7)
AGE STRATIFICATION
------------------
Age Groups: {', '.join(table_data['metadata']['age_labels'])}
Age Bins: {table_data['metadata']['age_bins']}
Age Centers (for plotting): {[f'{center:.1f}' for center in table_data['metadata']['age_centers']]}
The analysis stratifies data across these age groups to examine age-related changes
in lesion burden for both groups and genders.
METRICS ANALYZED
----------------
1. WMH Area (mm²):
- Absolute measurement of total white matter hyperintensity volume
- Column: {table_data['metadata']['metrics']['area']['column']}
- Units: Square millimeters (mm²)
- Clinical significance: Larger values indicate greater lesion burden
2. WMH Ratio (%):
- Normalized measurement relative to total brain area/volume
- Column: WMHRatio
- Units: Percentage (%)
- Clinical significance: Controls for individual brain size differences
CLINICAL CONTEXT
----------------
White Matter Hyperintensities (WMH):
- Bright signal areas on T2-weighted and FLAIR MRI sequences
- Associated with small vessel disease, aging, and neurodegeneration
- In MS: May represent demyelination, inflammation, or tissue damage
- Age-related increase is normal but accelerated in pathological conditions
Expected Patterns:
- MS patients typically show higher WMH burden than healthy controls
- Age-related increase in WMH burden in both groups
- MS may show accelerated age-related progression
- Gender differences may exist in lesion development patterns
ADVANCED STATISTICAL APPROACH
------------------------------
The analysis employs a sophisticated statistical testing framework:
Normality Assessment:
- Shapiro-Wilk tests performed on both groups (HC and MS) for each metric
- Significance threshold: p < 0.05 indicates non-normal distribution
- Results inform choice between parametric and non-parametric tests
Statistical Test Selection:
1. If both groups pass normality: Independent t-test (parametric)
2. If one or both groups fail normality: Mann-Whitney U test (non-parametric)
Effect Size Calculations:
- Parametric tests: Cohen's d
* Small effect: d = 0.2
* Medium effect: d = 0.5
* Large effect: d = 0.8
- Non-parametric tests: Rank-biserial correlation (r ≈ Z/√N, large-sample approximation)
Note: termed "Effect size r" in earlier versions; corrected per Reviewer 2, Minor Comment 4.
* Small effect: r = 0.1
* Medium effect: r = 0.3
* Large effect: r = 0.5
For each combination of:
- Group (HC vs MS)
- Age group ({len(table_data['metadata']['age_labels'])} categories)
- Gender (Male vs Female)
- Metric (Area vs Ratio)
The following statistics are calculated:
- Sample size (N)
- Mean ± Standard Deviation
- Minimum and Maximum values
- Median values
- Interquartile ranges (for non-parametric reporting)
STATISTICAL OUTPUT INTERPRETATION
---------------------------------
Parametric Results (t-test):
- Reports: Mean ± SD for both groups
- Test statistic: t-value and degrees of freedom
- p-value for significance testing
- Cohen's d for effect size magnitude
Non-parametric Results (Mann-Whitney U):
- Reports: Mean ± SD AND Median [IQR] for both groups
- Test statistic: U-value (or equivalent Z-score)
- p-value for significance testing
- Rank-biserial correlation (r) for effect magnitude (Mann-Whitney U; r ≈ Z/√N)
Normality Test Results:
- Shapiro-Wilk p-values reported for each group
- p < 0.05 indicates significant departure from normality
- Informs test selection rationale
INTERPRETATION GUIDELINES
-------------------------
Stacked Area Plot Interpretation:
- Height of bottom layer = Male mean value
- Height of top layer = Female mean value
- Total height = Combined mean (male + female means)
- Wider gaps between age points indicate larger differences
- Steeper slopes indicate rapid changes with age
Clinical Significance Thresholds:
- Minimal lesion burden: < 500 mm² (approximate)
- Mild lesion burden: 500-5000 mm²
- Moderate lesion burden: 5000-15000 mm²
- Severe lesion burden: > 15000 mm²
(Note: These are approximate guidelines and may vary by study protocol)
Expected Clinical Patterns:
- HC group: Low baseline with gradual age-related increase
- MS group: Higher baseline with potentially steeper age-related progression
- Gender differences: May reflect hormonal or lifestyle factors
- Age acceleration: MS may show earlier onset of lesion accumulation
Statistical Significance Levels:
- p < 0.001: Highly significant (strong evidence)
- p < 0.01: Very significant (moderate to strong evidence)
- p < 0.05: Significant (sufficient evidence)
- p ≥ 0.05: Non-significant (insufficient evidence)
Effect Size Interpretation:
- Cohen's d or r values indicate practical significance
- Large effect sizes may be clinically meaningful even if p > 0.05
- Small p-values with small effect sizes may lack clinical relevance
DATA QUALITY CONSIDERATIONS
----------------------------
- Zero values in plots indicate no subjects in that age/gender combination
- Small sample sizes may lead to unstable mean estimates
- Standard deviations provide insight into data variability within groups
- WMH measurements are sensitive to MRI acquisition parameters
- Manual/automated segmentation differences may affect absolute values
- Ratios help normalize for technical and anatomical variations
- Normality violations automatically trigger non-parametric alternatives
STATISTICAL TESTING DETAILS
----------------------------
Overall group comparisons (HC vs MS) performed separately for each metric:
1. WMH Area Analysis:
- Automatic normality assessment using Shapiro-Wilk test
- Test selection based on normality results
- Effect size calculation appropriate to test type
- Comprehensive reporting of descriptive statistics
2. WMH Ratio Analysis:
- Independent statistical analysis from area measurements
- Same rigorous normality assessment and test selection
- Separate effect size calculations
- Controls for multiple testing considerations
Test Assumptions:
- Parametric tests: Normality, independence, homogeneity of variance
- Non-parametric tests: Independence, similar distributions
- Both assume random sampling from populations of interest
OUTPUT FILES GENERATED
-----------------------
1. Figure: total_lesion_burden_analysis.png
- 2x2 subplot layout with stacked area plots showing proportional contributions
- High resolution (DPI: {getattr(self.config, 'DPI', 300)})
- White background for publication quality
2. Enhanced Statistics Tables (CSV):
- lesion_area_detailed_stats.csv (includes all descriptive statistics)
- lesion_ratio_detailed_stats.csv (includes all descriptive statistics)
- lesion_statistical_comparisons.csv (NEW: comprehensive test results)
3. Plot Data Tables (CSV):
- lesion_area_plot_data.csv (means and sample sizes)
- lesion_ratio_plot_data.csv (means and sample sizes)
4. Contribution Analysis (CSV):
- lesion_area_contributions.csv (NEW: proportional contributions)
- lesion_ratio_contributions.csv (NEW: proportional contributions)
5. Statistical Results Summary (CSV):
- lesion_normality_results.csv (NEW: normality test outcomes)
- lesion_effect_sizes.csv (NEW: effect size calculations)
6. Stacked Area Values (CSV):
- lesion_area_stacked_data.csv
- lesion_ratio_stacked_data.csv
7. This Documentation:
- lesion_burden_analysis_documentation.txt
TECHNICAL SPECIFICATIONS
-------------------------
Figure Specifications:
- Size: 16" x 12" (width x height)
- DPI: {getattr(self.config, 'DPI', 300)}
- Background: White
- Font sizes: Title=14pt (bold), Axis labels=12pt
- Grid: Enabled with 30% transparency
- Legend: Enabled for each subplot
Data Processing:
- Missing data handled by excluding from calculations
- Zero values used when no subjects available in category
- Robust statistics (median/IQR) preferred over mean/SD
- Sample size weighting for statistical calculations
Statistical Libraries:
- scipy.stats: Shapiro-Wilk, t-test, Mann-Whitney U
- numpy: Mathematical operations and statistical functions
- pandas: Data manipulation and summary statistics
Quality Control:
- Automatic handling of missing values
- Robust error handling for edge cases
- Comprehensive logging of statistical decisions
- Validation of statistical assumptions
LIMITATIONS AND CONSIDERATIONS
------------------------------
1. Sample Size Variations:
- Some age/gender combinations may have small sample sizes
- Unbalanced groups may affect statistical power
- Non-parametric tests more robust to unequal sample sizes
2. Multiple Comparisons:
- Two separate statistical tests performed (area and ratio)
- Consider Bonferroni correction: α = 0.05/2 = 0.025
- Family-wise error rate may be inflated without correction
3. Age Grouping Effects:
- Discretized age groups may mask continuous age effects
- Age bin boundaries are predetermined and may not reflect natural breakpoints
- Consider continuous age modeling for more detailed analysis
4. Lesion Measurement Considerations:
- WMH detection depends on MRI sequence parameters
- Segmentation methods (manual vs automated) may introduce variability
- Small lesions may be missed due to resolution limitations
- Partial volume effects at tissue boundaries
5. Stacked Area Representation:
- Visual emphasis on gender differences may overshadow group differences
- Direct comparison between groups requires careful interpretation
- Mean values used may not reflect distribution skewness
6. Statistical Assumptions:
- Automatic normality testing with Shapiro-Wilk (sensitive to large samples)
- Independence assumption may be violated in related subjects
- Equal variances assumed for t-tests (consider Welch's t-test alternative)
RECOMMENDED FOLLOW-UP ANALYSES
------------------------------
1. Advanced Statistical Approaches:
- Age as continuous variable (linear/polynomial regression)
- Two-way ANOVA with interaction terms (group × gender × age)
- Mixed-effects models for correlated data
- Bootstrap confidence intervals for robust inference
2. Multiple Comparisons Corrections:
- Bonferroni correction for family-wise error control
- False Discovery Rate (FDR) control for exploratory analyses
- Planned comparisons vs. post-hoc testing strategies
3. Effect Size and Power Analysis:
- Post-hoc power calculations for observed effects
- Sample size calculations for future studies
- Confidence intervals around effect size estimates
4. Alternative Statistical Approaches:
- Bayesian analysis for probabilistic interpretation
- Permutation tests for distribution-free inference
- Robust statistical methods for outlier resistance
5. Longitudinal Analysis:
- If temporal data available, analyze lesion progression rates
- Mixed-effects models for individual trajectories
- Survival analysis for time to lesion threshold
- Correlation with clinical severity measures
- Longitudinal tracking of lesion changes
- Predictive modeling for disease progression
QUALITY ASSURANCE CHECKLIST
----------------------------
✓ Normality testing performed automatically
✓ Appropriate statistical test selected based on data properties
✓ Effect sizes calculated and reported
✓ Both parametric and non-parametric results available
✓ Comprehensive descriptive statistics provided
✓ Multiple output formats for different use cases
✓ Documentation includes interpretation guidelines
✓ Limitations and assumptions clearly stated
✓ Recommendations for follow-up analyses provided
QUALITY CONTROL RECOMMENDATIONS
--------------------------------
1. Data Validation:
- Check for implausible values (negative areas, ratios > 100%)
- Identify and investigate outliers
- Verify age group assignments
2. Technical Validation:
- Compare manual vs automated segmentation on subset
- Inter-rater reliability for manual segmentations
- Phantom studies for scanner consistency
3. Clinical Validation:
- Correlation with clinical disability measures
- Agreement with radiological assessment
- Validation against established biomarkers
CONTACT AND METHODOLOGY
-----------------------
This analysis was generated using an automated pipeline for lesion burden assessment
with enhanced statistical testing capabilities.
For questions about methodology or interpretation, refer to:
- Original research protocol and statistical analysis plan
- Relevant neuroimaging analysis guidelines
- Statistical consulting resources for complex designs
Analysis Pipeline Version: Enhanced Statistical Testing v2.0
Statistical Methods: Automatic normality assessment with adaptive test selection
Last Updated: {datetime.now().strftime('%Y-%m-%d')}
END OF DOCUMENTATION
====================
"""
# Save documentation to file
doc_filename = os.path.join(config.OUTPUT_DIR, 'lesion_burden_analysis_documentation.txt')
with open(doc_filename, 'w', encoding='utf-8') as f:
f.write(doc_content)
print(f"\n{'=' * 80}")
print("COMPREHENSIVE LESION BURDEN DOCUMENTATION GENERATED")
print(f"{'=' * 80}")
print(f"Documentation saved to: lesion_burden_analysis_documentation.txt")
print(f"File contains detailed explanation of figure, methods, and clinical interpretation.")
print(f"- Enhanced statistical testing methodology")
print(f"- Normality assessment and test selection")
print(f"- Effect size calculations and interpretation")
print(f"- Figure interpretation and clinical relevance")
def _generate_lesion_burden_tables(self, table_data):
"""Generate comprehensive tables from the lesion burden analysis"""
# Table 1: Enhanced detailed statistics by group, age, and gender
print(f"\n{'=' * 80}")
print("TABLE 1: DETAILED STATISTICS BY GROUP, AGE, AND GENDER")
print(f"{'=' * 80}")
for metric_name in ['area', 'ratio', 'ratio_skull']:
unit = 'mm²' if metric_name == 'area' else '%'
if metric_name == 'area':
metric_title = 'WMH Area'
elif metric_name == 'ratio':
metric_title = 'WMH Ratio'
else:
metric_title = 'WMH Ratio-Skull'
print(f"\n{metric_title} ({unit}):")
print("-" * 70)
# Create DataFrame for this metric
rows = []
for group in ['HC', 'MS']:
for age_label in self.config.AGE_LABELS:
for gender in ['Male', 'Female']:
stats = table_data['detailed_stats'][group][metric_name][age_label][gender]
rows.append({
'Group': group,
'Age Group': age_label,
'Gender': gender,
'N': stats['count'],
'Mean': f"{stats['mean']:.2f}" if not np.isnan(stats['mean']) else 'N/A',
'SD': f"{stats['std']:.2f}" if not np.isnan(stats['std']) else 'N/A',
'Min': f"{stats['min']:.2f}" if not np.isnan(stats['min']) else 'N/A',
'Max': f"{stats['max']:.2f}" if not np.isnan(stats['max']) else 'N/A',
'Median': f"{stats['median']:.2f}" if not np.isnan(stats['median']) else 'N/A',
'IQR_25': f"{np.nan:.2f}" if np.isnan(
stats['median']) else f"{stats['median'] - stats['std'] / 2:.2f}",
'IQR_75': f"{np.nan:.2f}" if np.isnan(
stats['median']) else f"{stats['median'] + stats['std'] / 2:.2f}"
})
df = pd.DataFrame(rows)
print(df.to_string(index=False))
# Save to CSV
filename = f'lesion_{metric_name}_detailed_stats.csv'
df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Enhanced table saved to: {filename}")
# Table 2: Statistical Comparisons (NEW - Enhanced with normality and effect size info)
print(f"\n{'=' * 80}")
print("TABLE 2: STATISTICAL COMPARISONS (HC vs MS) - ENHANCED")
print(f"{'=' * 80}")
# Extract statistical results from the analysis
if hasattr(self, 'results') and 'lesion_burden' in self.results:
lesion_results = self.results['lesion_burden']
comparison_rows = []
for metric_type in ['area_analysis', 'ratio_analysis']:
metric_name = 'Area' if metric_type == 'area_analysis' else 'Ratio'
unit = 'mm²' if metric_type == 'area_analysis' else '%'
analysis = lesion_results[metric_type]
# Extract basic statistics
hc_mean, hc_std = analysis['hc_stats']
ms_mean, ms_std = analysis['ms_stats']
# Extract test results
comparison = analysis['comparison']
test_type = analysis['test_type']
normality_info = analysis.get('normality_info', {})
effect_size_info = analysis.get('effect_size', (None, None))
row = {
'Metric': f'{metric_name} ({unit})',
'HC_Mean': f"{hc_mean:.2f}" if not np.isnan(hc_mean) else 'N/A',
'HC_SD': f"{hc_std:.2f}" if not np.isnan(hc_std) else 'N/A',
'MS_Mean': f"{ms_mean:.2f}" if not np.isnan(ms_mean) else 'N/A',
'MS_SD': f"{ms_std:.2f}" if not np.isnan(ms_std) else 'N/A',
'Test_Type': test_type if test_type else 'N/A',
'Test_Statistic': f"{comparison.statistic:.3f}" if comparison else 'N/A',
'P_Value': f"{comparison.pvalue}" if comparison else 'N/A',
'Effect_Size': f"{effect_size_info[0]:.3f}" if effect_size_info[0] is not None else 'N/A',
'Effect_Size_Type': 'Cohen_d' if test_type == 'parametric' else 'r',
'HC_Normality_p': f"{normality_info.get('group1_shapiro_p', np.nan):.3f}" if 'group1_shapiro_p' in normality_info else 'N/A',
'MS_Normality_p': f"{normality_info.get('group2_shapiro_p', np.nan):.3f}" if 'group2_shapiro_p' in normality_info else 'N/A',
'Normality_Passed': 'Yes' if test_type == 'parametric' else 'No' if test_type == 'non_parametric' else 'N/A'
}
comparison_rows.append(row)
comparison_df = pd.DataFrame(comparison_rows)
print("\nStatistical Comparison Results:")
print("-" * 100)
print(comparison_df.to_string(index=False))
# Save to CSV
comparison_df.to_csv(os.path.join(config.OUTPUT_DIR, 'lesion_statistical_comparisons.csv'),
index=False)
print(f"\nStatistical comparisons saved to: lesion_statistical_comparisons.csv")
# Table 3: Plot data (means used for visualization)
print(f"\n{'=' * 80}")
print("TABLE 3: PLOT DATA (MEANS BY AGE GROUP)")
print(f"{'=' * 80}")
for metric_name in ['area', 'ratio', 'ratio_skull']:
unit = 'mm²' if metric_name == 'area' else '%'
if metric_name == 'area':
metric_title = 'WMH Area'
elif metric_name == 'ratio':
metric_title = 'WMH Ratio'
else:
metric_title = 'WMH Ratio-Skull'
print(f"\n{metric_title} - Mean Values and Contributions Used in Plot ({unit}):")
print("-" * 70)
# Create enhanced plot data table
plot_rows = []
for i, age_label in enumerate(self.config.AGE_LABELS):
age_center = table_data['plot_data']['HC'][metric_name]['age_centers'][i]
# HC data
hc_male_mean = table_data['plot_data']['HC'][metric_name]['male_means'][i]
hc_female_mean = table_data['plot_data']['HC'][metric_name]['female_means'][i]
hc_combined = table_data['plot_data']['HC'][metric_name]['combined_means'][i]
hc_male_contrib = table_data['plot_data']['HC'][metric_name]['male_contributions'][i]
hc_female_contrib = table_data['plot_data']['HC'][metric_name]['female_contributions'][i]
# MS data
ms_male_mean = table_data['plot_data']['MS'][metric_name]['male_means'][i]
ms_female_mean = table_data['plot_data']['MS'][metric_name]['female_means'][i]
ms_combined = table_data['plot_data']['MS'][metric_name]['combined_means'][i]
ms_male_contrib = table_data['plot_data']['MS'][metric_name]['male_contributions'][i]
ms_female_contrib = table_data['plot_data']['MS'][metric_name]['female_contributions'][i]
row = {
'Age_Group': age_label,
'Age_Center': f"{age_center:.1f}",
'HC_Male_Mean': f"{hc_male_mean:.2f}",
'HC_Female_Mean': f"{hc_female_mean:.2f}",
'HC_Combined_Mean': f"{hc_combined:.2f}",
'HC_Male_Contribution': f"{hc_male_contrib:.2f}",
'HC_Female_Contribution': f"{hc_female_contrib:.2f}",
'MS_Male_Mean': f"{ms_male_mean:.2f}",
'MS_Female_Mean': f"{ms_female_mean:.2f}",
'MS_Combined_Mean': f"{ms_combined:.2f}",
'MS_Male_Contribution': f"{ms_male_contrib:.2f}",
'MS_Female_Contribution': f"{ms_female_contrib:.2f}",
'HC_Male_N': table_data['plot_data']['HC'][metric_name]['male_counts'][i],
'HC_Female_N': table_data['plot_data']['HC'][metric_name]['female_counts'][i],
'MS_Male_N': table_data['plot_data']['MS'][metric_name]['male_counts'][i],
'MS_Female_N': table_data['plot_data']['MS'][metric_name]['female_counts'][i]
}
plot_rows.append(row)
plot_df = pd.DataFrame(plot_rows)
print(plot_df.to_string(index=False))
# Save to CSV
filename = f'lesion_{metric_name}_plot_data.csv'
plot_df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Enhanced plot data saved to: {filename}")
# Table 4: Contribution Analysis (NEW)
print(f"\n{'=' * 80}")
print("TABLE 4: PROPORTIONAL CONTRIBUTION ANALYSIS")
print(f"{'=' * 80}")
for metric_name in ['area', 'ratio']:
unit = 'mm²' if metric_name == 'area' else '%'
metric_title = 'Lesion Area' if metric_name == 'area' else 'Lesion Ratio'
print(f"\n{metric_title} - Proportional Contributions ({unit}):")
print("-" * 70)
# Create contribution analysis table
contrib_rows = []
for group in ['HC', 'MS']:
for i, age_label in enumerate(self.config.AGE_LABELS):
male_contrib = table_data['plot_data'][group][metric_name]['male_contributions'][i]
female_contrib = table_data['plot_data'][group][metric_name]['female_contributions'][i]
combined_mean = table_data['plot_data'][group][metric_name]['combined_means'][i]
male_count = table_data['plot_data'][group][metric_name]['male_counts'][i]
female_count = table_data['plot_data'][group][metric_name]['female_counts'][i]
total_count = male_count + female_count
# Calculate proportions
if combined_mean > 0:
male_prop = (male_contrib / combined_mean) * 100 if combined_mean > 0 else 0
female_prop = (female_contrib / combined_mean) * 100 if combined_mean > 0 else 0
else:
male_prop = 0
female_prop = 0
sample_male_prop = (male_count / total_count) * 100 if total_count > 0 else 0
sample_female_prop = (female_count / total_count) * 100 if total_count > 0 else 0
row = {
'Group': group,
'Age_Group': age_label,
'Combined_Mean': f"{combined_mean:.2f}",
'Male_Contribution': f"{male_contrib:.2f}",
'Female_Contribution': f"{female_contrib:.2f}",
'Male_Prop_of_Mean': f"{male_prop:.1f}%",
'Female_Prop_of_Mean': f"{female_prop:.1f}%",
'Male_Sample_Prop': f"{sample_male_prop:.1f}%",
'Female_Sample_Prop': f"{sample_female_prop:.1f}%",
'Total_N': total_count
}
contrib_rows.append(row)
contrib_df = pd.DataFrame(contrib_rows)
print(contrib_df.to_string(index=False))
# Save to CSV
filename = f'lesion_{metric_name}_contributions.csv'
contrib_df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Contribution analysis saved to: {filename}")
# Table 5: Effect Size Interpretation (NEW)
if hasattr(self, 'results') and 'lesion_burden' in self.results:
print(f"\n{'=' * 80}")
print("TABLE 5: EFFECT SIZE INTERPRETATION GUIDE")
print(f"{'=' * 80}")
effect_size_rows = []
lesion_results = self.results['lesion_burden']
for metric_type in ['area_analysis', 'ratio_analysis']:
metric_name = 'Area' if metric_type == 'area_analysis' else 'Ratio'
analysis = lesion_results[metric_type]
effect_size_info = analysis.get('effect_size', (None, None))
test_type = analysis['test_type']
if effect_size_info[0] is not None:
effect_size = effect_size_info[0]
if test_type == 'parametric':
# Cohen's d interpretation
if abs(effect_size) < 0.2:
magnitude = "Negligible"
elif abs(effect_size) < 0.5:
magnitude = "Small"
elif abs(effect_size) < 0.8:
magnitude = "Medium"
else:
magnitude = "Large"
else:
# Effect size r interpretation
if abs(effect_size) < 0.1:
magnitude = "Negligible"
elif abs(effect_size) < 0.3:
magnitude = "Small"
elif abs(effect_size) < 0.5:
magnitude = "Medium"
else:
magnitude = "Large"
row = {
'Metric': metric_name,
'Effect_Size_Value': f"{effect_size:.3f}",
'Effect_Size_Type': "Cohen's d" if test_type == 'parametric' else "Rank-biserial correlation (r)",
'Magnitude': magnitude,
'Interpretation': f"{magnitude} effect size indicating {'substantial' if magnitude in ['Medium', 'Large'] else 'minimal'} practical difference"
}
effect_size_rows.append(row)
if effect_size_rows:
effect_df = pd.DataFrame(effect_size_rows)
print("\nEffect Size Interpretations:")
print("-" * 60)
print(effect_df.to_string(index=False))
# Save to CSV
effect_df.to_csv(os.path.join(config.OUTPUT_DIR, 'lesion_effect_sizes.csv'), index=False)
print(f"\nEffect size interpretations saved to: lesion_effect_sizes.csv")
# Table 6: Stacked area values (cumulative for visualization)
print(f"\n{'=' * 80}")
print("TABLE 6: STACKED AREA VALUES (FOR AREA CHART)")
print(f"{'=' * 80}")
for metric_name in ['area', 'ratio', 'ratio_skull']:
unit = 'mm²' if metric_name == 'area' else '%'
if metric_name == 'area':
metric_title = 'WMH Area'
elif metric_name == 'ratio':
metric_title = 'WMH Ratio'
else:
metric_title = 'WMH Ratio-Skull'
print(f"\n{metric_title} - Stacked Values ({unit}):")
print("-" * 70)
# Create stacked data table
stacked_rows = []
for group in ['HC', 'MS']:
for i, age_label in enumerate(self.config.AGE_LABELS):
male_mean = table_data['plot_data'][group][metric_name]['male_means'][i]
female_mean = table_data['plot_data'][group][metric_name]['female_means'][i]
row = {
'Group': group,
'Age Group': age_label,
'Male Layer (0 to Male)': f"0.00 to {male_mean:.2f}",
'Female Layer (Male to Total)': f"{male_mean:.2f} to {male_mean + female_mean:.2f}",
'Total Height': f"{male_mean + female_mean:.2f}",
'Male Contribution': f"{male_mean:.2f}",
'Female Contribution': f"{female_mean:.2f}"
}
stacked_rows.append(row)
stacked_df = pd.DataFrame(stacked_rows)
print(stacked_df.to_string(index=False))
# Save to CSV
filename = f'lesion_{metric_name}_stacked_data.csv'
stacked_df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Stacked data saved to: {filename}")
print(f"\n{'=' * 80}")
print("All lesion burden tables have been generated and saved to CSV files!")
print(f"{'=' * 80}")
print("\nGenerated Files Summary:")
print("- lesion_area_detailed_stats.csv (Complete descriptive statistics)")
print("- lesion_ratio_detailed_stats.csv (Complete descriptive statistics)")
print("- lesion_area_plot_data.csv (Data used for visualization)")
print("- lesion_ratio_plot_data.csv (Data used for visualization)")
print("- lesion_area_stacked_data.csv (Stacked area chart values)")
print("- lesion_ratio_stacked_data.csv (Stacked area chart values)")
print("- lesion_area_robust_stats.csv (Non-parametric robust measures)")
print("- lesion_ratio_robust_stats.csv (Non-parametric robust measures)")
print("- lesion_burden_analysis_documentation.txt (Comprehensive documentation)")
print("ENHANCED STATISTICAL TABLES GENERATED!")
print("- total_lesion_burden_analysis.png (Figure)")
print("• Enhanced descriptive statistics with IQR")
print("• Comprehensive statistical comparison results")
print("• Detailed plot data with contributions")
print("• Proportional contribution analysis")
print("• Effect size interpretations")
print("• All tables saved as CSV files for further analysis")
def assess_normality(self, data, variable_name="", group_name="", alpha=0.05):
"""
Comprehensive normality assessment with multiple criteria
"""
import scipy.stats as stats
import numpy as np
if len(data) < 3:
return False, {"reason": "Insufficient data", "n": len(data)}
# Remove NaN values
clean_data = data.dropna() if hasattr(data, 'dropna') else data[~np.isnan(data)]
if len(clean_data) < 3:
return False, {"reason": "Insufficient valid data after removing NaN", "n": len(clean_data)}
# Shapiro-Wilk test
shapiro_stat, shapiro_p = stats.shapiro(clean_data)
# Skewness and kurtosis
skewness = stats.skew(clean_data)
kurtosis_val = stats.kurtosis(clean_data)
# Multiple criteria for normality
shapiro_normal = shapiro_p > alpha
skew_normal = abs(skewness) < 2
kurtosis_normal = abs(kurtosis_val) < 7
is_normal = shapiro_normal and skew_normal and kurtosis_normal
assessment = {
'shapiro_statistic': shapiro_stat,
'shapiro_p': shapiro_p,
'shapiro_normal': shapiro_normal,
'skewness': skewness,
'skew_normal': skew_normal,
'kurtosis': kurtosis_val,
'kurtosis_normal': kurtosis_normal,
'n_samples': len(clean_data),
'variable': variable_name,
'group': group_name
}
return is_normal, assessment
def standardized_group_comparison(self, group1_data, group2_data, group1_name="Group1", group2_name="Group2",
variable_name=""):
"""
Standardized approach for group comparisons with automatic test selection
"""
import scipy.stats as stats
import numpy as np
# Clean data
g1_clean = group1_data.dropna() if hasattr(group1_data, 'dropna') else group1_data[~np.isnan(group1_data)]
g2_clean = group2_data.dropna() if hasattr(group2_data, 'dropna') else group2_data[~np.isnan(group2_data)]
if len(g1_clean) < 2 or len(g2_clean) < 2:
return None, {"error": "Insufficient data for comparison"}
# Test normality for both groups
g1_normal, g1_assessment = self.assess_normality(g1_clean, variable_name, group1_name)
g2_normal, g2_assessment = self.assess_normality(g2_clean, variable_name, group2_name)
# Choose appropriate test
if g1_normal and g2_normal:
# Use parametric test
test_stat, p_value = stats.ttest_ind(g1_clean, g2_clean)
# Calculate Cohen's d
pooled_std = np.sqrt(((len(g1_clean) - 1) * np.var(g1_clean, ddof=1) +
(len(g2_clean) - 1) * np.var(g2_clean, ddof=1)) /
(len(g1_clean) + len(g2_clean) - 2))
cohens_d = (np.mean(g1_clean) - np.mean(g2_clean)) / pooled_std
result = {
'test_type': 'Independent t-test',
'test_statistic': test_stat,
'p_value': p_value,
'effect_size_type': "Cohen's d",
'effect_size': cohens_d,
'parametric': True,
'group1_stats': {
'mean': np.mean(g1_clean),
'std': np.std(g1_clean, ddof=1),
'n': len(g1_clean)
},
'group2_stats': {
'mean': np.mean(g2_clean),
'std': np.std(g2_clean, ddof=1),
'n': len(g2_clean)
}
}
else:
# Use non-parametric test
test_stat, p_value = stats.mannwhitneyu(g1_clean, g2_clean, alternative='two-sided')
# Rank-biserial correlation for Mann-Whitney U (Reviewer 2, Minor Comment 4).
# r ≈ Z / √N — large-sample approximation of rank-biserial correlation.
n_total = len(g1_clean) + len(g2_clean)
z_score = stats.norm.ppf(1 - p_value / 2) if p_value > 0 else 0
effect_size_r = abs(z_score) / np.sqrt(n_total)
result = {
'test_type': 'Mann-Whitney U test',
'test_statistic': test_stat,
'p_value': p_value,
'effect_size_type': 'Rank-biserial correlation (r)',
'effect_size': effect_size_r,
'parametric': False,
'group1_stats': {
'median': np.median(g1_clean),
'q25': np.percentile(g1_clean, 25),
'q75': np.percentile(g1_clean, 75),
'n': len(g1_clean)
},
'group2_stats': {
'median': np.median(g2_clean),
'q25': np.percentile(g2_clean, 25),
'q75': np.percentile(g2_clean, 75),
'n': len(g2_clean)
}
}
# Add normality assessment details
result['normality_assessment'] = {
'group1': g1_assessment,
'group2': g2_assessment,
'both_normal': g1_normal and g2_normal
}
return result, None
def ms_subgroup_analysis(self):
"""Analyze MS lesion subtypes with detailed stratification - both area and ratio"""
print("\n" + "=" * 60)
print("MS SUBGROUP LESION ANALYSIS")
print("=" * 60)
# Filter MS patients only
ms_data = self.data[self.data[self.config.COLUMNS['group']] == 'MS'].copy()
# Create 3x2 subplot layout
fig, axes = plt.subplots(3, 2, figsize=(16, 18))
fig.patch.set_facecolor('white') # Ensure white background
age_centers = [np.mean(age_range) for age_range in self.config.AGE_BINS]
# Analysis for: All MS, Female MS, Male MS
subgroups = [
('All MS Patients', ms_data),
('Female MS Patients', ms_data[ms_data['Gender'] == 'Female']),
('Male MS Patients', ms_data[ms_data['Gender'] == 'Male'])
]
# Metrics to plot: absolute area (left column) and normalized ratio (right column)
metrics = {
'area': {
'columns': {
'pewmh': self.config.COLUMNS.get('peri_wmh', 'peri_wmh'), # Absolute area columns
'dwmh': self.config.COLUMNS.get('deep_wmh', 'deep_wmh'),
'jcwmh': self.config.COLUMNS.get('juxta_wmh', 'juxta_wmh')
},
'ylabel': 'WMH Subtype Area (mm²)',
'title_suffix': 'WMH Subtype Areas'
},
'ratio': {
'columns': {
'pewmh': 'peri_wmh_ratio', # Ratio columns
'dwmh': 'deep_wmh_ratio',
'jcwmh': 'juxta_wmh_ratio'
},
'ylabel': 'WMH Subtype Ratio (%)',
'title_suffix': 'WMH Subtype Ratios'
}
}
# Dictionary to store all table data
table_data = {
'detailed_stats': {},
'plot_data': {},
'gender_comparisons': {},
'metadata': {
'age_bins': self.config.AGE_BINS,
'age_labels': self.config.AGE_LABELS,
'age_centers': age_centers,
'subgroups': [name for name, _ in subgroups],
'metrics': metrics,
'colors': self.config.COLORS,
'lesion_subtypes': ['PEWMH', 'DWMH', 'JCWMH']
}
}
for row_idx, (group_title, data) in enumerate(subgroups):
# Initialize group data in tables
table_data['detailed_stats'][group_title] = {}
table_data['plot_data'][group_title] = {}
for col_idx, (metric_type, metric_info) in enumerate(
[('area', metrics['area']), ('ratio', metrics['ratio'])]):
ax = axes[row_idx, col_idx]
# Initialize metric data in tables
table_data['detailed_stats'][group_title][metric_type] = {}
table_data['plot_data'][group_title][metric_type] = {
'age_centers': age_centers.copy(),
'age_labels': self.config.AGE_LABELS.copy(),
'pewmh_means': [],
'dwmh_means': [],
'jcwmh_means': [],
'pewmh_stds': [],
'dwmh_stds': [],
'jcwmh_stds': [],
'pewmh_counts': [],
'dwmh_counts': [],
'jcwmh_counts': []
}
# Prepare data for three-layer stacked area plot
pewmh_means = []
dwmh_means = []
jcwmh_means = []
pewmh_stds = []
dwmh_stds = []
jcwmh_stds = []
pewmh_counts = []
dwmh_counts = []
jcwmh_counts = []
for age_idx, age_label in enumerate(self.config.AGE_LABELS):
age_group_data = data[data['AgeGroup'] == age_label]
# Calculate statistics for each lesion subtype
pewmh_data = age_group_data[metric_info['columns']['pewmh']].dropna()
dwmh_data = age_group_data[metric_info['columns']['dwmh']].dropna()
jcwmh_data = age_group_data[metric_info['columns']['jcwmh']].dropna()
# Means for plotting
pewmh_mean = pewmh_data.mean() if len(pewmh_data) > 0 else 0
dwmh_mean = dwmh_data.mean() if len(dwmh_data) > 0 else 0
jcwmh_mean = jcwmh_data.mean() if len(jcwmh_data) > 0 else 0
# Standard deviations
pewmh_std = pewmh_data.std() if len(pewmh_data) > 0 else 0
dwmh_std = dwmh_data.std() if len(dwmh_data) > 0 else 0
jcwmh_std = jcwmh_data.std() if len(jcwmh_data) > 0 else 0
# Sample counts
pewmh_count = len(pewmh_data)
dwmh_count = len(dwmh_data)
jcwmh_count = len(jcwmh_data)
# Store for plotting
pewmh_means.append(pewmh_mean)
dwmh_means.append(dwmh_mean)
jcwmh_means.append(jcwmh_mean)
pewmh_stds.append(pewmh_std)
dwmh_stds.append(dwmh_std)
jcwmh_stds.append(jcwmh_std)
pewmh_counts.append(pewmh_count)
dwmh_counts.append(dwmh_count)
jcwmh_counts.append(jcwmh_count)
# Store detailed statistics for tables
if age_label not in table_data['detailed_stats'][group_title][metric_type]:
table_data['detailed_stats'][group_title][metric_type][age_label] = {}
for lesion_type, lesion_data in [('PEWMH', pewmh_data), ('DWMH', dwmh_data),
('JCWMH', jcwmh_data)]:
table_data['detailed_stats'][group_title][metric_type][age_label][lesion_type] = {
'count': len(lesion_data),
'mean': lesion_data.mean() if len(lesion_data) > 0 else np.nan,
'std': lesion_data.std() if len(lesion_data) > 0 else np.nan,
'min': lesion_data.min() if len(lesion_data) > 0 else np.nan,
'max': lesion_data.max() if len(lesion_data) > 0 else np.nan,
'median': lesion_data.median() if len(lesion_data) > 0 else np.nan,
'q25': lesion_data.quantile(0.25) if len(lesion_data) > 0 else np.nan,
'q75': lesion_data.quantile(0.75) if len(lesion_data) > 0 else np.nan
}
# Store plot data
table_data['plot_data'][group_title][metric_type]['pewmh_means'] = pewmh_means
table_data['plot_data'][group_title][metric_type]['dwmh_means'] = dwmh_means
table_data['plot_data'][group_title][metric_type]['jcwmh_means'] = jcwmh_means
table_data['plot_data'][group_title][metric_type]['pewmh_stds'] = pewmh_stds
table_data['plot_data'][group_title][metric_type]['dwmh_stds'] = dwmh_stds
table_data['plot_data'][group_title][metric_type]['jcwmh_stds'] = jcwmh_stds
table_data['plot_data'][group_title][metric_type]['pewmh_counts'] = pewmh_counts
table_data['plot_data'][group_title][metric_type]['dwmh_counts'] = dwmh_counts
table_data['plot_data'][group_title][metric_type]['jcwmh_counts'] = jcwmh_counts
# Create three-layer stacked area plot
ax.fill_between(age_centers, 0, pewmh_means,
color=self.config.COLORS['pewmh'], alpha=0.8, label='PEWMH')
ax.fill_between(age_centers, pewmh_means,
np.array(pewmh_means) + np.array(dwmh_means),
color=self.config.COLORS['dwmh'], alpha=0.8, label='DWMH')
ax.fill_between(age_centers, np.array(pewmh_means) + np.array(dwmh_means),
np.array(pewmh_means) + np.array(dwmh_means) + np.array(jcwmh_means),
color=self.config.COLORS['jcwmh'], alpha=0.8, label='JCWMH')
# Formatting
# Calculate panel letter (A through F for 3x2 layout)
panel_idx = row_idx * 2 + col_idx
panel_letter = chr(65 + panel_idx) # 65 is ASCII for 'A'
ax.set_title(f'{panel_letter}. {group_title} - {metric_info["title_suffix"]}', fontsize=18, fontweight='bold')
# ax.set_title(f'{group_title} - {metric_info["title_suffix"]}', fontsize=14, fontweight='bold')
ax.set_xlabel('Age (years)', fontsize=16)
ax.set_ylabel(metric_info['ylabel'], fontsize=16)
ax.legend(loc='upper right', fontsize=15) #, frameon=True, fancybox=True, shadow=True)
ax.grid(True, alpha=0.3)
ax.set_xticks(age_centers)
ax.set_xticklabels(self.config.AGE_LABELS, fontsize=16)
plt.tight_layout()
plt.savefig(os.path.join(config.OUTPUT_DIR, 'ms_subgroup_analysis.png'),
dpi=self.config.DPI, bbox_inches='tight', facecolor='white')
# Generate comprehensive documentation
self._generate_ms_subgroup_documentation(table_data)
# Generate and save tables
self._generate_ms_subgroup_tables(table_data)
# Statistical analysis for both area and ratio metrics
print(f"\n{'=' * 50}")
print("MS LESION SUBTYPE STATISTICS")
print(f"{'=' * 50}")
# Analysis for absolute areas
print(f"\nMS Lesion Subtype Areas (mm²):")
area_columns = [
('PEWMH', metrics['area']['columns']['pewmh']),
('DWMH', metrics['area']['columns']['dwmh']),
('JCWMH', metrics['area']['columns']['jcwmh'])
]
for subtype_name, column in area_columns:
if column in ms_data.columns:
subtype_data = ms_data[column].dropna()
if len(subtype_data) > 0:
print(
f"{subtype_name}: median [IQR] = {subtype_data.median():.2f} [{subtype_data.quantile(0.25):.2f}-{subtype_data.quantile(0.75):.2f}] mm²")
else:
print(f"{subtype_name}: No valid data available")
else:
print(f"{subtype_name}: Column '{column}' not found in data")
# Analysis for ratios
print(f"\nMS Lesion Subtype Ratios (%):")
ratio_columns = [
('PEWMH', 'peri_wmh_ratio'),
('DWMH', 'deep_wmh_ratio'),
('JCWMH', 'juxta_wmh_ratio')
]
for subtype_name, column in ratio_columns:
subtype_data = ms_data[column].dropna()
print(
f"{subtype_name}: median [IQR] = {subtype_data.median():.3f} [{subtype_data.quantile(0.25):.3f}-{subtype_data.quantile(0.75):.3f}]%")
print(f"\n{'=' * 50}")
print("GENDER COMPARISON WITH STANDARDIZED STATISTICAL APPROACH")
print(f"{'=' * 50}")
gender_comparisons = {}
# Area comparisons with standardized approach
print(f"\nArea Comparisons (mm²) - Standardized Statistical Testing:")
for subtype_name, column in area_columns:
if column in ms_data.columns:
male_data = ms_data[ms_data['Gender'] == 'Male'][column].dropna()
female_data = ms_data[ms_data['Gender'] == 'Female'][column].dropna()
if len(male_data) > 0 and len(female_data) > 0:
# Use standardized comparison
comparison_result, error = self.standardized_group_comparison(
male_data, female_data, "Male", "Female", f"{subtype_name} Area"
)
if comparison_result:
print(f"\n{subtype_name} - {comparison_result['test_type']}:")
if comparison_result['parametric']:
print(
f" Male mean ± SD: {comparison_result['group1_stats']['mean']:.2f} ± {comparison_result['group1_stats']['std']:.2f} mm² (n={comparison_result['group1_stats']['n']})")
print(
f" Female mean ± SD: {comparison_result['group2_stats']['mean']:.2f} ± {comparison_result['group2_stats']['std']:.2f} mm² (n={comparison_result['group2_stats']['n']})")
else:
print(
f" Male median [IQR]: {comparison_result['group1_stats']['median']:.2f} [{comparison_result['group1_stats']['q25']:.2f}-{comparison_result['group1_stats']['q75']:.2f}] mm² (n={comparison_result['group1_stats']['n']})")
print(
f" Female median [IQR]: {comparison_result['group2_stats']['median']:.2f} [{comparison_result['group2_stats']['q25']:.2f}-{comparison_result['group2_stats']['q75']:.2f}] mm² (n={comparison_result['group2_stats']['n']})")
print(f" Test statistic: {comparison_result['test_statistic']:.3f}")
print(f" P-value: {comparison_result['p_value']}")
print(f" {comparison_result['effect_size_type']}: {comparison_result['effect_size']:.3f}")
# Normality test results
g1_normal = comparison_result['normality_assessment']['group1']['shapiro_normal']
g2_normal = comparison_result['normality_assessment']['group2']['shapiro_normal']
print(
f" Normality: Male p={comparison_result['normality_assessment']['group1']['shapiro_p']:.3f} ({'Normal' if g1_normal else 'Non-normal'}), "
f"Female p={comparison_result['normality_assessment']['group2']['shapiro_p']:.3f} ({'Normal' if g2_normal else 'Non-normal'})")
if subtype_name not in gender_comparisons:
gender_comparisons[subtype_name] = {}
gender_comparisons[subtype_name]['area'] = comparison_result
else:
print(f"{subtype_name}: {error}")
# Ratio comparisons with standardized approach
print(f"\nRatio Comparisons (%) - Standardized Statistical Testing:")
for subtype_name, column in ratio_columns:
male_data = ms_data[ms_data['Gender'] == 'Male'][column].dropna()
female_data = ms_data[ms_data['Gender'] == 'Female'][column].dropna()
if len(male_data) > 0 and len(female_data) > 0:
# Use standardized comparison
comparison_result, error = self.standardized_group_comparison(
male_data, female_data, "Male", "Female", f"{subtype_name} Ratio"
)
if comparison_result:
print(f"\n{subtype_name} - {comparison_result['test_type']}:")
if comparison_result['parametric']:
print(
f" Male mean ± SD: {comparison_result['group1_stats']['mean']:.3f} ± {comparison_result['group1_stats']['std']:.3f}% (n={comparison_result['group1_stats']['n']})")
print(
f" Female mean ± SD: {comparison_result['group2_stats']['mean']:.3f} ± {comparison_result['group2_stats']['std']:.3f}% (n={comparison_result['group2_stats']['n']})")
else:
print(
f" Male median [IQR]: {comparison_result['group1_stats']['median']:.3f} [{comparison_result['group1_stats']['q25']:.3f}-{comparison_result['group1_stats']['q75']:.3f}]% (n={comparison_result['group1_stats']['n']})")
print(
f" Female median [IQR]: {comparison_result['group2_stats']['median']:.3f} [{comparison_result['group2_stats']['q25']:.3f}-{comparison_result['group2_stats']['q75']:.3f}]% (n={comparison_result['group2_stats']['n']})")
print(f" Test statistic: {comparison_result['test_statistic']:.3f}")
print(f" P-value: {comparison_result['p_value']}")
print(f" {comparison_result['effect_size_type']}: {comparison_result['effect_size']:.3f}")
# Normality test results
g1_normal = comparison_result['normality_assessment']['group1']['shapiro_normal']
g2_normal = comparison_result['normality_assessment']['group2']['shapiro_normal']
print(
f" Normality: Male p={comparison_result['normality_assessment']['group1']['shapiro_p']:.3f} ({'Normal' if g1_normal else 'Non-normal'}), "
f"Female p={comparison_result['normality_assessment']['group2']['shapiro_p']:.3f} ({'Normal' if g2_normal else 'Non-normal'})")
if subtype_name not in gender_comparisons:
gender_comparisons[subtype_name] = {}
gender_comparisons[subtype_name]['ratio'] = comparison_result
else:
print(f"{subtype_name}: {error}")
# Store gender comparison data in table_data
table_data['gender_comparisons'] = gender_comparisons
# Store comprehensive results
self.results['ms_subtypes'] = {
'area_analysis': {
subtype: {
'all_stats': (ms_data[column].dropna().median(),
ms_data[column].dropna().quantile(0.25),
ms_data[column].dropna().quantile(0.75)) if column in ms_data.columns and len(
ms_data[column].dropna()) > 0 else None
} for subtype, column in area_columns
},
'ratio_analysis': {
subtype: {
'all_stats': (ms_data[column].dropna().median(),
ms_data[column].dropna().quantile(0.25),
ms_data[column].dropna().quantile(0.75))
} for subtype, column in ratio_columns
},
'gender_comparison': gender_comparisons,
'table_data': table_data # Add table data to results
}
return self.results['ms_subtypes']
def _generate_ms_subgroup_tables(self, table_data):
"""Generate comprehensive tables from the MS subgroup analysis"""
# Table 1: Detailed statistics by subgroup, age, and lesion subtype
print(f"\n{'=' * 80}")
print("TABLE 1: DETAILED STATISTICS BY SUBGROUP, AGE, AND LESION SUBTYPE")
print(f"{'=' * 80}")
for metric_name in ['area', 'ratio']:
unit = 'mm²' if metric_name == 'area' else '%'
metric_title = 'WMH Subtype Area' if metric_name == 'area' else 'WMH Subtype Ratio'
print(f"\n{metric_title} ({unit}):")
print("-" * 80)
# Create DataFrame for this metric
rows = []
for subgroup in table_data['metadata']['subgroups']:
for age_label in self.config.AGE_LABELS:
for lesion_type in ['PEWMH', 'DWMH', 'JCWMH']:
if (subgroup in table_data['detailed_stats'] and
metric_name in table_data['detailed_stats'][subgroup] and
age_label in table_data['detailed_stats'][subgroup][metric_name] and
lesion_type in table_data['detailed_stats'][subgroup][metric_name][age_label]):
stats = table_data['detailed_stats'][subgroup][metric_name][age_label][lesion_type]
rows.append({
'Subgroup': subgroup,
'Age Group': age_label,
'Lesion Type': lesion_type,
'N': stats['count'],
'Mean': f"{stats['mean']:.2f}" if not np.isnan(stats['mean']) else 'N/A',
'SD': f"{stats['std']:.2f}" if not np.isnan(stats['std']) else 'N/A',
'Median': f"{stats['median']:.2f}" if not np.isnan(stats['median']) else 'N/A',
'Q25': f"{stats['q25']:.2f}" if not np.isnan(stats['q25']) else 'N/A',
'Q75': f"{stats['q75']:.2f}" if not np.isnan(stats['q75']) else 'N/A',
'Min': f"{stats['min']:.2f}" if not np.isnan(stats['min']) else 'N/A',
'Max': f"{stats['max']:.2f}" if not np.isnan(stats['max']) else 'N/A'
})
if rows:
df = pd.DataFrame(rows)
print(df.to_string(index=False))
# Save to CSV
filename = f'ms_subgroup_{metric_name}_detailed_stats.csv'
df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Table saved to: {filename}")
else:
print("No data available for detailed statistics table.")
# Table 2: Plot data (means used for visualization)
print(f"\n{'=' * 80}")
print("TABLE 2: PLOT DATA (MEANS BY AGE GROUP AND LESION SUBTYPE)")
print(f"{'=' * 80}")
for metric_name in ['area', 'ratio']:
unit = 'mm²' if metric_name == 'area' else '%'
metric_title = 'WMH Subtype Area' if metric_name == 'area' else 'WMH Subtype Ratio'
print(f"\n{metric_title} - Mean Values Used in Plot ({unit}):")
print("-" * 90)
# Create plot data table
plot_rows = []
for subgroup in table_data['metadata']['subgroups']:
if (subgroup in table_data['plot_data'] and
metric_name in table_data['plot_data'][subgroup]):
plot_data = table_data['plot_data'][subgroup][metric_name]
for i, age_label in enumerate(self.config.AGE_LABELS):
if i < len(plot_data['age_centers']):
age_center = plot_data['age_centers'][i]
row = {
'Subgroup': subgroup,
'Age Group': age_label,
'Age Center': f"{age_center:.1f}",
'PEWMH Mean': f"{plot_data['pewmh_means'][i]:.2f}",
'DWMH Mean': f"{plot_data['dwmh_means'][i]:.2f}",
'JCWMH Mean': f"{plot_data['jcwmh_means'][i]:.2f}",
'PEWMH N': plot_data['pewmh_counts'][i],
'DWMH N': plot_data['dwmh_counts'][i],
'JCWMH N': plot_data['jcwmh_counts'][i],
'Total Mean': f"{plot_data['pewmh_means'][i] + plot_data['dwmh_means'][i] + plot_data['jcwmh_means'][i]:.2f}"
}
plot_rows.append(row)
if plot_rows:
plot_df = pd.DataFrame(plot_rows)
print(plot_df.to_string(index=False))
# Save to CSV
filename = f'ms_subgroup_{metric_name}_plot_data.csv'
plot_df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Plot data saved to: {filename}")
else:
print("No data available for plot data table.")
# Table 3: Three-layer stacked area values
print(f"\n{'=' * 80}")
print("TABLE 3: THREE-LAYER STACKED AREA VALUES")
print(f"{'=' * 80}")
for metric_name in ['area', 'ratio']:
unit = 'mm²' if metric_name == 'area' else '%'
metric_title = 'WMH Subtype Area' if metric_name == 'area' else 'WMH Subtype Ratio'
print(f"\n{metric_title} - Stacked Layer Values ({unit}):")
print("-" * 100)
# Create stacked data table
stacked_rows = []
for subgroup in table_data['metadata']['subgroups']:
if (subgroup in table_data['plot_data'] and
metric_name in table_data['plot_data'][subgroup]):
plot_data = table_data['plot_data'][subgroup][metric_name]
for i, age_label in enumerate(self.config.AGE_LABELS):
if i < len(plot_data['pewmh_means']):
pewmh_mean = plot_data['pewmh_means'][i]
dwmh_mean = plot_data['dwmh_means'][i]
jcwmh_mean = plot_data['jcwmh_means'][i]
# Calculate cumulative layer boundaries
layer1_end = pewmh_mean
layer2_end = pewmh_mean + dwmh_mean
layer3_end = pewmh_mean + dwmh_mean + jcwmh_mean
row = {
'Subgroup': subgroup,
'Age Group': age_label,
'PEWMH Layer': f"0.00 to {layer1_end:.2f}",
'DWMH Layer': f"{layer1_end:.2f} to {layer2_end:.2f}",
'JCWMH Layer': f"{layer2_end:.2f} to {layer3_end:.2f}",
'Total Height': f"{layer3_end:.2f}",
'PEWMH Contribution': f"{pewmh_mean:.2f}",
'DWMH Contribution': f"{dwmh_mean:.2f}",
'JCWMH Contribution': f"{jcwmh_mean:.2f}",
'PEWMH %': f"{(pewmh_mean / layer3_end * 100):.1f}%" if layer3_end > 0 else "0.0%",
'DWMH %': f"{(dwmh_mean / layer3_end * 100):.1f}%" if layer3_end > 0 else "0.0%",
'JCWMH %': f"{(jcwmh_mean / layer3_end * 100):.1f}%" if layer3_end > 0 else "0.0%"
}
stacked_rows.append(row)
if stacked_rows:
stacked_df = pd.DataFrame(stacked_rows)
print(stacked_df.to_string(index=False))
# Save to CSV
filename = f'ms_subgroup_{metric_name}_stacked_data.csv'
stacked_df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Stacked data saved to: {filename}")
else:
print("No data available for stacked data table.")
# Table 4: Gender comparison results
print(f"\n{'=' * 80}")
print("TABLE 4: GENDER COMPARISON RESULTS")
print(f"{'=' * 80}")
if table_data['gender_comparisons']:
gender_rows = []
for lesion_type, comparisons in table_data['gender_comparisons'].items():
for metric_type, comparison_result in comparisons.items():
if isinstance(comparison_result, dict) and 'test_type' in comparison_result:
unit = 'mm²' if metric_type == 'area' else '%'
if comparison_result['parametric']:
male_stat = f"{comparison_result['group1_stats']['mean']:.3f} ± {comparison_result['group1_stats']['std']:.3f}"
female_stat = f"{comparison_result['group2_stats']['mean']:.3f} ± {comparison_result['group2_stats']['std']:.3f}"
stat_type = "Mean ± SD"
else:
male_stat = f"{comparison_result['group1_stats']['median']:.3f} [{comparison_result['group1_stats']['q25']:.3f}-{comparison_result['group1_stats']['q75']:.3f}]"
female_stat = f"{comparison_result['group2_stats']['median']:.3f} [{comparison_result['group2_stats']['q25']:.3f}-{comparison_result['group2_stats']['q75']:.3f}]"
stat_type = "Median [IQR]"
row = {
'Lesion Type': lesion_type,
'Metric': metric_type.upper(),
'Unit': unit,
'Test Used': comparison_result['test_type'],
'Statistic Type': stat_type,
'Male': male_stat,
'Female': female_stat,
'Male N': comparison_result['group1_stats']['n'],
'Female N': comparison_result['group2_stats']['n'],
'Test Statistic': f"{comparison_result['test_statistic']:.3f}",
'P-value': f"{comparison_result['p_value']}",
'Effect Size': f"{comparison_result['effect_size_type']}: {comparison_result['effect_size']:.3f}",
'Significant (α=0.05)': 'Yes' if comparison_result['p_value'] < 0.05 else 'No',
'Significant (Bonferroni α=0.0083)': 'Yes' if comparison_result[
'p_value'] < 0.0083 else 'No',
'Male Normality': 'Normal' if comparison_result['normality_assessment']['group1'][
'shapiro_normal'] else 'Non-normal',
'Female Normality': 'Normal' if comparison_result['normality_assessment']['group2'][
'shapiro_normal'] else 'Non-normal'
}
gender_rows.append(row)
if gender_rows:
gender_df = pd.DataFrame(gender_rows)
print(gender_df.to_string(index=False))
# Save to CSV
filename = 'ms_subgroup_gender_comparisons.csv'
gender_df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Gender comparison results saved to: {filename}")
else:
print("No gender comparison results available.")
else:
print("No gender comparison data available.")
# Table 5: Summary statistics for each subgroup and metric
print(f"\n{'=' * 80}")
print("TABLE 5: SUMMARY STATISTICS BY SUBGROUP AND METRIC")
print(f"{'=' * 80}")
for metric_name in ['area', 'ratio']:
unit = 'mm²' if metric_name == 'area' else '%'
metric_title = 'WMH Subtype Area' if metric_name == 'area' else 'WMH Subtype Ratio'
print(f"\n{metric_title} - Summary Across All Age Groups ({unit}):")
print("-" * 90)
# Create summary statistics table
summary_rows = []
for subgroup in table_data['metadata']['subgroups']:
for lesion_type in ['PEWMH', 'DWMH', 'JCWMH']:
if (subgroup in table_data['detailed_stats'] and
metric_name in table_data['detailed_stats'][subgroup]):
# Aggregate across all age groups for this subgroup/lesion type
all_means = []
all_medians = []
total_count = 0
for age_label in self.config.AGE_LABELS:
if (age_label in table_data['detailed_stats'][subgroup][metric_name] and
lesion_type in table_data['detailed_stats'][subgroup][metric_name][age_label]):
stats = table_data['detailed_stats'][subgroup][metric_name][age_label][lesion_type]
if stats['count'] > 0:
all_means.append(stats['mean'])
all_medians.append(stats['median'])
total_count += stats['count']
if all_means:
row = {
'Subgroup': subgroup,
'Lesion Type': lesion_type,
'Total N': total_count,
'Age Groups': len(all_means),
'Mean of Means': f"{np.mean(all_means):.2f}",
'Median of Medians': f"{np.median(all_medians):.2f}",
'Range of Means': f"{min(all_means):.2f} - {max(all_means):.2f}",
'Range of Medians': f"{min(all_medians):.2f} - {max(all_medians):.2f}"
}
summary_rows.append(row)
if summary_rows:
summary_df = pd.DataFrame(summary_rows)
print(summary_df.to_string(index=False))
# Save to CSV
filename = f'ms_subgroup_{metric_name}_summary_stats.csv'
summary_df.to_csv(os.path.join(config.OUTPUT_DIR, filename), index=False)
print(f"Summary statistics saved to: {filename}")
else:
print("No data available for summary statistics table.")
print(f"\n{'=' * 80}")
print("All MS subgroup tables have been generated and saved to CSV files!")
print(f"{'=' * 80}")
print("\nGenerated Files Summary:")
print("- ms_subgroup_area_detailed_stats.csv (Complete descriptive statistics)")
print("- ms_subgroup_ratio_detailed_stats.csv (Complete descriptive statistics)")
print("- ms_subgroup_area_plot_data.csv (Data used for visualization)")
print("- ms_subgroup_ratio_plot_data.csv (Data used for visualization)")
print("- ms_subgroup_area_stacked_data.csv (Three-layer stacked values)")
print("- ms_subgroup_ratio_stacked_data.csv (Three-layer stacked values)")
print("- ms_subgroup_gender_comparisons.csv (Statistical gender comparisons)")
print("- ms_subgroup_area_summary_stats.csv (Summary across age groups)")
print("- ms_subgroup_ratio_summary_stats.csv (Summary across age groups)")
print("- ms_subgroup_analysis_documentation.txt (Comprehensive documentation)")
print("- ms_subgroup_analysis.png (Figure)")
def _generate_ms_subgroup_documentation(self, table_data):
"""Generate comprehensive documentation explaining the MS subgroup analysis figure"""
from datetime import datetime
# Create comprehensive documentation
doc_content = f"""
MS SUBGROUP LESION ANALYSIS - COMPREHENSIVE DOCUMENTATION
=========================================================
Generated on: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}
OVERVIEW
--------
This analysis examines white matter hyperintensity (WMH) lesion subtypes specifically
in Multiple Sclerosis (MS) patients. The analysis stratifies lesions by anatomical
location and provides detailed comparisons across age groups, gender, and measurement
types (absolute area vs. normalized ratios).
FIGURE DESCRIPTION
------------------
The figure consists of a 3x2 subplot layout (16" width x 18" height):
Layout Structure:
- Row 1: All MS Patients (combined analysis)
- Row 2: Female MS Patients only
- Row 3: Male MS Patients only
- Left Column: Absolute WMH Subtype Areas (mm²)
- Right Column: WMH Subtype Ratios (%)
Subplot Details:
1. Top-Left: All MS - WMH Subtype Areas
2. Top-Right: All MS - WMH Subtype Ratios
3. Middle-Left: Female MS - WMH Subtype Areas
4. Middle-Right: Female MS - WMH Subtype Ratios
5. Bottom-Left: Male MS - WMH Subtype Areas
6. Bottom-Right: Male MS - WMH Subtype Ratios
VISUALIZATION METHOD
--------------------
Chart Type: Three-Layer Stacked Area Plot
- Each subplot uses three-layer stacked area charts showing lesion subtypes across age groups
- PEWMH layer (bottom): Fills from 0 to PEWMH mean value
- DWMH layer (middle): Fills from PEWMH to PEWMH + DWMH mean
- JCWMH layer (top): Fills from PEWMH + DWMH to total mean
- This visualization shows both individual subtype contributions and total lesion burden
Color Scheme:
- PEWMH: {table_data['metadata']['colors'].get('pewmh', 'Color defined in config')} (alpha=0.8)
- DWMH: {table_data['metadata']['colors'].get('dwmh', 'Color defined in config')} (alpha=0.8)
- JCWMH: {table_data['metadata']['colors'].get('jcwmh', 'Color defined in config')} (alpha=0.8)
AGE STRATIFICATION
------------------
Age Groups: {', '.join(table_data['metadata']['age_labels'])}
Age Bins: {table_data['metadata']['age_bins']}
Age Centers (for plotting): {[f'{center:.1f}' for center in table_data['metadata']['age_centers']]}
The analysis stratifies data across these age groups to examine age-related changes
in lesion subtype distribution for MS patients.
LESION SUBTYPES ANALYZED
------------------------
1. PEWMH (Periventricular White Matter Hyperintensities):
- Location: Adjacent to the ventricular system
- Clinical significance: Often associated with MS pathology and severity
- Area Column: {table_data['metadata']['metrics']['area']['columns']['pewmh']}
- Ratio Column: peri_wmh_ratio
2. DWMH (Deep White Matter Hyperintensities):
- Location: Near but not directly adjacent to ventricles
- Clinical significance: May represent different pathophysiological processes
- Area Column: {table_data['metadata']['metrics']['area']['columns']['dwmh']}
- Ratio Column: deep_wmh_ratio
3. JCWMH (Juxtacortical White Matter Hyperintensities):
- Location: At the interface between white and gray matter
- Clinical significance: Associated with cortical involvement in MS
- Area Column: {table_data['metadata']['metrics']['area']['columns']['jcwmh']}
- Ratio Column: juxta_wmh_ratio
CLINICAL CONTEXT
----------------
MS Lesion Distribution Patterns:
- MS lesions preferentially affect certain brain regions
- Periventricular regions are classically involved in MS
- Juxtacortical lesions may indicate disease progression
- Age-related changes may reflect disease evolution or natural aging
Expected Clinical Patterns:
- PEWMH typically most prominent in MS patients
- Age-related increase in all lesion subtypes
- Gender differences may reflect hormonal or genetic factors
- Individual variation in lesion distribution patterns
MEASUREMENT TYPES
-----------------
1. Absolute Area (mm²):
- Direct measurement of lesion area
- Units: Square millimeters (mm²)
- Clinical significance: Reflects total lesion load per subtype
2. Normalized Ratio (%):
- Lesion area relative to total brain area/volume
- Units: Percentage (%)
- Clinical significance: Controls for individual brain size differences
STATISTICAL APPROACH
--------------------
For each combination of:
- Subgroup (All MS, Female MS, Male MS)
- Age group ({len(table_data['metadata']['age_labels'])} categories)
- Lesion subtype (PEWMH, DWMH, JCWMH)
- Metric (Area vs Ratio)
The following statistics are calculated:
- Sample size (N)
- Mean ± Standard Deviation (for visualization)
- Median and Interquartile Range [Q25-Q75] (for robust statistics)
- Minimum and Maximum values
- 25th and 75th percentiles
Non-parametric Statistics:
- Mann-Whitney U tests for gender comparisons within each lesion subtype
- Median and IQR reported for robustness to outliers
- Appropriate for skewed lesion distribution data
INTERPRETATION GUIDELINES
-------------------------
Three-Layer Stacked Plot Interpretation:
- Bottom layer height = PEWMH mean contribution
- Middle layer height = DWMH mean contribution
- Top layer height = JCWMH mean contribution
- Total stack height = Combined lesion burden across all subtypes
- Layer thickness indicates relative contribution of each subtype
Clinical Pattern Recognition:
- Dominant lesion subtype can be identified by layer thickness
- Age-related changes visible as slope steepness
- Gender differences apparent by comparing male vs female rows
- Subtype-specific patterns may indicate different pathological processes
Expected Subtype Hierarchy:
- PEWMH often dominant in MS (thickest layer)
- DWMH and JCWMH may show age-dependent changes
- Individual variation in subtype distribution patterns
GENDER STRATIFICATION ANALYSIS
-------------------------------
The analysis includes separate visualizations for:
1. All MS patients (combined analysis)
2. Female MS patients only
3. Male MS patients only
Gender Comparison Features:
- Direct visual comparison between male and female patterns
- Statistical testing for gender differences in each lesion subtype
- Separate analysis for both area and ratio measurements
- Age-stratified patterns within each gender
DATA QUALITY CONSIDERATIONS
----------------------------
- Zero values indicate no subjects in that age/subtype combination
- Small sample sizes in gender-stratified analyses may reduce statistical power
- Lesion subtype classification depends on anatomical definition accuracy
- Manual segmentation variability may affect subtype boundaries
- Automated methods may have subtype-specific detection biases
STATISTICAL TESTING METHODOLOGY
--------------------------------
Gender Comparisons:
- Mann-Whitney U test for each lesion subtype
- Separate tests for area and ratio measurements
- Non-parametric approach suitable for skewed lesion data
- Two-sided alternative hypothesis
Multiple Testing Considerations:
- Multiple comparisons performed across lesion subtypes
- Consider Bonferroni correction: α = 0.05/6 = 0.0083 for significance
- False Discovery Rate (FDR) correction may be more appropriate
OUTPUT FILES GENERATED
-----------------------
1. Figure: ms_subgroup_analysis.png
- 3x2 subplot layout with three-layer stacked area plots
- High resolution (DPI: {getattr(self.config, 'DPI', 300)})
- White background for publication quality
2. Detailed Statistics Tables (CSV):
- ms_subgroup_area_detailed_stats.csv: Complete descriptive statistics for area measurements
- ms_subgroup_ratio_detailed_stats.csv: Complete descriptive statistics for ratio measurements
3. Plot Data Tables (CSV):
- ms_subgroup_area_plot_data.csv: Mean values and counts used for area visualization
- ms_subgroup_ratio_plot_data.csv: Mean values and counts used for ratio visualization
4. Stacked Area Values (CSV):
- ms_subgroup_area_stacked_data.csv: Layer boundaries and contributions for area plots
- ms_subgroup_ratio_stacked_data.csv: Layer boundaries and contributions for ratio plots
5. Gender Comparison Tables (CSV):
- ms_subgroup_gender_comparisons.csv: Statistical test results comparing males vs females
6. Summary Statistics Tables (CSV):
- ms_subgroup_area_summary_stats.csv: Aggregated statistics across age groups for areas
- ms_subgroup_ratio_summary_stats.csv: Aggregated statistics across age groups for ratios
7. This Documentation:
- ms_subgroup_analysis_documentation.txt: Complete explanation of analysis and interpretation
TECHNICAL SPECIFICATIONS
-------------------------
Figure Specifications:
- Size: 16" x 18" (width x height) - taller for 3-row layout
- DPI: {getattr(self.config, 'DPI', 300)}
- Background: White
- Font sizes: Title=14pt (bold), Axis labels=12pt
- Grid: Enabled with 30% transparency
- Legend: Three-layer legend for each subplot
Data Processing:
- MS patients only (HC excluded from this analysis)
- Missing data handled by excluding from calculations (dropna)
- Zero values used when no subjects available in category
- Robust statistics (median/IQR) preferred for group summaries
Plotting Library: matplotlib
Statistical Library: scipy.stats (Mann-Whitney U tests)
Data Processing: pandas, numpy
TABLE DESCRIPTIONS
------------------
Table 1 - Detailed Statistics:
Contains complete descriptive statistics (N, mean, SD, median, Q25, Q75, min, max)
for each combination of subgroup, age group, and lesion subtype.
Table 2 - Plot Data:
Contains the exact mean values and sample counts used to generate the stacked area plots,
organized by subgroup and age group.
Table 3 - Stacked Area Values:
Shows the layer boundaries and individual contributions for the three-layer stacked plots,
including percentage contributions of each lesion subtype.
Table 4 - Gender Comparisons:
Statistical test results (Mann-Whitney U) comparing male vs female patients for each
lesion subtype, with both uncorrected and Bonferroni-corrected significance levels.
Table 5 - Summary Statistics:
Aggregated statistics across all age groups for each subgroup and lesion subtype,
showing overall patterns and variability.
LIMITATIONS AND CONSIDERATIONS
------------------------------
1. Sample Size Limitations:
- Gender-stratified analyses have reduced sample sizes
- Some age groups may have insufficient subjects for reliable estimates
- Power analysis recommended for gender comparisons
2. Lesion Subtype Definition:
- Anatomical boundaries between subtypes may be arbitrary
- Different segmentation protocols may yield different results
- Spatial resolution limits affecting small lesion detection
3. Multiple Comparisons:
- Six statistical tests performed (3 subtypes × 2 metrics)
- Risk of Type I error inflation
- Consider correction for multiple testing
4. Age Group Effects:
- Discretized age groups may mask continuous relationships
- Unequal age distributions between genders possible
- Cross-sectional design limits inferences about progression
5. MS Disease Heterogeneity:
- MS subtypes (relapsing-remitting, progressive) not considered
- Disease duration effects not analyzed
- Treatment effects not controlled
RECOMMENDED FOLLOW-UP ANALYSES
------------------------------
1. Disease Subtype Stratification:
- Separate analysis for RRMS, SPMS, PPMS if sample size permits
- Include disease duration as covariate
2. Advanced Statistical Modeling:
- Multivariate analysis of lesion subtype interdependencies
- Machine learning approaches for subtype pattern classification
- Longitudinal analysis if follow-up data available
3. Clinical Correlation Studies:
- Correlation with disability scores (EDSS, MSFC)
- Cognitive function associations
- Treatment response predictions
4. Spatial Analysis:
- Lesion location heat maps
- Connectivity-based lesion impact analysis
- Atlas-based regional quantification
5. Comparative Studies:
- Comparison with other neurological conditions
- Validation in independent MS cohorts
- Cross-scanner reproducibility studies
QUALITY CONTROL RECOMMENDATIONS
--------------------------------
1. Segmentation Validation:
- Inter-rater reliability assessment for lesion subtype classification
- Comparison of automated vs manual segmentation methods
- Test-retest reliability studies
2. Clinical Validation:
- Correlation with established MS biomarkers
- Agreement with radiological assessment
- Validation against histopathological data if available
3. Statistical Validation:
- Power analysis for gender comparisons
- Bootstrap confidence intervals for robust statistics
- Cross-validation of predictive models
CONTACT AND METHODOLOGY
-----------------------
This analysis was generated using an automated pipeline for MS lesion subtype assessment.
For questions about clinical interpretation, statistical methods, or lesion classification
protocols, refer to the original research protocol and neuroimaging analysis guidelines.
Analysis Pipeline Version: [Version info if available]
Last Updated: {datetime.now().strftime('%Y-%m-%d')}
REFERENCES AND FURTHER READING
-------------------------------
1. Filippi et al. (2019). Assessment of lesions on magnetic resonance imaging
in multiple sclerosis: practical guidelines. Brain.
2. Geurts et al. (2012). Cortical lesions in multiple sclerosis: combined
postmortem MR imaging and histopathology. AJNR Am J Neuroradiol.
3. Brownell & Hughes (1962). The distribution of plaques in the cerebrum in
multiple sclerosis. Journal of Neurology, Neurosurgery & Psychiatry.
4. Barkhof & Filippi (2009). MRI in multiple sclerosis. Journal of Magnetic
Resonance Imaging.
5. Thompson et al. (2018). Diagnosis of multiple sclerosis: 2017 revisions of
the McDonald criteria. Lancet Neurology.
END OF DOCUMENTATION
====================
"""
# Save documentation to file
doc_filename = os.path.join(config.OUTPUT_DIR, 'ms_subgroup_analysis_documentation.txt')
with open(doc_filename, 'w', encoding='utf-8') as f:
f.write(doc_content)
print(f"\n{'=' * 80}")
print("COMPREHENSIVE MS SUBGROUP DOCUMENTATION GENERATED")
print(f"{'=' * 80}")
print(f"Documentation saved to: ms_subgroup_analysis_documentation.txt")
print(f"File contains detailed explanation of MS lesion subtype analysis and clinical interpretation.")
def correlation_analysis(self):
"""

p4_excel_analysis_developed.py at commit 3e9edd4, under MIT · at the source

Overview

Authors: Mahdi Bashiri Bawil¹, Mousa Shamsi¹, Abolhassan Shakeri Bavil²

ORCID iDs: Mahdi Bashiri Bawil, Mousa Shamsi, Abolhassan Shakeri Bavil

Institutions: Sahand University of Technology (Iran); Tabriz University of Medical Sciences (Iran)

Journal: BMC medical imaging, volume 26, issue 1, article 387

Dates: received 7 December 2025; accepted 27 May 2026; published online 1 June 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1186/s12880-026-02481-2 · PMID 42226149 · PMCID PMC13449394 · OpenAlex W7163025612

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: structural MRI / diffusion (modality), human (organism), multiple sclerosis (population), clinical / translational (subfield)

Methods: Connectivity, Statistics, Preprocessing, Machine learning

Keywords: Multiple sclerosis (MS), Deep learning, Neuroimaging, MRI, Population-specific biomarkers, FLAIR segmentation

MeSH: Brain*, Deep Learning*, Multiple Sclerosis*, Neuroimaging*, Adult, Female, Humans, Iran, Magnetic Resonance Imaging, Male, Middle Aged, Middle Eastern People, Retrospective Studies, Young Adult (* major topic)

Topic: Multiple Sclerosis Research Studies (Pathology and Forensic Medicine, Medicine), according to OpenAlex

Citations: not cited yet (Europe PMC); 59 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 23 matches between paragraphs and lines of code.

Mahdi-Bashiri/MS-DeepBrain-Study

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 3e9edd415b0beccf7ad19f036fc3dbfdcf1d8076, 22 April 2026

Languages: Python (29)

Size: 455 files, 29 scripts

Software Heritage: not archived

Found in: “Data availability”

Holds: README, license file, environment (requirements.txt)

Not found: CITATION.cff, tests, continuous integration, documentation

Tools: NumPy (20 files), TensorFlow (15 files), Matplotlib (10 files), Keras (9 files), scikit-image (9 files), NiBabel (8 files), OpenCV (8 files), scikit-learn (8 files), SciPy (6 files), pandas (5 files), seaborn (4 files), FSL (2 files), imageio (1 file), Pillow (1 file), Plotly (1 file), pydicom (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

32 files

Phase1_data_preprocessing/p4_preprocess.py, Python, 800 lines
Phase2_data_preparation_for_model_training/generating_3L_masks.py, Python, 351 lines
Phase3_model_training_and_inferencing_and_evaluation/for_GM/model_training_scripts/p1_compute_class_weights.py, Python, 336 lines
Phase3_model_training_and_inferencing_and_evaluation/for_GM/model_training_scripts/p1_data_loader.py, Python, 847 lines, 1 match
Phase3_model_training_and_inferencing_and_evaluation/for_GM/model_training_scripts/p1_pix2pix_var5.py, Python, 1,313 lines, 1 match
Phase3_model_training_and_inferencing_and_evaluation/for_GM/model_training_scripts/p1_predict_new_data_gm.py, Python, 477 lines
Phase3_model_training_and_inferencing_and_evaluation/for_GM/model_training_scripts/unet_model.py, Python, 87 lines
Phase3_model_training_and_inferencing_and_evaluation/for_GM/model_training_scripts/utility_functions.py, Python, 97 lines
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/data_splits/for_assignment.py, Python, 234 lines
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/attn_unet_model.py, Python, 85 lines
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/base_runner_all.py, Python, 23 lines
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/dlv3_unet_model.py, Python, 198 lines, 1 match
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/dlv3_unet_model_GN.py, Python, 247 lines, 2 matches
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_compute_class_weights.py, Python, 353 lines
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_data_loader.py, Python, 912 lines, 1 match
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_error_analysis.py, Python, 1,033 lines
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_folds_results_aggregator.py, Python, 611 lines
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_inference.py, Python, 1,146 lines, 2 matches
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_run_experiments_all.py, Python, 576 lines
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_unet_viz.py, Python, 640 lines
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/p4_variant_all_net.py, Python, 1,051 lines, 2 matches
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/trans_unet_model.py, Python, 125 lines, 2 matches
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/unet_model.py, Python, 87 lines
Phase3_model_training_and_inferencing_and_evaluation/for_WMH_Vent/model_training_scripts/utility_functions.py, Python, 96 lines
Phase4_data_processing/p4_brain_mri_extractor_excel.py, Python, 349 lines, 1 match
Phase4_data_processing/p4_core_processing.py, Python, 1,069 lines
Phase4_data_processing/p4_predict_new_data.py, Python, 516 lines, 1 match
Phase4_data_processing/p4_preprocess.py, Python, 800 lines
Phase5_statistical_analysis/p4_excel_analysis_developed.py, Python, 4,176 lines, 9 matches
LICENSE, License, 21 lines
LICENSE.txt, License, 21 lines
README.md, Text, 559 lines

At the source: github.com/Mahdi-Bashiri/MS-DeepBrain-Study

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
29 scripts, each with its path and the digest of its content;
23 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: Mahdi-Bashiri/MS-DeepBrain-Study

Read it in the paper: doi.org/10.1186/s12880-026-02481-2.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 3 authors, 6 keywords, 14 MeSH terms, 58 references.

Cite

This paper

Bawil, M. B., Shamsi, M., & Bavil, A. S. (2026). Deep learning-based neuroanatomical profiling reveals population-specific brain changes in multiple sclerosis: a large-scale Middle Eastern study. BMC medical imaging, 26(1), 387. https://doi.org/10.1186/s12880-026-02481-2

BibTeX

@article{bawil2026deep, author = {Bawil, Mahdi Bashiri and Shamsi, Mousa and Bavil, Abolhassan Shakeri}, title = {{Deep learning-based neuroanatomical profiling reveals population-specific brain changes in multiple sclerosis: a large-scale Middle Eastern study}}, journal = {BMC medical imaging}, year = {2026}, month = jun, volume = {26}, number = {1}, pages = {387}, publisher = {BMC}, issn = {1471-2342}, doi = {10.1186/s12880-026-02481-2}, url = {https://doi.org/10.1186/s12880-026-02481-2}, pmid = {42226149}, pmcid = {PMC13449394} }

RIS

TY - JOUR AU - Bawil, Mahdi Bashiri AU - Shamsi, Mousa AU - Bavil, Abolhassan Shakeri TI - Deep learning-based neuroanatomical profiling reveals population-specific brain changes in multiple sclerosis: a large-scale Middle Eastern study T2 - BMC medical imaging J2 - BMC Med Imaging PY - 2026 DA - 2026/06/01 VL - 26 IS - 1 SP - 387 SN - 1471-2342 PB - BMC DO - 10.1186/s12880-026-02481-2 UR - https://doi.org/10.1186/s12880-026-02481-2 LA - en ER -

CSL-JSON

{ "id": "10.1186/s12880-026-02481-2", "type": "article-journal", "title": "Deep learning-based neuroanatomical profiling reveals population-specific brain changes in multiple sclerosis: a large-scale Middle Eastern study", "container-title": "BMC medical imaging", "author": [ { "family": "Bawil", "given": "Mahdi Bashiri" }, { "family": "Shamsi", "given": "Mousa" }, { "family": "Bavil", "given": "Abolhassan Shakeri" } ], "container-title-short": "BMC Med Imaging", "volume": "26", "issue": "1", "page": "387", "DOI": "10.1186/s12880-026-02481-2", "PMID": "42226149", "PMCID": "PMC13449394", "ISSN": "1471-2342", "publisher": "BMC", "URL": "https://doi.org/10.1186/s12880-026-02481-2", "language": "en", "issued": { "date-parts": [ [ 2026, 6, 1 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41597-026-07184-5 [code]: A Multiple Sclerosis MRI Dataset with Tri-Mask Annotations for Lesion Segmentation.
Journal: Scientific data
In common: pydicom, Keras, TensorFlow, 9 other tools, multiple sclerosis, structural MRI / diffusion, 6 references, 3 authors
[2] doi:10.1186/s12938-026-01555-0 [code]: Incorporating normal periventricular changes for enhanced pathological white matter hyperintensity segmentation: on multiclass deep learning approaches.
Journal: Biomedical engineering online
In common: Keras, TensorFlow, OpenCV, 7 other tools, structural MRI / diffusion, 5 references, 3 authors
[3] doi:10.1002/epi.70296 [code]: Fully automated three-dimensional deep learning-based magnetic resonance imaging segmentation of brain cavities in epilepsy surgery.
Journal: Epilepsia
In common: imageio, Keras, TensorFlow, 10 other tools, structural MRI / diffusion, clinical / translational, 1 reference
[4] doi:10.3389/fnins.2026.1870124 [code]: An end-to-end pipeline for automated fetal brain segmentation and biometry from 3D SSFP MRI.
Journal: Frontiers in neuroscience
In common: imageio, FSL, OpenCV, 8 other tools, structural MRI / diffusion, 4 references
[5] doi:10.3389/frai.2026.1771088 [code]: Few-shot deployment of pretrained MRI transformers in brain imaging tasks.
Journal: Frontiers in artificial intelligence
In common: pydicom, imageio, OpenCV, 9 other tools, structural MRI / diffusion, 1 reference
[6] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: pydicom, Keras, TensorFlow, 10 other tools
[7] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: imageio, Keras, TensorFlow, 10 other tools
[8] doi:10.1371/journal.pdig.0001316 [code]: Reliability of a convolutional neural network in segmenting multiple sclerosis lesions from MRI: Impact of data augmentation, image modality and tolerance with U-Net architecture.
Journal: PLOS digital health
In common: scikit-image, NiBabel, scikit-learn, 3 other tools, multiple sclerosis, structural MRI / diffusion, 6 references
[9] doi:10.1080/07853890.2026.2685416 [code]: Pulmonary and cerebral damage in COVID-19 survivors: is there any association?
Journal: Annals of medicine
In common: pydicom, imageio, Keras, 8 other tools, structural MRI / diffusion, clinical / translational
[10] doi:10.1371/journal.pcbi.1014571 [code]: SynAPSeg: A novel dataset and image analysis framework for deep learning-based synapse detection and quantification.
Journal: PLoS computational biology
In common: imageio, Keras, TensorFlow, 9 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.