/* The macro omvssh2 is used by the macro CumIncV (below) */

%macro omvssh2(Data,Time,nStrata,obsnr);

 * udvælger data til beregning af cov2(t);
 data retur;     
  set &Data;
     n=_N_;
     if (n gt &obsnr) then delete;
    %do j=0 %to &nStrata; 
     omvp0&j=0;
    %end;
  run;

 %do j=1 %to &nStrata;
 data retur;
  set retur end=last;
   if last then call symput('p',p0&j);

 data retur;
  set retur;
   lagp0=lag(p00);
   if (lagp0 eq 0) then omvP0&j=.;
    if (lagp0 ne 0) then omvP0&j=(&p-p0&j)/lagp0;
   if (&Time eq 0) then omvP0&j=&p;
*   drop lagp0;
 %end;

 data retur;
  set retur;
  omvp00=1;
  %do j=1 %to &nStrata;
   omvp00=omvp00-omvp0&j;
  %end;
run;

%mend omvssh2;








%macro 
 CumIncV(StackDat,Data,Strata,Time,D,Surv,nCov,ResData,zbeta,EstData);


/* The number of strata (nStrata) */


proc sort data=&Data;
 by &Strata;

data nstrat;
 set &Data end=last;
  by &Strata;
  firstS=first.&Strata;
  retain nStrata 0;
   nStrata+firstS;
  if last then call symput('nStrata',nStrata);
  drop firstS;
run;


/* Observations in Data are deleted */


data newData;
 set &Data;
  start=(&Time eq 0);
  retain komb 0;
   komb+start;
  med=0;
 %do k=0 %to (&nStrata-1) %by 1;
   med=med+(komb eq (1+&k*(&nStrata+1)));
 %end;
 if (med eq 1);
 drop start med komb;




/* The number of jumps (definition of nFailure) */

proc sort data=&StackDat;
 by &Strata &Time;

proc means data=&StackDat sum noprint;
 var &D;
 by &Strata &Time;
 output out=nFData Sum=nFailure;

data nFData (keep=&Time &Strata nFailure);
 set nFData;
  if nFailure ne 0;
  
proc sort data=nFData;
 by &Strata &Time;


data NyData;
 merge newData nFdata;
 by &Strata &Time;
 if (nFailure eq .) then nFailure=0;



/* The number of observations per time point */

 proc freq data=NyData noprint;
  tables &Time / out=nspring;

 proc sort data=nspring;
  by count;


 data maxcount;
  set nspring end=last;
   by count;
   if first.count then do;
    n+count-n;
    output;
   end;  
  if last then call symput('antal',n);
run;



/* Include strata vector in NyData -> strata */

 proc sort data=NyData;
  by &Strata &Time;

 data strata;
  set NyData; by &Strata;
  retain stratum 0;
  stratum+first.&Strata;



/* Time vector */

 proc sort data=NyData;
  by &Time;

 data Time (keep=&Time); 
  set NyData;
 %do i=1 %to &antal %by 1;
  if &Time=lag(&Time) then delete;
 %end;


/***********************************************************************/
/************************ The probabilities ****************************/


/* Calculate Ai */

 %do i=1 %to &nStrata %by 1;
  data data&i;
   set strata;
   if stratum=&i;

  data data&i (keep = &Time A&i);
   merge data&i Time; by &Time;  
   retain Surv&i;
   if not (&Surv=.) then Surv&i=&Surv;
   A&i=-log(Surv&i);

  drop &Surv;

 %end;


/* A=sum(A_1,...,A_nStrata) and S */

data A_and_S;
  set data1;
  A=A1;
  %do i=2 %to &nStrata %by 1;
   merge data&i; by &Time;
   A+A&i;
 %end;

 dA=A-lag(A);
 if dA eq . then dA=0;

 retain S 1;
  S+S*(1-dA)-S;
 
 lagS=lag(S);

 drop A S dA;

run;


/* merge with data strata -> data  */

 proc sort data=strata; by &Time;

 data data;
  merge strata A_and_S; by &Time;

 proc sort data=data; 
  by &Strata &Time;


