/************************************************************************/
/*                                                                      */
/*                                                                      */
/*  This macro computes the direct adjusted cumulative incidence curves */
/*  for K treatment groups, based on a proportional subdistribution     */
/*  hazards model. Please note that the hazards ratio between any two   */
/*  treatments is assumed to be a constant.                             */
/*                                                                      */
/*  Macro parameters:                                                   */
/*     inputdata - the input sas data name;                             */
/*     time      - the survival time variable;                          */
/*     event     - 0 (censor) 1 (cause of interest) 2 (other causes);   */
/*     group     - the treatment group variable,                        */
/*                      which must take values 1,...,K for K groups.    */
/*     covlist   - a list of covariate names;                           */
/*     outdata   - the output sas data name.                            */
/*                                                                      */
/*  The output dataset contains:                                        */
/*     Time - the event times                                           */
/*     CIFi, i=1,...,K                                                  */
/*     SEi, i=1,...,K                                                   */
/*     SEi_j, 1<=i<j<=K                                                 */
/*                                                                      */
/************************************************************************/

%macro DACIF(inputdata,time, event, group, covlist, outdata);

proc means data=&inputdata noprint;
    var &group;
    output out=maxout max(&group)=numgroup;
run;


data _null_;
    set maxout;
    call symput('numgroup', numgroup);	
run;


proc iml;
    use &inputdata;
    read all var {&covlist} into x;
    close &inputdata;

    numcov=ncol(x);

    create ncovout from numcov[colname='numcov'];
    append from numcov;
    close ncovout;
quit;

data _null_;
    set ncovout;
    call symput('numcov', numcov);
run;

%let tcov=%eval(&numgroup-1+&numcov);

proc iml;
    use &inputdata;
    read all var {&time} into time;
    read all var {&event} into event;
    read all var {&group} into group;
    read all var {&covlist} into x;
    close &inputdata;    
 
    numobs=nrow(time);

    gmat=j(&numgroup, &numgroup-1, 0);
    do i=2 to &numgroup;
 	gmat[i, i-1]=1;
    end;

    zmat=j(numobs, &numgroup-1+&numcov, 0);
    do i=1 to numobs;
	zmat[i,]=gmat[group[i],]||x[i,];
    end;

    out=time||event||zmat;

    names={'time' 'cause' %do i=1 %to &numgroup-1+&numcov; "z&i" %end;};
    create indata from out[colname=names];
    append from out;
    close indata;
quit;
 
proc sort data=indata; by descending time descending cause; run;

data indata;
    set indata;
    atrisk+1;
    run;

proc sort data=indata; by time; run;

data cdata;
    set indata; where cause=0;
    keep time atrisk;
run;

proc sort data=cdata; by time descending atrisk; run;

proc freq data=cdata noprint;
	tables time /out=ccount;
run;

proc sort data=ccount; by time; run;

data cdata;
    merge cdata ccount; by time;
    if first.time;
    retain ctemp 1;
    ctemp=ctemp*(1-count/atrisk);
    csurv=ctemp;
    dcna=count/atrisk;
    drop percent ctemp;
    ctime=time;
    time=time+0.00001;
run;

data indata;
    merge indata (in=in1) cdata (keep=time csurv); by time;
    retain ctemp 1;
    if csurv^=. then ctemp=csurv;
    else csurv=ctemp;
    if in1;
run;

data e1time;
    set indata; where cause=1;
    keep time csurv;
run;

proc sort; by time; run;

data e1time;
    set e1time; by time;
    if first.time;
run;

proc phreg  data=indata out=best noprint;
    model time*cause(0,2)=%do i=1 %to &numgroup-1+&numcov; z&i %end; ;
run;	



%let tnum=%eval(&tcov*(&tcov+1)/2);

data obs10 obs2;
    set indata;
    if cause in (0,1) then output obs10;
    else output obs2;
run;

