---
title: "Tutorial 1: Getting Started with IJSE"
output:
  rmarkdown::html_vignette:
    css: ijse.css
vignette: >
  %\VignetteIndexEntry{Tutorial 1: Getting Started with IJSE}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r quick-setup, include=FALSE}
library(IJSE)
options(width = 62)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

## What you will learn

This tutorial runs without brms, Stan or a C++ compiler. You will construct
matched posterior draws and log likelihoods, distinguish posterior uncertainty
from sampling uncertainty, calculate a nonlinear quantity, and verify a streamed
calculation. All observations are simulated.

Install the supplied source release after installing its required dependency:

```{r quick-install, eval=FALSE, purl=FALSE}
install.packages("posterior")
install.packages("IJSE_0.3.3.tar.gz",
                 repos = NULL, type = "source")
library(IJSE)
```

The package supplies three interfaces. `ijse_draws()` accepts matched draws and
log likelihoods from any backend. `ijse()` extracts them from a brms fit.
`IJ_se()` preserves the original brms interface and returns only a named vector
of IJ standard errors. Start with `ijse_draws()` to see what the calculation uses.

## Construct an example with a known posterior

Generate 80 independent observations with mean 0.5 and SD 2. Fit a deliberately
narrow working likelihood, $y_i\mid\mu\sim N(\mu,1)$, and prior
$\mu\sim N(0,4)$. The posterior has variance $v=(N+1/4)^{-1}$ and mean
$m=v\sum_i y_i$. Sampling it directly separates this first lesson from MCMC.

```{r quick-data}
set.seed(431)
y <- rnorm(80, mean = 0.5, sd = 2)
n <- length(y)
v <- 1 / (n + 1 / 4)
m <- v * sum(y)
mu <- rnorm(20000, mean = m, sd = sqrt(v))
draws <- cbind(mu = mu)
rownames(draws) <- paste0("s", seq_along(mu))
ll <- vapply(y, function(yi) dnorm(yi, mu, log = TRUE),
             numeric(length(mu)))
rownames(ll) <- rownames(draws)
c(draw_rows = nrow(draws), likelihood_rows = nrow(ll),
  observations = ncol(ll))
```

Rows index posterior draws; likelihood columns index observations. Column $i$
contains $\log p(y_i\mid\mu^{(s)})$ for every retained draw $s$. The two inputs
must use exactly the same draw rows. Row names provide a useful check, but do
not prove that someone assigned correct labels to the data.

Use the likelihood, excluding the prior. A log posterior or one joint log
likelihood column cannot replace the required observation contributions.

```{r quick-fit}
result <- ijse_draws(draws, ll, diagnostics = "none")
data.frame(Estimate = unname(coef(result)),
  Posterior_SD = unname(result$posterior_sd),
  IJ_SE = unname(result$se))
```

`Estimate` is the posterior expectation, approximated by the draw mean.
`Posterior_SD` is the spread of the posterior distribution. `IJ_SE` estimates
how the posterior expectation varies across repeated samples of the observations.
In this example the working variance is too small, so the two uncertainty
measures have a substantive reason to differ. A larger IJ SE is a result to
interpret, rather than an automatic sign that a model was repaired.

## Read intervals and covariance correctly

```{r quick-intervals}
confint(result, level = 0.95)
vcov(result)
vcov(result, type = "posterior")
```

`confint()` returns a normal interval centered on the posterior expectation,
using the IJ SE. It is a sampling-uncertainty interval, rather than an equal-tailed
posterior credible interval. `vcov(result)` contains IJ sampling covariance;
`vcov(result, type = "posterior")` contains covariance of the posterior draws.
Both are matrices even with one quantity. The complete `summary(result)` also
reports the SE ratio and intervals.

An IJ interval can be poorly calibrated if the estimator is biased, if the
sampling units are wrong, or if the variance approximation is inaccurate.
Changing a standard error does not remove those problems by itself.

## Check against a calculation you know

