STAT 572 Survival Analysis HW7
Problem: Data set: WHAS, Covariates: AGE, SEX, MITYPE. Perform survival regression.

1. Data Description.
In this case, there is one continuous variable (AGE), and two categorical variables (SEX, MITYPE). We are studying how these three predicators affect patients’ survival time. For SEX, Male=0 and Female=1. For MI Type of patients, there are three values: 1 | Q-Wave | ST segment elevation myocardial infarction (STEMI) | 2 | Not Q-Wave | non-ST segment elevation myocardial infarction (NSTEMI) | 3 | Indeterminate | |
Table 1.1 MI type explanation

2. Variable Analysis
Since AGE is a continuous variable, to be convenience, we make it as categorical variable CAGE using 25%, 50% and 75% quantile at FSTAT=1. Using PROC UNIVARIATE to decide the cut point as follows: Variables | 25% Quantile | 50% Quantile | 75% Quantile | AGE | 64 | 72 | 80 |
Table 2.1 Quantiles Cut Point for AGE
Thus, all the variables are categorical, such that we can use the LIFETEST to test if different levels of covariates have the same effect on the survival curve.
We are doing hypothesis testing for each predictor:
H0:SSEX=0=SSEX=1 vs. H1:Not H0
H0:SMI=1=SMI=2=S(MI=2) vs. H1:Not H0
H0:SCAGE=0=SCAGE=1=SCAGE=2=S(CAGE=3) vs. H1:Not H0
First, we will test the equality over SEX, and compare the KM curve using LIFETEST and PHREG. Test | Chi-Square | DF | Pr > Chi-Square | Log-Rank | 9.4271 | 1 | 0.0021 | Wilcoxon | 9.2239 | 1 | 0.0024 | -2Log(LR) | 12.0677 | 1 | 0.0005 |
Table 2.2 Test of Equality over SEX using LIFETEST

Figure 2.1 KM Survival Curve for Male (Red) and Female (Blue)
From Table 2.2 and Figure 2.1, we can easily see that there is significant difference between male and female.
Then, we use PHREG to test SEX and draw the PHREG curve. Test | Chi-Square | DF | Pr>ChiSq | Likelihood Ratio | 9.1280 | 1 | 0.0025 | Score | 9.3800 | 1 | 0.0022 | Wald | 9.2647 | 1 | 0.0023 |
Table 2.3 Significance test for SEX effect using PHREG

Figure 2.2 Comparison of KM Curve (Blue) and COX Curve (Black) for SEX
For the MI Type, since there are only 6 observations with 1 censored when MI=3, we will test twice with MI=3 and without MI=3. Test | Chi-Square | DF | Pr > Chi-Square | | Test | Chi-Square | DF | Pr > Chi-Square | Log-Rank | 13.1026 | 2 | 0.0014 | | Log-Rank | 2.2637 | 1 | 0.1324 | Wilcoxon | 16.7644 | 2 | 0.0002 | | Wilcoxon | 6.5772 | 1 | 0.0103 | -2Log(LR) | 18.0153 | 2 | 0.0001 | | -2Log(LR) | 0.8075 | 1 | 0.3688 |
Table 2.3 Test for Equality over MI Type (Left: Include MI=3)

Figure 2.2 KM Survival Curve for MI Type (Red: Q-Wave; Blue: Not Q-Wave)
From table 2.3 and Figure 2.2, there is not a significant difference between Q-Wave and Not Q-Wave when the survival time is getting longer. However, MI=3 has a significant difference between the other two type, and people with MI=3 always have very short survival time. Thus, we consider to reclassify the MI Type as MI=0 if Q-Wave or Not-Q-Wave, and MI=1 if not determined.
In the Cox Model, the test results are different. Testing Global Null Hypothesis: BETA=0 | Test | Chi-Square | DF | Pr > ChiSq | Likelihood Ratio | 0.5332 | 1 | 0.4652 | Score | 0.5299 | 1 | 0.4666 | Wald | 0.5296 | 1 | 0.4668 |
Table 2.4 Significance test for MITYPE effect using PHREG