proc iml;

    use obs10;
    read all var{time} into time10;
    read all var{cause} into cause10;
    read all var{%do i=1 %to &numgroup-1+&numcov; z&i %end;} into z10;
    read all var{%do i=&numgroup %to &numgroup-1+&numcov; z&i %end;} into zzmat01;
    close obs10;

    use obs2;
    read all var{time} into time2;
    read all var{csurv} into ckm2;
    read all var{cause} into cause2;
    read all var{%do i=1 %to &numgroup-1+&numcov; z&i %end;} into z2;
    read all var{%do i=&numgroup %to &numgroup-1+&numcov; z&i %end;} into zzmat2;
    close obs2;

    use best;
    read all var{%do i=1 %to &numgroup-1+&numcov; z&i %end;} into best;
    close best;


    use e1time;
    read all var{time} into etime;
    read all var{csurv} into eckm;
    close e1time;    

    use cdata;
    read all var{ctime} into ctime;
    read all var{dcna} into dcna;
    read all var{atrisk} into pi;
    close cdata;

    numctime=nrow(ctime);    
            
    numtime=nrow(etime);

    numobs10=nrow(time10);
    numobs2=nrow(time2);
 
 
    cidx=j(numctime,1,0);
    do i=1 to numctime;
	if ctime[i]>etime[numtime] then cidx[i]=-1;
    	else do;
	    do j=1 to numtime until(etime[j]>=ctime[i]);
	    end;
	    cidx[i]=j;
	end;
    end;


    eidx=j(numtime,1,0);
    do i=1 to numtime;    	
	do j=1 to numobs10 until(time10[j]=etime[i]);
	end;
	eidx[i]=j;
    end;
     

    b=best; 

    numdeath=j(numtime,1,0);
    covsum=j(numtime,&tcov,0);    
    
    do i=1 to numtime;            
        do j=1 to numobs10;
	    if time10[j]=etime[i] & cause10[j]=1  then do;
                numdeath[i] = numdeath[i] + 1;
                covsum[i,] = covsum[i,] + z10[j,];            
	    end;
	end;
    end;


    incre=1;

do iter=1 to 20 until(incre<.0005);    

    s0_1=j(numtime,1,0);
    s1_1=j(numtime,&tcov,0);
    s2_1=j(numtime,&tnum,0);

    score=j(1,&tcov,0);
    fisher=j(&tcov,&tcov,0);
    loglike=0;        

    expbz10=j(numobs10,1,0);
    do i=1 to numobs10;
	expbz10[i]=b*t(z10[i,]);
    end;
    expbz10=exp(expbz10);
    expbz2=j(numobs2,1,0);
    do i=1 to numobs2;
	expbz2[i]=b*t(z2[i,]);
    end;
    expbz2=exp(expbz2);
        

    pres0=0;
    pres1=j(1,&tcov,0);
    pres2=j(1,&tnum,0);
    obsidx=numobs10;
    do i=numtime to 1 by -1;
	s0_1[i]=pres0;
	s1_1[i,]=pres1;
	s2_1[i,]=pres2;
	do j=eidx[i] to obsidx;
            s0_1[i] = s0_1[i] + expbz10[j];
            s1_1[i,] = s1_1[i,] + z10[j,]#expbz10[j];
            tonext=1;
            do m=1 to &tcov;
                do n=m to &tcov;
                    s2_1[i,tonext] = s2_1[i,tonext] + z10[j,m]#z10[j,n]#expbz10[j];
                    tonext = tonext + 1;
                end;
            end;
	end;
	pres0=s0_1[i];
    	pres1=s1_1[i,];
	pres2=s2_1[i,];
	obsidx=eidx[i]-1;
    end;

    s0_2=j(numtime,1,0);
    s1_2=j(numtime,&tcov,0);
    s2_2=j(numtime,&tnum,0);
    do i=1 to numtime;
	do j=1 to numobs2;
	    cweight=eckm[i]/ckm2[j];
	    weight=min(1,cweight);
            s0_2[i] = s0_2[i] + expbz2[j]#weight;
            s1_2[i,] = s1_2[i,] + z2[j,]#expbz2[j]#weight;
            tonext=1;
            do m=1 to &tcov;
                do n=m to &tcov;
                    s2_2[i,tonext] = s2_2[i,tonext] + z2[j,m]#z2[j,n]#expbz2[j]#weight;
                    tonext = tonext + 1;
                end;
            end;
	end;
    end;


    s0=s0_1+s0_2;
    s1=s1_1+s1_2;
    s2=s2_1+s2_2;

    do i=1 to numtime;

        score = score + covsum[i,] - s1[i,]#(numdeath[i]/s0[i]);

        tonext=1;
        do m=1 to &tcov;
            do n=m to &tcov;
                fisher[m,n] = fisher[m,n] + s2[i,tonext]#(numdeath[i]/s0[i])
                    - s1[i,m]#s1[i,n]#(numdeath[i]/s0[i]##2);
                tonext = tonext + 1;
            end;
        end;
	
        loglike = loglike + b*t(covsum[i,])-numdeath[i]*log(s0[i]); 
    end;       


	do m=2 to &tcov;
	    do n=1 to m-1;
		fisher[m,n]=fisher[n,m];
	    end;
	end;  



	oldb=b;
    b = b + score*inv(fisher);
    
	bchange=b-oldb;
    incre=max(abs(bchange));

