PCA on a single fUSI recording¶

This example shows how to use principal component analysis (PCA) to decompose a fUSI recording into principal axes of variance.

ConfUSIus supports two PCA orientations:

mode="temporal" (default): center across time for each voxel, then find orthogonal voxel-space directions that maximize variance of projected time samples.
mode="spatial": center across voxels for each time point (after transposition), then find orthogonal spatial maps that maximize variance across voxels.

Both return time courses with shape (time, component) via transform, but the statistical objective is applied along different dimensions.

PCA finds an ordered orthogonal set of axes in feature space¹ such that projection onto the first k axes captures the maximum possible variance among all k-dimensional linear projections. Equivalently, these axes are the dominant eigenvectors of the covariance matrix, or of the correlation matrix if variables are properly standardized. If you are interested in linear covariance structure in your fUSI data, PCA is a useful place to start.

Load a fUSI recording¶

To demonstrate the use of PCA on fUSI data, we use a spontaneous activity recording from the Nunez-Elizalde 2022 dataset. See the Datasets user guide for more details on how to download this dataset using ConfUSIus.

from pathlib import Path

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import confusius as cf
from confusius.datasets import fetch_nunez_elizalde_2022
from confusius.decomposition import PCA
from confusius.signal import standardize

# Adapt background color to the current Matplotlib style.
bg_color = mpl.colors.to_hex(mpl.rcParams["figure.facecolor"])

# Keep notebook output compact for large DataArray displays.
xr.set_options(display_expand_data=False)

bids_root = fetch_nunez_elizalde_2022(
    subjects="CR022",
    sessions="20201011",
    tasks="spontaneous",
    acqs="slice03",
)

pwd_path = (
    Path(bids_root)
    / "sub-CR022"
    / "ses-20201011"
    / "fusi"
    / "sub-CR022_ses-20201011_task-spontaneous_acq-slice03_pwd.nii.gz"
)
data = cf.load(pwd_path).compute()
data

Correct for brain motion¶

This recording contains some brain motion, which we can mitigate by performing a rigid translation correction with register_volumewise. This is a common preprocessing step for fUSI data, and it can help avoid spurious components driven by brain motion.

# The `learning_rate` controls the step size of the optimization. A value of 1e-2 is
# a common default that balances convergence speed and stability for typical fUSI data.
data = cf.registration.register_volumewise(data, learning_rate=1e-2)

/home/runner/work/confusius/confusius/.venv/lib/python3.14/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Temporal PCA (`mode="temporal"`)¶

Before performing PCA, we standardize the recording by centering and scaling each voxel's time series to zero mean and unit variance. The standardize function can be used for this purpose. This ensures that PCA captures patterns of correlation rather than patterns of covariance; otherwise, PCA may be dominated by high-variance voxels, such as those near large blood vessels, which may not be of primary interest.

data_std = standardize(data)

In ConfUSIus, the PCA model wraps the familiar scikit-learn PCA model while preserving the fUSI DataArray metadata and coordinates. PCA expects the same arguments as the scikit-learn model, such as n_components for the number of principal components to compute, and random_state for reproducibility (see the API documentation for more details).

By default, PCA uses mode="temporal" (fit on (time, voxels)). A mode="spatial" option is also available, analogous to spatial ICA.

Here, we fit a PCA model with all available components.

pca_t = PCA(random_state=0, mode="temporal")
signals_t = pca_t.fit_transform(data_std)
signals_t

xarray.DataArray

time: 751
component: 751

73.07 -8.361 -4.076 11.58 11.83 ... 0.01418 -0.8985 -0.297 1.711e-05

