Bayesian reanalysis of clinical trials

 # https://benjamin-andrew.shinyapps.io/andromeda_shock_bayesian/?fbclid=IwAR1qHXH4TdjT6h1APC0IS07tXcmRO-36UKvSMb_h8NRzH82TLw8zXXYqToo&mibextid=Zxz2cZ

plot_posterior <- function(est_type = 1, pt_est = 0.5, upper_ci = 0.9, ci_width = 0.95, hr_a = 10, hr_b = 30, or_a = 10, or_b = 30, trt_n = 50, ctrl_n = 50, theta_in = 1, sd_in = 0.42, hr_in = 0.5, pr_in = 0.05, ci_in = 89) {

library(tidyverse)

library(RColorBrewer)

  short_lab <- ifelse(est_type == 1, "HR", ifelse(est_type == 2, "OR", "RR"))

  long_lab <- ifelse(est_type == 1, "Hazard Ratio", ifelse(est_type == 2, "Odds Ratio", "Risk Ratio"))

  likelihood_p <- (1 - ((1 - ci_width) / 2))

  or_c <- trt_n - or_a

  or_d <- ctrl_n - or_b

  prior_theta <- log(theta_in)

  prior_sd <- sd_in

# Calculate Likelihood Parameters

  likelihood_theta <- ifelse(est_type == 1, log(pt_est),

                             ifelse(est_type == 2, log(((or_a + 0.5) * (or_d + 0.5)) / ((or_b + 0.5) * (or_c + 0.5))),

                                    log((or_a / (or_a + or_c)) / (or_b / (or_b + or_d)))))

  likelihood_sd <- ifelse(est_type == 1, sqrt((1 / hr_a) + (1 / hr_b)),

                          ifelse(est_type == 2, sqrt((1 / (or_a + 0.5)) + (1 / (or_b + 0.5)) + (1 / (or_c + 0.5)) + (1 / (or_d + 0.5))),

                                 sqrt((or_c / (or_a * (or_a + or_c))) + (or_d / (or_b * (or_b + or_d))))))

# Calculate Posterior Parameters

  post_theta <- ((prior_theta / (prior_sd)^2) + (likelihood_theta / likelihood_sd^2)) / 

    ((1 / (prior_sd)^2) + (1 / likelihood_sd^2))

  post_sd <- sqrt(1 / ((1 / (prior_sd)^2) + (1 / likelihood_sd^2)))

# Plot data

  x <- seq(-5, 3, by = 0.01)

  prior_plot <- dnorm(x, prior_theta, prior_sd)

  likelihood_plot <- dnorm(x, likelihood_theta, likelihood_sd)

  posterior_plot <- dnorm(x, post_theta, post_sd)

  plot_data <- tibble(

    x = rep(x, 3)

  ) %>%

    mutate(

      dist = rep(c("prior", "likelihood", "posterior"), each = nrow(.) / 3),

      y = c(prior_plot, likelihood_plot, posterior_plot),

      x = exp(x),

      y = exp(y)

    )

# Max xlim

  max_x <- max(likelihood_theta + (exp(likelihood_theta) * 1.5), 2)

# Credible interval

  lower_cred <- round(exp(qnorm((1 - (ci_in/100)) / 2, post_theta, post_sd)), 2)

  upper_cred <- round(exp(qnorm(1 - (1 - (ci_in/100)) / 2, post_theta, post_sd)), 2)

  mid_cred <- round(exp(qnorm(0.5, post_theta, post_sd)), 2)

# Plot

  ggplot(plot_data, aes(x = x, y = y, group = dist)) + 

    geom_vline(xintercept = 1, linetype = "dashed",

               color = "grey50", alpha = 0.75, size = 0.75) + 

    geom_line(aes(color = dist), size = 1.1) + 

    scale_color_brewer(name = NULL, type = "qual", palette = "Dark2",

                       breaks = c("likelihood", "prior", "posterior"),

                       labels = c("Likelihood", "Prior", "Posterior")) + 

    xlim(0, max_x) + 

    labs(

      x = long_lab,

      y = "Probability Density"

    ) + 

    annotate(geom = "text",

             label = paste("Posterior probability ", short_lab,

                           paste(" < 1: ", 

                                 round(pnorm(log(1), post_theta, post_sd, 

                                             lower.tail = TRUE), 3), sep = ""), sep = ""),

             x = max_x, y = max(plot_data$y), hjust = 1,

             fontface = "bold") + 

    annotate(geom = "text",

             label = paste("Posterior median (", ci_in,

                           paste("% credible interval): ",

                                 mid_cred,

                                 paste(" (", lower_cred, sep = ""),

                           paste(", ", upper_cred, sep = ""),

                           paste(")", sep = ""), sep = ""), sep = ""),

             x = max_x, y = max(plot_data$y) - (2 * max(plot_data$y)/25), hjust = 1,

             fontface = "bold") + 

    theme_classic() + 

    theme(

      legend.position = "bottom",

      text = element_text(family = "Gill Sans MT"),

      axis.ticks.y = element_blank(),

      axis.text.y = element_blank(),

      axis.title = element_text(size = 15),

      axis.text = element_text(size = 12),

      legend.text = element_text(size = 15)

    )

}


