#!/usr/bin/env Rscript
args<-commandArgs(TRUE); root<-if(length(args))normalizePath(args[1]) else normalizePath('.')
.libPaths(c(file.path(root,'code/Rlib'),.libPaths()))
library(survey);library(survival);library(splines);library(ggplot2)
options(survey.lonely.psu='adjust',stringsAsFactors=FALSE)
out<-Sys.getenv('NHANES_OUTDIR',file.path(root,'results'));dir.create(out,recursive=TRUE,showWarnings=FALSE)
input<-file.path(root,'data/analysis.csv')
d<-read.csv(input)
required<-c('tyg','age','sex','race','SDMVSTRA','SDMVPSU','weight_20yr','PERMTH_EXM','MORTSTAT','cvd_death','bmi','diabetes','hypertension','cvd_history','cancer')
missing<-setdiff(required,names(d));if(length(missing))stop('DATA CONTRACT INCOMPLETE: ',paste(missing,collapse=', '))
if(!nrow(d))stop('No analysis observations')
d$time<-d$PERMTH_EXM;d$allcause<-d$MORTSTAT;d$cvd<-d$cvd_death;d$wt<-d$weight_20yr
if(any(!is.finite(d$wt)|d$wt<=0))stop('Nonpositive/invalid fasting weight')
for(v in c('sex','race','diabetes','hypertension','cvd_history','cancer','education','smoking','alcohol'))if(v%in%names(d))d[[v]]<-factor(d[[v]])
d<-d[complete.cases(d[,c('tyg','time','allcause','cvd','SDMVSTRA','SDMVPSU','wt')])&d$time>0,]
extra<-intersect(c('education','pir','smoking','bmi','diabetes','hypertension','cvd_history','cancer'),names(d))
adj<-c('age','sex','race',extra);rhs<-c('1','age+sex+race',paste(adj,collapse='+'))
base_design<-svydesign(ids=~SDMVPSU,strata=~SDMVSTRA,weights=~wt,data=d,nest=TRUE)
q<-as.numeric(svyquantile(~tyg,base_design,quantiles=c(.25,.5,.75),ci=FALSE)[[1]])
d$quartile<-cut(d$tyg,c(-Inf,q,Inf),labels=paste0('Q',1:4));med<-sapply(1:4,function(j)median(d$tyg[d$quartile==paste0('Q',j)]));d$trend<-med[as.integer(d$quartile)]
des<-function(x){x<-x[is.finite(x$wt)&x$wt>0&!is.na(x$SDMVPSU)&!is.na(x$SDMVSTRA),];svydesign(ids=~SDMVPSU,strata=~SDMVSTRA,weights=~wt,data=droplevels(x),nest=TRUE)}
D<-des(d);write.csv(data.frame(quantile=c(.25,.5,.75),tyg=q),file.path(out,'quartile_cutpoints.csv'),row.names=FALSE)
extract<-function(m,event,model,analysis){b<-coef(m);se<-sqrt(diag(vcov(m)));data.frame(event,model,analysis,term=names(b),logHR=b,SE=se,HR=exp(b),lower=exp(b-1.96*se),upper=exp(b+1.96*se),P=2*pnorm(-abs(b/se)),N=m$n,events=m$nevent,row.names=NULL)}
fit<-function(e,exposure,adj='1',dd=D){aa<-strsplit(adj,'+',fixed=TRUE)[[1]];aa<-aa[aa=='1'|vapply(aa,function(v)if(v=='1')TRUE else length(unique(na.omit(dd$variables[[v]])))>1,logical(1))]; if(!length(aa))aa<-'1'; svycoxph(as.formula(paste0('Surv(time,',e,')~',exposure,'+',paste(aa,collapse='+'))),design=dd)}
save_table<-function(x,n)write.csv(x,file.path(out,paste0(n,'.csv')),row.names=FALSE)
results<-list();ph<-list();curves<-list();tests<-list();turns<-list()
for(e in c('allcause','cvd')){
for(k in 1:3)for(ex in c('tyg','quartile','trend')){m<-fit(e,ex,rhs[k]);results[[length(results)+1]]<-extract(m,e,paste0('M',k),ex)}
cc<-d[complete.cases(d[,adj]),];DD<-des(cc);m<-fit(e,'tyg',rhs[3],DD);ap<-coxph(formula(m),data=cc,weights=wt/mean(wt),x=TRUE);z<-cox.zph(ap);ph[[e]]<-data.frame(event=e,term=rownames(z$table),z$table,row.names=NULL)
png(file.path(out,paste0('ph_',e,'.png')),width=1800,height=1800,res=300);plot(z);dev.off();pdf(file.path(out,paste0('ph_',e,'.pdf')));plot(z);dev.off()
# Natural cubic spline with four fixed weighted quantile knots (5,35,65,95%).
kn<-as.numeric(svyquantile(~tyg,DD,c(.05,.35,.65,.95),ci=FALSE)[[1]]);B<-ns(cc$tyg,knots=kn[2:3],Boundary.knots=kn[c(1,4)]);for(j in 1:3)cc[[paste0('s',j)]]<-B[,j];DS<-des(cc);sm<-fit(e,'s1+s2+s3',rhs[3],DS)
# Nonlinearity by testing departure from a linear spline basis through residualized bases.
lmB<-lm(B~cc$tyg);resB<-residuals(lmB);QR<-qr(resB);keep<-QR$pivot[1:2];for(j in 1:2)cc[[paste0('n',j)]]<-resB[,keep[j]];DN<-des(cc);nm<-fit(e,'tyg+n1+n2',rhs[3],DN)
ov<-regTermTest(sm,~s1+s2+s3,method='Wald');nl<-regTermTest(nm,~n1+n2,method='Wald');tests[[e]]<-data.frame(event=e,overall_P=ov$p,nonlinear_P=nl$p,knot1=kn[1],knot2=kn[2],knot3=kn[3],knot4=kn[4])
grid<-seq(quantile(cc$tyg,.01),quantile(cc$tyg,.99),length.out=150);ref<-as.numeric(svyquantile(~tyg,DD,.5,ci=FALSE)[[1]]);gB<-predict(B,grid);rB<-predict(B,ref);delta<-sweep(gB,2,as.vector(rB));idx<-match(c('s1','s2','s3'),names(coef(sm)));bet<-coef(sm)[idx];vv<-vcov(sm)[idx,idx];eta<-as.vector(delta%*%bet);se<-sqrt(rowSums((delta%*%vv)*delta));rr<-data.frame(tyg=grid,HR=exp(eta),lower=exp(eta-1.96*se),upper=exp(eta+1.96*se),reference=ref);save_table(rr,paste0('rcs_',e));p<-ggplot(rr,aes(tyg,HR))+geom_ribbon(aes(ymin=lower,ymax=upper),alpha=.18)+geom_line()+geom_hline(yintercept=1,lty=2)+theme_bw()+labs(x='TyG',y='Hazard ratio (95% CI)',subtitle=sprintf('Overall P=%.3g; nonlinear P=%.3g',ov$p,nl$p));ggsave(file.path(out,paste0('rcs_',e,'.png')),p,width=6,height=4,dpi=300);ggsave(file.path(out,paste0('rcs_',e,'.pdf')),p,width=6,height=4)
# Exploratory change point selected by pseudo-partial-likelihood; Wald hinge test, not valid naive LRT.
candidates<-seq(quantile(cc$tyg,.15),quantile(cc$tyg,.85),length.out=31);ll<-rep(NA,length(candidates));for(i in seq_along(candidates)){cc$hinge<-pmax(0,cc$tyg-candidates[i]);fm<-fit(e,'tyg+hinge',rhs[3],des(cc));ll[i]<-fm$ll[2]};cpt<-candidates[which.max(ll)];cc$hinge<-pmax(0,cc$tyg-cpt);fm<-fit(e,'tyg+hinge',rhs[3],des(cc));bb<-coef(fm)[c('tyg','hinge')];vv<-vcov(fm)[c('tyg','hinge'),c('tyg','hinge')];turns[[e]]<-data.frame(event=e,breakpoint=cpt,left_HR=exp(bb[1]),left_lower=exp(bb[1]-1.96*sqrt(vv[1,1])),left_upper=exp(bb[1]+1.96*sqrt(vv[1,1])),right_HR=exp(sum(bb)),right_lower=exp(sum(bb)-1.96*sqrt(sum(vv))),right_upper=exp(sum(bb)+1.96*sqrt(sum(vv))),hinge_P=regTermTest(fm,~hinge)$p)
km<-svykm(as.formula(paste0('Surv(time,',e,')~quartile')),D,se=FALSE);for(ext in c('png','pdf')){if(ext=='png')png(file.path(out,paste0('km_',e,'.png')),1800,1400,res=300)else pdf(file.path(out,paste0('km_',e,'.pdf')));plot(km,xlab='Months',ylab='Survival',pars=list(col=c('#0072B2','#D55E00','#009E73','#CC79A7'),lty=1:4));legend('bottomleft',paste0('Q',1:4),col=c('#0072B2','#D55E00','#009E73','#CC79A7'),lty=1:4);dev.off()};capture.output(svylogrank(as.formula(paste0('Surv(time,',e,')~quartile')),D),file=file.path(out,paste0('logrank_',e,'.txt')))
}
save_table(do.call(rbind,results),'cox_results');save_table(do.call(rbind,ph),'ph_test');save_table(do.call(rbind,tests),'rcs_tests');save_table(do.call(rbind,turns),'breakpoints')
# Baseline weighted means/SE and categorical percentages with design tests.
tab<-list();for(v in adj){for(g in levels(d$quartile)){s<-subset(D,quartile==g);if(is.numeric(d[[v]])){a<-svymean(as.formula(paste0('~',v)),s,na.rm=TRUE);tab[[length(tab)+1]]<-data.frame(variable=v,level='mean',group=g,estimate=as.numeric(coef(a)),SE=as.numeric(SE(a)))}else{a<-svymean(as.formula(paste0('~',v)),s,na.rm=TRUE);tab[[length(tab)+1]]<-data.frame(variable=v,level=names(coef(a)),group=g,estimate=as.numeric(coef(a)),SE=as.numeric(SE(a)))}}};save_table(do.call(rbind,tab),'baseline')
# Subgroups: adjust for all covariates excluding constant subgroup term.
d$age_group<-factor(ifelse(d$age<60,'<60','>=60'));d$bmi_group<-cut(d$bmi,c(-Inf,25,30,Inf),right=FALSE);sg<-list();inter<-list()
for(e in c('allcause','cvd'))for(v in c('age_group','sex','race','bmi_group','diabetes','hypertension','cvd_history')){a<-setdiff(adj,v);form<-paste(a,collapse='+');for(g in levels(factor(d[[v]]))){dd<-d[!is.na(d[[v]])&d[[v]]==g,];aa<-a[sapply(dd[,a,drop=FALSE],function(x)length(unique(na.omit(x)))>1)];m<-try(fit(e,'tyg',paste(aa,collapse='+'),des(dd)),silent=TRUE);if(!inherits(m,'try-error'))sg[[length(sg)+1]]<-transform(extract(m,e,v,g),subgroup=g)};m<-try(fit(e,paste0('tyg*',v),form,des(d)),silent=TRUE);if(!inherits(m,'try-error')){bn<-names(coef(m));ii<-grep(':',bn);b<-coef(m)[ii];V<-vcov(m)[ii,ii,drop=FALSE];p<-pchisq(as.numeric(t(b)%*%solve(V,b)),length(ii),lower.tail=FALSE);inter[[length(inter)+1]]<-data.frame(event=e,subgroup=v,P_interaction=p)}}
save_table(do.call(rbind,sg),'subgroup_results');save_table(do.call(rbind,inter),'interaction_tests')
sens<-list();for(e in c('allcause','cvd')){for(s in c('exclude_early_deaths','exclude_CVD_cancer')){dd<-if(s=='exclude_early_deaths')d[!(d$allcause==1&d$time<=24),]else d[as.character(d$cvd_history)=='0'&as.character(d$cancer)=='0',];m<-fit(e,'tyg',rhs[3],des(dd));sens[[length(sens)+1]]<-extract(m,e,'M3',s)};m<-coxph(as.formula(paste0('Surv(time,',e,')~tyg+',rhs[3])),data=d);sens[[length(sens)+1]]<-extract(m,e,'M3','unweighted')};save_table(do.call(rbind,sens),'sensitivity_results');writeLines(capture.output(sessionInfo()),file.path(out,'sessionInfo.txt'))
# Multiple imputation: include outcomes, follow-up and sampling variables as predictors;
# exclude identifiers and exposure/outcome/design from imputation. Pool with Rubin rules.
if(anyNA(d[,adj,drop=FALSE])){
 if(!requireNamespace('mice',quietly=TRUE))stop('mice required for missing-covariate sensitivity')
 mi_data<-d[,unique(c('tyg','time','allcause','cvd','wt','SDMVSTRA','SDMVPSU',adj))];ini<-mice::mice(mi_data,maxit=0,printFlag=FALSE);meth<-ini$method;pred<-ini$predictorMatrix
 fixed<-c('tyg','time','allcause','cvd','wt','SDMVSTRA','SDMVPSU');meth[fixed]<-'';pred[,'SDMVPSU']<-0
 imp<-mice::mice(mi_data,m=5,maxit=10,method=meth,predictorMatrix=pred,seed=20261027,printFlag=FALSE);pooled<-list()
 for(e in c('allcause','cvd')){fits<-lapply(1:5,function(i)fit(e,'tyg',rhs[3],des(mice::complete(imp,i))));pool<-mitools::MIcombine(lapply(fits,coef),lapply(fits,vcov));b<-pool$coefficients;se<-sqrt(diag(pool$variance));df<-pool$df;pooled[[e]]<-data.frame(event=e,model='M3',analysis='multiple_imputation_m5',term=names(b),logHR=b,SE=se,HR=exp(b),lower=exp(b-qt(.975,df)*se),upper=exp(b+qt(.975,df)*se),P=2*pt(-abs(b/se),df));};save_table(do.call(rbind,pooled),'multiple_imputation');saveRDS(imp,file.path(out,'imputations.rds'))
}else save_table(data.frame(status='No missing covariates; imputation equals complete-data analysis'),'multiple_imputation_status')
# Forest plots preserve confidence intervals on a logarithmic HR scale.
sgt<-do.call(rbind,sg);sgt<-sgt[sgt$term=='tyg',];for(e in unique(sgt$event)){z<-sgt[sgt$event==e,];z$label<-paste(z$model,z$analysis,sep=': ');p<-ggplot(z,aes(HR,reorder(label,HR)))+geom_vline(xintercept=1,lty=2)+geom_errorbarh(aes(xmin=lower,xmax=upper),height=.18)+geom_point()+scale_x_log10()+theme_bw()+labs(x='Hazard ratio per 1 TyG (95% CI)',y='Subgroup');ggsave(file.path(out,paste0('forest_',e,'.png')),p,width=7,height=8,dpi=300);ggsave(file.path(out,paste0('forest_',e,'.pdf')),p,width=7,height=8)}
btest<-list();for(v in adj){if(is.numeric(d[[v]])){fm<-svyglm(as.formula(paste0(v,'~quartile')),D);pp<-regTermTest(fm,~quartile)$p}else pp<-svychisq(as.formula(paste0('~',v,'+quartile')),D,statistic='F')$p.value;btest[[v]]<-data.frame(variable=v,P=pp)};save_table(do.call(rbind,btest),'baseline_tests')
# Auxiliary unweighted LRT at the survey-selected change point; ignores breakpoint search.
aux<-list();for(e in c('allcause','cvd')){dd<-d[complete.cases(d[,adj]),];cp<-turns[[e]]$breakpoint;dd$hinge<-pmax(0,dd$tyg-cp);m0<-coxph(as.formula(paste0('Surv(time,',e,')~tyg+',rhs[3])),data=dd);m1<-coxph(as.formula(paste0('Surv(time,',e,')~tyg+hinge+',rhs[3])),data=dd);stat<-2*diff(c(m0$loglik[2],m1$loglik[2]));aux[[e]]<-data.frame(event=e,breakpoint=cp,LRT=stat,df=1,P=pchisq(stat,1,lower.tail=FALSE),limitation='Auxiliary unweighted LRT; does not account for survey design or breakpoint search')};save_table(do.call(rbind,aux),'auxiliary_LRT')

# Freeze corrected PH and expanded change-point analyses after the main pipeline.
for(supp in c('ph_supplement.R','breakpoint_supplement.R')) {status<-system2('Rscript',c(file.path(root,'code',supp),root));if(status!=0)stop('Supplement failed: ',supp)}

fs<-file.path(root,'code/forest_final.py');system2('python',fs)
if(exists('imp')){capture.output(imp$loggedEvents,file=file.path(out,'MI_loggedEvents.txt'));write.csv(as.data.frame.table(imp$chainMean,responseName='chain_mean'),file.path(out,'MI_chainMean.csv'),row.names=FALSE)}
writeLines(capture.output(sessionInfo()),file.path(out,'sessionInfo.txt'))