<jupyter_start><jupyter_text>Analyzing replicability of structure-based multivariate BWAS on the Human Connectome Project dataset Imports<jupyter_code>import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge
from sklearn.svm import SVR
from sklearn.model_selection import KFold, train_test_split, cross_val_predict, GridSearchCV
from sklearn.pipeline import Pipeline
from sklearn.decomposition import PCA
from joblib import Parallel, delayed
from mlxtend.evaluate import permutation_test
import fnmatch
sns.set(rc={"figure.figsize":(4, 2)})
sns.set_style("whitegrid")<jupyter_output><empty_output><jupyter_text>Load HCP dataWe load functional network matrices (netmats) from the HCP1200-release, as published on connectomeDB: https://db.humanconnectome.org/Due to licensoing issues, data is not supplied with the repository, but can be downloaded from the ConnectomeDB.See [hcp_data/readme.md](hcp_data/readme.md) for more details.<jupyter_code># HCP data can be obtained from the connectomeDB
# data is not part of this repository
behavior = pd.read_csv('hcp_data/hcp1200_behavioral_data.csv')
behavior = behavior.set_index('Subject', drop=True)

# convert age to numeric
age = []
for s in behavior['Age']:
    if s == '36+':
        age.append(36)
    else:
        split = s.split(sep='-')
        age.append(np.mean((float(split[0]), float(split[1]))))

behavior['age'] = age
behavior<jupyter_output><empty_output><jupyter_text>Function to prepare target variable<jupyter_code>def create_data(target='CogTotalComp_AgeAdj', feature_data="FS_*_Thck"):
    # it's a good practice to use pandas for merging, messing up subject order can be painful

    features = fnmatch.filter(behavior.columns.values, feature_data)

    df = behavior.dropna(subset = [target] + features)
    y = df[target].values
    X = df[features].values
    return X, y
features = fnmatch.filter(behavior.columns.values, "FS_*_Thck")
len(features)<jupyter_output><empty_output><jupyter_text>Function implementing a single bootstrap iterationWe define a workhorse function which:- randomly samples the discovery and the replication datasets,- creates cross-validated estimates of predictive performance within the discovery sample- finalizes the model by fitting it to the whole discovery sample (overfits the discovery but not the replication sample)- use it to predict the replication sample<jupyter_code>def bootstrap_workhorse(X, y, sample_size, model, random_state, shuffle_y=False):

    #create discovery and replication samples by random sampling from the whole dataset (without replacement)

    if shuffle_y:
        rng = np.random.default_rng(random_state)
        y = rng.permutation(y)

    X_discovery, X_replication, y_discovery, y_replication = train_test_split(X, y, train_size=sample_size, test_size=sample_size, shuffle=True, random_state=random_state)

    cv = KFold(10)
    # obtain cross-validated predictions in the discovery sample

    predicted_discovery_cv = np.zeros_like(y_discovery)
    cor_per_fold = np.zeros(cv.n_splits)
    i = 0
    for train, test in cv.split(X=X_discovery, y=y_discovery):
        model.fit(X=X_discovery[train], y=y_discovery[train])
        predicted_discovery_cv[test] = model.predict(X=X_discovery[test])
        cor_per_fold[i] = np.corrcoef(y_discovery[test], predicted_discovery_cv[test])[0,1]
        i += 1
    # correlation between the cross-validated predictions and observations in the discovery sample
    # this is the correct, unbiased estimate!
    # calculated as mean test performance across all folds
    r_disc_cv = np.mean(cor_per_fold)
    # finalize model by training it on the full discovery sample (without cross-validation)
    final_model = model.fit(X=X_discovery, y=y_discovery)
    # obtain predictions with the final model on the discovery sample, note that this model actually overfits this sample.
    # we do this only to demonstrate biased estimates
    predicted_discovery_overfit = final_model.predict(X=X_discovery)
    # here we obtain the biased effect size (r) estimates for demonstrational purposes
    r_disc_overfit = np.corrcoef(predicted_discovery_overfit, y_discovery)[0, 1]

    # We use the final model to predict the replication sample
    # This is correct (no overfitting here), the final model did not see this data during training
    predicted_replication = final_model.predict(X=X_replication)
    # we obtain the out-of-sample prediction performance estimates
    r_rep = np.corrcoef(predicted_replication, y_replication)[0, 1]

    # below we calculate permutation-based p-values for all three effect size estimates (in-sample unbiased, in-sample biased, out-of-sample)
    # (one sided tests, testing for positive correlation)
    p_disc_cv = permutation_test(predicted_discovery_cv, y_discovery, method='approximate', num_rounds=1000, func=lambda x, y: np.corrcoef(x, y)[1][0],seed=random_state)
    p_disc_overfit = permutation_test(predicted_discovery_overfit, y_discovery, method='approximate', num_rounds=1000, func=lambda x, y: np.corrcoef(x, y)[1][0],seed=random_state)
    p_rep = permutation_test(predicted_replication, y_replication, method='approximate', num_rounds=1000, func=lambda x, y: np.corrcoef(x, y)[1][0],seed=random_state)
    # return results
    return r_disc_cv, r_disc_overfit, r_rep, p_disc_cv, p_disc_overfit, p_rep<jupyter_output><empty_output><jupyter_text>All set, now we start the analysis.  Replicability with sample sizes n=50, 100, 200, 300 and maxHere we train a few different models on 100 bootstrap samples.We aggregate the results of our workhorse function in `n_bootstrap`=100 bootstrap cases (run in parallel).The whole process is repeated for all sample sizes, fetaure_sets and target variables. Here we test age and 5 cognitive variables, including 'cognitive ability' (the 'main' target variable in the target paper)- age: age group of the participants- CogTotalComp_AgeAdj: total cognitive ability- PMAT24_A_CR, : Fluid Intelligence (Penn Progressive Matrices)- CardSort_AgeAdj: Executive Function/Cognitive Flexibility (Dimensional Change Card Sort)- Flanker_AgeAdj: Executive Function/Inhibition (Flanker Task)- PicSeq_AgeAdj: Episodic Memory (Picture Sequence Memory)  Reproducing the PCA+SVR-based model from the target paper Like in the target paper:- Both PCA and SVR are done inside the cross-validation- PCA reatains the firts k principal components that together explain 50% of the variance- scikit-learn makes sure that PCA is only fit for the training samples- both for the test sets (in the cross-validation) and the replication sample PCA is not re-fit, bt features are simply transformed with the already fit PCA<jupyter_code>%%time

