/*------------------------------------------------------------------* | The documentation and code below is supplied by HSR CodeXchange. | *------------------------------------------------------------------*/ /*------------------------------------------------------------------* | MACRO NAME : schoen | SHORT DESC : Returns Schoenfeld residuals (raw and scaled) | from the Cox model. *------------------------------------------------------------------* | CREATED BY : Therneau, Terry (04/14/2004 11:50) | : Bergstralh, Erik *------------------------------------------------------------------* | PURPOSE | | 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 *------------------------------------------------------------------* | OPERATING SYSTEM COMPATIBILITY | | UNIX SAS v8 : YES | UNIX SAS v9 : | MVS SAS v8 : | MVS SAS v9 : | PC SAS v8 : | PC SAS v9 : *------------------------------------------------------------------* | MACRO CALL | | %schoen ( | time= , | event= , | xvars= , | data=_last_, | strata= , | outsch=schr, | outbt=schbt, | plot=r, | vref=yes, | points=yes, | df=4, | pvars= , | alpha=.05, | rug=no, | ties=efron | ); *------------------------------------------------------------------* | REQUIRED PARAMETERS | | Name : time | Default : | Type : Variable Name (Single) | Purpose : time variable for survival | | Name : event | Default : | Type : Variable Name (Single) | Purpose : event variable for survival, 1=event, 0=censored | | Name : xvars | Default : | Type : Variable Name (List) | Purpose : list of covariates for Cox model | *------------------------------------------------------------------* | OPTIONAL PARAMETERS | | Name : data | Default : _last_ | Type : Dataset Name | Purpose : name of analysis dataset. Default is last dataset | created. | | Name : strata | Default : | Type : Variable Name (Single) | Purpose : stratification variable(one only) for stratifed Cox | models. | | Name : outsch | Default : schr | Type : Dataset Name | Purpose : 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. | | Name : outbt | Default : schbt | Type : Dataset Name | Purpose : 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. | | Name : plot | Default : r | Type : Text | Purpose : 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. | | Name : vref | Default : yes | Type : Text | Purpose : indicator to control plotting of a vertical reference | line at y=0. Values are yes(default) and no. | | Name : points | Default : yes | Type : Text | Purpose : yes,no. Indicates whether to plot the actual data points. | | Name : df | Default : 4 | Type : Number (Single) | Purpose : degrees of freedom for smoothing process Possible | values are 3 - 7. | | Name : pvars | Default : | Type : Variable Name (List) | Purpose : variables to plot. Default is all xzvars. | | Name : alpha | Default : .05 | Type : Number (Single) | Purpose : confidence coefficient for plotting standard | error bars. Default is .05. A value of 0 means | do not plot SE bars. | | Name : rug | Default : no | Type : Text | Purpose : indicator to control plotting of rug of x values. | Values are yes and no(default). | | Name : ties | Default : efron | Type : Text | Purpose : method used to break ties in phreg. Values are | efron(default), breslow, discrete, exact. | *------------------------------------------------------------------* | RETURNED INFORMATION | | 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. *------------------------------------------------------------------* | Copyright 2004 Mayo Clinic College of Medicine. | | This program is free software; you can redistribute it and/or | modify it under the terms of the GNU General Public License as | published by the Free Software Foundation; either version 2 of | the License, or (at your option) any later version. | | This program is distributed in the hope that it will be useful, | but WITHOUT ANY WARRANTY; without even the implied warranty of | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | General Public License for more details. *------------------------------------------------------------------*/ %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;