/*********************************************************************************/ /* */ /* BAND FOR DIFFERENCES OF ADJUSTED SURVIVALS BASED ON STRATIFIED COX MODEL */ /* */ /* %SURVBAND aims at providing a confidence band for differences between */ /* adjusted survival of any two treatments. Direct adjustment method was */ /* implemented to account for unbalanced distribution of confounding */ /* variables in treatment groups. */ /* */ /* */ /* Macro parameters: */ /* inputdata - the input sas data name; */ /* time - the survival time variable; */ /* event - the event indicator; */ /* group - the treatment group variable, */ /* which must take values 1,...,K for K groups. */ /* covlist - a list of covariate names; */ /* timecov - the covariate that is modelled as time-dependent variables */ /* with piecewise constant HR's based on the time cutpoint */ /* provided in the next parameter; */ /* cutpoint - the time cutpoint so that the covariate at prior position */ /* has piecewise constant HR's; */ /* starttime - the start time point of the interval for simulation; */ /* stoptime - the stop time point of the interval for simulation; */ /* alpha - input alpha value to get a 100(1-alpha)% confidence interval; */ /* nsimu - number of realizations that will be generated;elected; */ /* seed - a large positive integer used as seed of pseudorandom number */ /* generator algorithm; */ /* outdata - the output sas data name. */ /* */ /* */ /*********************************************************************************/ %macro TIMEBAND(inputdata, time, event, group, covlist,timecov,cutpoint, starttime, stoptime, alpha, nsimu, seed, outdata); proc means data=&inputdata noprint; var &group; output out=maxout max(&group)=numgroup; run; data _null_; set maxout; call symput('numgroup', numgroup); run; proc iml; use &inputdata; read all var {&covlist} into x; close &inputdata; numcov=ncol(x); create ncovout from numcov[colname='numcov']; append from numcov; close ncovout; quit; data _null_; set ncovout; call symput('numcov', numcov); run; proc sort data=&inputdata (where=(&event=1)) out=timedata; by descending &time; run; data _null_; set timedata; if _n_=1 then call symput('lastt',time); run; %put ****** &starttime **** &stoptime ****; %if %sysevalf(&stoptime>&lastt) %then %do; %let stoptime=&lastt; %put The provided stop time was greater than the last event time. It was set to &lastt; %end; %let firstt=0; %do gp=1 %to &numgroup; proc sort data=&inputdata (where=(&event=1 & &group=&gp)) out=timedata; by &time; run; data _null_; set timedata; if _n_=1 then call symput('gptime',time); run; %if %sysevalf(&firstt<&gptime) %then %let firstt=&gptime; %end; %if %sysevalf(&starttime<&firstt) %then %do; %let starttime=&firstt; %put The provided start time was smaller than the maximum of earlest event times across groups. It was set to &firstt; %end; %put provided t &starttime &stoptime; %put data t &firstt &lastt; %if %sysevalf(&starttime>&stoptime) %then %do; %let starttime=&firstt; %let stoptime=&lastt; %put The provided times points are not valid.; %put The start time point was set to &firstt; %put The stop time point was set to &lastt; %end; /***************************************/ /* */ /* Model: a stratified Cox model */ /* */ /***************************************/ proc iml; use &inputdata; read all var {&time} into time; read all var {&event} into event; read all var {&group} into group; read all var {&covlist} into x; read all var {&timecov} into tx; close &inputdata; out=time||event||group||x||tx; names={'time' 'event' 'strata' %do i=1 %to &numcov; "z&i" %end; 'tz'}; create indata from out[colname=names]; append from out; close indata; quit; /**********************************************/ /* */ /* Get regression coefficient estimates from */ /* proc phreg, read in these estimates and */ /* calculate s0(b,t), s1(b,t). */ /* */ /**********************************************/ proc sort data=indata; by descending time descending event; run; data alltime; set indata (keep=time event); by descending time; if first.time; if event; drop event; run; proc sort data=alltime; by time; run; proc sort data=indata; by time; run; data event; merge indata (where=(event=1)) alltime (in=inall); by time; if inall; retain index 0; if first.time then index=index+1; run; proc phreg data=indata covout outest=best noprint; model time*event(0)=%do i=1 %to &numcov; z&i %end; timez1 timez2; timez1=0; timez2=0; if time<=&cutpoint then timez1=tz; if time>&cutpoint then timez2=tz; strata strata; run; data best sigma; set best; if _type_='PARMS' then output best; if _type_='COV' then output sigma; run; %let znum1=%eval(&numcov+1); %let znum2=%eval(&numcov+2); %do group=1 %to &numgroup; data gpone; set indata; where strata=&group; z&znum1=0; z&znum2=0; if time<=&cutpoint then z&znum1=tz; if time>&cutpoint then z&znum2=tz; run; proc iml; use gpone; read all var {time} into time; read all var {%do i=1 %to %eval(&numcov+2); z&i %end;} into zmat; close gpone; use best; read all var{%do i=1 %to &numcov; z&i %end; timez1 timez2} into b; read all var{%do i=1 %to &numcov; z&i %end; timez1 timez1} into b2; close best; numobs=nrow(time); zb=j(numobs,1,0); zb2=j(numobs,1,0); do i=1 to numobs; zb[i]=zmat[i,]*t(b); zb2[i]=zmat[i,]*t(b2); end; outmat=time||zmat||zb||zb2; names={'time' %do i=1 %to %eval(&numcov+2); "z&i" %end; 'zb' 'zb2'}; create gpone from outmat[colname=names]; append from outmat; close gpone; quit; proc sort; by descending time; run; data riskset&group; set gpone (keep=time %do i=1 %to %eval(&numcov+2); z&i %end; zb zb2); by descending time; s0t1+exp(zb); s0t2+exp(zb2); if time>&cutpoint then s0=s0t1; else s0=s0t2; %do i=1 %to &numcov; s1t1_&i+z&i * exp(zb); s1t2_&i+z&i * exp(zb2); if time>&cutpoint then s1_&i=s1t1_&i; else s1_&i=s1t2_&i; %end; %let vn1=%eval(&numcov+1); %let vn2=%eval(&numcov+2); temp1+z&vn2*exp(zb); if time>&cutpoint then temp2+z&vn2*exp(zb2); else temp2+z&vn1*exp(zb2); s1_&vn2=0; if time>&cutpoint then s1_&vn2=temp1; s1_&vn1=0; if time<=&cutpoint then s1_&vn1=temp2; if last.time; keep time s0 %do i=1 %to %eval(&numcov+2); s1_&i %end;; run; proc sort data=riskset&group; by time; run; ODS LISTING CLOSE; proc freq data=indata; where strata=&group & event=1; table time/out=dcount&group; run; ODS LISTING; data riskset&group; merge alltime (in=inalltime) riskset&group dcount&group (keep=time count); by time; if inalltime; if count=. then count=0; retain ts0 %do i=1 %to %eval(&numcov+2); ts1_&i %end; 999; if s0^=. then ts0=s0; else s0=ts0; %do i=1 %to %eval(&numcov+2); if s1_&i^=. then ts1_&i=s1_&i; else s1_&i=ts1_&i; %end; drop ts:; run; %end; data riskset; merge %do i=1 %to &numgroup; riskset&i (rename=(s0=s0&i count=count&i %do j=1 %to %eval(&numcov+2); s1_&j=s1&i._&j %end;)) %end;; run; proc iml; use riskset; read all var{time} into time; %do i=1 %to &numgroup; read all var{s0&i} into s0&i; read all var{%do j=1 %to %eval(&numcov+2); s1&i._&j %end;} into s1&i; read all var{count&i} into count&i; %end; close riskset; use best; read all var{%do i=1 %to &numcov; z&i %end; timez1 timez2} into b; close best; use sigma; read all var{%do i=1 %to &numcov; z&i %end; timez1 timez2} into sigma; close sigma; use indata; read all var{%do i=1 %to &numcov; z&i %end;} into zmat; read all var{tz} into tz; close indata; use event; read all var{time} into timeall; read all var{strata} into strata; read all var{index} into index; read all var{%do i=1 %to &numcov; z&i %end;} into zzall; read all var{tz} into tzall; close event; rinit=ranuni(&seed); numeobs=nrow(strata); numtime=nrow(time); numobs=nrow(zmat); numcov=ncol(s11); out=time; timeindex=j(numtime,1,0); tidx=numtime; do i=numeobs to 1 by -1; if index[i]=tidx then do; timeindex[tidx]=i; tidx=tidx-1; end; end; /************************************************************/ /* */ /* Find direct adjusted survival of K treatments and SE's */ /* */ /************************************************************/ do i=1 to numtime until (time[i]>&cutpoint); end; cutidx=i-1; do i=1 to numeobs until (timeall[i]>&cutpoint); end; ecutidx=i-1; %do trt=1 %to &numgroup; cumuhaz&trt=j(numtime,1,0); w&trt=j(numtime,1,0); wFc&trt=j(numtime,1,0); ctemp=0; wtemp=0; do i=1 to cutidx; ctemp=ctemp+count&trt[i]/s0&trt[i]; wtemp=wtemp+count&trt[i]/s0&trt[i]/s0&trt[i]; cumuhaz&trt[i]=ctemp; w&trt[i]=wtemp; end; wfctemp=0; do i=cutidx+1 to numtime; ctemp=ctemp+count&trt[i]/s0&trt[i]; wtemp=wtemp+count&trt[i]/s0&trt[i]/s0&trt[i]; wfctemp=wfctemp+count&trt[i]/s0&trt[i]/s0&trt[i]; cumuhaz&trt[i]=ctemp; w&trt[i]=wtemp; wfc&trt[i]=wfctemp; end; f&trt=j(numtime,1,0); fexpbz&trt._1=j(numtime,1,0); fexpbz&trt._2=j(numtime,1,0); fh&trt=j(numtime,numcov,0); do i=1 to numobs; ztemp1=zmat[i,]||tz[i]||0; ztemp2=zmat[i,]||0||tz[i]; expbz1=exp(ztemp1*t(b)); expbz2=exp(ztemp2*t(b)); surv=j(numtime,1,0); do j=1 to cutidx; surv[j]=exp(-cumuhaz&trt[j])##expbz1; end; do j=cutidx+1 to numtime; temp=cumuhaz&trt[j]#expbz2+cumuhaz&trt[cutidx]#(expbz1-expbz2); surv[j]=exp(-temp); end; f&trt=f&trt+surv; fexpbz&trt._1=fexpbz&trt._1+surv#expbz1; fexpbz&trt._2=fexpbz&trt._2+surv#expbz2; h=j(numtime,numcov,0); htemp=j(1, numcov, 0); do j=1 to cutidx; htemp=htemp + expbz1#count&trt[j]/s0&trt[j]*(ztemp1-s1&trt[j,]/s0&trt[j]); h[j,]=htemp; end; htemp2=j(1,numcov,0); do j=cutidx+1 to numtime; htemp2=htemp2 + expbz2#count&trt[j]/s0&trt[j]*(ztemp2-s1&trt[j,]/s0&trt[j]); h[j,]=htemp2+h[cutidx,]; end; fh&trt=fh&trt+surv#h; end; f&trt=f&trt/numobs; fexpbz&trt._1=fexpbz&trt._1/numobs; fexpbz&trt._2=fexpbz&trt._2/numobs; fh&trt=fh&trt/numobs; term1=j(numtime,1,0); do i=1 to cutidx; term1[i]=(fexpbz&trt._1[i]##2)#w&trt[i]; end; do i=cutidx+1 to numtime; term1[i]=(fexpbz&trt._1[i]##2)#w&trt[cutidx] + (fexpbz&trt._2[i]##2)#wFc&trt[i]; end; term2=j(numtime,1,0); do i=1 to numtime; term2[i]=fh&trt[i,]*sigma*t(fh&trt[i,]); end; var=term1+term2; se&trt=var##0.5; out=out||f&trt||se&trt; intt&trt=j(numeobs,1,0); do i=1 to numeobs; if strata[i]=&trt then intt&trt[i]=1/s0&trt[index[i]]; end; %end; /*****************************************************************/ /* */ /* Find SE's for differences between direct adjusted survivals */ /* */ /*****************************************************************/ %do trt1=1 %to &numgroup-1; %do trt2=%eval(&trt1+1) %to &numgroup; term1=j(numtime,1,0); do i=1 to cutidx; term1[i]=(fexpbz&trt1._1[i]##2)#w&trt1[i]; end; do i=cutidx+1 to numtime; term1[i]=(fexpbz&trt1._1[i]##2)#w&trt1[cutidx] + (fexpbz&trt1._2[i]##2)#wFc&trt1[i]; end; term2=j(numtime,1,0); do i=1 to cutidx; term2[i]=(fexpbz&trt2._1[i]##2)#w&trt2[i]; end; do i=cutidx+1 to numtime; term2[i]=(fexpbz&trt2._1[i]##2)#w&trt2[cutidx] + (fexpbz&trt2._2[i]##2)#wFc&trt2[i]; end; term3=j(numtime,1,0); do i=1 to numtime; term3[i]=(fh&trt1[i,]-fh&trt2[i,])*sigma*t(fh&trt1[i,]-fh&trt2[i,]); end; var=term1+term2+term3; se=var##0.5; out=out||se; se&trt1&trt2=se; do i=1 to numtime until(cumuhaz&trt1[i]>0.000001 & cumuhaz&trt2[i]>.0000001); end; start&trt1&trt2=i; %end; %end; names={"time" %do i=1 %to &numgroup; "surv&i" "se&i" %end; %do j=1 %to &numgroup; %do k=%eval(&j+1) %to &numgroup; "se&j._&k" %end; %end; }; create mout from out[colname=names]; append from out; close mout; /****************************************************************************/ /* */ /* Conduct simulation and find the maximum absolute value of each process */ /* */ /****************************************************************************/ eventzmat=j(numeobs, numcov,0); do j=1 to ecutidx; ztemp=zzall[j,]||tzall[j]||0; %do idx=1 %to &numgroup; if strata[j]=&idx then eventzmat[j,]=ztemp-s1&idx[index[j],]/s0&idx[index[j]]; %end; end; do j=ecutidx+1 to numeobs; ztemp=zzall[j,]||0||tzall[j]; %do idx=1 %to &numgroup; if strata[j]=&idx then eventzmat[j,]=ztemp-s1&idx[index[j],]/s0&idx[index[j]]; %end; end; do i=1 to &nsimu; term1=0; term2=0; term3=0; gvec=j(numeobs,1,0); do j=1 to numeobs; gvec[j]=rannor(0); end; allmart=j(1,numcov,0); do j=1 to numeobs; allmart=allmart + eventzmat[j,]#gvec[j]; end; %do trt1=1 %to &numgroup-1; %do trt2=%eval(&trt1+1) %to &numgroup; Gint&trt1=intt&trt1#gvec; Gint&trt2=intt&trt2#gvec; cumuG&trt1._1=j(numeobs,1,0); cumuG&trt1._2=j(numeobs,1,0); cumuG&trt2._1=j(numeobs,1,0); cumuG&trt2._2=j(numeobs,1,0); temp1=0; temp2=0; do j=1 to ecutidx; temp1=temp1+Gint&trt1[j]; temp2=temp2+Gint&trt2[j]; cumuG&trt1._1[j]=temp1; cumuG&trt2._1[j]=temp2; end; temp1=0; temp2=0; do j=ecutidx+1 to numeobs; temp1=temp1+Gint&trt1[j]; temp2=temp2+Gint&trt2[j]; cumuG&trt1._2[j]=temp1; cumuG&trt2._2[j]=temp2; end; max=0; do j=start&trt1&trt2 to cutidx; if time[j]>=&starttime & time[j]<=&stoptime then do; term1=cumuG&trt1._1[timeindex[j]]#fexpbz&trt1._1[j]; term2=cumuG&trt2._1[timeindex[j]]#fexpbz&trt2._1[j]; term3=(fh&trt1[j,]-fh&trt2[j,])*sigma*t(allmart); simu=(term1-term2+term3)/se&trt1&trt2[j]; if max=&starttime & time[j]<=&stoptime then do; term1=CumuG&trt1._1[ecutidx]#fexpbz&trt1._1[cutidx] + cumuG&trt1._2[timeindex[j]]#fexpbz&trt1._2[j]; term2=CumuG&trt2._1[ecutidx]#fexpbz&trt2._1[cutidx] + cumuG&trt2._2[timeindex[j]]#fexpbz&trt2._2[j]; term3=(fh&trt1[j,]-fh&trt2[j,])*sigma*t(allmart); simu=(term1-term2+term3)/se&trt1&trt2[j]; if max&stoptime then do; diff&trt1._&trt2._low=.; diff&trt1._&trt2._up=.; end; %end; %end; run; proc print; run; %mend;