Figure 2.3 Comparison of KM Curve (Blue) and COX Curve (Black) for MITYPE
Since there are 6 observations when MI=3, the KM curve for MI=3 is only 6 points at left side of the plot, while PHREG curves are parallel for each MITYPE. However, from table 2.4, the MITYPE is not significant.
Then, we will test the equality over categorical age. Test | Chi-Square | DF | Pr > Chi-Square | Log-Rank | 81.2589 | 3 | <.0001 | Wilcoxon | 74.1006 | 3 | <.0001 | -2Log(LR) | 122.7237 | 3 | <.0001 |
Table 2.4 Test of Equality over CAGE

Figure 2.4 KM Survival Curve for CAGE (0: Red; 1: Blue; 2: Green; 3: Black)
From table 2.4 and figure 2.3, we can easily see that there exist significant differences between each age categories. Also, from the plot, we can see with people getting older, their survival probability is reducing.

Figure 2.5 of KM Curve (Blue) and COX Curve (Black) for CAGE
From figure 2.5, Cox curve and KM curve are not very close.

3. Comparing Cox to Weibull Estimated Curve
In the last homework, we have the Weibull model which contains variable AGE SEX and their interaction. If we treat AGE and a categorical variable CAGE, then we can get four estimated Weibull curves and Cox curves, as well as we can get two curves for SEX.

Figure 3.1 Weibull (Red) and COX (Black) curve for SEX

Figure 3.2 Weibull (Red) and COX (Black) curve for CAGE
Since the estimated Weibull curve will extend to a very large time value, such that we cannot see the Cox curve if they are in the same figure. We just cut off a certain range of the curve. From figure 3.1 and 3.2, we can see the Weibull curve is close to Cox curve.

4. Checking Assumptions.
Since the data was collected individually, we assume that the observations are independent from each other. We need to check the assumption of linearity and proportional hazard. a) Check Linearity:
Only variable AGE is continuous, so we need to check the linearity assumption of AGE. h(t)=h0(t)expXβ => loght=logh0t+Xβ
We will use several methods to check the linearity of AGE, because it is continuous variable. 1) Categorical Method
Making AGE to five categories, then the output will give the estimated coefficients for Intercept, CAGE0, CAGE1, CAGE2 and CAGE3.

CAGE_0 | -1.51143 | CAGE_1 | -0.97181 | CAGE_2 | -0.47724 | CAGE_3 | -0.04285 |
Table 4.1 Parameter Estimation for CAGE

Then, we have the plot for categorical age:

Figure 4.1 CAGE vs. Log(h/h0)

From Figure 4.1, we can see it is a straight line. Thus, the linearity of age is appropriate.

Also, using TEST statement in PHREG to test Linearity, the result is: Linear Hypotheses Testing Results

Wald Label Chi-Square DF Pr > ChiSq

Test 1 20.5238 1 <.0001

Table 4.2 Linearity Test using Test statement

2) Restricted Cubic Spline
The RCS Macro will generate a new variable which is approximately a cubic term. In this problem, we use RCS Macro to generate variable AGE1. Then test significance.
H0: βAGE1=0 vs. H1: βAGE1≠0
SAS Output: Analysis of Maximum Likelihood Estimates | Parameter | DF | Parameter
Estimate | Standard
Error | Chi-Square | Pr > ChiSq | Hazard
Ratio | SEX | 1 | 0.08397 | 0.13243 | 0.4021 | 0.5260 | 1.088 | AGE | 1 | 0.05187 | 0.01097 | 22.3743 | <.0001 | 1.053 | AGE1 | 1 | -0.01297 | 0.01337 | 0.9403 | 0.3322 | 0.987 |
Table 4.3 Cubic Restricted Spline Test for Linearity

From the table above, the new variable AGE1 is obviously not significant. Thus, we accept null hypothesis, which is AGE is linearly entering the model.

b) Check Proportional Hazard
We have several methods to test the proportional hazard, such as LLS curve, Schoenfeld Residual plot and checking if there is any time dependent variable. 1) LLS Curve

