Paul Johnson, CRMDA, KU <pauljohn@ku.edu>
Charles Redmon, CRMDA, KU <redmon@ku.edu>
Please visit http://crmda.ku.edu/guides
Keywords: multilevel, mixed effects, lme4, BUGS, Stan
May. 31, 2018
Abstract
Author, please REMEMBER TO INCLUDE AN ABSTRACT BEFORE FINALIZING THIS DOCUMENT!
The following packages will be used in this session:
library(data.table)
library(lme4)Loading required package: MatrixLoading required package: methodslibrary(R2OpenBUGS)
library(rstan)Loading required package: ggplot2Loading required package: StanHeadersrstan (Version 2.17.3, GitRev: 2e1f913d3ca3)For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
Attaching package: 'rstan'The following object is masked from 'package:R2OpenBUGS':
monitorIn order to make confidence interval plots in the style of Gelman & Hill (2007), we load this external function:
rdir <- "R"
source(file.path(rdir, "drawCI.R"))## set directories
wdir <- "workingdata"
bugsdir <- "OpenBUGS"
standir <- "stan"
## read in data
df <- readRDS(file.path(wdir, "srrs2.mn.rds"))
## create aggregate data frame, one row per county
dt <- as.data.table(df)
setkey(dt, "county.numbers")
ctyDT <- dt[ , list(sample.size = length(radonlog),
cty.mns = mean(radonlog, na.rm=TRUE),
cty.vars = var(radonlog, na.rm=TRUE),
county.name = unique(county.name)),
by = c("county.numbers")]
ctyDT[ , cty.sds.sep := sqrt(cty.vars/sample.size)]
head(ctyDT) county.numbers sample.size cty.mns cty.vars county.name
1: 1 4 0.6604064 0.2107976 AITKIN
2: 2 52 0.8332496 0.5933375 ANOKA
3: 3 3 1.0483380 0.5630051 BECKER
4: 4 7 1.1409857 0.9363348 BELTRAMI
5: 5 4 1.2524357 0.1801393 BENTON
6: 6 3 1.5130101 0.2663447 BIG STONE
cty.sds.sep
1: 0.2295635
2: 0.1068192
3: 0.4332071
4: 0.3657350
5: 0.2122141
6: 0.2979624## J is the number of counties
J <- length(unique(df$county.name))
## G&H create their own jittered variant for plotting
ctyDT[ , ss.jittered := sample.size * exp(runif(J, -.1, .1))]m1 <- lm(radonlog ~ 0 + county.name, data=df)
m1.pred <- predict(m1, se.fit=TRUE, newdata=ctyDT)Next we put the fit and standard error estimates onto the aggregate data frame.
ctyDT$m1.fit <- m1.pred$fit
ctyDT$m1.se <- m1.pred$se.fitPlot results:
m1ci2 <- drawCI(ctyDT$ss.jittered, ctyDT$m1.fit, ctyDT$m1.se,
main = "LM Estimates of Intercepts")
Identify cases with high values
ctyDT[ctyDT$m1.fit > 2, c("county.numbers", "county.name", "m1.fit")] county.numbers county.name m1.fit
1: 31 JACKSON 2.020569
2: 33 KANDIYOHI 2.061874
3: 36 LAC QUI PARLE 2.598696
4: 40 LINCOLN 2.135670
5: 50 MURRAY 2.493205
6: 51 NICOLLET 2.165039
7: 81 WATONWAN 2.229174
8: 82 WILKIN 2.230014ctyDT[ctyDT$m1.fit < 0.5, c("county.numbers", "county.name", "m1.fit")] county.numbers county.name m1.fit
1: 35 KOOCHICHING 0.4074616
2: 37 LAKE 0.3220285
3: 79 WASECA 0.4347479m1ci2 <- drawCI(ctyDT$ss.jittered, ctyDT$m1.fit, ctyDT$m1.se,
main="LM Estimates of Intercepts",
focal=c(36, 50, 81, 35, 37, 79), labelfocal=TRUE)
lme1 <- lmer(radonlog ~ (1 | county.numbers), data = df)
summary(lme1)Linear mixed model fit by REML ['lmerMod']
Formula: radonlog ~ (1 | county.numbers)
Data: df
REML criterion at convergence: 2259.4
Scaled residuals:
Min 1Q Median 3Q Max
-4.4661 -0.5734 0.0441 0.6432 3.3516
Random effects:
Groups Name Variance Std.Dev.
county.numbers (Intercept) 0.09581 0.3095
Residual 0.63662 0.7979
Number of obs: 919, groups: county.numbers, 85
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.31258 0.04891 26.84These are the random effect estimates:
lme1.ranef <- ranef(lme1, condVar = TRUE)Get the random effects in an easy-to-use format
lme1.ranef.df <- as.data.frame(lme1.ranef)## Insert the estimates into the county-level data
## fit is beta0 + random effect
ctyDT[ , "lme1.fit"] <- lme1.ranef.df$condval + fixef(lme1)[["(Intercept)"]]
## Check that lme1.fit is same as coef output
ctyDT[ , "lme1.coef"] <- coef(lme1)$county.numbers["(Intercept)"]
## Insert std errs from lmer in county data frame
ctyDT[ , "lme1.se"] <- lme1.ranef.df$condsdpar(mfrow=c(1,2))
m1ci <- drawCI(ctyDT$ss.jittered, ctyDT$m1.fit, ctyDT$m1.se,
main="LM Estimates of Intercepts",
focal=c(1, 36, 40, 77), labelfocal=TRUE)
lme1ci <- drawCI(ctyDT$ss.jittered, ctyDT$lme1.fit, ctyDT$lme1.se,
ylim=m1ci$ylim, main="lme4 Intercept Estimates",
focal=c(1, 36, 40, 77), labelfocal=TRUE)
Set up data:
mData <- list(N = nrow(df),
J=J,
y = df$radonlog,
county = df$county.numbers)Initial values:
radon.inits <- function (){
list (a = rnorm(J),
mu.a = rnorm(1),
sigma.y = runif(1),
sigma.a = runif(1))
}Run model:
radon.parameters <- c ("a", "mu.a", "sigma.y", "sigma.a")
## Run scratch files in temporary directory
bayeswd <- tempdir()
## Copy the BUGS code file to the bayes scratch dir
file.copy(file.path(bugsdir, "radon.multilevel.nopred.txt"),
bayeswd, overwrite=TRUE)[1] TRUEmBugs1 <- bugs(model.file="radon.multilevel.nopred.txt",
data=mData, inits=radon.inits, n.iter=1000, n.chains=3,
parameters.to.save=radon.parameters,
OpenBUGS.pgm="/usr/bin/OpenBUGS",
useWINE=FALSE, working.directory=bayeswd, clearWD=FALSE)Plot:
plot(mBugs1)
Save results:
ctyDT$bugs.fit <- mBugs1$median$a
ctyDT$bugs.se <- mBugs1$sd$aPlot comparison:
par(mfrow=c(1,2))
lme1ci <- drawCI(ctyDT$ss.jittered, ctyDT$lme1.fit, ctyDT$lme1.se,
main="lme4 Intercept Estimates", ylim=c(0.3, 2.2),
focal=c(36, 50, 81, 35, 37, 79), labelfocal=TRUE)
mBugs1ci <- drawCI(ctyDT$ss.jittered, ctyDT$bugs.fit, ctyDT$bugs.se,
ylim=lme1ci$ylim, main="OpenBUGS Intercepts",
focal=c(36, 50, 81, 35, 37, 79), labelfocal=TRUE)
Correlation:
plot(lme1.fit ~ bugs.fit, ctyDT)
mStan1 <- stan(file.path(standir, "radon-1.stan"),
data=mData, warmup=1000, iter=2000, chains=4)Warning in readLines(file, warn = TRUE): incomplete final line found
on '/home/pauljohn/GIT/CRMDA/workshops-lfs/topics/sem/sem-5/sem-5-2-
bayes_regression/writeup/stan/radon-1.stan'
SAMPLING FOR MODEL 'radon-1' NOW (CHAIN 1).
Gradient evaluation took 0.000126 seconds
1000 transitions using 10 leapfrog steps per transition would take 1.26 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.10452 seconds (Warm-up)
0.879396 seconds (Sampling)
1.98392 seconds (Total)
SAMPLING FOR MODEL 'radon-1' NOW (CHAIN 2).
Gradient evaluation took 5.2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.52 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.10706 seconds (Warm-up)
1.71077 seconds (Sampling)
2.81783 seconds (Total)
SAMPLING FOR MODEL 'radon-1' NOW (CHAIN 3).
Gradient evaluation took 5.2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.52 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.06736 seconds (Warm-up)
0.877279 seconds (Sampling)
1.94464 seconds (Total)
SAMPLING FOR MODEL 'radon-1' NOW (CHAIN 4).
Gradient evaluation took 5.3e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.53 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.09689 seconds (Warm-up)
1.09558 seconds (Sampling)
2.19247 seconds (Total)Extract predicted intercepts:
post <- extract(mStan1)
ints <- apply(post$u, 2, function(x) x + post$b0)Add fit to results dataset:
ctyDT$stan.fit <- apply(ints, 2, median)
ctyDT$stan.se <- apply(ints, 2, sd)
ctyDT$stan.mean <- apply(ints, 2, mean)Plot comparison:
par(mfrow=c(1,3))
lme1ci <- drawCI(ctyDT$ss.jittered, ctyDT$lme1.fit, ctyDT$lme1.se,
main="lme4 Intercept Estimates", ylim=c(0.3, 2.2),
focal=c(36, 50, 81, 35, 37, 79), labelfocal=TRUE)
mBugs1ci <- drawCI(ctyDT$ss.jittered, ctyDT$bugs.fit, ctyDT$bugs.se,
ylim=lme1ci$ylim, main="OpenBUGS Intercepts",
focal=c(36, 50, 81, 35, 37, 79), labelfocal=TRUE)
mStan1ci <- drawCI(ctyDT$ss.jittered, ctyDT$stan.fit, ctyDT$stan.se,
ylim=lme1ci$ylim, main="Stan Intercepts",
focal=c(36, 50, 81, 35, 37, 79), labelfocal=TRUE)
Correlation matrix:
pairs(~ lme1.fit + bugs.fit + stan.fit, ctyDT)
print(mStan1, pars="sigmaU")Inference for Stan model: radon-1.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigmaU 0.32 0 0.05 0.23 0.28 0.31 0.35 0.42 1224 1
Samples were drawn using NUTS(diag_e) at Thu May 31 11:24:11 2018.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).mStan2 <- stan(file.path(standir, "radon-1b.stan"),
data=mData, warmup=1000, iter=2000, chains=4)
SAMPLING FOR MODEL 'radon-1b' NOW (CHAIN 1).
Gradient evaluation took 8.5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.85 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 0.96185 seconds (Warm-up)
0.691232 seconds (Sampling)
1.65308 seconds (Total)
SAMPLING FOR MODEL 'radon-1b' NOW (CHAIN 2).
Gradient evaluation took 5.9e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.59 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.15808 seconds (Warm-up)
0.448342 seconds (Sampling)
1.60642 seconds (Total)
SAMPLING FOR MODEL 'radon-1b' NOW (CHAIN 3).
Gradient evaluation took 5.7e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.57 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 0.912232 seconds (Warm-up)
0.579736 seconds (Sampling)
1.49197 seconds (Total)
SAMPLING FOR MODEL 'radon-1b' NOW (CHAIN 4).
Gradient evaluation took 5.2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.52 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.00278 seconds (Warm-up)
0.353058 seconds (Sampling)
1.35584 seconds (Total)Warning: There were 3917 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupWarning: Examine the pairs() plot to diagnose sampling problemsmStan3 <- stan(file.path(standir, "radon-1c.stan"),
data=mData, warmup=1000, iter=2000, chains=4)Warning in readLines(file, warn = TRUE): incomplete final line found
on '/home/pauljohn/GIT/CRMDA/workshops-lfs/topics/sem/sem-5/sem-5-2-
bayes_regression/writeup/stan/radon-1c.stan'
SAMPLING FOR MODEL 'radon-1c' NOW (CHAIN 1).
Gradient evaluation took 8.9e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.89 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.08651 seconds (Warm-up)
1.26972 seconds (Sampling)
2.35624 seconds (Total)
SAMPLING FOR MODEL 'radon-1c' NOW (CHAIN 2).
Gradient evaluation took 5.5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.55 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.15757 seconds (Warm-up)
0.68128 seconds (Sampling)
1.83885 seconds (Total)
SAMPLING FOR MODEL 'radon-1c' NOW (CHAIN 3).
Gradient evaluation took 5.8e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.58 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.05829 seconds (Warm-up)
1.85484 seconds (Sampling)
2.91313 seconds (Total)
SAMPLING FOR MODEL 'radon-1c' NOW (CHAIN 4).
Gradient evaluation took 5.7e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.57 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.14459 seconds (Warm-up)
0.449589 seconds (Sampling)
1.59418 seconds (Total)Warning: There were 3985 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Warning: Examine the pairs() plot to diagnose sampling problemspost2 <- extract(mStan2)
ints2 <- apply(post2$u, 2, function(x) x + post2$b0)
post3 <- extract(mStan3)
ints3 <- apply(post3$u, 2, function(x) x + post3$b0)Add fit to results dataset:
ctyDT$stan2.fit <- apply(ints2, 2, median)
ctyDT$stan2.se <- apply(ints2, 2, sd)
ctyDT$stan2.mean <- apply(ints2, 2, mean)
ctyDT$stan3.fit <- apply(ints3, 2, median)
ctyDT$stan3.se <- apply(ints3, 2, sd)
ctyDT$stan3.mean <- apply(ints3, 2, mean)Plot comparison:
par(mfrow=c(2, 2))
m1ci <- drawCI(ctyDT$ss.jittered, ctyDT$m1.fit, ctyDT$m1.se,
main="LM baseline", ylim=c(-0.5, 3.5),
focal=c(36, 50, 81, 35, 37, 79), labelfocal=TRUE)
stan1ci <- drawCI(ctyDT$ss.jittered, ctyDT$stan.fit, ctyDT$stan.se,
main="Stan baseline", ylim=c(-0.5, 3.5),
focal=c(36, 50, 81, 35, 37, 79), labelfocal=TRUE)
stan2ci <- drawCI(ctyDT$ss.jittered, ctyDT$stan2.fit, ctyDT$stan2.se,
ylim=c(-0.5, 3.5), main="SD on U Prior: 0.6",
focal=c(36, 50, 81, 35, 37, 79), labelfocal=TRUE)
stan3ci <- drawCI(ctyDT$ss.jittered, ctyDT$stan3.fit, ctyDT$stan3.se,
ylim=c(-0.5, 3.5), main="SD on U Prior: 1.0",
focal=c(36, 50, 81, 35, 37, 79), labelfocal=TRUE)
With uniform prior more iterations may be needed.
mStan4 <- stan(file.path(standir, "radon-1d.stan"),
data=mData, warmup=1000, iter=2000, chains=4)
SAMPLING FOR MODEL 'radon-1d' NOW (CHAIN 1).
Gradient evaluation took 8.2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.82 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.12818 seconds (Warm-up)
1.69523 seconds (Sampling)
2.82341 seconds (Total)
SAMPLING FOR MODEL 'radon-1d' NOW (CHAIN 2).
Gradient evaluation took 5.5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.55 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.43812 seconds (Warm-up)
0.531661 seconds (Sampling)
1.96978 seconds (Total)
SAMPLING FOR MODEL 'radon-1d' NOW (CHAIN 3).
Gradient evaluation took 7.8e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.78 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.2872 seconds (Warm-up)
0.552427 seconds (Sampling)
1.83963 seconds (Total)
SAMPLING FOR MODEL 'radon-1d' NOW (CHAIN 4).
Gradient evaluation took 5.7e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.57 seconds.
Adjust your expectations accordingly!
Iteration: 1 / 2000 [ 0%] (Warmup)
Iteration: 200 / 2000 [ 10%] (Warmup)
Iteration: 400 / 2000 [ 20%] (Warmup)
Iteration: 600 / 2000 [ 30%] (Warmup)
Iteration: 800 / 2000 [ 40%] (Warmup)
Iteration: 1000 / 2000 [ 50%] (Warmup)
Iteration: 1001 / 2000 [ 50%] (Sampling)
Iteration: 1200 / 2000 [ 60%] (Sampling)
Iteration: 1400 / 2000 [ 70%] (Sampling)
Iteration: 1600 / 2000 [ 80%] (Sampling)
Iteration: 1800 / 2000 [ 90%] (Sampling)
Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.23697 seconds (Warm-up)
0.904189 seconds (Sampling)
2.14116 seconds (Total)Warning: There were 3999 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupWarning: Examine the pairs() plot to diagnose sampling problemspost4 <- extract(mStan4)
ints4 <- apply(post4$u, 2, function(x) x + post4$b0)Add fit to results dataset:
ctyDT$stan4.fit <- apply(ints4, 2, median)
ctyDT$stan4.se <- apply(ints4, 2, sd)Plot comparison:
par(mfrow=c(1, 2))
m1ci <- drawCI(ctyDT$ss.jittered, ctyDT$m1.fit, ctyDT$m1.se,
main="LM baseline", ylim=c(-0.5, 3.5),
focal=c(36, 50, 81, 35, 37, 79), labelfocal=TRUE)
stan4ci <- drawCI(ctyDT$ss.jittered, ctyDT$stan.fit, ctyDT$stan.se,
main="Uniform Prior on U", ylim=c(-0.5, 3.5),
focal=c(36, 50, 81, 35, 37, 79), labelfocal=TRUE)
R version 3.4.4 (2018-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] methods stats graphics grDevices utils datasets
[7] base
other attached packages:
[1] rstan_2.17.3 StanHeaders_2.17.2 ggplot2_2.2.1
[4] R2OpenBUGS_3.2-3.2 lme4_1.1-17 Matrix_1.2-14
[7] data.table_1.10.4-3 stationery_0.80
loaded via a namespace (and not attached):
[1] lavaan_0.5-23.1097 splines_3.4.4 lattice_0.20-35
[4] colorspace_1.3-2 htmltools_0.3.6 stats4_3.4.4
[7] yaml_2.1.16 mgcv_1.8-23 rlang_0.1.6
[10] nloptr_1.0.4 pillar_1.1.0 foreign_0.8-69
[13] plyr_1.8.4 stringr_1.2.0 MatrixModels_0.4-1
[16] munsell_0.4.3 rockchalk_1.8.111 gtable_0.2.0
[19] codetools_0.2-15 kutils_1.45 coda_0.19-1
[22] evaluate_0.10.1 inline_0.3.14 knitr_1.19
[25] SparseM_1.77 quantreg_5.35 pbkrtest_0.4-7
[28] parallel_3.4.4 Rcpp_0.12.15 xtable_1.8-2
[31] backports_1.1.2 scales_0.5.0 gridExtra_2.3
[34] mnormt_1.5-5 digest_0.6.15 stringi_1.2.2
[37] openxlsx_4.0.17 grid_3.4.4 rprojroot_1.3-2
[40] quadprog_1.5-5 tools_3.4.4 magrittr_1.5
[43] lazyeval_0.2.1 tibble_1.4.2 car_2.1-6
[46] pbivnorm_0.6.0 MASS_7.3-49 minqa_1.2.4
[49] rmarkdown_1.8 boot_1.3-20 nnet_7.3-12
[52] nlme_3.1-137 compiler_3.4.4 Available under Created Commons license 3.0 