Illustrating the hazard function across values of an independent continuous varable

I'm running a Cox model in which my main interest is the effect of a continuous independent variable on the hazard of committing crime.

I want to illustrate the effect of this variable. For this purpose, I produced the hazard function based on 5 values of this variable, and got a graph with 5 lines, one for each value, where the x-axis is time and the y-axis is the hazard function. However, since I'm not interested in the change of hazard over time, but rather in the difference in hazard between the 5 values of the independent variable, I calculated the mean of the hazard for each line (i.e., the values across all time points). Just to illustrate, this gives me, for example, that the mean hazard for a value of 0 in the independent variable is 0.5, whereas it is 1.5 for a value of 1.

Does that make any sense to report my results this way? I'm open to other suggestions.

Thanks!

It's tricky to try to work with the hazard per se in a Cox model, as it isn't directly estimated. The Cox model's regression coefficient for a predictor is the log-hazard difference versus its reference value. After the model is fit you can then estimate the hazard function for a set of reference predictor values and adjust it for other predictor values, but that hazard will be exactly 0 between event times.

What you seem to be doing it to calculate an average hazard over some period of time (cumulative hazard at time $$t$$ divided by $$t$$) as a function of the continuous predictor. I'd be reluctant to do that, as that depends on the choice of time span, with cumulative hazards unbounded over time. See Wikipedia.

What's typically done in this situation is to start with reference values for all predictors in the model.* For example, choose the average value of your continuous predictor as its reference, and choose representative values of the other predictors.

Then make a set of predictions (with confidence intervals) from the model, relative to that set of reference values, over a range of values of your continuous predictor of interest. You can then plot the linear predictor (log-hazard) from the model or exponentiate it to get the hazard ratio with respect to the reference condition.

One warning: if you include only a single term for a continuous predictor in your model, the model assumes that the log-hazard is linearly associated with the predictor value. That's often not a good assumption. If you did that, then the plot I recommend will be a straight line on the log-hazard scale.

You should consider modeling a continuous predictor flexibly, for example via a regression spline. See Chapter 2 of Frank Harrell's Regression Modeling Strategies, for example.

*Software for fitting Cox models chooses a set of reference values, but the default choices might not necessarily make sense for your application.

• Dear @EdM, thank you for the detailed answer. Would it make sense to calculate the Nelson–Aalen cumulative hazard function for different values of the independent variable, and then report the cumulative hazard at the last time point? This way I can say, for example, that the cumulative hazard to commit crime is 0.20 for those with a value of X in the independent variable, but 0.30 for those with a value of Y. Hope I'm being clear enough.
– Eran
Commented Dec 11, 2023 at 18:02
• @Eran you could do that, but it would be easier for others to understand if you reported the event probability as a function of the continuous predictor value at that last time point. That comes from the Nelson-Aalen cumulative hazard, anyway. Cumulative hazards are hard to grasp for many people; event probabilities are much easier. You will have to choose and specify in your report the values of the other predictors, if you want to report absolute hazards or probabilities.
– EdM
Commented Dec 11, 2023 at 22:23
• Great, thanks @EdM. If I understand correctly, in STATA at least, these are calculated not from the Cox regression estimates. So if I report them, would it make more sense to do that before I show the Cox results, or after?
– Eran
Commented Dec 12, 2023 at 3:16
• @Eran survival probability predictions start with fitting a Cox model, using the coefficient estimates to get the baseline cumulative hazard (Nelson-Aalen-Breslow), adjusting the cumulative hazard for specific predictor values, and converting cumulative hazard to survival probability. See this page, for example. In that sense the probabilities aren't calculated "from the Cox regression estimates," but they come from calculations starting from the fitted model. Report the Cox model and then the probability estimates derived from it.
– EdM
Commented Dec 12, 2023 at 13:20

Instead of looking at the hazard, you could use the survival probability using the contsurvplot package (https://robindenz1.github.io/contsurvplot/). For example, you could use the plot_surv_area() function to create a survival area plot, which visualizes the survival probability as a function of time and the continuous variable. These plots look like this:

In the associated paper we go into more depth and showcase some alternatives. You can find the paper here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10392888/