Baseline hazard function for a stratified Cox PH-model

[John Payne]

Hello everyone,

Does anyone know if it is possible to get the baseline hazard function (graph) from a stratified Cox PH-model?

I have fit a Cox-model using stcox, but two of the dummy variables are in violation of the proportional hazard assumption. From what I understand, this problem can be solved by stratfying the model based on these two dummy variables. This did in fact solve the PH assumption violation. However, I cannot recover the baseline hazard function anymore using "stcurve", which is the whole point of my study.

Here is the code for my Cox model after having stset my data:

stcox X1 X2 X3 X4 X5 X6 X7, vce(cluster CLUSTERID) nohr strata(X8 X9)

Does anyone know how/if it is possible to get this graph?

Any help is greatly appreciated.

Best regards,

John

[Steve Samuels]

sts graph will do this. Be sure to read the Manual description of the adjust() option.

[Andrea Discacciati]

Steve answered your question. However, if the "whole point of [your] study" is to estimate the (baseline) hazard function, then I would consider using other statistical models than the Cox PH model. More specifically, I am thinking of parametric models (e.g.: Poisson or Royston–Parmar models).

[John Payne]

Thank you very much guys, problem solved. Really appreciate your helpfulness.

I am in fact also estimating a parametric model with a Weibull distribution in order to compare the results for these two models. This brings me on to another question..

Is it possible to merge these two functions in the same graph? That is, I would like both baseline hazard functions in the same XY axis. Not two seperate graphs next to each other.

[Steve Samuels]

You can do it using addplot. Here's an example from the leukemia data set analyzed in the ST manual The curves are at the values wbc3cat=2. The dashed lines are from streg. . As Andrea suggests, it would be worthwhile to try the Royston-Parmar-Lambert stratified proportional hazards model, fit by Paul Lambert's stpm2 (SSC). The great advantage is that that you get a flexible parametric equation for the hazard function.

use http://www.stata-press.com/data/r14/leukemia, clear
tempfile t1 t2
stset weeks, failure(relapse)
save `t1', replace

/* streg hazard */
streg  treatment1 wbc3cat ,distribution(weibull)
stcurve, hazard ///
at(treatment1 =0 wbc3cat=2)  at(treatment1=1 wbc3cat=2) ///
outfile(`t2', replace)

use `t2', clear
append using `t1'  // add streg/stcurve results
gen wbcc2 = wbc3cat-2  /* Center at 2 */

/* stcox stratified hazard with streg added */
#delim ;
sts graph, hazard noboundary
strata(treatment1) adjust(wbcc2)

/* add streg hazards */
addplot(scatter haz1 haz2 _t if haz1<0.2,
ms(i i) c(l l) lpattern(shortdash longdash))
;
#delim cr
graph export gplot.png, replace

[John Payne]

Wow, that looks great. Thank you very much, Steve. I will try to do that for my models.

I also have two last questions, if I may.

I am considering adding another variable to my model. The variable is actually the (natural) log transformation of a number (the size of a company).

But when I add this variable to my model, the baseline hazard function completely flips. It goes from increasing with time (which makes theoretically sense) to decreasing with time (doesn't make empirically sense). So I am wondering if this has anything to do with the log transformation? If so, is there a way around this problem?

Lastly, I would just like to ask a (potentially stupid) question relating to the stratification.

What (if any) is the difference between estimating the baseline hazard function for a stratified model (strata is a dummy variable), compared to simply creating a subsample for this variable and then estimating the baseline hazard function?

[Steve Samuels]

Without seeing the commands, results, and graphs for the the two models (with and without added variable), it's difficult to say much. So follow the instructions of FAQ 12 for showing the commands and results and mimic the graph export statement in my post. The "attachment" button is to the left of the "A" in the upper-right of the editor window.. To answer your second question: in a stratified model, the coefficients for the non-stratum variables are the same in both strata though the hazard functions differ. If you subsample, the coefficients for the same variable in one stratum are estimated independently of the variables. I suggest that you consult a survival text. Stata Press sells a good one (An Introduction to Survival Analysis using Stata).

[Steve Samuels]

John, start a new topic with your first question, as it is far afield of the original.

[John Payne]

Hi again, Steve.

I managed to add the two graphs together following your suggestion above. However, Stata won't allow me to produce the graphs for a stratified Weibull model in the same way as for the Cox-model.

Instead I had to use two different methods for each model:

stcox X1 X2 X3, strata(X4)
sts graph, strata(X4) adjustfor(X1 X2 X3)

streg X1 X2 X3 X4, dist(weibull)
stcurve, hazard at1(X1=0 X2=0 X3=0 X4=0) at2(X1=0 X2=0 X3=0 X4=1)

And from there add the two graphs together.

If my aim is to compare the Cox baseline hazard function with the Weibull baseline hazard function, this is not really the same is it? Because the Weibull model is not stratified..?

stcurve won't allow me to recover the baseline hazard function for a stratified model, and sts graph followed by "strata(X4) adjustfor(X1 X2 X3)" seems to produce a graph that doesn't look like a Weibull distribution (it looks like the Cox function).

Hope someone can help me with this.

Thanks,

John

[Steve Samuels]

Sorry, I didn't use the stratified Weibull in my previous example. Since I rarely use the standard parametric models, I wasn't aware of the strata option. To get the hazard function, you'll need to center all non-strata covariates to a meaningfull value. Then compute the estimated hazard from the model parameters. See page 247 of the ST Manual (Stata 14).
I've compared the results to the non-stratified Weibull above and found the curves to be nearly identical.

To answer your last concern: sts graph produces adjusted curves by fitting a Cox model internally. It doesn't have to follow a prior model command.

use http://www.stata-press.com/data/r14/leukemia,clear
gen wbcc2 = wbc3cat-2  /* Center covariate (at 2) */

stset weeks, failure(relapse)

streg   wbcc2   ,distribution(weibull) strata(treatment1)
scalar b1  =  _b[_t:_cons]
scalar b2 =   b1 + _b[_t:_Streatment_1]
scalar p1 =  exp(_b[ln_p:_cons])
scalar p2 =  p1*exp(_b[ln_p:_Streatment_1])

gen haz1 = p1*exp(b1)*_t^(p1-1)
gen haz2 = p2*exp(b2)*_t^(p2-1)

/* stcox stratified hazard  */
#delim ;
sts graph, hazard noboundary
strata(treatment1) adjust(wbcc2)
plotopts(lpattern(solid))
title("Stratified Weibull & Cox")

/* add streg hazards */
addplot(scatter haz1 haz2 _t if haz1<0.2,
ms(i i) c(l l) lpattern(longdash longdash))
;
#delim cr
graph export gplot.png, replace

[John Payne]

Once again, thank you very much Steve!

Problem solved. Your help is greatly appreciated.