RCT

虛無假設 H0: 效果≥ 臨界值 delta 對立假設 H1: 效果< delta 單尾檢定 p 值 <0.05 單尾 95% CI：上限是平均值+1.645 SE，下限是平均值-1.645 SE 治療意向分析 ITT 及依據計畫書分析 PP 都有效才是可靠的

 Sys.setlocale("LC_CTYPE", "C.UTF-8") library(dplyr) # https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(15)61505-7/fulltext # Set parameters n <- 17187  # total number of patients as per the ISIS-2 trial age_mean <- 60 age_sd <- 10 gender_distribution <- c(0.6, 0.4)  # Assuming 60% male, 40% female treatment_types <- c('Streptokinase', 'Aspirin', 'Both', 'Neither')  # Four treatment groups zodiac_signs <- c('Aries', 'Taurus', 'Gemini', 'Cancer', 'Leo', 'Capricorn', 'Libra',                    'Scorpio', 'Sagittarius', 'Pisces', 'Aquarius', 'Virgo') # Simulate data set.seed(1)  # for reproducibility data <- tibble(   patient_id = 1:n,   age = round(rnorm(n, age_mean, age_sd)),   gender = sample(c('Male', 'Female'), n, replace = TRUE, prob = gender_distribution),   treatment = sample(treatme...

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

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

 https://benjamin-andrew.shinyapps.io/bayesian_trials/ Bayesian analysis of RCT data https://discourse.datamethods.org/t/bayesian-analysis-of-rct-data/1997 ANDROMEDA-SHOCK (or, how to intepret HR 0.76, 95% CI 0.55-1.02, p=0.06) https://discourse.datamethods.org/t/andromeda-shock-or-how-to-intepret-hr-0-76-95-ci-0-55-1-02-p-0-06/1349?fbclid=IwZXh0bgNhZW0BMQABHZl-ZyjEjp0uC3AroV4tHD4uxQ2apk8KRUD0QPmjpoWeHhuYIhwLt1Khog_aem_AWdNbe-zxt7yGCq0VdpfxAgfv4SMMR_sj6uZ_aVPGDevXxQ-vuVziIK1EaG1XpvJo1E The statistical properties of RCTs and a proposal for shrinkage 21 https://onlinelibrary.wiley.com/doi/full/10.1002/sim.9173 David Spiegelhalter’s Gullible Skeptic, and a Bayesian “Hard-Nosed Skeptic” Reanalysis of the ANDROMEDA-SHOCK Trial https://www.bayesianspectacles.org/david-spiegelhalters-gullible-skeptic-and-a-bayesian-hard-nosed-skeptic-reanalysis-of-the-andromeda-shock-trial/?fbclid=IwZXh0bgNhZW0BMQABHYL4q4fmISoKBdOdG63AYYpW7Ctf1r4nDhSdM58APNPMb9_kS3rN1Em15Q_aem_AWcOKXB3novlr9AAbiAJPv2HqEKR...

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...