/* ---------------------------------------------------------
 program:       PBC0-v92.SAS
 input:         pbc.sas7bdat 
 output:
 does:          Cox regression with stratification
                Dummy vs. True
------------------------------------------------------------ */;

libname data 'H:\bios180\coursepack\SAS-unc-workshop-2008\';

**************************************************;
DATA A;
  set DATA.PBC;
  * create dummy vars for edema levels;
  edema0_5=(edema=0.5);
  edema1=(edema=1);
  trt=(trt1=1);
  run;

**************************************************;

title2 "Measure the treatment effect(Z1)";
proc phreg data=A;
  model time*delta(0)=trt;
run;

* The trial supports not using the drug in this disease;

***********************************************;
* The data on the 312 PBC randomized patients can be ;
* used to build a statistical model for the influence;
* of covariates on disease outcome;

title2 "A two-group comparison of the gender groups(Z3)";
* Score test is logrank test when there is no ties;

proc phreg data=A;
  model time*delta(0)=sex;
run;

*************************************************;
title2 'Look at the presence of edema(Z7)';

proc phreg data=A;
  model time*delta(0)=edema;
run;

***********************************************;
* Now factor edema, and perform the score and LR tests
* of heterogeneity. Is there strong evidence for  ;
* departure from linear trend?  ;

title2 'Edema factored';
proc phreg data=A;
  model time*delta(0)=edema0_5 edema1;
run;

title2 'Likelihood Ratio test of edema linear vs. factored';
* cut and past LR, compute p-value;
data LRtest;
   LR = (1223.840 - 1222.391);
   pvalue=1-probchi(LR,1);
run;
proc print data=LRtest; run;

***********************************************;
title2 'Look at the continuous covariate age(Z2)';

proc phreg data=A;
  model time*delta(0)=age;
run;

* What is the interpretation of the risk ratio?   ;

***********************************************;
title2 'Estimating Survival Function and Cumulative Hazard Function';
data D;
  input sex age edema trt bili;
  cards;
  0 0 0 0 0
        ;
run;

* Estimate survival function and cumulative hazard function; 
* at 0 of the covariates;
proc phreg data=A ;
  model time*delta(0)= sex age edema trt bili;
  baseline out=B1  covariates=D survival=s
                logsurv=ls loglogs=lls cumhaz=chaz;
run;


* Estimate survival function and cumulative hazard function;
* at the mean value of the each covarite;
proc phreg data=A ;
  model time*delta(0)= sex age edema trt bili;
  baseline out=B2 survival=s
                logsurv=ls loglogs=lls cumhaz=chaz;
run;


*combine data sets for the plots;
data B;
 set B1 B2;
  if age = 0 then type='at 0 covars';
            else type='at avg.covars';
run;


title1 h=1.5 "Cox regression for PBC data";
title2 'Estimated Cumulative Hazard Function';
proc sgplot data=B;
step y=chaz x=time /group=type;
run;
quit;

title1 h=1.5 'Cox regression for  PBC data';
title2 'Estimated Survival Function';
proc sgplot data=B;
step y=s x=time /group=type;
run;
quit;