RCT

# https://jamanetwork.com/journals/jama/fullarticle/2819082 Use pymc3 pystan for a bayesian reanalysis with weakly informative cauchy gaussian prior mean and half-cauchy gaussian sd with a smaller scale, 95% hdi, posterior probability of the effect being more than 0, and a sensitivity analysis with different scales, and run it without showing the code: Of 447 patients enrolled, 227 were randomized to acetaminophen and 220 to placebo. Days alive and free of organ support to day 28 were not meaningfully different for acetaminophen (20.2 days; 95% CI, 18.8 to 21.6) vs placebo (19.6 days; 95% CI, 18.2 to 21.0; P = .56; difference, 0.6; 95% CI, −1.4 to 2.6). OR 1.78 (95% CI, 0.98-3.22). Use pystan for a bayesian reanalysis of rct with weakly informative prior mean and sd, 95% hdi, posterior probability of the effect being more than 1, and a sensitivity analysis with different scales, and run it without showing the code

虛無假設 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...

import pymc3 as pm import numpy as np import arviz as az # Empagliflozin after Acute Myocardial Infarction 24 https://www.nejm.org/doi/full/10.1056/NEJMoa2314051 log_or = np.log(0.9) log_or_lower = np.log(0.76) log_or_upper = np.log(1.06) # The scale for the Cauchy distribution cauchy_scale = 0.707 # Model setup with pm.Model() as model:   # Prior distribution for the log-OR     log_or_prior = pm.Cauchy("log_OR", alpha=0, beta=cauchy_scale)     # Likelihood (normal approximation of the distribution of log-OR)     se_log_or = (log_or_upper - log_or_lower) / 3.92     likelihood = pm.Normal("likelihood", mu=log_or_prior, sigma=se_log_or, observed=np.array([log_or]))     # Sample from the posterior distribution     trace = pm.sample(500000, cores=1, tune=150000, chains=4) # Convert log-OR samples to HR samples hr_samples = np.exp(trace["log_OR"]) # Calculate the posterior probability that HR < 1 post_prob_hr_less_than_1 = n...

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

 # The ANDROMEDA-SHOCK Randomized Clinical Trial 19 https://jamanetwork.com/journals/jama/fullarticle/2724361 #Bayesian reanalysis https://www.atsjournals.org/doi/10.1164/rccm.201905-0968OC 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: Among 424 patients randomized (mean age, 63 years; 226 [53%] women), 416 (98%) completed the trial. By day 28, 74 patients (34.9%) in the peripheral perfusion group and 92 patients (43.4%) in the lactate group had died (hazard ratio, 0.75 [95% CI, 0.55 to 1.02]; P = .06; risk difference, −8.5% [95% CI, −18.2% to 1.2%]).

雙盲、隨機化、意向治療分析、沒有交叉或外部污染、沒有選擇偏差、沒有退出偏差 隨機分配： 表一謬誤：不要比較各組的基礎值，因為隨機分配保證各組干擾因子（無論已知、未知，有測量、未測量）的期望值一定是平衡的。 假設：可忽略性（沒有未觀察到的干擾因子，亦即在觀察到的共變數之下，治療跟潛在結局沒有相關）、Y(0), Y(1)  ⫫  T|X 1. SUTVA（潛在結局的穩定單位處理值假設）：讓我們能估計平均治療效果 ATE 和治療組的平均治療效果 ATT。集群 RCT 能減少單位之間的干擾、溢出效應、背景效應（也能用分層 RCT 減少）。 2. 排除限制：某些人不受治療影響。工具變項則是除了經由內在變項 x 以外不能影響 y 3. 非零性：每個可能的分組都不能是零 4. 一致性 5. 單位之間無干擾 6. 治療是穩定的 7. 沒有測量誤差 8. 沒有退出 9. 沒有不遵從性 10. 沒有設限資料 • 內在效度：消除選擇偏差和干擾因子（造成治療和結局的因子）。 • 外在效度：隨機分配保證因果關係是可移動的 transportable （各組是可交換的、治療效果可以應用在具有相同疾病的不同族群）。由於嚴格的收納/排除條件，RCT 不能像隨機抽樣（調查研究）一樣能保證效果可以推廣 generalizable 到整個目標族群。 可移動性的條件：可識別性（數據可供識別因果圖 DAG）、強交換性（各組的聯合分布在重新排列 permutation 之後維持不變，源域具有 DAG 全部的相關變數，源域和目標域的數據分布有重疊/相似性、因果關係 DAG 不變）、相同的干預 可交換性：弱交換性（重新排列後的聯合分布不變，但是邊際分布改變）、強交換性（重新排列後的聯合分布和邊際分布都不變）、部分交換性（只有子集合具有交換性） 潛在結局、反事實 事先決定要校正哪些跟結局相關的因子（ANCOVA）能增加精確性。 多重比較：FWER, FDR 統計檢定力 治療效果的異質性（HTE）、條件平均治療效果（CATE）：大數據能用機器學習分析，比交互作用、次群分析更好。 驚訝 S 值：S = −log2(p), p=(1/2)ˢ —————————————————————- p 值 0.5 0.2  0.1  0.05  0.01  0.005 ...

