/*********************************************************************************/ /* */ /* 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; */ /* 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 SURVBAND(inputdata, time, event, group, covlist, 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; close &inputdata; out=time||event||group||x; names={'time' 'event' 'strata' %do i=1 %to &numcov; "z&i" %end;}; 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;; strata strata; output out=coxout xbeta=zb; run; proc sort data=coxout; by strata descending time; run; data best sigma; set best; if _type_='PARMS' then output best; if _type_='COV' then output sigma; run; %do group=1 %to &numgroup; data riskset&group; set coxout (keep=time strata %do i=1 %to &numcov; z&i %end; zb); where strata=&group; by descending time; s0+exp(zb); %do i=1 %to &numcov; s1_&i+z&i * exp(zb); %end; if last.time; keep time s0 %do i=1 %to &numcov; 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 &numcov; ts1_&i %end; 999; if s0^=. then ts0=s0; else s0=ts0; %do i=1 %to &numcov; if s1_&i^=. then ts1_&i=s1_&i; else s1_&i=ts1_&i; %end; drop ts:; run; data dstrata&group; set indata; where strata=&group; run; %end; data riskset; merge %do i=1 %to &numgroup; riskset&i (rename=(s0=s0&i count=count&i %do j=1 %to &numcov; 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 &numcov; 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;} into b; close best; use sigma; read all var{%do i=1 %to &numcov; z&i %end;} into sigma; close sigma; use indata; read all var{%do i=1 %to &numcov; z&i %end;} into zmat; close indata; use event; 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; 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 trt=1 %to &numgroup; cumuhaz&trt=j(numtime,1,0); w&trt=j(numtime,1,0); ctemp=0; wtemp=0; do i=1 to numtime; 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; f&trt=j(numtime,1,0); fexpbz=j(numtime,1,0); fh=j(numtime,numcov,0); do i=1 to numobs; expbz=exp(zmat[i,]*t(b)); surv=exp(-cumuhaz&trt)##expbz; f&trt=f&trt+surv; fexpbz=fexpbz+surv#expbz; h=j(numtime,numcov,0); htemp=j(1, numcov, 0); do j=1 to numtime; htemp=htemp + count&trt[j]/s0&trt[j]*(zmat[i,]-s1&trt[j,]/s0&trt[j]); h[j,]=htemp; end; fh=fh+surv#h#expbz; end; f&trt=f&trt/numobs; fh&trt=fh/numobs; fexpbz&trt=fexpbz/numobs; term1=(fexpbz&trt##2)#w&trt; 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=(fexpbz&trt1##2)#w&trt1; term2=(fexpbz&trt2##2)#w&trt2; 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 numeobs; %do idx=1 %to &numgroup; if strata[j]=&idx then eventzmat[j,]=zzall[j,]-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=j(numeobs,1,0); cumuG&trt2=j(numeobs,1,0); temp1=0; temp2=0; do j=1 to numeobs; temp1=temp1+Gint&trt1[j]; temp2=temp2+Gint&trt2[j]; cumuG&trt1[j]=temp1; cumuG&trt2[j]=temp2; end; max=0; do j=start&trt1&trt2 to numtime; if time[j]>=&starttime & time[j]<=&stoptime then do; term1=cumuG&trt1[timeindex[j]]#fexpbz&trt1[j]; term2=cumuG&trt2[timeindex[j]]#fexpbz&trt2[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;