array([[ 7.3066635e+01, -8.3605938e+00, -4.0763493e+00, ...,
         1.0819744e+00,  1.6663108e+00,  4.4641160e-06],
       [ 5.2037067e+01,  2.2309378e+01, -1.4225268e+01, ...,
         8.1317788e-01, -4.4685724e-01, -5.8623887e-06],
       [ 4.5891502e+01,  2.2932816e+01, -2.1960232e+01, ...,
        -7.2219506e-02, -7.4409738e-02,  8.7668259e-06],
       ...,
       [ 7.4077156e+01, -5.9905930e+00,  4.1795883e+00, ...,
         1.2390769e+00,  8.1436560e-02,  5.4601655e-06],
       [ 7.5606865e+01, -6.2057152e+00,  1.7687122e+01, ...,
        -1.1357114e+00, -4.9661902e-01, -5.1508582e-06],
       [ 5.1808727e+01, -1.8981846e+01,  1.5910517e+01, ...,
        -8.9845270e-01, -2.9696307e-01,  1.7107752e-05]],
      shape=(751, 751), dtype=float32)

Coordinates: (2)
- time
  (time)
  float64
  10.61 10.91 11.21 ... 235.4 235.7
  units :
  s
  volume_acquisition_reference :
  start
  volume_acquisition_duration :
  0.3
```
array([ 10.608,  10.908,  11.208, ..., 235.095, 235.395, 235.695], shape=(751,))
```
- component
  (component)
  int64
  0 1 2 3 4 5 ... 746 747 748 749 750
```
array([  0,   1,   2, ..., 748, 749, 750], shape=(751,))
```
Attributes: (1)
long_name :
PCA signals

Explained variance¶

explained_variance_ratio_ gives the fraction of total variance captured along each selected principal axis. For centered data with shape (time, space), the rank is at most min(space, time - 1), so a final near-zero entry appears in the spectrum; we omit it here for clarity. The scree plot on the left highlights the rapid decay of successive components, while the cumulative curve on the right shows how quickly the total explained variance saturates. These plots can help guide the choice of how many components to retain for further analysis, for example by selecting the "elbow" of the screen plot or a threshold for cumulative variance.

variance_ratio = pca_t.explained_variance_ratio_.isel(component=slice(None, -1))
component_ids = variance_ratio.component.values + 1
cumulative_variance = np.cumsum(variance_ratio.values) * 100

fig, axes = plt.subplots(1, 2, figsize=(10, 3.4), constrained_layout=True)

axes[0].plot(
    component_ids, variance_ratio.values * 100, marker="o", ms=4, color="tab:blue"
)
axes[0].set_yscale("log")
axes[0].set_xlabel("Principal component")
axes[0].set_ylabel("Explained variance (%)")
axes[0].set_title("Scree plot")

axes[1].plot(component_ids, cumulative_variance, marker="o", ms=4, color="tab:blue")
axes[1].set_xlabel("Principal component")
axes[1].set_ylabel("Cumulative explained variance (%)")
_ = axes[1].set_title("Cumulative variance")

Temporal PCA maps and corresponding time courses¶

maps_ stores principal axes in voxel space as a (component, z, y, x) DataArray. transform returned the associated temporal scores in signals_t, a (time, component) DataArray.

Looking at both together is more informative than either alone. The maps show the dominant covariance structure in space, while the scores show how strongly each pattern is expressed over time.

n_show = 6
fig = plt.figure(figsize=(10.5, 8.5), constrained_layout=True)
fig.patch.set_facecolor(bg_color)
gs = fig.add_gridspec(n_show, 2, width_ratios=[1, 3])

axes_tc = [fig.add_subplot(gs[i, 1]) for i in range(n_show)]
for ax in axes_tc[1:]:
    ax.sharex(axes_tc[0])

for i, comp in enumerate(range(n_show)):
    var = float(pca_t.explained_variance_ratio_.sel(component=comp)) * 100

    component_map = pca_t.maps_.isel(component=[comp])
    vmax = float(np.abs(component_map).max())
    cf.plotting.plot_volume(
        component_map,
        axes=fig.add_subplot(gs[i, 0]),
        slice_mode="component",
        cmap="coolwarm",
        vmin=-vmax,
        vmax=vmax,
        show_axes=False,
        show_colorbar=False,
        show_titles=False,
        bg_color=bg_color,
    )

    signals_t.sel(component=comp).plot(ax=axes_tc[i], lw=1.1)
    axes_tc[i].set_title(f"PC {comp + 1} ({var:.1f} %)")
    axes_tc[i].set_xlabel("")
    axes_tc[i].set_ylabel("")

