Plot of Hazard Ratio against a Continuous Predictor

[SriniR]

I want to plot hazard ratio against time to treatment. As a first step, I wanted to replicate the procedure that plots survival probability against time to treatment, in the post https://communities.sas.com/t5/SAS-Procedures/Output-covariates-from-proc-phreg-not-working/m-p/2728..., but I am getting the following message and am not able to plot the survival probability against time to treatment.

 73        

 74         proc phreg data = all_os;

 75         class  gender_m (ref='1') /* bmi_cat (ref='2') rgn (ref='1') stg_grp  (ref='1') cci_grp (ref = '1') */;

 76         model os*os_censor(0)=time_to_trtmt gender_m /*age_index  rgn bmi_cat  stg_grp cci_grp */ / RISKLIMITS=PL ;

 77         *   model os*os_censor(0) = time_to_trtmt;

 78            store out = out_tti_survival;

 79         run;

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

 NOTE: The PHREG procedure generated the model item store WORK.OUT_TTI_SURVIVAL.

 NOTE: PROCEDURE PHREG used (Total process time):

       real time           0.17 seconds

       cpu time            0.17 seconds

 80        

 81         ods graphics on;

 82         proc plm restore = out_tti_survival;

 83            effectplot slicefit (x = time_to_trtmt sliceby = gender_m) ;

 84         ods graphics off;

 85         run;

 WARNING: The EFFECTPLOT statement cannot be used for the specified model.

 NOTE: PROCEDURE PLM used (Total process time):

       real time           0.02 seconds

       cpu time            0.03 seconds

I am using SAS Studio Release: 3.8 (Basic Edition) and SAS release: 9.04.01M6P11072018.

As a next step, I want to plot the hazard ratio against time to treatment. 

Finally, I want to be able to use Recursive Partitioning Analysis (RPA) to get a plot similar to the attached one. This plot is from Kompelli AR, LI H, Neskey DM. Impact of Delay in Treatment Initiation on Overall Survival in Laryngeal Cancers. Otolaryngology–Head and Neck Surgery 2019, Vol. 160(4) 651–657.

[Reeza]

Isn't the assumption for cox regression that the hazardratio is constant over time? Do you have a time dependent variable in your model?

[SriniR]

---

@Reeza wrote:
Isn't the assumption for cox regression that the hazardratio is constant over time? Do you have a time dependent variable in your model?

---

Please see the article http://support.sas.com/resources/papers/proceedings09/236-2009.pdf which is given in the post for which I have provided the link. As Rick has suggested in that post, I tried to use PROC PLM which does not seem to work.

[Reeza]

Survival is not the same as the hazardratio, so are you asking for the predicted survival or the hazardratio? Your question seems to consider them the same thing which isn't correct AFAIK.

Assuming you're looking at

[SriniR]

---

@Reeza wrote:
Survival is not the same as the hazardratio, so are you asking for the predicted survival or the hazardratio? Your question seems to consider them the same thing which isn't correct AFAIK.

Assuming you're looking at

---

No, I am not assuming that survival and HR are the same thing. As a first step, I wanted to plot the survival against time to treatment. I want to see what the survival probability versus time to treatment plot looks like. Then I want to figure out how to plot HR against time to treatment. Finally, I have to figure out how to use Recursive Partitioning Analysis to get the HR versus time to treatment plot and see what difference RPA makes to this plot. Hope I have cleared things up a little bit.

[Reeza]

I cannot seem to find a logical reason why it shouldn't work and it doesn't work for any variation so perhaps that option isn't supported for PROC PHREG though it doesn't say that anywhere as far as I can see. 

Maybe @Rick_SAS himself can comment 🙂

There's full code below to generate an example/test code.

Spoiler

proc format;
value yesno 0='no' 10='yes';
run;

