---
title: "Tutorial 2: brms, Clusters and Joint Contrasts"
output:
  rmarkdown::html_vignette:
    css: ijse.css
vignette: >
  %\VignetteIndexEntry{Tutorial 2: brms, Clusters and Joint Contrasts}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r brms-setup, include=FALSE}
library(IJSE)
stopifnot(requireNamespace("brms", quietly = TRUE))
options(width = 62)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
se_table <- function(x) {
  round(data.frame(Estimate = coef(x),
    Posterior_SD = x$posterior_sd, IJ_SE = x$se), 3)
}
```

## What you will learn

This tutorial uses two installed, synthetic brms fits. You can run all evaluated
code without compiling Stan or downloading data. Install brms separately if
needed. The optional refit at the end requires cmdstanr, CmdStan and a compiler.

Both examples fit a working asymmetric-Laplace median regression,
$y\sim 1+x$, with identity link, estimated scale, and fixed quantile 0.5. Data
have generating conditional median $0.5+0.7x$. The independent example has
heteroscedastic normal errors. The clustered example has 30 groups of ten, with
within-group correlation 0.6 in the generated errors. The fitted likelihood uses
working independence and no random effects. These are worked examples, rather
than an additional coverage study.

## Load a fit together with its sampling IDs

```{r brms-load}
path <- system.file("extdata", "tutorial",
  "median-examples.rds", package = "IJSE")
stopifnot(nzchar(path))
examples <- readRDS(path)
ind <- examples$independent
cl <- examples$clustered
stopifnot(inherits(ind$fit, "brmsfit"),
          inherits(cl$fit, "brmsfit"))
```

The bundle keeps `fit`, original fitting `data`, and `cluster` separately.
brms can discard a variable unused by its model formula, so a missing
`fit$data$group` must not silently become your independent-observation default.
Use the saved IDs and verify the fitting rows.

```{r brms-check-rows}
fit <- cl$fit
ids <- cl$cluster
stopifnot(length(ids) == nrow(cl$data),
  identical(ids, cl$data$group),
  isTRUE(all.equal(fit$data$x, cl$data$x)),
  isTRUE(all.equal(fit$data$y, cl$data$y)),
  !anyNA(ids), length(unique(ids)) == 30L)
```

For your own univariate regression, select complete fitting rows explicitly
before fitting and save their original row IDs and cluster IDs. A length check
alone cannot establish substantive row order. The checks above apply to this
ordinary regression; likelihoods with special observation orderings require
separate care.

## Begin with independent observations

```{r brms-independent}
independent <- ijse(ind$fit, diagnostics = "posterior")
se_table(independent)
confint(independent)
stopifnot(!independent$clustered,
          independent$n_units == 200L)
stopifnot(isTRUE(all.equal(
  IJ_se(ind$fit), independent$se)))
```

By default, `ijse()` selects parameters whose names start with `b_`.
`coef()` returns posterior expectations. The original `IJ_se()` returns only
the IJ SE vector; its pooled inference is the same and it skips discarded
diagnostics. Select exact parameter names when your transformation needs a
particular subset or a distributional parameter.

## Use clusters as the sampling units

Compare two calculations on the *same* clustered-data fit. This changes the
resampling units, not the posterior draws or point estimates.

```{r brms-clustered}
observation_ij <- ijse(fit, diagnostics = "none")
cluster_ij <- ijse(fit, cluster_var = ids,
  save_influence = TRUE, diagnostics = "all")
se_table(cluster_ij)
data.frame(Observation_IJ = observation_ij$se,
           Cluster_IJ = cluster_ij$se)
stopifnot(cluster_ij$clustered,
  cluster_ij$n_obs == 300L, cluster_ij$n_units == 30L,
  isTRUE(all.equal(coef(observation_ij), coef(cluster_ij))))
```

Cluster aggregation sums likelihood sensitivities within each group before
forming sampling covariance. It treats groups as independent sampling units.
It does not change the working likelihood, fit random effects, or automatically
validate arbitrary hierarchical models. Its normal intervals also do not supply
a universal correction for few or highly unequal clusters.

## Compute predictions and a joint contrast

For this identity-link median model, the predictor at fixed $x=1$ is the
intercept plus the slope. At $x=-1$, it is the intercept minus the slope. Compute
both on every draw, with their difference, so all covariance terms are retained.
These are posterior expectations of fixed-reference predictors; estimating a
reference distribution would add another source of sampling uncertainty.

```{r brms-transform}
predictions <- ijse(fit, cluster_var = ids,
  variables = c("b_Intercept", "b_x"),
  transform = function(d) {
    cbind(at_minus1 = d[, "b_Intercept"] - d[, "b_x"],
          at_plus1 = d[, "b_Intercept"] + d[, "b_x"],
          difference = 2 * d[, "b_x"])
  }, diagnostics = "posterior")
se_table(predictions)
confint(predictions, parm = "difference", level = 0.90)
```

The contrast has coefficient vector $(0,2)$ on the original intercept and slope.
Verify its SE using the original full covariance matrix.

```{r brms-contrast-check}
a <- c(b_Intercept = 0, b_x = 2)
v <- vcov(cluster_ij)[names(a), names(a), drop = FALSE]
contrast_se <- sqrt(drop(t(a) %*% v %*% a))
stopifnot(isTRUE(all.equal(
  unname(predictions$se["difference"]), contrast_se)))
