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(trace["diff"] < 0)

# Print the posterior probability

print("Posterior probability that the difference is negative:", post_prob_diff_negative)

# Calculate the 95% HDI

hdi_95 = az.hdi(trace, hdi_prob=0.95)

# Print the 95% HDI

print("95% HDI:")

print("Lower bound:", hdi_95["diff"][0])

print("Upper bound:", hdi_95["diff"][1])

# Calculate convergence diagnostics

summary = az.summary(trace, kind='diagnostics')

print(summary)

az.plot_trace(trace)