"""Local-only MIDUS prediction study. Emits aggregates, never participant records.""" import os os.environ.setdefault('OPENBLAS_NUM_THREADS','1') from pathlib import Path import json,hashlib import numpy as np import pandas as pd import pyreadstat OUT=Path('results');OUT.mkdir(exist_ok=True) AUDIT={} BASE=['na','na2','h_na','h_na2','stress','h_stress','day','next_weekend'] SLEEP=['h_sleep','sleep_dev','sleep_dev2'] EVENT=['h_events','events_dev'] PA=['pa','h_pa'] NIGHT=['night_dev','night_dev2'] def dump(name,obj): (OUT/name).write_text(json.dumps(obj,indent=2,allow_nan=False)+'\n') def clean(s,lo,hi):return pd.to_numeric(s,errors='coerce').where(lambda a:a.between(lo,hi)) def read_data(): meta=json.loads(Path('data/sources.json').read_text());res={} for wave,prefix in [('M2','B2'),('M3','C2')]: src=next(x for x in meta['files'] if Path(x['path']).name.startswith(wave+'_')) f=Path(src['path']);assert hashlib.sha256(f.read_bytes()).hexdigest()==src['sha256'] raw,_=pyreadstat.read_sav(f,apply_value_formats=False,user_missing=True) raw.columns=raw.columns.str.upper() d=pd.DataFrame({'id':raw.M2ID.astype(int),'family':raw.M2FAMNUM,'day':raw[prefix+'DDAY'].astype(int),'weekday':clean(raw[prefix+'DWEEKD'],1,7)}) for name,items,nreq in [('na',list(range(1,7))+list(range(13,21)),12),('pa',list(range(7,13))+list(range(21,28)),11)]: q=raw[[prefix+'DC'+str(i) for i in items]].apply(lambda s:clean(s,0,4)) d[name]=q.mean(axis=1).where(q.count(axis=1)>=nreq) supplied=prefix+('DNEGAV' if name=='na' else 'DPOSAV') if supplied in raw: m=d[name].notna()&raw[supplied].between(0,4) diff=np.abs(d.loc[m,name]-raw.loc[m,supplied]) full=q.count(axis=1).eq(len(items)) & m assert np.allclose(d.loc[full,name],raw.loc[full,supplied],atol=1e-6),'Fully observed item mapping conflict' bad=(diff>1e-6); nb=int(bad.sum()) AUDIT[wave+'_'+name+'_supplied_score_check']={'compared':int(m.sum()),'full_response_scores_match':True,'partial_response_discrepancies':nb if nb==0 or nb>=10 else 'suppressed_below_10','analysis_uses_preregistered_item_reconstruction':True} h=clean(raw[prefix+'DAW1CH'],0,24);m=clean(raw[prefix+'DAW1CM'],0,60) d['sleep']=(h+m/60).where(lambda x:(x>0)&(x<=24)) for name,items in [('stress',range(1,8)),('events',[8,9,10,12 if wave=='M2' else 13])]: q=raw[[prefix+'DF'+str(i) for i in items]] d[name]=q.eq(1).sum(axis=1).astype(float).where(q.isin([1,2]).all(axis=1)) assert not d.duplicated(['id','day']).any();assert d.day.between(1,8).all() d=d.sort_values(['id','day']).reset_index(drop=True) AUDIT[wave+'_input']={'participants':int(d.id.nunique()),'rows':len(d),'valid_na':int(d.na.notna().sum()),'valid_sleep':int(d.sleep.notna().sum()),'valid_events':int(d.events.notna().sum()),'valid_stress':int(d.stress.notna().sum())} res[wave]=d # Use public family identifiers consistently; conflicts are errors, not silent recoding. both=pd.concat([res['M2'][['id','family']],res['M3'][['id','family']]]) valid=both[both.family.between(1,999998)] assert valid.groupby('id').family.nunique().max()==1 fam=valid.drop_duplicates('id').set_index('id').family.to_dict() for d in res.values(): d['group']=[('F'+str(int(fam[i]))) if i in fam else ('P'+str(i)) for i in d.id] d['fold']=d.group.map(lambda g:int(hashlib.sha256(('sieve-finch-recovery-v1:'+g).encode()).hexdigest(),16)%5) AUDIT['cross_wave_repeated_participants']=len(set(res['M2'].id)&set(res['M3'].id)) return res def construct(d): d=d.copy();g=d.groupby('id',sort=False) for c in ['na','pa','sleep','events','stress']: prior_n=g[c].transform(lambda s:s.notna().cumsum()-s.notna().astype(int)) prior_sum=g[c].transform(lambda s:s.fillna(0).cumsum()-s.fillna(0)) d['n_'+c]=prior_n;d['h_'+c]=prior_sum/prior_n.replace(0,np.nan) for c in ['na','sleep','day','weekday']: d['next_'+c]=g[c].shift(-1) d['na2']=d.na**2;d['h_na2']=d.h_na**2 d['sleep_dev']=d.sleep-d.h_sleep;d['sleep_dev2']=d.sleep_dev**2 d['events_dev']=d.events-d.h_events d['night_dev']=d.next_sleep-d.h_sleep;d['night_dev2']=d.night_dev**2 d['next_weekend']=((d.weekday%7+1).isin([6,7])).astype(float) d['gap']=d.na-d.h_na filters=[('index_days',d.day.between(3,7)),('adjacent_target',d.next_day.eq(d.day+1)&d.next_weekday.eq(d.weekday%7+1)),('observed_current_and_target',d.na.notna()&d.next_na.notna()),('two_prior_na',d.n_na.ge(2)),('stress_and_elevation',d.stress.ge(1)&d.gap.gt(0)),('current_exposures',d[['sleep','events','weekday']].notna().all(axis=1)),('prior_exposures',d[['n_sleep','n_events','n_stress']].ge(2).all(axis=1))] eligible=pd.Series(True,index=d.index);counts={} for name,f in filters: eligible &=f counts[name]={'occasions':int(eligible.sum()),'participants':int(d.loc[eligible,'id'].nunique())} x=d.loc[eligible].copy().reset_index(drop=True) assert np.isfinite(x[BASE+SLEEP+EVENT+['next_na']]).all().all() return x,counts def weights(d): w=1/d.groupby('id').id.transform('size').to_numpy(float);return w*len(w)/w.sum() def fit_predict(train,test,cols): w=weights(train);X=train[cols].to_numpy(float);y=train.next_na.to_numpy(float) mu=np.average(X,axis=0,weights=w);sd=np.sqrt(np.average((X-mu)**2,axis=0,weights=w));sd[sd<1e-12]=1 Z=np.column_stack([np.ones(len(X)),(X-mu)/sd]);pen=np.eye(Z.shape[1]);pen[0,0]=0 coef=np.linalg.solve(Z.T@(w[:,None]*Z)+pen,Z.T@(w*y)) # Independent augmented least-squares solution catches weighting/penalty/intercept mistakes. aug=np.vstack([np.sqrt(w)[:,None]*Z,np.sqrt(pen)]) target=np.r_[np.sqrt(w)*y,np.zeros(Z.shape[1])] independent=np.linalg.lstsq(aug,target,rcond=None)[0] assert np.max(np.abs(coef-independent))<1e-8 T=np.column_stack([np.ones(len(test)),(test[cols].to_numpy(float)-mu)/sd]) pred=np.clip(T@coef,0,4) return pred def predictions(train,test,models): out=pd.DataFrame(index=test.index) out['persistence']=test.na.to_numpy();out['history']=test.h_na.to_numpy() for name in models:out[name]=np.nan checks=[] for fold in range(5): a=train[train.fold!=fold];b=test[test.fold==fold] assert len(a)>0 and len(b)>0 assert set(a.group).isdisjoint(set(b.group));assert set(a.id).isdisjoint(set(b.id)) for name,cols in models.items():out.loc[b.index,name]=fit_predict(a,b,cols) checks.append({'fold':fold,'training_people':int(a.id.nunique()),'test_people':int(b.id.nunique()) if b.id.nunique()>=10 else 'suppressed_below_10','family_overlap':0}) assert np.isfinite(out).all().all() return out,checks def evaluate(d,pred,day_weight=False): y=d.next_na.to_numpy();sq=(pred-y[:,None])**2;ab=np.abs(pred-y[:,None]) # Check target gap algebra independently. for n in pred: a=((pred[n]-d.h_na)-(d.next_na-d.h_na))**2 assert np.allclose(a,sq[n],atol=1e-12) group_info=d[['id','group']].drop_duplicates().set_index('id') person_sq=sq.assign(id=d.id).groupby('id').mean();person_ab=ab.assign(id=d.id).groupby('id').mean() if day_weight: stat=sq.assign(group=d.group).groupby('group').sum() cnt=d.groupby('group').size().reindex(stat.index).to_numpy(float) means=sq.mean();mae=ab.mean() else: stat=person_sq.join(group_info).groupby('group').sum() cnt=group_info.groupby('group').size().reindex(stat.index).to_numpy(float) means=person_sq.mean();mae=person_ab.mean() n=len(stat);rng=np.random.default_rng(20261007) mult=rng.multinomial(n,np.ones(n)/n,size=2000) boot=mult@stat.to_numpy()/(mult@cnt)[:,None] names=list(stat.columns) def contrast(a,b): vals=100*(boot[:,names.index(a)]-boot[:,names.index(b)])/boot[:,names.index(a)] ci=np.quantile(vals,[.025,.975]);point=100*(means[a]-means[b])/means[a] return {'percent_mse_reduction':float(point),'ci95_low':float(ci[0]),'ci95_high':float(ci[1]),'established_positive':bool(ci[0]>0),'above_one_percent_benchmark':bool(ci[0]>1),'upper_excludes_one_percent':bool(ci[1]<1)} result={'participants':int(d.id.nunique()),'families':int(d.group.nunique()),'occasions':len(d),'inference_size_threshold_met':bool(d.id.nunique()>=50 and len(d)>=200),'metrics':{name:{'mse':float(means[name]),'rmse':float(np.sqrt(means[name])),'mae':float(mae[name])} for name in pred}} result['comparisons']={name+'_vs_base':contrast('base',name) for name in pred if name not in ['base','persistence','history']} if 'sleep' in pred and 'events' in pred: vals=100*(boot[:,names.index('sleep')]-boot[:,names.index('events')])/boot[:,names.index('base')] ci=np.quantile(vals,[.025,.975]) result['event_minus_sleep_gain_percentage_points']={'estimate':float(100*(means['sleep']-means['events'])/means['base']),'ci95_low':float(ci[0]),'ci95_high':float(ci[1])} if 'overnight' in pred:result['comparisons']['overnight_vs_joint']=contrast('joint','overnight') gaps=d.assign(current_gap=d.na-d.h_na,next_gap=d.next_na-d.h_na).groupby('id')[['current_gap','next_gap']].mean().mean() result['descriptive_gaps']={k:float(v) for k,v in gaps.items()} folds=[] for f in range(5): m=d.fold==f if d.loc[m,'id'].nunique()<10:folds.append({'fold':f,'suppressed':True});continue vals=sq.loc[m].assign(id=d.loc[m,'id']).groupby('id').mean().mean() folds.append({'fold':f,'people':int(d.loc[m,'id'].nunique()),'occasions':int(m.sum()),'joint_percent_mse_reduction':float(100*(vals['base']-vals['joint'])/vals['base'])}) result['folds']=folds return result def validate_temporal(): # Entirely synthetic, no participant records. d=pd.DataFrame({'id':[1]*8,'family':[1]*8,'day':np.arange(1,9),'weekday':np.arange(8)%7+1,'na':[.1,.2,.8,.4,.6,.2,.8,.3],'pa':[2]*8,'stress':[1]*8,'sleep':[7,8,6,7,8,7,6,8],'events':[0,1,2,1,0,1,1,2],'group':['F1']*8,'fold':[0]*8}) a,_=construct(d);b=d.copy();b.loc[b.day==4,['na','sleep','events']]=[2.0,12,5];c,_=construct(b) x=a[a.day==3];z=c[c.day==3] assert len(x)==len(z)==1 assert np.allclose(x[BASE+SLEEP+EVENT],z[BASE+SLEEP+EVENT]) assert float(x.h_na.iloc[0])==np.mean([.1,.2]) assert x.next_na.iloc[0]!=z.next_na.iloc[0] # A missed day must never turn a two-day gap into a one-day prediction. e,_=construct(d[d.day!=4]);assert 3 not in set(e.day) AUDIT['synthetic_history_future_invariance']=True;AUDIT['synthetic_no_gap_bridging']=True def independent_history_check(raw,eligible): # Independent row loop; no expanding sums or group transforms. lookup={(int(r.id),int(r.day)):r for r in eligible.itertuples()} checked=0 for ident,days in raw.groupby('id'): rows=list(days.sort_values('day').itertuples()) for r in rows: key=(int(ident),int(r.day)) if key not in lookup:continue made=lookup[key] for name in ['na','sleep','events','stress']: values=[float(getattr(z,name)) for z in rows if z.day