Time to Reperfusion

Preamble

The following is a first exploration of

r markdown
the rethinking package written by Professor Richard McElreath

The data relates to time to reperfusion for ST elevation myocardial infarction.

Any visitor is urged to look at Professor McElreath’s site. The youtube videos are excellent. I have just received a copy of his book and cannot wait to dive in.

Data

Prospective consecutive STEMI undergoing percutaneous coronary intervention were collected. In this analysis time to reperfusion (T) was defined as time from first medical contact to time of first device use. The effect of time of day and day of week is to be assessed.

setwd("c:/users/mark/documents/work/stemibayes/")
stemi<-read.csv("atcm.csv",header=TRUE)

Bayesian Models (based on McElreath Rethinking)

library(rethinking)

Model 1

model1<-map2stan(
  alist(T~dlnorm(mu,sigma),
        mu<-a+day[DOW]+ b*TOD,
        a~dnorm(0,1),
        day[DOW]~dnorm(0,1),
        b~dnorm(0,1),
        sigma~dcauchy(0,1)),
  data=stemi)

## 
## SAMPLING FOR MODEL 'T ~ dlnorm(mu, sigma)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1, Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1, Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1, Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1, Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%]  (Sampling)
##  Elapsed Time: 0.254 seconds (Warm-up)
##                0.131 seconds (Sampling)
##                0.385 seconds (Total)
## 
## [1] "The following numerical problems occured the indicated number of times after warmup on chain 1"
##                                                                                         count
## Exception thrown at line 23: lognormal_log: Scale parameter is 0, but must be > 0!          2
## Exception thrown at line 23: lognormal_log: Scale parameter is inf, but must be finite!     2
## [1] "When a numerical problem occurs, the Hamiltonian proposal gets rejected."
## [1] "See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected"
## [1] "If the number in the 'count' column is small,  do not ask about this message on stan-users."
## 
## SAMPLING FOR MODEL 'T ~ dlnorm(mu, sigma)' 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)
## 
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

precis(model1,depth=2)

##         Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a       4.50   0.39       3.86       5.09   232 1.00
## day[1]  0.61   0.38       0.01       1.22   245 1.00
## day[2]  0.51   0.43      -0.16       1.18   262 1.00
## day[3]  0.28   0.40      -0.34       0.91   292 1.00
## day[4]  0.31   0.42      -0.41       0.96   318 1.00
## day[5]  0.77   0.39       0.18       1.42   233 1.01
## day[6]  0.58   0.47      -0.12       1.36   417 1.00
## day[7]  1.25   0.50       0.55       2.08   475 1.00
## b      -0.02   0.02      -0.05       0.01   615 1.00
## sigma   0.86   0.08       0.73       0.98   717 1.01

Diagnostics

plot(model1)
dashboard(model1)

Trace plots

Dashboard

Model 2

stemi$todsq<-stemi$TOD*stemi$TOD
model2<-map2stan(alist(
  T~dlnorm(mu,sigma),
  mu<-a+day[DOW]+ b1*TOD+b2*todsq,
  a~dnorm(0,1),
  day[DOW]~dnorm(0,1),
  b1~dnorm(0,1),
  b2~dnorm(0,1),
  sigma~dcauchy(0,1)),
  data=stemi)

## 
## SAMPLING FOR MODEL 'T ~ dlnorm(mu, sigma)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1, Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1, Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1, Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1, Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%]  (Sampling)
##  Elapsed Time: 1.041 seconds (Warm-up)
##                0.432 seconds (Sampling)
##                1.473 seconds (Total)
## 
## [1] "The following numerical problems occured the indicated number of times after warmup on chain 1"
##                                                                                         count
## Exception thrown at line 26: lognormal_log: Scale parameter is inf, but must be finite!    15
## [1] "When a numerical problem occurs, the Hamiltonian proposal gets rejected."
## [1] "See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected"
## [1] "If the number in the 'count' column is small,  do not ask about this message on stan-users."
## 
## SAMPLING FOR MODEL 'T ~ dlnorm(mu, sigma)' 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)
## 
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

precis(model2, depth = 2)

##         Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a       4.98   0.47       4.18       5.67   226    1
## day[1]  0.80   0.44       0.09       1.45   300    1
## day[2]  0.59   0.45      -0.10       1.35   305    1
## day[3]  0.48   0.45      -0.08       1.29   328    1
## day[4]  0.52   0.47      -0.20       1.30   346    1
## day[5]  0.85   0.42       0.19       1.56   293    1
## day[6]  0.60   0.52      -0.19       1.48   372    1
## day[7]  1.20   0.52       0.36       2.00   378    1
## b1     -0.17   0.06      -0.27      -0.08   227    1
## b2      0.01   0.00       0.00       0.01   249    1
## sigma   0.79   0.08       0.66       0.90   383    1

Diagnostics

plot(model2)
dashboard(model2)

Trace plots

Dashboard

Predictive Posterior

td<-seq(0,24,0.5)
new<-data.frame(DOW=1,TOD=td,todsq=td^2)
pred<-exp(link(model2,data=new,n=1000))

## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

pred.mean<-apply(pred,2,mean)
pred.HPDI<-apply(pred,2,HPDI)
plot(T[DOW==1]~TOD[DOW==1],data=stemi,
     xlab="Time of Day",ylab="Time to reperfusion (min)",
     main="Monday")
lines(td,pred.mean)
shade(pred.HPDI,td)

Function

dowfunm2<-function(x){
  ifelse(x==1,lab<-"Monday",
         ifelse(x==2,lab<-"Tuesday",
                ifelse(x==3,lab<-"Wednesday",
                       ifelse(x==4,lab<-"Thursday",
                              ifelse(x==5,lab<-"Friday",
                                     ifelse(x==6,lab<-"Saturday",lab<-"Sunday"
                                     ))))))
         tdn<-seq(0,24,0.5)
         dafr<-data.frame(DOW=x,TOD=tdn,todsq=tdn^2)
         lnk<-exp(link(model2,data=dafr,n=1000))
         lnk.mean<-apply(lnk,2,mean)
         lnk.HPDI<-apply(lnk,2,HPDI)
         par(mar=c(2,2,2,2))
         plot.new()
      plot(T[DOW==x]~TOD[DOW==x],data=stemi,xlim=c(0,24),ylim=c(0,600),
           xlab="Time of Day",
           ylab="Time to reperfusion (min)",
           main=lab)
      lines(tdn,lnk.mean)
      shade(lnk.HPDI,tdn)
    
  
         }

Model 3

model3<-map2stan(
  alist(
    T~dlnorm(mu, sigma),
    mu<-a+day[DOW]+gamma*TOD,
    gamma<-b1+b2*day[DOW],
    day[DOW]~dnorm(0,1),
    a~dnorm(0,1),
    b1~dnorm(0,1),
    b2~dnorm(0,1),
    sigma~dcauchy(0,1)
  ),
  data=stemi
  
)

## 
## SAMPLING FOR MODEL 'T ~ dlnorm(mu, sigma)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1, Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1, Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1, Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1, Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%]  (Sampling)
##  Elapsed Time: 0.516 seconds (Warm-up)
##                0.26 seconds (Sampling)
##                0.776 seconds (Total)
## 
## [1] "The following numerical problems occured the indicated number of times after warmup on chain 1"
##                                                                                         count
## Exception thrown at line 29: lognormal_log: Scale parameter is inf, but must be finite!    11
## [1] "When a numerical problem occurs, the Hamiltonian proposal gets rejected."
## [1] "See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected"
## [1] "If the number in the 'count' column is small,  do not ask about this message on stan-users."
## 
## SAMPLING FOR MODEL 'T ~ dlnorm(mu, sigma)' 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)
## 
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

precis(model3,depth=2)

##         Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## day[1]  0.49   0.53      -0.33       1.33   511 1.00
## day[2]  0.66   0.53      -0.20       1.49   545 1.00
## day[3] -0.06   0.62      -1.01       0.99   573 1.00
## day[4]  0.15   0.62      -0.77       1.22   766 1.00
## day[5]  1.65   0.49       0.98       2.51   477 1.00
## day[6] -0.03   0.77      -1.23       1.18   511 1.02
## day[7]  1.58   0.69       0.44       2.64   488 1.00
## a       4.22   0.41       3.55       4.83   388 1.00
## b1      0.04   0.03      -0.01       0.09   381 1.00
## b2     -0.07   0.01      -0.09      -0.05   449 1.01
## sigma   0.79   0.07       0.67       0.90   699 1.00

Diagnostics

plot(model3)
dashboard(model3)

Trace plots

Dashboard

Model Comparison

com<-compare(model1,model2,model3)
plot(com)

## Effect of Day of Week

wde<-function(x){
  s<-extract.samples(model2)
  ws<-exp(s$a+s$day[,7]+s$b1*x+s$b2*x^2)-exp(s$a+s$day[,3]+s$b1*x+s$b2*x^2)
  hist(ws,xlim=c(-100,800),ylim=c(0,500),xlab="Difference in time to reperfusion: Sunday-Wednesday",main=paste("Time: ",x,":00 hours"))
  text(600,400,paste("P(T<0):",length(ws[ws<0])/length(ws)))
  text(600,220,paste("median (IQR)",round(quantile(ws,c(0.5)),0),"(",round(quantile(ws,c(0.25)),0),",",round(quantile(ws,c(0.75)),0),")"))
}

The following counter-factual estimates the delay between a Sunday presentation and a Wednesay presentation with varying hours.

par(mfrow=c(3,2))
par(oma=c(0,0,0,0))
wde(2)
wde(9)
wde(12)
wde(16)
wde(20)
wde(23)

Time to Reperfusion

ubpdqn

October 5, 2016

Preamble

Data

Bayesian Models (based on McElreath Rethinking)

Model 1

Diagnostics

Model 2

Diagnostics

Predictive Posterior

Function

Model 3

Diagnostics

Model Comparison