Data

Egholm et al report from a Danish Registry on the 30 day cardiac event rates post non-cardiac surgery in patients with prior PCI compared with patients undergoing similar procedures without a history of ischemic heart disease. The cardiac events recorded were: myocardial infarction, cardiac death (and all cause mortality)

Reading Data

Data from the paper was tabulated into a spreadsheet.

library(rethinking)
da<-read.csv("e:/pandoc/data/surgerypci.csv",header=TRUE,sep=",")
names(da)
## [1] "Time"                       "IHD"                       
## [3] "MI.Number"                  "MI.Events"                 
## [5] "Cardiac.Death.Number"       "Cardiac.Death.Events"      
## [7] "All.Cause.Mortality.Number" "All.Cause.Mortality.Events"

Bayesian Models

Myocardial Infarction

m1<-map2stan(
  alist(
    MI.Events~dpois(lambda),
    log(lambda)<-a+ihd[IHD]+time[Time]+log(MI.Number),
    a~dnorm(0,1),
    ihd[IHD]~dnorm(0,1),
    time[Time]~dnorm(0,1)
      ),
  data=da,
  iter=1000
)
## 
## SAMPLING FOR MODEL 'MI.Events ~ dpois(lambda)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 1, Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 1, Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 1, Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 1, Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 1, Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 1, Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 1, Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 1, Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 1, Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 1, Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 1, Iteration: 1000 / 1000 [100%]  (Sampling)
##  Elapsed Time: 0.05 seconds (Warm-up)
##                0.045 seconds (Sampling)
##                0.095 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'MI.Events ~ dpois(lambda)' 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)
## 
## [ 50 / 500 ]
[ 100 / 500 ]
[ 150 / 500 ]
[ 200 / 500 ]
[ 250 / 500 ]
[ 300 / 500 ]
[ 350 / 500 ]
[ 400 / 500 ]
[ 450 / 500 ]
[ 500 / 500 ]
Precis
precis(m1,depth=2)
##          Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a       -2.63   0.63      -3.65      -1.66   135 1.00
## ihd[1]  -2.27   0.57      -3.20      -1.41   111 1.02
## ihd[2]  -0.39   0.56      -1.23       0.54    97 1.02
## time[1]  0.14   0.49      -0.63       0.90    80 1.00
## time[2] -0.81   0.53      -1.68      -0.01    80 1.01
## time[3] -1.96   0.50      -2.76      -1.21    78 1.00
plot(precis(m1,depth=2))

Simulation

sam<-extract.samples(m1)
evr<-exp(sam$a+sam$ihd[,2]+sam$time[,1])
hist(evr,xlab="30 day myocardial infarction rate",main="Myocardial Infarction Rate: Prior PCI < 1 month")

evnoihd<-exp(sam$a+sam$ihd[,1]+sam$time[,1])
evlate<-exp(sam$a+sam$ihd[,1]+sam$time[,3])
or1<-evr/evnoihd
hist(or1,xlab="Odds Ratio: Myocadial Infarction rate: PCI within 1 month: control",main="Myocardial Infarction")

Cardiac Death

m2<-map2stan(
  alist(
    Cardiac.Death.Events~dpois(lambda),
    log(lambda)<-a+ihd[IHD]+time[Time]+log(Cardiac.Death.Number),
    a~dnorm(0,1),
    ihd[IHD]~dnorm(0,1),
    time[Time]~dnorm(0,1)
      ),
  data=da,
  iter=1000
)
## 
## SAMPLING FOR MODEL 'Cardiac.Death.Events ~ dpois(lambda)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 1, Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 1, Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 1, Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 1, Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 1, Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 1, Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 1, Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 1, Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 1, Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 1, Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 1, Iteration: 1000 / 1000 [100%]  (Sampling)
##  Elapsed Time: 0.074 seconds (Warm-up)
##                0.04 seconds (Sampling)
##                0.114 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'Cardiac.Death.Events ~ dpois(lambda)' 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)
## 
## [ 50 / 500 ]
[ 100 / 500 ]
[ 150 / 500 ]
[ 200 / 500 ]
[ 250 / 500 ]
[ 300 / 500 ]
[ 350 / 500 ]
[ 400 / 500 ]
[ 450 / 500 ]
[ 500 / 500 ]

Precis

precis(m2,depth=2)
##          Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a       -2.88   0.59      -3.71      -1.94   137 1.00
## ihd[1]  -2.42   0.56      -3.22      -1.49   132 1.01
## ihd[2]  -0.56   0.56      -1.37       0.35   124 1.03
## time[1]  0.14   0.53      -0.68       0.98   114 1.01
## time[2] -1.23   0.60      -2.21      -0.28   113 1.01
## time[3] -1.75   0.53      -2.61      -0.99   131 1.01
plot(precis(m2,depth=2))

sam2<-extract.samples(m2)
evrcd<-exp(sam2$a+sam2$ihd[,2]+sam2$time[,1])
hist(evrcd,xlab="30 day cardiac death rate",main="Cardiac Death Rate: Prior PCI < 1 month")

