Generating and interpreting hazard ratios from Cox models using lincom and margins commands

[SL Navarro]

Hi all -

There are some similar posts about this topic but none that specifically answer my question. Using Stata v14.1, I ran a Cox model which included an interaction term for two continuous variables. This term was significant. I followed up with a lincom command to get the hazard ratio for the interaction, but this was not significant. I am not sure why one is significant and the other is not. I'm assuming I must be doing something wrong with my code. I also tried using the margins command as suggested in another post, but this yielded an error message, "factor variables may not contain noninteger values". My variables do not have any noninteger values. Any help sorting this out would be much appreciated.

stcox EPAplusDHA f60solfb age bmi smoking educ f60alc tothormonestat texpwk ///
> f60enrgy f60fldeq f60calc redmeat dmarm hrtarm cadarm colorel c.EPAplusDHA#c.f60solfb

failure _d: colorectal == 1
analysis time _t: time
id: id

Iteration 0: log likelihood = -22630.355
Iteration 1: log likelihood = -22385.838
Iteration 2: log likelihood = -22385.147
Iteration 3: log likelihood = -22385.093
Iteration 4: log likelihood = -22385.092
Iteration 5: log likelihood = -22385.092
Refining estimates:
Iteration 0: log likelihood = -22385.092

Cox regression -- Breslow method for ties

No. of subjects = 134,017 Number of obs = 134,017
No. of failures = 1,952
Time at risk = 1569353.803
LR chi2(18) = 490.53
Log likelihood = -22385.092 Prob > chi2 = 0.0000

-----------------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
------------------------+----------------------------------------------------------------
EPAplusDHA | .5044889 .2055583 -1.68 0.093 .2269981 1.121195
f60solfb | .9652983 .0208153 -1.64 0.101 .9253512 1.00697
age | 1.063833 .0035574 18.50 0.000 1.056883 1.070828
bmi | 1.024562 .0043645 5.70 0.000 1.016044 1.033152
smoking | 1.103438 .0370171 2.93 0.003 1.033219 1.178428
educ | 1.021955 .0308403 0.72 0.472 .9632625 1.084225
f60alc | .9986251 .0022477 -0.61 0.541 .9942293 1.00304
tothormonestat | .7813583 .0367145 -5.25 0.000 .7126133 .8567351
texpwk | .9965047 .0018905 -1.85 0.065 .9928063 1.000217
f60enrgy | 1.000179 .0000817 2.19 0.028 1.000019 1.000339
f60fldeq | .9997964 .0001836 -1.11 0.267 .9994367 1.000156
f60calc | .9997787 .0000787 -2.81 0.005 .9996245 .999933
redmeat | .9957916 .0571086 -0.07 0.941 .8899222 1.114256
dmarm | 1.015401 .0315733 0.49 0.623 .9553661 1.079208
hrtarm | 1.013579 .0219306 0.62 0.533 .971494 1.057486
cadarm | .9951817 .0420857 -0.11 0.909 .9160212 1.081183
colorel | 1.055499 .0367989 1.55 0.121 .9857837 1.130145
|
c.EPAplusDHA#c.f60solfb | 1.133301 .068782 2.06 0.039 1.0062 1.276457
-----------------------------------------------------------------------------------------

lincom EPAplusDHA + f60pect + c.EPAplusDHA#c.f60solfb, hr

( 1) EPAplusDHA + f60solfb + c.EPAplusDHA#c.f60solfb = 0

------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
(1) | .5518976 .199766 -1.64 0.101 .2714913 1.121918
------------------------------------------------------------------------------

margins EPAplusDHA#f60solfb

EPAplusDHA: factor variables may not contain noninteger values
r(452);

Sandi

[Clyde Schechter]

First, you have done some editing to your code and output, as they are not consistent with each other. Your -lincom- command mentions a variable f60pect which is not shown in your output either from the Cox regression or in the -lincom- output.

Next, assuming that you really meant f60solfb instead of f60pect, that command will give you the hazard ratio conditional on f60solfb == 1 & EPAplusDHA == 1. If that's what you're trying to get, that's fine. But I have to say that without even knowing anything about your research goals, I suspect it isn't really what you want because it's just pretty odd and arbitrary when the two variables involved are both continuous.

The -margins- command appropriately rejected your request because that syntax for -margins- is only available for discrete variables. If you would like to know the hazard ratio for specific designated values of f60solfb and EPAplusDHA, the correct syntax is:

margins, at(f60solfb=(list of interesting values for f60solfb) EPAplusDHA = (list of interesting values for EPAplusDHA))

and you will get the hazard ratios corresponding to each combination of the f60solfb and EPAplusDHA values you specify.

Now, the other thing to remember with continuous variable interactions is that you have a different hazard ratio for every combination of f60solfb and EPAplusDHA. No matter how strongly "significant" the interaction term is in the regression, there will, in general, be values of f60solfb and EPAplusDHA for which the hazard ratio at those values is not significant. So the problem is that your expectation about the outcome of -lincom- based on the coefficient of the interaction term is unfounded.

And what I said about interaction terms in the preceding paragraph applies to interactions of dichotomies as well. It is entirely possible to have a highly significant A#B coefficient, and yet the effects of A at B= 0 and at B=1 are both non-significant. (Or vice-versa switching the roles of A and B.) The difference between two non-significant things can be significant. The difference between two signifcant things can be non-significant. And the difference between a significant and a non-significant thing can be significant or not.

I recommend you not focus on significance. You end up in all kinds of paradoxes because the non-technical meaning of the word leads you down the garden path to conclusions that are incorrect. When I read regression outputs, I generally ignore the p-value column and focus on my coefficients (or hazard ratios or odds ratios, etc.) and their magnitudes, and then look at the confidence intervals or the standard errors to get a sense of the precision of my results. If I have some time later on when I have nothing better to do, I might look at the p-values. Remember that p-values test a null hypothesis (coef = 0, or hr = 1) that is almost always a straw man. Their just not very useful, and you have to be super-careful not to get tricked by them.

[SL Navarro]

Thanks for the quick reply. I did edit my output as I realized I had copied and pasted the wrong output, but the values are correct. Should be f60solfb.

So if I'm understanding correctly, lincom isn't the best option to get a hazard ratio for an interaction for two continuous variables. I don't have any "interesting values" for these variables in which to use in the margins command as they are continuous unless I arbitrarily specify an intake threshold (as these are dietary variables). Is there any way to get an overall hazard ratio for an interaction of two continuous variables? Obviously 1.13 doesn't make sense when each individual variable HR is <1. Is this directly interpretable and is the over all HR truly >1?

[Clyde Schechter]

Is there any way to get an overall hazard ratio for an interaction of two continuous variables?

No. Because no such thing exists. You won't find it in Stata, nor anywhere else.

Obviously 1.13 doesn't make sense when each individual variable HR is <1. Is this directly interpretable and is the over all HR truly >1?

Yes, it does make sense. You can have one hazard ratio that is, say 0.75 and another that is 0.8475. Both of those are less than 1, but their ratio is 1.13.

The point is that the "hazard ratio" for an interaction term shown in the -stcox- output is not actually a hazard ratio: it is a ratio of hazard ratios. When you use an interaction model, there is no such thing as "the hazard ratio" for the interaction. There is a ratio of hazard ratios given in the -stcox- output, though it is a rather odd construct that is a bit difficult for most people to grasp. It is the ratio by which the hazard ratio for f60solfb increases in association with a unit increase in EPAplusDHA (or vice versa).

[SL Navarro]

Very helpful, thank you!