/* Determine p's */

 data data; 
  set data; by &Strata;

   if first.&Strata then do;
    lagS=1; end;
   
   %do i=1 %to &nStrata %by 1;
    lagA&i=lag(A&i);
    if ((stratum eq &i) and (first.&Strata)) then lagA&i=0;
    if (stratum ne &i) then lagA&i=0;

    l&i=(stratum=&i)*(lagS*(A&i-lagA&i));
    retain p&i 0;
    p&i+l&i; 
    if (stratum ne &i) then p&i=0;
    drop l&i;
   %end;

    p=0;
  %do i=1 %to &nStrata %by 1;
    p=p+p&i; 
    drop p&i;
  %end;
   drop lagS;
run;

proc sort data=data;
 by &Strata &Time;
run;



/*** S^0, S^1, E and the residuals ***/


 /* Make strata vector in &Resdata */
  proc sort data=&Resdata;
   by &Strata &Time;

  data Resdata;
   set &Resdata; by &Strata;
   retain stratum 0;
   stratum+first.&Strata;
  run;


 %do i=1 %to &nStrata;

  data Resdat&i;
   set Resdata;
    if stratum=&i;
    drop stratum &Strata;
  run;


 /* Unique vector of time points */
 /* The number of observations for each timepoint */

  proc freq data=Resdat&i noprint;
   tables &Time / out=nspring;

  proc sort data=nspring;
   by count;

  data maxcount;
   set nspring end=last;
    by count;
    if first.count then do;
     n+count-n;
     output;
    end;  
   if last then call symput('resantal',n);
 run;


 /* merge Resdat&i and Time */
 data Resdat&i;
  set Resdat&i;
  by &Time;
   l=last.&Time;
  run;

 proc sort data=Resdat&i;
  by descending &Time;
 run;

 /* S^0 and E (when covariates are present) */
 data Resdat&i;
  set Resdat&i; 

    if (&nCov eq 0) then &zbeta=0;
    b0=exp(&zbeta);
    retain S0_&i 0;
    S0_&i+b0;
    drop b0;
 /* */
    %if (&nCov ge 1) %then %do;
     %do j=1 %to &nCov;
       b1_&j=cov&j * exp(&zbeta);
     %end;
     %do j=1 %to &nCov;
      retain S1_&i&j 0;
      S1_&i&j+b1_&j;
      drop b1_&j;
     %end;
     %do j=1 %to &nCov;
      if (S0_&i eq 0) then E_&i&j=.;
      if (S0_&i ne 0) then E_&i&j=S1_&i&j/S0_&i;
      drop S1_&i&j;
     %end;
    %end;
 run;
 /* */

 /* define time vector */
 data Resdat&i;
  set Resdat&i;
  if (l eq 1);
  %do j=1 %to &nCov;  
   drop Cov&j;
  %end;
  drop l;
 run;

 proc sort data=Resdat&i;
  by &Time;
 run;

 %end; * end of i-loop;


/* Determine p's, Ss, jump processes, covariates on strata */


%do i=1 %to &nStrata;
 data data&i;
  set data;
   if stratum=&i;
   p0&i=p;
   n0&i=nFailure;
    zbeta&i=&zbeta;
  drop nFailure stratum &Strata &zbeta;
 run;

 proc sort data=data&i;
  by &Time;
  run;




 data data&i;
  merge data&i Resdat&i;
   by &Time;  
 /* */
   %if (&nCov ge 1) %then %do;
    %do j=1 %to &nCov;
     Res_&i&j=Cov&j-E_&i&j;
    drop Cov&j E_&i&j;
   %end;
  %end;
 /* */
run;

