Subgroup analysis in ISIS-2

 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(treatment_types, n, replace = TRUE, prob = rep(1/4, 4)),

  zodiac_sign = sample(zodiac_signs, n, replace = TRUE)

)

data$treatment <- factor(data$treatment)

data$treatment <- relevel(data$treatment, ref = 'Neither')

# Function to assign death probability based on treatment

assign_death_probability <- function(treatment) {

  if (treatment == 'Streptokinase') {

    return(rbinom(1, 1, 0.097))

  } else if (treatment == 'Aspirin') {

    return(rbinom(1, 1, 0.097))

  } else if (treatment == 'Both') {

    return(rbinom(1, 1, 0.087))

  } else {

    return(rbinom(1, 1, 0.120))

  }

}

# Apply the function to each row

data <- data %>%

  rowwise() %>%

  mutate(vascular_death = assign_death_probability(treatment))

# Mapping of zodiac signs to Traditional Chinese

zodiac_chinese <- c('Aries' = '白羊座', 'Taurus' = '金牛座', 'Gemini' = '雙子座', 

                    'Cancer' = '巨蟹座', 'Leo' = '獅子座', 'Virgo' = '處女座', 

                    'Libra' = '天秤座', 'Scorpio' = '天蠍座', 'Sagittarius' = '射手座', 

                    'Pisces' = '雙魚座', 'Aquarius' = '水瓶座', 'Capricorn' = '摩羯座')

# Convert zodiac signs to Traditional Chinese

data <- data %>%

  mutate(zodiac_sign = zodiac_chinese[zodiac_sign])

head(data)

library(rms)

lrm(vascular_death ~ treatment+ zodiac_sign, data=data)

if(FALSE){

"死亡              係數    標準誤   Wald Z    P

截距          -2.0565 0.0994 -20.70 <0.0001 

Aspirin      -0.2653 0.0694  -3.82 0.0001  

Both          -0.6230 0.0766  -8.13 <0.0001 

Streptokinase -0.298 0.07 -4.26 <0.0001 

天蠍座      -0.0117 0.1302  -0.09 0.9285  

射手座      -0.0261 0.1314  -0.20 0.8428  

巨蟹座      -0.0350 0.1325  -0.26 0.7918  

摩羯座       0.2499 0.1256   1.99 0.0466  

水瓶座       0.1392 0.1285   1.08 0.2787  

獅子座       0.1904 0.1263   1.51 0.1317  

白羊座       0.1529 0.1283   1.19 0.2333  

處女座       0.1063 0.1277   0.83 0.4054  

金牛座       0.0292 0.1291   0.23 0.8212  

雙子座       0.1332 0.1267   1.05 0.2932  

雙魚座       0.0294 0.1313   0.22 0.8231 "

}

留言

張貼留言

這個網誌中的熱門文章

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