Neuroanatomy reflects individual variability in impulsivity in youth.

Code ↔ Paper

2 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 2 matches

[1] § Methods › Predictive modeling ↔ nbdc_abcd_bpm_smri_impulsivity.ipynb, lines 611 to 760 · score 0.59 · evaluated model, permuting, shuffled, variable, optimized, split · Report a mismatch
[2] § Methods › Predictive modeling ↔ nbdc_abcd_bpm_smri_impulsivity.ipynb, lines 443 to 489 · score 0.53 · ridge regression models, females, thickness, volume, brain · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: 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

Jupyter notebook · 2,873 lines · 116 KB · no license · 2 matches

# %% [markdown]
# ## Import Libraries and Packages
# %%
#general imports
import os
import sys
import warnings
import numpy as np
import pandas as pd
import pickle
#suppress warnings
warnings.filterwarnings('ignore')
#sklearn imports
from sklearn.metrics import explained_variance_score
from sklearn.metrics.pairwise import cosine_similarity
from sklearn.linear_model import Ridge
from sklearn.model_selection import (train_test_split, GridSearchCV, cross_val_score,
GroupKFold, GroupShuffleSplit, StratifiedKFold)
#tqdm import
from tqdm import tqdm
#stats import
from statsmodels.stats.multitest import fdrcorrection as fdr
#plotting import
import matplotlib.pyplot as plt
import seaborn as sns
# %% [markdown]
# ## Load in Data
# %%
#set directories
base_dir = '/Users/elvishadhamala/OneDrive_Northwell/NBDC_PHI/'
results_dir = '/Users/elvishadhamala/Documents/BPM/nbdc_analyses/impulsivity/unrelatedparticipants/'
#define data files to load
data_files = {
'gen': 'abcd_6.0/core/g.csv',
'qc': 'abcd_6.0/core/mri/mr_qc.csv',
'adm': 'abcd_6.0/core/mri/mr_adm.csv',
'smri': 'abcd_6.0/core/mri/mr_smri.csv',
'mh': 'abcd_6.0/core/mh.csv'
}
#load all data files
print("Loading ABCD 6.0 data files...")
data = {}
for name, filepath in data_files.items():
data[name] = pd.read_csv(os.path.join(base_dir, filepath), header=0)
print(f"{name}: {data[name].shape}")
print("Data loading complete!")
#extract to individual variables (if you prefer)
gen, qc, adm, smri, mh = data['gen'], data['qc'], data['adm'], data['smri'], data['mh']
# %% [markdown]
# ## Clean Up Data
# %%
#function to apply standard QC criteria for sMRI
def smri_qc_filter(qc_df):
"""Apply standard sMRI QC criteria"""
mask = (qc_df.mr_y_qc__incl__smri__t1_indicator == 1) & (qc_df.mr_y_qc__clfind_score < 3)
print(f"sMRI QC: {len(qc_df)} → {mask.sum()} subjects/timepoints")
return qc_df[mask]
#perform sMRI QC
qc_pass = smri_qc_filter(qc)
#get QC-passed participant/session combinations
qc_lookup = qc_pass[['participant_id', 'session_id']]
#filter all other dataframes based on imaging QC
gen_qc = gen.merge(qc_lookup, on=['participant_id', 'session_id'])
mh_qc = mh.merge(qc_lookup, on=['participant_id', 'session_id'])
smri_qc = smri.merge(qc_lookup, on=['participant_id', 'session_id'])
adm_qc = adm.merge(qc_lookup, on=['participant_id', 'session_id'])
print(f"General: {len(gen)} → {len(gen_qc)}")
print(f"Mental Health: {len(mh)} → {len(mh_qc)}")
print(f"Imaging: {len(smri)} → {len(smri_qc)}")
print(f"MR Admin: {len(adm)} → {len(adm_qc)}")
#validation check
def validate_qc(qc_df, other_df, df_name):
qc_set = set(zip(qc_df.participant_id, qc_df.session_id))
data_set = set(zip(other_df.participant_id, other_df.session_id))
invalid = data_set - qc_set
print(f"✓ {df_name}: All data passed QC" if len(invalid) == 0 else f"✗ {df_name}: {len(invalid)} entries didn't pass QC!")
validate_qc(qc_pass, gen_qc, "General")
validate_qc(qc_pass, mh_qc, "Mental Health")
validate_qc(qc_pass, smri_qc, "Imaging")
validate_qc(qc_pass, adm_qc, "MR Admin")
# %%
#find relevant mental health data
#cbcl_cols = [col for col in mh_qc.columns if col.startswith('mh_p_cbcl') and col.endswith('_sum')]
#pps_cols = ['mh_y_pps__severity_score', 'mh_y_pps_count']
#gbi_cols = ['mh_p_gbi_sum']
#bpm_cols = [col for col in mh_qc.columns if col.startswith('mh_t_bpm') and col.endswith('_sum')]
mh_cols = ['mh_y_bisbas__bis_sum',
'mh_y_bisbas__bas__rr_sum', 'mh_y_bisbas__bas__dr_sum', 'mh_y_bisbas__bas__fs_sum',
'mh_y_upps__nurg_sum', 'mh_y_upps__purg_sum', 'mh_y_upps__pers_sum', 'mh_y_upps__plan_sum', 'mh_y_upps__sens_sum']
mh_subset = mh_qc[['participant_id', 'session_id'] + mh_cols]
#keep only participants with complete MH data
mh_complete = mh_subset.dropna()
print(f"MH variables: {len(mh_cols)}")
print(f"Complete cases: {len(mh_qc)} → {len(mh_complete)}")
# %%
#get the participant/session combinations with complete MH data
mh_lookup = mh_complete[['participant_id', 'session_id']]
#filter all other dataframes to only include those with complete MH data
gen_mh = gen_qc.merge(mh_lookup, on=['participant_id', 'session_id'])
smri_mh = smri_qc.merge(mh_lookup, on=['participant_id', 'session_id'])
adm_mh = adm_qc.merge(mh_lookup, on=['participant_id', 'session_id'])
print("After filtering for complete MH data:")
print(f"General: {len(gen_qc)} → {len(gen_mh)}")
print(f"Imaging: {len(smri_qc)} → {len(smri_mh)}")
print(f"MR Admin: {len(adm_qc)} → {len(adm_mh)}")
print(f"MH: {len(mh_complete)} (complete cases)")
# ==========================================
# ENSURE COMPLETE IMAGING DATA
# ==========================================
#get participant/session combinations that have BOTH complete MH AND complete imaging
smri_mh = smri_mh.dropna()
complete_imaging_lookup = smri_mh[['participant_id', 'session_id']]
#filter all dataframes to only include those with BOTH complete MH and imaging
gen_complete = gen_mh.merge(complete_imaging_lookup, on=['participant_id', 'session_id'])
smri_complete = smri_mh.merge(complete_imaging_lookup, on=['participant_id', 'session_id'])
adm_complete = adm_mh.merge(complete_imaging_lookup, on=['participant_id', 'session_id'])
mh_complete = mh_complete.merge(complete_imaging_lookup, on=['participant_id', 'session_id'])
print("\nAfter filtering for complete MH AND imaging data:")
print(f"General: {len(gen_mh)} → {len(gen_complete)}")
print(f"Imaging: {len(smri_mh)} → {len(smri_complete)}")
print(f"MR Admin: {len(adm_mh)} → {len(adm_complete)}")
print(f"MH: {len(mh_complete)} → {len(mh_complete)}")
# %%
# ==========================================
# FAMILY FILTERING - ONE PARTICIPANT PER FAMILY
# ==========================================
if(results_dir=='/Users/elvishadhamala/Documents/BPM/nbdc_analyses/impulsivity/unrelatedparticipants/'):
print("\n" + "="*60)
print("FILTERING TO ONE PARTICIPANT PER FAMILY")
print("="*60)
def select_one_participant_per_family(data_dict, family_col='ab_g_stc__design_id__fam',
participant_col='participant_id', session_col='session_id',
random_seed=42):
"""
Select one participant per family, keeping all their timepoints
"""
print("Selecting one participant per family...")
# Get participant-family mapping from general data
participant_info = data_dict['gen_complete'][[participant_col, family_col, session_col]].copy()
participant_info = participant_info.dropna(subset=[family_col])
print(f"Total participant-session records: {len(participant_info)}")
print(f"Unique participants: {participant_info[participant_col].nunique()}")
print(f"Unique families: {participant_info[family_col].nunique()}")
# Get unique participant-family pairs
participant_family_map = participant_info[[participant_col, family_col]].drop_duplicates()
print(f"Participant-family combinations: {len(participant_family_map)}")
# Count participants per family
family_counts = participant_family_map[family_col].value_counts()
print(f"Family size distribution:")
for size, count in family_counts.value_counts().sort_index().items():
print(f" {size} participant(s): {count} families")
# For each family, randomly select one participant
np.random.seed(random_seed)
selected_participants = []
for family_id in participant_family_map[family_col].unique():
family_participants = participant_family_map[
participant_family_map[family_col] == family_id
][participant_col].tolist()
# Randomly select one participant from this family
selected_participant = np.random.choice(family_participants)
selected_participants.append(selected_participant)
print(f"Selected {len(selected_participants)} participants (one per family)")
# Create lookup for selected participants (all their sessions)
selected_participant_sessions = []
for participant_id in selected_participants:
participant_sessions = participant_info[
participant_info[participant_col] == participant_id
][[participant_col, session_col]]
selected_participant_sessions.append(participant_sessions)
if selected_participant_sessions:
selected_lookup_df = pd.concat(selected_participant_sessions, ignore_index=True)
else:
selected_lookup_df = pd.DataFrame(columns=[participant_col, session_col])
print(f"Total participant-session combinations: {len(selected_lookup_df)}")
# Filter all dataframes
filtered_data = {}
for name, df in data_dict.items():
if participant_col in df.columns and session_col in df.columns:
filtered_df = df.merge(selected_lookup_df, on=[participant_col, session_col])
filtered_data[name] = filtered_df
print(f"Filtered {name}: {len(df)} → {len(filtered_df)} records")
else:
filtered_data[name] = df
# Create summary DataFrame
selected_summary = pd.DataFrame({
participant_col: selected_participants
})
# Add family info
selected_summary = selected_summary.merge(
participant_family_map, on=participant_col
)
print(f"\nFamily filtering summary:")
print(f" Original participants: {participant_info[participant_col].nunique()}")
print(f" Original families: {participant_info[family_col].nunique()}")
print(f" Selected participants: {len(selected_participants)}")
print(f" Selected families: {selected_summary[family_col].nunique()}")
print(f" Records per dataset: {len(selected_lookup_df)}")
# Verify one participant per family
family_check = selected_summary[family_col].value_counts()
max_participants_per_family = family_check.max()
if max_participants_per_family > 1:
print(f" ⚠️ WARNING: Some families have more than one participant (max: {max_participants_per_family})")
else:
print(f" ✅ Confirmed: Each family has exactly one participant")
return filtered_data, selected_summary
# Apply family filtering
print("Applying family filtering to complete datasets...")
# Create dictionary of complete datasets
complete_data = {
'gen_complete': gen_complete,
'smri_complete': smri_complete,
'mh_complete': mh_complete,
'adm_complete': adm_complete
}
# Apply family filtering
filtered_family_data, selected_families = select_one_participant_per_family(
complete_data,
family_col='ab_g_stc__design_id__fam',
participant_col='participant_id',
session_col='session_id',
random_seed=42 # Set seed for reproducibility
)
# Update your dataframes
gen_complete = filtered_family_data['gen_complete']
smri_complete = filtered_family_data['smri_complete']
mh_complete = filtered_family_data['mh_complete']
adm_complete = filtered_family_data['adm_complete']
print(f"\nFinal sample sizes after family filtering:")
print(f"General: {len(gen_complete)}")
print(f"Imaging: {len(smri_complete)}")
print(f"MH: {len(mh_complete)}")
print(f"Admin: {len(adm_complete)}")
# Show sample sizes by timepoint
print(f"\nSample sizes by timepoint after family filtering:")
timepoints = sorted(gen_complete['session_id'].unique())
for tp in timepoints:
n_participants = len(gen_complete[gen_complete['session_id'] == tp])
n_families = gen_complete[gen_complete['session_id'] == tp]['ab_g_stc__design_id__fam'].nunique()
n_males = np.sum(gen_complete[gen_complete['session_id'] == tp]['ab_g_stc__cohort_sex'] == 1)
n_females = np.sum(gen_complete[gen_complete['session_id'] == tp]['ab_g_stc__cohort_sex'] == 2)
print(f" {tp}: {n_participants} participants from {n_families} families,{n_males} males, {n_females} females")
# Final verification
print(f"\nFinal verification:")
total_participants = gen_complete['participant_id'].nunique()
total_families = gen_complete['ab_g_stc__design_id__fam'].nunique()
print(f" Unique participants: {total_participants}")
print(f" Unique families: {total_families}")
if total_participants == total_families:
print(f" ✅ SUCCESS: Participants = Families (one participant per family)")
else:
print(f" ⚠️ WARNING: Participants ≠ Families")
# Show family information
print(f"\nFamily selection details:")
print(f" Random seed used: 42")
print(f" Selection method: Random choice per family")
print(f" All timepoints preserved for selected participants")
print("✅ Family filtering complete!")
else:
print("Family filtering skipped.")
# %%
#extract structural columns, excluding the summary totals
thickness_cols = [col for col in smri_complete.columns
if col.startswith('mr_y_smri__thk__dsk')
and col.endswith('mean')
and col != 'mr_y_smri__thk__dsk__lh_mean'
and col != 'mr_y_smri__thk__dsk__rh_mean'
and col != 'mr_y_smri__thk__dsk_mean']
area_cols = [col for col in smri_complete.columns
if col.startswith('mr_y_smri__area__dsk')
and col.endswith('sum')
and col != 'mr_y_smri__area__dsk__lh_sum'
and col != 'mr_y_smri__area__dsk__rh_sum'
and col != 'mr_y_smri__area__dsk_sum']
volume_cols = [col for col in smri_complete.columns
if col.startswith('mr_y_smri__vol__dsk')
and col.endswith('sum')
and col != 'mr_y_smri__vol__dsk__lh_sum'
and col != 'mr_y_smri__vol__dsk__rh_sum'
and col != 'mr_y_smri__vol__dsk_sum']
subc = ['_ab', '_ag', '_cbc', '_cd', '_hc', '_pl', '_pt', '_th', '_vdc', '_bs']
subcvolume_cols = [col for col in smri_complete.columns
if col.startswith('mr_y_smri__vol__aseg')
and col.endswith('sum')
and any(phrase in col for phrase in subc)]
icv_col = ['mr_y_smri__vol__aseg__icv_sum']
# all structural measures (individual regions only)
all_structural_cols = thickness_cols + area_cols + volume_cols + subcvolume_cols + icv_col
smri_all = smri_complete[['participant_id', 'session_id'] + all_structural_cols].copy()
smri_all_complete = smri_all.dropna()
print(f"All structural measures: {len(smri_all)} → {len(smri_all_complete)} complete cases")
print(f"Total measures: {len(all_structural_cols)} ({len(thickness_cols)} thickness + {len(area_cols)} area + {len(volume_cols)} volume + {len(subcvolume_cols)} subcortical volume + 1 ICV)")
# %%
#make a copy of the sMRI data to work with
smri_normalized = smri_all_complete.copy()
#get the ICV values for division
icv_values = smri_normalized['mr_y_smri__vol__aseg__icv_sum']
#apply ICV proportion correction to area columns
for col in area_cols:
smri_normalized[col] = smri_normalized[col] / icv_values
#apply ICV proportion correction to volume columns
for col in volume_cols:
smri_normalized[col] = smri_normalized[col] / icv_values
#apply ICV proportion correction to subcortical volume columns
for col in subcvolume_cols:
smri_normalized[col] = smri_normalized[col] / icv_values
#thickness columns remain unchanged
print(f"Normalized {len(area_cols)} area columns and {len(volume_cols)} volume columns and {len(subcvolume_cols)} subcortical volume columns by ICV")
print(f"Thickness columns ({len(thickness_cols)}) left unchanged")
print(f"Final dataset shape: {smri_normalized.shape}")
#check the normalization worked for area
print(f"\nExample area column before/after normalization:")
print(f"Original range: {smri_all_complete[area_cols[0]].min():.2f} - {smri_all_complete[area_cols[0]].max():.2f}")
print(f"Normalized range: {smri_normalized[area_cols[0]].min():.6f} - {smri_normalized[area_cols[0]].max():.6f}")
#check the normalization worked for volume
print(f"\nExample volume column before/after normalization:")
print(f"Original range: {smri_all_complete[volume_cols[0]].min():.2f} - {smri_all_complete[volume_cols[0]].max():.2f}")
print(f"Normalized range: {smri_normalized[volume_cols[0]].min():.6f} - {smri_normalized[volume_cols[0]].max():.6f}")
#check the normalization worked for subcortical volume
print(f"\nExample subcortical volume column before/after normalization:")
print(f"Original range: {smri_all_complete[subcvolume_cols[0]].min():.2f} - {smri_all_complete[subcvolume_cols[0]].max():.2f}")
print(f"Normalized range: {smri_normalized[subcvolume_cols[0]].min():.6f} - {smri_normalized[subcvolume_cols[0]].max():.6f}")
# %%
#define ML configuration
ML_CONFIG = {
'rep': 100, # number of repetitions
'k': 3, # number of CV folds
'train_size': 0.80, # training set proportion
'test_size': 0.20, # test set proportion
'max_iter': 1000000,
'solver': 'lsqr',
'alphas': [0, 0.00001, 0.0001, 0.001, 0.004, 0.007, 0.01, 0.04, 0.07, 0.1,
0.4, 0.7, 1, 1.5, 2, 2.5, 3, 3.5, 4, 5, 10, 15, 20, 30, 40, 50,
60, 70, 80, 100]
}
#define how many null models to do
n_null = 1000
print(f"ML Configuration loaded:")
print(f" Repetitions: {ML_CONFIG['rep']}")
print(f" CV folds: {ML_CONFIG['k']}")
print(f" Train/test split: {ML_CONFIG['train_size']:.1%}/{ML_CONFIG['test_size']:.1%}")
print(f" Alpha values: {len(ML_CONFIG['alphas'])} values")
# %%
def train_test_ridge_model(X, Y, site, beh_names, brain_type, group_name, timepoint,
all_results, results_dir, config=ML_CONFIG):
"""
Train ans Test ridge regression models with nested CV, apply Haufe transformation
Parameters:
-----------
X : array
Brain features (n_samples x n_features)
Y : array
Behavioral targets (n_samples x n_behaviors)
site : array
Site IDs for grouping
beh_names : list
Names of behavioral variables
brain_type : str
Type of brain data ('thickness', 'area', 'volume')
group_name : str
Group identifier ('all', 'male', 'female')
timepoint : str
Timepoint identifier
all_results : dict
Master results dictionary to update
results_dir : str
Directory to save results
config : dict
ML configuration parameters
Returns:
--------
dict : Results dictionary with all metrics
"""
rep = config['rep']
k = config['k']
train_size = config['train_size']
alphas = config['alphas']
param_grid = {'alpha': alphas}
#get dimensions of X and Y data
n_feat = X.shape[1]
n_beh = Y.shape[1]
print(f"\nTraining {brain_type} → MH model")
print(f"Group: {group_name}, Timepoint: {timepoint}")
print(f"Samples: {X.shape[0]}, Features: {n_feat}, Behaviors: {n_beh}")
#initialize result arrays
results = {
'r2': np.zeros([n_beh, rep]),
'var': np.zeros([n_beh, rep]),
'corr': np.zeros([n_beh, rep]),
'opt_alpha': np.zeros([n_beh, rep]),
'featimp': np.zeros([n_beh, rep, n_feat]),
'featimp_haufe': np.zeros([n_beh, rep, n_feat]),
'config': {
'brain_type': brain_type,
'group_name': group_name,
'timepoint': timepoint,
'n_samples': X.shape[0],
'n_features': n_feat,
'n_behaviors': n_beh,
'behavior_names': beh_names
}
}
#train model for each behavior
for b in range(n_beh):
beh_data = Y[:, b]
print(f"\n Training {beh_names[b]} ({b+1}/{n_beh})")
#progress bar for repetitions
pbar = tqdm(range(rep), desc=f" Reps",
bar_format='{l_bar}{bar}| {n_fmt}/{total_fmt} [{elapsed}<{remaining}]')
#repetitions loop with progress bar
for p in pbar:
# Train/test split
train_inds, test_inds = next(
GroupShuffleSplit(test_size=1-train_size, n_splits=1, random_state=p)
.split(X, groups=site)
)
#split data
x_train, x_test = X[train_inds], X[test_inds]
y_train, y_test = beh_data[train_inds], beh_data[test_inds]
site_train = site[train_inds]
#convert to double
y_train, y_test = np.double(y_train), np.double(y_test)
#nested CV for hyperparameter optimization
regr = Ridge(normalize=True, max_iter=config['max_iter'], solver=config['solver'])
cv_split = GroupKFold(n_splits=k)
grid_search = GridSearchCV(
estimator=regr,
param_grid=param_grid,
cv=cv_split,
scoring='explained_variance',
n_jobs=-1,
verbose=0
)
#fit and get best parameters
grid_search.fit(x_train, y_train, groups=site_train)
best_alpha = grid_search.best_params_['alpha']
results['opt_alpha'][b, p] = best_alpha
#train final model with optimized alpha
model = Ridge(alpha=best_alpha, normalize=True,
max_iter=config['max_iter'], solver=config['solver'])
model.fit(x_train, y_train)
#evaluate model
results['r2'][b, p] = model.score(x_test, y_test)
preds = model.predict(x_test).ravel()
results['var'][b, p] = explained_variance_score(y_test, preds)
corr_val = np.corrcoef(y_test.ravel(), preds)[1, 0]
results['corr'][b, p] = corr_val
#feature importance (original and Haufe-transformed)
results['featimp'][b, p, :] = model.coef_
#haufe transformation
y_train_preds = model.predict(x_train).ravel()
cov_x = np.cov(x_train.T)
cov_y = np.cov(y_train_preds)
results['featimp_haufe'][b, p, :] = np.matmul(cov_x, model.coef_) * (1/cov_y)
#update progress bar with current correlation
pbar.set_postfix({'corr': f'{corr_val:.3f}'})
pbar.close()
#print summary for this behavior
mean_corr = np.mean(results['corr'][b, :])
std_corr = np.std(results['corr'][b, :])
print(f" Final: {mean_corr:.3f} ± {std_corr:.3f}")
#create result key and add to master results
result_key = f"{timepoint}_{brain_type}_{group_name}"
all_results['models'][result_key] = results
#add to summary
all_results['summary'][result_key] = {
'mean_correlation': float(np.mean(results['corr'])),
'std_correlation': float(np.std(results['corr'])),
'mean_r2': float(np.mean(results['r2'])),
'std_r2': float(np.std(results['r2'])),
'mean_explained_var': float(np.mean(results['var'])),
'std_explained_var': float(np.std(results['var'])),
}
#save incremental backup
backup_filename = f"{results_dir}/all_results_backup.pkl"
with open(backup_filename, 'wb') as f:
pickle.dump(all_results, f)
print(f"✓ Saved to master results: {result_key}")
print(f"✓ Mean correlation: {np.mean(results['corr']):.3f}")
return results
# %%
def generate_null_models(X, Y, site, beh_names, brain_type, group_name, timepoint,
real_results, all_results, results_dir, n_null, config=ML_CONFIG):
"""
Generate null models with permuted behavioral data using random alphas from real models
Parameters:
-----------
X : array
Brain features (n_samples x n_features)
Y : array
Behavioral targets (n_samples x n_behaviors)
site : array
Site IDs for grouping
beh_names : list
Names of behavioral variables
brain_type : str
Type of brain data ('thickness', 'area', 'volume')
group_name : str
Group identifier ('all', 'male', 'female')
timepoint : str
Timepoint identifier
real_results : dict
Results from real models (contains optimized alphas)
all_results : dict
Master results dictionary to update
results_dir : str
Directory to save results
n_null : int
Number of null permutations to generate
config : dict
ML configuration parameters (default: ML_CONFIG)
Returns:
--------
dict : Null results dictionary
"""
train_size = config['train_size']
#get dimensions of X and Y data
n_feat = X.shape[1]
n_beh = Y.shape[1]
print(f"\nGenerating {n_null} NULL {brain_type} → MH models")
print(f"Group: {group_name}, Timepoint: {timepoint}")
print(f"Samples: {X.shape[0]}, Features: {n_feat}, Behaviors: {n_beh}")
#initialize result arrays
null_results = {
'r2': np.zeros([n_beh, n_null]),
'var': np.zeros([n_beh, n_null]),
'corr': np.zeros([n_beh, n_null]),
'alpha_used': np.zeros([n_beh, n_null]),
'config': {
'brain_type': brain_type,
'group_name': group_name,
'timepoint': timepoint,
'n_samples': X.shape[0],
'n_features': n_feat,
'n_behaviors': n_beh,
'behavior_names': beh_names,
'n_null': n_null,
'model_type': 'null'
}
}
#generate null models for each behavior
for b in range(n_beh):
beh_data = Y[:, b]
print(f"\n Generating NULL {beh_names[b]} ({b+1}/{n_beh})")
#get optimized alphas from real model for this behavior
if real_results is not None and 'opt_alpha' in real_results:
available_alphas = real_results['opt_alpha'][b, :]
else:
available_alphas = np.array([1.0])
print(f" Warning: No real model alphas available, using default")
#permute behavioral data while preserving site structure
unique_sites = np.unique(site)
permuted_beh = np.zeros_like(beh_data)
for site_id in unique_sites:
site_mask = site == site_id
site_indices = np.where(site_mask)[0]
permuted_indices = np.random.permutation(site_indices)
permuted_beh[site_indices] = beh_data[permuted_indices]
#null models loop
for p in range(n_null):
# Train/test split (same random state pattern as real models)
train_inds, test_inds = next(
GroupShuffleSplit(test_size=1-train_size, n_splits=1, random_state=p)
.split(X, groups=site)
)
#split data
x_train, x_test = X[train_inds], X[test_inds]
y_train, y_test = permuted_beh[train_inds], permuted_beh[test_inds]
#convert to double
y_train, y_test = np.double(y_train), np.double(y_test)
#randomly select alpha from real model's optimized alphas
random_alpha = np.random.choice(available_alphas)
null_results['alpha_used'][b, p] = random_alpha
#train null model with random alpha
model = Ridge(alpha=random_alpha, normalize=True,
max_iter=config['max_iter'], solver=config['solver'])
model.fit(x_train, y_train)
#evaluate model
null_results['r2'][b, p] = model.score(x_test, y_test)
preds = model.predict(x_test).ravel()
null_results['var'][b, p] = explained_variance_score(y_test, preds)
corr_val = np.corrcoef(y_test.ravel(), preds)[1, 0]
null_results['corr'][b, p] = corr_val
#print summary for this behavior
mean_corr = np.mean(null_results['corr'][b, :])
std_corr = np.std(null_results['corr'][b, :])
print(f" NULL Final: {mean_corr:.3f} ± {std_corr:.3f}")
#create result key and add to master results
null_key = f"{timepoint}_{brain_type}_{group_name}"
all_results['null_models'][null_key] = null_results
#add to null summary
all_results['null_summary'][null_key] = {
'mean_null_correlation': float(np.mean(null_results['corr'])),
'std_null_correlation': float(np.std(null_results['corr'])),
'mean_null_r2': float(np.mean(null_results['r2'])),
'std_null_r2': float(np.std(null_results['r2'])),
'mean_null_var': float(np.mean(null_results['var'])),
'std_null_var': float(np.std(null_results['var'])),
'n_null': n_null
}
#save incremental backup
backup_filename = f"{results_dir}/null_results_backup.pkl"
with open(backup_filename, 'wb') as f:
pickle.dump(all_results, f)
print(f"✓ Saved NULL to master results: {null_key}")
print(f"✓ Mean NULL correlation: {np.mean(null_results['corr']):.3f}")
return null_results
# %%
# gt unique timepoints and define brain types
timepoints = sorted(smri_normalized['session_id'].unique())
sexes = ['all', 'female', 'male']
brain_types = {
'thickness': thickness_cols,
'area': area_cols,
'volume': volume_cols,
'subcvolume': subcvolume_cols
}
final_filename = f"{results_dir}/all_results_complete.pkl"
if os.path.exists(final_filename):
print(f"Loading existing results from {final_filename}")
with open(final_filename, 'rb') as f:
all_results = pickle.load(f)
print(f"Loaded {len(all_results['models'])} existing models")
print(f"Loaded {len(all_results['null_models'])} existing null models")
else:
print("Creating new results dictionary")
#master dictionary to store all results (including null models)
all_results = {
'models': {}, #real model results
'null_models': {}, #null model results
'metadata': { #analysis info
'ml_config': ML_CONFIG,
'mh_behaviors': mh_cols,
'feature_columns': {
'thickness': thickness_cols,
'area': area_cols,
'volume': volume_cols,
'subcvolume': subcvolume_cols
},
'timepoints': [],
'brain_types': ['thickness', 'area', 'volume', 'subcvolume'],
'sex_groups': ['all', 'female', 'male']
},
'summary': {}, #real model summaries
'null_summary': {} #null model summaries
}
print("Starting ML Pipeline...")
print(f"Timepoints: {timepoints}")
print(f"Brain types: {list(brain_types.keys())}")
print(f"Groups: {sexes}")
print(f"Results will be saved to: {results_dir}")
print("="*60)
#main loops
for tp in timepoints:
print(f"\n{'='*20} TIMEPOINT {tp} {'='*20}")
#add timepoint to metadata if not already there
if tp not in all_results['metadata']['timepoints']:
all_results['metadata']['timepoints'].append(tp)
#get data for this timepoint
smri_tp = smri_normalized[smri_normalized['session_id'] == tp].copy()
mh_tp = mh_complete[mh_complete['session_id'] == tp].copy()
gen_tp = gen_mh[gen_mh['session_id'] == tp].copy() # Get site and demographics
#merge all datasets
merged_tp = smri_tp.merge(mh_tp, on=['participant_id', 'session_id'])
merged_tp = merged_tp.merge(gen_tp[['participant_id', 'session_id', 'ab_g_dyn__design_site', 'ab_g_stc__cohort_sex']],
on=['participant_id', 'session_id'])
if len(merged_tp) == 0:
print(f"No data for timepoint {tp}, skipping...")
continue
print(f"Total participants at timepoint {tp}: {len(merged_tp)}")
#get site information
site = merged_tp['ab_g_dyn__design_site'].values
#prepare behavioral data
Y = merged_tp[mh_cols].values
#loop through sexes
for sex in sexes:
print(f"\n--- GROUP: {sex.upper()} ---")
#filter by sex
if sex == 'all':
data_subset = merged_tp.copy()
Y_subset = Y.copy()
site_subset = site.copy()
elif sex == 'male':
sex_mask = merged_tp['ab_g_stc__cohort_sex'] == 1 #'1' for male
data_subset = merged_tp[sex_mask].copy()
Y_subset = Y[sex_mask]
site_subset = site[sex_mask]
elif sex == 'female':
sex_mask = merged_tp['ab_g_stc__cohort_sex'] == 2 #'2' for female
data_subset = merged_tp[sex_mask].copy()
Y_subset = Y[sex_mask]
site_subset = site[sex_mask]
if len(data_subset) < 50: # Minimum sample size
print(f"Skipping {sex} group (n={len(data_subset)}, too small)")
continue
print(f"{sex} group: {len(data_subset)} participants")
print(f"Number of sites: {len(np.unique(site_subset))}")
#loop through brain types
#loop through brain types
for brain_type, feature_cols in brain_types.items():
result_key = f"{tp}_{brain_type}_{sex}"
# Check if this model already exists
if result_key in all_results['models']:
print(f" {result_key} already exists, skipping...")
continue
print(f"\nTraining {brain_type} model for {sex} at timepoint {tp}...")
#prepare brain data for this type
X_brain = data_subset[feature_cols].values
# ==========================================
# REAL MODEL
# ==========================================
print(" Training REAL model...")
results = train_test_ridge_model(
X_brain, Y_subset, site_subset, mh_cols,
brain_type, sex, tp, all_results, results_dir
)
# ==========================================
# NULL MODEL (shuffled Y data)
# ==========================================
print(" Training NULL model...")
# Generate null models using same structure
null_results = generate_null_models(
X_brain, Y_subset, site_subset, mh_cols,
brain_type, sex, tp, results, all_results, results_dir,
config=ML_CONFIG, n_null=n_null
)
result_key = f"{tp}_{brain_type}_{sex}"
print(f" Completed: {result_key}")
print(f" Real mean correlation: {np.mean(results['corr']):.3f}")
print(f" Null mean correlation: {np.mean(null_results['corr']):.3f}")
print(f" Improvement: {np.mean(results['corr']) - np.mean(null_results['corr']):.3f}")
print("-" * 40)
#final comprehensive save
print(f"\n{'='*50}")
print("SAVING FINAL RESULTS...")
#add final metadata
all_results['metadata']['total_models'] = len(all_results['models'])
all_results['metadata']['total_null_models'] = len(all_results['null_models'])
#save final complete results
final_filename = f"{results_dir}/all_results_complete.pkl"
with open(final_filename, 'wb') as f:
pickle.dump(all_results, f)
#save just the summary for quick access
summary_filename = f"{results_dir}/summary_only.pkl"
with open(summary_filename, 'wb') as f:
pickle.dump(all_results['summary'], f)
#save null summary
null_summary_filename = f"{results_dir}/null_summary_only.pkl"
with open(null_summary_filename, 'wb') as f:
pickle.dump(all_results['null_summary'], f)
#save just metadata
metadata_filename = f"{results_dir}/metadata_only.pkl"
with open(metadata_filename, 'wb') as f:
pickle.dump(all_results['metadata'], f)
print("\nML Pipeline completed!")
print("🎉 FINAL RESULTS SAVED!")
print(f"📁 Directory: {results_dir}/")
print(f"📊 Total real models: {len(all_results['models'])}")
print(f"📊 Total null models: {len(all_results['null_models'])}")
print("📋 Files saved:")
print(" - all_results_complete.pkl")
print(" - summary_only.pkl")
print(" - null_summary_only.pkl")
print(" - metadata_only.pkl")
print(" - all_results_backup.pkl (latest backup)")
print(f"📊 Models in final results:")
for key in sorted(all_results['models'].keys()):
print(f" ✓ {key}")
# %%
def generate_brain_predictions_from_existing_models(all_results, merged_data_dict, config):
"""
Generate brain predictions using the same methodology as existing models
but store predictions and test indices for MLM analysis
"""
brain_predictions_dict = {}
# Get metadata
mh_cols = all_results['metadata']['mh_behaviors']
feature_columns = all_results['metadata']['feature_columns']
for result_key, brain_results in all_results['models'].items():
tp, brain_type, sex = result_key.split('_')
print(f"Generating predictions for {result_key}...")
if tp not in merged_data_dict:
continue
merged_tp = merged_data_dict[tp]
feature_cols = feature_columns[brain_type]
# Filter by sex
if sex == 'all':
data_subset = merged_tp.copy()
elif sex == 'male':
sex_mask = merged_tp['ab_g_stc__cohort_sex'] == 1
data_subset = merged_tp[sex_mask].copy()
elif sex == 'female':
sex_mask = merged_tp['ab_g_stc__cohort_sex'] == 2
data_subset = merged_tp[sex_mask].copy()
if len(data_subset) < 50:
continue
# Prepare data
X_brain = data_subset[feature_cols].values
Y = data_subset[mh_cols].values
site = data_subset['ab_g_dyn__design_site'].values
# Storage for all iterations
all_predictions = []
all_true_values = []
all_test_indices = []
all_test_data = []
# Generate predictions using same methodology as original models
for p in range(config['rep']):
# Same train/test split as original
train_inds, test_inds = next(
GroupShuffleSplit(test_size=config['test_size'], n_splits=1, random_state=p)
.split(X_brain, groups=site)
)
X_train, X_test = X_brain[train_inds], X_brain[test_inds]
Y_train, Y_test = Y[train_inds], Y[test_inds]
site_train = site[train_inds]
test_data = data_subset.iloc[test_inds].copy()
# Get optimized alphas from original results
opt_alphas = brain_results['opt_alpha'][:, p] # alphas for this iteration
iteration_predictions = []
# Train model for each behavior
for b, behavior in enumerate(mh_cols):
y_train_beh = Y_train[:, b]
best_alpha = opt_alphas[b]
# Train Ridge model with optimized alpha
model = Ridge(alpha=best_alpha, normalize=True,
max_iter=config['max_iter'], solver=config['solver'])
model.fit(X_train, y_train_beh)
# Generate predictions
predictions = model.predict(X_test)
iteration_predictions.append(predictions)
# Store results
iteration_predictions = np.column_stack(iteration_predictions)
all_predictions.append(iteration_predictions)
all_true_values.append(Y_test)
all_test_indices.append(test_inds)
all_test_data.append(test_data)
brain_predictions_dict[result_key] = {
'predictions': all_predictions, # List of arrays: [iter][sample, behavior]
'true_values': all_true_values, # List of arrays: [iter][sample, behavior]
'test_indices': all_test_indices, # List of arrays: [iter][sample]
'test_data': all_test_data, # List of DataFrames: [iter]
'behaviors': mh_cols
}
print(f" Generated {len(all_predictions)} iterations of predictions")
return brain_predictions_dict
# %%
def create_merged_data_dict():
"""
Create merged_data_dict using the same data preparation as your ML pipeline
"""
print("Creating merged_data_dict...")
# Get unique timepoints (same as your original code)
timepoints = sorted(smri_normalized['session_id'].unique())
merged_data_dict = {}
for tp in timepoints:
print(f"Processing timepoint {tp}...")
# Get data for this timepoint (same as your original code)
smri_tp = smri_normalized[smri_normalized['session_id'] == tp].copy()
mh_tp = mh_complete[mh_complete['session_id'] == tp].copy()
gen_tp = gen_mh[gen_mh['session_id'] == tp].copy()
# Merge all datasets (same as your original code)
merged_tp = smri_tp.merge(mh_tp, on=['participant_id', 'session_id'])
merged_tp = merged_tp.merge(gen_tp[['participant_id', 'session_id', 'ab_g_dyn__design_site', 'ab_g_stc__cohort_sex']],
on=['participant_id', 'session_id'])
if len(merged_tp) == 0:
print(f"No data for timepoint {tp}, skipping...")
continue
print(f" {tp}: {len(merged_tp)} participants")
# Store in dictionary
merged_data_dict[tp] = merged_tp
print(f"Created merged_data_dict with {len(merged_data_dict)} timepoints")
return merged_data_dict
# Create the merged_data_dict
merged_data_dict = create_merged_data_dict()
# Verify the structure
print("\nMerged data structure:")
for tp, data in merged_data_dict.items():
print(f"{tp}: {data.shape}")
print(f" Columns include: {data.columns.tolist()[:10]}...") # Show first 10 columns
# Now you can use it with your function
if 'all_results' in locals():
brain_predictions_dict = generate_brain_predictions_from_existing_models(
all_results, merged_data_dict, ML_CONFIG
)
# %%
#add the significance testing function
def get_null_p(x, null):
"""Compare performance between a model and its null"""
pval = 1 - np.mean(x - null >= 0)
return pval
# %%
# ==========================================
# SIGNIFICANCE TESTING WITH BISBAS AND UPPS SUBSET CORRECTION
# ==========================================
print(f"\n{'='*50}")
print("TESTING MODEL SIGNIFICANCE WITH BISBAS AND UPPS SUBSET CORRECTION...")
print(f"{'='*50}")
# Load results if not already in memory
if 'all_results' not in locals():
with open(f"{results_dir}/all_results_complete.pkl", 'rb') as f:
all_results = pickle.load(f)
# ==========================================
# DEFINE BEHAVIORAL SUBSETS FOR BISBAS AND UPPS
# ==========================================
# Define your behavioral subsets
behavioral_subsets = {
'bisbas': [
'mh_y_bisbas__bis_sum',
'mh_y_bisbas__bas__rr_sum',
'mh_y_bisbas__bas__dr_sum',
'mh_y_bisbas__bas__fs_sum'
], # 4 BIS/BAS behaviors
'upps': [
'mh_y_upps__nurg_sum',
'mh_y_upps__purg_sum',
'mh_y_upps__pers_sum',
'mh_y_upps__plan_sum',
'mh_y_upps__sens_sum'
] # 5 UPPS behaviors
}
# Print subset information
print(f"Behavioral subsets defined:")
for subset_name, behaviors in behavioral_subsets.items():
print(f" {subset_name.upper()}: {len(behaviors)} behaviors")
for beh in behaviors:
print(f" - {beh}")
print(f"\nTotal behaviors: {sum(len(behaviors) for behaviors in behavioral_subsets.values())}")
print(f"Expected total: {len(mh_cols)}")
# Verify all behaviors are included
all_subset_behaviors = []
for behaviors in behavioral_subsets.values():
all_subset_behaviors.extend(behaviors)
missing_behaviors = [beh for beh in mh_cols if beh not in all_subset_behaviors]
extra_behaviors = [beh for beh in all_subset_behaviors if beh not in mh_cols]
if missing_behaviors:
print(f"⚠️ Missing behaviors from subsets: {missing_behaviors}")
if extra_behaviors:
print(f"⚠️ Extra behaviors in subsets: {extra_behaviors}")
if not missing_behaviors and not extra_behaviors:
print(f"✅ All behaviors correctly assigned to subsets")
model_keys = list(all_results['models'].keys())
# Group model keys by timepoint and sex
timepoint_sex_groups = {}
for model_key in model_keys:
parts = model_key.split('_')
timepoint = parts[0] # e.g., 'ses-00A'
sex_group = parts[2] # e.g., 'all', 'male', 'female'
group_key = f"{timepoint}_{sex_group}"
if group_key not in timepoint_sex_groups:
timepoint_sex_groups[group_key] = []
timepoint_sex_groups[group_key].append(model_key)
print(f"\nTimepoint-Sex groups found: {list(timepoint_sex_groups.keys())}")
for group_key, keys in timepoint_sex_groups.items():
print(f" {group_key}: {len(keys)} models")
# Initialize dictionaries to store significance results
significance_results = {}
fdr_results = {}
group_summaries = {}
# ==========================================
# PROCESS EACH TIMEPOINT-SEX GROUP SEPARATELY
# ==========================================
for group_key, group_model_keys in timepoint_sex_groups.items():
timepoint, sex_group = group_key.split('_', 1)
print(f"\n{'='*40}")
print(f"PROCESSING: {timepoint.upper()} - {sex_group.upper()} GROUP")
print(f"{'='*40}")
# Initialize storage for this group
group_significance_results = {}
group_fdr_results = {}
# ==========================================
# PROCESS EACH BEHAVIORAL SUBSET SEPARATELY
# ==========================================
for subset_name, subset_behaviors in behavioral_subsets.items():
print(f"\n--- Processing behavioral subset: {subset_name.upper()} ---")
# Extract data for this timepoint-sex group and behavioral subset
subset_real_corrs = []
subset_null_corrs = []
subset_model_info = []
subset_significance_results = {}
for model_key in group_model_keys:
real_model = all_results['models'][model_key]
null_model = all_results['null_models'][model_key]
behavior_names = real_model['config']['behavior_names']
# Find indices of behaviors in this subset
subset_indices = []
subset_behavior_names = []
for i, behavior in enumerate(behavior_names):
if behavior in subset_behaviors:
subset_indices.append(i)
subset_behavior_names.append(behavior)
if len(subset_indices) == 0:
print(f" No behaviors from {subset_name.upper()} found in {model_key}")
continue
print(f" {model_key}: {len(subset_indices)} behaviors from {subset_name.upper()}")
# Extract correlations for subset behaviors only
real_corrs = np.mean(real_model['corr'], axis=1)[subset_indices] # Shape: (n_subset_behaviors,)
null_corrs = null_model['corr'][subset_indices, :] # Shape: (n_subset_behaviors, n_reps)
# Store for FDR correction within this subset
subset_real_corrs.extend(real_corrs)
subset_null_corrs.extend(null_corrs.T)
# Store model info for each behavior in subset
for i, behavior_idx in enumerate(subset_indices):
subset_model_info.append({
'model_key': model_key,
'behavior_idx': behavior_idx,
'behavior_name': behavior_names[behavior_idx],
'subset_name': subset_name
})
# Calculate p-values for this model's subset behaviors
p_vals = np.zeros(len(subset_indices))
for i, behavior_idx in enumerate(subset_indices):
real_mean = real_corrs[i]
null_dist = null_corrs[i, :]
p_vals[i] = get_null_p(real_mean, null_dist)
# Store results for this model and subset
subset_key = f"{model_key}_{subset_name}"
subset_significance_results[subset_key] = {
'p_values': p_vals,
'real_correlations': real_corrs,
'null_correlations': null_corrs,
'behavior_names': subset_behavior_names,
'behavior_indices': subset_indices,
'subset_name': subset_name
}
if len(subset_real_corrs) == 0:
print(f" No data found for subset {subset_name.upper()} in {group_key}")
continue
print(f" Collected {len(subset_real_corrs)} tests for {subset_name.upper()} in {group_key}")
# ==========================================
# FDR CORRECTION WITHIN THIS BEHAVIORAL SUBSET
# ==========================================
print(f" Applying FDR correction within {subset_name.upper()}...")
# Collect all p-values for this subset
subset_all_p_vals = []
for subset_key in subset_significance_results.keys():
if subset_name in subset_key:
subset_all_p_vals.extend(subset_significance_results[subset_key]['p_values'])
subset_all_p_vals = np.array(subset_all_p_vals)
if len(subset_all_p_vals) == 0:
continue
# Apply FDR correction within this subset
subset_fdr_rejected, subset_fdr_p_corrected = fdr(subset_all_p_vals, alpha=0.05, method='indep')
print(f" Tests in {subset_name.upper()}: {len(subset_all_p_vals)}")
print(f" Significant after FDR: {np.sum(subset_fdr_rejected)}")
print(f" Significance rate: {np.mean(subset_fdr_rejected):.3f}")
# ==========================================
# ORGANIZE FDR RESULTS BY MODEL WITHIN SUBSET
# ==========================================
idx = 0
for model_key in group_model_keys:
subset_key = f"{model_key}_{subset_name}"
if subset_key not in subset_significance_results:
continue
n_beh = len(subset_significance_results[subset_key]['p_values'])
# Extract FDR results for this model and subset
model_fdr_rejected = subset_fdr_rejected[idx:idx+n_beh]
model_fdr_p_corrected = subset_fdr_p_corrected[idx:idx+n_beh]
# Store FDR results
group_fdr_results[subset_key] = {
'fdr_rejected': model_fdr_rejected,
'fdr_p_corrected': model_fdr_p_corrected,
'n_significant': np.sum(model_fdr_rejected),
'real_correlations': subset_significance_results[subset_key]['real_correlations'],
'p_values_uncorrected': subset_significance_results[subset_key]['p_values'],
'behavior_names': subset_significance_results[subset_key]['behavior_names'],
'behavior_indices': subset_significance_results[subset_key]['behavior_indices'],
'subset_name': subset_name
}
idx += n_beh
# Store results for this timepoint-sex group
significance_results.update(subset_significance_results)
fdr_results.update(group_fdr_results)
# Create group summary across all subsets
group_total_tests = 0
group_total_significant = 0
subset_summaries = {}
for subset_name in behavioral_subsets.keys():
subset_tests = 0
subset_significant = 0
for key in group_fdr_results.keys():
if f"_{subset_name}" in key and group_key in key:
subset_tests += len(group_fdr_results[key]['fdr_rejected'])
subset_significant += group_fdr_results[key]['n_significant']
if subset_tests > 0:
subset_summaries[subset_name] = {
'tests': subset_tests,
'significant': subset_significant,
'rate': subset_significant / subset_tests
}
group_total_tests += subset_tests
group_total_significant += subset_significant
group_summaries[group_key] = {
'timepoint': timepoint,
'sex_group': sex_group,
'total_tests': group_total_tests,
'total_significant': group_total_significant,
'significance_rate': float(group_total_significant/group_total_tests) if group_total_tests > 0 else 0.0,
'subset_summaries': subset_summaries,
'models': group_model_keys,
'fdr_alpha': 0.05,
'fdr_method': 'indep'
}
# Print summary for this group
print(f"\n{group_key} SUMMARY:")
print(f" Total tests across subsets: {group_total_tests}")
print(f" Total significant: {group_total_significant}")
if group_total_tests > 0:
print(f" Overall rate: {group_total_significant/group_total_tests:.3f}")
else:
print(f" Overall rate: 0.000")
for subset_name, subset_summary in subset_summaries.items():
print(f" {subset_name.upper()}: {subset_summary['significant']}/{subset_summary['tests']} ({subset_summary['rate']:.3f})")
# Show significant models for this group
print(f" Significant models by subset:")
for subset_name in behavioral_subsets.keys():
print(f" {subset_name.upper()}:")
for model_key in group_model_keys:
subset_key = f"{model_key}_{subset_name}"
if subset_key in group_fdr_results:
n_sig = group_fdr_results[subset_key]['n_significant']
n_total = len(group_fdr_results[subset_key]['fdr_rejected'])
if n_sig > 0:
print(f" {model_key}: {n_sig}/{n_total} ✓")
else:
print(f" {model_key}: {n_sig}/{n_total}")
# ==========================================
# OVERALL SUMMARY BY TIMEPOINT, SEX, AND BEHAVIORAL SUBSET
# ==========================================
print(f"\n{'='*60}")
print("SIGNIFICANCE SUMMARY BY TIMEPOINT, SEX, AND BEHAVIORAL SUBSET:")
print("="*60)
# Create organized summary
timepoint_summary = {}
total_tests_all = 0
total_significant_all = 0
for group_key in sorted(group_summaries.keys()):
summary = group_summaries[group_key]
timepoint = summary['timepoint']
sex_group = summary['sex_group']
if timepoint not in timepoint_summary:
timepoint_summary[timepoint] = {}
timepoint_summary[timepoint][sex_group] = summary
total_tests_all += summary['total_tests']
total_significant_all += summary['total_significant']
# Print organized summary
for timepoint in sorted(timepoint_summary.keys()):
print(f"\n{timepoint}:")
for sex_group in ['all', 'female', 'male']:
if sex_group in timepoint_summary[timepoint]:
summary = timepoint_summary[timepoint][sex_group]
print(f" {sex_group:>6}: {summary['total_significant']:>3}/{summary['total_tests']:<3} ({summary['significance_rate']:.3f})")
# Show subset breakdown
for subset_name, subset_summary in summary['subset_summaries'].items():
print(f" {subset_name.upper():>8}: {subset_summary['significant']:>2}/{subset_summary['tests']:<2} ({subset_summary['rate']:.3f})")
else:
print(f" {sex_group:>6}: No data")
print(f"\nGRAND TOTAL:")
print(f" Total tests: {total_tests_all}")
print(f" Total significant: {total_significant_all}")
if total_tests_all > 0:
print(f" Overall rate: {total_significant_all/total_tests_all:.3f}")
else:
print(f" Overall rate: 0.000")
# Summary by subset across all timepoints and sexes
print(f"\nSUMMARY BY BEHAVIORAL SUBSET:")
subset_totals = {}
for subset_name in behavioral_subsets.keys():
subset_tests = 0
subset_significant = 0
for group_summary in group_summaries.values():
if subset_name in group_summary['subset_summaries']:
subset_tests += group_summary['subset_summaries'][subset_name]['tests']
subset_significant += group_summary['subset_summaries'][subset_name]['significant']
subset_totals[subset_name] = {
'tests': subset_tests,
'significant': subset_significant,
'rate': subset_significant / subset_tests if subset_tests > 0 else 0.0
}
print(f" {subset_name.upper():>8}: {subset_significant:>3}/{subset_tests:<3} ({subset_totals[subset_name]['rate']:.3f})")
# ==========================================
# ADD TO MASTER RESULTS AND SAVE
# ==========================================
print(f"\nAdding BIS/BAS and UPPS subset significance results to master results...")
# Add significance results to all_results
all_results['significance_subset'] = significance_results
all_results['fdr_results_subset'] = fdr_results
all_results['significance_summary_subset'] = {
'behavioral_subsets': behavioral_subsets,
'by_timepoint_and_sex': group_summaries,
'by_timepoint': timepoint_summary,
'by_subset': subset_totals,
'overall': {
'total_tests': total_tests_all,
'total_significant': total_significant_all,
'significance_rate': float(total_significant_all/total_tests_all) if total_tests_all > 0 else 0.0,
'fdr_alpha': 0.05,
'fdr_method': 'indep',
'correction_level': 'within_bisbas_and_upps_subsets'
}
}
# Save updated results
final_filename = f"{results_dir}/all_results_complete_with_bisbas_upps_significance.pkl"
with open(final_filename, 'wb') as f:
pickle.dump(all_results, f)
# Save just subset significance results
subset_significance_filename = f"{results_dir}/bisbas_upps_significance_results_only.pkl"
with open(subset_significance_filename, 'wb') as f:
pickle.dump({
'behavioral_subsets': behavioral_subsets,
'significance': significance_results,
'fdr_results': fdr_results,
'summary': all_results['significance_summary_subset']
}, f)
print(f"\n🎉 BIS/BAS AND UPPS SUBSET SIGNIFICANCE TESTING COMPLETED!")
print(f"📁 Results saved to:")
print(f" - all_results_complete_with_bisbas_upps_significance.pkl")
print(f" - bisbas_upps_significance_results_only.pkl")
# %%
#print out the significant models
def show_significant_models():
"""Show only models with significant results"""
print("Models with significant results:")
print("="*50)
for model_key, results in fdr_results.items():
if results['n_significant'] > 0:
print(f"{model_key}: {results['n_significant']} significant behaviors")
sig_mask = results['fdr_rejected']
for i, is_sig in enumerate(sig_mask):
if is_sig:
behavior = results['behavior_names'][i]
corr = results['real_correlations'][i]
p_corr = results['fdr_p_corrected'][i]
print(f" {behavior}: r={corr:.3f}, p_fdr={p_corr:.4f}")
print()
#run quick summaries
show_significant_models()
# %%
def analyze_model_performance_vs_age_correlation_by_behavior_significant_only(brain_predictions_dict, gen_complete, mh_cols, all_results):
"""
Analyze correlation between model performance and age-behavior correlations
- Separately for each behavior
- Only for models that predicted behaviors better than chance (FDR significant)
For each significant behavior-model combination:
a = correlation between observed and predicted values
b = correlation between observed values and age
Then compute correlation between a and b across iterations for that specific behavior
"""
print("\n🔍 ANALYZING MODEL PERFORMANCE vs AGE-BEHAVIOR CORRELATIONS")
print("BY BEHAVIOR - SIGNIFICANT MODELS ONLY")
print("=" * 80)
# Get significance results
if 'fdr_results_subset' not in all_results:
print("❌ No significance results found. Run significance testing first!")
return None
fdr_results = all_results['fdr_results_subset']
behavioral_subsets = all_results['significance_summary_subset']['behavioral_subsets']
age_analysis_results = {}
total_significant_behaviors = 0
total_analyzed_behaviors = 0
for result_key, pred_data in brain_predictions_dict.items():
print(f"\nAnalyzing {result_key}...")
# Parse model components - result_key format is "timepoint_braintype_sexgroup"
try:
parts = result_key.split('_')
timepoint = parts[0]
brain_type = parts[1]
sex_group = parts[2]
except:
print(f" ⚠️ Could not parse {result_key}, skipping...")
continue
model_key = result_key # This is the model key without behavior name
behaviors = pred_data['behaviors'] # List of all behaviors for this model
print(f" Model: {model_key}")
print(f" Behaviors in model: {len(behaviors)}")
# Process each behavior in this model
for behavior_idx, behavior_name in enumerate(behaviors):
print(f" Checking behavior: {behavior_name}")
# Find which subset this behavior belongs to
behavior_subset = None
for subset_name, subset_behaviors in behavioral_subsets.items():
if behavior_name in subset_behaviors:
behavior_subset = subset_name
break
if behavior_subset is None:
print(f" ⚠️ Behavior {behavior_name} not found in any subset, skipping...")
continue
# Check significance for this specific behavior
subset_key = f"{model_key}_{behavior_subset}"
if subset_key not in fdr_results:
print(f" ⚠️ No significance results for {subset_key}, skipping...")
continue
fdr_result = fdr_results[subset_key]
fdr_behavior_names = fdr_result['behavior_names']
fdr_rejected = fdr_result['fdr_rejected']
# Find index of this behavior in the significance results
try:
fdr_behavior_idx = fdr_behavior_names.index(behavior_name)
is_significant = fdr_rejected[fdr_behavior_idx]
except ValueError:
print(f" ⚠️ Behavior {behavior_name} not found in significance results, skipping...")
continue
if not is_significant:
print(f" ❌ {behavior_name} not significant, skipping...")
continue
print(f" ✅ {behavior_name} is significant, analyzing...")
total_significant_behaviors += 1
n_iterations = len(pred_data['predictions'])
# Storage for correlations across iterations for this specific behavior
behavior_model_performance_corrs = [] # a values
behavior_age_behavior_corrs = [] # b values
# Process each iteration
for iteration in range(n_iterations):
predictions = pred_data['predictions'][iteration] # Shape: (n_samples, n_behaviors)
true_values = pred_data['true_values'][iteration] # Shape: (n_samples, n_behaviors)
test_indices = pred_data['test_indices'][iteration] # Shape: (n_samples,)
# Extract age data for these participants
try:
ages = gen_complete.iloc[test_indices]['ab_g_dyn__visit_age'].values
except:
continue
if len(ages) != len(true_values):
continue
# Extract data for this specific behavior
obs = true_values[:, behavior_idx] # Use behavior_idx from the behaviors list
pred = predictions[:, behavior_idx]
# Remove NaN values
valid_mask = ~(np.isnan(obs) | np.isnan(pred) | np.isnan(ages))
if valid_mask.sum() < 10: # Need minimum samples
continue
obs_clean = obs[valid_mask]
pred_clean = pred[valid_mask]
ages_clean = ages[valid_mask]
# a: model performance correlation for this behavior
try:
model_perf_corr = np.corrcoef(obs_clean, pred_clean)[0, 1]
if np.isnan(model_perf_corr):
continue
except:
continue
# b: age-behavior correlation for this behavior
try:
age_behavior_corr = np.corrcoef(obs_clean, ages_clean)[0, 1]
if np.isnan(age_behavior_corr):
continue
except:
continue
behavior_model_performance_corrs.append(model_perf_corr)
behavior_age_behavior_corrs.append(age_behavior_corr)
# Convert to arrays
behavior_model_performance_corrs = np.array(behavior_model_performance_corrs)
behavior_age_behavior_corrs = np.array(behavior_age_behavior_corrs)
print(f" Collected {len(behavior_model_performance_corrs)} valid iterations")
if len(behavior_model_performance_corrs) < 5:
print(f" ⚠️ Insufficient data for correlation analysis")
continue
# Compute correlation between a and b for this specific behavior
try:
from scipy.stats import pearsonr
performance_age_corr, p_value = pearsonr(behavior_model_performance_corrs, behavior_age_behavior_corrs)
except:
performance_age_corr = np.nan
p_value = np.nan
# Create unique result key for this behavior-model combination
behavior_result_key = f"{model_key}_{behavior_name.replace('mh_y_', '').replace('_', '')}"
# Store results for this specific behavior
age_analysis_results[behavior_result_key] = {
'behavior_name': behavior_name,
'model_key': model_key,
'timepoint': timepoint,
'brain_type': brain_type,
'sex_group': sex_group,
'behavior_subset': behavior_subset,
'behavior_idx': behavior_idx,
'performance_age_correlation': performance_age_corr,
'p_value': p_value,
'n_datapoints': len(behavior_model_performance_corrs),
'model_perf_mean': np.mean(behavior_model_performance_corrs),
'model_perf_std': np.std(behavior_model_performance_corrs),
'age_corr_mean': np.mean(behavior_age_behavior_corrs),
'age_corr_std': np.std(behavior_age_behavior_corrs),
'is_significant_predictor': True # All included behaviors are significant
}
print(f" Performance-age correlation: {performance_age_corr:.3f} (p={p_value:.3f})")
total_analyzed_behaviors += 1
print(f"\n📊 ANALYSIS SUMMARY:")
print(f" Total significant behaviors found: {total_significant_behaviors}")
print(f" Total behaviors analyzed: {total_analyzed_behaviors}")
return age_analysis_results
# %%
def print_significant_age_correlations_by_behavior_fdr_corrected(age_analysis_results):
"""
Print correlation values for significant model-behavior combinations,
with FDR correction across all tests
"""
print(f"\n🎯 SIGNIFICANT PERFORMANCE vs AGE-BEHAVIOR CORRELATIONS")
print("BY INDIVIDUAL BEHAVIOR (SIGNIFICANT PREDICTORS ONLY)")
print("WITH FDR CORRECTION ACROSS ALL TESTS")
print("=" * 80)
if not age_analysis_results:
print("No results to display.")
return
# Collect all results with valid p-values
all_results = []
all_p_values = []
for result_key, results in age_analysis_results.items():
corr = results['performance_age_correlation']
p_val = results['p_value']
# Only include results with valid p-values
if not np.isnan(p_val):
all_results.append({
'result_key': result_key,
'behavior_name': results['behavior_name'],
'timepoint': results['timepoint'],
'brain_type': results['brain_type'],
'sex_group': results['sex_group'],
'behavior_subset': results['behavior_subset'],
'correlation': corr,
'p_value': p_val,
'n_datapoints': results['n_datapoints'],
'model_perf_mean': results['model_perf_mean'],
'age_corr_mean': results['age_corr_mean']
})
all_p_values.append(p_val)
if len(all_results) == 0:
print("No valid results to analyze.")
return []
print(f"Total tests performed: {len(all_results)}")
# Apply FDR correction
from statsmodels.stats.multitest import fdrcorrection
all_p_values = np.array(all_p_values)
fdr_rejected, fdr_p_corrected = fdrcorrection(all_p_values, alpha=0.05, method='indep')
# Add FDR results to each result
for i, result in enumerate(all_results):
result['fdr_rejected'] = fdr_rejected[i]
result['fdr_p_corrected'] = fdr_p_corrected[i]
# Determine significance level for uncorrected p-value
if result['p_value'] < 0.001:
result['uncorrected_significance'] = "***"
elif result['p_value'] < 0.01:
result['uncorrected_significance'] = "**"
elif result['p_value'] < 0.05:
result['uncorrected_significance'] = "*"
else:
result['uncorrected_significance'] = ""
# Determine significance level for FDR corrected p-value
if result['fdr_p_corrected'] < 0.001:
result['fdr_significance'] = "***"
elif result['fdr_p_corrected'] < 0.01:
result['fdr_significance'] = "**"
elif result['fdr_p_corrected'] < 0.05:
result['fdr_significance'] = "*"
else:
result['fdr_significance'] = ""
# Separate significant results (uncorrected vs FDR corrected)
uncorrected_significant = [r for r in all_results if r['p_value'] < 0.05]
fdr_significant = [r for r in all_results if r['fdr_rejected']]
print(f"Significant before FDR correction: {len(uncorrected_significant)}")
print(f"Significant after FDR correction: {len(fdr_significant)}")
print(f"FDR correction alpha: 0.05")
# Print FDR significant results
if len(fdr_significant) > 0:
print(f"\n📊 FDR-CORRECTED SIGNIFICANT RESULTS:")
print("=" * 80)
# Sort by FDR corrected p-value
fdr_significant.sort(key=lambda x: x['fdr_p_corrected'])
# Print header
print(f"{'Behavior':<25} {'Model':<25} {'r':<8} {'p_uncorr':<10} {'p_FDR':<10} {'FDR_sig':<7}")
print("-" * 85)
# Print each FDR significant result
for result in fdr_significant:
behavior_clean = result['behavior_name'].replace('mh_y_', '').replace('_', '')[:24]
timepoint_label = result['timepoint'].replace('ses-', 'Y').replace('A', '')
model_desc = f"{result['brain_type'][:4]}-{result['sex_group'][:3]}-{timepoint_label}"[:24]
print(f"{behavior_clean:<25} {model_desc:<25} {result['correlation']:>7.3f} "
f"{result['p_value']:>9.4f} {result['fdr_p_corrected']:>9.4f} "
f"{result['fdr_significance']:>6s}")
else:
print(f"\n❌ No results survive FDR correction at α = 0.05")
# Print comparison table of uncorrected significant results
if len(uncorrected_significant) > 0:
print(f"\n📋 ALL UNCORRECTED SIGNIFICANT RESULTS (with FDR status):")
print("=" * 100)
# Sort by uncorrected p-value
uncorrected_significant.sort(key=lambda x: x['p_value'])
# Print header
print(f"{'Behavior':<25} {'Model':<25} {'r':<8} {'p_uncorr':<10} {'p_FDR':<10} {'FDR_pass':<8}")
print("-" * 95)
# Print each uncorrected significant result
for result in uncorrected_significant:
behavior_clean = result['behavior_name'].replace('mh_y_', '').replace('_', '')[:24]
timepoint_label = result['timepoint'].replace('ses-', 'Y').replace('A', '')
model_desc = f"{result['brain_type'][:4]}-{result['sex_group'][:3]}-{timepoint_label}"[:24]
fdr_pass = "✓" if result['fdr_rejected'] else "✗"
print(f"{behavior_clean:<25} {model_desc:<25} {result['correlation']:>7.3f} "
f"{result['p_value']:>9.4f} {result['fdr_p_corrected']:>9.4f} "
f"{fdr_pass:>7s}")
# Summary statistics
print(f"\n📊 SUMMARY STATISTICS:")
print("=" * 50)
# Before FDR correction
if len(uncorrected_significant) > 0:
uncorr_corr_values = [r['correlation'] for r in uncorrected_significant]
print(f"\nUNCORRECTED SIGNIFICANT (p < 0.05):")
print(f" Count: {len(uncorrected_significant)}")
print(f" Mean correlation: {np.mean(uncorr_corr_values):6.3f}")
print(f" Range: {np.min(uncorr_corr_values):6.3f} to {np.max(uncorr_corr_values):6.3f}")
# Count by significance level (uncorrected)
uncorr_highly_sig = len([r for r in uncorrected_significant if r['p_value'] < 0.001])
uncorr_very_sig = len([r for r in uncorrected_significant if 0.001 <= r['p_value'] < 0.01])
uncorr_sig = len([r for r in uncorrected_significant if 0.01 <= r['p_value'] < 0.05])
print(f" By uncorrected significance level:")
print(f" p < 0.001: {uncorr_highly_sig}")
print(f" p < 0.01: {uncorr_very_sig}")
print(f" p < 0.05: {uncorr_sig}")
# After FDR correction
if len(fdr_significant) > 0:
fdr_corr_values = [r['correlation'] for r in fdr_significant]
print(f"\nFDR-CORRECTED SIGNIFICANT:")
print(f" Count: {len(fdr_significant)}")
print(f" Mean correlation: {np.mean(fdr_corr_values):6.3f}")
print(f" Range: {np.min(fdr_corr_values):6.3f} to {np.max(fdr_corr_values):6.3f}")
# Count by significance level (FDR corrected)
fdr_highly_sig = len([r for r in fdr_significant if r['fdr_p_corrected'] < 0.001])
fdr_very_sig = len([r for r in fdr_significant if 0.001 <= r['fdr_p_corrected'] < 0.01])
fdr_sig = len([r for r in fdr_significant if 0.01 <= r['fdr_p_corrected'] < 0.05])
print(f" By FDR-corrected significance level:")
print(f" p_FDR < 0.001: {fdr_highly_sig}")
print(f" p_FDR < 0.01: {fdr_very_sig}")
print(f" p_FDR < 0.05: {fdr_sig}")
# Group by categories (FDR significant only)
print(f"\n FDR significant by brain type:")
brain_counts = {}
for result in fdr_significant:
bt = result['brain_type']
brain_counts[bt] = brain_counts.get(bt, 0) + 1
for bt, count in sorted(brain_counts.items()):
print(f" {bt}: {count}")
print(f"\n FDR significant by behavioral subset:")
subset_counts = {}
for result in fdr_significant:
subset = result['behavior_subset']
subset_counts[subset] = subset_counts.get(subset, 0) + 1
for subset, count in sorted(subset_counts.items()):
print(f" {subset.upper()}: {count}")
# Interpretation
pos_corrs = len([r for r in fdr_significant if r['correlation'] > 0])
neg_corrs = len([r for r in fdr_significant if r['correlation'] < 0])
print(f"\n FDR significant interpretation:")
print(f" Positive correlations: {pos_corrs} (better prediction of age-related behaviors)")
print(f" Negative correlations: {neg_corrs} (better prediction of age-independent behaviors)")
else:
print(f"\nNo results survive FDR correction.")
# False discovery rate
if len(all_results) > 0:
estimated_false_discoveries = len(uncorrected_significant) * 0.05 # Expected under null
print(f"\n📈 FALSE DISCOVERY ESTIMATION:")
print(f" Expected false discoveries (5% of {len(all_results)} tests): {estimated_false_discoveries:.1f}")
print(f" Observed uncorrected significant: {len(uncorrected_significant)}")
print(f" FDR-corrected significant: {len(fdr_significant)}")
if len(uncorrected_significant) > 0:
print(f" Estimated FDR in uncorrected results: {estimated_false_discoveries/len(uncorrected_significant):.3f}")
return {
'all_results': all_results,
'uncorrected_significant': uncorrected_significant,
'fdr_significant': fdr_significant,
'fdr_alpha': 0.05,
'n_tests': len(all_results)
}
# Update the main analysis section
if 'brain_predictions_dict' in locals() and brain_predictions_dict is not None:
print("🚀 RUNNING BEHAVIOR-SPECIFIC MODEL PERFORMANCE vs AGE CORRELATION ANALYSIS")
print("SIGNIFICANT PREDICTORS ONLY WITH FDR CORRECTION")
# Run analysis
age_analysis_by_behavior = analyze_model_performance_vs_age_correlation_by_behavior_significant_only(
brain_predictions_dict, gen_complete, mh_cols, all_results
)
if age_analysis_by_behavior:
# Print significant results with FDR correction
fdr_corrected_results = print_significant_age_correlations_by_behavior_fdr_corrected(age_analysis_by_behavior)
print(f"\n🎉 BEHAVIOR-SPECIFIC AGE CORRELATION ANALYSIS COMPLETED!")
print(f"📊 Results available in age_analysis_by_behavior dictionary")
print(f"📊 FDR-corrected results available in fdr_corrected_results dictionary")
else:
print("❌ No behavior-specific results generated")
else:
print("❌ No brain_predictions_dict found. Generate predictions first!")
# %%
# ==========================================
# FIND LONGITUDINAL PARTICIPANTS USING EXISTING VARIABLES
# ==========================================
def find_longitudinal_participants_from_existing(gen_complete, smri_normalized, mh_complete):
"""Find participants with both baseline and 2-year data using existing complete datasets"""
print("🔍 FINDING LONGITUDINAL PARTICIPANTS FROM EXISTING DATASETS")
print("=" * 60)
# Get participants with baseline data
baseline_gen = set(gen_complete[gen_complete['session_id'] == 'ses-00A']['participant_id'])
baseline_smri = set(smri_normalized[smri_normalized['session_id'] == 'ses-00A']['participant_id'])
baseline_mh = set(mh_complete[mh_complete['session_id'] == 'ses-00A']['participant_id'])
print(f"Baseline participants:")
print(f" Demographics: {len(baseline_gen)}")
print(f" Brain (SMRI): {len(baseline_smri)}")
print(f" Behavior (MH): {len(baseline_mh)}")
# Get participants with 4-year data
sixyear_gen = set(gen_complete[gen_complete['session_id'] == 'ses-02A']['participant_id'])
sixyear_smri = set(smri_normalized[smri_normalized['session_id'] == 'ses-02A']['participant_id'])
sixyear_mh = set(mh_complete[mh_complete['session_id'] == 'ses-02A']['participant_id'])
print(f"\n2-year participants:")
print(f" Demographics: {len(sixyear_gen)}")
print(f" Brain (SMRI): {len(sixyear_smri)}")
print(f" Behavior (MH): {len(sixyear_mh)}")
# Find longitudinal participants (have both baseline and 4-year)
longitudinal_gen = baseline_gen.intersection(sixyear_gen)
longitudinal_smri = baseline_smri.intersection(sixyear_smri)
longitudinal_mh = baseline_mh.intersection(sixyear_mh)
print(f"\nLongitudinal participants:")
print(f" Demographics: {len(longitudinal_gen)}")
print(f" Brain (SMRI): {len(longitudinal_smri)}")
print(f" Behavior (MH): {len(longitudinal_mh)}")
# Find complete longitudinal participants (have both brain and behavior longitudinal data)
complete_longitudinal = longitudinal_smri.intersection(longitudinal_mh).intersection(longitudinal_gen)
print(f" Complete (Brain + Behavior + Demographics): {len(complete_longitudinal)}")
return {
'longitudinal_gen': list(longitudinal_gen),
'longitudinal_smri': list(longitudinal_smri),
'longitudinal_mh': list(longitudinal_mh),
'complete_longitudinal': list(complete_longitudinal)
}
# Find longitudinal participants
longitudinal_participants = find_longitudinal_participants_from_existing(
gen_complete, smri_normalized, mh_complete
)
# ==========================================
# CREATE SMRI_DIFF (BRAIN DIFFERENCES)
# ==========================================
def create_smri_diff(smri_normalized, longitudinal_smri_participants):
"""Create smri_diff dataset with baseline to 2-year differences"""
print(f"\n🧠 CREATING SMRI_DIFF DATASET")
print("=" * 40)
# Filter for longitudinal participants only
smri_longitudinal = smri_normalized[
smri_normalized['participant_id'].isin(longitudinal_smri_participants)
].copy()
print(f"SMRI longitudinal data: {len(smri_longitudinal)} rows")
# Separate baseline and 4-year data
smri_baseline = smri_longitudinal[smri_longitudinal['session_id'] == 'ses-00A'].copy()
smri_sixYear = smri_longitudinal[smri_longitudinal['session_id'] == 'ses-02A'].copy()
print(f" Baseline: {len(smri_baseline)} participants")
print(f" 2-year: {len(smri_sixYear)} participants")
# Merge baseline and 4-year data
smri_merged = smri_baseline.merge(
smri_sixYear,
on='participant_id',
suffixes=('_baseline', '_6year')
)
print(f" Merged: {len(smri_merged)} participants")
# Get brain measure columns (exclude ID and session columns)
id_cols = ['participant_id', 'session_id_baseline', 'session_id_6year']
brain_cols = [col for col in smri_baseline.columns
if col not in ['participant_id', 'session_id']]
print(f" Brain measure columns: {len(brain_cols)}")
# Calculate differences (4-year - baseline)
smri_diff_data = {'participant_id': smri_merged['participant_id']}
for col in brain_cols:
baseline_col = f"{col}_baseline"
sixYear_col = f"{col}_6year"
if baseline_col in smri_merged.columns and sixYear_col in smri_merged.columns:
smri_diff_data[col] = smri_merged[sixYear_col] - smri_merged[baseline_col]
else:
print(f" ⚠️ Missing columns for {col}")
# Create smri_diff DataFrame
smri_diff = pd.DataFrame(smri_diff_data)
print(f"✅ Created smri_diff: {smri_diff.shape}")
print(f" Participants: {len(smri_diff)}")
print(f" Brain measures: {len(smri_diff.columns) - 1}") # -1 for participant_id
return smri_diff
# Create smri_diff
smri_diff = create_smri_diff(smri_normalized, longitudinal_participants['longitudinal_smri'])
# ==========================================
# CREATE MH_DIFF (BEHAVIOR DIFFERENCES)
# ==========================================
def create_mh_diff(mh_complete, longitudinal_mh_participants):
"""Create mh_diff dataset with baseline to 2-year differences"""
print(f"\n📈 CREATING MH_DIFF DATASET")
print("=" * 40)
# Filter for longitudinal participants only
mh_longitudinal = mh_complete[
mh_complete['participant_id'].isin(longitudinal_mh_participants)
].copy()
print(f"MH longitudinal data: {len(mh_longitudinal)} rows")
# Separate baseline and 4-year data
mh_baseline = mh_longitudinal[mh_longitudinal['session_id'] == 'ses-00A'].copy()
mh_sixYear = mh_longitudinal[mh_longitudinal['session_id'] == 'ses-02A'].copy()
print(f" Baseline: {len(mh_baseline)} participants")
print(f" 4-year: {len(mh_sixYear)} participants")
# Merge baseline and 4-year data
mh_merged = mh_baseline.merge(
mh_sixYear,
on='participant_id',
suffixes=('_baseline', '_6year')
)
print(f" Merged: {len(mh_merged)} participants")
# Get behavior columns (exclude ID and session columns)
id_cols = ['participant_id', 'session_id_baseline', 'session_id_6year']
behavior_cols = [col for col in mh_baseline.columns
if col not in ['participant_id', 'session_id', 'site_id_l']]
print(f" Behavior columns: {len(behavior_cols)}")
# Calculate differences (4-year - baseline)
mh_diff_data = {'participant_id': mh_merged['participant_id']}
for col in behavior_cols:
baseline_col = f"{col}_baseline"
sixYear_col = f"{col}_6year"
if baseline_col in mh_merged.columns and sixYear_col in mh_merged.columns:
mh_diff_data[col] = mh_merged[sixYear_col] - mh_merged[baseline_col]
else:
print(f" ⚠️ Missing columns for {col}")
# Create mh_diff DataFrame
mh_diff = pd.DataFrame(mh_diff_data)
print(f"✅ Created mh_diff: {mh_diff.shape}")
print(f" Participants: {len(mh_diff)}")
print(f" Behavior measures: {len(mh_diff.columns) - 1}") # -1 for participant_id
return mh_diff
# Create mh_diff
mh_diff = create_mh_diff(mh_complete, longitudinal_participants['longitudinal_mh'])
# ==========================================
# CREATE FILTERED DATASETS FOR COMPLETE LONGITUDINAL PARTICIPANTS
# ==========================================
def create_complete_longitudinal_datasets(smri_diff, mh_diff, gen_complete, complete_longitudinal_participants):
"""Create datasets filtered to participants with complete longitudinal data"""
print(f"\n🔄 CREATING COMPLETE LONGITUDINAL DATASETS")
print("=" * 50)
# Filter smri_diff to complete longitudinal participants
smri_diff_complete = smri_diff[
smri_diff['participant_id'].isin(complete_longitudinal_participants)
].copy()
# Filter mh_diff to complete longitudinal participants
mh_diff_complete = mh_diff[
mh_diff['participant_id'].isin(complete_longitudinal_participants)
].copy()
# Get demographics for complete longitudinal participants
gen_longitudinal = gen_complete[
gen_complete['participant_id'].isin(complete_longitudinal_participants)
].copy()
print(f"Complete longitudinal datasets:")
print(f" smri_diff_complete: {smri_diff_complete.shape}")
print(f" mh_diff_complete: {mh_diff_complete.shape}")
print(f" gen_longitudinal: {gen_longitudinal.shape}")
return smri_diff_complete, mh_diff_complete, gen_longitudinal
# Create complete longitudinal datasets
smri_diff_complete, mh_diff_complete, gen_longitudinal = create_complete_longitudinal_datasets(
smri_diff, mh_diff, gen_complete, longitudinal_participants['complete_longitudinal']
)
# ==========================================
# SUMMARY AND VERIFICATION
# ==========================================
print(f"\n📋 LONGITUDINAL DATA EXTRACTION SUMMARY")
print("=" * 60)
print(f"Participant counts:")
print(f" Brain longitudinal (smri_diff): {len(smri_diff)}")
print(f" Behavior longitudinal (mh_diff): {len(mh_diff)}")
print(f" Complete longitudinal: {len(smri_diff_complete)}")
print(f"\nDataset shapes:")
print(f" smri_diff: {smri_diff.shape}")
print(f" mh_diff: {mh_diff.shape}")
print(f" smri_diff_complete: {smri_diff_complete.shape}")
print(f" mh_diff_complete: {mh_diff_complete.shape}")
print(f" gen_longitudinal: {gen_longitudinal.shape}")
print(f"\nData verification:")
# Check that all participants in complete datasets are the same
smri_diff_ids = set(smri_diff_complete['participant_id'])
mh_complete_ids = set(mh_diff_complete['participant_id'])
gen_complete_ids = set(gen_longitudinal['participant_id'])
if smri_diff_ids == mh_complete_ids == gen_complete_ids:
print(f" ✅ All complete datasets have matching participants")
else:
print(f" ⚠️ Participant mismatch in complete datasets")
print(f" SMRI: {len(smri_diff_ids)}, MH: {len(mh_complete_ids)}, GEN: {len(gen_complete_ids)}")
# Show sample of changes
print(f"\nSample brain changes (first participant, first 5 measures):")
brain_measures = [col for col in smri_diff.columns if col != 'participant_id'][:5]
for measure in brain_measures:
sample_change = smri_diff[measure].iloc[0]
print(f" {measure}: {sample_change:.4f}")
print(f"\nSample behavior changes (first participant, first 5 behaviors):")
behavior_measures = [col for col in mh_diff.columns if col != 'participant_id'][:5]
for behavior in behavior_measures:
sample_change = mh_diff[behavior].iloc[0]
print(f" {behavior}: {sample_change:.4f}")
print(f"\n✨ LONGITUDINAL DIFFERENCE DATASETS CREATED!")
print(f"📊 Available for analysis:")
print(f" - smri_diff: Brain changes (2yr - baseline) for {len(smri_diff)} participants")
print(f" - mh_diff: Behavior changes (2yr - baseline) for {len(mh_diff)} participants")
print(f" - smri_diff_complete: Brain changes for complete longitudinal sample")
print(f" - mh_diff_complete: Behavior changes for complete longitudinal sample")
print(f" - gen_longitudinal: Demographics for complete longitudinal sample")
# %%
# ==========================================
# LONGITUDINAL DIFFERENCE PREDICTION PIPELINE
# ==========================================
def run_longitudinal_diff_predictions(smri_diff, mh_diff, gen_longitudinal,
results_dir, ml_config, n_null=100):
"""
Train NEW models optimized for longitudinal difference data using existing functions
"""
print("🧠 TRAINING NEW MODELS FOR LONGITUDINAL DIFFERENCE DATA")
print("="*70)
# Create separate results directory for longitudinal analysis
longitudinal_results_dir = os.path.join(results_dir, 'longitudinal_diff_results_2y')
os.makedirs(longitudinal_results_dir, exist_ok=True)
# Check if results already exist
longitudinal_final_filename = f"{longitudinal_results_dir}/longitudinal_diff_results_complete.pkl"
if os.path.exists(longitudinal_final_filename):
print(f"Loading existing longitudinal results from {longitudinal_final_filename}")
with open(longitudinal_final_filename, 'rb') as f:
longitudinal_results = pickle.load(f)
print(f"Loaded {len(longitudinal_results['models'])} existing longitudinal models")
print(f"Loaded {len(longitudinal_results['null_models'])} existing longitudinal null models")
else:
print("Creating new longitudinal results dictionary...")
# Initialize SEPARATE results dictionary for longitudinal analysis
longitudinal_results = {
'models': {}, # NEW trained models for differences
'null_models': {}, # Null model results
'metadata': { # Analysis info
'analysis_type': 'longitudinal_differences_new_models',
'ml_config': ml_config,
'mh_behaviors': [col for col in mh_diff.columns if col != 'participant_id'],
'feature_columns': {}, # Will be populated
'timepoint': 'longitudinal_difference_2yr_minus_baseline',
'brain_types': ['thickness', 'area', 'volume', 'subcvolume'], # Will be populated
'sex_groups': ['all', 'female', 'male'],
'creation_date': pd.Timestamp.now().isoformat()
},
'summary': {}, # Real model summaries
'null_summary': {} # Null model summaries
}
# Get behavior columns from mh_diff
mh_diff_cols = [col for col in mh_diff.columns if col != 'participant_id']
# Define brain types based on what's available in smri_diff
print("🔍 Identifying brain measures in difference data...")
# Get all brain columns from smri_diff (excluding participant_id)
all_brain_cols = [col for col in smri_diff.columns if col != 'participant_id']
# Categorize columns by brain type (match your original column patterns)
longitudinal_brain_types = {}
# Thickness columns - look for thickness patterns in difference data
thickness_diff_cols = [col for col in all_brain_cols if 'thick' in col.lower()]
if len(thickness_diff_cols) > 0:
longitudinal_brain_types['thickness'] = thickness_diff_cols
# Area columns - look for area patterns in difference data
area_diff_cols = [col for col in all_brain_cols if 'area' in col.lower()]
if len(area_diff_cols) > 0:
longitudinal_brain_types['area'] = area_diff_cols
# Volume columns - look for volume patterns in difference data
volume_diff_cols = [col for col in all_brain_cols if 'vol' in col.lower() and 'aseg' not in col.lower()]
if len(volume_diff_cols) > 0:
longitudinal_brain_types['volume'] = volume_diff_cols
# Subcortical volume columns - look for subcortical patterns
subcvolume_diff_cols = [col for col in all_brain_cols if 'vol' in col.lower() and 'aseg' in col.lower()]
if len(subcvolume_diff_cols) > 0:
longitudinal_brain_types['subcvolume'] = subcvolume_diff_cols
print(f"Brain types found in difference data:")
for brain_type, cols in longitudinal_brain_types.items():
print(f" {brain_type}: {len(cols)} features")
# Update metadata
longitudinal_results['metadata']['feature_columns'] = longitudinal_brain_types
longitudinal_results['metadata']['brain_types'] = list(longitudinal_brain_types.keys())
longitudinal_results['metadata']['mh_behaviors'] = mh_diff_cols
print(f"\nLongitudinal participants:")
print(f" Brain differences: {len(smri_diff)}")
print(f" Behavior differences: {len(mh_diff)}")
print(f" Behavior measures: {len(mh_diff_cols)}")
print(f"Results will be saved to: {longitudinal_results_dir}")
print("="*70)
# Merge datasets for analysis
print("\n📊 MERGING LONGITUDINAL DATA...")
# Merge smri_diff with mh_diff
merged_diff = smri_diff.merge(mh_diff, on='participant_id', how='inner')
print(f"After merging brain and behavior differences: {len(merged_diff)} participants")
# Merge with demographics (need sex and site info)
# Use baseline demographics since we need sex/site info
gen_baseline = gen_longitudinal[gen_longitudinal['session_id'] == 'ses-00A'].copy()
merged_diff = merged_diff.merge(
gen_baseline[['participant_id', 'ab_g_dyn__design_site', 'ab_g_stc__cohort_sex']],
on='participant_id',
how='inner'
)
print(f"After merging with demographics: {len(merged_diff)} participants")
if len(merged_diff) == 0:
print("❌ No merged longitudinal data available!")
return None
# Update participant count
longitudinal_results['metadata']['n_participants'] = len(merged_diff)
# Prepare data
Y_diff = merged_diff[mh_diff_cols].values # Behavior differences
site = merged_diff['ab_g_dyn__design_site'].values
print(f"Final dataset shape: {merged_diff.shape}")
print(f"Y (behavior differences) shape: {Y_diff.shape}")
print(f"Number of sites: {len(np.unique(site))}")
# Sex groups to analyze
sexes = ['all', 'female', 'male']
# Loop through sexes
for sex in sexes:
print(f"\n{'='*20} SEX GROUP: {sex.upper()} {'='*20}")
# Filter by sex
if sex == 'all':
data_subset = merged_diff.copy()
Y_subset = Y_diff.copy()
site_subset = site.copy()
elif sex == 'male':
sex_mask = merged_diff['ab_g_stc__cohort_sex'] == 1
data_subset = merged_diff[sex_mask].copy()
Y_subset = Y_diff[sex_mask]
site_subset = site[sex_mask]
elif sex == 'female':
sex_mask = merged_diff['ab_g_stc__cohort_sex'] == 2
data_subset = merged_diff[sex_mask].copy()
Y_subset = Y_diff[sex_mask]
site_subset = site[sex_mask]
if len(data_subset) < 50: # Minimum sample size
print(f"Skipping {sex} group (n={len(data_subset)}, too small)")
continue
print(f"{sex} group: {len(data_subset)} participants")
print(f"Number of sites: {len(np.unique(site_subset))}")
# Loop through brain types
for brain_type, feature_cols in longitudinal_brain_types.items():
print(f"\n--- BRAIN TYPE: {brain_type.upper()} ---")
# Check if this model already exists
result_key = f"longitudinal_diff_{brain_type}_{sex}"
if result_key in longitudinal_results['models']:
print(f" {result_key} already exists, skipping...")
continue
# Prepare brain difference data
X_brain_diff = data_subset[feature_cols].values
print(f"X (brain differences) shape: {X_brain_diff.shape}")
print(f"Feature columns: {len(feature_cols)}")
# Check for valid data
if X_brain_diff.shape[1] == 0:
print(f"No valid brain features for {brain_type}, skipping...")
continue
print(f"\nTraining NEW {brain_type} model for {sex} (longitudinal differences)...")
# ==========================================
# TRAIN NEW REAL MODEL using your existing function
# ==========================================
print(" Training NEW REAL model on difference data...")
results = train_test_ridge_model(
X_brain_diff, Y_subset, site_subset, mh_diff_cols,
f"{brain_type}_diff", sex, "longitudinal_diff",
longitudinal_results, longitudinal_results_dir,
config=ml_config
)
# ==========================================
# TRAIN NEW NULL MODEL using your existing function
# ==========================================
print(" Training NEW NULL model on difference data...")
null_results = generate_null_models(
X_brain_diff, Y_subset, site_subset, mh_diff_cols,
f"{brain_type}_diff", sex, "longitudinal_diff", results,
longitudinal_results, longitudinal_results_dir,
config=ml_config, n_null=n_null
)
print(f" Completed: {result_key}")
print(f" Real mean correlation: {np.mean(results['corr']):.3f}")
print(f" Null mean correlation: {np.mean(null_results['corr']):.3f}")
print(f" Improvement over null: {np.mean(results['corr']) - np.mean(null_results['corr']):.3f}")
print("-" * 50)
# Save NEW longitudinal results
print(f"\n{'='*50}")
print("SAVING NEW LONGITUDINAL DIFFERENCE RESULTS...")
# Update metadata
longitudinal_results['metadata']['total_models'] = len(longitudinal_results['models'])
longitudinal_results['metadata']['total_null_models'] = len(longitudinal_results['null_models'])
longitudinal_results['metadata']['last_updated'] = pd.Timestamp.now().isoformat()
# Save complete longitudinal results
with open(longitudinal_final_filename, 'wb') as f:
pickle.dump(longitudinal_results, f)
# Save summary only
longitudinal_summary_filename = f"{longitudinal_results_dir}/longitudinal_diff_summary.pkl"
with open(longitudinal_summary_filename, 'wb') as f:
pickle.dump(longitudinal_results['summary'], f)
# Save null summary
longitudinal_null_summary_filename = f"{longitudinal_results_dir}/longitudinal_diff_null_summary.pkl"
with open(longitudinal_null_summary_filename, 'wb') as f:
pickle.dump(longitudinal_results['null_summary'], f)
# Save metadata
longitudinal_metadata_filename = f"{longitudinal_results_dir}/longitudinal_diff_metadata.pkl"
with open(longitudinal_metadata_filename, 'wb') as f:
pickle.dump(longitudinal_results['metadata'], f)
print("\nNEW Longitudinal Difference ML Pipeline completed!")
print("🎉 NEW LONGITUDINAL MODELS TRAINED AND SAVED!")
print(f"📁 Directory: {longitudinal_results_dir}/")
print(f"📊 Total NEW longitudinal models: {len(longitudinal_results['models'])}")
print(f"📊 Total NEW longitudinal null models: {len(longitudinal_results['null_models'])}")
print("📋 Files saved:")
print(" - longitudinal_diff_results_complete.pkl")
print(" - longitudinal_diff_summary.pkl")
print(" - longitudinal_diff_null_summary.pkl")
print(" - longitudinal_diff_metadata.pkl")
print(f"\n📊 NEW longitudinal models trained:")
for key in sorted(longitudinal_results['models'].keys()):
if key in longitudinal_results['summary']:
mean_corr = longitudinal_results['summary'][key]['mean_correlation']
print(f" ✓ {key}: {mean_corr:.3f}")
return longitudinal_results
# ==========================================
# RUN NEW LONGITUDINAL DIFFERENCE ANALYSIS
# ==========================================
print("🚀 TRAINING NEW MODELS FOR LONGITUDINAL DIFFERENCE DATA...")
# Train NEW models optimized for longitudinal differences
longitudinal_results = run_longitudinal_diff_predictions(
smri_diff=smri_diff,
mh_diff=mh_diff,
gen_longitudinal=gen_longitudinal,
results_dir=results_dir,
ml_config=ML_CONFIG,
n_null=100 # Adjust as needed
)
if longitudinal_results:
print(f"\n✨ LONGITUDINAL DIFFERENCE ANALYSIS COMPLETED!")
print(f"📊 NEW models trained specifically for longitudinal changes")
print(f"📁 Original cross-sectional results: {results_dir}/all_results_complete.pkl")
print(f"📁 NEW longitudinal results: {results_dir}/longitudinal_diff_results/longitudinal_diff_results_complete.pkl")
# Show summary
print(f"\n📊 Longitudinal model performance:")
for key, summary in longitudinal_results['summary'].items():
print(f" {key}: {summary['mean_correlation']:.3f} ± {summary['std_correlation']:.3f}")
else:
print("❌ Longitudinal analysis failed")
# %%
# ==========================================
# SIGNIFICANCE TESTING FOR LONGITUDINAL DIFFERENCE MODELS
# ==========================================
def test_longitudinal_significance(longitudinal_results, results_dir):
"""
Test significance for longitudinal difference models with FDR correction
"""
print(f"\n{'='*60}")
print("TESTING LONGITUDINAL MODEL SIGNIFICANCE...")
print(f"{'='*60}")
# Load longitudinal results if not provided
if longitudinal_results is None:
longitudinal_results_dir = os.path.join(results_dir, 'longitudinal_diff_results_2y')
longitudinal_final_filename = f"{longitudinal_results_dir}/longitudinal_diff_results_complete.pkl"
if os.path.exists(longitudinal_final_filename):
with open(longitudinal_final_filename, 'rb') as f:
longitudinal_results = pickle.load(f)
else:
print("❌ No longitudinal results found!")
return None
model_keys = list(longitudinal_results['models'].keys())
print(f"Longitudinal models found: {len(model_keys)}")
for key in model_keys:
print(f" - {key}")
# Group model keys by sex group (since all are same timepoint: longitudinal_diff)
sex_groups = {}
for model_key in model_keys:
# Extract sex group from model key (e.g., 'longitudinal_diff_thickness_all' -> 'all')
parts = model_key.split('_')
sex_group = parts[-1] # Last part is sex group
if sex_group not in sex_groups:
sex_groups[sex_group] = []
sex_groups[sex_group].append(model_key)
print(f"\nSex groups found: {list(sex_groups.keys())}")
for sex_group, keys in sex_groups.items():
print(f" {sex_group}: {len(keys)} models")
# Initialize dictionaries to store significance results
longitudinal_significance_results = {}
longitudinal_fdr_results = {}
longitudinal_group_summaries = {}
# Process each sex group separately
for sex_group, group_model_keys in sex_groups.items():
print(f"\n{'='*40}")
print(f"PROCESSING SEX GROUP: {sex_group.upper()}")
print(f"{'='*40}")
# Extract data for this sex group only
group_real_corrs = []
group_null_corrs = []
group_model_info = []
group_significance_results = {}
print(f"Extracting correlation data for {sex_group} group...")
for model_key in group_model_keys:
real_model = longitudinal_results['models'][model_key]
null_model = longitudinal_results['null_models'][model_key]
n_beh = real_model['config']['n_behaviors']
# Get mean correlations across repetitions for each behavior
real_corrs = np.mean(real_model['corr'], axis=1) # Shape: (n_behaviors,)
null_corrs = null_model['corr'] # Shape: (n_behaviors, n_reps)
# Store for FDR correction within this sex group
group_real_corrs.extend(real_corrs)
group_null_corrs.extend(null_corrs.T) # Transpose to get (n_reps, n_behaviors), then flatten
# Store model info for each behavior
for b in range(n_beh):
group_model_info.append({
'model_key': model_key,
'behavior_idx': b,
'behavior_name': real_model['config']['behavior_names'][b]
})
# Calculate p-values for this model
p_vals = np.zeros(n_beh)
for b in range(n_beh):
real_mean = real_corrs[b]
null_dist = null_corrs[b, :]
p_vals[b] = get_null_p(real_mean, null_dist)
# Store results for this model
group_significance_results[model_key] = {
'p_values': p_vals,
'real_correlations': real_corrs,
'null_correlations': null_corrs,
'behavior_names': real_model['config']['behavior_names']
}
print(f"Collected data for {len(group_real_corrs)} tests in {sex_group} group")
# ==========================================
# FDR CORRECTION WITHIN THIS SEX GROUP
# ==========================================
print(f"Applying FDR correction within {sex_group} group...")
# Collect all p-values for this sex group
group_all_p_vals = []
for model_key in group_model_keys:
group_all_p_vals.extend(group_significance_results[model_key]['p_values'])
group_all_p_vals = np.array(group_all_p_vals)
# Apply FDR correction within this sex group
group_fdr_rejected, group_fdr_p_corrected = fdr(group_all_p_vals, alpha=0.05, method='indep')
print(f" Tests in {sex_group}: {len(group_all_p_vals)}")
print(f" Significant after FDR: {np.sum(group_fdr_rejected)}")
print(f" Significance rate: {np.mean(group_fdr_rejected):.3f}")
# ==========================================
# ORGANIZE FDR RESULTS BY MODEL WITHIN SEX GROUP
# ==========================================
print(f"Organizing FDR results for {sex_group}...")
idx = 0
group_fdr_results = {}
for model_key in group_model_keys:
n_beh = len(group_significance_results[model_key]['p_values'])
# Extract FDR results for this model
model_fdr_rejected = group_fdr_rejected[idx:idx+n_beh]
model_fdr_p_corrected = group_fdr_p_corrected[idx:idx+n_beh]
# Store FDR results
group_fdr_results[model_key] = {
'fdr_rejected': model_fdr_rejected,
'fdr_p_corrected': model_fdr_p_corrected,
'n_significant': np.sum(model_fdr_rejected),
'real_correlations': group_significance_results[model_key]['real_correlations'],
'p_values_uncorrected': group_significance_results[model_key]['p_values'],
'behavior_names': group_significance_results[model_key]['behavior_names']
}
idx += n_beh
# Store results for this sex group
longitudinal_significance_results.update(group_significance_results)
longitudinal_fdr_results.update(group_fdr_results)
# Create group summary
group_total_significant = sum([group_fdr_results[key]['n_significant'] for key in group_model_keys])
longitudinal_group_summaries[sex_group] = {
'sex_group': sex_group,
'total_tests': len(group_all_p_vals),
'total_significant': int(group_total_significant),
'significance_rate': float(group_total_significant/len(group_all_p_vals)) if len(group_all_p_vals) > 0 else 0.0,
'models': group_model_keys,
'fdr_alpha': 0.05,
'fdr_method': 'indep'
}
# Print summary for this group
print(f"\n{sex_group} GROUP SUMMARY:")
print(f" Total tests: {len(group_all_p_vals)}")
print(f" Significant: {group_total_significant}")
print(f" Rate: {group_total_significant/len(group_all_p_vals):.3f}" if len(group_all_p_vals) > 0 else " Rate: 0.000")
# Show significant models for this group
print(f" Model results:")
for model_key in group_model_keys:
n_sig = group_fdr_results[model_key]['n_significant']
n_total = len(group_fdr_results[model_key]['fdr_rejected'])
brain_type = model_key.split('_')[2] # Extract brain type from model key
print(f" {brain_type}: {n_sig}/{n_total} significant")
# ==========================================
# OVERALL SUMMARY FOR LONGITUDINAL MODELS
# ==========================================
print(f"\n{'='*60}")
print("LONGITUDINAL MODEL SIGNIFICANCE SUMMARY:")
print("="*60)
total_tests_all = 0
total_significant_all = 0
for sex_group in sorted(longitudinal_group_summaries.keys()):
summary = longitudinal_group_summaries[sex_group]
print(f"{sex_group:>8}: {summary['total_significant']:>3}/{summary['total_tests']:<3} ({summary['significance_rate']:.3f})")
total_tests_all += summary['total_tests']
total_significant_all += summary['total_significant']
print(f"\nGRAND TOTAL:")
print(f" Total tests: {total_tests_all}")
print(f" Total significant: {total_significant_all}")
print(f" Overall rate: {total_significant_all/total_tests_all:.3f}")
# ==========================================
# DETAILED RESULTS BY BRAIN TYPE AND BEHAVIOR
# ==========================================
print(f"\n{'='*60}")
print("DETAILED RESULTS BY BRAIN TYPE:")
print("="*60)
# Organize results by brain type
brain_type_results = {}
for model_key, fdr_result in longitudinal_fdr_results.items():
brain_type = model_key.split('_')[2] # Extract brain type
sex_group = model_key.split('_')[3] # Extract sex group
if brain_type not in brain_type_results:
brain_type_results[brain_type] = {}
brain_type_results[brain_type][sex_group] = fdr_result
# Print detailed results
for brain_type in sorted(brain_type_results.keys()):
print(f"\n{brain_type.upper()}:")
for sex_group in ['all', 'female', 'male']:
if sex_group in brain_type_results[brain_type]:
result = brain_type_results[brain_type][sex_group]
print(f" {sex_group:>6}: {result['n_significant']}/{len(result['fdr_rejected'])} significant")
# Show top significant behaviors
if result['n_significant'] > 0:
sig_indices = np.where(result['fdr_rejected'])[0]
sig_corrs = result['real_correlations'][sig_indices]
sig_behaviors = [result['behavior_names'][i] for i in sig_indices]
# Sort by correlation strength
sorted_indices = np.argsort(np.abs(sig_corrs))[::-1]
print(f" Top significant behaviors:")
for i in sorted_indices[:3]: # Show top 3
behavior = sig_behaviors[i]
corr = sig_corrs[i]
print(f" {behavior}: r={corr:.3f}")
# ==========================================
# ADD TO LONGITUDINAL RESULTS AND SAVE
# ==========================================
print(f"\nAdding significance results to longitudinal results...")
# Add significance results to longitudinal_results
longitudinal_results['significance'] = longitudinal_significance_results
longitudinal_results['fdr_results'] = longitudinal_fdr_results
longitudinal_results['significance_summary'] = {
'by_sex_group': longitudinal_group_summaries,
'overall': {
'total_tests': total_tests_all,
'total_significant': total_significant_all,
'significance_rate': float(total_significant_all/total_tests_all) if total_tests_all > 0 else 0.0,
'fdr_alpha': 0.05,
'fdr_method': 'indep',
'correction_level': 'within_sex_group'
}
}
# Save updated longitudinal results
longitudinal_results_dir = os.path.join(results_dir, 'longitudinal_diff_results_2y')
longitudinal_final_filename = f"{longitudinal_results_dir}/longitudinal_diff_results_complete_with_significance.pkl"
with open(longitudinal_final_filename, 'wb') as f:
pickle.dump(longitudinal_results, f)
# Save just significance results
longitudinal_significance_filename = f"{longitudinal_results_dir}/longitudinal_significance_results_only.pkl"
with open(longitudinal_significance_filename, 'wb') as f:
pickle.dump({
'significance': longitudinal_significance_results,
'fdr_results': longitudinal_fdr_results,
'summary': longitudinal_results['significance_summary']
}, f)
print(f"\n🎉 LONGITUDINAL SIGNIFICANCE TESTING COMPLETED!")
print(f"📁 Results saved to:")
print(f" - longitudinal_diff_results_complete_with_significance.pkl")
print(f" - longitudinal_significance_results_only.pkl")
return longitudinal_results
# ==========================================
# RUN LONGITUDINAL SIGNIFICANCE TESTING
# ==========================================
print("🧠 TESTING SIGNIFICANCE FOR LONGITUDINAL MODELS...")
# Test significance for longitudinal models
longitudinal_results_with_significance = test_longitudinal_significance(
longitudinal_results,
results_dir
)
# ==========================================
# COMPARE CROSS-SECTIONAL VS LONGITUDINAL SIGNIFICANCE
# ==========================================
def compare_significance_rates(results_dir):
"""Compare significance rates between cross-sectional and longitudinal models"""
print(f"\n{'='*60}")
print("COMPARING CROSS-SECTIONAL VS LONGITUDINAL SIGNIFICANCE:")
print("="*60)
# Load cross-sectional results
cross_sectional_file = f"{results_dir}/all_results_complete_with_significance_by_timepoint_sex.pkl"
if os.path.exists(cross_sectional_file):
with open(cross_sectional_file, 'rb') as f:
cross_sectional_results = pickle.load(f)
cs_overall = cross_sectional_results['significance_summary']['overall']
print(f"Cross-sectional significance rate: {cs_overall['significance_rate']:.3f}")
print(f" ({cs_overall['total_significant']}/{cs_overall['total_tests']} tests)")
else:
print("Cross-sectional significance results not found")
# Load longitudinal results
longitudinal_dir = os.path.join(results_dir, 'longitudinal_diff_results_2y')
longitudinal_file = f"{longitudinal_dir}/longitudinal_diff_results_complete_with_significance.pkl"
if os.path.exists(longitudinal_file):
with open(longitudinal_file, 'rb') as f:
longitudinal_results = pickle.load(f)
long_overall = longitudinal_results['significance_summary']['overall']
print(f"Longitudinal significance rate: {long_overall['significance_rate']:.3f}")
print(f" ({long_overall['total_significant']}/{long_overall['total_tests']} tests)")
else:
print("Longitudinal significance results not found")
# Run comparison
compare_significance_rates(results_dir)
# %%

