Quickstart¶

nltools keeps things simple by making use of a few key concepts that can be flexibly combined to perform a wide variety of analyses.

Everything on this page runs in your browser: a real Python interpreter, nothing installed. Press Run to recompute a cell; the first Run downloads Python and nltools, and each cell needs the ones above it to have run first.

Basics¶

Working with neuroimaging data¶

Start by choosing the grid every image lives on. set_brainspace sets it for the session, and 3 mm is what the example dataset below uses: 71,020 voxels against 238,955 at the 2 mm default, small enough to keep this page in a tab.

0.6.0.dev1
BrainSpaceConfig(template='default', resolution=3mm)
  mask: 3mm-MNI152-2009fsl-mask.nii.gz
  brain: 3mm-MNI152-2009fsl-brain.nii.gz
  plot: 3mm-MNI152-2009fsl-T1.nii.gz

import nltools
from nltools import set_brainspace

print(nltools.__version__)
print(set_brainspace(resolution=3))

load_haxby_example simulates a whole experiment and downloads nothing but the MNI template: the eight object conditions of the Haxby task in a randomized block design, each condition driving an 8 mm sphere where that category responds in the real data. It hands back one BrainData per run, images by voxels with one row per TR, and the DesignMatrix that generated it. The fetchers in nltools.datasets load your own data.

nltools.data.braindata.BrainData(data=(72, 71020), resolution=3.0mm, space=mni, mask=3mm-MNI152-2009fsl-mask.nii.gz)
shape: (9, 2)
┌──────────────┬───────┐
│ condition    ┆ count │
│ ---          ┆ ---   │
│ str          ┆ u32   │
╞══════════════╪═══════╡
│ bottle       ┆ 6     │
│ cat          ┆ 6     │
│ chair        ┆ 6     │
│ face         ┆ 6     │
│ house        ┆ 6     │
│ rest         ┆ 24    │
│ scissors     ┆ 6     │
│ scrambledpix ┆ 6     │
│ shoe         ┆ 6     │
└──────────────┴───────┘

from nltools.datasets import load_haxby_example

brains, design_matrices = load_haxby_example()
data, design = brains[0], design_matrices[0]

print(data)
print(data.Y["condition"].value_counts().sort("condition"))

.Y carries the condition of every TR, so indexing pulls one condition out and mean() averages it; smooth first, because six face TRs against one subject's noise make a speckled map. Face blocks minus rest peaks in the right fusiform, and threshold="99.5%" keeps the strongest half percent of voxels — the blue is the patches the other categories drive, which the rest TRs carry a little of. iplot() draws it in an interactive niivue viewer: drag the sliders to rewindow the map, scroll a panel to move through slices.

Interactive viewer (needs a live kernel)

smoothed = data.smooth(fwhm=6)
faces = smoothed[smoothed.Y["condition"] == "face"].mean()
baseline = smoothed[smoothed.Y["condition"] == "rest"].mean()

(faces - baseline).iplot(threshold="99.5%")

More in Working with BrainData, which covers loading, indexing, masks and ROIs, plotting and saving images.

Working with an experimental design¶

A DesignMatrix is a dataframe that knows it describes a timeseries: it carries a sampling frequency, and convolve applies a hemodynamic response to each task regressor, renaming it <column>_c0. The example's design arrives convolved: eight condition regressors and an intercept over 72 TRs:

DesignMatrix(sampling_freq=0.4, shape=(72, 9))
  convolved (8): ['bottle_c0', 'cat_c0', 'chair_c0', 'face_c0', 'house_c0', 'scissors_c0', 'scrambledpix_c0', 'shoe_c0']
  confounds (1): ['.nl_poly_0']

print(design)

design.plot()

More in Working with DesignMatrix, which covers confounds, polynomial terms, multi-run designs and events files.

Working with similarities & distances¶

Adjacency holds a square matrix over a set of nodes: a correlation matrix, a distance matrix, a network. BrainData.distance makes one out of images, and it is worth doing inside a region, since over 71,020 voxels two patterns mostly differ by noise. create_sphere draws regions in MNI millimetres and apply_mask keeps the voxels inside them — here the example's eight response spheres. Correlation distance ignores how strongly a pattern responds overall and compares only its shape across those voxels:

nltools.data.braindata.BrainData(data=(8, 616), resolution=3.0mm, space=mni, mask=None)
nltools.data.adjacency.Adjacency(shape=(8, 8), Y=(0, 0), is_symmetric=True, matrix_type=distance)

from nltools import concatenate
from nltools.data import Adjacency
from nltools.mask import create_sphere

ventral_stream = {
    "face": [40, -50, -20],  # right fusiform face area
    "cat": [-40, -50, -20],  # left fusiform
    "bottle": [46, -78, -6],  # right lateral occipital
    "scissors": [-46, -78, -6],  # left lateral occipital
    "shoe": [36, -66, -16],  # right posterior fusiform
    "chair": [-36, -66, -16],  # left posterior fusiform
    "house": [-26, -44, -10],  # left parahippocampal place area
    "scrambledpix": [0, -88, 2],  # early visual cortex
}
ventral_temporal = create_sphere(
    list(ventral_stream.values()), radius=8, mask=data.mask
)
patterns = concatenate(
    [data[data.Y["condition"] == name].mean() for name in ventral_stream]
).apply_mask(ventral_temporal)
neural = Adjacency(
    patterns.distance(metric="correlation"), labels=list(ventral_stream)
)

print(patterns)
print(neural)

neural.plot()

