options nocenter nodate ls=80;

title 'rhDNase Data';


/* 1. Read data from an existing ASCII file                               */
data rhDNase;
  infile "c:/WORKSHOP/rhDNase.dat";
  input id trt fev fev2 time1 time2 status etype enum enum1 enum2;
  fevc  = fev - 61.07783;  /* mean(fev) = 61.07783 */
  gtime = time2 - time1;

  status1 = 0.0;
  if (etype = 1) then status1 = status;
run;

proc sort data=rhDNase; by id time1 time2; run;

proc print data=rhDNase(obs=10); 
  var id enum etype time1 time2 gtime status status1 trt fev fevc;
run;


/* 1. Nonparametric Estimates for Control Arm */

/* 1(A). Data includes patients in risk set when they have exacerbations  */
proc phreg data=rhDNase(where=(trt = 0)) covs(aggregate);
  model (time1, time2)*status1(0) = / ties=breslow;
  id id;
  baseline out=CMFDATA_ALL
           cmf=cmf stdcumhaz=cmf_se_poisson stdcmf=cmf_se_robust;
run;

proc print data=CMFDATA_ALL(obs=10);
  var time2 cmf cmf_se_poisson cmf_se_robust;
run;


/* 1(B). Data excludes patients from risk set when they have exacerbations */
proc phreg data=rhDNase(where=((etype = 1) & (trt = 0))) covs(aggregate);
  model (time1, time2)*status1(0) = / ties=breslow;
  id id;
  baseline out=CMFDATA_etype1
           cmf=cmf stdcumhaz=cmf_se_poisson stdcmf=cmf_se_robust;
run;

proc print data=CMFDATA_etype1(obs=10);
  var time2 cmf cmf_se_poisson cmf_se_robust;
run;


/* 2. Andersen-Gill Analysis */

/* 2(A). Data includes patients in risk set when they have exacerbations  */
proc phreg data=rhDNase;
  model gtime*status1(0) = trt fevc / ties=breslow;
run;


/* 2(B). Data excludes patients from risk set when they have exacerbations */
proc phreg data=rhDNase(where=(etype = 1));
  model gtime*status1(0) = trt fevc / ties=breslow;
run;