---
title: "Introduction to combreg"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to combreg}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## The model

`combreg` fits Bayesian regression models whose response data are binary vectors constrained to an integral polytope,

$$\mathcal{Y} = \{\, y \in \{0,1\}^d : A y \le b \,\}.$$

Examples include matchings, survey with skip-logic, network flow, and other combinatorial structures. The proposed model considers each observation having a latent Gaussian utility

$$\zeta_i = \beta^\top x_i + \varepsilon_i, \qquad
\varepsilon_i \sim N(0, I_d),$$

and the observed response is the feasible utility maximizer

$$y_i = \arg\max_{y \in \mathcal{Y}} \; \zeta_i^\top y.$$

Posterior sampling uses a Metropolis–Hastings-within-Gibbs scheme on an augmented posterior in which each observation carries a dual variable $u_i \ge 0$ of the optimality of $y_i$; the dual updates use hit-and-run sampling implemented in C++. For more details, see the paper "Statistical Modeling of Combinatorial Response Data". 

## Defining constraints

Constraints are declared once and validated (including an exhaustive total unimodularity check for small systems):

```{r}
library(combreg)

A <- rbind(
  c(1, 1, 0),
  c(0, 1, 1)
)
con <- crr_constraints(A, b = c(1, 1))
con
```

Here `A` and `b` describe the polytope $\{\, y \in \{0,1\}^d : A y \le b \,\}$ directly: each row of `A y <= b` is one linear inequality carving the feasible set out of the binary cube. The two rows above encode "at most one of coordinates 1--2" and "at most one of coordinates 2--3", corresponding to $y_{i1}+y_{i2}\le 1$ and $y_{i2}+y_{i3}\le 1$, respectively. The feasible set is thus $\{(0,0,0), (1,0,0), (0,1,0), (0,0,1),(1,0,1)\}$.

`is_tum()` and `is_feasible()` are available as standalone utilities. The method requires the polytope $\{y\in[0,1]^d: Ay\le b\}$ to be integral, for which the total unimodularity of `A` is a sufficient condition, so it is checked at construction; you can also test a matrix directly:

```{r}
is_tum(A)                                     # TRUE: this A is TUM
is_tum(rbind(c(1, 1, 0), c(1, 0, 1), c(0, 1, 1)))  # FALSE (odd cycle)
```

`is_feasible()` checks response rows against the polytope --- each returned
logical answers "is this row in $\mathcal{Y}$?":

```{r}
candidates <- rbind(
  c(1, 0, 1),   # A y = (1, 1) <= (1, 1): feasible
  c(1, 1, 0)    # A y = (2, 1): violates row 1
)
is_feasible(con, candidates)
```

The same pattern generalizes to any TUM system. A single "at most one active coordinate" simplex, for instance, is one dense row:

```{r}
d <- 3
A_simplex <- matrix(1, nrow = 1, ncol = d)
con_simplex <- crr_constraints(A_simplex, b = 1)
con_simplex
is_feasible(con_simplex, rbind(diag(d)[1, ], rep(1, d)))  # e_1 ok, all-ones not
```

## Simulating data and fitting

`simulate_crr()` draws covariates, coefficients, and utilities, and maps each utility vector to its constrained maximizer by integer programming (this requires the `lpSolve` package):

```{r, eval = has_lpsolve}
sim <- simulate_crr(n = 60, p = 2, constraints = con, seed = 1)
head(sim$Y)
```

```{r, eval = has_lpsolve}
fit <- crr(sim$Y, sim$X, con,
           kernel = "exponential",
           n_iter = 400, warmup = 200, seed = 1)
fit
```

Posterior summaries and diagnostics:

```{r, eval = has_lpsolve}
head(summary(fit))
plot(fit, pars = c("beta[1,1]", "beta[2,1]"))
```

Draws convert directly to the `coda` ecosystem via `coda::as.mcmc(fit)`, and to `posterior::draws_array` via `posterior::as_draws(fit)`. For the full diagnostics workflow, including the structured `crr_diagnostics()` report, violin / ESS / efficiency / residual plots, and benchmarking with `crr_benchmark()`, see `vignette("diagnostics")`.

## Why constraints matter

