Skip to contents

This helper function constructs a ggplot2 geom to plot Kaplan-Meier curves from a network containing survival or time-to-event outcomes. This is useful for overlaying the "raw" survival data on the estimated survival functions created with plotted with plot.surv_nma_summary(), but can also be used standalone to plot Kaplan-Meier curves before fitting a model.

Usage

geom_km(
  network,
  ...,
  transform = c("identity", "cloglog", "log", "cumhaz"),
  curve_args = list(),
  cens_args = list()
)

Arguments

network

A nma_data network object containing survival outcomes

...

Additional arguments passed to survival::survfit()

transform

Character string giving the transformation to apply to the KM curves before plotting. The default is "identity" for no transformation; other options are "cloglog" for \(\log(-\log(S))\), "log" for \(\log(S)\), or "cumhaz" for the cumulative hazard \(-\log(S)\).

curve_args

Optional list of arguments to customise the curves plotted with ggplot2::geom_step()

cens_args

Optional list of arguments to customise the censoring marks plotted with ggplot2::geom_point()

Value

A ggplot2 geom list that can be added to a ggplot2 plot object

Examples

# Set up newly-diagnosed multiple myeloma network

head(ndmm_ipd)
#>          study trt       studyf trtf      age iss_stage3 response_cr_vgpr male
#> 1 McCarthy2012 Pbo McCarthy2012  Pbo 50.81625          0                1    0
#> 2 McCarthy2012 Pbo McCarthy2012  Pbo 62.18165          0                0    0
#> 3 McCarthy2012 Pbo McCarthy2012  Pbo 51.53762          1                1    1
#> 4 McCarthy2012 Pbo McCarthy2012  Pbo 46.74128          0                1    1
#> 5 McCarthy2012 Pbo McCarthy2012  Pbo 62.62561          0                1    1
#> 6 McCarthy2012 Pbo McCarthy2012  Pbo 49.24520          1                1    0
#>   eventtime status
#> 1 31.106516      1
#> 2  3.299623      0
#> 3 57.400000      0
#> 4 57.400000      0
#> 5 57.400000      0
#> 6 30.714460      0
head(ndmm_agd)
#>        study     studyf trt trtf eventtime status
#> 1 Morgan2012 Morgan2012 Pbo  Pbo  18.72575      1
#> 2 Morgan2012 Morgan2012 Pbo  Pbo  63.36000      0
#> 3 Morgan2012 Morgan2012 Pbo  Pbo  34.35726      1
#> 4 Morgan2012 Morgan2012 Pbo  Pbo  10.77826      1
#> 5 Morgan2012 Morgan2012 Pbo  Pbo  63.36000      0
#> 6 Morgan2012 Morgan2012 Pbo  Pbo  14.52966      1

ndmm_net <- combine_network(
  set_ipd(ndmm_ipd,
          study, trt,
          Surv = Surv(eventtime / 12, status)),
  set_agd_surv(ndmm_agd,
               study, trt,
               Surv = Surv(eventtime / 12, status),
               covariates = ndmm_agd_covs))
# Plot KM curves using ggplot2
library(ggplot2)

# We need to create an empty ggplot object to add the curves to
ggplot() + geom_km(ndmm_net)


# Adding plotting options, facets, axis labels, and a plot theme
ggplot() +
  geom_km(ndmm_net,
          curve_args = list(linewidth = 0.5),
          cens_args = list(size = 3, shape = 124)) +
  facet_wrap(vars(Study)) +
  labs(xlab = "Time", ylab = "Survival Probability") +
  theme_multinma()


# Using the transform argument to produce log-log plots (e.g. to assess the
# proportional hazards assumption)
ggplot() +
  geom_km(ndmm_net, transform = "cloglog") +
  facet_wrap(vars(Study)) +
  theme_multinma()


# Using the transform argument to produce cumulative hazard plots
ggplot() +
  geom_km(ndmm_net, transform = "cumhaz") +
  facet_wrap(vars(Study)) +
  theme_multinma()


