<jupyter_start><jupyter_text>Fit a ridge model with wordnet featuresIn this example, we model the fMRI responses with semantic "wordnet" features,manually annotated on each frame of the movie stimulus. The model is aregularized linear regression model, known as ridge regression. Since thismodel is used to predict brain activity from the stimulus, it is called a(voxelwise) encoding model.This example reproduces part of the analysis described in Huth et al (2012)[1]_. See this publication for more details about the experiment, the wordnetfeatures, along with more results and more discussions.*Wordnet features:* The features used in this example are semantic labelsmanually annotated on each frame of the movie stimulus. The semantic labelsinclude nouns (such as "woman", "car", or "building") and verbs (such as"talking", "touching", or "walking"), for a total of 1705 distinct categorylabels. To interpret our model, labels can be organized in a graph of semanticrelashionship based on the `Wordnet `_ dataset.*Summary:* We first concatenate the features with multiple temporal delays toaccount for the slow hemodynamic response. We then use linear regression to fita predictive model of brain activity. The linear regression is regularized toimprove robustness to correlated features and to improve generalizationperformance. The optimal regularization hyperparameter is selected over agrid-search with cross-validation. Finally, the model generalizationperformance is evaluated on a held-out test set, comparing the modelpredictions to the corresponding ground-truth fMRI responses.  Path of the data directory<jupyter_code>import os
from voxelwise_tutorials.io import get_data_home
directory = os.path.join(get_data_home(), "vim-5")
print(directory)
# modify to use another subject
subject = "S01"<jupyter_output><empty_output><jupyter_text>Load the dataWe first load the fMRI responses. These responses have been preprocessed asdecribed in [1]_. The data is separated into a training set ``Y_train`` and atesting set ``Y_test``. The training set is used for fitting models, andselecting the best models and hyperparameters. The test set is later usedto estimate the generalization performance of the selected model. Thetest set contains multiple repetitions of the same experiment to estimatean upper bound of the model prediction accuracy (cf. previous example).<jupyter_code>import numpy as np
from voxelwise_tutorials.io import load_hdf5_array

file_name = os.path.join(directory, "responses", f"{subject}_responses.hdf")
Y_train = load_hdf5_array(file_name, key="Y_train")
Y_test = load_hdf5_array(file_name, key="Y_test")

print("(n_samples_train, n_voxels) =", Y_train.shape)
print("(n_repeats, n_samples_test, n_voxels) =", Y_test.shape)<jupyter_output><empty_output><jupyter_text>If we repeat an experiment multiple times, part of the fMRI responses mightchange. However the modeling features do not change over the repeats, so thevoxelwise encoding model will predict the same signal for each repeat. Tohave an upper bound of the model prediction accuracy, we keep only therepeatable part of the signal by averaging the test repeats.<jupyter_code>Y_test = Y_test.mean(0)

print("(n_samples_test, n_voxels) =", Y_test.shape)<jupyter_output><empty_output><jupyter_text>We fill potential NaN (not-a-number) values with zeros.<jupyter_code>Y_train = np.nan_to_num(Y_train)
Y_test = np.nan_to_num(Y_test)<jupyter_output><empty_output><jupyter_text>Then, we load the semantic "wordnet" features, extracted from the stimulus ateach time point. The features corresponding to the training set are noted``X_train``, and the features corresponding to the test set are noted``X_test``.<jupyter_code>feature_space = "wordnet"

file_name = os.path.join(directory, "features", f"{feature_space}.hdf")
X_train = load_hdf5_array(file_name, key="X_train")
X_test = load_hdf5_array(file_name, key="X_test")

print("(n_samples_train, n_features) =", X_train.shape)
print("(n_samples_test, n_features) =", X_test.shape)<jupyter_output><empty_output><jupyter_text>Define the cross-validation schemeTo select the best hyperparameter through cross-validation, we must define across-validation splitting scheme. Since fMRI time-series are autocorrelatedin time, we should preserve as much as possible the temporal correlation.In other words, since consecutive time samples are correlated, we should notput one time sample in the training set and the immediately following timesample in the validation set. Thus, we define here a leave-one-run-outcross-validation split that keeps each recording run intact.<jupyter_code>from sklearn.model_selection import check_cv
from voxelwise_tutorials.utils import generate_leave_one_run_out

