params <-
list(EVAL = TRUE)

## ----SETTINGS-knitr, include=FALSE--------------------------------------------
stopifnot(require(knitr))
opts_chunk$set(
  comment=NA,
  eval = if (isTRUE(exists("params"))) params$EVAL else FALSE,
  dev = "png",
  dpi = 150,
  fig.asp = 0.618,
  fig.width = 5,
  out.width = "60%",
  fig.align = "center"
)

## ----eval=FALSE---------------------------------------------------------------
# library("rstan")
# 
# # Prepare data
# url <- "http://stat.columbia.edu/~gelman/arm/examples/arsenic/wells.dat"
# wells <- read.table(url)
# wells$dist100 <- with(wells, dist / 100)
# X <- model.matrix(~ dist100 + arsenic, wells)
# standata <- list(y = wells$switch, X = X, N = nrow(X), P = ncol(X))
# 
# # Fit model
# fit_1 <- stan("logistic.stan", data = standata)
# print(fit_1, pars = "beta")

## ----eval=FALSE---------------------------------------------------------------
# library("loo")
# 
# # Extract pointwise log-likelihood
# # using merge_chains=FALSE returns an array, which is easier to
# # use with relative_eff()
# log_lik_1 <- extract_log_lik(fit_1, merge_chains = FALSE)
# 
# # as of loo v2.0.0 we can optionally provide relative effective sample sizes
# # when calling loo, which allows for better estimates of the PSIS effective
# # sample sizes and Monte Carlo error
# r_eff <- relative_eff(exp(log_lik_1), cores = 2)
# 
# # preferably use more than 2 cores (as many cores as possible)
# # will use value of 'mc.cores' option if cores is not specified
# loo_1 <- loo(log_lik_1, r_eff = r_eff, cores = 2)
# print(loo_1)

## ----eval=FALSE---------------------------------------------------------------
# standata$X[, "arsenic"] <- log(standata$X[, "arsenic"])
# fit_2 <- stan(fit = fit_1, data = standata)
# 
# log_lik_2 <- extract_log_lik(fit_2, merge_chains = FALSE)
# r_eff_2 <- relative_eff(exp(log_lik_2))
# loo_2 <- loo(log_lik_2, r_eff = r_eff_2, cores = 2)
# print(loo_2)

## ----eval=FALSE---------------------------------------------------------------
# # Compare
# comp <- loo_compare(loo_1, loo_2)

## ----eval=FALSE---------------------------------------------------------------
# print(comp) # can set simplify=FALSE for more detailed print output