In giulianonetto/bayesdca: Bayesian Decision Curve Analysis

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dpi = 300, fig.width = 7
)

library(tidyverse)
library(bayesDCA)
library(rms)
theme_set(theme_bw())
data("Docking2021Data")
d <- Docking2021Data; rm(Docking2021Data)

aml_pmp <- d %>% filter(fig_id == "Fig3B")
d$aps_lasso_class2 <- ifelse(
  d$aps_lasso_score > median(aml_pmp$aps_lasso_score),
  "High", "Low"
) %>% factor(levels = c("Low", "High"))
dd <- datadist(aml_pmp)
options(datadist = "dd")
fit <- cph(Surv(survival_weeks, survival_event) ~ aps_lasso_score,
           data = aml_pmp, x = TRUE, y = TRUE, surv = TRUE)

plot_preds <- function(.t) {
  .pred <- Predict(fit, aps_lasso_score, fun = \(i) 1-i, time = .t)
  .xint <- quantile(
    aml_pmp$aps_lasso_score,
    c(1/2)
  )
  .yint <- Predict(fit, aps_lasso_score = .xint, fun = \(i) 1-i, time = .t)
  .yint_pretty <- paste0(round(.yint$yhat*100), "%")
  years <- c("one", "two", "three", "four", "five", "six")[.t/52]
  years <- paste0(years, ifelse(years == "one", " year", " years"))
  ggplot(.pred) +
    geom_rug(
      data = aml_pmp %>% filter(survival_weeks <= 52 & survival_event == 1),
      aes(x = aps_lasso_score,
          color = aps_lasso_class),
      inherit.aes = F, sides = 't'
    ) +
    geom_rug(
      data = aml_pmp %>% filter(survival_weeks > 52 | survival_event == 1),
      aes(x = aps_lasso_score,
          color = aps_lasso_class),
      inherit.aes = F,sides = 'b'
    ) +
    geom_vline(
      xintercept = .xint,
      linetype = "longdash",
      alpha = .3
    ) +
    geom_hline(
      yintercept = .yint$yhat,
      linetype = "longdash",
      alpha = .3
    ) +
    annotate(
      "text", x = -.25, y = .85, 
      label = str_glue(
        '"Low APS" group has up to\n{.yint_pretty} risk of death within {years}'
      )
    ) +
    geom_curve(
      aes(
        y = 0.7, yend = 0.7,
        x = min(aml_pmp$aps_lasso_score),
        xend = .xint
      ),
      curvature = -.1,
      color = "blue"
    ) +
    labs(
      y = str_glue("{.t} week death probability")
    ) +
    scale_y_continuous(
      breaks = scales::pretty_breaks(7),
      labels = scales::percent
    ) +
    scale_color_manual(
      values = list(
        "Low" = "blue",
        "High" = "red"
      )
    )
}
plot_preds(52)

external validation

LAML cohort

laml <- d %>% filter(fig_id == "Fig3C")
dd2 <- datadist(laml)
options(datadist = "dd2")
s <- Surv(laml$survival_weeks, laml$survival_event)
v <- val.surv(fit, newdata = laml, u = 52,
              S = s)
plot(v)
estimates <- survest(fit, newdata = laml, times = 52)
plot(laml$aps_lasso_score, 1 - estimates$surv,
     col = c("black", "red")[
       as.numeric(laml$survival_event==1) + 1
     ])
rcorr.cens(x = estimates$surv, S = s)

fit2 <- cph(Surv(survival_weeks, survival_event) ~ aps_lasso_class2,
           data = laml, x = TRUE, y = TRUE, surv = TRUE)
s <- Surv(laml$survival_weeks, laml$survival_event)
v <- val.surv(fit2, newdata = laml, u = 52,
              S = s)
plot(v)
estimates <- survest(fit2, newdata = laml, times = 52)
plot(laml$aps_lasso_score, 1 - estimates$surv,
     col = c("black", "red")[
       as.numeric(laml$survival_event==1) + 1
     ])
rcorr.cens(x = estimates$surv, S = s)
rms::Predict(fit2, aps_lasso_class2, time = 52)

fit2 <- cph(Surv(survival_weeks, survival_event) ~ aps_lasso_class,
           data = laml, x = TRUE, y = TRUE, surv = TRUE)
s <- Surv(laml$survival_weeks, laml$survival_event)
v <- val.surv(fit2, newdata = laml, u = 52,
              S = s)
plot(v)
estimates <- survest(fit2, newdata = laml, times = 52)
plot(laml$aps_lasso_score, 1 - estimates$surv,
     col = c("black", "red")[
       as.numeric(laml$survival_event==1) + 1
     ])
rcorr.cens(x = estimates$surv, S = s)
rms::Predict(fit2, aps_lasso_class, time = 52)