---
title: "Diagnostics and benchmarking"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Diagnostics and benchmarking}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 7, fig.height = 4.5)
has_lpsolve <- requireNamespace("lpSolve", quietly = TRUE)
knitr::opts_chunk$set(eval = has_lpsolve)
```

This vignette shows the diagnostic and benchmarking workflow for `crr()`
fits: the structured diagnostics report, the plotting methods, and the
one-call benchmark against the unconstrained baseline. It mirrors, at a
CRAN-friendly scale, the simulation comparisons in *Statistical Modeling of
Combinatorial Response Data*.

## Setup: simulate and fit

```{r}
library(combreg)

A <- rbind(
  c(1, 1, 0),
  c(0, 1, 1)
)
con <- crr_constraints(A, b = c(1, 1))
sim <- simulate_crr(n = 80, p = 2, constraints = con, seed = 42)

fit <- crr(sim$Y, sim$X, con,
           n_iter = 800, warmup = 400, chains = 2, seed = 42)
fit
```

## The diagnostics report

`crr_diagnostics()` assembles everything needed to judge a run in one
object: per-coefficient posterior summaries, effective sample sizes (total
and per second of sampling), split-Rhat, MH acceptance, in-sample regression
fit, and posterior predictive goodness-of-fit checks. In simulation the true
coefficients are known, so passing `beta` additionally reports estimation
error (RMSE) and posterior-interval coverage. Potential problems (Rhat above
1.05, ESS below 100, very low acceptance, extreme predictive p-values) are
collected as warnings.

```{r}
report <- crr_diagnostics(fit, beta = sim$beta)
report
```

The underlying per-parameter table is available as `report$table`, and the
predictive-check object as `report$ppc`.

### Reading the goodness-of-fit checks

The raw exact-match and accuracy rates are *not* by themselves evidence for
or against the model: the response is the maximizer of a noisy utility, so
even the true $\beta$ cannot predict $y_i$ perfectly — point-prediction
accuracy is bounded by the intrinsic noise. The posterior predictive checks
calibrate this. Each replicate re-simulates responses from the fitted model
and recomputes the fit statistics; the model is adequate when the observed
statistics are typical of that distribution (two-sided p-values away from
0). Here the observed exact-match rate sits in the middle of its predictive
distribution — the seemingly modest rate is exactly what a *correct* model
predicts for this noise level, and the RMSE and coverage against the known
truth confirm accurate estimation. `crr_ppc()` runs the checks standalone
with more replicates if needed.

### Standalone ESS and split-Rhat

The two MCMC diagnostics that feed the report table are also exported
directly. `crr_ess()` returns the effective sample size per coefficient —
the number of independent draws the autocorrelated chain is worth. As a
rough rule, 100–400 effective draws per quantity of interest is enough for
stable posterior means and intervals; below ~100 the Monte Carlo error is no
longer negligible (this is the threshold the report warns on).

```{r}
crr_ess(fit)
```

`crr_rhat()` returns the split-Rhat potential scale reduction factor, which
compares within- to between-half-chain variance. Values at 1 indicate the
chains have mixed to a common distribution; the usual cut-offs are 1.01
(strict) to 1.05 (lenient), above which the run has not converged and should
be lengthened.

```{r}
crr_rhat(fit)
```

Here every coefficient clears both bars, consistent with the empty warnings
block in the report above.

### Standalone posterior predictive checks

`crr_ppc()` runs the goodness-of-fit checks on their own, and `n_rep`
controls the number of replicates — more replicates sharpen the predictive
distribution and stabilise the p-values at the cost of extra integer-program
solves per replicate. The report uses 50 by default; a hundred or so is a
reasonable standalone budget:

```{r}
ppc <- crr_ppc(fit, n_rep = 100, seed = 1)
ppc
```

The p-values are two-sided: each asks how extreme the observed statistic is
within its own predictive distribution, so values far from 0 mean the
observed data are *typical* of the fitted model. This is the formal version
of the point made above — a modest raw exact-match rate is not evidence
against the model when it sits comfortably inside the band the model itself
predicts, whereas a p-value near 0 would flag genuine misfit. The observed
statistic and its predictive mean and 5–95% band are in `ppc$observed`,
`ppc$replicated`, and `ppc$p_value` for further inspection.

### Reading the report table

The per-parameter table behind the printed summary is a plain data frame in
`report$table`, with columns `mean`, `sd`, the interval bounds, `ess`,
`ess_per_sec`, and `rhat`. Pulling it out lets you find the parameters that
gate the run — sort by ascending ESS to see the slowest-mixing coordinates,
or by descending Rhat for the least-converged:

```{r}
tab <- report$table
head(tab[order(tab$ess), c("parameter", "ess", "ess_per_sec", "rhat")])
head(tab[order(-tab$rhat), c("parameter", "rhat", "ess")])
```

When a report carries warnings (as in the harder matching example below),
this is the practical way to see *which* coefficients are the bottleneck and
by how much, rather than reading the truncated printout. The lightweight
accessors `summary(fit)`, `coef(fit)`, `fitted(fit)`, and `residuals(fit)`
return the posterior summary data frame, the posterior-mean coefficient
matrix, and the fitted / residual response matrices for the same kind of
direct manipulation.

## Plotting methods

All plots operate on the post-warmup, thinned draws.

Trace and autocorrelation:

```{r}
plot(fit, type = "trace", pars = c("beta[1,1]", "beta[2,1]"))
plot(fit, type = "acf", pars = c("beta[1,1]", "beta[2,1]"))
```

Marginal posteriors of all coefficients as violins (median dot, central 95%
interval bar):

```{r}
plot(fit, type = "violin")
```

Effective sample sizes and sampling efficiency (ESS per second of sampling
time, the efficiency measure used for method comparison in the paper):

```{r}
plot(fit, type = "ess")
plot(fit, type = "ess_time")
```

Residual heat map of the posterior point-estimate fit
(`Y - fitted(fit)`, entries in $\{-1, 0, 1\}$):

```{r, fig.height = 5.5}
plot(fit, type = "residual")
```

When the truth is known, the coefficient estimation error
$\hat\beta - \beta$ as a heat map:

```{r}
plot(fit, type = "beta_diff", beta = sim$beta)
```

For anything beyond these built-ins, convert the draws with
`coda::as.mcmc(fit)` or `posterior::as_draws(fit)` and use those ecosystems
directly.

## Benchmarking against the unconstrained baseline

`crr_benchmark()` fits several methods to the same data with identical
settings and reports timing, efficiency, in-sample fit, and — when the true
coefficients are known, as in simulation — RMSE and posterior-interval
coverage:

```{r}
bm <- crr_benchmark(sim$Y, sim$X, con, beta = sim$beta,
                    n_iter = 500, warmup = 250, seed = 42)