Eight conditions make 28 distinct pairs, and that is all an Adjacency stores; the plot fills the diagonal back in with the zero a pattern has with itself. The dark cells are the within-category pairs — face with cat, the four man-made objects with each other — while houses and scrambled pictures have no close partner.

More in Working with Adjacency, which covers thresholds, Fisher z, stacking subjects, regression and graphs.

Common analysis workflows¶

Those three objects are all you need. Every permutation and bootstrap test below uses 200 resamples instead of the default 5,000, so each cell takes a second.

Mapping neural responses¶

Fitting a GLM is fit(model="glm", X=design), and compute_contrasts asks it a question — here, how much more each voxel responds to faces than to houses:

data.fit(model="glm", X=design)

data.compute_contrasts("face_c0 - house_c0").plot(
    title="faces - houses, one subject"
)

One subject is not a result: the map a paper reports is a test across subjects. n_runs=5 returns five fresh draws of the experiment, each with its own block order and noise, to stand in for five subjects; stack their contrast maps with concatenate, ttest gives the voxelwise one-sample test, and threshold zeroes every voxel whose p-value misses a cutoff:

from nltools.algorithms import threshold

subjects, subject_designs = load_haxby_example(n_runs=5)
group = concatenate(
    [
        subject.fit(model="glm", X=subject_design).compute_contrasts(
            "face_c0 - house_c0"
        )
        for subject, subject_design in zip(subjects, subject_designs)
    ]
)
group_t = group.ttest()

threshold(group_t["t"], group_t["p"], thr=0.001).plot(title="group t, p < 0.001")

ttest returns the mean, t, z and p maps as separate images.

Predicting neural responses¶

An encoding model turns the same equation around: features go in, and the model is judged on data it never saw. The features here are every convolved condition regressor at twelve delays — 96 columns over 72 TRs, so least squares has no unique solution and the ridge penalty is what makes the fit possible at all. Two runs, one to train on and one held out, both z-scored because ridge fits no intercept and the data sit on a baseline of 100. ridge_cv picks the penalty by cross-validation, one per voxel, and the fit lands on .model: a weight map per feature in betas, the chosen penalty in alpha, and the variance explained on the training run in r2:

nltools.data.braindata.BrainData(data=(96, 71020), resolution=3.0mm, space=mni, mask=3mm-MNI152-2009fsl-mask.nii.gz)
penalty per voxel: [   1.   10.  100. 1000.]
best training r2: 0.82

import numpy as np

def delayed(design):
    """The eight convolved regressors at twelve delays, side by side."""
    names = [name for name in design.columns if name.endswith("_c0")]
    regressors = np.column_stack([np.asarray(design[name]) for name in names])
    shifted = [np.pad(regressors, ((lag, 0), (0, 0))) for lag in range(12)]
    return np.column_stack([block[: len(regressors)] for block in shifted])

runs, run_designs = load_haxby_example(n_runs=2)
train = runs[0].standardize(method="zscore")
held_out = runs[1].standardize(method="zscore")
train.fit(
    model="ridge",
    X=delayed(run_designs[0]),
    ridge_alpha=[1, 10, 100, 1000],
    ridge_cv=5,
    random_state=0,
)

print(train.model.betas)
print(f"penalty per voxel: {np.unique(train.model.alpha.data)}")
print(f"best training r2: {train.model.r2.data.max():.2f}")

Applying those weights to the held-out run's features predicts its timecourse. Correlating that prediction with what the run actually did, voxel by voxel, gives a performance map: not how much a voxel responds, but how well the model accounts for it:

from nltools.data import BrainData

def voxel_correlation(observed, predicted):
    """Correlate two arrays over the same voxels, one voxel at a time."""
    observed = (observed - observed.mean(axis=0)) / observed.std(axis=0)
    predicted = (predicted - predicted.mean(axis=0)) / predicted.std(axis=0)
    return (observed * predicted).mean(axis=0)

encoding_scores = voxel_correlation(
    held_out.data, train.predict(X=delayed(run_designs[1])).data
)

BrainData(encoding_scores[None, :], mask=held_out.mask).plot(
    title="Held-out run, r per voxel"
)

Analyzing neural patterns¶

Representational similarity analysis compares geometries instead of voxels: how far apart the conditions are in the brain, against how far apart a model says they should be. The neural side is the distance matrix from the basics section, and the model side puts two conditions far apart when their categories differ. Both sides are an Adjacency, and similarity correlates them, Spearman by default, with a permutation test that shuffles rows and columns together:

rho = 0.75, p = 0.010 (200 permutations)

categories = np.array(["animate"] * 2 + ["object"] * 4 + ["scene", "control"])
model_rdm = Adjacency(
    (categories[:, None] != categories[None, :]).astype(float),
    matrix_type="distance",
    labels=neural.labels,
)
rsa = neural.similarity(model_rdm, n_permute=200, random_state=0, n_jobs=1)

print(f"rho = {rsa['correlation']:.2f}, p = {rsa['p']:.3f} (200 permutations)")

neural.plot_stacked(
    model_rdm, upper_title="Correlation distance", lower_title="Category model"
)

plot_stacked puts the two matrices in one square, the measured geometry above the diagonal and the model below it.

Decoding asks the same region the opposite question: given the pattern, which condition was it? predict cross-validates a classifier over images and reports the accuracy per fold. Twelve face and house TRs, three folds:

accuracy 0.92, per fold [0.75 1.   1.  ]

decoded = (
    data[data.Y["condition"].is_in(["face", "house"])]
    .apply_mask(ventral_temporal)
    .predict(y="condition", estimator="linear_svc", cv=3)
)

print(f"accuracy {decoded.mean_score:.2f}, per fold {decoded.scores}")