# indice of first sample of each run
run_onsets = load_hdf5_array(file_name, key="run_onsets")
print(run_onsets)<jupyter_output><empty_output><jupyter_text>We define a cross-validation splitter, compatible with ``scikit-learn`` API.<jupyter_code>n_samples_train = X_train.shape[0]
cv = generate_leave_one_run_out(n_samples_train, run_onsets)
cv = check_cv(cv)  # copy the cross-validation splitter into a reusable list<jupyter_output><empty_output><jupyter_text>Define the modelNow, let's define the model pipeline.We first center the features, since we will not use an intercept. The meanvalue in fMRI recording is non-informative, so each run is detrended anddemeaned independently, and we do not need to predict an intercept value inthe linear model.However, we prefer to avoid normalizing by the standard deviation of eachfeature. If the features are extracted in a consistent way from the stimulus,their relative scale is meaningful. Normalizing them independently from eachother would remove this information. Moreover, the wordnet features areone-hot-encoded, which means that each feature is either present (1) or notpresent (0) in each sample. Normalizing one-hot-encoded features is notrecommended, since it would scale disproportionately the infrequent features.<jupyter_code>from sklearn.preprocessing import StandardScaler
scaler = StandardScaler(with_mean=True, with_std=False)<jupyter_output><empty_output><jupyter_text>Then we concatenate the features with multiple delays to account for thehemodynamic response. Due to neurovascular coupling, the recorded BOLD signalis delayed in time with respect to the stimulus onset. With different delayedversions of the features, the linear regression model will weigh each delayedfeature with a different weight to maximize the predictions. With a sampleevery 2 seconds, we typically use 4 delays [1, 2, 3, 4] to cover thehemodynamic response peak. In the next example, we further describe thishemodynamic response estimation.<jupyter_code>from voxelwise_tutorials.delayer import Delayer
delayer = Delayer(delays=[1, 2, 3, 4])<jupyter_output><empty_output><jupyter_text>Finally, we use a ridge regression model. Ridge regression is a linearregression with L2 regularization. The L2 regularization improves robustnessto correlated features and improves generalization performance. However, theL2 regularization is controlled by a hyperparameter ``alpha`` that needs tobe tuned for each dataset. This regularization hyperparameter is usuallyselected over a grid search with cross-validation, selecting thehyperparameter that maximizes the predictive performances on the validationset. More details about cross-validation can be found in the `scikit-learndocumentation`_.For computational reasons, when the number of features is larger than thenumber of samples, it is more efficient to solve ridge regression using the(equivalent) dual formulation [2]_. This dual formulation is equivalent tokernel ridge regression with a linear kernel. Here, we have 3600 trainingsamples, and 1705 * 4 = 6820 features (we multiply by 4 since we use 4 timedelays), therefore it is more efficient to use kernel ridge regression.With one target, we could directly use the pipeline in ``scikit-learn``'s``GridSearchCV``, to select the optimal regularization hyperparameter(``alpha``) over cross-validation. However, ``GridSearchCV`` can onlyoptimize a single score across all voxels (targets). Thus, in themultiple-target case, ``GridSearchCV`` can only optimize (for example) themean score over targets. Here, we want to find a different optimalhyperparameter per target/voxel, so we use the package `himalaya`_ which implements a``scikit-learn`` compatible estimator ``KernelRidgeCV``, with hyperparameterselection independently on each target.<jupyter_code>from himalaya.kernel_ridge import KernelRidgeCV<jupyter_output><empty_output><jupyter_text>``himalaya`` implements different computational backends,including two backends that use GPU for faster computations. The twoavailable GPU backends are "torch_cuda" and "cupy". (Each backend is onlyavailable if you installed the corresponding package with CUDA enabled. Checkthe ``pytorch``/``cupy`` documentation for install instructions.)Here we use the "torch_cuda" backend, but if the import fails we continuewith the default "numpy" backend. The "numpy" backend is expected to beslower since it only uses the CPU.<jupyter_code>from himalaya.backend import set_backend
backend = set_backend("torch_cuda", on_error="warn")
print(backend)<jupyter_output><empty_output><jupyter_text>To speed up model fitting on GPU, we use single precision float numbers.(This step probably does not change significantly the performances on non-GPUbackends.)<jupyter_code>X_train = X_train.astype("float32")
X_test = X_test.astype("float32")<jupyter_output><empty_output><jupyter_text>Since the scale of the regularization hyperparameter ``alpha`` is unknown, weuse a large logarithmic range, and we will check after the fit that besthyperparameters are not all on one range edge.<jupyter_code>alphas = np.logspace(1, 20, 20)<jupyter_output><empty_output><jupyter_text>We also indicate some batch sizes to limit the GPU memory.<jupyter_code>kernel_ridge_cv = KernelRidgeCV(
    alphas=alphas, cv=cv,
    solver_params=dict(n_targets_batch=500, n_alphas_batch=5,
                       n_targets_batch_refit=100))<jupyter_output><empty_output><jupyter_text>Finally, we use a ``scikit-learn`` ``Pipeline`` to link the different stepstogether. A ``Pipeline`` can be used as a regular estimator, calling``pipeline.fit``, ``pipeline.predict``, etc. Using a ``Pipeline`` can beuseful to clarify the different steps, avoid cross-validation mistakes, orautomatically cache intermediate results. See the ``scikit-learn```documentation `_ formore information.<jupyter_code>from sklearn.pipeline import make_pipeline
