# DASS Analysis Clinic Case 009; 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))

import platform, sys, hashlib
arm=np.repeat([0.,1.],4); y=np.array([1.,2.,3.,4.,3.,4.,5.,10.]); z=np.array([1.,2.,3.,4.,1.,2.,3.,np.nan])
X=np.column_stack((np.ones(8),arm)); full=ols(X,y); keep=~np.isnan(z)
same=ols(X[keep],y[keep]); adjusted=ols(np.column_stack((X[keep],z[keep])),y[keep]); reverse=ols(np.column_stack((np.ones(8),1-arm)),y)
print('All rows treatment-control:',full[0][1]); print('Complete rows treatment-control:',same[0][1]); print('Adjusted complete rows:',adjusted[0][1]); print('Reverse contrast:',reverse[0][1])
print('Included IDs:',np.arange(1,9)); print('Python:',sys.version,'platform:',platform.platform(),'NumPy:',np.__version__)
print('SHA256 of script:',hashlib.sha256(Path(__file__).read_bytes()).hexdigest())
assert np.allclose([full[0][1],same[0][1],adjusted[0][1],reverse[0][1]],[3,1.5,2,-3])
assert np.allclose(X@full[0],np.column_stack((np.ones(8),1-arm))@reverse[0])