CODEWITHSUNDEEP
rsurvival
Lily Clements
Lily Clements

Reputation: 125

Plotting the Hazard Function in R using survminer?

We can use survminer to plot the survival function or cumulative hazard function, but I cannot see a way to use it to plot the hazard function.

For example,

library(survival)
library(tidyverse)
library(survminer)

data(lung)

# Run Kaplan-Meier on the data
mod.lung <- survfit(Surv(time, status) ~ 1, data = lung)

# Kaplan-Meier Survival Curve
ggsurvplot(mod.lung)

# Cumulative Hazard
ggsurvplot(mod.lung, fun = function(y) -log(y))

Since the cumulative hazard function is H(t) = -log(S(t)) then I just need to add in fun = function(y) -log(y) to get the cumulative hazard plot.

The hazard function is h(t) = -d/dt log(S(t)), and so I am unsure how to use this to get the hazard function in a survminer plot.

An alternative definition of the hazard function is h(t) = f(t)/S(t), however, I'm unsure how to use this to get the plot.

I have found ways to get a hazard plot using ggplot2, for example

survival.table1 <- broom::tidy(mod.lung) %>% filter(n.event > 0)
survival.table1 <- survival.table1 %>% mutate(hazard = n.event / (n.risk * (lead(time) - time)))
ggplot() +
  geom_step(data = survival.table1, aes(x = time, y = hazard)) +
  labs(x = "Time", y = "Hazard")

However, I mainly wish to find a way with the survminer package, partly to have some consistency.

Thanks

Upvotes: 6

Views: 11602

Answers (1)

Brent
Brent

Reputation: 4906

A couple of years late, but here for others.

survplot can only be used to plot the hazard if the estimate was created by the psm function. The psm function of the rms library fits the accelerated failure time family of parametric survival models (defaulting to a Weibull distribution). Other available distributions are in the documentation for the survreg package:

These include "weibull", "exponential", "gaussian", "logistic","lognormal" and "loglogistic". Otherwise, it is assumed to be a user defined list conforming to the format described in survreg.distributions

library(rms)
mod.lung <- psm(Surv(time, status) ~ 1, data = lung)
survplot(mod.lung, what="hazard")

For non parametric survival models, the muhaz library might be more useful. This example uses the default "epanechnikov" boundary kernel function. You may wish to explore different bandwidth options - see the muhaz package documentation.

library(muhaz)
mod.lung <- muhaz(lung$time, lung$status - 1) # status must be 1 for failure and 0 for censored
plot(mod.lung)

Alternatively, to apply B-splines instead of kernel density smoothing to the hazard, have a look at the bshazard library

library(bshazard)
mod.lung <- bshazard(Surv(time, status) ~ 1, data = lung)
plot(mod.lung)

Upvotes: 5

Related Questions