bm
```

Two effects visible here reproduce the paper's findings: the constrained
sampler estimates $\beta$ with smaller error and better-calibrated
intervals, because the unconstrained probit baseline conditions on each
coordinate independently and is biased when responses are in fact
constrained maximizers.

The fitted models are kept in `bm$fits`, so any of the diagnostics above can
be applied to each method afterwards, e.g.
`plot(bm$fits$unconstrained, type = "violin")`.

## A non-trivial example: bipartite matching

The paper's motivating application is matching data. Take a $3 \times 3$
bipartite matching: $y \in \{0,1\}^9$ indicates which of the nine possible
edges $(l, r)$ are used, and each node may be matched at most once. The
constraint matrix is the node–edge incidence matrix of $K_{3,3}$ — totally
unimodular by construction:

```{r}
L <- 3; R <- 3
d <- L * R
A <- matrix(0, L + R, d)
for (l in seq_len(L)) {
  for (r in seq_len(R)) {
    edge <- (l - 1) * R + r
    A[l, edge] <- 1        # left node l matched at most once
    A[L + r, edge] <- 1    # right node r matched at most once
  }
}
con_match <- crr_constraints(A, b = rep(1, L + R))
con_match
```

With $d = 9$ coordinates, $m = 6$ constraints, and $p = 3$ covariates the
posterior has 27 coefficients and the feasible set is the set of partial
matchings — a substantially harder target than the toy problem above:

```{r}
sim_m <- simulate_crr(n = 80, p = 3, constraints = con_match, seed = 7)
fit_m <- crr(sim_m$Y, sim_m$X, con_match,
             n_iter = 700, warmup = 350, chains = 2, seed = 7)
report_m <- crr_diagnostics(fit_m, beta = sim_m$beta)
print(report_m, max_rows = 8)
```

This report is doing its job: at this chain length the sampler has **not**
converged — split-Rhat well above 1.05 and effective sample sizes in the
tens are flagged, the MH acceptance rate is much lower than in the small
problem, and the coefficient RMSE and interval coverage against the known
truth are correspondingly degraded. Note the contrast with the posterior
predictive checks, which still pass: the *model* is correct for these data
(the package's test suite verifies the sampler against exact enumerable
posteriors), the *chains* are simply too short. This is the expected
behaviour of the MH-within-Gibbs sampler as $d$ grows — the paper's
production runs use long, heavily thinned chains (e.g. thinning by 25)
precisely because of this autocorrelation, and the ESS-per-second plot is
the tool for budgeting such runs:

```{r}
plot(fit_m, type = "ess_time")
```

The coefficient-error heat map shows the estimation error concentrated in a
few coordinates rather than spread uniformly — the signature of incomplete
mixing rather than model bias:

```{r, fig.height = 4}
plot(fit_m, type = "beta_diff", beta = sim_m$beta)
```

Even with short chains, ignoring the matching constraints costs far more
than the remaining Monte Carlo error. The unconstrained baseline mixes fast
(it is a conjugate Gibbs sampler) but converges to the wrong posterior:

```{r}
bm_m <- crr_benchmark(sim_m$Y, sim_m$X, con_match, beta = sim_m$beta,
                      n_iter = 600, warmup = 300, seed = 7)
bm_m
```

For production use on problems of this size, budget chain length from the
ESS-per-second numbers above and rerun with `n_iter` and `thin` increased
until the report's Rhat and ESS warnings disappear — convergence at this
dimension takes tens of thousands of iterations, which is why the report's
warnings, rather than the raw fit statistics, should drive that decision.