evcdnoihd<-exp(sam2$a+sam2$ihd[,1]+sam2$time[,1])
evcdlate<-exp(sam2$a+sam2$ihd[,1]+sam2$time[,3])
or1<-evrcd/evcdnoihd
hist(or1,xlab="Odds Ratio: Cardiac Death rate: PCI within 1 month: control",main="Cardiac Death")

All Cause Mortality

m3<-map2stan(
  alist(
    All.Cause.Mortality.Events~dpois(lambda),
    log(lambda)<-a+ihd[IHD]+time[Time]+log(All.Cause.Mortality.Number),
    a~dnorm(0,1),
    ihd[IHD]~dnorm(0,1),
    time[Time]~dnorm(0,1)
      ),
  data=da,
  iter=1000
)
## 
## SAMPLING FOR MODEL 'All.Cause.Mortality.Events ~ dpois(lambda)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 1, Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 1, Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 1, Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 1, Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 1, Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 1, Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 1, Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 1, Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 1, Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 1, Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 1, Iteration: 1000 / 1000 [100%]  (Sampling)
##  Elapsed Time: 0.086 seconds (Warm-up)
##                0.081 seconds (Sampling)
##                0.167 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'All.Cause.Mortality.Events ~ dpois(lambda)' 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.001 seconds (Sampling)
##                0.001 seconds (Total)
## 
## [ 50 / 500 ]
[ 100 / 500 ]
[ 150 / 500 ]
[ 200 / 500 ]
[ 250 / 500 ]
[ 300 / 500 ]
[ 350 / 500 ]
[ 400 / 500 ]
[ 450 / 500 ]
[ 500 / 500 ]

Precis

precis(m3,depth = 2)
##          Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a       -1.79   0.62      -2.84      -0.92   103 1.03
## ihd[1]  -1.08   0.59      -1.93      -0.12    96 1.00
## ihd[2]  -0.96   0.59      -1.75       0.01    95 1.00
## time[1]  0.23   0.50      -0.53       1.02   129 1.02
## time[2] -0.63   0.50      -1.48       0.12   126 1.02
## time[3] -1.11   0.49      -1.86      -0.30   124 1.02
plot(precis(m3,depth=2))

Cardiac Death Model with Interaction Term

mint<-map2stan(
  alist(
    Cardiac.Death.Events~dpois(lambda),
    log(lambda)<-a+ihd[IHD]+time[Time]+log(Cardiac.Death.Number)+b*time[Time]*ihd[IHD],
    a~dnorm(0,1),
    ihd[IHD]~dnorm(0,1),
    time[Time]~dnorm(0,1),
    b~dnorm(0,1)
      ),
  data=da,
  iter=1000
)
## 
## SAMPLING FOR MODEL 'Cardiac.Death.Events ~ dpois(lambda)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 1, Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 1, Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 1, Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 1, Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 1, Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 1, Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 1, Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 1, Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 1, Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 1, Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 1, Iteration: 1000 / 1000 [100%]  (Sampling)
##  Elapsed Time: 0.077 seconds (Warm-up)
##                0.059 seconds (Sampling)
##                0.136 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'Cardiac.Death.Events ~ dpois(lambda)' 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)
## 
## [ 50 / 500 ]
[ 100 / 500 ]
[ 150 / 500 ]
[ 200 / 500 ]
[ 250 / 500 ]
[ 300 / 500 ]
[ 350 / 500 ]
[ 400 / 500 ]
[ 450 / 500 ]
[ 500 / 500 ]

Precis

precis(mint,depth=2)
##          Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a       -3.58   0.53      -4.48      -2.77   142 1.01
## ihd[1]  -2.24   0.53      -3.09      -1.47   157 1.01
## ihd[2]  -0.10   0.63      -1.19       0.70   152 1.00
## time[1]  0.66   0.69      -0.46       1.71   150 1.00
## time[2] -1.18   0.66      -2.27      -0.19   202 1.00
## time[3] -1.84   0.53      -2.72      -1.08   169 1.00
## b        0.24   0.08       0.12       0.37   253 1.00
plot(precis(mint,depth=2))

smp<-extract.samples(mint)
evcdmnt<-exp(smp$a+smp$ihd[,2]+smp$time[,1]+smp$b*smp$ihd[,2]*smp$time[,1])
evcdmnt2<-exp(smp$a+smp$ihd[,1]+smp$time[,1]+smp$b*smp$ihd[,1]*smp$time[,1])
df1<-data.frame(eventrate=evcdmnt)
df2<-data.frame(eventrate=evcdmnt2)
df1$group<-"PCI<1 month"
df2$group<-"No IHD"

dfc<-rbind(df1,df2)
ggplot(dfc,aes(eventrate,fill=group))+geom_histogram(aes(y=..density..))+ggtitle("30 Day Cardiac Death Rate:\n PCI<1 month versus no IHD ") + 
     theme(plot.title = element_text(lineheight=.8, face="bold"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.