In bashlee/test: Provides a suite of extension functions for NMA using the `gemtc` package

knitr::opts_chunk$set(echo = FALSE)

library(dplyr)
library(reshape2)
library(ggplot2)
library(gemtc)
library(pander)


# Settings for survival extrapolations (could for example be provided in a global _settings.R file, which is sourced in)
ref.trt <- "B"
ref.std <- "STUDY2"
horizon.months <- c(24, 60, 120)

General information

The NMA was run via the gemtc package.

Input settings

## read in results from Bayesian fit
files_all <- dir(system.file("extdata", "results", package = "gemtcPlus"))
files_rda <- grep("fit-pwe-[0-9].RData", files_all, value = TRUE)
out_all <- vector(mode = "list", length = length(files_rda))
for(i in seq_along(files_rda)){
  load(system.file("extdata",
                   "results",
                   files_rda[i],
                   package = "gemtcPlus"))
  out_all[[i]] <- out
  rm(out)
}

Input data: study overview

knitr::kable(out_all[[1]]$data.arms, caption = "__Table__ Input data (studies and treatment arms)")

Model comparison

knitr::kable(get_pwe_comparison(out_all), caption = "__Table__ Model comparison")

Hazard ratio estimates

# loop through fits
for(i in seq_along(out_all)){
  out <- out_all[[i]]
  title <- paste(out$descr, ", cut point(s) at", paste(unlist(out$model.pars), collapse = ", "), "months", collapse =" ")
  cat("## ", title, "  \n")

  ## Tables: Hazard ratio estimates for each segment
  treatments <- unique(out$data.arms$treatment[order(out$data.arms$treatmentn)])
  # logHR <- get_pwe_contrasts(fit = out$fit, treatments = treatments, ref = ref.trt, 
  #                            cut.pts = out$model.pars$cut.pts, digits = 3,
  #                            exponentiate = FALSE)
  # HR    <- get_pwe_contrasts(fit = out$fit, treatments = treatments, ref = ref.trt, 
  #                            cut.pts = out$model.pars$cut.pts, digits = 3,
  #                            exponentiate = TRUE)
  HR_rev<- get_pwe_contrasts(fit = out$fit, treatments = treatments, ref = ref.trt, 
                             cut.pts = out$model.pars$cut.pts, digits = 3,
                             exponentiate = TRUE, reverse = TRUE)

  # print(knitr::kable(logHR %>% select(-x, -xend), caption = "__Table__ Log hazard ratio estimates of other treatments vs A"))
  # cat("\n\n")
  # 
  # print(knitr::kable(HR %>% select(-x, -xend), caption = "__Table__ Hazard ratio estimates of other treatments vs A"))
  # cat("\n\n")
  print(knitr::kable(HR_rev %>% select(-x, -xend), caption = "__Table__ Hazard ratio estimates of A vs other treatments"))
  cat("\n\n")


  ymax <- 10
  dhr <- HR_rev
  drib <- data.frame(x = as.vector(t(dhr[c("x", "xend")])),           # data structure for ribbons
                     ylo = rep(dhr$lCrI, each = 2),                   #  cap ribbon at ymax
                     yup = rep(dhr$uCrI, each = 2),
                     Comparison = rep(dhr$Comparison, each = 2)) %>%
    mutate(yup = ifelse(yup > ymax, ymax, yup))

  fig <- ggplot(data = dhr) + 
    geom_ribbon(data = drib, aes(x = x, ymin = ylo, ymax = yup), fill = "lightblue", alpha = 0.8) +
    geom_hline(aes(yintercept = 1), col = "darkgrey") +
    geom_segment(aes(x = x, xend = xend, y = Median, yend = Median)) +
    facet_wrap(~Comparison, ncol = 2) +
    scale_y_log10(breaks = c(0.1, 0.5, 1, 2, 10)) +
    coord_cartesian(ylim = c(0.1, ymax)) +
    xlab("Month") + ylab("Hazard ratio") +
    theme_bw()

  cat("__Figure__ Hazard ratio estimates A vs other treatments\n")
  plot(fig)
  cat("\n\n")

  rm(out)
}

Survivor function estimates