For this normal model, the exact IJ SE under the package's finite-sample
convention is $v\sqrt{N s_y^2}$.

```{r quick-oracle}
exact_se <- v * sqrt(n * var(y))
data.frame(Analytic_IJ = exact_se,
           Draw_based_IJ = unname(result$se))
```

The difference here comes from a finite number of posterior draws. There is
also an exact identity for this particular model: estimated IJ SE divided by
exact IJ SE equals the draws' sample variance divided by $v$. The following
check tests numerical implementation, without assuming that a random draw
estimate equals its population target exactly.

```{r quick-identity}
stopifnot(isTRUE(all.equal(
  unname(result$se) / exact_se, var(mu) / v,
  tolerance = 1e-10)))
```

That identity is specific to this example. For MCMC, inspect the selected
quantities and the stability of the covariance calculation; a large draw count
alone does not establish precise IJ inference.

## Transform every draw

To estimate $E(\mu^2\mid y)$, square each draw and then average. Squaring the
posterior mean instead gives a different target. Select input parameters before
applying a transformation, and return one row per input draw.

```{r quick-transform}
nonlinear <- ijse_draws(draws, ll,
  transform = function(d) {
    cbind(mu = d[, "mu"], mu_squared = d[, "mu"]^2)
  }, diagnostics = "none")
coef(nonlinear)
c(Exact_E_mu_squared = m^2 + v, Squared_E_mu = m^2)
vcov(nonlinear)
```

The off-diagonal covariance matters for joint contrasts. For example, obtain the
SE of $E(\mu^2-2\mu\mid y)$ either directly from transformed draws or by the
linear combination of the full IJ covariance matrix.

```{r quick-contrast}
contrast <- ijse_draws(draws, ll,
  transform = function(d) d[, "mu"]^2 - 2 * d[, "mu"],
  diagnostics = "none")
a <- c(mu = -2, mu_squared = 1)
contrast_se <- sqrt(drop(t(a) %*% vcov(nonlinear) %*% a))
stopifnot(isTRUE(all.equal(
  unname(contrast$se), contrast_se)))
```

The same principle applies to products of jointly sampled parameters. Pairing
rows from separate fitted models does not generally create a joint posterior.

## Evaluate likelihood blocks

A callback accepts observation indices and returns a draws-by-observations
matrix for those indices, in the requested order. It must preserve a matrix
when there is only one likelihood column.

```{r quick-stream}
read_likelihood <- function(ii) {
  block <- vapply(y[ii],
    function(yi) dnorm(yi, mu, log = TRUE),
    numeric(length(mu)))
  rownames(block) <- rownames(draws)
  block
}
streamed <- ijse_draws(draws, read_likelihood,
  n_obs = n, block_size = 7, diagnostics = "none")
stopifnot(isTRUE(all.equal(vcov(streamed), vcov(result))))
```

This tutorial keeps the dense matrix to verify equality. In production, a
callback can generate blocks without ever allocating that full matrix. Blocking
limits likelihood storage; its runtime advantage depends on the backend. It is
not a general speed guarantee. If subsetting a stored matrix in your callback,
use `ll[, ii, drop = FALSE]`.

## Keep a useful result record

```{r quick-record}
record <- list(
  estimate = coef(result),
  ij_covariance = vcov(result),
  posterior_covariance = vcov(result, type = "posterior"),
  interval = confint(result),
  draws = result$n_draws, observations = result$n_obs,
  units = result$n_units, clustered = result$clustered,
  version = as.character(packageVersion("IJSE")))
stopifnot(record$units == n, !record$clustered)
```

When reporting an analysis, state the posterior expectation, working likelihood,
sampling units, interval reference, and posterior simulation settings. Keep IDs,
transformation code and row-selection decisions with the result. Continue with
`vignette("ijse-brms", package = "IJSE")` for fitted models and clustering, or
`vignette("ijse-diagnostics", package = "IJSE")` for diagnostics and error cases.