# 2.4. Searchlight : finding voxels containing information¶

## 2.4.1. Principle of the Searchlight¶

Searchlight was introduced in Information-based functional brain mapping, Nikolaus Kriegeskorte,
Rainer Goebel and Peter Bandettini (PNAS 2006) and consists in scanning the
images volume with a *searchlight*. Briefly, a ball of given radius is
scanned across the brain volume and the prediction accuracy of a
classifier trained on the corresponding voxels is measured.

## 2.4.2. Preparing the data¶

### 2.4.2.1. Loading¶

Fetching the data from internet and loading it can be done with the provided functions (see Loading data):

```
import numpy as np
import nibabel
from nilearn import datasets
from nilearn.image import new_img_like
haxby_dataset = datasets.fetch_haxby_simple()
# print basic information on the dataset
print('Anatomical nifti image (3D) is located at: %s' % haxby_dataset.mask)
print('Functional nifti image (4D) is located at: %s' % haxby_dataset.func)
fmri_filename = haxby_dataset.func
fmri_img = nibabel.load(fmri_filename)
y, session = np.loadtxt(haxby_dataset.session_target).astype('int').T
conditions = np.recfromtxt(haxby_dataset.conditions_target)['f0']
```

### 2.4.2.2. Reshaping the data¶

For this example we need:

- to put X in the form
*n_samples*x*n_features* - compute a mean image for visualization background
- limit our analysis to the face and house conditions (like in the decoding tutorial)

```
from nilearn.image import index_img
condition_mask = np.logical_or(conditions == b'face', conditions == b'house')
fmri_img = index_img(fmri_img, condition_mask)
y, session = y[condition_mask], session[condition_mask]
conditions = conditions[condition_mask]
```

### 2.4.2.3. Masking¶

One of the main elements that distinguish Searchlight from other algorithms is the notion of structuring element that scans the entire volume. If this seems rather intuitive, it has in fact an impact on the masking procedure.

Most of the time, fMRI data is masked and then given to the algorithm. This is not possible in the case of Searchlight because, to compute the score of non-masked voxels, some masked voxels may be needed. This is why two masks will be used here :

*mask_img*is the anatomical mask*process_mask_img*is a subset of mask and contains voxels to be processed.

*process_mask_img* will then be used to restrain computation to one slice, in the
back of the brain. *mask_img* will ensure that no value outside the brain is
taken into account when iterating with the sphere.

```
mask_img = nibabel.load(haxby_dataset.mask)
# .astype() makes a copy.
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)
```

## 2.4.3. Third Step: Setting up the searchlight¶

### 2.4.3.1. Classifier¶

The classifier used by default by `SearchLight`

is LinearSVC with C=1 but
this can be customed easily by passing an estimator parameter to the
cross validation. See scikit-learn documentation for other classifiers.

### 2.4.3.2. Score function¶

Here we use precision as metrics to measure the proportion of true
positives among all positive results for one class. Others metrics can be
specified by the “scoring” argument to the `SearchLight`

, as
detailed in the scikit-learn documentation

### 2.4.3.3. Cross validation¶

`SearchLight`

will iterate on the volume and give a score to each voxel. This
score is computed by running a classifier on selected voxels. In order to make
this score as accurate as possible (and avoid overfitting), a cross validation
is made.

As `SearchLight`

is computationally costly, we have chosen a cross
validation method that does not take too much time.
*K*-Fold along with *K* = 4 is a
good compromise between running time and quality.

```
from sklearn.cross_validation import KFold
cv = KFold(y.size, n_folds=4)
```

## 2.4.4. Running Searchlight¶

Running `SearchLight`

is straightforward now that everything is set.
The only
parameter left is the radius of the ball that will run through the data.
Kriegskorte et al. use a 4mm radius because it yielded the best detection
performance in their simulation.

```
searchlight = nilearn.decoding.SearchLight(
mask_img,
process_mask_img=process_mask_img,
radius=5.6, n_jobs=n_jobs,
verbose=1, cv=cv)
searchlight.fit(fmri_img, y)
```

## 2.4.5. Visualization¶

### 2.4.5.1. Searchlight¶

As the activation map is cropped, we use the mean image of all scans as a background. We can see here that voxels in the visual cortex contains information to distinguish pictures showed to the volunteers, which was the expected result.

```
# Use the fmri mean image as a surrogate of anatomical data
from nilearn import image
from nilearn.plotting import plot_stat_map, show
mean_fmri = image.mean_img(fmri_img)
plot_stat_map(new_img_like(mean_fmri, searchlight.scores_), mean_fmri,
title="Searchlight", display_mode="z", cut_coords=[-16],
colorbar=False)
### F_score results
p_ma = np.ma.array(p_unmasked, mask=np.logical_not(process_mask))
plot_stat_map(new_img_like(mean_fmri, p_ma), mean_fmri,
title="F-scores", display_mode="z",
cut_coords=[-16], colorbar=False)
show()
```

