%macro compcif(ds=, time=, cens=, group=,val1=,val2=);

******************************************************;
* It compares the CIF between two groups             *;
* It is based on the test described in               *;
* Pepe and Mori paper Statistics in Medicine 1993    *;
******************************************************;
* ds=data set                                        *;
* time= time to first event observed                 *;
* cens= censor variable coded as:                    *;
*           1= event of interest                     *;
*           2= competing risks type of event         *;
*           0= censored (no failure observed)        *;
* group=group variable coded as val1 and val2        *;
******************************************************;

data d1;set &ds;
if &cens in (0,2) then c02=1;else c02=0;
if &cens ^ in (0,1) then dcr=1;else dcr=0;
if &cens=1 then dev=1;else dev=0;

data dg1;set d1;if &group=&val1;keep &time c02 dev dcr;
proc means data=dg1 sum noprint;
var c02 dev dcr ;class &time;
output out=dg11 sum=c021 dev1 dcr1 ;
data dg11;set dg11;if _type_=1;keep &time c021 dev1 dcr1 ;

data dg2;set d1;if &group=&val2; keep &time c02 dev dcr;
proc means data=dg2 sum noprint;
var c02 dev dcr ;class &time;
output out=dg21 sum=c022 dev2 dcr2 ;
data dg21;set dg21;if _type_=1;keep &time c022 dev2 dcr2 ;

proc means data=d1 noprint;var &time;class &group;
output out=d2 max=tau;
data d21;set d2;if &group=&val1;rename _freq_=n1 tau=tau1;keep _freq_ tau;
data d22;set d2;if &group=&val2;rename _freq_=n2 tau=tau2;keep _freq_ tau;
data d2;merge d21 d22;tau=min(tau1,tau2);keep n1 n2 tau;
data dg;merge dg11 dg21;by &time;
data dg;set dg;
if dev1=. then dev1=0;if dcr1=. then dcr1=0;if c021=. then c021=0;
if dev2=. then dev2=0;if dcr2=. then dcr2=0;if c022=. then c022=0;
data dg;set dg;if _n_=1 then set d2;if &time<=tau;
data dg;set dg;retain atrisk1 s1 c1 f1 fcr1 atrisk2 s2 c2 f2 fcr2 ;
c021m=lag(c021);dev1m=lag(dev1);c022m=lag(c022);dev2m=lag(dev2);
if _n_=1 then do;
atrisk1=n1;
s1=(atrisk1-dev1-dcr1)/atrisk1;
c1=(atrisk1-c021)/atrisk1;
f1=dev1/atrisk1;
fcr1=dcr1/atrisk1;
atrisk2=n2;
s2=(atrisk2-dev2-dcr2)/atrisk2;
c2=(atrisk2-c022)/atrisk2;
f2=dev2/atrisk2;
fcr2=dcr2/atrisk2;
end;
else do;
atrisk1=atrisk1-c021m-dev1m;
s1=s1*(atrisk1-dev1-dcr1)/atrisk1;
c1=c1*(atrisk1-c021)/atrisk1;
atrisk2=atrisk2-c022m-dev2m;
s2=s2*(atrisk2-dev2-dcr2)/atrisk2;
c2=c2*(atrisk2-c022)/atrisk2;
end;
s1m=lag(s1);s2m=lag(s2); c1m=lag(c1);c2m=lag(c2);
if _n_=1 then do;
s1m=1;s2m=1; c1m=1;c2m=1;
end;

if _n_>1 then do;
f1=f1+dev1/atrisk1*s1m;
fcr1=fcr1+dcr1/atrisk1*s1m;
f2=f2+dev2/atrisk2*s2m;
fcr2=fcr2+dcr2/atrisk2*s2m;
end;

w=c1m*c2m*(n1+n2)/(n1*c1m+n2*c2m);
data temp;set dg;keep &time;
proc sort data=temp;by descending &time ;
data temp;set temp;timem=lag(&time);
proc sort data=temp;by &time ;

data dg;merge temp dg;by &time;
data dg;set dg;
retain v1part v21part v22part;
if timem=. then delete;
deltat=timem-&time;
spart=w*(f1-f2)*deltat;

** variance;
if _n_=1 then do;
v1i=w*deltat;
v21i=v1i*f1;
v22i=v1i*f2;
v1part=w*deltat;
v21part=v1i*f1;
v22part=v1i*f2;
end;
else do;
v1i=w*deltat;
v21i=v1i*f1;
v22i=v1i*f2;
v1part=v1part+w*deltat;
v21part=v21part+v1i*f1;
v22part=v22part+v1i*f2;
end;

data dlast;set dg end=last;if last;
keep v1part v21part v22part;rename v1part=v1tot v21part=v21tot v22part=v22tot;
data dg;set dg;if _n_=1 then set dlast;
data dg;set dg;
v1=v1tot-v1part+v1i;
v21=v21tot-v21part+v21i;
v22=v22tot-v22part+v22i;
v11=(1-fcr1)*v1-v21;
v12=(1-fcr2)*v1-v22;
sigma1=(v11**2*dev1+v21**2*dcr1)/(atrisk1*(atrisk1-1));
sigma2=(v12**2*dev2+v22**2*dcr2)/(atrisk2*(atrisk2-1));

proc means data=dg noprint;var spart sigma1 sigma2;
output out=var sum=s sigma1 sigma2;
data var;set var;keep s sigma1 sigma2;

data result;merge var d2;
s=s*sqrt(n1*n2/(n1+n2));
sigma=n1*n2/(n1+n2)*(sigma1+sigma2);
chisq=s**2/sigma;
pvalue=1-probchi(chisq,1);


proc report data=dg nowd colwidth=12 spacing=2 headline headskip split="*";
column &time atrisk1 dev1 dcr1 f1 fcr1 ;
define &time/display "Time" format=7.2;
define atrisk1/display "Number*at risk" format=7.;
define dev1/display "Number*events of*interest" format=9.;
define dcr1/display "Number*competing*risks" format=9.;
define f1/display "CIF*event of*interest" format=8.6;
define fcr1/display "CIF*competing*risk" format=9.6;
title1 "Estimates for group &val1";
run;

proc report data=dg nowd colwidth=12 spacing=2 headline headskip split="*";
column &time atrisk2 dev2 dcr2 f2 fcr2 ;
define &time/display "Time" format=7.2;
define atrisk2/display "Number*at risk" format=7.;
define dev2/display "Number*events of*interest" format=9.;
define dcr2/display "Number*competing*risks" format=9.;
define f2/display "CIF*event of*interest" format=8.6;
define fcr2/display "CIF*competing*risk" format=9.6;
title1 "Estimates for group &val2";
run;



proc report data=result nowd colwidth=12 spacing=2 headline headskip split="*";
column n1 n2 s chisq pvalue;
define n1/display "N*Group*&val1" format=7.;
define n2/display "N*Group*&val2" format=7.;
define s/display "Score" format=13.4;
define chisq/display "Chi-square" format=13.4;
define pvalue/display "p-value" format=7.5;
%mend;