/* PubH 7450
   Exampel 9.4b:  Using an alternative specification in PHREG to
                  fit a PHM for the left-truncated and right-censored 
                  Channing House data on p.313; see also Fig 4.11 on p.125
   Reason: SAS manual says that "Although left truncation can be 
           accommodated in PROC PHREG through the counting process style 
           of input, such specification does not allow survival estimates 
           to be output. Using the ENTRY= option in PROC PHREG for left 
           truncation does not suppress computing the survival estimates.", 
           which however is probably incorrect, BUT ...
*/

options ls=80 center pagesize=60 number label;
goptions device = PS;
         /* graphs are saved in a file names sasgraph.ps */

data channing;
infile '/home/merganser/weip/public_html/course/7450/data/channing.dat' 
       firstobs=2;
input Obs status Age1 Age2 Diff Gender;

/* 1) Obtain a nonparametric (default) K-M estimate of the survival 
      for males: */
data males;
  set channing;
  if Gender=1;

proc phreg data=males ;
model Age2*status(0) = / entry=Age1;
output out=a survival=S num_left=Y ;

/*
proc print data=a;
*/
proc sort data=a;
  by Age2;
data Pred0;
  set a;
  Pattern=1;
  keep Age2 S Pattern;

/* 2) Obtain a semi-parametric PHM-based estimate of the survival for males: */
data malecov;
  input Gender;
  datalines;
1
;

proc phreg data=channing ;
model Age2*status(0) = Gender / entry=Age1;
output out=a survival=surv num_left=Y ;
baseline covariates=malecov out =Pred1 survival=S 
         / nomean;

/*
proc print data=Pred1;
*/
data Pred1;
  set Pred1;
  Pattern=2;
  keep Age2 S Pattern;

/* merge the two datasets to faciliatte drawing a plot */
data both;
  set Pred0 Pred1;
run;

   legend1 label=none shape=symbol(4, .8)
      value=(f=swiss h=.8 'Nonparametric' 'PHM-based');
   axis1 label=(h=1 f=swiss a=90) minor=(n=1);
   axis2 label=(h=1 f=swiss 'Age in Months') minor=(n=4);

   proc gplot data=both;
      plot S*Age2=Pattern / legend=legend1 vaxis=axis1
                            haxis=axis2  ;
      symbol1 interpol=stepLJ h=1 v=square  c=blue;
      symbol2 interpol=stepLJ h=1 v=diamond c=red;
      note f=swiss h=1.5 j=c 'Survival of male patients';
   run;

/* 3): check the PH assumption. */
proc phreg data=channing ;
model Age2*status(0) = Gender gt / entry=Age1;
gt = Gender*log(Age2);
run;
************************************1):

                              -2 LOG L = 278.2255

************************************2):
                             Model Fit Statistics
 
                                     Without           With
                    Criterion     Covariates     Covariates

                    -2 LOG L        1607.597       1604.433
                    AIC             1607.597       1606.433
                    SBC             1607.597       1609.604

                    Analysis of Maximum Likelihood Estimates
 
                   Parameter     Standard                               Hazard
  Variable   DF     Estimate        Error   Chi-Square   Pr > ChiSq      Ratio

  Gender      1     -0.31596      0.17313       3.3305       0.0680      0.729

************************************3):
                             Model Fit Statistics
 
                                     Without           With
                    Criterion     Covariates     Covariates

                    -2 LOG L        1607.597       1601.727
                    AIC             1607.597       1605.727
                    SBC             1607.597       1612.068

                    Analysis of Maximum Likelihood Estimates
 
                   Parameter     Standard                               Hazard
  Variable   DF     Estimate        Error   Chi-Square   Pr > ChiSq      Ratio

  Gender      1    -26.90095     15.98891       2.8307       0.0925      0.000
  gt          1      3.85515      2.32008       2.7611       0.0966     47.236