params <-
list()

## ----setup, include=FALSE-----------------------------------
library(IJSE)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", dpi = 300,
  fig.align = "center", out.width = "100%",
  dev = if (knitr::is_latex_output()) "cairo_pdf" else "png",
  dev.args = if (knitr::is_latex_output()) list() else list(type = "cairo"))
options(knitr.kable.NA = "-")
data_dir <- params$data_dir
if (is.null(data_dir)) data_dir <- system.file("extdata", "study", package = "IJSE")
stopifnot(dir.exists(data_dir))
read_study <- function(name) read.csv(file.path(data_dir, name), stringsAsFactors = FALSE)
normal <- read_study("normal-summary.csv")
design <- read_study("normal-design.csv")
audit <- read_study("draw-audit-summary.csv")
quantile <- read_study("quantile-summary.csv")
accounting <- read_study("quantile-accounting.csv")
paired <- read_study("quantile-paired-coverage.csv")
source(system.file("study", "study-figures.R", package = "IJSE"), local = TRUE)
table_out <- function(x, caption, digits = 3) {
  knitr::kable(x, caption = caption, digits = digits, row.names = FALSE,
    format = if (knitr::is_latex_output()) "latex" else "html", booktabs = TRUE,
    longtable = FALSE, position = "H")
}
percent <- function(x, digits = 1) sprintf(paste0("%.", digits, "f"), 100 * x)
value <- function(scenario, method, column = "coverage") {
  normal[normal$scenario == scenario & normal$method == method, column]
}
qvalue <- function(scenario, parameter, method, column = "coverage") {
  quantile[quantile$scenario == scenario & quantile$parameter == parameter &
             quantile$method == method, column]
}
short_names <- c("Matched (50)", "Matched (200)", "Variance x4 (50)",
 "Variance x4 (200)", "t5 variance x4 (200)", "8 equal clusters",
 "30 equal clusters", "100 equal clusters", "30 unequal clusters", "Strong prior (20)")

## ----short-example------------------------------------------
set.seed(417)
y <- rnorm(200, sd = 2)
v <- 1 / (length(y) + 0.01)
m <- v * sum(y)
mu <- rnorm(6000, m, sqrt(v))
draws <- cbind(mu = mu)
ll <- vapply(y, function(yi) dnorm(yi, mu, log = TRUE),
             numeric(length(mu)))
result <- ijse_draws(draws, ll, diagnostics = "none")
summary(result)$coefficients[, 1:3]
vcov(result)
vcov(result, type = "posterior")

## ----design-table, results='asis', echo=FALSE---------------
d <- data.frame(Scenario = short_names, N = design$n,
                J = ifelse(design$clusters == 0, "-", design$clusters),
                Variance = design$true_variance, `Prior SD` = design$prior_sd,
                check.names = FALSE)
table_out(d, "Part A design. Every row uses 4,000 independently generated datasets.")

## ----normal-coverage-table, results='asis', echo=FALSE------
rows <- lapply(seq_len(nrow(design)), function(i) {
  s <- design$scenario[i]
  entry <- function(method) {
    z <- normal[normal$scenario == s & normal$method == method, ]
    if (!nrow(z)) return("-")
    paste0(percent(z$coverage), " (", percent(z$coverage_mcse, 2), ")")
  }
  data.frame(Scenario = short_names[i], Posterior = entry("Posterior"),
               `IJ (ind.)` = entry("IJ independent"), `IJ (cl.)` = entry("IJ cluster"),
               check.names = FALSE)
})
table_out(do.call(rbind, rows),
 "Coverage percentages, with Monte Carlo SE in percentage points in parentheses. Cluster IJ uses normal intervals.")

## ----normal-coverage-figure, echo=FALSE, fig.width=7, fig.height=6, fig.cap="Part A coverage. Points summarize 4,000 replications per scenario; horizontal bars are Wilson 95% intervals for simulation coverage. The dashed vertical line marks 95%. The t reference is a sensitivity analysis and is not the package default."----
normal_coverage_plot(data_dir)

## ----calibration-table, echo=FALSE, results='asis'----------
rows <- lapply(seq_len(nrow(design)), function(i) {
  s <- design$scenario[i]
  base <- normal[normal$scenario == s & normal$method == "Posterior", ]
  data.frame(Scenario = short_names[i], Bias = base$bias, `Emp. SD` = base$empirical_sd,
    `Post/SD` = base$se_ratio, `IJ ind/SD` = value(s, "IJ independent", "se_ratio"),
    `IJ cl/SD` = if (design$clusters[i]) value(s, "IJ cluster", "se_ratio") else NA_real_,
    check.names = FALSE)
})
table_out(do.call(rbind, rows), "Bias, empirical sampling SD, and mean estimated uncertainty divided by empirical SD.")

## ----draw-audit-table, echo=FALSE, results='asis'-----------
d <- data.frame(Draws = audit$draws, Ensembles = audit$ensembles,
  `Mean error (%)` = 100 * audit$mean_relative_error,
  `RMSE (%)` = 100 * audit$relative_rmse,
  `Theory SD (%)` = 100 * audit$theoretical_relative_sd, check.names = FALSE)
