# Independent model calculations using the author's declared measure construction.
import sys,json
import numpy as np
import pandas as pd
from scipy import stats
sys.path.insert(0,'code')
import nhanes as N
a=N.build(N.load('data/nhanes'));a['phq']=a[N.PHQ].sum(axis=1);a['height']=a.BMXHT
out={}
for col in ['struct','allostatic','bmi','height']:
 d=a.dropna(subset=[col,'RIDAGEYR','male']);y=d[col].to_numpy();sd=y.std();x=d.phq.to_numpy();X=np.column_stack([np.ones(len(d)),d.RIDAGEYR,d.RIDAGEYR**2,d.male,pd.get_dummies(d.cycle,drop_first=True,dtype=float),x]);Q,_=np.linalg.qr(X);ry=y-Q@(Q.T@y)
 rows=[]
 for c in range(8,23):
  h=np.maximum(0,x-c);rh=h-Q@(Q.T@h);b=rh@ry/(rh@rh);res=ry-b*rh;sse=res@res;se=np.sqrt(sse/(len(y)-X.shape[1]-1)/(rh@rh));t=b/se
  rows.append({'c':c,'change_sd':b/sd,'t':t,'pointwise_ci95_sd':[(b-1.96*se)/sd,(b+1.96*se)/sd],'sse':sse,'n_above':int((x>c).sum())})
 out[col]={'n':len(d),'max_signed_t':max(rows,key=lambda r:r['t']),'min_sse':min(rows,key=lambda r:r['sse']),'all_grid':rows}
h=a.dropna(subset=['struct','DPQ100']);h=h[h.DPQ100.isin([0,1,2,3])]
def pc(d):
 C=np.column_stack([np.ones(len(d))]+[stats.rankdata(d[c]) for c in ['cognitive','RIDAGEYR','male']]);Q,_=np.linalg.qr(C)
 x=stats.rankdata(d.struct);y=stats.rankdata(d.somatic);return float(np.corrcoef(x-Q@(Q.T@x),y-Q@(Q.T@y))[0,1])
out['strata']={'high':pc(h[h.DPQ100>=2]),'low':pc(h[h.DPQ100<=1])}
json.dump(out,open('diagnostic.json','w'),indent=2)
print(json.dumps({k:{kk:vv for kk,vv in v.items() if kk!='all_grid'} for k,v in out.items()},indent=2))