random_state = 42
n_bootstrap = 100

features = {
    'FS_CT': "FS_*_Thck"
}

models = {
    'PCA_SVR': Pipeline([('pca', PCA(n_components=0.5)),
                         ('svr', SVR())])

}

# We aggregate all results here:
df = pd.DataFrame(columns=['connectivity','model','target','n','r_discovery_cv','r_discovery_overfit','r_replication','p_discovery_cv','p_discovery_overfit','p_replication'])

for feature_set in features:
    for model in models:
        for target_var in ['age', 'CogTotalComp_AgeAdj', 'PMAT24_A_CR', 'Flanker_AgeAdj', 'CardSort_AgeAdj', 'PicSeq_AgeAdj']:
            for sample_size in [50, 100, 200, 300, 'max']:

                print('*****************************************************************')
                print(feature_set, model, target_var, sample_size)

                X, y = create_data(target=target_var, feature_data=features[feature_set])

                if sample_size=='max':
                    sample_size = int(len(y)/2)

                # create random seeds for each bootstrap iteration for reproducibility
                rng = np.random.default_rng(random_state)
                random_sates = rng.integers(np.iinfo(np.int32).max, size=n_bootstrap)

                # run bootstrap iterations in parallel
                r_discovery_cv, r_discovery_overfit, r_replication, p_discovery_cv, p_discovery_overfit, p_replication = zip(
                    *Parallel(n_jobs=-1)(
                    delayed(bootstrap_workhorse)(X, y, sample_size, models[model], seed) for seed in random_sates))

                tmp_data_frame = pd.DataFrame({
                    'connectivity' : feature_set,
                    'model' : model,
                    'target' : target_var,
                    'n' : sample_size,
                    'r_discovery_cv': r_discovery_cv,
                    'r_discovery_overfit': r_discovery_overfit,
                    'r_replication': r_replication,
                    'p_discovery_cv': p_discovery_cv,
                    'p_discovery_overfit': p_discovery_overfit,
                    'p_replication': p_replication
                })
                #sns.scatterplot(x='r_replication', y='r_discovery_cv', data=tmp_data_frame)
                #plt.ylabel('in-sample (r)')
                #plt.xlabel('out-of-sample (r_pred)')
                #plt.show()
                print(tmp_data_frame.r_discovery_cv.mean(), tmp_data_frame.r_replication.mean())

                for alpha in [0.05, 0.01, 0.005, 0.001]:
                    print('Replicability at alpha =', alpha, ':',
                          (tmp_data_frame.loc[tmp_data_frame['p_discovery_cv']<alpha,'p_replication']<alpha).sum() / (tmp_data_frame['p_discovery_cv']<0.05).sum() * 100, '%')

                df = pd.concat((df, tmp_data_frame))
                df.reset_index(drop=True, inplace=True)
                df.to_csv('res/results_thck_PCA_SVR.csv')

