Introduction to combreg
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):
library(combreg)
A <- rbind(
c(1, 1, 0),
c(0, 1, 1)
)
con <- crr_constraints(A, b = c(1, 1))
con
#> <crr_constraints>: 2 constraints on {0,1}^3 (stored as A y <= b)
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:
is_tum(A) # TRUE: this A is TUM
#> [1] TRUE
is_tum(rbind(c(1, 1, 0), c(1, 0, 1), c(0, 1, 1))) # FALSE (odd cycle)
#> [1] FALSE
#> attr(,"witness")
#> attr(,"witness")$rows
#> [1] 1 2 3
#>
#> attr(,"witness")$cols
#> [1] 1 2 3
#>
#> attr(,"witness")$determinant
#> [1] -2
is_feasible() checks response rows against the polytope
— each returned logical answers “is this row in \(\mathcal{Y}\)?”:
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)
#> [1] TRUE FALSE
The same pattern generalizes to any TUM system. A single “at most one
active coordinate” simplex, for instance, is one dense row:
d <- 3
A_simplex <- matrix(1, nrow = 1, ncol = d)
con_simplex <- crr_constraints(A_simplex, b = 1)
con_simplex
#> <crr_constraints>: 1 constraints on {0,1}^3 (stored as A y <= b)
is_feasible(con_simplex, rbind(diag(d)[1, ], rep(1, d))) # e_1 ok, all-ones not
#> [1] TRUE FALSE
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):
sim <- simulate_crr(n = 60, p = 2, constraints = con, seed = 1)
head(sim$Y)
#> [,1] [,2] [,3]
#> [1,] 1 0 1
#> [2,] 0 0 0
#> [3,] 1 0 1
#> [4,] 0 0 1
#> [5,] 0 0 0
#> [6,] 1 0 1
fit <- crr(sim$Y, sim$X, con,
kernel = "exponential",
n_iter = 400, warmup = 200, seed = 1)
fit
#> <crr_fit> method=mhwg kernel=exponential: n=60, p=2, d=3, chains=1, iterations=400 (warmup 200, thin 1)
#> Mean MH acceptance rate: 0.618
#> Adaptive zeta_block (target accept 0.6): 3
#> Total sampling time: 0.438s
Posterior summaries and diagnostics:
head(summary(fit))
#> parameter mean sd q2.5 q50 q97.5
#> 1 beta[1,1] -0.521653211 0.2206710 -0.9754898 -0.50747053 -0.1181152
#> 2 beta[2,1] 1.132347440 0.2842555 0.6341381 1.11014231 1.7030841
#> 3 beta[1,2] -0.007552769 0.2329463 -0.4749901 0.01156193 0.3876347
#> 4 beta[2,2] -0.060822241 0.2302698 -0.5022534 -0.06933810 0.4615147
#> 5 beta[1,3] -0.169853105 0.2283719 -0.6449864 -0.16139128 0.2525212
#> 6 beta[2,3] 0.817046029 0.2268734 0.4568968 0.79297514 1.3140923
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:
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))
#> constrained unconstrained
#> 0.1398470 0.1692571
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:
fit$zeta_block_tuned
#> [1] 3
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):
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
#> [1] 3
Customizing the sampler
Every knob is set through crr() and
crr_control(). A fully customized fit:
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
#> <crr_fit> method=mhwg kernel=half_gaussian: n=60, p=2, d=3, chains=2, iterations=400 (warmup 200, thin 2)
#> Mean MH acceptance rate: 0.61
#> Adaptive zeta_block (target accept 0.6): 3
#> Total sampling time: 0.424s
Running multiple chains lets you assess convergence with the usual
between-chain diagnostics (see
vignette("diagnostics")).
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.