Time to Reperfusion 2
ubpdqn
October 7, 2016
Preamble
This is a second analysis using a previous dataset.
Data
The times collected and definitions used in this dataset are different from the previous analysis. * Door to Lab* time and Lab to PCI time are outcomes. The weekend and time of day relationships will again be explored.
Door to PCI Time
t1<-read.csv("old01.csv",header=TRUE)
library(rethinking)## Loading required package: rstan## Loading required package: ggplot2## Loading required package: StanHeaders## rstan (Version 2.12.1, packaged: 2016-09-11 13:07:50 UTC, GitRev: 85f7a56811da)## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())## Loading required package: parallel## rethinking (Version 1.59)Visualizing the Data
In the following 1-Monday,2-Tuesday,…,Sunday-7.
ggplot(t1, aes(x =TOD, y =DoorToPCI, colour=DOW)) +
geom_point()+
facet_wrap( ~DOW)
Door to PCI time by time of data and gropuped by day of week
Histograms
ggplot(t1, aes(x=DoorToPCI)) +
geom_histogram(aes(y=..density..),fill="red")+geom_density(alpha=0.5)+
facet_wrap( ~DOW)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Histograms of Door To PCI time by day of week.
Model 1
The model:
\[ t \sim \text{LogNormal}(\mu,\sigma)\\ \mu=\alpha + \text{day}_\text{DOW}+\beta_1 \ \text{TOD} + \beta_2\ \text{TOD}^2 \\ \alpha\sim \text{Normal}(0,1)\\ \text{day}_\text{DOW}\sim \text{Normal}(0,1)\\ \beta_1\sim \text{Norma}(0,1)\\ \beta_2\sim\text{Normal}(0,1)\\ \sigma\sim\text{Cauchy}(0,1)\\ \]
d2pci1<-map2stan(
alist(DoorToPCI~dlnorm(mu,sigma),
mu<-a+day[DOW]+b1*TOD+b2*TOD*TOD,
a~dnorm(0,1),
day[DOW]~dnorm(0,1),
b1~dnorm(0,1),
b2~dnorm(0,1),
sigma~dcauchy(0,1)
),
data=t1,iter=1e4
)##
## SAMPLING FOR MODEL 'DoorToPCI ~ dlnorm(mu, sigma)' 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: 11.23 seconds (Warm-up)
## 5.908 seconds (Sampling)
## 17.138 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times after warmup on chain 1"
## count
## Exception thrown at line 25: lognormal_log: Scale parameter is inf, but must be finite! 10
## [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 'DoorToPCI ~ 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)## 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 ]precis(d2pci1,depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 4.07 0.35 3.54 4.67 616 1
## day[1] 0.59 0.35 0.04 1.14 623 1
## day[2] 0.46 0.35 -0.10 1.01 647 1
## day[3] 0.65 0.35 0.09 1.20 627 1
## day[4] 0.50 0.35 -0.03 1.08 622 1
## day[5] 0.62 0.35 0.09 1.21 628 1
## day[6] 0.74 0.35 0.18 1.30 630 1
## day[7] 0.66 0.35 0.11 1.24 638 1
## b1 -0.08 0.02 -0.12 -0.04 1614 1
## b2 0.00 0.00 0.00 0.00 1667 1
## sigma 0.50 0.02 0.46 0.54 1583 1Diagnostics
plot(d2pci1)
dashboard(d2pci1)
Model 2
The model
d2pci2<-map2stan(
alist(DoorToPCI~dlnorm(mu,sigma),
mu<-a+day[DOW]+gamma*TOD,
gamma<-b1+b2*TOD+b3*day[DOW],
a~dnorm(0,1),
day[DOW]~dnorm(0,1),
b1~dnorm(0,1),
b2~dnorm(0,1),
b3~dnorm(0,1),
sigma~dcauchy(0,1)
),
data=t1,iter=1e4
)##
## SAMPLING FOR MODEL 'DoorToPCI ~ dlnorm(mu, sigma)' 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: 14.758 seconds (Warm-up)
## 11.81 seconds (Sampling)
## 26.568 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times after warmup on chain 1"
## count
## Exception thrown at line 31: lognormal_log: Scale parameter is 0, but must be > 0! 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 'DoorToPCI ~ 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.001 seconds (Sampling)
## 0.001 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 ]precis(d2pci2,depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 4.12 0.39 3.51 4.71 766 1
## day[1] 0.54 0.39 -0.05 1.19 860 1
## day[2] 0.65 0.44 0.01 1.41 769 1
## day[3] 0.58 0.39 -0.01 1.23 838 1
## day[4] 0.62 0.41 -0.04 1.25 812 1
## day[5] 0.50 0.40 -0.16 1.13 852 1
## day[6] 0.62 0.41 -0.03 1.26 807 1
## day[7] 0.60 0.40 -0.04 1.23 833 1
## b1 -0.05 0.04 -0.11 0.02 793 1
## b2 0.00 0.00 0.00 0.00 2607 1
## b3 -0.07 0.03 -0.12 -0.02 527 1
## sigma 0.51 0.03 0.47 0.55 2302 1plot(d2pci2)
dashboard(d2pci2)
Model 3
d2pci3<-map2stan(
alist(DoorToPCI~dlnorm(mu,sigma),
mu<-a+day[DOW]+gamma*TOD,
gamma<-b1+b2*day[DOW],
a~dnorm(0,1),
day[DOW]~dnorm(0,1),
b1~dnorm(0,1),
b2~dnorm(0,1),
sigma~dcauchy(0,1)
),
data=t1,iter=1e4
)##
## SAMPLING FOR MODEL 'DoorToPCI ~ dlnorm(mu, sigma)' 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: 5.294 seconds (Warm-up)
## 6.063 seconds (Sampling)
## 11.357 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 'DoorToPCI ~ 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)## 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 ]precis(d2pci3)## 7 vector or matrix parameters omitted in display. Use depth=2 to show them.## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 3.75 0.35 3.16 4.29 554 1
## b1 0.02 0.03 -0.02 0.06 571 1
## b2 -0.05 0.03 -0.11 -0.01 499 1
## sigma 0.52 0.03 0.48 0.56 2037 1plot(d2pci3)
dashboard(d2pci3)
Compare Models
cp<-compare(d2pci1,d2pci2,d2pci3)
cp## WAIC pWAIC dWAIC weight SE dSE
## d2pci1 2110.4 9.4 0.0 0.96 26.41 NA
## d2pci2 2116.6 12.0 6.1 0.04 27.11 3.98
## d2pci3 2128.4 10.8 18.0 0.00 27.91 6.98plot(cp)
Predictions
mf<-function(x,y,t)
{
s<-extract.samples(d2pci1,n=5000)
diff<-exp(s$a+s$day[,x]+s$b1*t+s$b2*t*t)-
exp(s$a+s$day[,y]+s$b1*t+s$b2*t*t)
df<-data.frame(diff,names=c("Difference"))
ggplot(df,aes(x=diff))+geom_histogram(aes(y=..density..) , fill="red",color="black")+
geom_text(mapping=aes(x=45,y=0.04,
label=paste("P(T<0): ",round(length(diff[diff<0])/5000,2))),
color="blue",size=4)+
geom_text(mapping=aes(x=45,y=0.03,label=paste("Median:",round(quantile(diff,c(0.5))))),
color="blue",size=4)+
scale_x_continuous(limits = c(-60, 60))+
scale_y_continuous(limits = c(0,0.05))+
ggtitle(paste("Difference between day", x, "and day ",y,"at ",t,":00"))+
xlab("Time (minutes)")+theme(plot.title=element_text(size=12))
}Prediction from Model 1
todf<-function(t){
library(gridExtra)
grid.arrange(mf(7,1,t),
mf(7,2,t),
mf(7,3,t),
mf(7,4,t),
mf(7,5,t),
mf(7,6,t),
ncol=2,nrow=3)
}2 am
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 2 rows containing non-finite values (stat_bin).## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 3 rows containing non-finite values (stat_bin).## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 1 rows containing non-finite values (stat_bin).
8 am
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Midday
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
4 pm
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
8 pm
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
11 pm
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.