<jupyter_start><jupyter_text>Crash course on Nilearn and MVPA (Multi-Variate Pattern Analysis)In this lesson, we will cover some popular methods in MVPA in the context of brain imaging. Namely,* Searchlight [Nikolaus et al. 2006 (PNAS)]* Space-Net [Dohmatob et al. 2014 (PRNI)] Load Haxby 2001 data<jupyter_code>from nilearn.datasets import fetch_haxby
data_files = fetch_haxby()
# Load Target labels
import numpy as np
labels = np.recfromcsv(data_files.session_target[0], delimiter=" ")

labels<jupyter_output><empty_output><jupyter_text>Restrict to face and house conditions<jupyter_code>target = labels['labels']
condition_mask = np.logical_or(target == b"face", target == b"house")<jupyter_output><empty_output><jupyter_text>Split data into train and test pools<jupyter_code>condition_mask_train = np.logical_and(condition_mask, labels['chunks'] <= 6)
condition_mask_test = np.logical_and(condition_mask, labels['chunks'] > 6)
from nilearn.image import index_img
from nilearn.masking import compute_epi_mask
func_filenames = data_files.func[0]
X_train = index_img(func_filenames, condition_mask_train)
X_test = index_img(func_filenames, condition_mask_test)
y_train = target[condition_mask_train]
y_test = target[condition_mask_test]<jupyter_output><empty_output><jupyter_text>Compute brain mask and background EPI image<jupyter_code>from nilearn.image import mean_img
mask_img = compute_epi_mask(func_filenames)
background_img = mean_img(func_filenames)
background_img.to_filename("bg_img.nii.gz")<jupyter_output><empty_output><jupyter_text>Searchlight============ Prepare masks* mask_img is the original mask* process_mask_img is a subset of mask_img, it contains the voxels that should be processed (we only keep the slice z = 26 and the back of the brain to speed up computation)<jupyter_code>from nilearn.decoding import SearchLight
from nilearn.image import new_img_like
process_mask = mask_img.get_data().astype(np.int)
picked_slice = 27
process_mask[..., (picked_slice + 1):] = 0
process_mask[..., :picked_slice] = 0
process_mask[:, 30:] = 0
process_mask_img = new_img_like(mask_img, process_mask)<jupyter_output><empty_output><jupyter_text>Now, fit a Searchlight model on your data<jupyter_code>searchlight = SearchLight(mask_img, process_mask_img=process_mask_img, radius=5.6, verbose=1)
searchlight.fit(X_train, y_train)<jupyter_output>[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    8.4s finished<jupyter_text>Plot the model coefficients to see whats really going on.<jupyter_code>%matplotlib inline
from nilearn.plotting import plot_stat_map, show
coef_img = new_img_like(background_img, searchlight.scores_)
plot_stat_map(coef_img, background_img,
              title="Search Light",
              cut_coords=(-34, -16), display_mode="yz")<jupyter_output><empty_output><jupyter_text>SpaceNet========* This is a family of priors which simultaneously impose sparsity and structure on the maps of regressor coefficients.* It includes GraphNet [Grosenick et al 2013], TV-L1 [Gramfort et al. 2013 (PRNI)], etc.![caption](haxby_comp.png)Formally, the model is\begin{equation}\text{minimize}_{\bf{w} \in \mathbb R^p}\mathcal L(y, X\bf{w}) + \alpha (\rho \|\bf{w}\|_1 + (1-\rho) \Omega(\bf{w}))\end{equation}where* $\bf{w} \in \mathbb R^p$ is the vector of model coefficients and $p$ is the number of voxels in the brain mask (say, $p \sim 10^6$ in a full brain mask!)* $\mathcal L(y, X\bf{w})$ is the loss term. For example $\mathcal L(y, X\bf{w}) = \frac{1}{2}\|X\bf{w} - y\|_2^2$ in regression settings, equals the logistic loss in classification settings, etc.* $\alpha > 0$ is the regularization parameter and controls and controls the relative importance of the penalty and loss terms. $\rho \in [0, 1]$ is a mixing parameter and controls the tradeoff between sparsity and structure. For example  $\rho = 1$ corresponds to the classical LASSO.* The penalty term $\|\bf{w}\|_1$ induces sparsity on the model coefficients $\bf{w}$: all coefficients are forced to zero, except the most "predictive" ones!* $\Omega(\bf{w})$ is a structure inducing penalty on the model coefficients:For example \begin{equation}\Omega(\bf{w}) = \frac{1}{2}\|\nabla \bf{w}\|_F^2 := \sum_{j=1}^p|(\nabla_x \bf{w})_j|^2 + |(\nabla_y \bf{w})_j|^2 + |(\nabla_z \bf{w})_j|^2\end{equation} in GraphNet / S-Lasso and for TV-L1, we have\begin{equation}\Omega(\bf{w}) = \|\bf{w}\|_{TV} := \sum_{j=1}^p\sqrt{|(\nabla_x \bf{w})_j|^2 + |(\nabla_y \bf{w})_j|^2 + |(\nabla_z \bf{w})_j|^2}\end{equation}* GraphNet is Tikhonov regularization on the spatial gradient of the model coefficients $\bf{w}$. It will produce smooth blobs, but no sharp edges.* TV semi-norm promotes smooth blobs with sharp edges (i.e has a segmentation property, which is a very good thing). However, the resulting model is much harder to optimize than GraphNet. **GraphNet / S-Lasso**<jupyter_code>from nilearn.decoding import SpaceNetClassifier

decoder = SpaceNetClassifier(memory="nilearn_cache", memory_level=2, penalty='graph-net', verbose=1, loss="mse")
decoder.fit(X_train, y_train)
y_pred = decoder.predict(X_test)
accuracy = (y_pred == y_test).mean() * 100.
print("Graph-net classification accuracy : %g%%" % accuracy)
from nilearn.plotting import plot_stat_map, show
coef_img = decoder.coef_img_
plot_stat_map(coef_img, background_img,
              title="graph-net: accuracy %g%%" % accuracy,
              cut_coords=(-34, -16), display_mode="yz")<jupyter_output><empty_output><jupyter_text>**TV-L1**<jupyter_code>decoder = SpaceNetClassifier(memory="nilearn_cache", memory_level=2, penalty='tv-l1', loss="mse", verbose=1)
decoder.fit(X_train, y_train)
y_pred = decoder.predict(X_test)
accuracy = (y_pred == y_test).mean() * 100.
print("Graph-net classification accuracy : %g%%" % accuracy)
from nilearn.plotting import plot_stat_map, show
coef_img = decoder.coef_img_
plot_stat_map(coef_img, background_img,
              title="graph-net: accuracy %g%%" % accuracy,
              cut_coords=(-34, -16), display_mode="yz")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1292302587923026, 'fraction_numerical': 0.015428002654280026, 'mean_word_length': 5.986095017381229, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '391465', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2124, 'n_tokens_neox': 2018, 'n_words': 551, 'path': '4_2_spacenet.ipynb', 'repo_name': 'dohmatob/nilearn_tutorials', 'size': '6028'}