Chapter Contents
Chapter Contents
Previous
Previous
Next
Next
The PHREG Procedure

Example 49.8: Multiple Failure Outcomes

For survival data with multiple failure outcomes for each subject, the failures may be repetitions of the same kind of event or they may be events of different nature.

The Andersen-Gill (AG) model and the Prentice, Williams, and Peterson (1981) models, also referred to as the PWP models, can be used in the analysis of repeated failure outcomes of the same kind, while the marginal analysis approach of Wei, Lin, and Weissfeld (1989), also referred to as the WLW analysis, can be applied to both multiple events of the same types and multiple events of different types. The bladder cancer data listed in Wei, Lin, and Weissfeld (1989) are used to illustrate these methods of analyses.

The data consist of 86 patients with superficial bladder tumors, which were removed when they entered the study. Of these patients, 48 were randomized into the placebo group, and 38 were randomized into the thiotepa group. Many patients had multiple recurrences of tumors during the study, and new tumors were removed at each visit. The data set contains the first four recurrences of the tumor for each patient, and each recurrence time was measured from the patient's entry time into the study.

The input data consist of the following eight variables:

In the data set Bladder, four observations are created for each patient, each corresponding to one of the four tumor recurrences. In addition to values of Trt, Number, and Size for the patient, each observation contains the following variables:

For instance, a patient with only one recurrence time at month 6, who was followed until month 10, will have values for Visit, TStart, TStop, and Status of (1,0,6,1), (2,6,10,0), (3,10,10,0), and (4,10,10,0). If the follow-up time of a patient is beyond the time of the fourth tumor recurrence, an extra observation is created to represent this last follow-up period. For instance, a patient with recurrences at months 2, 15, 34, and 50, who was followed until month 61, will have five observations with values for Visit, TStart, TStop, and Status of (1,0,2,1), (2,2,15,1), (3,15,34,1), (4,34,50,1), and (5,50,61,0). In the former situation, the last two observations are redundant for the AG model, but they are important for the WLW analysis. In the latter situation, the fifth observation is not needed for the WLW model, but it is indispensable to the AG analysis.

The following SAS statements create the data set Bladder:

   data Bladder;
      keep ID TStart TStop Status Trt Number Size Visit;
      retain ID TStart 0;
      array tt T1-T4;
      infile datalines missover;
      input Trt Time Number Size T1-T4;
      ID + 1;
      TStart=0;
      do over tt;
         Visit=_i_;
         if tt = . then do;
            TStop=Time;
            Status=0;
         end;
         else do;
            TStop=tt;
            Status=1;
         end;
         output;
         TStart=TStop;
      end;
      if (TStart < Time) then do;
         TStop= Time;
         Status=0;
         Visit=5;
         output;
      end;
      datalines;
   1       0       1     1
   1       1       1     3
   1       4       2     1
   1       7       1     1
   1       10      5     1
   1       10      4     1   6
   1       14      1     1
   1       18      1     1
   1       18      1     3   5
   1       18      1     1   12  16
   1       23      3     3
   1       23      1     3   10  15
   1       23      1     1   3   16  23
   1       23      3     1   3   9   21
   1       24      2     3   7   10  16  24
   1       25      1     1   3   15  25
   1       26      1     2
   1       26      8     1   1
   1       26      1     4   2   26
   1       28      1     2   25
   1       29      1     4
   1       29      1     2
   1       29      4     1
   1       30      1     6   28  30
   1       30      1     5   2   17  22
   1       30      2     1   3   6   8   12
   1       31      1     3   12  15  24
   1       32      1     2
   1       34      2     1
   1       36      2     1
   1       36      3     1   29
   1       37      1     2
   1       40      4     1   9   17  22  24
   1       40      5     1   16  19  23  29
   1       41      1     2
   1       43      1     1   3
   1       43      2     6   6
   1       44      2     1   3   6   9
   1       45      1     1   9   11  20  26
   1       48      1     1   18
   1       49      1     3
   1       51      3     1   35
   1       53      1     7   17
   1       53      3     1   3   15  46  51
   1       59      1     1
   1       61      3     2   2   15  24  30
   1       64      1     3   5   14  19  27
   1       64      2     3   2   8   12  13
   2       1       1     3
   2       1       1     1
   2       5       8     1   5
   2       9       1     2
   2       10      1     1
   2       13      1     1
   2       14      2     6   3
   2       17      5     3   1   3   5   7
   2       18      5     1
   2       18      1     3   17
   2       19      5     1   2
   2       21      1     1   17  19
   2       22      1     1
   2       25      1     3
   2       25      1     5
   2       25      1     1
   2       26      1     1   6   12  13
   2       27      1     1   6
   2       29      2     1   2
   2       36      8     3   26  35
   2       38      1     1
   2       39      1     1   22  23  27  32
   2       39      6     1   4   16  23  27
   2       40      3     1   24  26  29  40
   2       41      3     2
   2       41      1     1
   2       43      1     1   1   27
   2       44      1     1
   2       44      6     1   2   20  23  27
   2       45      1     2
   2       46      1     4   2
   2       46      1     4
   2       49      3     3
   2       50      1     1
   2       50      4     1   4   24  47
   2       54      3     4
   2       54      2     1   38
   2       59      1     3
   ;

