Registrar Training
ubpdqn
October 20, 2016
Data
DAP for diagnostic coronary angiographywas collected for first and second year advanced trainees in Cardiology after similar periods of training. This explores “learning effect”.
setwd("e:/pandoc/rmd")
d<-read.csv("reg.csv")
d<-data.frame(d)
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)Visualization
ggplot(d,aes(x=DAP))+geom_histogram(aes(y=..density..),fill="red",color="black")+facet_grid(Experience~Access)+ggtitle("Dose Area Product Diagnostic Coronary Angiograpy\nby access type and experience")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Models
Model 1
m1<-map2stan(
alist(
DAP~dlnorm(mu,sigma),
mu<-a+experience[Experience]+access[Access],
a~dnorm(1,1),
experience[Experience]~dnorm(0,1),
access[Access]~dnorm(0,1),
sigma~dcauchy(0,1)
),
data=d,
iter=1e4
)##
## SAMPLING FOR MODEL 'DAP ~ 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: 15.989 seconds (Warm-up)
## 19.78 seconds (Sampling)
## 35.769 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times after warmup on chain 1"
## count
## Exception thrown at line 24: lognormal_log: Scale parameter is inf, but must be finite! 6
## [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 'DAP ~ 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
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 4.69 0.67 3.59 5.73 1316 1
## experience[1] 1.90 0.62 0.96 2.91 1231 1
## experience[2] 1.75 0.62 0.82 2.78 1239 1
## access[1] 1.74 0.60 0.89 2.79 1342 1
## access[2] 1.95 0.60 1.12 3.02 1351 1
## sigma 0.55 0.02 0.52 0.57 1838 1
Model 2
m2<-map2stan(
alist(
DAP~dlnorm(mu,sigma),
mu<-a+experience[Experience]+access[Access]+experience[Experience]*access[Access],
a~dnorm(1,1),
experience[Experience]~dnorm(0,1),
access[Access]~dnorm(0,1),
sigma~dcauchy(0,1)
),
data=d,
iter=1e4
)##
## SAMPLING FOR MODEL 'DAP ~ 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: 39.644 seconds (Warm-up)
## 47.882 seconds (Sampling)
## 87.526 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times after warmup on chain 1"
## count
## Exception thrown at line 24: lognormal_log: Scale parameter is inf, but must be finite! 8
## [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 'DAP ~ 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
p2<-precis(m2,depth = 2)
p2## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 2.45 0.94 0.98 3.98 1645 1
## experience[1] 1.70 0.43 1.02 2.38 1047 1
## experience[2] 1.64 0.42 0.98 2.31 1050 1
## access[1] 1.59 0.41 0.92 2.24 1068 1
## access[2] 1.67 0.42 1.01 2.36 1060 1
## sigma 0.55 0.02 0.52 0.57 1713 1plot(p2)
plot(m2)
dashboard(m2)
Compare Models
cm<-compare(m1,m2)
cm## WAIC pWAIC dWAIC weight SE dSE
## m2 10927.9 4 0 0.5 40.73 NA
## m1 10927.9 4 0 0.5 40.82 0.24plot(cm)
Predictions
Model 1 is used in the following. The learning effect is supported for both femoral and radial access.
sam<-extract.samples(m1,n=5000)
fy1<-exp(sam$a+sam$access[,1]+sam$experience[,1])
fy2<-exp(sam$a+sam$access[,1]+sam$experience[,2])
lef<-fy1-fy2
hist(lef,main="Learning Effect: Femoral Access ",xlab="DAP: Year1 -Year2")
text(1100,500,paste("P(change in DAP< 0): ",round(length(lef[lef<0])/5000,3)))
sam<-extract.samples(m1,n=5000)
ry1<-exp(sam$a+sam$access[,2]+sam$experience[,1])
ry2<-exp(sam$a+sam$access[,2]+sam$experience[,2])
ler<-ry1-ry2
hist(ler,main="Learning Effect: Radial Access ",xlab="DAP: Year1 -Year2")
text(1600,600,paste("P(change in DAP< 0): ",round(length(ler[ler<0])/5000,3)))
a1<-ry1-fy1
hist(a1,main="Learning Effect: Radial v Femoral (First year) ",xlab="DAP: Radial - Femoral")
text(1500,700,paste("P(change in DAP< 0): ",round(length(a1[a1<0])/5000,3)))
a2<-ry2-fy2
hist(a2,main="Learning Effect: Radial v Femoral (Second year) ",xlab="DAP: Radial - Femoral")
text(1250,700,paste("P(change in DAP< 0): ",round(length(a2[a2<0])/5000,3)))
## Mean 2.5% 25% 50% 75%
## Learning effect: femoral 619.2913 122.5528 430.7119 609.5256 794.1793
## Learning effect: radial 766.8394 150.7434 534.7962 762.5665 978.5176
## Femoral v radial (year 1) 999.6646 602.5915 854.7235 996.4276 1140.4839
## Femoral v radial (year2) 852.1165 519.3226 730.6903 849.0985 974.9072
## 97.5%
## Learning effect: femoral 1170.894
## Learning effect: radial 1441.337
## Femoral v radial (year 1) 1421.792
## Femoral v radial (year2) 1213.560
Summary
This Bayesian analysis has findings consistent with classical analysis including quantile regression
Cross sectional analysis of experienced v novice trainees are consistent with a learning effect for both femoral and radial approaches. (Longitudinal analysis is required)
Radial procedures have higher dose than femoral procedures. The gap is narroiwer in the second year trainees.
Other factors (number of procedures, patient characteristics ) were not available and limits.