data VALung;
drop check m;
retain Therapy Cell;
infile cards column=column;
length Check $ 1;
label Time='time to death in days'
Kps='Karnofsky performance scale'
Duration='months from diagnosis to randomization'
Age='age in years'
Prior='prior therapy'
Cell='cell type'
Therapy='type of treatment';
format Prior yesno.;
M=Column;
input Check $ @@;
if M>Column then M=1;
if Check='s'|Check='t' then do;
input @M Therapy $ Cell $;
delete;
end;
else do;
input @M Time Kps Duration Age Prior @@;
Status=(Time>0);
Time=abs(Time);
end;
datalines;
standard squamous
72 60 7 69 0 411 70 5 64 10 228 60 3 38 0 126 60 9 63 10
118 70 11 65 10 10 20 5 49 0 82 40 10 69 10 110 80 29 68 0
314 50 18 43 0 -100 70 6 70 0 42 60 4 81 0 8 40 58 63 10
144 30 4 63 0 -25 80 9 52 10 11 70 11 48 10
standard small
30 60 3 61 0 384 60 9 42 0 4 40 2 35 0 54 80 4 63 10
13 60 4 56 0 -123 40 3 55 0 -97 60 5 67 0 153 60 14 63 10
59 30 2 65 0 117 80 3 46 0 16 30 4 53 10 151 50 12 69 0
22 60 4 68 0 56 80 12 43 10 21 40 2 55 10 18 20 15 42 0
139 80 2 64 0 20 30 5 65 0 31 75 3 65 0 52 70 2 55 0
287 60 25 66 10 18 30 4 60 0 51 60 1 67 0 122 80 28 53 0
27 60 8 62 0 54 70 1 67 0 7 50 7 72 0 63 50 11 48 0
392 40 4 68 0 10 40 23 67 10
standard adeno
8 20 19 61 10 92 70 10 60 0 35 40 6 62 0 117 80 2 38 0
132 80 5 50 0 12 50 4 63 10 162 80 5 64 0 3 30 3 43 0
95 80 4 34 0
standard large
177 50 16 66 10 162 80 5 62 0 216 50 15 52 0 553 70 2 47 0
278 60 12 63 0 12 40 12 68 10 260 80 5 45 0 200 80 12 41 10
156 70 2 66 0 -182 90 2 62 0 143 90 8 60 0 105 80 11 66 0
103 80 5 38 0 250 70 8 53 10 100 60 13 37 10
test squamous
999 90 12 54 10 112 80 6 60 0 -87 80 3 48 0 -231 50 8 52 10
242 50 1 70 0 991 70 7 50 10 111 70 3 62 0 1 20 21 65 10
587 60 3 58 0 389 90 2 62 0 33 30 6 64 0 25 20 36 63 0
357 70 13 58 0 467 90 2 64 0 201 80 28 52 10 1 50 7 35 0
30 70 11 63 0 44 60 13 70 10 283 90 2 51 0 15 50 13 40 10
test small
25 30 2 69 0 -103 70 22 36 10 21 20 4 71 0 13 30 2 62 0
87 60 2 60 0 2 40 36 44 10 20 30 9 54 10 7 20 11 66 0
24 60 8 49 0 99 70 3 72 0 8 80 2 68 0 99 85 4 62 0
61 70 2 71 0 25 70 2 70 0 95 70 1 61 0 80 50 17 71 0
51 30 87 59 10 29 40 8 67 0
test adeno
24 40 2 60 0 18 40 5 69 10 -83 99 3 57 0 31 80 3 39 0
51 60 5 62 0 90 60 22 50 10 52 60 3 43 0 73 60 3 70 0
8 50 5 66 0 36 70 8 61 0 48 10 4 81 0 7 40 4 58 0
140 70 3 63 0 186 90 3 60 0 84 80 4 62 10 19 50 10 42 0
45 40 3 69 0 80 40 4 63 0
test large
52 60 4 45 0 164 70 15 68 10 19 30 4 39 10 53 60 12 66 0
15 30 5 63 0 43 60 11 49 10 340 80 10 64 10 133 75 1 65 0
111 60 5 64 0 231 70 18 67 10 378 80 4 65 0 49 30 3 37 0
;

proc phreg data=VALung;
class Prior(ref='no') Cell(ref='large') Therapy(ref='standard');
model Time*Status(0) = Kps Duration Age Cell Therapy;
store out=myModel;
run;

proc plm source = myModel;
effectplot fit(x=KPS);
run;

[SriniR]

---

@Reeza wrote:

I cannot seem to find a logical reason why it shouldn't work and it doesn't work for any variation so perhaps that option isn't supported for PROC PHREG though it doesn't say that anywhere as far as I can see. 

Maybe @Rick_SAS himself can comment 🙂

There's full code below to generate an example/test code.

Spoiler

proc format;
value yesno 0='no' 10='yes';
run;

