Scale of y-axis of a plot showing the predicted outcome risk for a continious predictor in a logistic regression model

Question

I was previously informed in another thread that my plot on the y-axis is showing odds and not odds ratios. The following example code and figure depict this.

library(survival)
data(cancer, package="survival")
d <- colon[, c("age", "status")]
table(d$status)

dd <- datadist(d); options(datadist='dd')
fit <- lrm(status ~ rcs(age, 4), data=d)
predictions <- Predict(fit, age, ref.zero=TRUE, fun=exp)
ggplot(predictions) + ylab("Predicted odds")

In my understanding: The odds are defined as the probability that the event will occur divided by the probability that the event will not occur. With regard to my example, the overall odds would be 0.98.

> table(d$status)

  0   1 
938 920 
> round(920/938, 2)
[1] 0.98

Based on this, I do understand that a single point in the plot above shows the odds for a specific age.

From this understanding, I would assume that the scale on the y-axis would be hazard and not hazard ratio for a time to event outcome. However, the two following screenshots show figures of published papers with hazard ratios on the y-axis. What do I get wrong?

Johannesen CDL, et al. BMJ. 2020.

Rawshani A, et al. N Engl J Med. 2018.

Update after Björn's answer

Thanks for your detailed answer! Let me rephrase it in my own words to see if I understood it correctly.

1. When the plot shows odds or probability, the value on the Y-axis must be interpreted separately for each value of x, and the 95% CI should not be zero at any point (case #1).
2. When the plot shows odds ratio (or possibly risk ratio), the value on the Y-axis must be interpreted in relation to the reference x value; the y value should be equal to 1 for the reference x value, and the 95% CI should be zero at this point (case #2).

I tried to demonstrate this using examples with the rms package and the spline package. I believe that I did it correct using the rms package, but I am not sure how to do it (case #2) with the spline package.

### Case #1 using rms
library(survival)
data(cancer, package="survival")
d <- colon[, c("age", "status")]
table(d$status)

dd <- datadist(d); options(datadist='dd')
# dd$limits["Adjust to","age"] <- 60
fit <- lrm(status ~ rcs(age, 4), data=d)
predictions <- Predict(fit, age, ref.zero=FALSE, fun=exp)
ggplot(predictions) + ylab("Predicted Odds")

### Case #2 using rms
library(survival)
data(cancer, package="survival")
d <- colon[, c("age", "status")]
table(d$status)

dd <- datadist(d); options(datadist='dd')
dd$limits["Adjust to","age"] <- 60
fit <- lrm(status ~ rcs(age, 4), data=d)
predictions <- Predict(fit, age, ref.zero=TRUE, fun=exp)

age_range <- data.frame(age = seq(from = 27, to = 81, length = 100))

ggplot(predictions) + ylab("Predicted Odds Ratio")

### Case #1 using spline
library(survival)
data(cancer, package="survival")
d <- colon[, c("age", "status")]

fit <- glm("status ~ ns(age, 3)", family=binomial(link="logit"), data=d)
age_range <- data.frame(age = seq(from = 27, to = 81, length = 100))


predictions <- predict(fit, newdata = age_range, se.fit=TRUE, type = "link", level = 0.95)
x <- data.frame(age = age_range$age,
                 effect = exp(predictions$fit),
                 lower = exp(predictions$fit - 1.96 * predictions$se.fit),
                 upper = exp(predictions$fit + 1.96 * predictions$se.fit),
                 model = "ns")
ggplot(data = x, aes(x = age, y = effect)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
  xlab("Age") +
  ylab("Predicted Odds")

Answer 1

There are two different choices on how to display these kind of data, which I will first illustrate with you first odds / odds-ratio example

1. Give the predicted odds ($\text{odds at age X}:=p(\text{event}|\text{age}=X)/p(\text{no event}|\text{age}=X)$) for each predictor value.
2. Give the odds ratio compared with a reference value (e.g. compared with the odds at age 50), which would be defined as $\text{odds ratio for age } X \text{ compared with age }50 := \text{odds at age X} / \text{odds at age 50}$.

A a binary outcome you can relatively easily produce both displays. I would also be easy to do this for an exponential time-to-event model, where a single hazard rate describes the distribution, and a hazard ratio can be calculated as the ratio of the hazard rate over the hazard rate for a reference category. However, many time-to-even analyses use Cox regression (but similar considerations as below apply to other time-to-event models like Weibull regression, too), where the the baseline hazard function over time is a non-parametrically estimated function that cannot easily be summarized into a single number. I.e. when using Cox regression you usually have \begin{equation} \text{hazard rate}(t | age = X) = \text{baseline hazard function}(t) \times e^{ \text{intercept} + \beta \times \text{age}}. \end{equation} There might also be other terms in the regression equation (other than intercept and age), which would complicate things further (especially if those covariates are correlated with age). At best you could plot the hazard at a fixed time $t$ you need to choose. However, once you look at a hazard ratio compared vs. age 50, you get ride of the baseline hazard function \begin{equation} \text{hazard ratio for age X compared with age 50} := \text{hazard rate}(t | age = X) / \text{hazard rate}(t | age = 50) = e^{ \beta \times (\text{age} - 50)}. \end{equation}

That's why plotting a hazard ratio relative to some reference value is a logical choice in a time-to-event setting.

You notice that this is what was done by noting that there is some value where the curves intersect the reference line for a hazard ratio of $1.0$ and where the confidence interval has width zero. That would be the reference value that was chosen.