%macro cuminc(ds= , time= , cenvble= ,interest= ,group= )     ;
     **********************************************************;  
     *  this macro calculates:                                *;
     *        the CIF for the event of interest               *;     
     *        the variance for CIF based on the delta method  *;        
     *        the conditional probability  (CP)               *;    
     *        the variance for conditional probability based  *;
     *              on Pepe & Mori, stat in med, 1993         *;
     *  this macro creates the graph for                      *;
     *        the CIF for event of interest                   *;
     *        and calculates the conditional probability.     *;
     **********************************************************;
     *  macro parameters:                                     *;
     *                                                        *;
     *  1. ds                                                 *;
     *  the name of the input dataset.                        *;
     *                                                        *;
     *  2. time                                               *;
     *  the name of the time variable                         *;
     *                                                        *;
     *  3. cenvble.                                           *;
     *  the name of the variable with the censoring/event codes;  
     *  in order to use this macro, censored observations must*;
     *  be coded as 0 (zero) and events must have numeric     *;
     *  codes,although they need not be consecutive.          *;
     *                                                        *;
     *  4. interest                                           *;
     *  denotes the value of the censoring variable (cenvble) *;
     *  that is of primary interest.the macro produces printed*; 
     *  results including a listing of the cif for each group *;
     *  and a graph of the cif for the event of interest      *;
     *  for each group. it does not calculate or produce      *;
     *  output for the events other than the one of interest. *;
     *                                                        *;
     *  5. group                                              *;
     *  represents the grouping variable.                     *;
     *  the group codes can be either character or numeric.   *;           
     *                                                        *;
     *  the variance of CIF is calcualted based on            *;       
     *  the delta method                                      *;
     *  the variance of CP is calculated based on             *;         
     *   the formula given in the Pepe and Mori 1993 paper    *;           
     *                                                        *;
     **********************************************************;

title " "    ;
footnote " " ;

proc sort data = &ds (keep = &time &cenvble &group) 
out = sorted ;
by &group &time    ;
run                     ;

data sorted summary        ;
set sorted end = eof  ;
by &group &time       ;
retain grpcat 0 grpn freq totevt totcen totint tototh ;
     add1 = first.&group   ;
	 grpcat + add1       ;
