/* ***************************************************************************; Macro name: %schoen Function: Returns Schoenfeld residuals and scaled Schoenfeld residuals from the Cox model. The Schoenfeld residuals(r) are defined only at event times. There is one residual for each covariate in the Cox model. The scaled residuals(Bt) are defined as: Bt= B + D*cov(B)*r', B=cox model beta, D=total #events and r=Schoenfeld residuals. Bt is an estimate of the hazard function at follow-up time t. Therefore a plot of Bt vs time can be used to assess whether the proportional hazards assumption is valid. As time is frequently skewed, plots of Bt vs rank time or '1-Kaplan- Meier' for the entire dataset(which is a monotone function of time) may be preferred. The correlation of Bt with time provides a test of the proportional hazards assumption. Developers: E. Bergstralh and T. Therneau Section of Biostatistics Mayo Clinic Rochester, Mn 55905 e-mail: bergstra@mayo.edu therneau@mayo.edu Date: May 6, 1993 Call statement: %schoen(data=,time=,event=,xvars=, strata=,outsch=,outbt=,plot=, vref=,points=,df=,pvars,alpha=,rug=,ties=); Parameters: Required: time=time variable for survival event=event variable for survival, 1=event, 0=censored xvars=list of covariates for Cox model Optional: data= name of analysis dataset. Default is last dataset created. strata=stratification variable(one only) for stratifed Cox models. outsch=name of output data set containing Schoenfeld residuals. One obs per each event time. The variables containing the residuals have the same name as the covariates (xvars). The data set also includes the time variable and the strata variable. Default data set name is schr. outbt =name of output data set containing the scaled Schoenfeld residuals(Bt). One obs per each event time. The variables containing the scaled residuals have the same name as the covariates (xvars). The dataset also includes the time variable, the rank of the time(rtime) and a variable(probevt) which is equal to '1 - overall Kaplan-Meier' at the given time. Default data set name is schbt. plot=t,r,k,n. Indicates that SAS/Graph plots of Bt vs time(t), rank of time(r) or '1 - overall Kaplan-Meier' (k) are to be done. Default is r. For no plots use n. The name of the graphics catalog is gschbt. vref= indicator to control plotting of a vertical reference line at y=0. Values are yes(default) and no. points=yes,no. Indicates whether to plot the actual data points. Default is yes. df= degrees of freedom for smoothing process Possible values are 3 - 7. Default is 4. pvars= variables to plot. Default is all xzvars. alpha= confidence coefficient for plotting standard error bars. Default is .05. A value of 0 means do not plot SE bars. rug= indicator to control plotting of rug of x values. Values are yes and no(default). ties= method used to break ties in phreg. Values are efron(default), breslow, discrete, exact. Output: The macro prints the PHREG output used to fit the Cox model, correlation coefficients of Bt vs time and SAS/Graph plots of Bt vs time. References: Schoenfeld, D. (1982). Partial residuals for the proportional hazards regression model. Biometrika 69, 239-41. Grambsch PM, Therneau TM (1993). Proportional hazards tests and diagnostics based on weighted residuals. University of Minnesota, Division of Biostatistics. Research Report 93-002. ***************************************************************************; */ %macro schoen(data=_last_,time=,event=,xvars=,vref=yes, strata=,outsch=schr,outbt=schbt,plot=r,points=yes, df=4,pvars=,alpha=.05,rug=no,ties=efron); footnote"schoen macro: event=&event time=&time strata=&strata"; footnote2"Xvars= &xvars"; %let bad=0; %if &time= %then %do; %put ERROR: NO TIME VARIABLE GIVEN.; %let bad=1; %end; %if &event= %then %do; %put ERROR: NO EVENT VARIABLE GIVEN.; %let bad=1; %end; %if &xvars= %then %do; %put ERROR: NO XVARS GIVEN.; %let bad=1; %end; %if &df<3 | &df>7 %then %do; %put ERROR: DF MUST BE BETWEEN 3 AND 7.; %let bad=1; %end; %if %upcase(&points)^=YES & %upcase(&points)^=NO %then %do; %put ERROR: POINTS MUST BE YES OR NO.; %let bad=1; %end; %if %upcase(&rug)^=YES & %upcase(&rug)^=NO %then %do; %put ERROR: RUG MUST BE YES OR NO.; %let bad=1; %end; %if %upcase(&vref)^=YES & %upcase(&vref)^=NO %then %do; %put ERROR: VREF MUST BE YES OR NO.; %let bad=1; %end; %if %upcase(&plot)^=T & %upcase(&plot)^=R & %upcase(&plot)^=K & %upcase(&plot)^=N %then %do; %put ERROR: PLOT MUST BE T, R, K, OR N.; %let bad=1; %end; %if %upcase(&ties)^=EFRON & %upcase(&ties)^=BRESLOW & %upcase(&ties)^=DISCRETE & %upcase(&ties)^=EXACT %then %do; %put ERROR: TIES MUST BE EFRON, BRESLOW, DISCRETE, EXACT.; %let bad=1; %end; %let badalpha=0; data _null_; a=symget('alpha'); if a<0 | a>=1 then call symput('badalpha',1); if a=0 then doalpha=0; else doalpha=1; call symput('doalpha',doalpha); run; %if &badalpha=1 %then %do; %put ERROR: ALPHA MUST BE BETWEEN 0 AND 1; %let bad=1; %end; %let nxs=0; %*count the number of x vars; %do %until(%scan(&xvars,&nxs+1,' ')= ); %let nxs=%eval(&nxs+1); %end; %if &pvars^= %then %let plotvars=&pvars; %else %let plotvars=&xvars; %let nplotvar=0; %do %until(%scan(&plotvars,&nplotvar+1,' ')= ); %let nplotvar=%eval(&nplotvar+1); %end; %do i=1 %to &nplotvar; %let pv=%scan(&plotvars,&i,' '); %let found=0; %do j=1 %to &nxs; %if &pv=%scan(&xvars,&j,' ') %then %let found=1; %end; %if &found=0 %then %do; %put ERROR: PLOT VARIABLE &PV NOT ON XVARS LIST.; %let bad=1; %end; %end; %macro xlist(prefix); %*set-up var list code; %do j=1 %to &nxs; &prefix&j %end; %mend xlist; %if &bad=0 %then %do; data _setup; set &data; %*delete obs with mising values; xmiss=0; %do i=1 %to &nxs; xx=%scan(&xvars,&i); if xx=. then xmiss=1; %end; if &event=. or &time=. or xmiss=1 then delete; %* run phreg and grab output datasets; proc phreg data=_setup covout outest=_est ; model &time*&event(0)= &xvars / ties=&ties ; %if &strata ^= %then %do; strata &strata; %end; output out=_sch xbeta=xbeta ressch=rsch1-rsch&nxs wtressch=wrsch1-wrsch&nxs; data _sch1; set _sch(keep=&strata &time &event rsch1-rsch&nxs); drop &event; if &event=1; rename %do i=1 %to &nxs; %let thisx=%scan(&xvars,&i,' '); rsch&i=&thisx %end; ; proc sort; by &time; data &outsch; set _sch1; proc sort; by &strata &time; ** calculate overall Kaplan-Meier; proc lifetest noprint data=_setup outs=_km; time &time*&event(0); data _km2; set _km; keep &time probevt; if _censor_=0; **keep event times; probevt=1-survival; label probevt='1-Overall Kaplan-Meier'; data _null_; set _est; %if %substr(&sysver,1,1)=6 %then %do; if _type_='PARMS' & _name_='ESTIMATE'; %end; %if %substr(&sysver,1,1)=8 %then %do; if _type_='PARMS' & upcase(_name_)=upcase("&time") ; %end; %do i=1 %to &nxs; %let thisx=%scan(&xvars,&i,' '); call symput("est&i",&thisx); %end; data _sch2; set _sch(keep=&strata &time &event wrsch1-wrsch&nxs); keep &strata &time &xvars; if &event=1; %do i=1 %to &nxs; %let thisx=%scan(&xvars,&i,' '); &thisx=&&est&i + wrsch&i; %end; proc sort; by &time; data _bt; merge _sch2(in=ins) _km2(keep=&time probevt); by &time; if ins; proc sort; by &strata &time; run; symbol1 i=join l=1 c=black; symbol2 i=join l=2 c=black; symbol3 i=join l=2 c=black; symbol4 i=none v=dot h=.1 cm c=black; %macro dovars(tvar); %if %length(&tvar)=8 %then %let t=%substr(&tvar,1,7); %else %let t=&tvar; %daspline(&tvar,nk=&df,data=&outbt) %do i=1 %to &nplotvar; %let x=%scan(&plotvars,&i,' '); data __temp1; set &outbt; &&_&t %global _print_; %let _print_=off; %let ndummies=%eval(&df-2); %if %substr(&sysver,1,1)=8 %then %do; ods listing close; %end; proc genmod; model &x=&tvar &t.1-&t.&ndummies/dist=normal obstats %if &doalpha=1 %then %do; alpha=&alpha %end; ; make 'OBSTATS' out=__temp2; run; %if %substr(&sysver,1,1)=8 %then %do; ods listing; %end; data __temp3; merge __temp1 __temp2; /* lower=xbeta-2*std; upper=xbeta+2*std; */ %if %upcase(&rug)=YES %then %do; data __anno; set __temp3; x=&tvar; y=0; xsys='2'; ysys='1'; function='move'; output; y=5; function='draw'; output; %end; proc sort data=__temp3; by &tvar; proc gplot data=__temp3 gout=gschbt %if %upcase(&rug)=YES %then %do; annotate=__anno %end; ; %if &doalpha=1 %then %do; plot xbeta*&tvar=1 lower*&tvar=2 upper*&tvar=3 %if %upcase(&points)=YES %then %do; &x*&tvar=4 %end; /overlay vaxis=axis1 haxis=axis2 %if %upcase(&vref)=YES %then %do; vref=0 lvref=3 %end; ; %end; %if &doalpha=0 %then %do; plot xbeta*&tvar=1 %if %upcase(&points)=YES %then %do; &x*&tvar=4 %end; /overlay vaxis=axis1 haxis=axis2 %if %upcase(&vref)=YES %then %do; vref=0 lvref=3 %end; ; %end; axis1 label=(r=0 a=90 "smooth(&x) df=&df"); %if &tvar=&time %then %do; axis2 label=("&tvar"); %end; %if &tvar=rtime %then %do; axis2 label=("Rank for Variable &time"); %end; %if &tvar=probevt %then %do; axis2 label=('1-Overall Kaplan-Meier'); %end; run; quit; %end; %mend dovars; title2'Scaled residuals(Bt) as a fcn of time.'; proc rank data=_bt out=&outbt; var &time; ranks rtime; proc corr pearson data=&outbt; with &xvars; var &time rtime probevt; label probevt='1 - Overall Kaplan-Meier'; %if %upcase(&plot)=T %then %do; %dovars(&time) %end; %if %upcase(&plot)=R %then %do; %dovars(rtime) %end; %if %upcase(&plot)=K %then %do; %dovars(probevt) %end; run; quit; footnote; symbol; title2; proc datasets nolist; delete _setup _est _sch _sch1 _sch2 _km _km2 _bt %if %upcase(&rug)=YES %then %do; __anno %end; __temp1 __temp2 __temp3; run; quit; %end; %mend schoen;