nbdc_abcd_bpm_smri_impulsivity.ipynb at commit 85451d8, no license · at the source

Overview

Authors: Elvisha Dhamala^1,2, Erynn Christensen¹, Jamie L Hanson^3,4, Jocelyn A Ricard^5,6, Noelle Arcaro¹, Simran Bhola¹, Lisa Wiersch¹, Katharina Brosch¹, B T Thomas Yeo^{7,8,9,10,11,12}, Avram J Holmes¹³, Sarah W Yip^14,15,16,17

ORCID iDs: Elvisha Dhamala, Jocelyn A Ricard, Noelle Arcaro, Katharina Brosch, Avram J Holmes, Sarah W Yip

17 affiliations

Institute of Behavioral Sciences, Feinstein Institutes for Medical Research, Manhasset, USA
Department of Psychiatry, Donald and Barbara Zucker School of Medicine at Hofstra/Northwell, Uniondale, Hempstead, USA
Learning Research and Development Center, University of Pittsburgh, Pittsburgh, USA
Department of Psychology, University of Pittsburgh, Pittsburgh, USA
Stanford Neurosciences Interdepartmental Program, Stanford University School of Medicine, Stanford, USA
Department of Psychology, Stanford University, Stanford, USA
Centre for Sleep & Cognition & Centre for Translational Magnetic Resonance Research, Yong Loo Lin School of Medicine, Singapore, National University of Singapore, Singapore, Singapore
Department of Electrical and Computer Engineering, National University of Singapore, Singapore, Singapore
N.1 Institute for Health, National University of Singapore, Singapore, Singapore
Integrative Sciences and Engineering Programme (ISEP), National University of Singapore, Singapore, Singapore
Department of Medicine, Healthy Longevity Research Programme, Human Potential Translational Research Programme & Institute for Digital Medicine (WisDM), Yong Loo Lin School of Medicine, National University of Singapore, Singapore, Singapore
Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Charlestown, USA
Department of Psychiatry, Brain Health Institute, Rutgers University, Piscataway, USA
Department of Psychiatry, Yale School of Medicine, New Haven, USA
Child Study Center, Yale School of Medicine, New Haven, USA
Center for Brain and Mind Health, Yale School of Medicine, New Haven, USA
Wu Tsai Institute, Yale University, New Haven, USA

