Data from EXCEL and NOBLE Trials

This is a superficial exploration of the plausibility of the percutaneous coronary intervention being associated with clinically meaningful higher rate of composite major cardiovascular complications compared to coronary artery bypass surgery for left main disease. The data used is from the published EXCEL and NOBLE trials.

Generalized Linear Model of Common Event Rate

nep<-ceiling(3*121/5)
nec<-ceiling(3*81/5)
daten<-data.frame(group=c(1,2,1,2),number=c(as.integer(948),as.integer(957),as.integer(592),as.integer((592))),events=c(208,174,nep,nec))
model0<-map2stan(
  alist(
    events~dbinom(number,p),
   logit(p)<-a,
    a~dnorm(0,5)
    ),
  data=daten,
  iter=1e4
)
## 
## SAMPLING FOR MODEL 'events ~ dbinom(number, p)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 1, Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 1, Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 1, Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 1, Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 1, Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 1, Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 1, Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 1, Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 1, Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 1, Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 1, Iteration: 10000 / 10000 [100%]  (Sampling)
##  Elapsed Time: 0.078 seconds (Warm-up)
##                0.073 seconds (Sampling)
##                0.151 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'events ~ dbinom(number, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 0 seconds (Warm-up)
##                0 seconds (Sampling)
##                0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 500 / 5000 ]
[ 1000 / 5000 ]
[ 1500 / 5000 ]
[ 2000 / 5000 ]
[ 2500 / 5000 ]
[ 3000 / 5000 ]
[ 3500 / 5000 ]
[ 4000 / 5000 ]
[ 4500 / 5000 ]
[ 5000 / 5000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
precis(model0,depth=2)
##    Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a -1.64   0.05      -1.72      -1.56  3663    1
plot(precis(model0,depth=2))

Posterior Probability of observed difference given common event rate
s1<-sim(model0,n=5000)
## [ 500 / 5000 ]
[ 1000 / 5000 ]
[ 1500 / 5000 ]
[ 2000 / 5000 ]
[ 2500 / 5000 ]
[ 3000 / 5000 ]
[ 3500 / 5000 ]
[ 4000 / 5000 ]
[ 4500 / 5000 ]
[ 5000 / 5000 ]
es1<-(s1[,1]-s1[,2])/950
es2<-(s1[,3]-s1[,4])/592
oe1<-daten$events[1]/daten$number[1]-daten$events[2]/daten$number[2]
po1<-length(es1[es1>oe1])/5000
hist(es1,main="Effect Size (EXCEL)",xlab="Effect Size")
abline(v=oe1,col="red",lwd=2)
text(-0.05,900,paste("P(Event Rate>Observed): ",round(po1,3)))

oe2<-daten$events[3]/daten$number[3]-daten$events[4]/daten$number[4]
po2<-length(es2[es2>oe2])/5000
hist(es2,main="Effect Size (NOBLE)",xlab="Effect Size")
abline(v=oe2,col="red",lwd=2)
text(-0.05,900,paste("P(Event Rate>Observed): ",round(po2,3)))

Generalized Linear Model: Intervention Predictor

nep<-ceiling(3*121/5)
nec<-ceiling(3*81/5)
daten<-data.frame(group=c(1,2,1,2),number=c(as.integer(948),as.integer(957),as.integer(592),as.integer((592))),events=c(208,174,nep,nec))
model1<-map2stan(
  alist(
    events~dbinom(950,p),
   logit(p)<-a_group[group],
    a_group[group]~dnorm(0,5)
    ),
  data=daten,
  iter=1e4
)
## 
## SAMPLING FOR MODEL 'events ~ dbinom(950, p)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 1, Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 1, Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 1, Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 1, Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 1, Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 1, Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 1, Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 1, Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 1, Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 1, Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 1, Iteration: 10000 / 10000 [100%]  (Sampling)
##  Elapsed Time: 0.11 seconds (Warm-up)
##                0.144 seconds (Sampling)
##                0.254 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'events ~ dbinom(950, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 0 seconds (Warm-up)
##                0 seconds (Sampling)
##                0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 500 / 5000 ]
[ 1000 / 5000 ]
[ 1500 / 5000 ]
[ 2000 / 5000 ]
[ 2500 / 5000 ]
[ 3000 / 5000 ]
[ 3500 / 5000 ]
[ 4000 / 5000 ]
[ 4500 / 5000 ]
[ 5000 / 5000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
precis(model1,depth=2)
##             Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_group[1] -1.75   0.06      -1.85      -1.64  3427    1
## a_group[2] -2.02   0.07      -2.14      -1.91  3218    1
plot(precis(model1,depth=2))

Composite End-point: PCI versus CABG

sane<-extract.samples(model1,n=5000)
pci<-logistic(sane$a_group[,1])
cabg<-logistic(sane$a_group[,2])
diff<-pci-cabg
hp<-HPDI(diff,prob=0.95)


hist(diff,main="Posterior Distribution",xlab="Event rate: PCI-CABG",ylim=c(0,1100))
text(0.01,1000,paste("mean (95%HPDI): ",round(mean(diff),3),"(",round(hp[1],3),",",round(hp[2],3),")"))

Odds Ratio

oddsr<-exp(sane$a_group[,1])/exp(sane$a_group[,2])
ohp<-HPDI(oddsr)
hist(oddsr,main="Posterior Distribution Odds Ratio for Composite End-point",
     xlab="Odds Ratio: PCI/CABG",
     ylim=c(0,1500))
text(1.5, 1000,paste("odds ratio (95% HPDI): ",round(mean(ohp),3),"(",round(ohp[1],3) ,",",round(ohp[2],3),")"))