```

Do not add marginal SEs to obtain a contrast SE. Use the covariance matrix or
transform the joint draws. A contrast between predictors is an associational
quantity here; the simulation and model do not make it a causal effect.

## Select draws without losing alignment

`draw_ids` uses post-warmup draw indices, ordered by iteration within chain,
then chain. IJSE applies the same indices to parameter and likelihood extraction.
For this example there are 1,000 retained draws in each of two chains. The
following selects half of each chain for a numerical stability comparison.

```{r brms-subset}
arr <- posterior::as_draws_array(fit,
  variable = c("b_Intercept", "b_x"), regex = FALSE)
n_iter <- posterior::niterations(arr)
n_chain <- posterior::nchains(arr)
keep <- unlist(lapply(seq_len(n_chain), function(j) {
  (j - 1L) * n_iter + seq.int(1L, n_iter, by = 2L)
}))
subset_result <- ijse(fit, cluster_var = ids,
  draw_ids = keep, diagnostics = "posterior")
stopifnot(subset_result$n_draws == length(keep))
data.frame(All_draws = cluster_ij$se,
           Selected_draws = subset_result$se)
```

The estimates can differ because fewer posterior draws estimate the covariance
less precisely. This comparison does not imply thinning is a precision remedy.
Inspect convergence and collect sufficient effective posterior information.

## Verify streamed likelihood extraction

```{r brms-stream}
streamed <- ijse(fit, cluster_var = ids,
  block_size = 40, diagnostics = "none")
stopifnot(isTRUE(all.equal(vcov(streamed),
                         vcov(cluster_ij))))
```

This univariate likelihood supports streaming. The tested brms implementation
does not support IJSE streaming for regression ARMA models, and multivariate
streaming is rejected explicitly. Use `block_size = NULL` in those cases.
Covariance ARMA models require explicit joint series likelihoods through
`ijse_draws()`: their native brms output contains LOO conditional densities.
Discrete truncated models also require explicit likelihood matrices because
backend lower-bound normalization is not validated. Backend prediction objects
can still consume memory when likelihoods are blocked.

For measurement-error predictors, fit with
`save_pars = brms::save_pars(latent = TRUE)` to retain matching latent draws.
Without them, IJSE rejects prediction-time resampling. The native likelihood
still conditions on those latents; an observed-data analysis needs a likelihood
that integrates them. Missing-response models can return NA likelihood columns;
these require a separately derived observed-data likelihood, rather than
dropping incomplete cases or substituting zero for NA entries indiscriminately.

## Review diagnostics and save the analysis

```{r brms-diagnostics}
data.frame(Rhat = cluster_ij$diagnostics$rhat,
  Bulk_ESS = cluster_ij$diagnostics$ess_bulk)
cluster_ij$diagnostics$chain_relative_range
```

These R-hat and ESS values refer to the selected quantities. Check the original
fit's diagnostics too, including varying nuisance parameters and divergent
transitions. Chain-specific IJ ranges are stability checks, not calibrated
Monte Carlo SEs of IJ covariance. The two tutorial fits have their generation
seeds, fit seeds, fitting settings and diagnostic record in the installed data
bundle's companion files.

```{r brms-save, eval=FALSE, purl=FALSE}
saveRDS(list(result = cluster_ij, cluster = ids,
  data = cl$data, transform = predictions$call,
  session = sessionInfo()), "my-ijse-analysis.rds")
```

For reporting, identify the estimand, working likelihood, 30 independent clusters,
normal interval reference, posterior settings and diagnostic limitations.
A single example's interval width does not demonstrate repeated-sampling
coverage. The separate `ijse-study` vignette evaluates that question.

## Optional: fit the model yourself

The following code is displayed but not run while building the tutorial. It
refits this example's data and preserves the IDs separately. Using another
backend, seed or toolchain can change finite-draw results.

```{r brms-refit, eval=FALSE, purl=FALSE}
model_data <- cl$data
model_ids <- model_data$group
refit <- brms::brm(
  brms::bf(y ~ x, quantile = 0.5), data = model_data,
  family = brms::asym_laplace(link_quantile = "identity"),
  prior = c(
    brms::set_prior("normal(0, 2.5)", class = "b"),
    brms::set_prior("normal(0, 2.5)", class = "Intercept"),
    brms::set_prior("exponential(1)", class = "sigma")),
  backend = "cmdstanr", chains = 2, cores = 2,
  iter = 2000, warmup = 1000, seed = 41302,
  control = list(adapt_delta = 0.99, max_treedepth = 12))
refit_ij <- ijse(refit, cluster_var = model_ids)
stopifnot(refit_ij$clustered, refit_ij$n_units == 30L)
```

To regenerate both installed examples, run
`tools/generate-tutorial-data.R` from the source release with an output directory
and model-cache directory. Refit execution is separate from package checks and
from the coverage simulation in the study.