import pandas as pd
from scipy import stats
import statsmodels.formula.api as smf

d = pd.read_csv("cells.csv")            # one row per cell: group, mouse, length

# Route A: one mean per mouse, then Welch's t-test on the means
m = d.groupby(["group", "mouse"], as_index=False)["length"].mean()
print(stats.ttest_ind(m.loc[m.group == "treated", "length"],
                      m.loc[m.group == "control", "length"], equal_var=False))

# Route B: random intercept per mouse (statsmodels reports a z-test)
fit = smf.mixedlm("length ~ group", d, groups=d["mouse"]).fit(reml=True)
print(fit.summary())                   # the p-value column is a z-test: do not report it with few mice

# Small-sample p-value and CI for Route B: compare estimate/SE with t on (mice - groups) df
est, se = fit.params["group[T.treated]"], fit.bse["group[T.treated]"]
df = d["mouse"].nunique() - 2                       # 10 mice, 2 groups -> 8
t = est / se
half = stats.t.ppf(0.975, df) * se
print(f"treated - control = {est:.2f}, 95% CI {est - half:.2f} to {est + half:.2f}, "
      f"t({df}) = {t:.2f}, p = {2 * stats.t.sf(abs(t), df):.3f}")
icc = fit.cov_re.iloc[0, 0] / (fit.cov_re.iloc[0, 0] + fit.scale)
print(f"ICC = {icc:.2f}")