--- 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.