Institutions: Feinstein Institute for Medical Research (United States); Hofstra University (United States); Donald & Barbara Zucker School of Medicine at Hofstra/Northwell (United States); University of Pittsburgh (United States); Stanford Medicine (United States); Stanford University (United States); National University of Singapore (Singapore); Massachusetts General Hospital (United States); Athinoula A. Martinos Center for Biomedical Imaging (United States); Rutgers, The State University of New Jersey (United States); Yale University (United States)

Journal: Molecular psychiatry, volume 31, issue 8, pages 4282-4296

Dates: received 16 June 2025; accepted 10 March 2026; published online 24 March 2026; in print 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1038/s41380-026-03526-2 · PMID 41876706 · PMCID PMC13364653 · OpenAlex W7140211364

Open access: hybrid, a free copy (OpenAlex)

Status: code verified (how the code is found and verified)

Categories: structural MRI / diffusion (modality), human (organism)

Methods: Statistics, Machine learning, Preprocessing

Keywords: Neuroscience, Biological techniques, Psychiatric disorders

MeSH: Impulsive Behavior*, Adolescent, Brain, Brain Mapping, Decision Making, Female, Gray Matter, Humans, Individuality, Machine Learning, Magnetic Resonance Imaging, Male, Neuroanatomy (* major topic)