The NMA baseline estimate from the r ref.trt arm from r ref.std is used. The contrast estimates from the NMA are then added to obtain the survivor functions for the other interventions.

# loop through fits
for(i in seq_along(out_all)){
  out <- out_all[[i]]
  title <- paste(out$descr, ", cut point(s) at", paste(unlist(out$model.pars), collapse = ", "), "months", collapse =" ")
  cat("## ", title, "  \n")

  ## Plots of survivor functions over time ("NMA result"), ref study/arm and timehorizons specified in settings function
  id.ref.std <- out$data.arms$studyn[which(out$data.arms$study == ref.std & out$data.arms$treatment == ref.trt)]
  id.ref.arm <- out$data.arms$arm[   which(out$data.arms$study == ref.std & out$data.arms$treatment == ref.trt)]
  treatments <- unique(out$data.arms$treatment[order(out$data.arms$treatmentn)])

  for(hm in horizon.months){
    S_extrap <- get_pwe_S(fit = out$fit, cut.pts = out$model.pars$cut.pts, ref.std = id.ref.std, ref.arm = id.ref.arm, treatments = treatments, 
                   time = seq(0, hm, 0.1))

    fig <- ggplot(data = S_extrap) +        
      geom_line(aes(x = time, y = S, col = treatment, linetype = treatment)) +
      ylim(0, 1) +
      xlab("Month") + ylab("Survival probability") +
      theme_bw() +
      theme(legend.title = element_blank())
    cat("__Figure__ Survivor function estimates (time horizon:", hm, "months) \n")
    plot(fig)
    cat("\n\n")

    fig <- ggplot(data = S_extrap) + 
      facet_wrap(~treatment) +
      #geom_line(aes(x = time, y = lCrI), linetype = "dashed") +
      #geom_line(aes(x = time, y = uCrI), linetype = "dashed") +
      geom_ribbon(aes(x = time, ymin = lCrI, ymax = uCrI), fill = "lightblue", alpha = 0.8) +
      geom_line(aes(x = time, y = S)) +
      ylim(0, 1) +
      xlab("Month") + ylab("Survival probability") +
      theme_bw()
    cat("__Figure__ Survivor function estimates by treatment (time horizon:", hm, "months) \n")
    plot(fig)
    cat("\n\n")

    rm(list = c("S_extrap", "fig"))
  }
  rm(out)
}

Model fit: observed KM data vs estimated S(t)

# loop through fits
for(i in seq_along(out_all)){
  out <- out_all[[i]]
  title <- paste(out$descr, ", cut point(s) at", paste(unlist(out$model.pars), collapse = ", "), "months", collapse =" ")
  cat("## ", title, "  \n")

  gof <- get_pwe_GoF(fit = out$fit, 
                     cut.pts = out$model.pars$cut.pts, 
                     data.arms = out$data.arms,
                     data.jg = out$data.jg)

  fig <- ggplot() + 
    geom_line(data = gof %>% filter(type == "nma"), aes(x = time, y = S, col = treatment)) +
    #geom_point(data = gof %>% filter(type == "obs"), aes(x = time, y = S, col = treatment)) + # too many studies (=panels) for good visibility
    geom_line(data = gof %>% filter(type == "obs"), aes(x = time, y = S, col = treatment), linetype = "dashed") +
    facet_wrap(~study, ncol = 2) +
    ylim(0, 1) + xlim(0, 36) +
    xlab("Month") + ylab("Survival probability") +
    theme_bw() +
    theme(legend.position = "top", legend.title = element_blank())

  cat("__Figure__ Goodness-of-fit: estimated (solid lines) and observed (dashed) survivor functions for each study\n")
  plot(fig)
  cat("\n\n")

  rm(list = c("gof", "fig"))
  rm(out)
}

Appendix

# loop through fits
for(i in seq_along(out_all)){
  out <- out_all[[i]]
  title <- paste(out$descr, ", cut point(s) at", paste(unlist(out$model.pars), collapse = ", "), "months", collapse =" ")
  cat("## ", title, "  \n\n")

  jginfo <- get_jags_info(out$fit, include.comments = TRUE)
  cat("```\n", jginfo, "\n```\n\n")

  rm(jginfo)
  rm(out)
}