Code part 2
# Define the key correlations we want to test
corr_pairs = [
('Width_proxy', 'IPR', 'width_proxy', 'ipr'),
('Width_proxy', 'Localization', 'width_proxy', 'localization'),
('λ_max', 'IPR', 'lambda_max', 'ipr'),
('Width_proxy', 'λ_max', 'width_proxy', 'lambda_max'),
]
r_obs = {}
for name1, name2, key1, key2 in corr_pairs:
r, _ = pearsonr(metrics_obs[key1], metrics_obs[key2])
r_obs[(name1, name2)] = r
print(f"Observed r({name1}, {name2}) = {r: .4f}")
# === Surrogate testing ===
n_surrogates = 300 # good balance between power and runtime (increase if needed)
print(f"\nRunning {n_surrogates} phase-randomized surrogates...")
surrogate_rs = {(name1, name2): [] for name1, name2, key1, key2 in corr_pairs}
for s in tqdm(range(n_surrogates), desc="Surrogates"):
# Phase randomize each channel independently
data_surr = np.array([phase_randomize(ch) for ch in data])
# Recompute full pipeline on surrogate
fc_surr = compute_dynamic_plv(data_surr, sfreq, verbose=False)
metrics_surr = compute_topc_metrics(fc_surr)
for name1, name2, key1, key2 in corr_pairs:
r_s, _ = pearsonr(metrics_surr[key1], metrics_surr[key2])
surrogate_rs[(name1, name2)].append(r_s)
# === Empirical p-values summary ===
print("\n" "="*70)
print("SURROGATE TEST RESULTS (Phase randomization null)")
print("="*70)
results_table = []
for name1, name2, key1, key2 in corr_pairs:
r_o = r_obs[(name1, name2)]
r_surr = np.array(surrogate_rs[(name1, name2)])
# Two-sided empirical p-value
count_extreme = np.sum(np.abs(r_surr) >= np.abs(r_o))
p_emp = (count_extreme 1) / (n_surrogates 1)
mean_surr = np.mean(r_surr)
std_surr = np.std(r_surr)
sig = "***" if p_emp < 0.001 else "**" if p_emp < 0.01 else "*" if p_emp < 0.05 else "ns"
print(f"\n{name1} ↔ {name2}")
print(f" Observed r = {r_o: .4f}")
print(f" Surrogate mean = {mean_surr: .4f} ± {std_surr:.4f}")
print(f" Empirical p = {p_emp:.4f} {sig}")
results_table.append({
'pair': f"{name1}–{name2}",
'r_obs': r_o,
'r_surr_mean': mean_surr,
'r_surr_std': std_surr,
'p_emp': p_emp,
'sig': sig
})
# Save results
np.save(os.path.join(out_dir, 'surrogate_r_distributions.npy'),
{str(k): np.array(v) for k,v in surrogate_rs.items()})
print(f"\nSaved surrogate correlation distributions to surrogate_r_distributions.npy")
# === Visualization: histograms ===
fig, axes = plt.subplots(2, 2, figsize=(12, 9))
axes = axes.flatten()
for idx, (name1, name2, key1, key2) in enumerate(corr_pairs):
ax = axes[idx]
r_surr = np.array(surrogate_rs[(name1, name2)])
r_o = r_obs[(name1, name2)]
ax.hist(r_surr, bins=40, color='
#5DADE2', edgecolor='white', alpha=0.85, density=True)
ax.axvline(r_o, color='
#E74C3C', linewidth=2.5, label=f'Observed r = {r_o: .3f}')
ax.axvline(np.mean(r_surr), color='gray', linestyle='--', linewidth=1.5,
label=f'Surr. mean = {np.mean(r_surr): .3f}')
ax.set_title(f'{name1} ↔ {name2}\nEmpirical p = {results_table[idx]["p_emp"]:.4f}', fontsize=11)
ax.set_xlabel('Pearson r (surrogates)', fontsize=10)
ax.set_ylabel('Density', fontsize=10)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.suptitle('Surrogate