pipeline = make_pipeline(
    scaler,
    delayer,
    kernel_ridge_cv,
)<jupyter_output><empty_output><jupyter_text>We can display the ``scikit-learn`` pipeline with an HTML diagram.<jupyter_code>from sklearn import set_config
set_config(display='diagram')  # requires scikit-learn 0.23
pipeline<jupyter_output><empty_output><jupyter_text>Fit the modelWe fit on the training set..<jupyter_code>_ = pipeline.fit(X_train, Y_train)<jupyter_output><empty_output><jupyter_text>..and score on the test set. Here the scores are the $R^2$ scores, withvalues in $]-\infty, 1]$. A value of $1$ means the predictionsare perfect.Note that since ``himalaya`` is implementing multiple-targetsmodels, the ``score`` method differs from ``scikit-learn`` API and returnsone score per target/voxel.<jupyter_code>scores = pipeline.score(X_test, Y_test)
print("(n_voxels,) =", scores.shape)<jupyter_output><empty_output><jupyter_text>If we fit the model on GPU, scores are returned on GPU using an array objectspecfic to the backend we used (such as a ``torch.Tensor``). Thus, we need tomove them into ``numpy`` arrays on CPU, to be able to use them for example ina ``matplotlib`` figure.<jupyter_code>scores = backend.to_numpy(scores)<jupyter_output><empty_output><jupyter_text>Plot the model prediction accuracyTo visualize the model prediction accuracy, we can plot it for each voxel ona flattened surface of the brain. To do so, we use a mapper that is specificto the each subject's brain. (Check previous example to see how to use themapper to Freesurfer average surface.)<jupyter_code>import matplotlib.pyplot as plt
from voxelwise_tutorials.viz import plot_flatmap_from_mapper