end;

    do i=1 to numobs10;
	expbz10[i]=b*t(z10[i,]);
    end;
    expbz10=exp(expbz10);
    do i=1 to numobs2;
	expbz2[i]=b*t(z2[i,]);
    end;
    expbz2=exp(expbz2);



    lambda1=j(numtime,1,0);    
    dlambda1=j(numtime,1,0);
    ltemp=0;
    do i=1 to numtime;
        dlambda1[i]=numdeath[i]/s0[i];
        ltemp=ltemp+dlambda1[i];
        lambda1[i]=ltemp;
    end;


    
    eta_1=j(numobs10,&tcov,0);
    dmart1_1=j(numobs10,numtime,0);
    do i=1 to numobs10;
        do j=1 to numtime;
            if time10[i]>=etime[j] then dmart1_1[i,j]=-dlambda1[j]#expbz10[i];
            if cause10[i]=1 & time10[i]=etime[j] then dmart1_1[i,j]=dmart1_1[i,j]+1;
            eta_1[i,] = eta_1[i,] + (z10[i,]-s1[j,]/s0[j])*dmart1_1[i,j];
 	end;
    end;

    eta_2=j(numobs2,&tcov,0);
    dmart1_2=j(numobs2,numtime,0);
    do i=1 to numobs2;
        do j=1 to numtime;
	    cweight=eckm[j]/ckm2[i];
	    weight=min(1,cweight);
            dmart1_2[i,j]=-dlambda1[j]#expbz2[i]#weight;
            eta_2[i,] = eta_2[i,] + (z2[i,]-s1[j,]/s0[j])*dmart1_2[i,j];
        end;
    end;
    
    numobs=numobs10+numobs2;
    dmart1=dmart1_1//dmart1_2;
    time=time10//time2;
    cause=cause10//cause2;
    zmat=z10//z2;
    eta=eta_1//eta_2;

    psi=j(numobs,&tcov,0);


    q=j(numctime,&tcov,0);
    do i=1 to numctime while (cidx[i]>0);
	do j=1 to numobs2 while (ctime[i]>time2[j]);
	    do k=cidx[i] to numtime;
                q[i,]=q[i,]+(z2[j,]-s1[k,]/s0[k])#dmart1_2[j,k];
	    end;
	end;
	q[i,]=q[i,]/pi[i];
    end;
   


    dcmart=j(numobs,numctime,0);
    do i=1 to numobs;
	do j=1 to numctime;
	    if time[i]>=ctime[j] then dcmart[i,j]=-dcna[j];
	    if time[i]=ctime[j] & cause[i]=0 then dcmart[i,j]=dcmart[i,j]+1;
	end;
	psi[i,]=-dcmart[i,]*q;
    end; 



    sigma=j(&tcov,&tcov,0); 
    do i=1 to numobs;
        sigma = sigma + t(eta[i,]+psi[i,])*(eta[i,]+psi[i,]);
    end;


            
    naivev = inv(fisher);
            
    robustv = inv(fisher)*sigma*inv(fisher);

	se=j(&tcov,1,0);
	do i=1 to &tcov;
		se[i]=sqrt(robustv[i,i]);
	end;
    
	bout=t(b)||se;
	hazout=etime||lambda1;

	create best from bout[colname={'Estimate' 'SE'}];
	append from bout;
	close best; 

	names={%do i=2 %to &numgroup; "G&i" %end; %do i=1 %to &numcov; "Z&i" %end;};
	create covest from robustv[colname=names];
	append from robustv;
	close covest;

	create basehaz from hazout[colname={'Time' 'Base_haz'}];
	append from hazout;
	close basehaz; 


	vut=j(numctime, numtime,0);
	do i=1 to numctime while (cidx[i]>0);
	    vtemp=0;
	    do j=cidx[i] to numtime;
		do k=1 to numobs2;
		    if time2[k]<ctime[i] then vtemp=vtemp+1/s0[j]#dmart1_2[k,j];
		end;
		vut[i,j]=-vtemp;
	    end;
	    vut[i,]=vut[i,]/pi[i];
	end;