Figure 4.2 LLS Curve for Categorical Covariates
Left side is LLS curve generated from LIFETEST, and right side is generated from PHREG. Since in Cox model, the proportional hazard is assumed, the LLS curve is always parallel no matter what the data is. Thus, the three plots in right side of Figure 4.2 cannot be used to check the PH assumption.
By looking at the left side of figure 4.2, the LLS curve for SEX and CAGE are roughly parallel, while the LLS curve for MI type is not.
For MI Type, since there are only six observations out of 481 total observations when MI=3, thus we will ignore these 6 observations to test if this covariate is significant, and based on the result of LIFETEST, MI=1 and MI=2 do not have significant effect for the survival probability. Thus, we will not use this covariate in our model, and still use the 481 total observations.

2) Shoenfeld Residual Plot
Since the Schoenfeld residuals are, in principle, independent of time, a plot that shows a non-random pattern against time is evidence of violation of the PH assumption.
We plot the Shoenfeld Residual of AGE and SEX as follows:

Figure 4.3 Schoenfeld Residual Plot for AGE and SEX
From the above plots, we can the Schoenfeld residuals are symmetrically distributed around 0, and they are independent from time. Thus the PH assumption is satisfied for the two covariates.

3) Time Dependent Method
We know that if any covariates are time dependent, then the Proportional Hazard assumption for the covariate is violated. Using the interaction between a covariate and time as a time dependent variable, and adding it in PROC PHREG procedure, we can get the p-value of significant test for the time dependent variable. Thus, we can see if the PH assumption is violated.

Analysis of Maximum Likelihood Estimates | Parameter | DF | Parameter
Estimate | Standard
Error | Chi-Square | Pr > ChiSq | Hazard
Ratio | AGE | 1 | 0.04687 | 0.00672 | 48.7004 | <.0001 | 1.048 | SEX | 1 | 0.10863 | 0.15995 | 0.4613 | 0.4970 | 1.115 | TSEX | 1 | -0.0000218 | 0.0001262 | 0.0299 | 0.8626 | 1.000 | TAGE | 1 | -5.4977E-6 | 5.31812E-6 | 1.0687 | 0.3012 | 1.000 |
Table 4.4 Checking PH using Time dependent variable
From the output above, we can see variable TSEX and TAGE, which are time dependent variables, are not significant. No evidence was shown to reject null hypothesis, which is they are time independent. It indicates that PH assumption is satisfied.
For the variable MITYPE, which we have decided not to use, when we check its PH assumption, it is not satisfied. The results are not shown, because it will not be used for our model. If some other variables, which we may use in future, does not satisfy the PH assumption, we will use BY or STRAT statement to fit separate Cox model. 4) Some other method in SAS
Using ASSESS statement, we can test the PH assumption. The output is as follows: Supremum Test for Proportionals Hazards
Assumption | Variable | Maximum Absolute
Value | Replications | Seed | Pr >
MaxAbsVal | AGE | 0.6458 | 1000 | 159346001 | 0.7500 | SEX | 0.8925 | 1000 | 159346001 | 0.3670 |
Table 4.5 Checking PH using ASSESS Statement in SAS
The p-values are large enough to accept null hypothesis, which is PH assumptions are satisfied. In this statement, SAS also gives some plots. Figure 4.4 Testing PH using ASSESS
The ASSESS statement produces a graphical display of empirical score process, which is based on the martingale residuals. [See Paul D. Allison, Survival Analysis Using SAS, page 173-174].

5. Checking Interactions.
When we add the interaction of AGE and SEX, the result shows all of the three terms are significant. Analysis of Maximum Likelihood Estimates | Parameter | DF | Parameter
Estimate | Standard
Error | Chi-Square | Pr > ChiSq | Hazard
Ratio | Label | AGE | 1 | 0.05175 | 0.00730 | 50.2241 | <.0001 | . | | SEX | 1 | 1.55289 | 0.79231 | 3.8414 | 0.0500 | . | | AGE*SEX | 1 | -0.02040 | 0.01088 | 3.5147 | 0.0608 | . | AGE * SEX |
Table 5.1 Significant test for interaction
As one can see, however, the Hazard Ratio is left blank if we add any interaction to the model. It is because that there will not be an exact hazard ratio to interpret. We use Likelihood Ratio Test to test if it is necessary to add this interaction. -2 Log L | Reduced Model (without interaction) | Full Model (with interaction) | 2770.751 | 2767.265 |
Table 5.2 -2 Log L for full model and reduced model
Obviously, the test statistic, the difference between the two values (3.5) is less than Chi (0.95, 1), which is 3.84, but it is bigger than Chi (0.90, 1) = 2.705. Given α=0.1, we need to add this interaction.

