Run permutation tests and correct the full family

from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
from scipy.stats import false_discovery_control
from IPython.display import display as _display
def display(value):
    if isinstance(value,pd.DataFrame):
        table=value.to_html(max_rows=12,escape=True)
        style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
        return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
    return _display(value)
RECIPE_ID='ma06-permutation-tests'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
                      "axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")

Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6
Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads
Core data mode: embedded snapshot; no network requests

def load_snapshot(recipe_id, expected_sha256, snapshot_key):
    """Read embedded data from this notebook, including a renamed copy."""
    for path in sorted(Path.cwd().glob("*.ipynb")):
        try:
            notebook = json.loads(path.read_text(encoding="utf-8"))
            metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
            blob = metadata.get("snapshots", {}).get(snapshot_key)
            if metadata.get("recipe_id") != recipe_id or not blob:
                continue
            if blob.get("sha256") != expected_sha256:
                continue
            raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
            if hashlib.sha256(raw).hexdigest() != expected_sha256:
                raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
            return json.loads(raw)
        except (OSError, json.JSONDecodeError):
            continue
    raise FileNotFoundError(
        "Save/open the complete .ipynb in the active working folder. "
        "The embedded snapshot could not be located for " + recipe_id)

SNAPSHOT_HASHES={'iris': 'cb2512f34911979eef3a42b7a9e6fbc6ca566fd2c0eaa7bd86daefd7ac1cf5ad'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

Embedded snapshots: ['iris']

def audit_numeric_table(frame,columns):
    """Audit observed rows; never silently impute or remove them."""
    if frame.empty or frame.index.has_duplicates:
        raise ValueError('Expected a nonempty table with unique row IDs')
    missing=[name for name in columns if name not in frame.columns]
    if missing:raise ValueError('Missing required columns: '+str(missing))
    values=frame[columns].to_numpy(dtype=float)
    if not np.isfinite(values).all():raise ValueError('Non-finite measurements; inspect missing values explicitly')
    return {'rows':len(frame),'features':len(columns),'finite':True,'duplicate_measurement_rows':int(frame[columns].duplicated().sum())}

print("Method ready: audit_numeric_table")

Method ready: audit_numeric_table

def permutation_mean_test(a,b,repetitions=1999,seed=42):
    """Two-sided independent-group mean permutation test, with the +1 correction."""
    a=np.asarray(a,dtype=float);b=np.asarray(b,dtype=float)
    if a.ndim!=1 or b.ndim!=1 or min(len(a),len(b))<2 or not np.isfinite(np.r_[a,b]).all() or repetitions<1:
        raise ValueError('Need two finite groups and positive permutation count')
    observed=float(a.mean()-b.mean());pool=np.r_[a,b];rng=np.random.default_rng(seed)
    null=np.empty(repetitions)
    for j in range(repetitions):
        shuffled=rng.permutation(pool);null[j]=shuffled[:len(a)].mean()-shuffled[len(a):].mean()
    exceedances=int(np.sum(np.abs(null)>=abs(observed)))
    return dict(difference=observed,p_value=(exceedances+1)/(repetitions+1),null_differences=null,exceedances=exceedances)

print("Method ready: permutation_mean_test")

Method ready: permutation_mean_test

def bh_adjust(pvalues):
    p = np.asarray(pvalues,dtype=float)
    if not len(p):
        return p.copy()
    if not np.isfinite(p).all() or (p<0).any() or (p>1).any():
        raise ValueError("p-values must be finite and in [0,1]")
    order = np.argsort(p,kind="stable")
    adjusted = np.minimum.accumulate((p[order]*len(p)/np.arange(1,len(p)+1))[::-1])[::-1]
    result = np.empty_like(p); result[order] = np.minimum(adjusted,1)
    return result

print("Method ready: bh_adjust")

Method ready: bh_adjust

iris=snapshot("iris")
frame=pd.DataFrame(iris["rows"],columns=iris["columns"])
features=list(frame.columns[:4])
audit=audit_numeric_table(frame,features)
assert len(frame)==150 and frame["species"].nunique()==3
assert (frame[features]>0).all().all()
print("Audit:",audit,"; measurement unit: cm")
print("Source variant:",iris["variant"])
display(frame.head(8))
groups=sorted(frame.species.unique());pairs=[(groups[i],groups[j]) for i in range(3) for j in range(i+1,3)]
records=[];first=None
for feature in features:
 for a,b in pairs:
  result=permutation_mean_test(frame.loc[frame.species.eq(a),feature],frame.loc[frame.species.eq(b),feature],repetitions=1999,seed=42)
  records.append({"feature":feature,"group_a":a,"group_b":b,"mean_difference_cm":result["difference"],"p_value":result["p_value"],"exceedances":result["exceedances"]})
  if first is None:first=result
results=pd.DataFrame(records);assert len(results)==12
display(results.head(8))

Audit: {'rows': 150, 'features': 4, 'finite': True, 'duplicate_measurement_rows': 3} ; measurement unit: cm
Source variant: Original UCI iris.data; no silent correction of the two documented transcription discrepancies

fig,ax=plt.subplots(figsize=(7,4));ax.hist(first["null_differences"],bins=40,edgecolor="white")
ax.axvline(first["difference"],color="red",label="Observed difference")
ax.set(xlabel="Permuted mean difference (cm)",ylabel="Permutations",title="First predeclared comparison under exchangeability");ax.legend();plt.tight_layout();plt.show()

results["BH_q"]=bh_adjust(results.p_value.to_numpy())
results["Bonferroni"]=np.minimum(results.p_value*len(results),1)
labels=[f"{r.feature.replace('_cm','')} {r.group_a.replace('Iris-','')} / {r.group_b.replace('Iris-','')}" for r in results.itertuples()]
fig,ax=plt.subplots(figsize=(9,6));ax.plot(results.p_value,np.arange(12),"o",label="Raw p");ax.plot(results.BH_q,np.arange(12),"x",label="BH q")
ax.axvline(.05,color="gray",linestyle="--");ax.set_yticks(np.arange(12),labels,fontsize=8)
ax.set(xscale="log",xlabel="Probability / adjusted value",ylabel="Predeclared test",title="All comparisons share one adjustment family");ax.legend();plt.tight_layout();plt.show()

effect=results.mean_difference_cm.to_numpy().reshape(4,3)
limit=np.max(np.abs(effect));fig,ax=plt.subplots(figsize=(7,4));im=ax.imshow(effect,cmap="coolwarm",vmin=-limit,vmax=limit)
ax.set_yticks(range(4),features);ax.set_xticks(range(3),[a.replace("Iris-","")+" / "+b.replace("Iris-","") for a,b in pairs],rotation=20,ha="right")
for i in range(4):
 for j in range(3):ax.text(j,i,f"{effect[i,j]:.2f}",ha="center",va="center")
fig.colorbar(im,ax=ax,label="Mean A − mean B (cm)");ax.set_title("Observed measurement differences");plt.tight_layout();plt.show()

assert np.allclose(results.p_value,(results.exceedances+1)/2000)
assert results.p_value.between(1/2000,1).all()
assert np.allclose(results.BH_q,false_discovery_control(results.p_value.to_numpy(),method="bh"))
assert (results.BH_q+1e-15>=results.p_value).all()
OUTPUT_DIR.mkdir(parents=True,exist_ok=True);results.to_csv(OUTPUT_DIR/"all-permutation-tests.csv",index=False)
pd.DataFrame({"null_mean_difference_cm":first["null_differences"]}).to_csv(OUTPUT_DIR/"first-null-distribution.csv",index=False)
print("Joint family:",len(results),"; Monte Carlo p floor:",1/2000,"; independent BH check passed")

Joint family: 12 ; Monte Carlo p floor: 0.0005 ; independent BH check passed

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://archive.ics.uci.edu/dataset/53/iris',timeout=20)
        response.raise_for_status()
        print("Source:",response.url)
        print(response.text[:700])
    except requests.RequestException as error:
        print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
    print("Online extension skipped; embedded core data unchanged")

Online extension skipped; embedded core data unchanged