table_out(d, "Numerical error relative to the exact IJ SE. Each row contains 80 independent posterior-draw ensembles.", digits = 2)

## ----draw-precision-figure, echo=FALSE, fig.width=7, fig.height=3.8, fig.cap="Part B numerical precision. Each dot is one of 80 independent posterior-draw ensembles. Boxes show the median and interquartile range; whiskers reach the most extreme observations within 1.5 interquartile ranges of the box. Dashed curves mark the exact IID-normal 2.5th and 97.5th percentiles of relative SE error, not confidence limits for its mean."----
draw_precision_plot(data_dir)

## ----quantile-accounting, echo=FALSE, results='asis'--------
d <- accounting[c("scenario", "datasets", "completed", "adequate", "retries")]
names(d) <- c("Scenario", "Datasets", "Completed", "Adequate", "Retried")
table_out(d, "Actual brms fit accounting. Retried datasets contribute their final fit.", digits = 0)

## ----quantile-calibration, echo=FALSE, results='asis'-------
keys <- unique(quantile[c("scenario", "parameter")])
rows <- lapply(seq_len(nrow(keys)), function(i) {
  key <- keys[i, ]
  z <- quantile[quantile$scenario == key$scenario & quantile$parameter == key$parameter, ]
  take <- function(method, column) {
    x <- z[z$method == method, column]
    if (!length(x)) NA_real_ else x
  }
  data.frame(Scenario = key$scenario, Term = key$parameter,
    `Emp. SD` = z$empirical_sd[1], `Post/SD` = take("Posterior", "se_ratio"),
    `IJ ind/SD` = take("IJ independent", "se_ratio"),
    `IJ cl/SD` = take("IJ cluster", "se_ratio"), check.names = FALSE)
})
table_out(do.call(rbind, rows), "Actual brms fits: empirical SD and mean SE-to-SD ratios for both regression coefficients.")

## ----quantile-coverage-table, echo=FALSE, results='asis'----
rows <- lapply(seq_len(nrow(keys)), function(i) {
  key <- keys[i, ]
  z <- quantile[quantile$scenario == key$scenario & quantile$parameter == key$parameter, ]
  entry <- function(method) {
    x <- z[z$method == method, ]
    if (!nrow(x)) return("-")
    paste0(percent(x$coverage), " (", percent(x$coverage_mcse, 2), ")")
  }
  data.frame(Scenario = key$scenario, Term = key$parameter,
    Posterior = entry("Posterior"), `IJ (ind.)` = entry("IJ independent"),
    `IJ (cl.)` = entry("IJ cluster"), check.names = FALSE)
})
table_out(do.call(rbind, rows), "Actual brms coverage percentages, with Monte Carlo SE in percentage points in parentheses.")

## ----quantile-coverage-figure, echo=FALSE, fig.width=7, fig.height=3.8, fig.cap="Part C: actual brms median-regression fits, with 200 datasets per scenario. Points show empirical coverage and bars show Wilson 95% intervals for coverage. The dashed line marks 95%. All IJ intervals use a normal reference. These error bars quantify simulation uncertainty, not MCMC uncertainty within a fit."----
quantile_coverage_plot(data_dir)

## ----brms-example-code, eval=FALSE--------------------------
# fit <- readRDS("study-output/quantile-clustered-example.rds")
# stopifnot("group" %in% names(fit$data))
# result <- ijse(fit, cluster_var = fit$data$group,
#                save_influence = TRUE, diagnostics = "all")
# stopifnot(result$clustered, result$n_units == 30L)
# summary(result)
# vcov(result)
# vcov(result, type = "posterior")
# case_weight_derivatives <- influence(result) / result$n_units

## ----worked-example-table, echo=FALSE, results='asis'-------
d <- read_study("worked-example.csv")
d <- d[c("parameter", "Estimate", "Posterior_SD", "IJ_SE")]
names(d) <- c("Parameter", "Estimate", "Posterior SD", "IJ SE")
table_out(d,
 "The first clustered brms dataset: posterior expectations, posterior SDs and clustered IJ SEs.")

## ----transformed-example-code, eval=FALSE-------------------
# prediction <- ijse(
#   fit, cluster_var = fit$data$group,
#   transform = function(d) {
#     cbind(median_at_x1 = d[, "b_Intercept"] + d[, "b_x"])
#   }, diagnostics = "posterior"
# )
# coef(prediction)
# confint(prediction)

## ----worked-transform-table, echo=FALSE, results='asis'-----
d <- read_study("worked-transform.csv")
names(d) <- c("Quantity", "Estimate", "Posterior SD", "IJ SE")
table_out(d,
 "Posterior expectation of the median predictor at the fixed reference value x=1.")

## ----reproduce-commands, eval=FALSE-------------------------
# # From the source directory, after installing IJSE 0.3.0:
# dir.create("study-output", showWarnings = FALSE)
# dir.create("model-cache", showWarnings = FALSE)
# system2("Rscript", c(
#   "tools/run-study.R", ".", "study-output", "4000", "80"))
# system2("Rscript", c(
#   "tools/run-quantile-study.R", ".", "study-output",
#   "model-cache", "200"))
# system2("Rscript", c(
#   "tools/summarise-study.R", ".", "study-output", "figures"))