data VALung;
drop check m;
retain Therapy Cell;
infile cards column=column;
length Check $ 1;
label Time='time to death in days'
Kps='Karnofsky performance scale'
Duration='months from diagnosis to randomization'
Age='age in years'
Prior='prior therapy'
Cell='cell type'
Therapy='type of treatment';
format Prior yesno.;
M=Column;
input Check $ @@;
if M>Column then M=1;
if Check='s'|Check='t' then do;
input @M Therapy $ Cell $;
delete;
end;
else do;
input @M Time Kps Duration Age Prior @@;
Status=(Time>0);
Time=abs(Time);
end;
datalines;
standard squamous
72 60 7 69 0 411 70 5 64 10 228 60 3 38 0 126 60 9 63 10
118 70 11 65 10 10 20 5 49 0 82 40 10 69 10 110 80 29 68 0
314 50 18 43 0 -100 70 6 70 0 42 60 4 81 0 8 40 58 63 10
144 30 4 63 0 -25 80 9 52 10 11 70 11 48 10
standard small
30 60 3 61 0 384 60 9 42 0 4 40 2 35 0 54 80 4 63 10
13 60 4 56 0 -123 40 3 55 0 -97 60 5 67 0 153 60 14 63 10
59 30 2 65 0 117 80 3 46 0 16 30 4 53 10 151 50 12 69 0
22 60 4 68 0 56 80 12 43 10 21 40 2 55 10 18 20 15 42 0
139 80 2 64 0 20 30 5 65 0 31 75 3 65 0 52 70 2 55 0
287 60 25 66 10 18 30 4 60 0 51 60 1 67 0 122 80 28 53 0
27 60 8 62 0 54 70 1 67 0 7 50 7 72 0 63 50 11 48 0
392 40 4 68 0 10 40 23 67 10
standard adeno
8 20 19 61 10 92 70 10 60 0 35 40 6 62 0 117 80 2 38 0
132 80 5 50 0 12 50 4 63 10 162 80 5 64 0 3 30 3 43 0
95 80 4 34 0
standard large
177 50 16 66 10 162 80 5 62 0 216 50 15 52 0 553 70 2 47 0
278 60 12 63 0 12 40 12 68 10 260 80 5 45 0 200 80 12 41 10
156 70 2 66 0 -182 90 2 62 0 143 90 8 60 0 105 80 11 66 0
103 80 5 38 0 250 70 8 53 10 100 60 13 37 10
test squamous
999 90 12 54 10 112 80 6 60 0 -87 80 3 48 0 -231 50 8 52 10
242 50 1 70 0 991 70 7 50 10 111 70 3 62 0 1 20 21 65 10
587 60 3 58 0 389 90 2 62 0 33 30 6 64 0 25 20 36 63 0
357 70 13 58 0 467 90 2 64 0 201 80 28 52 10 1 50 7 35 0
30 70 11 63 0 44 60 13 70 10 283 90 2 51 0 15 50 13 40 10
test small
25 30 2 69 0 -103 70 22 36 10 21 20 4 71 0 13 30 2 62 0
87 60 2 60 0 2 40 36 44 10 20 30 9 54 10 7 20 11 66 0
24 60 8 49 0 99 70 3 72 0 8 80 2 68 0 99 85 4 62 0
61 70 2 71 0 25 70 2 70 0 95 70 1 61 0 80 50 17 71 0
51 30 87 59 10 29 40 8 67 0
test adeno
24 40 2 60 0 18 40 5 69 10 -83 99 3 57 0 31 80 3 39 0
51 60 5 62 0 90 60 22 50 10 52 60 3 43 0 73 60 3 70 0
8 50 5 66 0 36 70 8 61 0 48 10 4 81 0 7 40 4 58 0
140 70 3 63 0 186 90 3 60 0 84 80 4 62 10 19 50 10 42 0
45 40 3 69 0 80 40 4 63 0
test large
52 60 4 45 0 164 70 15 68 10 19 30 4 39 10 53 60 12 66 0
15 30 5 63 0 43 60 11 49 10 340 80 10 64 10 133 75 1 65 0
111 60 5 64 0 231 70 18 67 10 378 80 4 65 0 49 30 3 37 0
;

proc phreg data=VALung;
class Prior(ref='no') Cell(ref='large') Therapy(ref='standard');
model Time*Status(0) = Kps Duration Age Cell Therapy;
store out=myModel;
run;

proc plm source = myModel;
effectplot fit(x=KPS);
run;

---

Thank you so much. I will try out the code. 

[SriniR]

I get the same message. Please see below:

1 OPTIONS NONOTES NOSTIMER NOSOURCE NOSYNTAXCHECK;

72

73 proc format;

74 value yesno 0='no' 10='yes';