The counting process MODEL specification is used to carry out the analysis of the AG model. Note that some of the observations in the data set Bladder have a degenerated interval of risk. The presence of these observations does not affect the results of the analysis since none of these observations are included in any of the risk sets. However, the procedure will run more efficiently without these observations; consequently, in the following SAS statements, the WHERE clause is used to eliminate these redundant observations.

   *title 'Andersen-Gill Multiplicative Hazards Model';
   proc phreg data=Bladder;
      model (TStart, TStop) * Status(0) = Trt Number Size;
      where TStart < TStop;
   run;

Results of fitting the AG model are shown in Output 49.8.1.

Output 49.8.1: Fitting Andersen-Gill Multiplicative Hazards Model

The PHREG Procedure

Model Information
Data Set WORK.BLADDER
Dependent Variable TStart
Dependent Variable TStop
Censoring Variable Status
Censoring Value(s) 0
Ties Handling BRESLOW

Summary of the Number of Event and Censored
Values
Total Event Censored Percent
Censored
190 112 78 41.05

Convergence Status
Convergence criterion (GCONV=1E-8) satisfied.

Model Fit Statistics
Criterion Without
Covariates
With
Covariates
-2 LOG L 934.210 920.159
AIC 934.210 926.159
SBC 934.210 934.315

Testing Global Null Hypothesis: BETA=0
Test Chi-Square DF Pr > ChiSq
Likelihood Ratio 14.0509 3 0.0028
Score 15.4173 3 0.0015
Wald 15.1736 3 0.0017

Analysis of Maximum Likelihood Estimates
Variable DF Parameter
Estimate
Standard
Error
Chi-Square Pr > ChiSq Hazard
Ratio
Trt 1 -0.40710 0.20007 4.1402 0.0419 0.666
Number 1 0.16065 0.04801 11.1980 0.0008 1.174
Size 1 -0.04009 0.07026 0.3256 0.5683 0.961

The WLW analysis regards the recurrence times as multivariate failure times, and it models the marginal distribution of each component time with the Cox proportional hazards model. No specific correlation structure is imposed on the multiple failure times. For the kth marginal model, let {\beta}_k denote the row vector of regression parameters, let \hat{{\beta}}_k denote the maximum likelihood estimate of {\beta}_k,let \hat{A}_k denote the covariance matrix obtained by inverting the information matrix, and let Ri denote the matrix of score residuals. Wei, Lin, and Weissfeld (1989) showed that the joint distribution of (\hat{{\beta}}_1, ... ,\hat{{\beta}}_4)'can be approximated by a multivariate normal distribution with mean vector ({\beta}_1, ... ,{\beta}_4)' and robust covariance matrix