plot_posterior(est_type = 1, pt_est = 0.5, upper_ci = 0.9, ci_width = 0.95,

               hr_a = 10, hr_b = 30, or_a = 10, or_b = 30, trt_n = 50, ctrl_n = 50,

               theta_in = 1, sd_in = 0.42, hr_in = 0.5, pr_in = 0.05, ci_in = 89)

留言

張貼留言

這個網誌中的熱門文章

Bayesian reanalysis of ORBITA

import pymc3 as pm import numpy as np import arviz as az # Percutaneous coronary intervention in stable angina (ORBITA): a double-blind, randomised controlled trial 18 https://www.sciencedirect.com/science/article/abs/pii/S0140673617327149 diff = 16.6 diff_lower = -8.9 diff_upper = 42 # The scale for the Cauchy distribution cauchy_scale = 5 # Model setup with pm.Model() as model:     # Prior distribution for the difference     diff_prior = pm.Cauchy("diff", alpha=0, beta=cauchy_scale)      # Likelihood (normal approximation of the distribution of the difference)     se_diff = (diff_upper - diff_lower) / 3.92     likelihood = pm.Normal("likelihood", mu=diff_prior, sigma=se_diff, observed=np.array([diff]))     # Sample from the posterior distribution     trace = pm.sample(500000, cores=1, tune=200000, chains=4) # Calculate the posterior probability that the difference is negative post_prob_diff_negative = np.mean(t...
閱讀完整內容

Bayesian reanalysis with pymc3

 # Empagliflozin after Acute Myocardial Infarction 24 https://www.nejm.org/doi/full/10.1056/NEJMoa2314051 Do a Bayesian reanalysis of a RCT for log(hazard ratio) with 500000 samplings with a weakly informative prior, calculate the probability of hazard ratio less than 1 and a sensitivity analysis using pymc3 and only show the results in traditional Chinese instead of the code: A total of 3260 patients were assigned to receive empagliflozin and 3262 to receive placebo. During a median follow-up of 17.9 months, a first hospitalization for heart failure or death from any cause occurred in 267 patients (8.2%) in the empagliflozin group and in 298 patients (9.1%) in the placebo group, with incidence rates of 5.9 and 6.6 events, respectively, per 100 patient-years (hazard ratio, 0.90; 95% confidence interval [CI], 0.76 to 1.06; P=0.21). 使用 pymc3 對一項 RCT 進行貝氏再分析，在弱信息先驗條件下抽樣 500000 次，計算危害比小於 1 的概率，並進行敏感性分析，只用繁體中文顯示結果，而不顯示代碼:  Do a Bayesian reanalysis of a RCT with 500000 samplings wit...
閱讀完整內容

Bayesian reanalysis of Apolipoprotein A1 Infusions

Apolipoprotein A1 Infusions and Cardiovascular Outcomes after Acute Myocardial Infarction 24 https://www.nejm.org/doi/full/10.1056/NEJMoa2400969 P(θ | D) = (P(D | θ) * P(θ)) / P(D) marginal likelihood P(D) = ∫ P(D | θ) * P(θ) dθ Likelihood: P(D | θ) ~ N(μ_l, σ_l^2) Likelihood ratio: (P(D | θ1))/(P(D | θ2)*) Bayes factor: (P(D | θ1)* P(θ1))/(P(D | θ2)* P(θ2)) Prior: P(θ) ~ N(μ_p, σ_p^2) P(θ | D) ∝ P(D | θ) * P(θ) When we normalize the product of the two PDFs by dividing by the marginal likelihood p(D), the constant term cancels out, and we are left with a normal distribution with mean μ_new and variance σ_new^2. Prior mean 0, sd 0.5 log(HR) = log(0.93) = -0.0726 SE = (log(1.05) - log(0.81)) / (2 * 1.96) = 0.0644 Prior precision = 1 / (0.35^2) = 8.16 Likelihood precision = 1 / (0.0644^2) = 241.18 Posterior precision = 8.16 + 241.18 = 249.34 Posterior mean = (8.16 * 0 + 241.18 * -0.0726) / 249.34 = -0.0701 Posterior SD = sqrt(1 / 249.34) = 0.0634 Posterior 95% credible interval for the lo...
閱讀完整內容