5

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

Lily Clements
  • 115
  • 1
  • 6
  • I don't see anything in `survminer` that addresses this issue. Have you considered using a package that does have plotting functions for hazard estimates? Look at output from `??"hazard"` as a start. – IRTFM Jul 01 '19 at 03:42
  • There's currently no packages that plots the hazard function (I looked). Your code is seems to work though. – Kristine Sep 03 '19 at 07:52

2 Answers2

5

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)
Brent
  • 4,611
  • 4
  • 38
  • 55
3

In the rms package, you can using the survplot function with "what" parameter specified as "hazard" to plot the hazard function verse time.

https://rdrr.io/cran/rms/man/survplot.html

Yellow Bridge
  • 31
  • 3