if eof then call symput ("grpcat",grpcat) ;
if first.&group then grpn = 1 ;
else grpn = grpn+1  ;
if last.&group then call symput ("grpn"||trim(put(grpcat,1.)), grpn ) ;
if &cenvble = 0 then cenxyz = 0 ;
else if &cenvble = &interest then cenxyz = 1 ;
else if &cenvble not in (0, &interest) then cenxyz = 2 ;
if first.&time then do  ;
freq = 1 ;
totint = (cenxyz = 1)  ;
tototh = (cenxyz = 2)  ;
end ;
else do ;
freq = freq + 1 ;
totint = totint + (cenxyz = 1)  ;
tototh = tototh + (cenxyz = 2)  ;
end ;
totevt = totint + tototh  ;
output sorted  ;
if last.&time then output summary  ;
drop add1 grpn  ;
run  ;
data summary1;
set summary ;
by &group &time ;
if first.&group then do i = 1 to 1 ;
output ;
end ;
output ;
run ;
data summary2 ;
set summary1 ;
by &group   ;
retain surv atrisk sprev cif f2 t2 p4 p5 p1 p2 p3 v;
if i = 1 then do ;
&time=0 ; &cenvble=. ;
freq=0 ; totevt=0 ; 
totint=0 ; tototh = 0 ;
end ;
if first.&group then do ;
atrisk = input(compress(trim(symget("grpn"||trim(put(grpcat,1.))))),8.);
sprev = 1  ;
surv = 1-totevt/atrisk  ;
cif = totint / atrisk ;
f2 = tototh / atrisk ; 
t2 = (totint/atrisk**2)*(1-totint/atrisk) ;
p1 = totevt/(atrisk*(atrisk-totevt)) ;
p2 = cif*totevt/(atrisk*(atrisk-totevt)) ;
p3 = (cif**2)*totevt/(atrisk*(atrisk-totevt)) ;
p4 = totint/(atrisk**2) ;
p5 = (cif*totint)/(atrisk**2)  ;
v = (totint*(1-f2)**2+tototh*cif**2)/(atrisk*(atrisk-1));
atrisk = atrisk - freq        ;
end ;
else do ;
sprev = surv ;
surv = surv*(1-totevt/atrisk) ;
cif = cif + sprev*(totint/atrisk) ;
f2 = f2 + sprev*(tototh/atrisk) ;
t2 = t2 + (sprev**2)*(totint/atrisk**2)*(1-totint/atrisk) ;
p1 = p1 + totevt/(atrisk*(atrisk-totevt)) ;
p2 = p2 + cif*totevt/(atrisk*(atrisk-totevt)) ;
p3 = p3 + (cif**2)*totevt/(atrisk*(atrisk-totevt)) ;
p4 = p4 + sprev*totint/(atrisk**2) ;
p5 = p5 + (sprev*cif*totint)/(atrisk**2)  ;
if atrisk>1 then
v=v+(totint*(1-f2)**2+tototh*cif**2)/(atrisk*(atrisk-1));
else v=v;
atrisk = atrisk - freq        ;
end ;
drop &cenvble cenxyz i grpcat ;
run ;
data summary3 ;
set summary2 ;
vardelta = (cif**2)*p1 - 2*cif*p2 + p3 + t2 -2*cif*p4 + 2*p5 ;
cp      = cif /(1-f2) ;
varcp = (sprev**2/(1-f2)**4)*v;
run ;
proc report data = summary3 nowd colwidth=12
spacing = 2 headline headskip split="*"    ;
column &time atrisk totint totevt cif vardelta cp varcp;
define cif  / display "CIF for*type &interest *events" format=8.6 ;
define vardelta  / display "Variance*for CIF*type &interest *events" 
format=8.6; 
define totint  /display "Number*of*type &interest*events" format=6.   ;
define totevt  / display "Total*number*of events"  format=7.  ;
define atrisk  / display "Number*left" format=6. ;
define cp  / display "CP for*type &interest *events" format=8.6 ;
define varcp  / display "Variance*for CP*type &interest *events" format=8.6 ;
by &group ;
where totevt>0;
run ;
    symbol1 c=default i=stepslj l=1	w=2	v=none ;
    symbol2 c=default i=stepslj l=2	w=2	v=none ;
    symbol3 c=default i=stepslj l=3	w=2	v=none ;
    symbol4 c=default i=stepslj l=4	w=2	v=none ;
    symbol5 c=default i=stepslj l=5	w=2	v=none ;
    symbol6 c=default i=stepslj l=6	w=2	v=none ;
    symbol7 c=default i=stepslj l=7	w=2	v=none ;

axis2 color=black
label=("CIF" font=swissl h=1 c=black)
width=2.0
value=(font=swissl)
style=1 ;
legend1 position=(top left inside)
down = &grpcat
label=(font=swissl h=0.7 c=black)
value=(font=swissl h=0.7 c=black)
mode=protect ;

proc gplot data = summary3 ;
plot cif * &time =&group / 
haxis = axis1 vaxis = axis2 frame 
legend = legend1 ;
title h=1.5 c=black f=swissl "CIF for event &&interest by &&group" ;
run ;
axis2 color=black
label=("cp" font=swissl h=1 c=black)
width=2.0
value=(font=swissl)
style=1 ;
proc gplot data = summary3 ;
plot cp * &time =&group / 
haxis = axis1 vaxis = axis2 frame 
legend = legend1 ;
title h=1.5 c=black f=swissl "CP for event &&interest by &&group" ;
run ;
%mend ;