Data
## 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)
dat<-read.csv("old02.csv",header=TRUE)
datf<-data.frame(dat$DOW,dat$TOD,dat$LabToPCI)
colnames(datf)<-c("DOW","TOD","LabToPCI")
lbl<-c("1"="Monday",
"2"="Tuesday",
"3"="Wednesday",
"4"="Thursday",
"5"="Friday",
"6"="Saturday",
"7"="Sunday"
)
Visualization
ggplot(datf,aes(x=TOD,y=LabToPCI))+
geom_point()+
facet_wrap(~DOW,labeller = labeller(DOW=lbl))
ggplot(datf,aes(x=LabToPCI))+
geom_histogram(aes(y=..density..),
color="black",fill="red")+
facet_wrap(~DOW,labeller = labeller(DOW=lbl))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Model 1
m1<-map2stan(
alist(
LabToPCI~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=datf,
iter=1e4
)
##
## SAMPLING FOR MODEL 'LabToPCI ~ 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: 2.527 seconds (Warm-up)
## 2.888 seconds (Sampling)
## 5.415 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! 1
## [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 'LabToPCI ~ 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 ]
Diagnostics
precis(m1,depth=2)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 3.10 0.34 2.57 3.67 486 1
## day[1] 0.30 0.34 -0.23 0.87 525 1
## day[2] 0.46 0.34 -0.05 1.05 484 1
## day[3] 0.54 0.34 -0.01 1.09 478 1
## day[4] 0.43 0.34 -0.13 0.98 470 1
## day[5] 0.49 0.34 -0.08 1.01 470 1
## day[6] 0.49 0.34 -0.04 1.05 458 1
## day[7] 0.46 0.34 -0.07 1.03 495 1
## b -0.01 0.00 -0.02 0.00 2219 1
## sigma 0.40 0.02 0.36 0.43 2148 1
plot(m1)
dashboard(m1)
Prediction
pfun<-function(t){
ma<-matrix(data=NA,nrow=1000,ncol=7)
for(i in 1:1000){
for(j in 1:7){
ma[i,j]<-j
}
}
v1<-factor(as.vector(ma))
nd<-data.frame(DOW=1:7,TOD=t )
lf<-link(m1,data=nd)
v2<-exp(as.vector(lf))
gd<-data.frame(v1,v2)
colnames(gd)<-c("Day","LabToPCI")
library(plyr)
cdat<-ddply(gd, "Day", summarise, day.median=median(LabToPCI))
ggplot(gd,aes(x=LabToPCI))+
geom_histogram(aes(y=..density..),fill="red")+
geom_density(alpha=0.2,fill="#FF6666")+
xlim(0,75)+
xlab("Time (minutes)")+
geom_vline(data=cdat,aes(xintercept=day.median))+
geom_text(data=cdat,aes(x=60,y=0.1,
label=paste("Median: ",round(day.median))))+
facet_wrap(~Day,labeller = labeller(Day=lbl),ncol = 1)+
ggtitle(paste("Time: ",t,":00 hours"))
}
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Model 2
m2<-map2stan(
alist(
LabToPCI~dlnorm(mu,sigma),
mu<-a+day[DOW]+b1*TOD+b2*TOD*day[DOW],
a~dnorm(0,1),
day[DOW]~dnorm(0,1),
b1~dnorm(0,1),
b2~dnorm(0,1),
sigma~dcauchy(0,1)
),
data=datf,
iter=1e4
)
##
## SAMPLING FOR MODEL 'LabToPCI ~ 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: 6.453 seconds (Warm-up)
## 7.226 seconds (Sampling)
## 13.679 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! 5
## [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 'LabToPCI ~ 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 ]
Diagnostics
precis(m2,depth=2)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 3.10 0.35 2.56 3.65 424 1
## day[1] 0.42 0.41 -0.22 1.06 503 1
## day[2] 0.47 0.38 -0.12 1.06 490 1
## day[3] 0.45 0.37 -0.16 1.03 488 1
## day[4] 0.49 0.37 -0.08 1.07 474 1
## day[5] 0.41 0.37 -0.24 0.94 494 1
## day[6] 0.51 0.37 -0.03 1.13 464 1
## day[7] 0.37 0.37 -0.18 0.97 470 1
## b1 0.02 0.03 -0.02 0.07 512 1
## b2 -0.07 0.03 -0.11 -0.02 816 1
## sigma 0.40 0.02 0.37 0.43 1718 1
plot(m2)
dashboard(m2)
Compare Models
cm<-compare(m1,m2)
cm
## WAIC pWAIC dWAIC weight SE dSE
## m1 1600.2 8.1 0.0 0.89 26.99 NA
## m2 1604.3 8.8 4.1 0.11 27.09 3.61
plot(cm)