mapper_file = os.path.join(directory, "mappers", f"{subject}_mappers.hdf")
ax = plot_flatmap_from_mapper(scores, mapper_file, vmin=0, vmax=0.4)
plt.show()<jupyter_output><empty_output><jupyter_text>We can see that the "wordnet" features successfully predict part of themeasured brain activity, with $R^2$ scores as high as 0.4. Note thatthese scores are generalization scores, since they are computed on a test setthat was not used during model fitting. Since we fitted a model independentlyin each voxel, we can inspect the generalization performances at the bestavailable spatial resolution: individual voxels.The best-predicted voxels are located in visual semantic areas like EBA, orFFA. This is expected since the wordnet features encode semantic informationabout the visual stimulus. For more discussions about these results, we referthe reader to the original publication [1]_.  Plot the selected hyperparametersSince the scale of alphas is unknown, we plot the optimal alphas selected bythe solver over cross-validation. This plot is helpful to refine the alphagrid if the range is too small or too large.Note that some voxels might be at the maximum regularization value in thegrid search. These are voxels where the model has no predictive power, thusthe optimal regularization parameter is large to lead to a prediction equalto zero. We do not need to extend the alpha range for these voxels.<jupyter_code>from himalaya.viz import plot_alphas_diagnostic
best_alphas = backend.to_numpy(pipeline[-1].best_alphas_)
plot_alphas_diagnostic(best_alphas=best_alphas, alphas=alphas)
plt.show()<jupyter_output><empty_output><jupyter_text>Visualize the regression coefficientsHere, we go back to the main model on all voxels. Since our model is linear,we can use the (primal) regression coefficients to interpret the model. Thebasic intuition is that the model will use larger coefficients on featuresthat have more predictive power.Since we know the meaning of each feature, we can interpret the largeregression coefficients. In the case of wordnet features, we can even build agraph that represents the features that are linked by a semanticrelationship. We first get the (primal) ridge regression coefficients from the fittedmodel.<jupyter_code>primal_coef = pipeline[-1].get_primal_coef()
primal_coef = backend.to_numpy(primal_coef)
print("(n_delays * n_features, n_voxels) =", primal_coef.shape)<jupyter_output><empty_output><jupyter_text>Here, we are only interested in the voxels with good generalizationperformances. We select an arbitrary threshold of 0.05 (R^2 score).<jupyter_code>primal_coef_selection = primal_coef[:, scores > 0.05]<jupyter_output><empty_output><jupyter_text>Then, we aggregate the coefficients across the different delays.<jupyter_code># split the ridge coefficients per delays
delayer = pipeline.named_steps['delayer']
primal_coef_per_delay = delayer.reshape_by_delays(primal_coef_selection,
                                                  axis=0)
print("(n_delays, n_features, n_voxels) =", primal_coef_per_delay.shape)

# average over delays
average_coef = np.mean(primal_coef_per_delay, axis=0)
print("(n_features, n_voxels) =", average_coef.shape)<jupyter_output><empty_output><jupyter_text>Even after averaging over delays, the coefficient matrix is still too largeto interpret it. Therefore, we use principal component analysis (PCA) toreduce the dimensionality of the matrix.<jupyter_code>from sklearn.decomposition import PCA

pca = PCA(n_components=4)
pca.fit(average_coef.T)
components = pca.components_
print("(n_components, n_features) =", components.shape)<jupyter_output><empty_output><jupyter_text>We can check the ratio of explained variance by each principal component.We see that the first four components already explain a large part of thecoefficients variance.<jupyter_code>print("PCA explained variance =", pca.explained_variance_ratio_)<jupyter_output><empty_output><jupyter_text>Similarly to [1]_, we correct the coefficients of features linked by asemantic relationship. When building the wordnet features, if a frame waslabeled with `wolf`, the authors automatically added the semantically linkedcategories `canine`, `carnivore`, `placental mammal`, `mamma`, `vertebrate`,`chordate`, `organism`, and `whole`. The authors thus argue that the samecorrection needs to be done on the coefficients.<jupyter_code>from voxelwise_tutorials.wordnet import load_wordnet
from voxelwise_tutorials.wordnet import correct_coefficients
_, wordnet_categories = load_wordnet(directory=directory)
components = correct_coefficients(components.T, wordnet_categories).T
components -= components.mean(axis=1)[:, None]
components /= components.std(axis=1)[:, None]<jupyter_output><empty_output><jupyter_text>Finally, we plot the first principal component on the wordnet graph. In suchgraph, edges indicate "is a" relationships (e.g. an `athlete` "is a"`person`). Each marker represents a single noun (circle) or verb (square).The area of each marker indicates the principal component magnitude, and thecolor indicates the sign (red is positive, blue is negative).<jupyter_code>from voxelwise_tutorials.wordnet import plot_wordnet_graph
from voxelwise_tutorials.wordnet import apply_cmap