用次群分層設計 事先決定要分析的少數次群 人數至少需要主要終點的 4-16 倍 校正多重比較 校正事先決定的共變數 治療效果異質性分析 混合（多層次）模型 貝氏分析

利用觀察數據來模仿隨機對照試驗(RCT)結果 優點: 1. 成本效益高,因為利用現有數據,無需招募參與者或進行新的試驗。  2. 可以比傳統 RCT 更快得出結果。 3. 樣本量更大,可提高結果的精確度。 4. 估計的治療效果更具普遍性，因為不像 RCT 一樣具有狹窄的收納/排除條件。 5. 驗證或補充RCT結果,提供進一步的治療有效性證據。 6. 減少：零時（前導時間/不死時間）偏差、錯誤分類偏差、選擇偏差（收納/排除條件）、適應症干擾、干擾因子。 缺點: 1. 對數據和治療效果的強假設可能並不成立。 2. 選擇偏差。   3. 無法控制干擾變數,導致治療效果估計偏差。 4. 不涉及隨機化,導致治療效果估計偏差。 5. 高度依賴統計模型。 6.假設和方法不完全透明或可重現。  7. 沒有標準化的模擬方法,不同研究結果難以比較。 8. 監管機構可能不接受模擬結果作為治療有效性的證據。

樣本不是隨機的： • 表ㄧ謬誤：各組的標準差沒有意義。表二謬誤是觀察性研究的。 • 只能估計樣本平均治療效果 SATE，不能估計族群平均治療效果 PATE。 分配是隨機的： • 各組的基礎期望值是一樣的，亦即沒有偏差（系統誤差）、各組的干擾因子（有測量、未測量）是平衡的。 • 各組的基礎差別是隨機的，不會造成估計治療效果的偏差，雖然用事先決定會影響結局的共變數校正來估計條件平均治療效果 CATE 能增加效率。 • 必須用治療意向分析 ITT，不能用符合計畫書分析。 次群分析： • 需要的樣本數至少是 4-16 倍，亦即需要大數據。 • 次群事先決定。 • 次群的數目愈少愈好。 • 用治療-次群交互作用分析，而非各別分析次群：各次群之間的干擾因子是不平衡的。 • 需要校正多重比較，例如用 FDR。 • 用機器學習分析 CATE 或治療效果的異質性 HTE 更好：連續變項、非線性/多重交互作用，但是需要大數據。 統計檢定力： 由於人力、物力、財力、時間的限制，大部分都會高估預期的效果量與事件率，也會低估標準差。實務上的統計檢定力大部分都是遠低於 0.8 的，因此容易有假陽性、假陰性的結果，同時會高估效果量，亦即贏家的詛咒，因此需要用收縮估計來校正。 p 值的意義： p <0.05: 用收縮估計校正治療效果的「贏家的詛咒」。 p ≥0.05: 用 pymc3 估計治療有效的機率。