/****************************************************************/ /* S A S A U T O C A L L L I B R A R Y */ /* */ /* NAME: CumIncid */ /* TITLE: Estimate Cumulative Incidence */ /* PRODUCT: STAT */ /* SYSTEM: ALL */ /* KEYS: survival data, failure probability, competing risks */ /* PROCS: LIFETEST */ /* DATA: */ /* SUPPORT: sasycs UPDATE: 14Nov2002 */ /* REF: Cooley, T.A., Leisenring, W., Crowley, J. and */ /* Storer, B.E. (1999) "Estimation of Failure */ /* Probabilities in the Presence of Competing Risks: */ /* New Representation of Old Estimators." */ /* Statistics in Medicine, 18, 695-706 */ /* */ /* Kalbfleish, J.D. and Prentice, R.L. (1980) */ /* The Statistical Analyisis of Failure Time Data, */ /* John Wiley and Sons, New York */ /* */ /* Marubini, E. and Valsecchi, M.G. (1995) */ /* Analyzing Survival Data from Clinical Trials and */ /* Observational Studies, John Wiley & Sons, New York */ /* */ /* MISC: The following SAS products are required to run */ /* this macro: BASE and STAT. */ /****************************************************************/ /*------------------------------------------------------------------ The %CumIncid macro computes the crude cumulative incidence function estimates for homogenous (no covariates) survival data whose endpoints are subjected to competing risks (Kalbfleish and Prentice, 1980). Standard errors and pointwise confidence limits are also computed. The estimated crude cumulative incidence curve is displayed as a step function using ODS graphics. To display the graph, you specify odsgraphics=ON and the htmfile= or rtffile option. If you specify the STRATA= argument, the %CumIncid macro computes a separate crude cumulative incidence curve for each stratum and displays all the curves in the same graph. Optionally, you can request the cause-specific failure cumulative incidence estimates (Marubini and Valsecchi, 1995) by specifying CSE in the OPTIONS= argument. Cooley et al (1999) refers to it as 1-KM, which is interpreted as a hypothetical probability that assumes the probability of failure from the cause of interest would not change if competing risks were removed. Since there is no way of checking such assumption, 1-KM is deemed unintepretable when there are competing risks. Example ------- data a (drop=_i_); array type1(10) _temporary_ (1 13 17 30 34 41 78 100 119 169); array type2(20) _temporary_ (1 6 8 13 13 15 33 37 44 45 63 80 89 89 91 132 144 171 183 240); array type3(5) _temporary_ (34 60 63 149 207); do _i_=1 to 10; status=1; t=type1{_i_}; output; end; do _i_=1 to 20; status=2; t=type2{_i_}; output; end; do _i_=1 to 5; status=0; t=type3{_i_}; output; end; run; data b (drop=_i_); array type1(10) _temporary_ (7 16 16 20 39 49 56 73 93 113); array type2(20) _temporary_ (1 2 4 6 8 9 10 13 17 17 17 18 18 27 29 39 50 69 76 110); array type3(5) _temporary_ (34 60 63 78 149); do _i_=1 to 10; status=1; t=type1{_i_}; output; end; do _i_=1 to 20; status=2; t=type2{_i_}; output; end; do _i_=1 to 5; status=0; t=type3{_i_}; output; end; run; data ab; set a(in=ina) b; if ina then group=1; else group=2; run; %CumIncid(data=ab, time=t, status=status, event=1, compete=2, censored=0, strata=group, alpha=.05, htmfile= foo); Reference --------- Cooley, T.A., Leisenring, W., Crowley, J., and Storer, B.E. (1999). "Estimation of failure probabilities in the presence of competing risks: new representation of old estimators," Statistics in Medicine, 18, 695-706. Kalbfleish, J.D. and Prentice, R.L. (1980). The Statistical Analysis of Failure Time Data. Wiley. Marubini, E. and Valsecchi, M.G. (1995). Analysising Survival Data from Clinical Trials and Observational Studies. Wiley. ---------------------------------------------------------------------- DISCLAIMER: THIS INFORMATION IS PROVIDED BY SAS INSTITUTE INC. AS A SERVICE TO ITS USERS. IT IS PROVIDED "AS IS". THERE ARE NO WARRANTIES, EXPRESSED OR IMPLIED, AS TO MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE REGARDING THE ACCURACY OF THE MATERIALS OR CODE CONTAINED HEREIN. ----------------------------------------------------------------------*/ *===================== Macro Argument List Start ======================; %macro CumIncid( data=, /* specifies the name of the input data set */ out=, /* specifies the name of the output data set that */ /* contains the cumulative incidence estimates */ time=, /* specifies the time variable */ status=, /* specifies the variable whose values indicate whether */ /* an observation corresponds to the event of interest, */ /* the competing risk event, or is censored */ event=, /* specifies the values of status= that correspond to */ /* the event of interest */ compete=, /* specifies the values of status= that correspond to */ /* the competing risk events */ censored=, /* specifies the values of status= that correspond to */ /* censored observations */ strata=, /* specifies the list of strata variables */ alpha=, /* specifies the complement of the confidence level for */ /* pointwise confidence limits */ odsgraphics=, /* on or off, on for turning on the ODS graphics on, */ /* and the cumulative incidence curve will be displayed. */ /* odsgraphics=on is the default. */ htmfile=, /* output htm file without the .htm extension */ rtffile=, /* output rtf file without the .rtf extension */ options=, /* Specify one or more of following values after the */ /* OPTIONS=. */ /* */ /* CSE produces the cause-specific failure cumulative */ /* incidence estimates. This estimate is */ /* precisely the complement of the Kaplan-Meier */ /* estimate of the survivor function with */ /* competing risk events being treated as */ /* censored values. When this option is */ /* specified, standard errors and lower and */ /* upper confidence limits are also produced. */ /* */ /* PLOTCL */ /* plots the pointwise confidence limits for */ /* individual cumulative incidence curves. */ /* NOPRINT */ /* suppress print of the cumulative incidence */ /* estimates */ ); %global _debug_; %local cseflag noprintflag plotclflag Confidence ne n conflictsymbol svalueleng startclock nostimer savenote lapsedtime; %let startclock = %sysfunc(datetime()); %let savenote = %sysfunc(getoption(notes)); %if &_debug_ ne 1 %then options nonotes;; %if &odsgraphics = %then %do; %let odsgraphics=on; %end; %if &htmfile ne %then %do; ods html file="&htmfile..htm"; %end; %if &rtffile ne %then %do; ods rtf file="&rtffile..rtf"; %end; data _null_; /*----- check OPTIONS= -----*/ list= lowcase(symget('options')); i = index(list, 'cse'); if i then substr(list, i, 3) = ' '; call symput('cseflag', put(i gt 0, 1.)); i = index(list, 'noprint'); if i then substr(list, i, 7) = ' '; call symput('noprintflag', put(i gt 0, 1.)); i = index(list, 'plotcl'); if i then substr(list, i, 6) = ' '; call symput('plotclflag', put(i gt 0, 1.)); /*----- system options -----*/ call symput('nostimer', put(getoption('stimer' ) eq 'NOSTIMER', 1.)); /*----- Confidence Level -------*/ alpha=.05; %if &alpha ne %then %do; alpha=α %end; call symput('Confidence',trim(put(100*(1-alpha),best8. -l))); /*---- No of event= and compete= values ----*/ s=trim(left(translate(symget('event'),","," "))); ne= length(s) - length(compress(s," ")) + 1; call symput('ne',put(ne,best8. -l)); s=trim(left(translate(symget('compete'),","," "))); nc= length(s) - length(compress(s," ")) + 1; call symput('nc',put(nc,best8. -l)); run; data _null_; array event[&ne] _temporary_(&event); array compete[&nc] _temporary_(&compete); do i= 1 to ≠ do j= 1 to &nc; if compete[j]=event[i] then do; call symput('conflictsymbol',put(compete[j],best8. -l)); goto skip1; end; end; end; skip1:; run; %if &_debug_ ne 1 %then options &savenote;; %if &conflictsymbol ne %then %do; %let _xrc_= %str(%sysfunc(trim(&conflictsymbol)) is both an EVENT= and COMPETE= value); %put ERROR: &_xrc_..; %goto exit; %end; %if &_debug_ ne 1 %then options nonotes;; ods listing close; %if &htmfile ne %then %do; ods html select none; %end; %if &rtffile ne %then %do; ods rtf select none; %end; proc lifetest data=&data; ods output ProductLimitEstimates=_km1; %if &compete= and &censored= %then %do; time &time; %end; %else %do; time &time * &status(&compete &censored); %end; strata &strata; run; proc lifetest data=&data; ods output ProductLimitEstimates=_km2; %if &censored= %then %do; time &time; %end; %else %do; time &time * &status(&censored); %end; strata &strata; run; *proc print data=_km1; *proc print data=_km2; data _d1 (keep= &time Censor dj nj Sj_1 wj); retain m e nj lastSurv; set _km2; by notsorted &time; if &time=0 then do; lastSurv= 1.; nj=left; m=0; end; else do; if first.&time then do; nj=nj-m; m=0; e=0; end; m + 1; if censor=0 then e+1; if Survival= . then do; dj=0; wj=nj; end; else do; dj=e; wj=nj-dj; end; end; Sj_1= lastSurv; if (Survival ~= .) then lastSurv=Survival; run; *proc print data=_d1; *title "Data _1"; *run; data _d2; merge _km1 _d1(rename=(censor=censor2)); run; data _d2; label Censor='0=event 1=censored' failure= '1 - KM' CumInc= 'Crude Cumulative Incidence Estimate'; retain Sj_1 atrisk cj mj ej ci2 0; set _d2; by notsorted &time; if &time=0 then do; CumInc= 0; StdErrCumInc= 0; cj= 0; ci2= 0; hj= 0; atrisk=left; mj=0; *delete; end; else do; if first.&time then do; atrisk= atrisk - mj; ej=0; mj=0; end; mj + 1; if censor=0 then ej + 1; if survival ~= . then do; d1j= ej; hj= ej / atrisk; cj= cj + Sj_1 * hj; CumInc= cj; ci2=CumInc; end; else do; CumInc= .; d1j=0; end; end; run; *proc print data=_d2 label; *title "Data _d2"; *run; ods listing; %if &htmfile ne %then %do; ods html select all; %end; %if &rtffile ne %then %do; ods rtf select all; %end; data _d3 (keep= stratum &strata &time censor CumInc StdErrCumInc Lower&Confidence.PctCumInc Upper&Confidence.PctCumInc %if &cseflag %then %do; CSE StdErrCSE Lower&Confidence.PctCSE Upper&Confidence.PctCSE %end; ); retain nstart _z; %if &strata ne %then %do; retain &strata; %end; length &time 8 censor 8 CumInc 8 StdErrCumInc 8 Lower&Confidence.PctCumInc 8 Upper&Confidence.PctCumInc 8; %if &cseflag %then %do; length CSE 8 StdErrCSE 8 Lower&Confidence.PctCSE 8 Upper&Confidence.PctCSE 8; %end; alpha= .05; %if &alpha ne %then %do; alpha=α%end; _z= probit(.5*alpha); label CumInc='Crude Cumulative Incidence Estimate' StdErrCumInc='Standard Error of Crude Cumulative Incidence Estimator' Lower&Confidence.PctCumInc= "Lower &Confidence% Pointwise Confidence Limit for CumInc" Upper&Confidence.PctCumInc= "Upper &Confidence% Pointwise Confidence Limit for CumInc"; set _d2(rename=(FAILURE=CSE STDERR=StdErrCSE)); by stratum; if first.stratum then nstart=_n_+1; if d1j>0 then do; x=ci2; sigma2=0; nn= _n_; do i=nstart to nn; set _d2 point=i; b1= x - ci2; if b1=0 then c1= 0; else c1= b1 * b1 * dj / (nj * wj); c2= Sj_1 * Sj_1 * (nj-d1j) * d1j / (nj **3); c3= 2 * b1 * Sj_1 * d1j / (nj * nj); sigma2 = sigma2 + c1 + c2 - c3; end; StdErrCumInc= sqrt(sigma2); if CumInc= . or CumInc= 0. then do; Lower&Confidence.PctCumInc= .; Upper&Confidence.PctCumInc= .; end; else do; _d= exp(_z * StdErrCumInc / CumInc); Lower&Confidence.PctCumInc= CumInc * _d; Upper&Confidence.PctCumInc= CumInc / _d; end; %if &cseflag %then %do; if CSE= . or CSE= 0. then do; Lower&Confidence.PctCSE= .; Upper&Confidence.PctCSE= .; end; else do; _d= exp(_z * StdErrCSE / CSE); Lower&Confidence.PctCSE= CSE * _d; Upper&Confidence.PctCSE= CSE / _d; end; %end; end; else do; StdErrCumInc= .;end; run; %if &noprintflag = 0 %then %do; proc print data=_d3(drop=stratum); title "Cumulative Incidence Estimates"; run; title; %end; %if &out ne %then %do; %if &_debug_ ne 1 %then options &savenote;; %if not &nostimer %then options nostimer;; data &out; set _d3(drop=stratum); run; %if &_debug_ ne 1 %then options nonotes;; %if not &nostimer %then options stimer;; %end; %if &odsgraphics ne on %then %goto exit2; %if &strata ne %then %do; data _null_; set _d3 end=eof; retain __l __k 0; length __v $ 32; if _n_ = 1 then do; call vname(&strata, __v); __l = length(__v); end; __k = max(__k, length(trim(&strata))); if eof then do; call symput('svalueleng', put(__k + __l + 1, 2.)); call symput('nstrata', put(stratum, 2.)); end; run; %end; %else %let nstrata=1; data _d4; %if &strata ne %then %do; length __strata $ &svalueleng; %end; retain _last 0; set _d3; if CumInc= . then CumInc= _last; else _last=CumInc; %if &strata ne %then %do; *__strata=cat("&strata=",&strata); __strata="&strata="||trim(left(&strata)); %end; format &time best8.; run; ods path (prepend) work._foo; %if &_debug_ = 1 %then %do; ods path show; %end; %if &plotclflag %then %do; %do i=1 %to &nstrata; data _d5; set _d4(where=(stratum=&i)) end=eof; _lower= Lower&Confidence.PctCumInc; _upper= Upper&Confidence.PctCumInc; if eof then do; %if &strata= %then %do; %let svalue=; %end; %else %do; call symput('svalue',put(__strata,&svalueleng..)); %end; end; run; proc template; define statgraph _cicl; layout lattice / columns=2 hrange=unionall vrange=unionall; %if &strata ne %then %do; sidebar / align=top; entrytitle "&svalue"; endsidebar; %end; cell; cellheader; entrytitle "Cumulative Incidence Plot"; endcellheader; layout overlay / yaxisopts=(ticks=(0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1.0)); stepplot y=CumInc x=&time; endlayout; endcell; cell; cellheader; entrytitle "&Confidence% Pointwise CI's"; endcellheader; layout overlay / yaxisopts=(ticks=(0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1.0)); Band yLimitUpper=_upper yLimitLower=_lower x=&time / modelname="ci" lines=false fill=true; stepplot y=CumInc x=&time / name="ci"; endlayout; endcell; end; run; data _null_; set _d5; file print ods=(template="_cicl"); put _ods_; output; run; %end; %end; %if &strata= %then %do; proc template; define statgraph _ci; layout gridded; entrytitle 'Cumulative Incidence Plot'; layout overlay / yaxisopts=(ticks=(0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1.0)); stepplot y=CumInc x=&time; endlayout; endlayout; end; run; %end; %else %do; proc template; define statgraph _ci; layout gridded / width=640 height=480; entrytitle 'Cumulative Incidence Plot'; layout overlay / yaxisopts=(ticks=(0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1.0)); stepplot y=CumInc x=&time / group=__strata name="ci"; discretelegend "ci" / hAlign=right vAlign=bottom across=1 border=true; endlayout; endlayout; end; run; %end; data _null_; set _d4; file print ods=(template="_ci"); put _ods_; output; run; ods path (remove) work._foo; %if &_debug_ = 1 %then %do; ods path show; %end; %if &odsgraphics=on %then %do; ods graphics off; %end; %exit2: %if &htmfile ne %then %do; ods html close; %end; %if &rtffile ne %then %do; ods rtf close; %end; %if &_debug_ ne 1 %then %do; proc datasets nolist; delete _d1 _d2 _d3 %if &odsgraphics = on %then %do; _d4 _foo(mt=itemstor) %if &plotclflag %then _d5; %end; ; quit; %end; %exit: options &savenote; %if not &nostimer %then %do; %let lapsedtime = %sysfunc(round(%sysevalf(%sysfunc(datetime()) - &startclock), 0.01)); %if %sysevalf(&lapsedtime >= 60) %then %do; %let lapsedtime = %sysfunc(round(%sysevalf(&lapsedtime / 60), 0.01)) minutes.; %end; %else %let lapsedtime = &lapsedtime seconds.; %put NOTE: The CumIncid macro used &lapsedtime; %end; %mend; %mend CumIncid; /****************/ /* end of macro */ /****************/