# DASS Analysis Clinic Case 001; Python dependencies: numpy, scipy.
from pathlib import Path
import numpy as np
from scipy import stats, optimize
HERE = Path(__file__).resolve().parent

def ols(X, y):
    b = np.linalg.lstsq(X, y, rcond=None)[0]
    residual = y - X @ b
    df = len(y) - X.shape[1]
    cov = (residual @ residual / df) * np.linalg.inv(X.T @ X)
    se = np.sqrt(np.diag(cov))
    ci = np.column_stack((b-stats.t.ppf(.975, df)*se,b+stats.t.ppf(.975,df)*se))
    return b, cov, ci

def power(n, d):
    # Matches R power.t.test(strict=FALSE): rejection tail in effect direction.
    return stats.nct.sf(stats.t.ppf(.975,2*n-2),2*n-2,d*np.sqrt(n/2))

data=np.genfromtxt(HERE/'case-001-shared-data.csv',delimiter=',',names=True,dtype=None,encoding='utf8')
group=np.char.strip(data['group'], '"')
a=data['y'][group=='treatment']; b=data['y'][group=='control']
n=len(a); df=2*n-2; effect=a.mean()-b.mean(); se=np.sqrt((a.var(ddof=1)+b.var(ddof=1))/n)
p=2*stats.t.sf(abs(effect/se),df)
for raw in [.01,.03,.05,.10,.20,.40,.80]:
 z=stats.norm.isf(raw/2); critical=stats.norm.isf(.025)
 print(f'p={raw:.2f}: observed z-power={stats.norm.cdf(z-critical)+stats.norm.cdf(-z-critical):.3f}')
ci=effect+np.array([-1,1])*stats.t.ppf(.975,df)*se
ci90=effect+np.array([-1,1])*stats.t.ppf(.95,df)*se
tost=max(stats.t.sf((effect+6)/se,df),stats.t.cdf((effect-6)/se,df))
sensitivity=optimize.brentq(lambda d:power(n,d)-.8,.01,3)
print(f'Difference={effect:.6f}; p={p:.6f}; 95% CI={ci}; 90% CI={ci90}; TOST p={tost:.6f}; sensitivity d={sensitivity:.6f}')
assert len(a)==len(b)==40 and ci90[0]>-6 and ci90[1]<6