See also

### 2.4.5.2. Comparing to massively univariate analysis: F_score or SPM¶

The standard approach to brain mapping is performed using *Statistical
Parametric Mapping* (SPM), using ANOVA (analysis of variance), and
parametric tests (F-tests ot t-tests).
Here we compute the *p-values* of the voxels [1].
To display the results, we use the negative log of the p-value.

```
p_ma = np.ma.array(p_unmasked, mask=np.logical_not(process_mask))
plot_stat_map(new_img_like(mean_fmri, p_ma), mean_fmri,
title="F-scores", display_mode="z",
cut_coords=[-16], colorbar=False)
show()
```

Parametric scores can be converted into p-values using a reference
theoretical distribution, which is known under specific assumptions
(hence the name *parametric*). In practice, neuroimaging signal has a
complex structure that might not match these assumptions. An exact,
non-parametric *permutation test* can be performed as an alternative
to the parametric test: the residuals of the model are permuted so as
to break any effect and the corresponding decision statistic is
recomputed. One thus builds the distribution of the decision statistic
under the hypothesis that there is no relationship between the tested
variates and the target variates. In neuroimaging, this is generally
done by swapping the signal values of all voxels while the tested
variables remain unchanged [2]. A voxel-wise analysis is then
performed on the permuted data. The relationships between the image
descriptors and the tested variates are broken while the value of the
signal in each particular voxel can be observed with the same
probability than the original value associated to that voxel. Note
that it is hereby assumed that the signal distribution is the same in
every voxel. Several data permutations are performed (typically
10,000) while the scores for every voxel and every data permutation
is stored. The empirical distribution of the scores is thus
constructed (under the hypothesis that there is no relationship
between the tested variates and the neuroimaging signal, the so-called
*null-hypothesis*) and we can compare the original scores to that
distribution: The higher the rank of the original score, the smaller
is its associated p-value. The
`nilearn.mass_univariate.permuted_ols`

function returns the
p-values computed with a permutation test.

```
from nilearn.mass_univariate import permuted_ols
```

```
# We use a two-sided t-test to compute p-values, but we keep trace of the
# effect sign to add it back at the end and thus observe the signed effect
neg_log_pvals, t_scores_original_data, _ = permuted_ols(
grouped_conditions_encoded, grouped_fmri_masked,
# + intercept as a covariate by default
n_perm=10000, two_sided_test=True,
n_jobs=1) # can be changed to use more CPUs
signed_neg_log_pvals = neg_log_pvals * np.sign(t_scores_original_data)
```

The number of tests performed is generally large when full-brain
analysis is performed (> 50,000 voxels). This increases the
probability of finding a significant activation by chance, a
phenomenon that is known to statisticians as the *multiple comparisons
problem*. It is therefore recommended to correct the p-values to take
into account the multiple tests. *Bonferroni correction* consists of
multiplying the p-values by the number of tests (while making sure the
p-values remain smaller than 1). Thus, we control the occurrence of one
false detection *at most*, the so-called *family-wise error control*.
A similar control can be performed when performing a permutation test:
For each permutation, only the maximum value of the F-statistic across
voxels is considered and is used to build the null distribution. It is
crucial to assume that the distribution of the signal is the same in
every voxel so that the F-statistics are comparable. This correction
strategy is applied in nilearn
`nilearn.mass_univariate.permuted_ols`

function.

We observe that the results obtained with a permutation test are less conservative than the ones obtained with a Bonferroni correction strategy.

In nilearn `nilearn.mass_univariate.permuted_ols`

function, we
permute a parametric t-test. Unlike F-test, a t-test can be signed
(*one-sided test*), that is both the absolute value and the sign of an
effect are considered. Thus, only positive effects
can be focused on. It is still possible to perform a two-sided test
equivalent to a permuted F-test by setting the argument
two_sided_test to True. In the example above, we do perform a two-sided
test but add back the sign of the effect at the end using the t-scores obtained
on the original (non-permuted) data. Thus, we can perform two one-sided tests
(a given contrast and its opposite) for the price of one single run.
The example results can be interpreted as follows: viewing faces significantly
activates the Fusiform Face Area as compared to viewing houses, while viewing
houses does not reveals significant supplementary activations as compared to
viewing faces.

[1] | The p-value is the probability of getting the observed values
assuming that nothing happens (i.e. under the null hypothesis).
Therefore, a small p-value indicates that there is a small chance
of getting this data if no real difference existed, so the observed
voxel must be significant. |

[2] | When the variate tested is a scalar (test of the intercept)
–which corresponds to a one sample test–, no swapping can be
performed but one can estimate the null distribution by assuming
symmetry about some reference value. When this value is zero, one can
randomly swap the sign of the target variates (the imaging
signal). nilearn
`nilearn.mass_univariate.permuted_ols` function automatically
adopts the suitable strategy according to the input data. |