Decomposition¶

Decomposition looks for structure with no labels at all: factor a dataset into a small number of components and see what they turn out to be. It is the tool for "what is in this data" rather than "does this data predict that".

BrainData.decompose wraps scikit-learn's decompositions, so the choice of algorithm is one keyword. Here it runs on the pain dataset — 28 subjects at three intensities — and the components turn out to track pain.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from joblib import Memory

from nltools.datasets import fetch_pain
from nltools import concatenate

memory = Memory(".tutorial-cache", verbose=0)

data = memory.cache(fetch_pain)()
data

nltools.data.braindata.BrainData(data=(84, 238955), resolution=2.0mm, space=mni, mask=2mm-MNI152-2009fsl-mask.nii.gz)

Center within subject¶

Every subject contributes three images, and subjects differ from each other far more than intensities differ within a subject. Left alone, the first components would describe who each image belongs to. Subtracting each subject's own mean removes that, leaving the part of each image that is about the intensity.

centered = concatenate(
    [
        data[data.X["SubjectID"] == subject].standardize()
        for subject in data.X["SubjectID"].unique(maintain_order=True)
    ]
)
centered

nltools.data.braindata.BrainData(data=(84, 238955), resolution=2.0mm, space=mni, mask=2mm-MNI152-2009fsl-mask.nii.gz)

Factor the data¶

method selects the algorithm: 'pca', 'ica', 'nnmf', 'fa', 'dictionary' or 'kernelpca'. axis decides which way round the factorization runs. axis='images' treats voxels as observations and images as features, so each component is a brain map and each image gets a loading on it. axis='voxels' does the opposite, which is what you want for resting-state networks over a timeseries.

The result carries decomposition_object (the fitted scikit-learn estimator, with its parameters and explained variance), components (a BrainData of one map per component), and weights (an images-by-components array).

N_COMPONENTS = 5

factors = centered.decompose(method="fa", axis="images", n_components=N_COMPONENTS)

print(sorted(factors))
print(f"components: {factors['components'].shape}")
print(f"weights:    {factors['weights'].shape}")

['components', 'decomposition_object', 'weights']
components: (5, 238955)
weights:    (84, 5)

The weight matrix says how much each image loads on each component, and the component maps say which voxels each one is made of.

weight_figure, weight_axis = plt.subplots(figsize=(5, 8))
sns.heatmap(factors["weights"], ax=weight_axis, center=0, cmap="RdBu_r")
weight_axis.set_xlabel("component")
weight_axis.set_ylabel("image")
weight_axis.set_title("Image loadings")
weight_figure.tight_layout()

# Bound rather than left as the cell's last expression, so the page gets the
# figures and not the list of Figure objects that `plot` returns.
component_figures = factors["components"].plot(limit=N_COMPONENTS)