Permutation Testing
The two examples below both use encoding models evaluated on PPA and FFA voxels, but test different hypotheses and therefore permute different things. See the permutation testing guide for the full step-by-step method.
Example 1
Is encoding performance above chance?
python
import numpy as np
from laion_fmri.subject import load_subject
sub = load_subject("sub-03")
# Load betas for PPA and FFA voxels (check sub.get_available_rois() for exact names)
betas_ppa = sub.get_betas(session="ses-01", roi="PPA") # (n_images, n_ppa)
betas_ffa = sub.get_betas(session="ses-01", roi="FFA") # (n_images, n_ffa)
betas = np.concatenate([betas_ppa, betas_ffa], axis=1) # (n_images, n_voxels)
# model_predictions: your encoding model's predicted responses
# shape (n_images, n_voxels), same image ordering as betas
def pearsonr_cols(A, B):
# Vectorised Pearson r for every column pair; returns shape (n_voxels,)
A = A - A.mean(0); B = B - B.mean(0)
return (A * B).sum(0) / np.sqrt((A**2).sum(0) * (B**2).sum(0))
def mean_r(pred, resp):
return np.nanmean(pearsonr_cols(pred, resp))
# Step 1: Observed test statistic
t_obs = mean_r(model_predictions, betas)
# Steps 2-3: Permute image indices in model predictions -> null distribution
# This breaks the correspondence between which image the model is
# predicting and which brain response it is compared against.
rng = np.random.default_rng(seed=0)
t_perm = np.array([
mean_r(model_predictions[rng.permutation(len(betas))], betas)
for _ in range(10000)
])
# Step 4: One-sided p-value
p = (1 + np.sum(t_perm >= t_obs)) / (1 + len(t_perm))
print(f"Observed mean r = {t_obs:.4f}, p = {p:.4f}")Permutes image indices in the model predictions. This breaks the model's correspondence with any specific image's brain response, generating a one-sided null distribution.
Example 2
Does CLIP outperform AlexNet as an encoding model?
python
import numpy as np
from laion_fmri.subject import load_subject
sub = load_subject("sub-03")
# Load betas for PPA and FFA voxels (check sub.get_available_rois() for exact names)
betas_ppa = sub.get_betas(session="ses-01", roi="PPA") # (n_images, n_ppa)
betas_ffa = sub.get_betas(session="ses-01", roi="FFA") # (n_images, n_ffa)
betas = np.concatenate([betas_ppa, betas_ffa], axis=1) # (n_images, n_voxels)
# clip_predictions: CLIP predictions, shape (n_images, n_voxels)
# alexnet_predictions: AlexNet predictions, shape (n_images, n_voxels)
def pearsonr_cols(A, B):
# Vectorised Pearson r for every column pair; returns shape (n_voxels,)
A = A - A.mean(0); B = B - B.mean(0)
return (A * B).sum(0) / np.sqrt((A**2).sum(0) * (B**2).sum(0))
def mean_r(pred, resp):
return np.nanmean(pearsonr_cols(pred, resp))
# Step 1: Observed test statistic — difference in mean Pearson r
r_clip = mean_r(clip_predictions, betas)
r_alexnet = mean_r(alexnet_predictions, betas)
t_obs = r_clip - r_alexnet
# Steps 2-3: For each permutation, randomly swap CLIP and AlexNet predictions per image,
# then recompute the correlation difference.
# Under H0 (no difference), which model produced which prediction is arbitrary per image.
rng = np.random.default_rng(seed=0)
n = len(betas)
t_perm = np.empty(10000)
for k in range(10000):
swap = rng.integers(0, 2, size=n, dtype=bool)
perm_clip = np.where(swap[:, None], alexnet_predictions, clip_predictions)
perm_alexnet = np.where(swap[:, None], clip_predictions, alexnet_predictions)
t_perm[k] = mean_r(perm_clip, betas) - mean_r(perm_alexnet, betas)
# Step 4: Two-sided p-value
p = (1 + np.sum(np.abs(t_perm) >= np.abs(t_obs))) / (1 + len(t_perm))
print(f"CLIP r = {r_clip:.4f}, AlexNet r = {r_alexnet:.4f}")
print(f"Observed difference = {t_obs:.4f}, p = {p:.4f}")Uses mean Pearson r as the outcome measure. For each permutation, randomly swaps CLIP and AlexNet predictions per image and recomputes the correlation difference — under H₀, which model produced which prediction is arbitrary, so this swap is equally plausible. Two-sided test.