df<jupyter_output>*****************************************************************
FS_CT PCA_SVR age 50
0.08589268853895356 0.11814430558345144
Replicability at alpha = 0.05 : 30.0 %
Replicability at alpha = 0.01 : 0.0 %
Replicability at alpha = 0.005 : 0.0 %
Replicability at alpha = 0.001 : 0.0 %
*****************************************************************
FS_CT PCA_SVR age 100
0.09976372318981333 0.1321490781638016
Replicability at alpha = 0.05 : 41.17647058823529 %
Replicability at alpha = 0.01 : 5.88235294117647 %
Replicability at alpha = 0.005 : 0.0 %
Replicability at alpha = 0.001 : 0.0 %
*****************************************************************
FS_CT PCA_SVR age 200
0.1319123546390201 0.12993985943621775
Replicability at alpha = 0.05 : 59.45945945945946 %
Replicability at alpha = 0.01 : 24.324324324324326 %
Replicability at alpha = 0.005 : 16.216216216216218 %
Replicability at alpha = 0.001 : 2.7027027027027026 %
*****************************************************************
FS_C[...]<jupyter_text>Now we fit a simple Ridge regression(no feature selection, no hyperparameter optimization)This can be expected to perform better on low samples than SVR.<jupyter_code>%%time

random_state = 42
n_bootstrap = 100

features = {
    'FS_CT': "FS_*_Thck"
}

models = {
    'ridge': Ridge()
}

# We aggregate all results here:
df = pd.DataFrame(columns=['connectivity','model','target','n','r_discovery_cv','r_discovery_overfit','r_replication','p_discovery_cv','p_discovery_overfit','p_replication'])

for feature_set in features:
    for model in models:
        for target_var in ['age', 'CogTotalComp_AgeAdj', 'PMAT24_A_CR', 'Flanker_AgeAdj', 'CardSort_AgeAdj', 'PicSeq_AgeAdj']:
            for sample_size in [50, 100, 200, 300, 'max']:

                print('*****************************************************************')
                print(feature_set, model, target_var, sample_size)

                X, y = create_data(target=target_var, feature_data=features[feature_set])

                if sample_size=='max':
                    sample_size = int(len(y)/2)

                # create random seeds for each bootstrap iteration for reproducibility
                rng = np.random.default_rng(random_state)
                random_sates = rng.integers(np.iinfo(np.int32).max, size=n_bootstrap)

                # run bootstrap iterations in parallel
                r_discovery_cv, r_discovery_overfit, r_replication, p_discovery_cv, p_discovery_overfit, p_replication = zip(
                    *Parallel(n_jobs=-1)(
                    delayed(bootstrap_workhorse)(X, y, sample_size, models[model], seed) for seed in random_sates))

                tmp_data_frame = pd.DataFrame({
                    'connectivity' : feature_set,
                    'model' : model,
                    'target' : target_var,
                    'n' : sample_size,
                    'r_discovery_cv': r_discovery_cv,
                    'r_discovery_overfit': r_discovery_overfit,
                    'r_replication': r_replication,
                    'p_discovery_cv': p_discovery_cv,
                    'p_discovery_overfit': p_discovery_overfit,
                    'p_replication': p_replication
                })
                print(tmp_data_frame.r_discovery_cv.mean(), tmp_data_frame.r_replication.mean())

                for alpha in [0.05, 0.01, 0.005, 0.001]:
                    print('Replicability at alpha =', alpha, ':',
                          (tmp_data_frame.loc[tmp_data_frame['p_discovery_cv']<alpha,'p_replication']<alpha).sum() / (tmp_data_frame['p_discovery_cv']<0.05).sum() * 100, '%')

                df = pd.concat((df, tmp_data_frame))
                df.reset_index(drop=True, inplace=True)
                df.to_csv('res/results_thck_Ridge.csv')

df<jupyter_output>*****************************************************************
FS_CT ridge age 50
0.15322562774488466 0.16690836840736034
Replicability at alpha = 0.05 : 55.00000000000001 %
Replicability at alpha = 0.01 : 10.0 %
Replicability at alpha = 0.005 : 5.0 %
Replicability at alpha = 0.001 : 0.0 %
*****************************************************************
FS_CT ridge age 100
0.15399594921982113 0.1668875020529473
Replicability at alpha = 0.05 : 58.97435897435898 %
Replicability at alpha = 0.01 : 12.82051282051282 %
Replicability at alpha = 0.005 : 2.564102564102564 %
Replicability at alpha = 0.001 : 2.564102564102564 %
*****************************************************************
FS_CT ridge age 200
0.2227858933204806 0.22609433664584924
Replicability at alpha = 0.05 : 94.31818181818183 %
Replicability at alpha = 0.01 : 65.9090909090909 %
Replicability at alpha = 0.005 : 52.27272727272727 %
Replicability at alpha = 0.001 : 26.136363636363637 %
************************************[...]<jupyter_text>Null scenario with random targetTo evaluate false positives with biased estimates<jupyter_code>%%time

