fMRI Multiverse

fMRI Multiverse#

Note: You can download this individual file as a Jupyter Notebook (.ipynb) file by clicking the download button at the top.

Pre-load the data#

For the fMRI multiverse analysis example, we will use data from 100 participants of the Autism Brain Imaging Data Exchange (ABIDE) preprocessed connectomes dataset.

  • Run the following code cell to download the data from nilearn into the data/ folder. To speed things up, you can also use one of our USB sticks and copy the entire data/ folder into the directory of this script.

  • Inspect the following code cell and visit the website of the dataset to familiarize yourself with the data. What do the individual options mean, and are they reasonable decisions for a multiverse analysis?

from tqdm import tqdm
import itertools
from nilearn import datasets

pipeline = ["cpac", "ccs", "dparsf", "niak"]                    # Preprocessing pipelines
band_pass = [True, False]                                       # Band-pass filtering   
global_signal = [True, False]                                   # Global signal regression 
parcellation = ["rois_aal", "rois_cc200", "rois_dosenbach160"]  # Parcellated time series data

# Download the ABIDE dataset with all combinations of the decision points
abide_dataset = {} # store in a dict (only needed to print some information at the end of this cell)
for pipe, bp, gsr, parc in tqdm(itertools.product(pipeline, band_pass, global_signal, parcellation)):
    bunch = datasets.fetch_abide_pcp(
        data_dir="./data", n_subjects=100, quality_checked=True, verbose=0,
        pipeline=pipe, derivatives=parc, band_pass_filtering=bp, global_signal_regression=gsr)
    abide_dataset[(pipe, bp, gsr, parc)] = bunch

print(f"Available pipelines: {list(abide_dataset.keys())}")
print(f"Number of subjects:  {len(abide_dataset[('cpac', True, True, 'rois_aal')].phenotypic)}")
print(f"Class distribution:  {abide_dataset[('cpac', True, True, 'rois_aal')].phenotypic['DX_GROUP'].value_counts()}")

48it [00:25,  1.91it/s]

Available keys:     [('cpac', True, True, 'rois_aal'), ('cpac', True, True, 'rois_cc200'), ('cpac', True, True, 'rois_dosenbach160'), ('cpac', True, False, 'rois_aal'), ('cpac', True, False, 'rois_cc200'), ('cpac', True, False, 'rois_dosenbach160'), ('cpac', False, True, 'rois_aal'), ('cpac', False, True, 'rois_cc200'), ('cpac', False, True, 'rois_dosenbach160'), ('cpac', False, False, 'rois_aal'), ('cpac', False, False, 'rois_cc200'), ('cpac', False, False, 'rois_dosenbach160'), ('ccs', True, True, 'rois_aal'), ('ccs', True, True, 'rois_cc200'), ('ccs', True, True, 'rois_dosenbach160'), ('ccs', True, False, 'rois_aal'), ('ccs', True, False, 'rois_cc200'), ('ccs', True, False, 'rois_dosenbach160'), ('ccs', False, True, 'rois_aal'), ('ccs', False, True, 'rois_cc200'), ('ccs', False, True, 'rois_dosenbach160'), ('ccs', False, False, 'rois_aal'), ('ccs', False, False, 'rois_cc200'), ('ccs', False, False, 'rois_dosenbach160'), ('dparsf', True, True, 'rois_aal'), ('dparsf', True, True, 'rois_cc200'), ('dparsf', True, True, 'rois_dosenbach160'), ('dparsf', True, False, 'rois_aal'), ('dparsf', True, False, 'rois_cc200'), ('dparsf', True, False, 'rois_dosenbach160'), ('dparsf', False, True, 'rois_aal'), ('dparsf', False, True, 'rois_cc200'), ('dparsf', False, True, 'rois_dosenbach160'), ('dparsf', False, False, 'rois_aal'), ('dparsf', False, False, 'rois_cc200'), ('dparsf', False, False, 'rois_dosenbach160'), ('niak', True, True, 'rois_aal'), ('niak', True, True, 'rois_cc200'), ('niak', True, True, 'rois_dosenbach160'), ('niak', True, False, 'rois_aal'), ('niak', True, False, 'rois_cc200'), ('niak', True, False, 'rois_dosenbach160'), ('niak', False, True, 'rois_aal'), ('niak', False, True, 'rois_cc200'), ('niak', False, True, 'rois_dosenbach160'), ('niak', False, False, 'rois_aal'), ('niak', False, False, 'rois_cc200'), ('niak', False, False, 'rois_dosenbach160')]
Number of subjects: 100
Class distribution: DX_GROUP
1    50
2    50
Name: count, dtype: int64

Exercise: Create and run the multiverse#

The multiverse which we will implement is is similar to the one perfomed by Dafflon et al. 2022. The main difference is that we will only apply two connectivity methods and also only use data for 100 participants to reduce memory and cpu load.

In short, we will predict an autism diagnosis based on static functional connectivity estimates. Available decision points for the preprocessed fMRI time series data are the following:

  • Preprocessing pipeline ('cpac', 'ccs', 'dparsf', 'niak')

  • Band pass filtering (True or False)

  • Global signal regression (True or False) -> If false, standard motion regression was performed

  • Parcellation atlas ('rois_aal', 'rois_cc200', 'rois_dosenbach160')

For the connectivity measure, the two methods from the comet toolbox are already included:

  • Pearson correlation (comet.connectivity.Static_Pearson)

  • Partial correlation (comet.connectivity.Static_Partial)

from comet import multiverse

forking_paths = {
    "pipeline": [...],          # Preprocessing pipelines
    "parcellation": [...],      # Parcellated time series data
    "band_pass": [...],         # Band-pass filtering   
    "global_signal": [...],     # Global signal regression     
    "connectivity":[            # Functional connectivity method
        {"name": "pearson", "func": "comet.connectivity.Static_Pearson(ts).estimate()"},
        {"name": "partial", "func": "comet.connectivity.Static_Partial(ts).estimate()"}]
}

def analysis_template():
    import comet
    import numpy as np
    from nilearn import datasets
    from sklearn.pipeline import Pipeline
    from sklearn.preprocessing import StandardScaler
    from sklearn.linear_model import LogisticRegression
    from sklearn.model_selection import StratifiedKFold, cross_val_score

    # Get data (if it is preloaded this will skip the download)
    data = datasets.fetch_abide_pcp(data_dir="./data", n_subjects=100, quality_checked=True, verbose=0,
                                    pipeline= ...,
                                    derivatives= ...,
                                    band_pass_filtering= ...,
                                    global_signal_regression= ...)

    time_series = data[...]
    diagnosis = data["phenotypic"]["DX_GROUP"]

    # Calculate FC
    tri_ix = None
    features = []

    for ts in time_series:
        FC = ...

        if tri_ix == None:
            tri_ix = np.triu_indices_from(FC, k=1)
        
        feat_vec = FC[tri_ix]
        features.append(feat_vec)

    # Prepare features (FC estimates) and target (autism/control)
    X = np.vstack(features)
    X[np.isnan(X)] = 0.0
    y = np.array(diagnosis)

    # Classification model
    model = Pipeline([('scaler', StandardScaler()), ('reg', LogisticRegression(penalty='l2'))])
    cv = StratifiedKFold(n_splits=10)
    accuracies = cross_val_score(model, X, y, cv=cv, scoring='accuracy')

    # Save the results
    comet.utils.save_universe_results({"accuracy": accuracies})

# Create and run the multiverse analysis
...

Discussion#

  • Discuss the specification curve. How would you interpret the results? What could be the next steps?

  • Do you see any statistical issues with the classification model?

Contents