Proportionality in Cox model

[Reza Haririan]

I have a cohort and looking for association of a continuous variable of interest with the outcome. When I include a binary confounder along the predictor of interest in the model, proportionality is not violated by Schoenfeld test, and the predictor of interest is significantly associated with outcome. However, when running the Cox model separately the association is significant for only one of the binary confounders (albeit the smaller subcohort of the two). And the pattern of smoothed hazard estimates over time are not similar. Stratifying the model by the binary confounder, the model still shows a significant association with the outcome. I cannot reconcile these contradictory relationships, and I'd appreciate any input what to do in this situation and why.

[Clyde Schechter]

Well, where to begin? At the highest level, you have subjected your data to hypothesis testing using three radically different models, so you should not be surprised that their results are not as similar as you had hoped.

That said, they may or may not actually be as dissimilar as you fear: the difference between statistically significant and not statistically significant results is, itself, not statistically significant. In fact, it isn't even meaningful. When you break up a single data set into two subsets and even run the same analysis (let alone when you run a radically different one), it is very easy for a "significant" result in the whole set to evaporate just due to sample size considerations. So comparing statistical significance of a result in the whole with that of a subset is something you should do only if you have previously supported doing that with a careful statistical power analysis. In any case, you should also look specifically at the coefficients you got in the two subset analyses and see whether they are actually very different from each other numerically: they may or may not be.

But let me return to the issue of the three very different models you have used, as I think this is a more important point. Let's call your continuous predictor of interest C, and let's call the dichotomous confounder d. In your original analysis where d is included as a covariate in the model, the underlying model is:

h(C, d) = h0*(bcC + bdd)  (Model 0)

where h is the hazard function for your outcome, b_c and b_d are the coefficients of C and d, and h₀ is the Baseline Hazard function. By baseline hazard function I mean the hazard function conditional on C = 0 and d = 0. Key observation: this analysis entails one coefficient for C and one baseline hazard function.

Now, when you break this up into two separate analyses, one for d = 0 and the other for d = 1, the two models are:

h(C | d = 0) = h00*(bc0C)  (Model 1)

h(C | d = 1) = h01*(bc1C)  (Model 2)

Notice that in addition to getting coefficients b_c1 and b_c1, you also get two separate baseline hazard functions, h_{00 and} h₀₁, which are the hazard functions when C = 0 and d = 0 in the first model, and the hazard functions when C = 0 and d = 1 in the second. It is the fact that this pair of models allows the baseline C = 0 hazard to differ depending on the value of d, in contrast to the single hazard function h₀ that applies in Model 0 that accounts for the fact that the smoothed hazards you get in Models 1 and 2 do not resemble what you saw in Model 0. The difference between Model 0 and Models 1 and 2, is very much like the difference between having d as a covariate in a simple linear regression (Model 0) and having a d X time interaction in the model (Models 1 and 2). These are very different models. Key observation: these models entail two separate baseline hazard functions and two separate coefficients for C.

Now let's look at what the model underlying the stratified analysis is.

h(C | d) = cond(d == 0, h#00, h#01)*(bcC)  (Model 3)

In model 3, just as in Models 1 and 2, we have two baseline hazards in play. (I use the superscripted # only so that the notation reflects the fact that the two baseline hazards in Model 3 are, in general, not the same as the ones in Models 1 and 2.) But, unlike Models 1 and 2 which gave two separate coefficients of C, one for each value of d, here we have a single coefficient of C that is applied regardless of the value of d. Key observation: this model entails two separate baseline hazard functions and a single coefficient for C.

So as you can see, these models are very, very different from each other, and there is no reason to expect them to give you similar results.

This leaves you with the dilemma of which set of results you should rely on. In the absence of good prior information to support the choice of one approach over the others, your best bet here is to check the goodness of fit of these models with the -estat gofplot- command.

[Reza Haririan]

Thank you very much for your detailed explanation. You can see the GOF plots for the 3 models here and next post. They look similar to me!! The N=51 for the d=1 group (HR 0.90, P=0.02) and N=219 for d==0 (HR 0.97, P=0.14). What would you recommend? Thanks again

[Reza Haririan]

It seems that the pasted graphs do not show.

[Clyde Schechter]

On the edit window's toolbar, in the second row, the third icon from the left is a button for inserting an image. Don't use it with .gph files. Export the graphs to another format such as .png and then insert them.

[Reza Haririan]

Using "Paste from Word" option:

[Reza Haririan]

Stratified

[Reza Haririan]

1

[Reza Haririan]

Please see the graphs in attached word document. Thank you

[Stephen Jenkins]

Reza Haririan Clyde asked you to post the graphs in png format here. That is what the forum FAQ refers to. (Use the "image" button to upload the files.) Please read the FAQ; please don't post Word documents.

[Reza Haririan]

gof 1

[Reza Haririan]

Gof 1

[Reza Haririan]

gof3

[Reza Haririan]

gof0/2

[Reza Haririan]

gof 1/2