random_state = 42
n_bootstrap = 100

features = {
    'FS_CT': "FS_*_Thck"
}

models = {
    'PCA_SVR': Pipeline([('pca', PCA(n_components=0.5)),
                         ('svr', SVR())])

}

# We aggregate all results here:
df = pd.DataFrame(columns=['connectivity','model','target','n','r_discovery_cv','r_discovery_overfit','r_replication','p_discovery_cv','p_discovery_overfit','p_replication'])

for feature_set in features:
    for model in models:
        for target_var in ['age', 'CogTotalComp_AgeAdj', 'PMAT24_A_CR', 'Flanker_AgeAdj', 'CardSort_AgeAdj', 'PicSeq_AgeAdj']:
            for sample_size in [50, 100, 200, 300, 'max']:

                print('*****************************************************************')
                print(feature_set, model, target_var, sample_size)

                X, y = create_data(target=target_var, feature_data=features[feature_set]) # gives a random y when target is None

                if sample_size=='max':
                    sample_size = int(len(y)/2)

                # create random seeds for each bootstrap iteration for reproducibility
                rng = np.random.default_rng(random_state)
                random_sates = rng.integers(np.iinfo(np.int32).max, size=n_bootstrap)

                # run bootstrap iterations in parallel, with shuffle_y=True
                r_discovery_cv, r_discovery_overfit, r_replication, p_discovery_cv, p_discovery_overfit, p_replication = zip(
                    *Parallel(n_jobs=-1)(
                    delayed(bootstrap_workhorse)(X, y, sample_size, models[model], seed, shuffle_y=True) for seed in random_sates))

                tmp_data_frame = pd.DataFrame({
                    'connectivity' : feature_set,
                    'model' : model,
                    'target' : target_var,
                    'n' : sample_size,
                    'r_discovery_cv': r_discovery_cv,
                    'r_discovery_overfit': r_discovery_overfit,
                    'r_replication': r_replication,
                    'p_discovery_cv': p_discovery_cv,
                    'p_discovery_overfit': p_discovery_overfit,
                    'p_replication': p_replication
                })

                #sns.scatterplot(x='r_replication', y='r_discovery_cv', data=tmp_data_frame)
                #plt.ylabel('in-sample (r)')
                #plt.xlabel('out-of-sample (r_pred)')
                #plt.show()
                print(tmp_data_frame.r_discovery_cv.mean(), tmp_data_frame.r_replication.mean())

                for alpha in [0.05, 0.01, 0.005, 0.001]:
                    print('Replicability at alpha =', alpha, ':',
                          (tmp_data_frame.loc[tmp_data_frame['p_discovery_cv']<alpha,'p_replication']<alpha).sum() / (tmp_data_frame['p_discovery_cv']<0.05).sum() * 100, '%')

                df = pd.concat((df, tmp_data_frame))
                df.reset_index(drop=True, inplace=True)
                df.to_csv('res/results_thck_null_PCA_SVR.csv')

df<jupyter_output>*****************************************************************
FS_CT PCA_SVR age 50
-0.006672762525657837 -0.007750463437365226
Replicability at alpha = 0.05 : 0.0 %
Replicability at alpha = 0.01 : 0.0 %
Replicability at alpha = 0.005 : 0.0 %
Replicability at alpha = 0.001 : 0.0 %
*****************************************************************
FS_CT PCA_SVR age 100
-0.0054298235725445355 0.007488203892819826
Replicability at alpha = 0.05 : 0.0 %
Replicability at alpha = 0.01 : 0.0 %
Replicability at alpha = 0.005 : 0.0 %
Replicability at alpha = 0.001 : 0.0 %
*****************************************************************
FS_CT PCA_SVR age 200
0.0006834925110597672 0.003009677488480728
Replicability at alpha = 0.05 : 20.0 %
Replicability at alpha = 0.01 : 0.0 %
Replicability at alpha = 0.005 : 0.0 %
Replicability at alpha = 0.001 : 0.0 %
*****************************************************************
FS_CT PCA_SVR age 300
-0.0021544089614625975 -0.006674668177288535
Replicabil[...]
{'fraction_non_alphanumeric': 0.1450635490963694, 'fraction_numerical': 0.05378881321392181, 'mean_word_length': 3.4452920143027415, 'pattern_counts': {'":': 1, '<': 35, '<?xml version=': 0, '>': 26, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '453699', 'license': 'MIT', 'n_tokens_mistral': 6575, 'n_tokens_neox': 5524, 'n_words': 1385, 'path': 'multivariate_BWAS_replicability_analysis_CT.ipynb', 'repo_name': 'spisakt/BWAS_comment', 'size': '18647'}