<jupyter_start><jupyter_text>MVPA and Searchlight with `nilearn`In this section we will show how you can use `nilearn` to perform multivariate pattern analysis (MVPA), either in a specific region of interest (ROI) or directly in the whole brain, for example with a Searchlight analysis. `nilearn`Although nilearn's visualizations are quite nice, its primary purpose was to facilitate machine learning in neuroimaging. It's in some sense the bridge between [nibabel](http://nipy.org/nibabel/) and [scikit-learn](http://scikit-learn.org/stable/). On the one hand, it reformats images to be easily passed to scikit-learn, and on the other, it reformats the results to produce valid nibabel images.DisclaimerFirst: This section is heavily based on the nilearn decoding tutorial.    Second: This section is not intended to teach machine learning, but to demonstrate a simple nilearn pipeline.  Setup<jupyter_code>import numpy as np
import nibabel as nb
import pandas as pd
from nilearn.image import resample_to_img, math_img, mean_img

%matplotlib inline
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Locating and verifying the machine learning datasetFirst, let's verify the location and size of the dataset we prepared in the previous notebook:<jupyter_code>func = '/home/neuro/workshop/notebooks/data/dataset_ML.nii.gz'
!nib-ls $func
# For quicker access, let's also create a functional mean image
img_mean = mean_img(func)
img_mean.shape<jupyter_output><empty_output><jupyter_text>Creating masksNow, let's create multiple masks to better restrict some of our analysis. What we'll need is a mask for the brain (and in this case exceptionally also the eyes), as well as a mask for three different regions of interests (ROI): The primary visual cortex (V1), the primary audiotry cortex (A1) and the eyes (EYE).<jupyter_code>from nilearn.plotting import plot_roi
def create_mask(img_roi, img_mean):
    
    # Resample mask to functional space
    img_resampled = resample_to_img(img_roi, img_mean)

    # Binarize ROI template
    data_binary = np.array(img_resampled.get_fdata()!=0, dtype=np.int8)

    # Dilate binary mask once
    data_dilated = binary_dilation(data_binary, iterations=1).astype(np.int8)

    # Save binary mask in NIfTI image
    mask = nb.Nifti1Image(data_dilated, img_resampled.affine, img_resampled.header)
    
    # Restrict mask only to voxels with functional values
    mask = math_img("(mean * mask)>0", mean=img_mean, mask=mask)
    mask.set_data_dtype('i1')
    
    return mask<jupyter_output><empty_output><jupyter_text>Now that we're ready, let's create some masks!<jupyter_code># Create global mask containing brain and eyes
brain = '/home/neuro/workshop/notebooks/data/templates/MNI152_T1_1mm_brain.nii.gz'
eyes = '/home/neuro/workshop/notebooks/data/templates/MNI152_T1_1mm_eye.nii.gz'
img_roi = math_img("(img1 + img2)>0", img1=brain, img2=eyes)
img_resampled = resample_to_img(img_roi, img_mean)
mask_global = math_img("((mean!=0)*img)>0.5", mean=img_mean, img=img_resampled)

# Plotting global mask
anat = '/home/neuro/workshop/notebooks/data/templates/MNI152_T1_1mm.nii.gz'
plot_roi(mask_global, img_mean, cmap='Paired', dim=-.5, draw_cross=False, annotate=False);
# Create eye mask containing both eyes
img_roi = math_img("img>0", img=eyes)
img_resampled = resample_to_img(img_roi, img_mean)
mask_eye = math_img("((mean!=0)*img)>0.5", mean=img_mean, img=img_resampled)

# Plotting global mask
anat = '/home/neuro/workshop/notebooks/data/templates/MNI152_T1_1mm.nii.gz'
plot_roi(mask_eye, img_mean, cmap='Paired', dim=-.5, draw_cross=False, annotate=False);<jupyter_output><empty_output><jupyter_text>To get region specific masks for V1 and A1, we will use [FSL's Harvard-Oxford atlas](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Atlases) as provided by [AtlasReader](https://github.com/miykael/atlasreader).<jupyter_code># Load labels from AtlasReader
atlas_path = '/opt/miniconda-latest/envs/neuro/lib/python3.7/site-packages/atlasreader/data/atlases/'
atlas_labels = pd.read_csv(atlas_path + 'labels_harvard_oxford.csv', index_col=0)
atlas_labels.head()
# Load atlas from AtlasReader
img_atlas = nb.load(atlas_path + 'atlas_harvard_oxford.nii.gz')
img_atlas.shape
# Load V1 ROI labels
idx_v1 = [92, 93, 94, 95]
display(atlas_labels.iloc[idx_v1])
# Create V1 mask
roi_v1 = mean_img([img_atlas.slicer[..., v] for v in idx_v1])
img_roi = math_img("img>0", img=roi_v1)
img_resampled = resample_to_img(img_roi, img_mean)
mask_v1 = math_img("((mean!=0)*img)>0.5", mean=img_mean, img=img_resampled)

# Plotting V1 mask
plot_roi(mask_v1, img_mean, cmap='Paired', dim=-.5, draw_cross=False, annotate=False);
# Load A1 ROI labels
idx_a1 = [88, 89, 90, 91]
display(atlas_labels.iloc[idx_a1])
# Create A1 mask
roi_a1 = mean_img([img_atlas.slicer[..., a] for a in idx_a1])
img_roi = math_img("img>0", img=roi_a1)
img_resampled = resample_to_img(img_roi, img_mean)
mask_a1 = math_img("((mean!=0)*img)>0.5", mean=img_mean, img=img_resampled)

# Plotting V1 mask
plot_roi(mask_a1, img_mean, cmap='Paired', dim=-.5, draw_cross=False, annotate=False);<jupyter_output><empty_output><jupyter_text>Labels and chunksNow that we have the dataset and the masks at hand, we're almost ready for some machine learning. Only thing still missing are the **target labels** (i.e. which samples / time-points are "eyes open" condition and which ones aren't, and **chunks/subject_identifier** to prevent data leakage between model training and validation.From the last section of the [Machine Learning Preparation](machine_learning_preparation.ipynb) notebook, we know that we have a total of 384 volumes in our `dataset_ML.nii.gz` file and that it's always 4 volumes of the condition `eyes closed`, followed by 4 volumes of the condition `eyes open`, etc. Therefore our labels should be as follows:<jupyter_code>labels = np.ravel([[['closed'] * 4, ['open'] * 4] for i in range(48)])
labels[:20]<jupyter_output><empty_output><jupyter_text>As mentioned above, the `chunks` variable is important if we want to do for example something like cross-validation. In our case we would ideally create 48 chunks, one for each subject. But because a cross-validation of 48 chunks takes very long, let's just create 6 chunks, containing always 8 subjects, i.e. 64 volumes in one group:<jupyter_code>chunks = np.ravel([[i] * 64 for i in range(6)])
chunks[:100]<jupyter_output><empty_output><jupyter_text>One way to do cross-validation is the so called **Leave-one-out cross-validation**. This approach trains on `(n - 1)` chunks, and classifies the remaining chunk, and repeats this for every chunk, also called **fold**. Therefore, a 6-fold cross-validation is one that divides the whole data into 6 different chunks.  ROI-based decoding analysisNow that the masks are ready, let's look at a ROI based decoding analysis. In other words, we'll use the response patterns of voxels in the mask to predict if the eyes were **closed** or **open** during a resting-state fMRI recording.To setup the procedure let's perform the analysis first just for the **ROI of V1**.<jupyter_code>mask2use = mask_v1<jupyter_output><empty_output><jupyter_text>Masking and Un-masking dataFor the classification with `nilearn`, we need our functional data in a 2D, sample-by-voxel matrix. To get that, we'll select all the voxels within our `mask`. To do this, we can use the `NiftiMasker` object, which applies a mask to a dataset and returns the masked voxels flattened to a 2D matrix.<jupyter_code># Mask ML dataset with V1 ROI mask
from nilearn.input_data import NiftiMasker
masker = NiftiMasker(mask_img=mask2use, standardize=False, detrend=False,
                     memory="nilearn_cache", memory_level=2)
samples = masker.fit_transform(func)
print(samples)<jupyter_output><empty_output><jupyter_text>Its shape corresponds to the number of time-points times the number of voxels in the mask.<jupyter_code>print(samples.shape)<jupyter_output><empty_output><jupyter_text>To recover the original data shape (giving us a masked and z-scored BOLD series), we simply use the masker's inverse transform:<jupyter_code>masked_epi = masker.inverse_transform(samples)<jupyter_output><empty_output><jupyter_text>Let's now visualize the masked epi.<jupyter_code>from nilearn.plotting import plot_stat_map

max_zscores = math_img("np.abs(img).max(axis=3)", img=masked_epi)
plot_stat_map(max_zscores, bg_img=anat, dim=-.5,
              draw_cross=False, annotate=False, colorbar=False,
              title='Maximum Amplitude per Voxel in Mask')<jupyter_output><empty_output><jupyter_text>Train/test split of datasetTo be able to test how well our model performs, we need to a test set, i.e. a dataset that wasn't used to train the model. The easiest solution to this is to use scikit-learns `train_test_split`.<jupyter_code>from sklearn.model_selection import train_test_split

# Split the data into training and test sets
X_train_v1, X_test_v1, y_train_v1, y_test_v1, c_train, c_test = train_test_split(
    samples, labels=='open', chunks, test_size=64, random_state=0, shuffle=False)

print('Shapes of X:', X_train_v1.shape, X_test_v1.shape)
print('Shapes of y:', y_train_v1.shape, y_test_v1.shape)<jupyter_output><empty_output><jupyter_text>Training and testing a linear support vector machine (linear SVC) model for V1 ROI<jupyter_code>from sklearn.svm import LinearSVC
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.model_selection import GridSearchCV
# Create pipeline
pipe = Pipeline([
    ('scaler', StandardScaler()),
    ('svc', LinearSVC(multi_class='ovr', penalty='l2',
                      loss='squared_hinge', max_iter=500))
])
# Create parameter grid to be explored during grid search
grid = {
    'svc__C': np.logspace(-6, 2, num=40),
}
# Create grid search object with a cross-validation with LeaveOneGroupOut
from sklearn.model_selection import LeaveOneGroupOut
grid_cv_v1 = GridSearchCV(pipe, grid, cv=LeaveOneGroupOut(),
                          return_train_score=True, refit=True, n_jobs=-1, verbose=1)
# Train the model and find optimal hyperparameter
grid_cv_v1.fit(X=X_train_v1, y=y_train_v1, groups=c_train);<jupyter_output><empty_output><jupyter_text>To better understand the model performance during the hyperparamter fine tuning, let's write a supportive plotting function.<jupyter_code>def plot_hyperparam_fitting(cv_results):

    # Store grid search parameters and outcomes in dataframe
    df_pred = pd.DataFrame(cv_results)
    columns = [c for c in df_pred.columns if 'time' not in c
               and 'split' not in c
               and 'rank' not in c
               and c!='params']
    df_pred = df_pred[columns].sort_values('mean_test_score', ascending=False)
    display(df_pred.head())
    
    # Plot the model fit information
    df_plot = df_pred.sort_values('param_svc__C')

    # Exsract relevant modelling metrics
    train_scores = df_plot['mean_train_score']
    valid_scores = df_plot['mean_test_score']
    std_tr = df_plot['std_train_score']
    std_va = df_plot['std_test_score']

    plt.figure(figsize=(12, 4))
    Cs = df_plot['param_svc__C']
    plt.semilogx(Cs, train_scores, label='Training Set')
    plt.semilogx(Cs, valid_scores, label='Validation Set')

    # Add marker and text for best score
    max_id = np.argmax(valid_scores)
    x_pos = Cs.iloc[max_id]
    y_pos = valid_scores.iloc[max_id]
    txt = '{:0.4f}'.format(y_pos)
    plt.scatter(x_pos, y_pos, marker='x', c='red', zorder=10)
    plt.text(x_pos, y_pos, txt, fontdict={'size': 18})

    # Quantify variance with ±std curves
    plt.fill_between(Cs.astype(float), train_scores-std_tr, train_scores+std_tr, alpha=0.3)
    plt.fill_between(Cs.astype(float), valid_scores-std_va, valid_scores+std_va, alpha=0.3)
    plt.ylabel('Performance metric')
    plt.xlabel('Model parameter')

    # Adjust x-lim, y-lim, add legend and adjust layout
    plt.legend()
    plt.show()
# Plot hyperparameter fine tuning outcome
plot_hyperparam_fitting(grid_cv_v1.cv_results_)<jupyter_output><empty_output><jupyter_text>Now that the model is trained and we've verified that the fit seems reasonable, let's compute the models performance on the withheld test set.<jupyter_code>score_test_v1 = grid_cv_v1.score(X_test_v1, y_test_v1)
print('Test score based on V1:', score_test_v1)<jupyter_output><empty_output><jupyter_text>Given that the chance level is 50%, having ~73% accuracy on predicting if a person has eyes opened or closed during a resting state recording, only based on V1 information is rather good!  Training and testing a linear support vector machine (linear SVC) model for A1 ROIWhat about the **regions A1**?<jupyter_code># Mask ML dataset with A1 ROI mask
mask2use = mask_a1
masker = NiftiMasker(mask_img=mask2use, standardize=False, detrend=False,
                     memory="nilearn_cache", memory_level=2)
samples = masker.fit_transform(func)
# Split the data into training and test sets
X_train_a1, X_test_a1, y_train_a1, y_test_a1, c_train, c_test = train_test_split(
    samples, labels=='open', chunks, test_size=64, random_state=0, shuffle=False)
# Create grid search object with a cross-validation with LeaveOneGroupOut
grid_cv_a1 = GridSearchCV(pipe, grid, cv=LeaveOneGroupOut(),
                          return_train_score=True, refit=True, n_jobs=-1, verbose=1)
# Train the model and find optimal hyperparameter
grid_cv_a1.fit(X=X_train_a1, y=y_train_a1, groups=c_train);
# Plot hyperparameter fine tuning outcome
plot_hyperparam_fitting(grid_cv_a1.cv_results_)
score_test_a1 = grid_cv_a1.score(X_test_a1, y_test_a1)
print('Test score based on A1:', score_test_a1)<jupyter_output><empty_output><jupyter_text>Below chance level,... It makes sense that A1 would be somewhere around chance level, but this value is nonetheless a bit low. Let's come back to this later.  Training and testing a linear support vector machine (linear SVC) model for eye ROILet's also check an obvious region to detect if eyes were opened or closed, the **eyes themselves**.<jupyter_code># Mask ML dataset with eye ROI mask
mask2use = mask_eye
masker = NiftiMasker(mask_img=mask2use, standardize=False, detrend=False,
                     memory="nilearn_cache", memory_level=2)
samples = masker.fit_transform(func)
# Split the data into training and test sets
X_train_eye, X_test_eye, y_train_eye, y_test_eye, c_train, c_test = train_test_split(
    samples, labels=='open', chunks, test_size=64, random_state=0, shuffle=False)
# Create grid search object with a cross-validation with LeaveOneGroupOut
grid_cv_eye = GridSearchCV(pipe, grid, cv=LeaveOneGroupOut(),
                           return_train_score=True, refit=True, n_jobs=-1, verbose=1)
# Train the model and find optimal hyperparameter
grid_cv_eye.fit(X=X_train_eye, y=y_train_eye, groups=c_train);
# Plot hyperparameter fine tuning outcome
plot_hyperparam_fitting(grid_cv_eye.cv_results_)
score_test_eye = grid_cv_eye.score(X_test_eye, y_test_eye)
print('Test score based on A1:', score_test_eye)<jupyter_output><empty_output><jupyter_text>That's a very high score. Then again, not surprising for the eyes themselves.  Conclusion of the ROI-based decoding analysis<jupyter_code>df_res = pd.DataFrame([0.5, score_test_v1, score_test_a1, score_test_eye],
                      columns=['Accuracy'], index=['chance', 'V1', 'A1', 'eyes']).T.round(2)*100
df_res
plt.figure(figsize=(7, 4))
plt.bar(df_res.columns, df_res.values.T.ravel())
plt.title('Decoding accuarcy for certain ROI');
print()<jupyter_output><empty_output><jupyter_text>Picking one of the two classes "eyes-closed" or "eyes-open" would lead to a prediction accuracy of 50%, i.e. chance-level. Not surprisingly, the signal recorded within the eyes lead to the highest prediction accuracy, followed by the primary visual cortex, while the hopefully unrelated primary auditory cortex doesn't provide better prediction accuracy then the chance-level.  Permutation testingThe approach above, with using a train and test set is important in machine learning to establish an optimized model and to judge how well it can generalize to a new dataset (here the test set).One way to test the quality of the prediction accuracy is to run the cross-validation multiple times, but permutate the labels of the volumes randomly. Afterward we can compare the accuracy value of the correct labels to the ones with the random / false labels. Luckily `Scikit-learn` already has a function that does this for us. So let's do it.**Note**: As a classifier, we will chose the best classifier/estimator from our grid search above. Additionally, we chose the number of iterations under `n_permutations` for the permutation testing rather low, to reduce computation time as well. This value should ideally be higher.  Permuation testing for V1-ROI<jupyter_code>from sklearn.model_selection import permutation_test_score

# Mask ML dataset with v1 ROI mask
mask2use = mask_v1
masker = NiftiMasker(mask_img=mask2use, standardize=False, detrend=False,
                     memory="nilearn_cache", memory_level=2)
samples = masker.fit_transform(func)

# Get best classifier from gridsearch
clf = grid_cv_v1.best_estimator_

# Run the permuation cross-validation
null_cv_scores = permutation_test_score(estimator=clf,
                                        X=samples,
                                        y=labels,
                                        groups=chunks,
                                        cv=LeaveOneGroupOut(),
                                        n_permutations=200,
                                        n_jobs=-1,
                                        verbose=1)<jupyter_output><empty_output><jupyter_text>So, let's take a look at the results:<jupyter_code># Plotting the null distribution
plt.title('Null-hypothesis prediction accuracy')
plt.hist(null_cv_scores[1], bins=25, density=True);
plt.vlines(null_cv_scores[0], 0, 25, color='r')
print('Red line indicates accuracy of actual labels.')

print('Prediction accuracy: %.02f' % (null_cv_scores[0] * 100),
      'p-value: %.04f' % (null_cv_scores[2]),
      sep='\n')<jupyter_output><empty_output><jupyter_text>Great! This means... Using resting-state fMRI images, we can predict if a person had their eyes open or closed with an accuracy significantly above chance level, when we only look at the activation from V1!  Permuation testing for A1-ROI<jupyter_code>from sklearn.model_selection import permutation_test_score

# Mask ML dataset with a1 ROI mask
mask2use = mask_a1
masker = NiftiMasker(mask_img=mask2use, standardize=False, detrend=False,
                     memory="nilearn_cache", memory_level=2)
samples = masker.fit_transform(func)

# Get best classifier from gridsearch
clf = grid_cv_a1.best_estimator_

# Run the permuation cross-validation
null_cv_scores = permutation_test_score(estimator=clf,
                                        X=samples,
                                        y=labels,
                                        groups=chunks,
                                        cv=LeaveOneGroupOut(),
                                        n_permutations=200,
                                        n_jobs=-1,
                                        verbose=1)<jupyter_output><empty_output><jupyter_text>So, let's take a look at the results:<jupyter_code># Plotting the null distribution
plt.title('Null-hypothesis prediction accuracy')
plt.hist(null_cv_scores[1], bins=25, density=True);
plt.vlines(null_cv_scores[0], 0, 25, color='r')
print('Red line indicates accuracy of actual labels.')

print('Prediction accuracy: %.02f' % (null_cv_scores[0] * 100),
      'p-value: %.04f' % (null_cv_scores[2]),
      sep='\n')<jupyter_output><empty_output><jupyter_text>As already seen before, using A1 to decode if somebody had eyes open or closed during the resting state recording doesn't work.  Permuation testing for EYE-ROI<jupyter_code>from sklearn.model_selection import permutation_test_score

# Mask ML dataset with eye ROI mask
mask2use = mask_eye
masker = NiftiMasker(mask_img=mask2use, standardize=False, detrend=False,
                     memory="nilearn_cache", memory_level=2)
samples = masker.fit_transform(func)

# Get best classifier from gridsearch
clf = grid_cv_eye.best_estimator_

# Run the permuation cross-validation
null_cv_scores = permutation_test_score(estimator=clf,
                                        X=samples,
                                        y=labels,
                                        groups=chunks,
                                        cv=LeaveOneGroupOut(),
                                        n_permutations=200,
                                        n_jobs=-1,
                                        verbose=1)<jupyter_output><empty_output><jupyter_text>So, let's take a look at the results:<jupyter_code># Plotting the null distribution
plt.title('Null-hypothesis prediction accuracy')
plt.hist(null_cv_scores[1], bins=25, density=True);
plt.vlines(null_cv_scores[0], 0, 25, color='r')
print('Red line indicates accuracy of actual labels.')

print('Prediction accuracy: %.02f' % (null_cv_scores[0] * 100),
      'p-value: %.04f' % (null_cv_scores[2]),
      sep='\n')<jupyter_output><empty_output><jupyter_text>And also here, the previous assumption gets supported. The eyes are of course a very good location to predict if a participant had eyes open or closed during a resting state recording.  Which region is driving the classification when we don't have prior knowledge?In the previous example, we had some intuition/hypothesis which regions to look into. But this might not always be the case. So how can we get some ideas about which regions might be driving the classification accuracy, when one is looking at the whole brain?There are many different ways to figure out which region is important for classification, but let us introduce you to two different approaches that you can use in `nilearn`: `SpaceNet` and  `Searchlight`**Note:** For both of these examples we will take the full brain mask prepared above.<jupyter_code>plot_roi(mask_global, img_mean, cmap='Paired', dim=-.5, draw_cross=False, annotate=False);<jupyter_output><empty_output><jupyter_text>SpaceNet: decoding with spatial structure for better mapsSpaceNet implements spatial penalties which improve brain decoding power as well as decoder maps. The results are brain maps which are both sparse (i.e regression coefficients are zero everywhere, except at predictive voxels) and structured (blobby). For more detail, check out `nilearn`'s section about [SpaceNet](http://nilearn.github.io/decoding/space_net.html).To train a SpaceNet on our data, let's first split the data into a training set (chunk 0-4) and a test set (chunk 5).<jupyter_code># Split the data into training and test sets
idx_train, idx_test, y_train, y_test, c_train, c_test = train_test_split(
    range(len(samples)), labels=='open', chunks, test_size=64, random_state=0, shuffle=False)
# Split the NIfTI image into train and test set
from nilearn.image import index_img
X_train = index_img(func, idx_train)
X_test = index_img(func, idx_test)<jupyter_output><empty_output><jupyter_text>Now we can fit the SpaceNet to our data with a TV-l1 penalty. ***Note*** again, that we reduced the number of `max_iter` to have a quick computation. In a realistic case this value should be around 1000.<jupyter_code>from nilearn.decoding import SpaceNetClassifier

# Fit model on train data and predict on test data
decoder = SpaceNetClassifier(penalty='tv-l1',
                             mask=mask_global,
                             max_iter=10,
                             cv=5,
                             standardize=True,
                             memory="nilearn_cache",
                             memory_level=2,
                             verbose=0)

decoder.fit(X_train, y_train)<jupyter_output><empty_output><jupyter_text>Now that the `SpaceNet` is fitted to the training data. Let's see how well it does in predicting the test data.<jupyter_code># Predict the labels of the test data
y_pred = decoder.predict(X_test)
# Retrun average accuracy
accuracy = (y_pred == y_test).mean() * 100.
print("\nTV-l1  classification accuracy : %g%%" % accuracy)<jupyter_output><empty_output><jupyter_text>Again above 80% prediction accuracy? But we wanted to know what's driving this prediction. So let's take a look at the fitting coefficients.<jupyter_code>from nilearn.plotting import plot_stat_map, show
coef_img = decoder.coef_img_
# Plotting the searchlight results on the glass brain
from nilearn.plotting import plot_glass_brain
plot_glass_brain(coef_img, black_bg=True, colorbar=True, display_mode='lyrz', symmetric_cbar=False,
                 cmap='magma', title='graph-net: accuracy %g%%' % accuracy)<jupyter_output><empty_output><jupyter_text>Cool! As expected the visual cortex (in the back of the head) and the eyes are driving the classification!  Searchlight: finding voxels containing informationNow the next question is: How high would the prediction accuracy be if we only take one small region to do the classification?To answer this question we can use something that is called a **Searchlight** approach. The searchlight approach was first proposed by [Kriegeskorte et al., 2006](https://pdfs.semanticscholar.org/985c/ceaca8606443f9129616a26bbbbf952f2d7f.pdf). It is a widely used approach for the study of the fine-grained patterns of information in fMRI analysis. Its principle is relatively simple: a small group of neighboring features is extracted from the data, and the prediction function is instantiated on these features only. The resulting prediction accuracy is thus associated with all the features within the group, or only with the feature on the center. This yields a map of local fine-grained information, that can be used for assessing hypothesis on the local spatial layout of the neural code under investigation.You can do a searchlight analysis in `nilearn` as follows:<jupyter_code>from nilearn.decoding import SearchLight
# Specify the mask in which the searchlight should be performed
mask = mask_global
# Specify the classifier to use
# For presentation purpose, let's use a GaussainNB classifier to have rather small computation time
from sklearn.naive_bayes import GaussianNB
clf = GaussianNB()
# Specify the radius of the searchlight sphere  that will scan the volume
# (the bigger the longer the computation)
sphere_radius = 8  # in mm<jupyter_output><empty_output><jupyter_text>Now we're ready to create the searchlight object.<jupyter_code># Create searchlight object
sl = SearchLight(mask_global,
                 process_mask_img=mask,
                 radius=sphere_radius,
                 estimator=clf,
                 cv=LeaveOneGroupOut(),
                 n_jobs=-1,
                 verbose=1)
# Run the searchlight algorithm
sl.fit(nb.load(func), labels, groups=chunks)<jupyter_output><empty_output><jupyter_text>That took a while. So let's take a look at the results.<jupyter_code># First we need to put the searchlight output back into an MRI image
from nilearn.image import new_img_like
searchlight_img = new_img_like(func, sl.scores_)<jupyter_output><empty_output><jupyter_text>Now we can plot the results. Let's plot it once on the glass brain and once from the side. For better interpretation on where the peaks are, let's set a minimum accuracy threshold of 60%.<jupyter_code>from nilearn.plotting import plot_glass_brain
plot_glass_brain(searchlight_img, black_bg=True, colorbar=True, display_mode='lyrz',
                 threshold=0.6, cmap='magma', title='Searchlight Prediction Accuracy')
from nilearn.plotting import plot_stat_map
plot_stat_map(searchlight_img, cmap='magma', bg_img=anat, colorbar=True,
              display_mode='x', threshold=0.6, cut_coords=[0, 6, 12, 18],
              title='Searchlight Prediction Accuracy');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08685689735419426, 'fraction_numerical': 0.01371236260785507, 'mean_word_length': 4.381888246628131, 'pattern_counts': {'":': 0, '<': 137, '<?xml version=': 0, '>': 146, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '34', 'dir': 'jupyter-structured-clean-dedup', 'id': '171022', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 8948, 'n_tokens_neox': 8445, 'n_words': 2756, 'path': 'workshop/notebooks/05b_machine_learning_nilearn.ipynb', 'repo_name': 'jooh/workshop_pybrain', 'size': '27931'}