6. Interpretation for parameters.
In Cox model, the parameters for covariates are easy to interpret. Their exponential terms are hazard ratio. However, to interpret them, the model should not include any interactions. Analysis of Maximum Likelihood Estimates | Parameter | DF | Parameter
Estimate | Standard
Error | Chi-Square | Pr > ChiSq | Hazard
Ratio | AGE | 1 | 0.04281 | 0.00553 | 59.8881 | <.0001 | 1.044 | SEX | 1 | 0.08404 | 0.13269 | 0.4011 | 0.5265 | 1.088 |
Table 6.1 Model without interaction
The hazard ratio in above table is the exponential function of the parameters.
A hazard ratio of two, for example, means a patient in one treatment group who has not died (or progressed, or whatever end point is tracked) at a certain time point has twice the probability of having died (or progressed...) by the next time point compared to a patient in the other treatment group.
In our model, the hazard ratio of SEX is 1.088, means female (SEX=1) who has not died at a certain time point has 1.088 times the probability of having died by next time point compared to male.
Using the same way, we can interpret the hazard ratio of AGE. The hazard ratio is between patient who is one more year older and patient at a certain age. 7. Prediction
In the end, we will make a prediction for the one who has the mean value of age (67), the worst value of sex (Female=1). Since our best model does not contain MI Type, we will skip this covariate. We will use the model in Table 6.1 (without interaction) and Table 5.1 (with interaction) to predict the patient’s survival probability respectively.
From SAS data step, by looking at the data set generated from SAS, we get the estimated probability. With interaction | Without interaction | 0.682—0.684 | 0.702—0.704 |
Table 7.1 Prediction of Survival Probability using two Cox models
There is not a significant difference between these two models for the predicted survival probability. Also, comparing to 0.678, the estimated survival probability using Weibull Model in HW5, we find they are very close.
Also, we can get the Cox curve comparing to Weibull Curve we generated in HW5.

Figure 7.1 Comparing Estimated Weibull (Red) and Cox (Black Step) Curve for Sixty-Seven Years Old Female Patient
Figure 7.1 shows both Weibull curve and Cox curve under the same predictor variables. They are very close. Further study is needed to assessing which model is more appropriate.

APPENDIX: MAIN SAS CODE
OPTIONS LS=80 PS=400 FORMDLIM='=';
DATA WHAS;
INFILE 'I:\Documents\My Wealth\572_Survival_Analysis\HW5\WHAS.txt'; INPUT ID AGE SEX CPK SHO CHF MIORD MITYPE YEAR YRGRP LENSTAY DSTAT LENFOL FSTAT;
KEEP ID AGE SEX MITYPE LENFOL FSTAT;
RUN;

PROC UNIVARIATE;
VAR AGE; WHERE FSTAT=1;
RUN;

DATA WHAS; SET WHAS; CAGE1=(64<=AGE<72); CAGE2=(72<=AGE<80); CAGE3=(AGE>=80); CAGE=1*CAGE1+2*CAGE2+3*CAGE3; RUN;

ods graphics on;
*Plot KM Curve and LLS curve using lifetest;
DATA AA; INPUT SEX @@; CARDS;
0 1
;run;
PROC LIFETEST DATA=WHAS OUTS=KM1 NOTABLE; STRATA SEX; *Change to SEX and CAGE to test separately; TIME LENFOL*FSTAT(0); RUN; DATA KM1; SET KM1; IF SEX=0 THEN SURV0=SURVIVAL; IF SEX=1 THEN SURV1=SURVIVAL; run;
PROC PHREG DATA=WHAS; MODEL LENFOL*FSTAT(0)=MITYPE;run; BASELINE OUT=BB1 COVARIATES=AA SURVIVAL=SURVIVAL;run;

DATA CC1; SET BB1; IF SEX=0 OR SEX=1; IF SEX=0 THEN SUR0=SURVIVAL; IF SEX=1 THEN SUR1=SURVIVAL;run;