first_component = components[0]
node_sizes = np.abs(first_component)
node_colors = apply_cmap(first_component, vmin=-2, vmax=2, cmap='coolwarm',
                         n_colors=2)

plot_wordnet_graph(node_colors=node_colors, node_sizes=node_sizes)
plt.show()<jupyter_output><empty_output><jupyter_text>According to the authors of [1]_, "this principal component distinguishesbetween categories with high stimulus energy (e.g. moving objects like`person` and `vehicle`) and those with low stimulus energy (e.g. stationaryobjects like `sky` and `city`)".In this example, because we use only a single subject and we perform adifferent voxel selection, our result is slightly different than in [1]_. Wealso use a different regularization parameter in each voxel, while in [1]_all voxels had the same regularization parameter. We do not aim atreproducing exactly the results in [1]_, but we rather describe the generalapproach. To project the principal component on the cortical surface, we first need touse the fitted PCA to transform the primal weights of all voxels.<jupyter_code># split the ridge coefficients per delays
primal_coef_per_delay = delayer.reshape_by_delays(primal_coef, axis=0)
print("(n_delays, n_features, n_voxels) =", primal_coef_per_delay.shape)
del primal_coef

# average over delays
average_coef = np.mean(primal_coef_per_delay, axis=0)
print("(n_features, n_voxels) =", average_coef.shape)
del primal_coef_per_delay

# transform with the fitted PCA
average_coef_transformed = pca.transform(average_coef.T).T
print("(n_components, n_voxels) =", average_coef_transformed.shape)
del average_coef

# We make sure vmin = -vmax, so that the colormap is centered on 0.
vmax = np.percentile(np.abs(average_coef_transformed), 99.9)

# plot the primal weights projected on the first principal component.
ax = plot_flatmap_from_mapper(average_coef_transformed[0], mapper_file,
                              vmin=-vmax, vmax=vmax, cmap='coolwarm')
plt.show()<jupyter_output><empty_output><jupyter_text>This flatmap shows in which brain regions the model has the largestprojection on the first component. Again, this result is different from theone in [1]_, and should only be considered as reproducing the generalapproach. Following [1]_, we also plot the next three principal components on thewordnet graph, mapping the three vectors to RGB colors.<jupyter_code>from voxelwise_tutorials.wordnet import scale_to_rgb_cube

next_three_components = components[1:4].T
node_sizes = np.linalg.norm(next_three_components, axis=1)
node_colors = scale_to_rgb_cube(next_three_components)
print("(n_nodes, n_channels) =", node_colors.shape)

plot_wordnet_graph(node_colors=node_colors, node_sizes=node_sizes)
plt.show()<jupyter_output><empty_output><jupyter_text>According to the authors of [1]_, "this graph shows that categories thoughtto be semantically related (e.g. athletes and walking) are representedsimilarly in the brain". Finally, we project these principal components on the cortical surface.<jupyter_code>from voxelwise_tutorials.viz import plot_3d_flatmap_from_mapper

voxel_colors = scale_to_rgb_cube(average_coef_transformed[1:4].T, clip=3).T
print("(n_channels, n_voxels) =", voxel_colors.shape)

ax = plot_3d_flatmap_from_mapper(voxel_colors[0], voxel_colors[1],
                                 voxel_colors[2], mapper_file=mapper_file,
                                 vmin=0, vmax=1, vmin2=0, vmax2=1, vmin3=0,
                                 vmax3=1)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07134697762970015, 'fraction_numerical': 0.006425511661113755, 'mean_word_length': 6.308173913043478, 'pattern_counts': {'":': 0, '<': 125, '<?xml version=': 0, '>': 126, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '31', 'dir': 'jupyter-structured-clean-dedup', 'id': '653856', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6173, 'n_tokens_neox': 5840, 'n_words': 2231, 'path': 'tutorials/notebooks/movies_3T/02_plot_wordnet_model.ipynb', 'repo_name': 'gallantlab/voxelwise_tutorials', 'size': '21010'}