/*	variance of direct adjusted CIF  */


    g=j(&numgroup, &numgroup-1,0);
    do i=2 to &numgroup;
        g[i,i-1]=1;
    end;

    outmat=etime;

    zzmat=zzmat01//zzmat2;

    zbar=j(numtime,&tcov,0);
    do i=1 to numtime;
	zbar[i,]=s1[i,]/s0[i];
    end;

        
%do gp=1 %to &numgroup;    

    zz=j(numobs, &numgroup-1+&numcov, 0);
    do 	i=1 to numobs;
	zz[i,]=g[&gp,]||zzmat[i,];
    end;

    dacif=j(numtime,1,0);
    zfn1=j(numtime,1,0);
    zfn2=j(numtime,&numgroup-1+&numcov,0);
    do i=1 to numtime;
	do j=1 to numobs;
	    dacif[i]=dacif[i]+1-exp(-lambda1[i]#exp(b*t(zz[j,])));
	    zfn1[i]=zfn1[i]+exp(-lambda1[i]#exp(b*t(zz[j,]))+b*t(zz[j,]));
	    zfn2[i,]=zfn2[i,]+zz[j,]#exp(-lambda1[i]#exp(b*t(zz[j,]))+b*t(zz[j,]));
	end;
    end;
    dacif=dacif/numobs;
    zfn1=zfn1/numobs;
    zfn2=zfn2/numobs;


    htz=j(numtime,&tcov,0);
    do i=1 to numtime;
	htemp=j(1,&tcov,0);
	do j=1 to i;
	    htemp=htemp+(zfn2[i,]-zfn1[i]#zbar[j,])#dlambda1[j];
 	end;
	htz[i,]=htemp;
    end;



    cifv=j(numtime,1,0);
    cifse=j(numtime,1,0);
    w1temp=j(numobs,1,0);

    do i=1 to numtime;
	w1=j(numobs,1,0);
	w2=j(numobs,1,0);
	w3=j(numobs,1,0);
	do j=1 to numobs;
	    w1temp[j]=w1temp[j]+1/s0[i]#dmart1[j,i];
	    w1[j]=w1temp[j]#zfn1[i];
	    w2[j]=htz[i,]*naivev*t(eta[j,]+psi[j,]);
	    w3[j]=w3[j]+dcmart[j,]*vut[,i]#zfn1[i];

	    tempw=w1[j]+w2[j]+w3[j];
	    cifv[i]=cifv[i]+tempw#tempw;
	end;
    end;

    cifse=cifv##0.5;

    outmat=outmat||dacif||cifse;

%end;


/* 	variance of the difference between direct adjusted CIF's  */

%do gp1=1 %to &numgroup-1;
    %do gp2=&gp1+1 %to &numgroup;

    zz1=j(numobs, &numgroup-1+&numcov, 0);
    zz2=j(numobs, &numgroup-1+&numcov, 0);
    do 	i=1 to numobs;
	zz1[i,]=g[&gp1,]||zzmat[i,];
	zz2[i,]=g[&gp2,]||zzmat[i,];
    end;

    zfn1=j(numtime,1,0);
    zfn2=j(numtime,&numgroup-1+&numcov,0);
    do i=1 to numtime;
	do j=1 to numobs;
	    zfn1[i]=zfn1[i]+exp(-lambda1[i]#exp(b*t(zz1[j,]))+b*t(zz1[j,]))
		-exp(-lambda1[i]#exp(b*t(zz2[j,]))+b*t(zz2[j,]));
	    zfn2[i,]=zfn2[i,]+zz1[j,]#exp(-lambda1[i]#exp(b*t(zz1[j,]))+b*t(zz1[j,]))
		-zz2[j,]#exp(-lambda1[i]#exp(b*t(zz2[j,]))+b*t(zz2[j,]));
	end;
    end;
    zfn1=zfn1/numobs;
    zfn2=zfn2/numobs;


    htz=j(numtime,&tcov,0);
    do i=1 to numtime;
	htemp=j(1,&tcov,0);
	do j=1 to i;
	    htemp=htemp+(zfn2[i,]-zfn1[i]#zbar[j,])#dlambda1[j];
	end;
	htz[i,]=htemp;
    end;	
   

    cif_diff_v=j(numtime,1,0);
    cif_diff_se=j(numtime,1,0);
    w1temp=j(numobs,1,0);
    do i=1 to numtime;
	w1=j(numobs,1,0);
	w2=j(numobs,1,0);
	w3=j(numobs,1,0);
	do j=1 to numobs;
	    w1temp[j]=w1temp[j]+1/s0[i]#dmart1[j,i];
	    w1[j]=w1temp[j]#zfn1[i];
	    w2[j]=htz[i,]*naivev*t(eta[j,]+psi[j,]);
	    w3[j]=w3[j]+dcmart[j,]*vut[,i]#zfn1[i];
	    tempw=w1[j]+w2[j]+w3[j];
	    cif_diff_v[i]=cif_diff_v[i]+tempw#tempw;
	end;
    end;
    cif_diff_se=cif_diff_v##0.5;

    outmat=outmat||cif_diff_se;

    %end;
%end;


    names={'Time' %do i=1 %to &numgroup; "CIF&i" "SE&i" %end;
		%do i=1 %to &numgroup; %do j=&i+1 %to &numgroup; "SE&i._&j" %end; %end;};

    create &outdata from outmat[colname=names];
    append from outmat;
    close &outdata;

quit;

data name;
    %do i=2 %to &numgroup;
    Variable="G&i"; output;
    %end;
    %do i=1 %to &numcov;
    Variable="Z&i"; output;
    %end;
run;

data best;
    merge name best;
    zstat=abs(Estimate/SE);
    Prob=2*probnorm(-zstat);
    drop zstat;
run;

data covest;
    merge name covest;
run;


title 'Cox model for a subdistribution function';
title2 'Estimates of the regression parameters';
proc print data=best; run;

title2 'Estimated variance-covariance matrix';
proc print data=covest; run;


title2 'Estimated baseline cumulative hazard function';
proc print data=basehaz; run;


title2 'Direct adjusted cumulative incidence functions';
proc print data=&outdata; run;


%mend;