/* C.2 Code for rhDNase Data Analyses */

options nocenter nodate ls=80;


/* Create data for fist and second gap time analysis */
data rhDNase_etype1;
  infile "rhDNase.dat";
  input id trt fev fev2 time1 time2 status etype enum enum1 enum2;
  gtime = time2 - time1;
  fevc  = fev - 61.07783;
  if (etype = 1);
run;

proc sort data=rhDNase_etype1; by id enum1; run;

data rhDNase1; 
  set rhDNase_etype1(where =(enum1 = 1)); 
run;

data rhDNase2; 
  merge rhDNase_etype1(where =(enum1 = 2) in=inrhDNase2)
        rhDNase1(keep=id gtime status
                 rename=(gtime=gtime1 status=status1)
                 where=(status1 = 1));
  by id;
  gtime1c = gtime1 - 70.96916;
  lgtime1 = log(gtime1);
  drop status1;
  if (inrhDNase2);
run;


/* Fit a Cox model to first gap times */
proc phreg data=rhDNase1;
  model gtime*status(0) = trt fevc / ties=breslow;
run;


/* Test for proportion hazards assumption */
/* Similar to cox.zph(, transform="log")  */
proc phreg data=rhDNase1;
  model gtime*status(0) = trt fevc trtt fevct / ties=breslow;
  trtt  = trt*log(gtime);
  fevct = fevc*log(gtime);
  proportionality_test_trt: test trtt;
  proportionality_test_fevc: test fevct;
  proportionality_test_global: test trtt, fevct;
run;


/*  Fit a Cox model to second gap times */
proc phreg data=rhDNase2;
  model gtime*status(0) = trt fevc gtime1c / ties=breslow;
run;


/* Fit a parametric log-normal model to first gap time */
proc lifereg data=rhDNase1;
  model gtime*status(0) = trt fevc / dist=lnormal;
run;


/* Fit a parametric log-normal model to second gap time */
proc lifereg data=rhDNase2;
  model gtime*status(0) = trt fevc lgtime1 / dist=lnormal;
run;