# Myelin g-ratio, knockout vs wild-type: axons as n, a per-axon regression of g-ratio on axon diameter with axons as n,
# the mean g-ratio per mouse, and each mouse's g-ratio at a common axon diameter (1 um) from a line of fibre on axon
# diameter (a line of g-ratio on axon diameter, tried first, did not remove the axon-size effect; see simulation.py). Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats
d = pd.read_csv("axons.csv")                         # one row per myelinated axon: genotype, mouse, axon, axon_um, fibre_um
d["g"] = d.axon_um / d.fibre_um
D0 = 1.0                                             # common axon diameter for the adjusted comparison
def at_d0(x):                                        # per-mouse line of fibre on axon diameter; g-ratio at D0 = D0 / fitted fibre
    b = np.polyfit(x.axon_um, x.fibre_um, 1); return D0 / (b[0] * D0 + b[1])
m = d.groupby(["genotype", "mouse"]).apply(lambda x: pd.Series({"g_mean": x.g.mean(), "g_at_1um": at_d0(x), "diam": x.axon_um.median(), "n": len(x)}), include_groups=False).reset_index()
for g, x in m.groupby("genotype"):
    print(f"{g}: {len(x)} mice, {int(x.n.sum())} axons; mean g-ratio {x.g_mean.mean():.3f}; median axon diameter {x.diam.mean():.2f} um")
k, w = d.g[d.genotype == "knockout"], d.g[d.genotype == "wild-type"]; t = stats.ttest_ind(k, w, equal_var=False)
print(f"A. Axons as n ({len(k)} vs {len(w)}), Welch t-test on g-ratio: difference {k.mean() - w.mean():.3f}, p = {t.pvalue:.2g}")
X = np.column_stack([np.ones(len(d)), d.genotype == "knockout", d.axon_um]).astype(float); y = d.g.to_numpy()
b, *_ = np.linalg.lstsq(X, y, rcond=None); res = y - X @ b; dfr = len(y) - 3
se = np.sqrt(res @ res / dfr * np.linalg.inv(X.T @ X)[1, 1])
print(f"B. Regression of g-ratio on genotype and axon diameter, axons as n: genotype difference {b[1]:.3f}, p = {2 * stats.t.sf(abs(b[1] / se), dfr):.2g}")
for v, lab in (("g_mean", "C. Mice as n, mean g-ratio per mouse"), ("g_at_1um", "D. Mice as n, each mouse's g-ratio at 1 um axon diameter")):
    k, w = m[v][m.genotype == "knockout"], m[v][m.genotype == "wild-type"]; t = stats.ttest_ind(k, w, equal_var=False); ci = t.confidence_interval()
    print(f"{lab} ({len(k)} vs {len(w)}), Welch t-test: difference {k.mean() - w.mean():.3f} (95% CI {ci.low:.3f} to {ci.high:.3f}), p = {t.pvalue:.2g}")