GOPTIONS RESET=ALL;
DATA FINAL1; SET KM1 CC1;
PROC GPLOT; PLOT(SURV0 SURV1 )*LENFOL/OVERLAY; SYMBOL1 V=NONE I=STEPL C=BLACK L=1; SYMBOL2 V=NONE I=STEPL C=BLACK L=1;
*this will plot the 2 KM survival curves ;
* notice the curves have different steps locations;
RUN;
PROC GPLOT; PLOT( SUR0 SUR1)*LENFOL/OVERLAY; SYMBOL1 V=NONE I=STEPL C=BLACK L=1; SYMBOL2 V=NONE I=STEPL C=BLACK L=1;
*this will plot the 2 PHREG curves;
* notice the curves have the same steps locations;
RUN;
PROC GPLOT; PLOT(SURV0 SURV1 SUR0 SUR1)*LENFOL/OVERLAY; SYMBOL1 V=NONE I=STEPL C=BLUE L=2 W=2; SYMBOL2 V=NONE I=STEPL C=BLUE L=2 W=2; SYMBOL3 V=NONE I=STEPL C=BLACK L=1 W=2; SYMBOL4 V=NONE I=STEPL C=BLACK L=1 W=2;
*this will plot the 2 KM survival curves and the 2 PHREG curves;
RUN;

DATA EST; INPUT SEX; CONTROL=1; LENFOL=0; FSTAT=0;
CARDS;
0
1
;
DATA LAST; SET WHAS EST;
PROC LIFEREG DATA=LAST; MODEL LENFOL*FSTAT(0)=SEX /D=WEIBULL; OUTPUT OUT=WEIB QUANTILES=.01 .05 .10 .15 .20 .25 .30 .35 .40 .45 .50 .55 .60 .65 .70 .75 .80 .85 .90 .95 .99 P=PRED CONTROL=CONTROL ;run;
DATA WEIB; SET WEIB; SURV=1-_PROB_; LENFOL=PRED; IF SEX=0 THEN SURVV0=SURV; IF SEX=1 THEN SURVV1=SURV;
DATA FINAL11; SET CC1 WEIB;
RUN;
PROC GPLOT; PLOT(SUR0 SUR1 SURVV0 SURVV1)*LENFOL /OVERLAY HAXIS=1 TO 8000 BY 400;; SYMBOL1 V=NONE I=STEPL C=BLACK L=1 W=2; SYMBOL2 V=NONE I=STEPL C=BLACK L=1 W=2; SYMBOL3 V=NONE I=JOIN C=RED L=3 W=2; SYMBOL4 V=NONE I=JOIN C=RED L=3 W=2; TITLE 'WEIBULL and COX CURVES for SEX';
RUN;
*Check Linearity;
PROC UNIVARIATE DATA=WHAS;
VAR AGE; WHERE FSTAT=1;
OUTPUT OUT=AAA PCTLPTS=0 20 40 60 80 100 PCTLPRE=AGE_;
RUN;
DATA WHAS1; SET WHAS; CAGE1=(62<=AGE<70); CAGE2=(70<=AGE<76); CAGE3=(76<=AGE<82); CAGE4=(AGE>=82); CAGE=1*CAGE1+2*CAGE2+3*CAGE3+4*CAGE4; AGE2=AGE*AGE; RUN;
PROC PHREG DATA=WHAS1; CLASS CAGE; MODEL LENFOL*FSTAT(0)=CAGE; RUN;

*Schoenfeld Residual;
PROC PHREG DATA=WHAS ; MODEL LENFOL*FSTAT(0)=AGE SEX/RISKLIMITS; OUTPUT OUT=C RESDEV=DEV WTRESSCH=SCH_AGE SCH_SEX RESSCH=AGE1 SEX1 ;RUN;

PROC GPLOT DATA=C; PLOT AGE1*LENFOL; PLOT SEX1*LENFOL; RUN;