%end; * end of i-loop;


 /* collect the residuals and Es */
 data vardata;
  set data1;
 run;
 %do i=2 %to &nStrata;
 data vardata;
   merge vardata data&i; by &Time;
 %end;
 run;


 /* complete the probabilities */
 %do i=1 %to &nStrata;

  data antmgl&i;
   set vardata;
    mgl=(p0&i eq .);
    if (mgl eq 0) then antmgl=0; 
      else antmgl+1; 
 
  proc sort data=antmgl&i;
   by antmgl;

  data antmgl&i;
   set antmgl&i end=last;
    by antmgl;
     if first.antmgl then do;
      n+1;
      output;
     end;  
     if last then call symput('max',n);
  run;

  data vardata;
   set vardata;
    %do j=1 %to &max;
     dummy=lag(p0&i);
     if (p0&i eq .) then p0&i=dummy;
    %end;
     drop dummy;
  run;

 %end;

 /* make p00 */
 data vardata;
  set vardata;
   p00=1;
   %do i=1 %to &nStrata;
    p00=p00-p0&i; 
   %end; 

  proc sort data=vardata;
   by &Time;
  run;

 /* delete redundant variables */
  data vardata; 
   set vardata;
   %do i=1 %to &nStrata;
     drop A&i lagA&i;
   %end;
run;




/*************************************************************************/
/************************* THE VARIANCES *********************************/



%if (&nCov ge 1) %then %do; 

/* define Ws */

data vardata;
 set vardata;
  %do i=1 %to &nStrata %by 1;      * define terms for W0i;
   %do j=1 %to &nCov %by 1;         
     lagb=exp(zbeta&i)*Res_&i&j/S0_&i;   * Res and S0 are left cont.;
     b0&i&j=lagb*N0&i;
     if (b0&i&j eq .) then b0&i&j=0;
    %end;
  %end;

    %do i=1 %to &nStrata;   * sum the terms -> W0ij;
   %do j=1 %to &nCov;
     retain W0&i&j 0;
     W0&i&j+b0&i&j;
   %end;
  %end;

  %do j=1 %to &nCov;     * W00j;
   W00&j=0;
   %do i=1 %to &nStrata;
     W00&j=W00&j-W0&i&j;
   %end;
  %end;

%do i=1 %to &nStrata;   * delete the redundant variables;
    %do j=1 %to &nCov;
     drop b0&i&j Res_&i&j;
    %end;    
%end;

 drop lagb;
run;


%end;




/****** make cov1 and cov2 ******/

/* define the terms of cov1 and cov2, int&j&i */
data vardata;
 set vardata;
  %do j=0 %to &nStrata;
   cov2_0&j=0;
   cov1_0&j=0;   * unchanged when no covariates are present;
  %end;

 %if (&nCov ge 1) %then %do;
   %do j=0 %to &nStrata;
    %do i=1 %to &nCov;
     int&j&i=0;
    %end;
   %end;
  %end; 

run;



/* delete the time points for which only censorings are observed */
data vardata;
 set vardata;
 ind=0;
 %do i=1 %to &nStrata;
    ind=ind+((n0&i eq 0) or (n0&i eq .));
 %end;
 if (&Time eq 0) then ind=0;   * to make sure time point 0 is not deleted;
 if ind=&nStrata then delete;
 drop ind;
run;


/* find the number of observations in vardata */
proc sort data=vardata;
 by &Time;

data nobs;
 set vardata end=last;
  by &Time;
   if first.&Time then do;
    n+1;
    output;
   end;
  if last then call symput('nobs',n);
run;

%do t=2 %to &nobs;


%omvssh2(vardata,&Time,&nStrata,&t);


/***  cov2 is made  ***/

data retur;
 set retur;
 %do j=0 %to &nStrata; 
  %do l=1 %to &nStrata;
    lagb=(lag(p00)**2) * (((&l eq &j)-lag(omvp0&j))**2) * 
	  (exp(zbeta&l) / S0_&l)**2;
    b&l&j=lagb*N0&l;
     if (b&l&j eq .) then b&l&j=0;
  %end;
 %end;

  drop lagb;

* sum over strata l;
%do j=0 %to &nStrata;
b&j=0; 
   %do l=1 %to &nStrata;
     b&j+b&l&j;
   %end;    
