options nocenter nodate ls=80;

title 'rhDNase Data';


/* 1. Here we read data from an existing ASCII file and create          */
/*    data sets for the first and second gap time analysis              */
data rhDNase_etype1;
  infile "c:/WORKSHOP/rhDNase.dat";
  input id trt fev fev2 time1 time2 status etype enum enum1 enum2;
  gtime = time2 - time1;
  fevc  = fev - 61.07783;  /* mean(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;  /* mean(gtime1) = 70.96916 */
  lgtime1 = log(gtime1);
  drop status1;
  if (inrhDNase2);
run;


/* 2. Cox Analysis of First and Second Gap Times                        */
/* 2(A). Here we fit a Cox model to first gap times                     */
proc phreg data=rhDNase1;
  model gtime*status(0) = trt fevc / ties=breslow;
run;

/* Here we test the proportion hazards assumption                        */
/* This is similar to the cox.zph(, transform="log") functionin R/S-Plus */
/* - The test statement is used to test all 'time-dependent' covariates. */
proc phreg data=rhDNase1;
  model gtime*status(0) = trt fevc trtt fevct / ties=breslow;
  trtt  = trt*log(gtime);
  fevct = fevc*log(gtime);
  proportionality_test_global: test trtt, fevct;
run;


/* 2(B). Here we fit a Cox model to the second gap times                 */
proc phreg data=rhDNase2;
  model gtime*status(0) = trt fevc gtime1c / ties=breslow;
run;
  

/* 3. Parametric Log-normal Model for First and Second Gap Times         */ 

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

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




/* 4. Fitting Modulated Markov Models                                    */

/* Model 1 */
proc phreg data=rhDNase_etype1;
  model (time1,time2)*status(0) = trt fevc
                                  / ties=breslow fconv=1e-04;
run;


/* Model 4 */
data rhDNaseTD;
  infile "c:/WORKSHOP/rhDNasetd80.dat";
  input id time1 time2 status enum trt fevc It Nt;

  trt_It0 = 0.0;
  trt_It1 = 0.0;
  if (It = 1) then trt_It0 = trt;
  if (It = 0) then trt_It1 = trt;
run;

proc sort data=rhDNaseTD; by id enum; run;
proc print data=rhDNaseTD(obs=18); run;

proc phreg data=rhDNaseTD;
  model (time1,time2)*status(0) = trt_It0 trt_It1 fevc Nt
                                  / ties=breslow fconv=1e-04;
run;