*options WTRESSCH,DFBETA, and LD are also useful in assessing model fit;
DATA D; SET C; LLENFOL=LOG(LENFOL); SLENFOL=SQRT(LENFOL);
PROC CORR DATA=D; VAR LENFOL LLENFOL SLENFOL; WITH SCH_AGE SCH_SEX; RUN;
PROC GPLOT DATA=C; *plot the deviance residuals; PLOT DEV*AGE; SYMBOL1 V=DOT C=BLACK; RUN;
PROC GPLOT DATA=C; PLOT SCH_AGE*LENFOL / CFRAME=white OVERLAY VAXIS=axis1 HAXIS=axis2 FRAME VREF=0 VMINOR=0 HMINOR=0 CAXIS = blue NAME='plot3'; symbol value=dot i=sm60s h=1.2 w=3; TITLE 'Schoenfeld Residuals for AGE'; axis1 label =(a=90 r=0 f='Arial''Schoenfeld Residuals') value=(f='Arial' ); axis2 label=( f='Arial' 'Time')value=(f='Arial') ; RUN; QUIT; *plot the Schoenfeld residuals for AGE;

SYMBOL1 V=DOT C=BLACK; RUN;
PROC GPLOT DATA=C; PLOT SCH_SEX*LENFOL / CFRAME=white OVERLAY VAXIS=axis1 HAXIS=axis2 FRAME VREF=0 VMINOR=0 HMINOR=0 CAXIS = blue NAME='plot3'; symbol value=dot i=sm60s h=1.2 w=3; TITLE 'Schoenfeld Residuals for SEX'; axis1 label =(a=90 r=0 f='Arial''Schoenfeld Residuals') value=(f='Arial' ); axis2 label=( f='Arial' 'Time')value=(f='Arial') ; RUN; QUIT; *plot the Schoenfeld residuals for SEX (not very informative); TITLE 'Schoenfeld Residuals for SEX'; SYMBOL1 V=DOT C=BLACK; RUN;

PROC PRINT DATA=C; RUN;

*ASSESS Statement;
PROC PHREG DATA=WHAS PLOTS(CL)=(SURV SUR SURVIVAL S); MODEL LENFOL*FSTAT(0)=AGE SEX; ASSESS PH / RESAMPLE;
RUN;

*Time dependent Method;
PROC PHREG DATA=WHAS; MODEL LENFOL*FSTAT(0)=AGE SEX TSEX TAGE; TSEX=SEX*LENFOL; TAGE=AGE*LENFOL; PROPORTIONALITY_TEST: TEST TSEX, TAGE;
*this is a time-dependent variable that tests for PH;
RUN;
*Estimated Survival Curve for mean age people;
DATA AAA; INPUT AGE SEX; CARDS;
67 1
;
PROC PHREG DATA=WHAS; MODEL LENFOL*FSTAT(0)=AGE SEX; BASELINE OUT=BBB COVARIATES=AAA SURVIVAL=SURVIVAL LOGLOGS=LOGLOGS;
DATA CCC; SET BBB; IF SEX=1 AND AGE=67 THEN SURV1=SURVIVAL; LLENFOL=LOG(LENFOL); proc sort; by LENFOL;
PROC GPLOT; PLOT SURV1 *LENFOL/OVERLAY; TITLE 'Cox Curve for 67 years old Female Patients'; SYMBOL1 V=dot I=STEPL C=RED ; run; DATA EST1; INPUT CONTROL SEX AGE LENFOL FSTAT;
CARDS;
1 1 67 0 0
;
DATA LAST; SET WHAS EST1;run;
PROC LIFEREG DATA=LAST; MODEL LENFOL*FSTAT(0)=AGE SEX/D=WEIBULL; OUTPUT OUT=WEIB1 QUANTILES=.01 .05 .10 .15 .20 .25 .30 .35 .40 .45 .50 .55 .60 .65 .70 .75 .80 .85 .90 .95 .99 P=PRED CONTROL=CONTROL ; RUN;
DATA WEIB2; SET WEIB1; SURV2=1-_PROB_; LENFOL=PRED;
PROC SORT DATA=WEIB2; BY LENFOL;
DATA FINAL; SET CCC WEIB2;
PROC GPLOT; PLOT (SURV2 SURV1)*LENFOL/OVERLAY HAXIS=1 TO 5000 BY 400 VAXIS=0.3 TO 1 BY 0.05; SYMBOL1 V=NONE I=JOIN C=RED; SYMBOL2 V=NONE I=STEPL C=BLACK; TITLE 'Weibull and Cox Curve for 67 years old Female Patients'; *this plots the Kaplan-Meier curve together with the estimated Weibull survival curve;
RUN;