Working through Chapter 8 of “Doing Bayesian Data Analysis” by John K. Kruschke.
A Complete Model
Load data
myData <- structure(list(y = c(0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L,
1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L)), .Names = "y", class = "data.frame", row.names = c(NA,
-50L))
y <- myData$y
Ntotal <- length(y)
dataList <- list(
y = myData$y,
Ntotal = length(y)
)
dplyr::glimpse(dataList)
## List of 2
## $ y : int [1:50] 0 1 0 0 0 0 0 0 0 0 ...
## $ Ntotal: int 50
Specify model
modelString <- "
model {
for (i in 1:Ntotal) {
y[i] ~ dbern(theta)
}
theta ~ dbeta(1, 1)
}
"
modelConnection <- textConnection(modelString)
Initialize chains
thetaInit <- sum(y)/length(y)
initsList <- list(theta = thetaInit)
#initsList <- function() {
# resampledY <- sample(y, replace = TRUE)
# thetaInit <- sum(resampledY)/length(resampledY)
# thetaInit <- 0.001 + 0.998 * thetaInit
# return(list(theta = thetaInit))
#}
Generate chains
library(rjags)
jagsModel <- jags.model(
file = modelConnection,
data = dataList,
inits = initsList,
n.chains = 3,
n.adapt = 500,
quiet = TRUE
)
close(modelConnection)
jagsModel
## JAGS model:
##
##
## model {
## for (i in 1:Ntotal) {
## y[i] ~ dbern(theta)
## }
## theta ~ dbeta(1, 1)
## }
##
## Fully observed variables:
## Ntotal y
update(jagsModel, n.iter = 500)
codaSamples <- coda.samples(
jagsModel,
variable.names = c("theta"),
n.iter = 3334
)
Examine chains
In using the functions in DBDA2E-utilities.R we have to edit the source of openGraph first and make its body empty. Can we do that here programmatically for reproducibility?
library(here)
setwd(here("DBDA2Eprograms"))
capture.output(source("DBDA2E-utilities.R")) -> x
diagMCMC(codaObject = codaSamples, parName = "theta")
Plot posterior
plotPost(codaSamples[, "theta"], main = "theta", xlab = bquote(theta))
## ESS mean median mode hdiMass hdiLow hdiHigh
## theta 10002 0.3084651 0.305687 0.3044928 0.95 0.1856722 0.4327306
## compVal pGtCompVal ROPElow ROPEhigh pLtROPE pInROPE pGtROPE
## theta NA NA NA NA NA NA NA
plotPost(
codaSamples[, "theta"],
main = "theta",
xlab = bquote(theta),
cenTend = "median",
compVal = 0.5,
ROPE = c(0.45, 0.55),
credMass = 0.9
)
## ESS mean median mode hdiMass hdiLow hdiHigh
## theta 10002 0.3084651 0.305687 0.3044928 0.9 0.2008087 0.4101135
## compVal pGtCompVal ROPElow ROPEhigh pLtROPE pInROPE pGtROPE
## theta 0.5 0.00259948 0.45 0.55 0.9815037 0.01839632 9.998e-05
Simplified scripts for frequently used analyses
setwd(here("DBDA2Eprograms"))
capture.output(source("Jags-Ydich-Xnom1subj-MbernBeta.R")) -> x
mcmcCoda <- genMCMC(data = myData, numSavedSteps = 10000)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 1
## Total graph size: 54
##
## Initializing model
##
## Burning in the MCMC chain...
## Sampling final MCMC chain...
diagMCMC(mcmcCoda, parName = "theta")
smryMCMC(mcmcCoda)
## Mean Median Mode ESS HDImass HDIlow HDIhigh
## theta 0.3067373 0.3042315 0.2957885 10000 0.95 0.1901645 0.4380788
## CompVal PcntGtCompVal ROPElow ROPEhigh PcntLtROPE PcntInROPE
## theta NA NA NA NA NA NA
## PcntGtROPE
## theta NA
## Mean Median Mode ESS HDImass HDIlow HDIhigh
## theta 0.3067373 0.3042315 0.2957885 10000 0.95 0.1901645 0.4380788
## CompVal PcntGtCompVal ROPElow ROPEhigh PcntLtROPE PcntInROPE
## theta NA NA NA NA NA NA
## PcntGtROPE
## theta NA
plotMCMC(mcmcCoda, data = myData)
Example: Difference of Biases
setwd(here("DBDA2Eprograms"))
myData <- read.csv("z6N8z2N7.csv")
s <- as.numeric(myData$s)
Ntotal <- length(y)
Nsubj <- length(unique(s))
dataList <- list(
y = y,
s = s,
Ntotal = Ntotal,
Nsubj = Nsubj
)
model {
for ( i in 1:Ntotal ) {
y[i] ~ dbern( theta[s[i]] )
}
for ( s in 1:Nsubj ) {
theta[s] ~ dbeta(2, 2)
}
}
setwd(here("DBDA2Eprograms"))
capture.output(source("Jags-Ydich-XnomSsubj-MbernBeta.R")) -> x
mcmcCoda <- genMCMC(data = myData, numSavedSteps = 10000)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 15
## Unobserved stochastic nodes: 2
## Total graph size: 38
##
## Initializing model
##
## Burning in the MCMC chain...
## Sampling final MCMC chain...
diagMCMC(mcmcCoda, parName = "theta[1]")
smryMCMC(mcmcCoda, compVal = NULL, compValDiff = 0.0)
## Mean Median Mode ESS HDImass HDIlow
## theta[1] 0.6671045 0.6759641 0.6864179 10000 0.95 0.4219037
## theta[2] 0.3631761 0.3548670 0.3267515 10000 0.95 0.1129291
## theta[1]-theta[2] 0.3039284 0.3138316 0.3260221 10000 0.95 -0.0601375
## HDIhigh CompVal PcntGtCompVal ROPElow ROPEhigh
## theta[1] 0.9118710 NA NA NA NA
## theta[2] 0.6414132 NA NA NA NA
## theta[1]-theta[2] 0.6713553 0 93.91 NA NA
## PcntLtROPE PcntInROPE PcntGtROPE
## theta[1] NA NA NA
## theta[2] NA NA NA
## theta[1]-theta[2] NA NA NA
## Mean Median Mode ESS HDImass HDIlow
## theta[1] 0.6671045 0.6759641 0.6864179 10000 0.95 0.4219037
## theta[2] 0.3631761 0.3548670 0.3267515 10000 0.95 0.1129291
## theta[1]-theta[2] 0.3039284 0.3138316 0.3260221 10000 0.95 -0.0601375
## HDIhigh CompVal PcntGtCompVal ROPElow ROPEhigh
## theta[1] 0.9118710 NA NA NA NA
## theta[2] 0.6414132 NA NA NA NA
## theta[1]-theta[2] 0.6713553 0 93.91 NA NA
## PcntLtROPE PcntInROPE PcntGtROPE
## theta[1] NA NA NA
## theta[2] NA NA NA
## theta[1]-theta[2] NA NA NA
plotMCMC(mcmcCoda, data = myData, compVal = NULL, compValDiff = 0.0)
Faster Sampling with Parallel Processing in runjags
library(runjags)
nChains <- 3
nAdaptSteps <- 1000
nBurninSteps <- 500
nUseSteps <- 10000
nThinSteps <- 2
library(runjags)
runJagsOut <- jagsModel <- run.jags(
method = "parallel",
model = here("model.txt"),
monitor = c("theta"),
data = dataList,
inits = initsList,
n.chains = nChains,
adapt = nAdaptSteps,
burnin = nBurninSteps,
sample = ceiling(nUseSteps/nChains),
thin = nThinSteps,
summarise = FALSE,
plots = FALSE
)
## Calling 3 simulations using the parallel method...
## Following the progress of chain 1 (the program will wait for all
## chains to finish before continuing):
## Welcome to JAGS 4.0.0 on Wed Sep 6 12:55:36 2017
## JAGS is free software and comes with ABSOLUTELY NO WARRANTY
## Loading module: basemod: ok
## Loading module: bugs: ok
## . . Reading data file data.txt
## . Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 1
## Total graph size: 70
##
## WARNING: Unused variable(s) in data table:
## Nsubj
## s
##
## . Reading parameter file inits1.txt
## . Initializing model
## . Adaptation skipped: model is not in adaptive mode.
## . Updating 500
## -------------------------------------------------| 500
## ************************************************** 100%
## . . Updating 6668
## -------------------------------------------------| 6650
## ************************************************** 100%
## * 100%
## . . . . Updating 0
## . Deleting model
## .
## All chains have finished
## Note: the model did not require adaptation
## Simulation complete. Reading coda files...
## Coda files loaded successfully
## Finished running the simulation
codaSamples <- as.mcmc.list(runJagsOut)
plot(codaSamples)