# This function can also be used to add KM data to plots of estimated survival
# curves from a fitted model, in a similar manner
# \donttest{
# Run newly-diagnosed multiple myeloma example if not already available
if (!exists("ndmm_fit")) example("example_ndmm", run.donttest = TRUE)
#> 
#> exmpl_> # Set up newly-diagnosed multiple myeloma network
#> exmpl_> 
#> exmpl_> head(ndmm_ipd)
#>          study trt       studyf trtf      age iss_stage3 response_cr_vgpr male
#> 1 McCarthy2012 Pbo McCarthy2012  Pbo 50.81625          0                1    0
#> 2 McCarthy2012 Pbo McCarthy2012  Pbo 62.18165          0                0    0
#> 3 McCarthy2012 Pbo McCarthy2012  Pbo 51.53762          1                1    1
#> 4 McCarthy2012 Pbo McCarthy2012  Pbo 46.74128          0                1    1
#> 5 McCarthy2012 Pbo McCarthy2012  Pbo 62.62561          0                1    1
#> 6 McCarthy2012 Pbo McCarthy2012  Pbo 49.24520          1                1    0
#>   eventtime status
#> 1 31.106516      1
#> 2  3.299623      0
#> 3 57.400000      0
#> 4 57.400000      0
#> 5 57.400000      0
#> 6 30.714460      0
#> 
#> exmpl_> head(ndmm_agd)
#>        study     studyf trt trtf eventtime status
#> 1 Morgan2012 Morgan2012 Pbo  Pbo  18.72575      1
#> 2 Morgan2012 Morgan2012 Pbo  Pbo  63.36000      0
#> 3 Morgan2012 Morgan2012 Pbo  Pbo  34.35726      1
#> 4 Morgan2012 Morgan2012 Pbo  Pbo  10.77826      1
#> 5 Morgan2012 Morgan2012 Pbo  Pbo  63.36000      0
#> 6 Morgan2012 Morgan2012 Pbo  Pbo  14.52966      1
#> 
#> exmpl_> ndmm_net <- combine_network(
#> exmpl_+   set_ipd(ndmm_ipd,
#> exmpl_+           study, trt,
#> exmpl_+           Surv = Surv(eventtime / 12, status)),
#> exmpl_+   set_agd_surv(ndmm_agd,
#> exmpl_+                study, trt,
#> exmpl_+                Surv = Surv(eventtime / 12, status),
#> exmpl_+                covariates = ndmm_agd_covs))
#> 
#> exmpl_> ## No test: 
#> exmpl_> # Fit Weibull (PH) model
#> exmpl_> ndmm_fit <- nma(ndmm_net, ## Don't show: 
#> exmpl_+ refresh = if (interactive()) 200 else 0,
#> exmpl_+ ## End(Don't show)
#> exmpl_+                 likelihood = "weibull",
#> exmpl_+                 prior_intercept = normal(scale = 100),
#> exmpl_+                 prior_trt = normal(scale = 10),
#> exmpl_+                 prior_aux = half_normal(scale = 10))
#> Note: Setting "Pbo" as the network reference treatment.
#> 
#> exmpl_> ndmm_fit
#> A fixed effects NMA with a weibull likelihood (log link).
#> Inference for Stan model: survival_param.
#> 4 chains, each with iter=2000; warmup=1000; thin=1; 
#> post-warmup draws per chain=1000, total post-warmup draws=4000.
#> 
#>                         mean se_mean   sd     2.5%      25%      50%      75%
#> d[Len]                 -0.54    0.00 0.05    -0.63    -0.57    -0.54    -0.51
#> d[Thal]                -0.11    0.00 0.09    -0.29    -0.17    -0.11    -0.05
#> lp__                -6229.99    0.06 2.49 -6235.81 -6231.42 -6229.61 -6228.18
#> shape[Attal2012]        1.30    0.00 0.06     1.18     1.25     1.30     1.34
#> shape[Jackson2019]      0.93    0.00 0.02     0.89     0.92     0.93     0.95
#> shape[McCarthy2012]     1.30    0.00 0.07     1.17     1.25     1.29     1.34
#> shape[Morgan2012]       0.88    0.00 0.03     0.82     0.86     0.88     0.90
#> shape[Palumbo2014]      1.02    0.00 0.07     0.88     0.97     1.02     1.06
#>                        97.5% n_eff Rhat
#> d[Len]                 -0.45  5031    1
#> d[Thal]                 0.06  4915    1
#> lp__                -6226.16  1548    1
#> shape[Attal2012]        1.42  4275    1
#> shape[Jackson2019]      0.98  4679    1
#> shape[McCarthy2012]     1.43  4275    1
#> shape[Morgan2012]       0.94  5574    1
#> shape[Palumbo2014]      1.16  4304    1
#> 
#> Samples were drawn using NUTS(diag_e) at Mon Sep 16 13:27:23 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at 
#> convergence, Rhat=1).
#> 
#> exmpl_> ## End(No test)
#> exmpl_> ## Don't show: 
#> exmpl_> if (requireNamespace("pkgdown", quietly = TRUE) && pkgdown::in_pkgdown()) {
#> exmpl_+   assign("ndmm_net", ndmm_net, .GlobalEnv)
#> exmpl_+   assign("ndmm_fit", ndmm_fit, .GlobalEnv)
#> exmpl_+ }
#> 
#> exmpl_> ## End(Don't show)
#> exmpl_> 
#> exmpl_> 
#> exmpl_> 
# }
# Plot estimated survival curves, and overlay the KM data
# \donttest{
plot(predict(ndmm_fit, type = "survival")) + geom_km(ndmm_net)

# }