with the submatrix Vij given by
V_{ij}= \hat{A}_i(R_i'R_j)\hat{A}_j

The PHREG procedure computes the DFBETA statistics, which are precisely the products R_k\hat{A}_k. By outputting the DFBETA statistics into an OUTPUT data set, you can subsequently use SAS/IML software to compute the robust covariance matrix in a straightforward manner.

In this example, there are four marginal proportional hazards models, one for each recurrence time. Instead of fitting one model at a time, you can fit all four marginal models in one analysis by using the STRATA statement. A new input data set (named Bladder2) is created from the data set Bladder by eliminating observations that have a Visit value of 5. These observations contain the follow-up information beyond the fourth recurrence time and are irrelevant in the fitting of the four marginal models. In addition, four treatment variables, Trt1, Trt2, Trt3, and Trt4, are created from variables Trt and Visit, representing the treatment group for each of the four values of Visit. For instance, Trt1=Trt when the Visit value is 1, and Trt1=0 when the Visit value is not 1. Likewise, variables Number1, Number2, Number3, and Number4 are created from the variables Number and Visit; variables Size1, Size2, Size3, and Size4 are created from the variables Size and Visit.

The following SAS statements create the data set Bladder2:

   data Bladder2;
      set Bladder;
      if Visit < 5;
      Trt1= Trt * (Visit=1);
      Trt2= Trt * (Visit=2);
      Trt3= Trt * (Visit=3);
      Trt4= Trt * (Visit=4);
      Number1= Number * (Visit=1);
      Number2= Number * (Visit=2);
      Number3= Number * (Visit=3);
      Number4= Number * (Visit=4);
      Size1= Size * (Visit=1);
      Size2= Size * (Visit=2);
      Size3= Size * (Visit=3);
      Size4= Size * (Visit=4);
   run;

The following SAS statements fit the marginal models. The parameter estimates are output to a data set named Est1; the DFBETA statistics for the treatment variables are output to a data set named Out1.

   *title 'Fitting Marginal Proportional Hazards Models';
   proc phreg data=Bladder2 outest=Est1;
      model TStop*Status(0)=Trt1-Trt4 Number1-Number4 Size1-Size4;
      output out=Out1 dfbeta=dt1-dt4 / order=data;
      strata Visit;
      id ID;
   run;

The output of this analysis is shown in Output 49.8.2

Output 49.8.2: Fitting Marginal Proportional Hazards Models

The PHREG Procedure

Model Information
Data Set WORK.BLADDER2
Dependent Variable TStop
Censoring Variable Status
Censoring Value(s) 0
Ties Handling BRESLOW

Summary of the Number of Event and Censored Values
Stratum Visit Total Event Censored Percent
Censored
1 1 86 47 39 45.35
2 2 86 29 57 66.28
3 3 86 22 64 74.42
4 4 86 14 72 83.72
Total   344 112 232 67.44

Convergence Status
Convergence criterion (GCONV=1E-8) satisfied.

Model Fit Statistics
Criterion Without
Covariates
With
Covariates
-2 LOG L 880.828 851.435
AIC 880.828 875.435
SBC 880.828 908.057

Testing Global Null Hypothesis: BETA=0
Test Chi-Square DF Pr > ChiSq
Likelihood Ratio 29.3932 12 0.0034
Score 33.0747 12 0.0009
Wald 31.0544 12 0.0019

Analysis of Maximum Likelihood Estimates
Variable DF Parameter
Estimate
Standard
Error
Chi-Square Pr > ChiSq Hazard
Ratio
Trt1 1 -0.51762 0.31576 2.6873 0.1012 0.596
Trt2 1 -0.61944 0.39318 2.4821 0.1151 0.538
Trt3 1 -0.69988 0.45994 2.3155 0.1281 0.497
Trt4 1 -0.65079 0.57744 1.2702 0.2597 0.522
Number1 1 0.23599 0.07608 9.6231 0.0019 1.266
Number2 1 0.13756 0.09190 2.2403 0.1345 1.147
Number3 1 0.16984 0.10521 2.6061 0.1065 1.185
Number4 1 0.32880 0.12528 6.8880 0.0087 1.389
Size1 1 0.06789 0.10125 0.4496 0.5025 1.070
Size2 1 -0.07612 0.13406 0.3224 0.5702 0.927
Size3 1 -0.21131 0.18240 1.3421 0.2467 0.810
Size4 1 -0.20317 0.23018 0.7791 0.3774 0.816


The following SAS statements calculate the robust covariance matrix for the treatment coefficients. The MEANS procedure sums up the DFBETA statistics for each subject and outputs the results to a SAS data set named Out2. The IML procedure then reads the DFBETA statistics from the data set Out2 and computes the robust variance, which is output to a SAS data set called RCov.

   proc means data=Out1 noprint;
      by ID;
      var dt1-dt4;
      output out=Out2 sum=dt1-dt4;

   proc iml;
      use Out2;
      read all var{dt1 dt2 dt3 dt4} into x;
      v=x` * x;
      reset noname;
      vname={"Trt1","Trt2","Trt3","Trt4"};
      print,"Estimated Covariance Matrix",, v[colname=vname
        rowname=vname format=10.5];
      create RCov from v[colname=vname rowname=vname];
      append from v[rowname=vname];

The estimated robust covariance matrix is displayed in Output 49.8.3.

Output 49.8.3: Robust Covariance Matrix for the Treatment Coefficients

Estimated Covariance Matrix

  Trt1 Trt2 Trt3 Trt4
Trt1 0.09456 0.06018 0.05677 0.04378
Trt2 0.06018 0.13243 0.13012 0.11604
Trt3 0.05677 0.13012 0.17236 0.15909
Trt4 0.04378 0.11604 0.15909 0.23981


The approximate multivariate normal distribution of the parameter estimators provides a basis for simultaneous inferences about the parameters. For example, you can test jointly the null hypothesis of no treatment effect for each tumor recurrence, or you can estimate the coefficient for the common treatment effect. A detailed IML program for the analysis is given by the following SAS statements (results not shown):

   proc iml;
      use Est1;
      read all var{Trt1 Trt2 Trt3 Trt4} into Trt;
      b= Trt`;
      use Out2;
      read all var{dt1 dt2 dt3 dt4} into x;
      v=x` * x;
      nparm= nrow(b);
      se=sqrt(vecdiag(v));
      reset noname;
      stitle={"Estimate", "   Std Error"};
      vname={"Trt1","Trt2","Trt3","Trt4"};
      tmpprt= b || se;
      print,tmpprt[colname=stitle rowname=vname format=10.5];
      print,"Estimated Covariance Matrix",,
             v[colname=vname rowname=vname format=10.5];

      /* H0: beta11=beta12=beta13=beta14=0 */
      chisq= b` * inv(v) * b;
      df= nrow(b);
      p= 1-probchi(chisq,df);
      print ,,"Testing H0: no treatment effects", ,
            "Wald Chi-Square = " chisq, "DF = " df,
            "p-value = "p[format=5.4],;

      /* Assume beta11=beta12=beta13=beta14 and
         estimate the common value */
      c= {1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1};
      cb= c * b;
      si= c * v * t(c);
      e= j(4,1,1);
      isi=inv(si);
      h= inv(e` * isi * e) * isi * e;
      b1= t(h) * cb;
      se= sqrt(t(h) * si * h);
      zscore= b1 / se;
      p= 1- probchi( zscore # zscore, 1);
      print ,"Estimation of the Common Parameter for Treatment",,
            "Optimal Weights = "h,
            "Estimate = " b1,
            "Standard Error = " se,
            "z-score =" zscore,
            "2-sided p-value = " p[format=5.4];
   quit;

The PHREG procedure can also be used to fit the PWP model. In the PWP model, the risk set for the (k+1) recurrence is restricted to those patients who have experienced the first k recurrences. For example, a patient who experienced only one recurrence is an event observation for the first recurrence; this patient is a censored observation for the second recurrence and should not be included in the risk set for the third or fourth recurrence. The following DATA step eliminates those observations that should not be in the risk sets, forming a new input data set (named Bladder3) for fitting the PWP models. The gap times between successive recurrences are also calculated.

The following SAS statements create the data set Bladder3:

   data Bladder3(drop=lstatus);
      retain lstatus;
      set Bladder2;
      by ID;
      if first.ID then lstatus=1;
      if (Status=0 and lstatus=0) then delete;
      lstatus=Status;
      GapTime=TStop-TStart;
   run;

The following statements fit the PWP total time model with noncommon effects:

   *title2 'PWP Total Time Model with Noncommon Effects';
   proc phreg data=Bladder3;
      model TStop*Status(0)=Trt1-Trt4 Number1-Number4
                                     Size1-Size4;
      strata Visit;
   run;

The following statements fit the PWP gap time model with noncommon effects:

   *title2 'PWP Gap Time Model with Noncommon Effects';
   proc phreg data=Bladder3;
      model GapTime*Status(0)=Trt1-Trt4 Number1-Number4
                              Size1-Size4;
      strata Visit;
   run;

Results of these two analyses are shown in Output 49.8.4 and Output 49.8.5, respectively.

Output 49.8.4: Fitting PWP Total Time Model with Noncommon Effects

The PHREG Procedure

Model Information
Data Set WORK.BLADDER3
Dependent Variable TStop
Censoring Variable Status
Censoring Value(s) 0
Ties Handling BRESLOW

Summary of the Number of Event and Censored Values
Stratum Visit Total Event Censored Percent
Censored
1 1 86 47 39 45.35
2 2 47 29 18 38.30
3 3 29 22 7 24.14
4 4 22 14 8 36.36
Total   184 112 72 39.13

Convergence Status
Convergence criterion (GCONV=1E-8) satisfied.


The PHREG Procedure

Model Fit Statistics
Criterion Without
Covariates
With
Covariates
-2 LOG L 743.098 725.677
AIC 743.098 749.677
SBC 743.098 782.299

Testing Global Null Hypothesis: BETA=0
Test Chi-Square DF Pr > ChiSq
Likelihood Ratio 17.4211 12 0.1344
Score 18.5546 12 0.0999
Wald 17.7388 12 0.1239

Analysis of Maximum Likelihood Estimates
Variable DF Parameter
Estimate
Standard
Error
Chi-Square Pr > ChiSq Hazard
Ratio
Trt1 1 -0.51757 0.31576 2.6868 0.1012 0.596
Trt2 1 -0.42584 0.40227 1.1206 0.2898 0.653
Trt3 1 -0.89894 0.53956 2.7757 0.0957 0.407
Trt4 1 -0.23739 0.68274 0.1209 0.7281 0.789
Number1 1 0.23605 0.07607 9.6287 0.0019 1.266
Number2 1 0.00117 0.09372 0.0002 0.9900 1.001
Number3 1 0.01468 0.13253 0.0123 0.9118 1.015
Number4 1 0.29306 0.22067 1.7637 0.1842 1.341
Size1 1 0.06790 0.10125 0.4498 0.5024 1.070
Size2 1 -0.12515 0.11708 1.1425 0.2851 0.882
Size3 1 -0.21520 0.17801 1.4615 0.2267 0.806
Size4 1 0.25135 0.29077 0.7472 0.3874 1.286

Output 49.8.5: Fitting PWP Gap Time Model with Noncommon Effects

The PHREG Procedure

Model Information
Data Set WORK.BLADDER3
Dependent Variable GapTime
Censoring Variable Status
Censoring Value(s) 0
Ties Handling BRESLOW

Summary of the Number of Event and Censored Values
Stratum Visit Total Event Censored Percent
Censored
1 1 86 47 39 45.35
2 2 47 29 18 38.30
3 3 29 22 7 24.14
4 4 22 14 8 36.36
Total   184 112 72 39.13

Convergence Status
Convergence criterion (GCONV=1E-8) satisfied.

Model Fit Statistics
Criterion Without
Covariates
With
Covariates
-2 LOG L 735.076 717.268
AIC 735.076 741.268
SBC 735.076 773.890


The PHREG Procedure

Testing Global Null Hypothesis: BETA=0
Test Chi-Square DF Pr > ChiSq
Likelihood Ratio 17.8089 12 0.1216
Score 19.6097 12 0.0748
Wald 18.4759 12 0.1020

Analysis of Maximum Likelihood Estimates
Variable DF Parameter
Estimate
Standard
Error
Chi-Square Pr > ChiSq Hazard
Ratio
Trt1 1 -0.51757 0.31576 2.6868 0.1012 0.596
Trt2 1 -0.25911 0.40511 0.4091 0.5224 0.772
Trt3 1 0.22105 0.54909 0.1621 0.6873 1.247
Trt4 1 -0.19498 0.64184 0.0923 0.7613 0.823
Number1 1 0.23605 0.07607 9.6287 0.0019 1.266
Number2 1 -0.00571 0.09667 0.0035 0.9529 0.994
Number3 1 0.12935 0.15970 0.6561 0.4180 1.138
Number4 1 0.42079 0.19816 4.5091 0.0337 1.523
Size1 1 0.06790 0.10125 0.4498 0.5024 1.070
Size2 1 -0.11636 0.11924 0.9524 0.3291 0.890
Size3 1 0.24995 0.23113 1.1695 0.2795 1.284
Size4 1 0.03557 0.29043 0.0150 0.9025 1.036


The following statements fit the PWP total time model with common effects:

   *title2 'PWP Total Time Model with Common Effects';
   proc phreg data=Bladder3;
      model TStop*Status(0)=Trt Number Size;
      strata Visit;
   run;

The following statements fit the PWP gap time model with common effects:

   *title2 'PWP Gap Time Model with Common Effects';
   proc phreg data=Bladder3;
      model GapTime*Status(0)=Trt Number Size;;
      strata Visit;
   run;

Results of these two analyses are not shown in this document.

Chapter Contents
Chapter Contents
Previous
Previous
Next
Next
Top
Top

Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.