Permutation Testing
The example below uses a CLIP encoding model evaluated on PPA and FFA voxels. 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
Is encoding performance higher in PPA than FFA?
python
import numpy as np
from laion_fmri.subject import load_subject
sub = load_subject("sub-03")
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)
# model_preds_ppa: predictions for PPA voxels, shape (n_images, n_ppa)
# model_preds_ffa: predictions for FFA voxels, shape (n_images, n_ffa)
def zscore(X):
return (X - X.mean(0)) / X.std(0)
# Per-image encoding score: mean prediction-response product across voxels
score_ppa = (zscore(model_preds_ppa) * zscore(betas_ppa)).mean(axis=1) # (n_images,)
score_ffa = (zscore(model_preds_ffa) * zscore(betas_ffa)).mean(axis=1) # (n_images,)
# Imagewise PPA-minus-FFA advantage
d = score_ppa - score_ffa # (n_images,)
# Step 1: Observed test statistic
t_obs = d.mean()
# Steps 2-3: Randomly flip sign of each image's difference -> null distribution
# Under the null, positive and negative differences are equally likely.
rng = np.random.default_rng(seed=0)
t_perm = np.array([
(rng.choice([-1, 1], size=len(d)) * d).mean()
for _ in range(10000)
])
# Step 4: Two-sided p-value
p = (1 + np.sum(np.abs(t_perm) >= np.abs(t_obs))) / (1 + len(t_perm))
print(f"Observed mean imagewise PPA advantage = {t_obs:.4f}, p = {p:.4f}")Computes a per-image PPA-minus-FFA encoding advantage, then randomly flips the sign of each image's difference. Under the null, positive and negative differences are equally likely. Two-sided test.