Bayesian reanalysis of clinical trials

 # Adapted from code by Dan Lane, pr_in prob of interest

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

# https://discourse.datamethods.org/t/bayesian-analysis-of-rct-data/1997/2?fbclid=IwZXh0bgNhZW0BMQABHVaVaDWAQl1XtzqVr4byQFje4WB3C_SbZC2u7b-rnv8OAJyaK1agsDF3fQ_aem_AWodvwieeBCpXTYGtS2vQ-psFnk549y7FLjaM9SCDRRLJJ6A1T-kRsYS2Hg-p9pEPyY

library(tidyverse)

library(RColorBrewer)

# Estimate Type Labels

est_type <- 1 # 1: HR, 2: OR, 3: RR

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"))

# hr_in HR of interest or MCID, pr_in prob of MCID

pt_est <- 0.5

upper_ci <- 0.9

ci_width <- 0.95

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

hr_a <- 10

hr_b <- 30

se_hr <- NA

se_or <- NA

or_a <- 10

or_b <- 30

trt_n <- 50

ctrl_n <- 50

or_c <- trt_n - or_a

or_d <- ctrl_n - or_b

theta_in <- 1

sd_in <- 0.42

hr_in <- 0.5

pr_in <- 0.05

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

ci_in <- 89

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)

# Colors

like_col <- brewer.pal(3, "Dark2")[1]

post_col <- brewer.pal(3, "Dark2")[2]

prior_col <- brewer.pal(3, "Dark2")[3]

# 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)

  )

留言

張貼留言

這個網誌中的熱門文章

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...
閱讀完整內容