# Dendritic spine density, knockout vs wild-type: segments, neurons or mice as n, and a Poisson model of all segments
# with log(length) as an offset. Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats
d = pd.read_csv("spines.csv")                        # one row per segment: genotype, mouse, neuron, segment, length_um, spines
d["density10"] = 10 * d.spines / d.length_um         # spines per 10 um
nr = d.groupby(["genotype", "mouse", "neuron"])[["spines", "length_um"]].sum().reset_index(); nr["density10"] = 10 * nr.spines / nr.length_um
ms = d.groupby(["genotype", "mouse"])[["spines", "length_um"]].sum().reset_index(); ms["density10"] = 10 * ms.spines / ms.length_um
for g, x in ms.groupby("genotype"):
    print(f"{g}: {len(x)} mice, {x.spines.sum()} spines on {x.length_um.sum():.0f} um; mouse densities {x.density10.min():.1f} to {x.density10.max():.1f} per 10 um")
def welch(df, label):
    k, w = df.density10[df.genotype == "knockout"], df.density10[df.genotype == "wild-type"]; t = stats.ttest_ind(k, w, equal_var=False)
    print(f"{label} ({len(k)} vs {len(w)}), Welch t-test: knockout / wild-type {k.mean() / w.mean():.2f}, p = {t.pvalue:.2g}")
    return t
welch(d, "A. Segments as n"); welch(nr, "B. Neurons as n")
k, w = ms.density10[ms.genotype == "knockout"], ms.density10[ms.genotype == "wild-type"]; t = stats.ttest_ind(k, w, equal_var=False); ci = t.confidence_interval()
print(f"C. Mice as n ({len(k)} vs {len(w)}), Welch t-test on spines per 10 um: difference {k.mean() - w.mean():.2f} (95% CI {ci.low:.2f} to {ci.high:.2f}), p = {t.pvalue:.2g}")
y = d.spines.to_numpy(float); off = np.log(d.length_um.to_numpy()); x = (d.genotype == "knockout").to_numpy(float)
def fit(X):                                          # Poisson regression by Fisher scoring
    b = np.linalg.lstsq(X, np.log(y + 0.5) - off, rcond=None)[0]
    for _ in range(100):
        mu = np.exp(off + X @ b); step = np.linalg.solve(X.T @ (X * mu[:, None]), X.T @ (y - mu)); b = b + step
        if np.max(np.abs(step)) < 1e-10: break
    mu = np.exp(off + X @ b); return b, 2 * np.sum(np.where(y > 0, y * np.log(np.where(y > 0, y, 1) / mu), 0) - (y - mu))
b1, dev1 = fit(np.column_stack([np.ones_like(x), x])); _, dev0 = fit(np.ones((len(y), 1))); lr = dev0 - dev1
print(f"D. Poisson model of all segments, log(length) offset, mice ignored: rate ratio {np.exp(b1[1]):.2f}, likelihood-ratio chi2(1) = {lr:.2f}, p = {stats.chi2.sf(lr, 1):.2g}")