Ignoring the constraint structure biases estimation. The `"unconstrained"` method fits independent Albert–Chib probit regressions per coordinate, a natural but wrong model for constrained responses:

```{r, eval = has_lpsolve}
fit_un <- crr(sim$Y, sim$X, method = "unconstrained",
              n_iter = 400, warmup = 200, seed = 1)

rmse <- function(est, truth) sqrt(mean((est - truth)^2))
c(constrained   = rmse(coef(fit), sim$beta),
  unconstrained = rmse(coef(fit_un), sim$beta))
```

The constrained sampler recovers $\beta$ markedly better because the likelihood correctly conditions on $y_i$ being the feasible maximizer rather than a vector of independent signs.

## The adaptive `zeta_block` controller

The latent-utility MH update does not refresh all $d$ coordinates at once. Each sweep resamples a random sub-block of the coordinates, and the block size --- `zeta_block` --- controls the trade-off between per-step acceptance and how far the utilities move: small blocks accept often but mix slowly, large blocks propose bolder moves that are rejected more often. The optimal size depends on the constraint geometry, not just on $d$, so a good default cannot be a fixed number.

The default `zeta_block = "adaptive"` handles this automatically. During warmup a dual-averaging controller tunes the block size to hit `zeta_target_accept` (default `0.6`), then freezes it for the sampling phase so the proposal is fixed (valid MCMC). The tuned value is reported by `print()` and stored on the fit; the `fit` above already used it:

```{r, eval = has_lpsolve}
fit$zeta_block_tuned
```

You therefore need not tune this by hand unless you witness a low acceptance rate during MH-Within-Gibbs sampling. You can also manually set `zeta_block` as a positive integer to turn off the adaptive controller. In this case, the block size is the minimum of $d$ and the specified integer. For example, you can set `zeta_block = 100` to cap the block size at 100 (note that since $d=3$ in this example, the block size will be capped at $3$, and $100$ will be ignored):

```{r, eval = has_lpsolve}
fit_fixed <- crr(sim$Y, sim$X, con,
                 n_iter = 300, warmup = 150, seed = 1,
                 control = crr_control(zeta_block = 100))
fit_fixed$zeta_block_tuned   # block sizes are capped at d
```

## Customizing the sampler

Every knob is set through `crr()` and `crr_control()`. A fully customized fit:

```{r, eval = has_lpsolve}
fit_custom <- crr(sim$Y, sim$X, con,
                  kernel = "half_gaussian",         # half-Gaussian dual kernel
                  prior  = crr_prior(sd = 10),      # wider N(0, 100) prior
                  n_iter = 400, warmup = 200,       # total / discarded sweeps
                  thin   = 2,                       # keep every 2nd draw
                  chains = 2,                       # independent chains (serial)
                  seed   = 1,                        # reproducible RNG stream
                  control = crr_control(
                    n_iter_hit_and_run = 20,        # inner hit-and-run steps
                    zeta_block = "adaptive",        # auto-tune the block size
                    n_threads = 1))                 # OpenMP threads (result-invariant)
fit_custom
```

Running multiple chains lets you assess convergence with the usual between-chain diagnostics (see `vignette("diagnostics")`).

## Formula interface

In the function `crr()`, the second argument `X` may be a one-sided formula evaluated against a data frame. Because `simulate_crr()` uses a no-intercept design, drop the intercept with `~ 0 + ...` so the design matrix matches the simulated one:

```{r, eval = has_lpsolve}
df <- data.frame(x1 = sim$X[, 1], x2 = sim$X[, 2])
fit_formula <- crr(sim$Y, ~ 0 + x1 + x2, con, data = df,
                   n_iter = 300, warmup = 150, seed = 1)
coef(fit_formula)
```

The formula path uses `model.matrix()`, so factor expansion and interactions work as in any R modelling function.

## Building custom models from primitives

Models with structured (non-linear-in-$X$) means can reuse the exported primitives (`init_dual()`, `sample_dual()`, `draw_utility()`, `sample_utility()`, `dual_feasible()`, `coef_precompute()`, `update_coef()`) and supply their own coefficient updates, so that no constraint logic needs to be reimplemented. See their reference pages for the exact conditional distributions.