Topic: Neural and Behavioral Psychology Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: Yale University; Northwell Health; National Research Foundation; National Research Foundation Singapore; Feinstein Institutes for Medical Research; Canadian Institutes of Health Research; Medical Research Council; National Medical Research Council (STaR20nov-0003, OFLCG19MAY-0035, OFIRG24jul-0049, CTGIIT23jan-0001, OFIRG24jan-0006, CG21APR1009); National Institute on Alcohol Abuse and Alcoholism (R01AA027553); National Institute of Mental Health; National Institute on Drug Abuse (R01 DA053301)

Citations: not cited yet (Europe PMC); 125 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 2 matches between paragraphs and lines of code.

elvisha/neuroanat_impulsivity

License: none: the authors keep all their rights

State: the link answers, verified on 30 September 2026

Evidence: files inventoried

Commit: 85451d816fb0de87cbbaa4bbcc04285f6fc54225, 14 October 2025

Languages: Jupyter (1)

Size: 2 files, 1 script

Software Heritage: not archived

Found in: “Code availability”

Holds: README, 1 notebook

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file), scikit-learn (1 file), SciPy (1 file), seaborn (1 file), statsmodels (1 file)

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

30 September 2026: the link answers

2 files

nbdc_abcd_bpm_smri_impulsivity.ipynb, Jupyter, 2,873 lines, 2 matches
README.md, Text, 3 lines