for ax in axes_tc[:-1]:
    ax.tick_params(labelbottom=False)
axes_tc[-1].set_xlabel("Time (s)")
_ = fig.suptitle("Temporal PCA maps and time courses (first 6 components)", fontsize=21)

Spatial PCA (`mode="spatial"`)¶

Spatial PCA transposes the matrix to (voxels, time) before decomposition. In this orientation, variance maximization is performed across voxels, yielding orthogonal spatial maps as principal components.

pca_s = PCA(random_state=0, mode="spatial")
signals_s = pca_s.fit_transform(data_std)
signals_s

xarray.DataArray

time: 751
component: 751

-3.187e+04 -1.548e+03 883.6 -3.861e+03 ... 1.22 39.7 15.21 0.0007041

array([[-3.1866018e+04, -1.5477708e+03,  8.8362970e+02, ...,
        -4.1422390e+01, -7.7635834e+01, -1.7511325e-04],
       [-7.8375811e+03, -7.9947441e+03,  1.9465931e+03, ...,
        -3.7357880e+01,  1.9214939e+01, -5.8000605e-04],
       [-7.6417998e+03, -1.1311187e+04,  3.1012839e+03, ...,
         7.6817799e+00,  5.1758451e+00, -2.5194747e-04],
       ...,
       [-3.3448461e+04,  1.1590634e+03,  5.9401001e+03, ...,
        -5.0094185e+01, -3.7692327e+00, -5.7497656e-04],
       [-3.1455812e+04,  6.7948394e+03,  4.7766562e+03, ...,
         4.7841347e+01,  2.4876919e+01, -2.5275841e-04],
       [-2.8597014e+04,  7.5429858e+03,  2.4295175e+02, ...,
         3.9704731e+01,  1.5206790e+01,  7.0413668e-04]],
      shape=(751, 751), dtype=float32)

Coordinates: (2)
- time
  (time)
  float64
  10.61 10.91 11.21 ... 235.4 235.7
  units :
  s
  volume_acquisition_reference :
  start
  volume_acquisition_duration :
  0.3
```
array([ 10.608,  10.908,  11.208, ..., 235.095, 235.395, 235.695], shape=(751,))
```
- component
  (component)
  int64
  0 1 2 3 4 5 ... 746 747 748 749 750
```
array([  0,   1,   2, ..., 748, 749, 750], shape=(751,))
```
Attributes: (1)
long_name :
PCA signals

As in temporal mode, we inspect maps and corresponding time courses jointly.

n_show = 6
fig = plt.figure(figsize=(10.5, 8.5), constrained_layout=True)
fig.patch.set_facecolor(bg_color)
gs = fig.add_gridspec(n_show, 2, width_ratios=[1, 3])

axes_tc = [fig.add_subplot(gs[i, 1]) for i in range(n_show)]
for ax in axes_tc[1:]:
    ax.sharex(axes_tc[0])

for i, comp in enumerate(range(n_show)):
    var = float(pca_s.explained_variance_ratio_.sel(component=comp)) * 100

    component_map = pca_s.maps_.isel(component=[comp])
    vmax = float(np.abs(component_map).max())
    cf.plotting.plot_volume(
        component_map,
        axes=fig.add_subplot(gs[i, 0]),
        slice_mode="component",
        cmap="coolwarm",
        vmin=-vmax,
        vmax=vmax,
        show_axes=False,
        show_colorbar=False,
        show_titles=False,
        bg_color=bg_color,
    )

    signals_s.sel(component=comp).plot(ax=axes_tc[i], lw=1.1)
    axes_tc[i].set_title(f"PC {comp + 1} ({var:.1f} %)")
    axes_tc[i].set_xlabel("")
    axes_tc[i].set_ylabel("")

for ax in axes_tc[:-1]:
    ax.tick_params(labelbottom=False)
axes_tc[-1].set_xlabel("Time (s)")
_ = fig.suptitle("Spatial PCA maps and time courses (first 6 components)", fontsize=21)

Total running time: 225.7 s

Launch in Binder Download .py Download .ipynb

We usually consider voxels as features and time points as samples, but PCA can be applied in either orientation. ↩