%end;

* sum over time;
%do j=0 %to &nStrata;
     retain cov2_&j 0;
     cov2_&j+b&j;
%end;

%do j=0 %to &nStrata;
    %do l=1 %to &nStrata;
     drop b&l&j b&j;
    %end;    
%end;

run;



/* cov2 at time point t */
%do j=0 %to &nStrata;
data retur;
 set retur end=last;
 %global cov02&j;
 if last then call symput("cov02&j",cov2_&j);
%end;


/* */
%if (&nCov ge 1) %then %do; 
/*** the integral for cov1 ***/

%do j=0 %to &nStrata;			* lag(P00) * d(W00) * P0j;
 %do i=1 %to &nCov;       
  lagW00&i=lag(W00&i);
  dW00&i=W00&i-lagW00&i;
  if dW00&i eq . then dW00&i=0;

  lagp00=lag(p00);
   if lagp00 eq . then lagp00=0;

  b1=lagp00*dW00&i*omvp0&j;

  retain int1&j&i 0;			* integral of type 1;
  int1&j&i+b1;



  %if (&j ne 0) %then %do;
   lagW0&j&i=lag(W0&j&i);
   dW0&j&i=W0&j&i-lagW0&j&i;

   b2&j&i=lagp00*dW0&j&i;
 
   retain int2&j&i 0;			* integral of type 2;
   int2&j&i+b2&j&i;
   drop b2&j&i;
  %end;

  %if (&j ne 0) %then %do; int&j&i=int1&j&i+int2&j&i; %end; *the integral;
  %if (&j eq 0) %then %do; int&j&i=int1&j&i; %end;

 %end; 
%end;

 drop b1;


/* int&j&i at time point t */

%do j=0 %to &nStrata;
 %do i=1 %to &nCov;
  %global int0&j&i;
  if last then call symput("int0&j&i",int&j&i);
 %end;
%end;

 %end;
/* */

/* cov2 at time point t is written in to vardata */

proc iml;
 edit vardata;
  %do j=0 %to &nStrata;
    cov2_0&j=&&cov02&j;
    replace point &t var {cov2_0&j};
  %end;
/* */
  %if (&nCov ge 1) %then %do;
  %do j=0 %to &nStrata;
   %do i=1 %to &nCov;
    int&j&i=&&int0&j&i;
    replace point &t var {int&j&i};;
   %end;
  %end;
  %end;
/* */
 close vardata;
quit;





%end; * end of t-loop;

/* */
%if (&nCov ge 1) %then %do;
/****** cov1 is made******/

/* cov1 */

proc iml;
 
 use &Estdata; 		* F_T_b;
 read all into var1;
 F_t_b=var1[2:(&ncov+1),1:&ncov];

 use vardata; 		* vector of integrals;
 %do i=0 %to &nStrata;
   %do j=1 %to &nCov;
    read all var{int&i&j};
   %end;

    int&i=int&i.1; run;

   %do j=2 %to &nCov;
    int&i=int&i||int&i&j;
   %end;
 %end;


%do i=0 %to &nStrata;   * cov1 is calculated;
 c1_0&i=repeat(0,nrow(int0),1); 
  do n=1 to nrow(int0);  
   vec=int&i[n,];
   c1_0&i[n,1]=vec * F_T_b * vec`;  
  end;
%end;

edit vardata;		* cov1 is saved;
 %do i=0 %to &nStrata;
   do n=1 to nrow(int0);
     cov1_0&i=c1_0&i[n,1];
     replace point n var {cov1_0&i};
   end;
 %end;
close vardata;

quit; 
%end;
/* */

/* sum cov1 and cov2 and delete redundant variables */

data data;
 set vardata;
 %do i=1 %to &nStrata;
  Var0&i=cov1_0&i+cov2_0&i;
 %end;
  Var00=cov1_00+cov2_00;
  keep &Time;
 %do i=0 %to &nStrata;
  keep P0&i var0&i;
 %end;
run;


%mend CumIncV;


*options nonotes;
run;