At the source: github.com/elvisha/neuroanat_impulsivity

Code availability statement

The paper has a code 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: elvisha/neuroanat_impulsivity

Read it in the paper: doi.org/10.1038/s41380-026-03526-2.

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;
1 script, each with its path and the digest of its content;
2 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.

Data availability statement

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

no repository, dataset or request procedure was recognized in it

Read it in the paper: doi.org/10.1038/s41380-026-03526-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, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 11 authors, 3 keywords, 13 MeSH terms, 11 funders, 112 references.

Cite

This paper

Dhamala, E., Christensen, E., Hanson, J. L., Ricard, J. A., Arcaro, N., Bhola, S., Wiersch, L., Brosch, K., Yeo, B. T. T., Holmes, A. J., & Yip, S. W. (2026). Neuroanatomy reflects individual variability in impulsivity in youth. Molecular psychiatry, 31(8), 4282-4296. https://doi.org/10.1038/s41380-026-03526-2

BibTeX

@article{dhamala2026neuroanatomy, author = {Dhamala, Elvisha and Christensen, Erynn and Hanson, Jamie L and Ricard, Jocelyn A and Arcaro, Noelle and Bhola, Simran and Wiersch, Lisa and Brosch, Katharina and Yeo, B T Thomas and Holmes, Avram J and Yip, Sarah W}, title = {{Neuroanatomy reflects individual variability in impulsivity in youth}}, journal = {Molecular psychiatry}, year = {2026}, month = mar, volume = {31}, number = {8}, pages = {4282--4296}, publisher = {Springer Nature}, issn = {1359-4184}, doi = {10.1038/s41380-026-03526-2}, url = {https://doi.org/10.1038/s41380-026-03526-2}, pmid = {41876706}, pmcid = {PMC13364653} }

RIS

TY - JOUR AU - Dhamala, Elvisha AU - Christensen, Erynn AU - Hanson, Jamie L AU - Ricard, Jocelyn A AU - Arcaro, Noelle AU - Bhola, Simran AU - Wiersch, Lisa AU - Brosch, Katharina AU - Yeo, B T Thomas AU - Holmes, Avram J AU - Yip, Sarah W TI - Neuroanatomy reflects individual variability in impulsivity in youth T2 - Molecular psychiatry J2 - Mol Psychiatry PY - 2026 DA - 2026/03/24 VL - 31 IS - 8 SP - 4282 EP - 4296 SN - 1359-4184 PB - Springer Nature DO - 10.1038/s41380-026-03526-2 UR - https://doi.org/10.1038/s41380-026-03526-2 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41380-026-03526-2", "type": "article-journal", "title": "Neuroanatomy reflects individual variability in impulsivity in youth", "container-title": "Molecular psychiatry", "author": [ { "family": "Dhamala", "given": "Elvisha" }, { "family": "Christensen", "given": "Erynn" }, { "family": "Hanson", "given": "Jamie L" }, { "family": "Ricard", "given": "Jocelyn A" }, { "family": "Arcaro", "given": "Noelle" }, { "family": "Bhola", "given": "Simran" }, { "family": "Wiersch", "given": "Lisa" }, { "family": "Brosch", "given": "Katharina" }, { "family": "Yeo", "given": "B T Thomas" }, { "family": "Holmes", "given": "Avram J" }, { "family": "Yip", "given": "Sarah W" } ], "container-title-short": "Mol Psychiatry", "volume": "31", "issue": "8", "page": "4282-4296", "DOI": "10.1038/s41380-026-03526-2", "PMID": "41876706", "PMCID": "PMC13364653", "ISSN": "1359-4184", "publisher": "Springer Nature", "URL": "https://doi.org/10.1038/s41380-026-03526-2", "language": "en", "issued": { "date-parts": [ [ 2026, 3, 24 ] ] } }

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.1126/sciadv.adn4202 [code]: Functional brain networks are associated with both sex and gender in children.
Journal: Science advances
In common: statsmodels, seaborn, scikit-learn, 4 other tools, 18 references, 2 authors
[2] doi:10.1126/sciadv.adn1862 [code]: Generalizable and replicable brain-based predictions of cognitive functioning across common psychiatric illness.
Journal: Science advances
In common: scikit-learn, pandas, SciPy, 2 other tools, 12 references, 3 authors
[3] doi:10.1038/s41467-026-73668-y [code]: Convergent and divergent brain-cognition development in early adolescence.
Journal: Nature communications
In common: seaborn, scikit-learn, pandas, 3 other tools, 13 references, author Elvisha Dhamala
[4] doi:10.1038/s42003-024-06836-9 [code]: A shared spatial topography links the functional connectome correlates of cocaine use disorder and dopamine D<sub>2/3</sub> receptor densities.
Journal: Communications biology
In common: seaborn, pandas, SciPy, 2 other tools, 1 reference, 3 authors
[5] doi:10.1038/s41514-026-00456-9 [code]: Exploring the link between body physiology and cognition: the role of the brain and aging.
Journal: npj aging
In common: statsmodels, seaborn, scikit-learn, 4 other tools, 10 references
[6] doi:10.1038/s44220-025-00523-2 [code]: Sex-specific differences in brain activity dynamics of youth with a family history of substance use disorder.
Journal: Nature. Mental health
In common: SciPy, Matplotlib, NumPy, 11 references
[7] doi:10.1038/s41380-025-03196-6 [code]: Impulsivity and neuroticism share distinct functional connectivity signatures with alcohol-use risk in youth.
Journal: Molecular psychiatry
In common: 6 references, author Sarah W. Yip
[8] doi:10.1016/j.dcn.2025.101624 [code]: Brain functional connectivity, but not neuroanatomy, captures the interrelationship between sex and gender in preadolescents.
Journal: Developmental cognitive neuroscience
In common: statsmodels, seaborn, scikit-learn, 4 other tools, structural MRI / diffusion, 7 references
[9] doi:10.1038/s41586-025-09250-1 [code]: Longer scans boost prediction and cut costs in brain-wide association studies.
Journal: Nature
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, 9 references
[10] doi:10.1162/imag_a_00233 [code]: Multilayer meta-matching: Translating phenotypic prediction models from multiple datasets to small data.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, 9 references

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.

Add the badge to the code's README

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown

<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1038/s41380-026-03526-2.svg)](https://openscicode.org/badge/doi/10.1038/s41380-026-03526-2)
<!-- oscr-badge:end -->

reStructuredText

.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1038/s41380-026-03526-2.svg
   :target: https://openscicode.org/badge/doi/10.1038/s41380-026-03526-2
   :alt: OSCR

.. oscr-badge:end

HTML

<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1038/s41380-026-03526-2"><img src="https://openscicode.org/badge/doi/10.1038/s41380-026-03526-2.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

elvisha/neuroanat_impulsivity: README.md on GitHub

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge 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: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

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