Group-label permutation: false positive rate
permute(what="groups") tests group-level differences (e.g. group_comparison()’s null) by shuffling subject labels and recomputing the aggregate effect-size map (cohen(a,b), hedges(a,b), …) from the real, still-spatially-smooth subject data on every draw. This is structurally different from the spatial null methods (moran/cornblath/…) benchmarked in bench01, which instead generate spatially constrained surrogate maps to permute
against a fixed observed map.
Because group-label permutation preserves each subject’s full spatial pattern rather than reshuffling parcel values, there’s a reasonable prior that it is already well calibrated – unlike naive parcel-value shuffling (bench01’s random baseline), it never destroys the spatial autocorrelation (SA) that inflates FPR in the first place. That prior has never been checked empirically. This notebook does that, using the same synthetic GRF (Gaussian Random Field) methodology as bench01 (Markello &
Misic, 2021).
Reference: Markello RD & Misic B (2021). Comparing spatial null models for brain maps. NeuroImage, 236, 118052. https://doi.org/10.1016/j.neuroimage.2021.118052
Reference: Dukart et al. (2021). JuSpace: A tool for spatial correlation analyses of magnetic resonance imaging data with nuclear imaging derived neurotransmitter maps. Human Brain Mapping. https://doi.org/10.1002/hbm.25244
Note: this notebook is compute-intensive. It is designed to be pre-run. For a quick smoke-test, set
N_COHORTS=5, N_X_PER_REP=5, N_PERM=100.
Scope: this notebook stays strictly on false positive rate (type I error). It does not inject true positives / measure power or FDR control across multiple candidate maps – that is deferred to a later, separate benchmark.
Benchmark design
permute(what="groups")’s null-generating unit is a whole subject cohort – one permutation draw relabels subjects and recomputes one aggregate row (e.g. cohen(a,b)) from the entire Y matrix. That’s different from bench01’s per-map null, so the “one big (N_PAIRS x N_PAIRS) matrix, take the diagonal” trick doesn’t directly apply (there’s no way to get N_PAIRS independent cohorts’ worth of output from one permute() call). Instead, for each condition:
Draw ``N_COHORTS`` independent cohorts: fresh GRF-map draws with arbitrary group labels. H0 (no true group difference) holds by construction, since group membership has nothing to do with which GRF map a subject happens to get.
For each cohort, draw ``N_X_PER_REP`` independent GRF maps as candidate X reference maps (disjoint from the cohort’s own draw and from each other).
Run
transform_y()→colocalize("spearman")→permute(what="groups")→get_p_values(), which returns one(1, N_X_PER_REP)row per cohort (multi-row transforms collapse to 1 row because"groups"mode force-poolspooled_p).Pool all
N_COHORTS * N_X_PER_REPp-values per condition;FPR = fraction(p < P_THRESH).
Compute cost is dominated by the n_perm loop, not by how many X columns are tested per call – so N_COHORTS=10 x N_X_PER_REP=50 needs ~5x fewer permute() calls than a N_COHORTS=50 x N_X_PER_REP=10 split, for the same statistical design.
Condition grid
axis |
unpaired |
paired |
|---|---|---|
alpha |
0.0, 1.0, 2.0, 3.0 |
same |
cohort size |
|
|
|
|
|
transform |
|
|
[1]:
# ── Configuration ─────────────────────────────────────────────────────────────
PARC = "Yan200"
ALPHAS = [0.0, 1.0, 2.0, 3.0]
N_COHORTS = 20 # independent cohorts per condition
N_X_PER_REP = 50 # independent GRF X maps tested per cohort (→ pooled n = 500/condition)
N_PERM = 1000 # null permutations per permute() call
P_THRESH = 0.05
SEED = 42
N_PROC = -1
# unpaired conditions: (n_a, n_b) cohort sizes × transforms
UNPAIRED_COHORT_SIZES = [(10, 10), (15, 15), (20, 20), (30, 30), (10, 15), (10, 20), (10, 30)]
UNPAIRED_TRANSFORMS = ["zscore(a,b)", "centile(a,b)", "cohen(a,b)", "hedges(a,b)"]
# paired conditions: n_subj (n_a = n_b = n_subj, forced) × transforms
PAIRED_SUBJ_SIZES = [10, 15, 20, 30]
PAIRED_TRANSFORMS = ["elemdiff(a,b)", "prc(a,b)", "pairedcohen(a,b)"]
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sn
from nispace import NiSpace
from nispace.datasets import fetch_reference
def sample_disjoint(grf_alpha, n_needed, rng):
"""n_needed rows sampled without replacement from grf_alpha."""
n = len(grf_alpha)
assert n >= n_needed, f"Need >={n_needed} maps at this alpha; have {n}"
idx = rng.choice(n, size=n_needed, replace=False)
return grf_alpha.iloc[idx]
[22]:
def run_unpaired(grf_alpha, n_a, n_b, transform, rep_seed):
"""One (alpha, cohort_size, transform) condition -> pooled p-value array."""
p_pool = []
for rep in range(N_COHORTS):
rng = np.random.default_rng(rep_seed + rep)
n_y = n_a + n_b
draw = sample_disjoint(grf_alpha, n_y + N_X_PER_REP, rng)
Y, X = draw.iloc[:n_y], draw.iloc[n_y:]
groups = np.array(["a"] * n_a + ["b"] * n_b)
nsp = NiSpace(x=X, y=Y, parcellation=PARC, verbose=False, n_proc=N_PROC)
nsp.fit()
nsp.transform_y(transform, groups=groups, store=True, verbose=False)
nsp.colocalize("spearman")
nsp.permute(
what="groups", Y_transform=transform, groups=groups,
# next param, group strategy is fixed for now; but explicit so a future
# sweep over "shuffle"/"draw" is a one-line hange, not a silent default-dependency
groups_strategy="shuffle",
n_perm=N_PERM, seed=rep_seed + rep, verbose=False,
)
p_row = nsp.get_p_values().values.flatten()
assert p_row.shape == (N_X_PER_REP,), f"Expected ({N_X_PER_REP},), got {p_row.shape}"
p_pool.extend(p_row.tolist())
return np.array(p_pool)
def run_paired(grf_alpha, n_subj, transform, rep_seed):
"""One (alpha, n_subj, transform) paired condition -> pooled p-value array."""
p_pool = []
for rep in range(N_COHORTS):
rng = np.random.default_rng(rep_seed + rep)
n_y = 2 * n_subj
draw = sample_disjoint(grf_alpha, n_y + N_X_PER_REP, rng)
Y, X = draw.iloc[:n_y], draw.iloc[n_y:]
groups = np.array(["a"] * n_subj + ["b"] * n_subj)
subjects = np.tile(np.arange(n_subj), 2)
nsp = NiSpace(x=X, y=Y, parcellation=PARC, verbose=False, n_proc=N_PROC)
nsp.fit()
nsp.transform_y(transform, groups=groups, subjects=subjects, store=True, verbose=False)
nsp.colocalize("spearman")
nsp.permute(
what="groups", Y_transform=transform, groups=groups, subjects=subjects,
# next param, group strategy is fixed for now; but explicit so a future
# sweep over "shuffle"/"draw" is a one-line hange, not a silent default-dependency
groups_strategy="shuffle",
n_perm=N_PERM, seed=rep_seed + rep, verbose=False,
)
p_row = nsp.get_p_values().values.flatten()
assert p_row.shape == (N_X_PER_REP,), f"Expected ({N_X_PER_REP},), got {p_row.shape}"
p_pool.extend(p_row.tolist())
return np.array(p_pool)
[23]:
n_pooled = N_COHORTS * N_X_PER_REP
ci_lo = P_THRESH - 1.96 * np.sqrt(P_THRESH * (1 - P_THRESH) / n_pooled)
ci_hi = P_THRESH + 1.96 * np.sqrt(P_THRESH * (1 - P_THRESH) / n_pooled)
print(f"95% CI for n_pooled={n_pooled}: [{ci_lo:.3f}, {ci_hi:.3f}]")
grf = fetch_reference("grf", parcellation=PARC, collection="ByAlpha", verbose=False)
records = []
cond_id = 0
t_start = time.time()
for alpha in ALPHAS:
grf_a = grf.loc[f"alpha-{alpha:.1f}"]
print(f"\n{'='*60}\nALPHA: {alpha}")
for (n_a, n_b) in UNPAIRED_COHORT_SIZES:
for transform in UNPAIRED_TRANSFORMS:
cond_id += 1
t0 = time.time()
p_pool = run_unpaired(grf_a, n_a, n_b, transform, SEED + cond_id * 1000)
fpr = (p_pool < P_THRESH).mean()
status = "OK" if ci_lo <= fpr <= ci_hi else "!!"
print(f" {status} unpaired N=({n_a},{n_b}) {transform:16s} FPR={fpr:.3f} "
f"({time.time()-t0:.0f}s, total {time.time()-t_start:.0f}s)")
records.append({
"alpha": alpha, "pairing": "unpaired", "cohort_size": f"({n_a},{n_b})",
"transform": transform, "fpr": fpr, "n_pooled": len(p_pool),
})
for n_subj in PAIRED_SUBJ_SIZES:
for transform in PAIRED_TRANSFORMS:
cond_id += 1
t0 = time.time()
p_pool = run_paired(grf_a, n_subj, transform, SEED + cond_id * 1000)
fpr = (p_pool < P_THRESH).mean()
status = "OK" if ci_lo <= fpr <= ci_hi else "!!"
print(f" {status} paired N={n_subj:<9d} {transform:16s} FPR={fpr:.3f} "
f"({time.time()-t0:.0f}s, total {time.time()-t_start:.0f}s)")
records.append({
"alpha": alpha, "pairing": "paired", "cohort_size": str(n_subj),
"transform": transform, "fpr": fpr, "n_pooled": len(p_pool),
})
results = pd.DataFrame(records)
print(f"\n{'='*60}\n{len(results)} conditions run in {time.time()-t_start:.0f}s total.")
95% CI for n_pooled=1000: [0.036, 0.064]
============================================================
ALPHA: 0.0
OK unpaired N=(10,10) zscore(a,b) FPR=0.051 (54s, total 54s)
OK unpaired N=(10,10) centile(a,b) FPR=0.062 (49s, total 103s)
OK unpaired N=(10,10) cohen(a,b) FPR=0.046 (39s, total 142s)
OK unpaired N=(10,10) hedges(a,b) FPR=0.051 (38s, total 180s)
OK unpaired N=(15,15) zscore(a,b) FPR=0.049 (42s, total 222s)
OK unpaired N=(15,15) centile(a,b) FPR=0.046 (55s, total 277s)
OK unpaired N=(15,15) cohen(a,b) FPR=0.051 (39s, total 316s)
OK unpaired N=(15,15) hedges(a,b) FPR=0.058 (37s, total 353s)
OK unpaired N=(20,20) zscore(a,b) FPR=0.037 (44s, total 397s)
OK unpaired N=(20,20) centile(a,b) FPR=0.046 (61s, total 458s)
OK unpaired N=(20,20) cohen(a,b) FPR=0.050 (37s, total 495s)
OK unpaired N=(20,20) hedges(a,b) FPR=0.038 (37s, total 532s)
OK unpaired N=(30,30) zscore(a,b) FPR=0.054 (52s, total 584s)
OK unpaired N=(30,30) centile(a,b) FPR=0.049 (83s, total 667s)
OK unpaired N=(30,30) cohen(a,b) FPR=0.056 (39s, total 706s)
OK unpaired N=(30,30) hedges(a,b) FPR=0.045 (38s, total 744s)
OK unpaired N=(10,15) zscore(a,b) FPR=0.054 (38s, total 782s)
OK unpaired N=(10,15) centile(a,b) FPR=0.057 (45s, total 828s)
OK unpaired N=(10,15) cohen(a,b) FPR=0.054 (36s, total 864s)
OK unpaired N=(10,15) hedges(a,b) FPR=0.061 (36s, total 900s)
!! unpaired N=(10,20) zscore(a,b) FPR=0.064 (38s, total 938s)
OK unpaired N=(10,20) centile(a,b) FPR=0.041 (47s, total 986s)
OK unpaired N=(10,20) cohen(a,b) FPR=0.054 (37s, total 1022s)
OK unpaired N=(10,20) hedges(a,b) FPR=0.049 (37s, total 1060s)
OK unpaired N=(10,30) zscore(a,b) FPR=0.051 (38s, total 1098s)
OK unpaired N=(10,30) centile(a,b) FPR=0.053 (46s, total 1145s)
OK unpaired N=(10,30) cohen(a,b) FPR=0.037 (37s, total 1182s)
OK unpaired N=(10,30) hedges(a,b) FPR=0.060 (37s, total 1219s)
OK paired N=10 elemdiff(a,b) FPR=0.060 (38s, total 1257s)
!! paired N=10 prc(a,b) FPR=0.036 (39s, total 1296s)
OK paired N=10 pairedcohen(a,b) FPR=0.058 (37s, total 1333s)
OK paired N=15 elemdiff(a,b) FPR=0.050 (41s, total 1374s)
OK paired N=15 prc(a,b) FPR=0.059 (42s, total 1416s)
OK paired N=15 pairedcohen(a,b) FPR=0.038 (36s, total 1452s)
OK paired N=20 elemdiff(a,b) FPR=0.053 (44s, total 1496s)
OK paired N=20 prc(a,b) FPR=0.054 (45s, total 1541s)
OK paired N=20 pairedcohen(a,b) FPR=0.045 (36s, total 1577s)
OK paired N=30 elemdiff(a,b) FPR=0.054 (50s, total 1627s)
OK paired N=30 prc(a,b) FPR=0.053 (51s, total 1678s)
OK paired N=30 pairedcohen(a,b) FPR=0.039 (36s, total 1714s)
============================================================
ALPHA: 1.0
!! unpaired N=(10,10) zscore(a,b) FPR=0.031 (38s, total 1753s)
OK unpaired N=(10,10) centile(a,b) FPR=0.049 (46s, total 1799s)
!! unpaired N=(10,10) cohen(a,b) FPR=0.067 (37s, total 1836s)
OK unpaired N=(10,10) hedges(a,b) FPR=0.042 (37s, total 1873s)
OK unpaired N=(15,15) zscore(a,b) FPR=0.045 (42s, total 1915s)
OK unpaired N=(15,15) centile(a,b) FPR=0.055 (54s, total 1970s)
OK unpaired N=(15,15) cohen(a,b) FPR=0.039 (37s, total 2007s)
OK unpaired N=(15,15) hedges(a,b) FPR=0.045 (37s, total 2044s)
OK unpaired N=(20,20) zscore(a,b) FPR=0.047 (44s, total 2088s)
OK unpaired N=(20,20) centile(a,b) FPR=0.057 (60s, total 2148s)
OK unpaired N=(20,20) cohen(a,b) FPR=0.037 (38s, total 2186s)
OK unpaired N=(20,20) hedges(a,b) FPR=0.049 (37s, total 2222s)
OK unpaired N=(30,30) zscore(a,b) FPR=0.061 (50s, total 2272s)
!! unpaired N=(30,30) centile(a,b) FPR=0.067 (75s, total 2348s)
OK unpaired N=(30,30) cohen(a,b) FPR=0.061 (37s, total 2385s)
OK unpaired N=(30,30) hedges(a,b) FPR=0.047 (39s, total 2424s)
OK unpaired N=(10,15) zscore(a,b) FPR=0.048 (43s, total 2467s)
OK unpaired N=(10,15) centile(a,b) FPR=0.037 (50s, total 2516s)
!! unpaired N=(10,15) cohen(a,b) FPR=0.031 (37s, total 2554s)
OK unpaired N=(10,15) hedges(a,b) FPR=0.052 (37s, total 2591s)
OK unpaired N=(10,20) zscore(a,b) FPR=0.063 (43s, total 2634s)
OK unpaired N=(10,20) centile(a,b) FPR=0.048 (48s, total 2682s)
!! unpaired N=(10,20) cohen(a,b) FPR=0.065 (38s, total 2720s)
OK unpaired N=(10,20) hedges(a,b) FPR=0.037 (37s, total 2757s)
OK unpaired N=(10,30) zscore(a,b) FPR=0.043 (38s, total 2795s)
OK unpaired N=(10,30) centile(a,b) FPR=0.060 (46s, total 2841s)
OK unpaired N=(10,30) cohen(a,b) FPR=0.046 (37s, total 2878s)
OK unpaired N=(10,30) hedges(a,b) FPR=0.048 (37s, total 2915s)
OK paired N=10 elemdiff(a,b) FPR=0.052 (38s, total 2953s)
OK paired N=10 prc(a,b) FPR=0.061 (40s, total 2993s)
OK paired N=10 pairedcohen(a,b) FPR=0.042 (36s, total 3029s)
OK paired N=15 elemdiff(a,b) FPR=0.049 (41s, total 3069s)
OK paired N=15 prc(a,b) FPR=0.052 (42s, total 3111s)
OK paired N=15 pairedcohen(a,b) FPR=0.055 (36s, total 3147s)
!! paired N=20 elemdiff(a,b) FPR=0.035 (44s, total 3191s)
OK paired N=20 prc(a,b) FPR=0.059 (45s, total 3236s)
OK paired N=20 pairedcohen(a,b) FPR=0.049 (36s, total 3272s)
OK paired N=30 elemdiff(a,b) FPR=0.042 (50s, total 3321s)
OK paired N=30 prc(a,b) FPR=0.056 (51s, total 3372s)
OK paired N=30 pairedcohen(a,b) FPR=0.043 (36s, total 3408s)
============================================================
ALPHA: 2.0
!! unpaired N=(10,10) zscore(a,b) FPR=0.093 (38s, total 3447s)
OK unpaired N=(10,10) centile(a,b) FPR=0.057 (46s, total 3492s)
OK unpaired N=(10,10) cohen(a,b) FPR=0.040 (37s, total 3529s)
OK unpaired N=(10,10) hedges(a,b) FPR=0.062 (36s, total 3566s)
OK unpaired N=(15,15) zscore(a,b) FPR=0.057 (41s, total 3606s)
OK unpaired N=(15,15) centile(a,b) FPR=0.055 (52s, total 3659s)
!! unpaired N=(15,15) cohen(a,b) FPR=0.067 (37s, total 3695s)
OK unpaired N=(15,15) hedges(a,b) FPR=0.047 (38s, total 3733s)
OK unpaired N=(20,20) zscore(a,b) FPR=0.047 (44s, total 3777s)
!! unpaired N=(20,20) centile(a,b) FPR=0.068 (65s, total 3843s)
OK unpaired N=(20,20) cohen(a,b) FPR=0.061 (39s, total 3882s)
!! unpaired N=(20,20) hedges(a,b) FPR=0.077 (37s, total 3919s)
OK unpaired N=(30,30) zscore(a,b) FPR=0.040 (50s, total 3969s)
OK unpaired N=(30,30) centile(a,b) FPR=0.039 (75s, total 4043s)
OK unpaired N=(30,30) cohen(a,b) FPR=0.061 (37s, total 4080s)
OK unpaired N=(30,30) hedges(a,b) FPR=0.052 (37s, total 4117s)
OK unpaired N=(10,15) zscore(a,b) FPR=0.057 (38s, total 4155s)
OK unpaired N=(10,15) centile(a,b) FPR=0.039 (49s, total 4205s)
OK unpaired N=(10,15) cohen(a,b) FPR=0.047 (39s, total 4243s)
OK unpaired N=(10,15) hedges(a,b) FPR=0.057 (38s, total 4281s)
!! unpaired N=(10,20) zscore(a,b) FPR=0.036 (39s, total 4321s)
!! unpaired N=(10,20) centile(a,b) FPR=0.029 (49s, total 4370s)
!! unpaired N=(10,20) cohen(a,b) FPR=0.035 (39s, total 4409s)
OK unpaired N=(10,20) hedges(a,b) FPR=0.050 (38s, total 4447s)
OK unpaired N=(10,30) zscore(a,b) FPR=0.043 (39s, total 4486s)
OK unpaired N=(10,30) centile(a,b) FPR=0.042 (47s, total 4533s)
OK unpaired N=(10,30) cohen(a,b) FPR=0.058 (37s, total 4570s)
OK unpaired N=(10,30) hedges(a,b) FPR=0.047 (37s, total 4607s)
OK paired N=10 elemdiff(a,b) FPR=0.038 (39s, total 4646s)
OK paired N=10 prc(a,b) FPR=0.042 (40s, total 4686s)
OK paired N=10 pairedcohen(a,b) FPR=0.051 (36s, total 4722s)
OK paired N=15 elemdiff(a,b) FPR=0.041 (42s, total 4764s)
OK paired N=15 prc(a,b) FPR=0.057 (43s, total 4807s)
OK paired N=15 pairedcohen(a,b) FPR=0.037 (38s, total 4845s)
OK paired N=20 elemdiff(a,b) FPR=0.053 (46s, total 4890s)
OK paired N=20 prc(a,b) FPR=0.055 (53s, total 4943s)
OK paired N=20 pairedcohen(a,b) FPR=0.052 (39s, total 4982s)
!! paired N=30 elemdiff(a,b) FPR=0.067 (62s, total 5044s)
!! paired N=30 prc(a,b) FPR=0.064 (63s, total 5107s)
OK paired N=30 pairedcohen(a,b) FPR=0.060 (42s, total 5149s)
============================================================
ALPHA: 3.0
OK unpaired N=(10,10) zscore(a,b) FPR=0.053 (45s, total 5195s)
OK unpaired N=(10,10) centile(a,b) FPR=0.044 (57s, total 5251s)
OK unpaired N=(10,10) cohen(a,b) FPR=0.048 (42s, total 5293s)
OK unpaired N=(10,10) hedges(a,b) FPR=0.050 (48s, total 5341s)
!! unpaired N=(15,15) zscore(a,b) FPR=0.029 (47s, total 5387s)
!! unpaired N=(15,15) centile(a,b) FPR=0.084 (64s, total 5452s)
OK unpaired N=(15,15) cohen(a,b) FPR=0.046 (42s, total 5493s)
OK unpaired N=(15,15) hedges(a,b) FPR=0.044 (42s, total 5535s)
OK unpaired N=(20,20) zscore(a,b) FPR=0.053 (55s, total 5590s)
!! unpaired N=(20,20) centile(a,b) FPR=0.024 (75s, total 5665s)
OK unpaired N=(20,20) cohen(a,b) FPR=0.045 (42s, total 5708s)
!! unpaired N=(20,20) hedges(a,b) FPR=0.064 (41s, total 5749s)
OK unpaired N=(30,30) zscore(a,b) FPR=0.048 (58s, total 5807s)
OK unpaired N=(30,30) centile(a,b) FPR=0.042 (86s, total 5893s)
OK unpaired N=(30,30) cohen(a,b) FPR=0.056 (42s, total 5935s)
OK unpaired N=(30,30) hedges(a,b) FPR=0.038 (41s, total 5976s)
!! unpaired N=(10,15) zscore(a,b) FPR=0.087 (44s, total 6019s)
!! unpaired N=(10,15) centile(a,b) FPR=0.009 (58s, total 6077s)
OK unpaired N=(10,15) cohen(a,b) FPR=0.054 (43s, total 6120s)
OK unpaired N=(10,15) hedges(a,b) FPR=0.063 (43s, total 6163s)
OK unpaired N=(10,20) zscore(a,b) FPR=0.048 (46s, total 6209s)
!! unpaired N=(10,20) centile(a,b) FPR=0.082 (56s, total 6265s)
OK unpaired N=(10,20) cohen(a,b) FPR=0.063 (42s, total 6307s)
OK unpaired N=(10,20) hedges(a,b) FPR=0.044 (41s, total 6348s)
OK unpaired N=(10,30) zscore(a,b) FPR=0.056 (46s, total 6394s)
OK unpaired N=(10,30) centile(a,b) FPR=0.037 (51s, total 6444s)
OK unpaired N=(10,30) cohen(a,b) FPR=0.055 (39s, total 6484s)
OK unpaired N=(10,30) hedges(a,b) FPR=0.045 (41s, total 6525s)
!! paired N=10 elemdiff(a,b) FPR=0.073 (43s, total 6567s)
!! paired N=10 prc(a,b) FPR=0.030 (45s, total 6612s)
OK paired N=10 pairedcohen(a,b) FPR=0.050 (39s, total 6651s)
OK paired N=15 elemdiff(a,b) FPR=0.047 (45s, total 6697s)
OK paired N=15 prc(a,b) FPR=0.046 (49s, total 6746s)
OK paired N=15 pairedcohen(a,b) FPR=0.038 (39s, total 6785s)
!! paired N=20 elemdiff(a,b) FPR=0.019 (52s, total 6837s)
!! paired N=20 prc(a,b) FPR=0.072 (53s, total 6890s)
OK paired N=20 pairedcohen(a,b) FPR=0.040 (41s, total 6930s)
!! paired N=30 elemdiff(a,b) FPR=0.028 (61s, total 6992s)
OK paired N=30 prc(a,b) FPR=0.060 (61s, total 7052s)
OK paired N=30 pairedcohen(a,b) FPR=0.063 (41s, total 7093s)
============================================================
160 conditions run in 7093s total.
Results
[24]:
# Heatmap: FPR per (transform × alpha), separately for unpaired cohort sizes and paired N_SUBJ
for pairing, cohort_sizes in [("unpaired", UNPAIRED_COHORT_SIZES),
("paired", PAIRED_SUBJ_SIZES)]:
ncol = 4
nrow = np.ceil(len(cohort_sizes) / ncol).astype(int)
fig, axes = plt.subplots(nrow, ncol, figsize=(3*ncol, 3*nrow))
if pairing=="unpaired":
cohort_sizes = [f"({a},{b})" for a, b in UNPAIRED_COHORT_SIZES]
for i, ax in enumerate(axes.ravel()):
if i >= len(cohort_sizes):
ax.set_axis_off()
continue
cohort_size = cohort_sizes[i]
sub = results[(results.pairing == pairing) & (results.cohort_size == str(cohort_size))]
piv = sub.pivot(index="transform", columns="alpha", values="fpr")
sn.heatmap(
piv, ax=ax, annot=True, fmt=".3f", annot_kws={"size": 9}, cmap="RdBu_r", vmin=0, vmax=0.25,
center=P_THRESH, cbar_kws={"label": "FPR"}, linewidths=0.5
)
ax.set_title(f"{pairing} N={cohort_size}", fontsize=10)
ax.set_ylabel("")
if not (i / ncol == i // ncol):
ax.set_yticklabels([])
if i >= (nrow * ncol - ncol):
ax.set_xlabel("GRF alpha")
else:
ax.set_xlabel("")
fig.suptitle(f"False positive rate, permute(what=\"groups\") "
f"(target={P_THRESH}, 95% CI [{ci_lo:.3f}–{ci_hi:.3f}])", fontsize=11)
plt.tight_layout(pad=2.0)
[25]:
# Line plot: FPR vs alpha, averaged across cohort sizes, one line per transform, split by pairing
fig, axes = plt.subplots(1, 2, figsize=(9, 4), sharey=True)
for ax, pairing in zip(axes, ["unpaired", "paired"]):
sub = results[results.pairing == pairing]
for transform, grp in sub.groupby("transform"):
mean_fpr = grp.groupby("alpha")["fpr"].mean()
ax.plot(mean_fpr.index, mean_fpr.values, label=transform, lw=2, marker="o")
ax.axhline(P_THRESH, color="k", lw=1, zorder=0)
ax.axhspan(ci_lo, ci_hi, alpha=0.12, color="green", zorder=0, label="95% CI")
ax.set_xlabel("GRF alpha (spatial autocorrelation)")
ax.set_title(pairing)
ax.set_xticks(ALPHAS)
ax.legend(fontsize=8, loc="upper left")
ax.set_ylim(0,0.25)
axes[0].set_ylabel("False positive rate")
fig.suptitle("FPR vs SA level, averaged across cohort sizes", fontsize=12)
plt.tight_layout(w_pad=2.5)
Interpretation
"shuffle" achieves near-nominal FPR across the entire grid
Across all 160 conditions (4 alphas x 7 unpaired cohort sizes x 4 unpaired transforms, and 4 alphas x 4 paired N_subj x 3 paired transforms):
mean FPR |
std |
min |
max |
|
|---|---|---|---|---|
unpaired (112 conditions) |
0.051 |
0.012 |
0.009 |
0.093 |
paired (48 conditions) |
0.050 |
0.011 |
0.019 |
0.073 |
Both are indistinguishable from the nominal 0.05 target, and critically, flat across every axis tested – no trend by alpha (unpaired 0.049-0.052, paired 0.047-0.051 across alpha=0..3), by cohort size (unpaired 0.049-0.053 across all 7 size combinations; paired 0.047-0.052 across N_subj=10..30), by transform (unpaired 0.049-0.052 across zscore/centile/cohen/hedges; paired 0.048-0.054 across elemdiff/prc/pairedcohen), or by balanced vs. unbalanced unpaired cohorts
(0.051 vs. 0.050). This directly resolves the question this notebook set out to answer: permute(what="groups")’s subject-label-permutation null (with the current "shuffle" default) is well calibrated across the SA range, cohort sizes, and transforms tested here, for both paired and unpaired designs.
Comparison to bench01
Unlike bench01’s spatial null methods – where FPR clearly depends on alpha (e.g. moran ranges from ~0.02 at alpha=1 to ~0.10 at alpha=3, and needs a tuned n_components/K to stay within range – "groups" mode’s calibration here shows no alpha-dependence at all. This makes sense mechanistically: permute(what="groups") never needs a spatial null model in the first place – it permutes subject labels and recomputes the aggregate map from real
subject data each draw, so whatever spatial autocorrelation is present in the data is automatically preserved in every permutation, with no tuning parameter (K, spin procedure, variogram kernel) to get wrong.
Caveat: the displayed 95% CI band is likely too narrow for this design
~18% of individual conditions are flagged !! (outside the plotted CI band), noticeably more than the ~5% you’d expect if each condition’s n_pooled=1000 p-values behaved like 1000 independent Bernoulli(0.05) draws. They don’t quite: the N_X_PER_REP=50 p-values sharing one cohort all derive from the same n_perm=1000 permutation null and the same (arbitrary but fixed) Y group-difference map, so they’re correlated with each other even though each is individually a valid test. The
effective sample size behind each condition’s FPR estimate is therefore somewhere between N_COHORTS=20 and n_pooled=1000, not the full 1000 the plotted CI assumes – so the CI band is an optimistic lower bound on the true sampling uncertainty, and the !! rate above should be read as “somewhat noisier than the band suggests,” not as evidence of miscalibration (the means, which aggregate across many independent cohorts, are the number to trust).