/* C.1 Tumorgenicity Data Analysis of Chapter 3 */

options nocenter nodate ls=80;

data RATS;
  infile "rats.dat";
  input id start stop status enum trt;
run;

proc sort data=RATS; by id enum; run;


/* C.1.2 Poisson Analysis with Piecewise Constant Rate */

data RATS_PW_TMP;
  set RATS;
  by id;
  stop = floor(stop);

  retain count1 .;
  retain count2 .;
  retain count3 .;
  retain count4 .;
  if (first.id) then do;
     count1 = 0; count2 = 0; count3 = 0; count4 = 0;
     if ( (stop >  0) & (stop <=  30) & (status = 1) ) then count1 = 1;
     if ( (stop > 30) & (stop <=  60) & (status = 1) ) then count2 = 1;
     if ( (stop > 60) & (stop <=  90) & (status = 1) ) then count3 = 1;
     if ( (stop > 90) & (stop <= 122) & (status = 1) ) then count4 = 1;
  end;
  else do;
     if ( (stop >  0) & (stop <=  30) & (status = 1) ) then count1 = count1 + 1;
     if ( (stop > 30) & (stop <=  60) & (status = 1) ) then count2 = count2 + 1;
     if ( (stop > 60) & (stop <=  90) & (status = 1) ) then count3 = count3 + 1;
     if ( (stop > 90) & (stop <= 122) & (status = 1) ) then count4 = count4 + 1;
  end;

  len1 = 30; len2 = 30; len3 = 30; len4 = 32;
  interval1 = 1; interval2 = 2; interval3 = 3; interval4 = 4;
  if (last.id);
run;

data RATS_PW;
  set RATS_PW_TMP(keep = id interval1 count1 len1 trt
                  rename=(interval1=interval count1=count len1=len))
      RATS_PW_TMP(keep = id interval2 count2 len2 trt
                  rename=(interval2=interval count2=count len2=len))
      RATS_PW_TMP(keep = id interval3 count3 len3 trt
                  rename=(interval3=interval count3=count len3=len))
      RATS_PW_TMP(keep = id interval4 count4 len4 trt
                  rename=(interval4=interval count4=count len4=len));
  loglen = log(len);
run;

proc sort data=RATS_PW; by id interval; run;


/* Fit a piecewise constant model for control rats */
proc genmod data=RATS_PW(where = (trt = 0));
  class interval / param=REF ref=FIRST;
  model count = interval / offset=loglen dist=poisson link=log covb;
run;

/* Fit a multiplicative model to estimate the treatment effect */
proc genmod data=RATS_PW;
  class interval trt / param=REF ref=FIRST;
  model count = interval trt / offset=loglen dist=poisson link=log covb;
run;


/* C.1.3 Nonparametric Analysis */

/* Nelson-Aalen estimates and SE for control rats */
proc phreg data=RATS(where = (trt = 0));
  model (start,stop)*status(0)= / ties=breslow;
  baseline out=NADATA cumhaz=na_mf stdcumhaz=se_na_mf;
run;

proc print data=NADATA; 
  var stop na_mf se_na_mf;  
run;


/* Cumulative mean function estimates, and Poisson-based and robust SE for control rats */
proc phreg data=RATS(where = (trt = 0)) covs(aggregate);
  model (start,stop)*status(0)= / ties=breslow;
  id id;
  baseline out=CMFDATA cmf=cmf stdcumhaz=cmf_se_poisson stdcmf=cmf_se_robust;
run;

proc print data=CMFDATA; 
  var stop cmf cmf_se_poisson cmf_se_robust;
run;
 

/* C.1.3 Semiparametric Poisson Analysis */

proc phreg data=RATS;
  model (start,stop)*status(0)= trt / ties=breslow;
run;


/* C.1.5 Robust Semiparamtric Analysis */ 

proc phreg data=RATS covs(aggregate);
  model (start,stop)*status(0)=trt / ties=breslow;
  id id;
run;