# Anxiety tests, knockout vs wild-type: eight measures from the elevated plus maze (EPM) and open field (OF), and four
# ways to draw a conclusion about anxiety. Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats
d = pd.read_csv("anxiety.csv")                       # one row per mouse
ANX = ["pct_open_time", "pct_open_entries", "open_entries", "pct_centre_time", "centre_entries"]   # higher = less anxious
ACT = ["closed_entries", "distance_m", "rearing"]                                                     # activity
ko, wt = d[d.genotype == "knockout"], d[d.genotype == "wild-type"]
p = {}
for c in ANX + ACT:
    p[c] = stats.ttest_ind(ko[c], wt[c], equal_var=False).pvalue
    print(f"{c}: wild-type {wt[c].mean():.1f}, knockout {ko[c].mean():.1f}, Welch p = {p[c]:.2g}")
pa = np.array([p[c] for c in ANX])
print(f"A. Any of the {len(ANX)} anxiety measures with p < 0.05: {'yes' if (pa < 0.05).any() else 'no'} ({np.sum(pa < 0.05)} of {len(ANX)})")
o = np.argsort(pa); holm = np.maximum.accumulate(np.minimum(1, pa[o] * (len(pa) - np.arange(len(pa)))))
print(f"B. Holm-adjusted over the {len(ANX)} anxiety measures: smallest adjusted p = {holm.min():.2g}")
t = stats.ttest_ind(ko.pct_open_time, wt.pct_open_time, equal_var=False); ci = t.confidence_interval()
print(f"C. Primary outcome chosen in advance, % time in open arms: difference {ko.pct_open_time.mean() - wt.pct_open_time.mean():.1f} points "
      f"(95% CI {ci.low:.1f} to {ci.high:.1f}), p = {t.pvalue:.2g}")
Z = ["pct_open_time", "pct_open_entries", "pct_centre_time"]                  # composite of the three percentage measures
z = sum((d[c] - wt[c].mean()) / wt[c].std() for c in Z) / len(Z)            # z-scores against the wild-type mean and SD
t = stats.ttest_ind(z[d.genotype == "knockout"], z[d.genotype == "wild-type"], equal_var=False); ci = t.confidence_interval()
print(f"D. Composite anxiety score (mean z of {', '.join(Z)}): difference {z[d.genotype == 'knockout'].mean() - z[d.genotype == 'wild-type'].mean():.2f} SD "
      f"(95% CI {ci.low:.2f} to {ci.high:.2f}), p = {t.pvalue:.2g}")