NOTE: Format YESNO has been output.

75 run;

NOTE: PROCEDURE FORMAT used (Total process time):

real time 0.00 seconds

cpu time 0.01 seconds

76 data VALung;

77 drop check m;

78 retain Therapy Cell;

79 infile cards column=column;

80 length Check $ 1;

81 label Time='time to death in days'

82 Kps='Karnofsky performance scale'

83 Duration='months from diagnosis to randomization'

84 Age='age in years'

85 Prior='prior therapy'

86 Cell='cell type'

87 Therapy='type of treatment';

88 format Prior yesno.;

89 M=Column;

90 input Check $ @@;

91 if M>Column then M=1;

92 if Check='s'|Check='t' then do;

93 input @M Therapy $ Cell $;

94 delete;

95 end;

96 else do;

97 input @M Time Kps Duration Age Prior @@;

98 Status=(Time>0);

99 Time=abs(Time);

100 end;

101 datalines;

NOTE: SAS went to a new line when INPUT statement reached past the end of a line.

NOTE: The data set WORK.VALUNG has 137 observations and 8 variables.

NOTE: DATA statement used (Total process time):

real time 0.00 seconds

cpu time 0.01 seconds

147 ;

148

149 proc phreg data=VALung;

150 class Prior(ref='no') Cell(ref='large') Therapy(ref='standard');

151 model Time*Status(0) = Kps Duration Age Cell Therapy;

152 store out=myModel;

153 run;

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

NOTE: The PHREG procedure generated the model item store WORK.MYMODEL.

NOTE: PROCEDURE PHREG used (Total process time):

real time 0.14 seconds

cpu time 0.12 seconds

NOTE: Format YESNO is already on the library WORK.FORMATS.

NOTE: Format YESNO has been output.

NOTE: There were 2 observations read from the data set WORK.DATA1.

NOTE: PROCEDURE FORMAT used (Total process time):

real time 0.00 seconds

cpu time 0.00 seconds

154 proc plm source = myModel;

155 effectplot fit(x=KPS);

NOTE: The format YESNO was loaded from item store WORK.MYMODEL.

NOTE: The data set WORK.DATA1 has 2 observations and 21 variables.

156 run;

WARNING: The EFFECTPLOT statement cannot be used for the specified model.

NOTE: PROCEDURE PLM used (Total process time):

real time 0.03 seconds

cpu time 0.04 seconds

157

158 OPTIONS NONOTES NOSTIMER NOSOURCE NOSYNTAXCHECK;

170

[Reeza]

I said it didn't work. The code is for tech support or Rick to test to make it easier for them to diagnose the issue. It's beyond me at this point in time.

[SriniR]

---

@Reeza wrote:
I said it didn't work. The code is for tech support or Rick to test to make it easier for them to diagnose the issue. It's beyond me at this point in time.

---

My apologies. I should have read post properly. 

[Rick_SAS]

The message

WARNING: The EFFECTPLOT statement cannot be used for the specified model.

means that the EFFECTPLOT statement does not support the PHREG model. Create an output data set from PHREG

and use PROC SGPLOT to create the graph you want.

[Reeza]

Can the documentation be updated to reflect this?

Given that PROC PHREG is supposed to be supported in PLM but not for specific statement when there's no indication in the documentation otherwiset is confusing.

The following procedures support the STORE statement for the generation of item stores that can be processed with the PLM procedure: GEE, GENMOD, GLIMMIX, GLM, GLMSELECT, LIFEREG, LOGISTIC, MIXED, ORTHOREG, PHREG, PROBIT, SURVEYLOGISTIC, SURVEYPHREG, and SURVEYREG. The RELIABILITY procedure in SAS/QC software also supports the STORE statement. For details about the STORE statement, see the section STORE Statement of Chapter 19: Shared Concepts and Topics.

[Rick_SAS]

Yes, I have asked the developer to update the doc for the EFFECTPLOT statement.

In general, the EFFECTPLOT statement supports fixed-effects for generalized linear models.

After some basic testing, my colleagues and I think that the following procedures are NOT supported by the EFFECTPLOT statement:

• LIFEREG, PHREG, RMSTREG, and SURVEYPHREG.
• PROC LOGISTIC with a STRATA statement (conditional logistic regression)
• The effect plot for the MIXED procedure only supports the fixed effects.
• The RELIABILITY procedure in SAS/QC (but you can use the EFFECTPLOT statement in the procedure).

I do not know whether this list is comprehensive.