Package {combreg}

combreg: Bayesian Regression for Combinatorial Response Data

Description

Bayesian regression whose response data are integer-valued vectors subject to combinatorial constraints in the form of Ay<=b. Implements the Metropolis-Hastings-within-Gibbs sampler of Zheng et al. (2026+) doi:10.48550/arXiv.2504.11630, an unconstrained probit baseline for comparison, benchmarking helpers, MCMC and regression diagnostics (effective sample sizes, split-Rhat, structured reports), plotting methods (trace, autocorrelation, violin, effective-sample-size, sampling-efficiency, residual heat map), and utilities for constraint validation (total unimodularity, feasibility) and data simulation.

Author(s)

Maintainer: Hugh Zheng hugh.stats@gmail.com [copyright holder]

See Also

Useful links:

• https://github.com/YuZh98/combreg

• Report bugs at https://github.com/YuZh98/combreg/issues

Convert to coda mcmc.list

Description

Post-warmup, thinned draws per chain.

Usage

## S3 method for class 'crr_fit'
as.mcmc(x, ...)

Arguments

Value

A coda::mcmc.list.

Convert to posterior draws_array

Description

Convert to posterior draws_array

Usage

## S3 method for class 'crr_fit'
as_draws(x, ...)

Arguments

Value

A posterior::draws_array (requires the 'posterior' package).

Posterior mean coefficients

Description

Posterior mean coefficients

Usage

## S3 method for class 'crr_fit'
coef(object, ...)

Arguments

Value

Posterior mean coefficient matrix (p x d).

Examples

if (requireNamespace("lpSolve", quietly = TRUE)) {
  con <- crr_constraints(rbind(c(1, 1, 0), c(0, 1, 1)), b = c(1, 1))
  sim <- simulate_crr(n = 50, p = 2, constraints = con, seed = 1)
  fit <- crr(sim$Y, sim$X, con, n_iter = 200, warmup = 100, seed = 1)
  coef(fit)   # posterior mean of beta, a p x d matrix (row k = covariate k,
              # column j = response coordinate j)
}

Precompute quantities for the conjugate coefficient update

Description

Precompute quantities for the conjugate coefficient update

Usage

coef_precompute(X, prior = crr_prior())

Arguments

Value

List with posterior covariance V and its lower Cholesky factor L, for use in update_coef().

Bayesian combinatorial response regression

Description

Fits the latent-utility regression model

\zeta_i = \beta^\top x_i + \varepsilon_i, \quad \varepsilon_i \sim N(0, I_d), \quad y_i = \arg\max_{y \in \mathcal{Y}} \zeta_i^\top y,

where \mathcal{Y} = \{y \in \{0,1\}^d : Ay \le b\} with totally unimodular A, by MH-within-Gibbs sampling on the dual-certificate augmented posterior (method = "mhwg"). method = "unconstrained" fits independent Albert-Chib probit regressions per coordinate, ignoring the constraints — useful as a baseline demonstrating constraint-ignoring bias.

Usage

crr(
  Y,
  X,
  constraints,
  data = NULL,
  method = c("mhwg", "unconstrained"),
  kernel = c("exponential", "half_gaussian"),
  prior = crr_prior(),
  n_iter = 5000,
  warmup = floor(n_iter/2),
  thin = 1,
  chains = 1,
  seed = NULL,
  control = crr_control(),
  verbose = FALSE
)

Arguments

Value

An object of class crr_fit; see summary.crr_fit().

Examples

if (requireNamespace("lpSolve", quietly = TRUE)) {
  con <- crr_constraints(rbind(c(1, 1, 0), c(0, 1, 1)), b = c(1, 1))
  sim <- simulate_crr(n = 50, p = 2, constraints = con, seed = 1)
  fit <- crr(sim$Y, sim$X, con, n_iter = 200, warmup = 100,
             chains = 1, seed = 1)
  summary(fit)
}

Benchmark samplers on a common data set

Description

Fits several methods to the same data and collects a comparison table: wall-clock time, MH acceptance, effective sample sizes and sampling efficiency (min ESS per second), in-sample fit, and — when the true coefficients are supplied, e.g. from simulate_crr() — estimation error (RMSE) and coverage of the central posterior intervals. This is the package's canonical way to compare the paper's constrained sampler with the unconstrained probit baseline (or any subset of registered methods) without writing bespoke benchmarking code.

Usage

crr_benchmark(
  Y,
  X,
  constraints,
  methods = c("mhwg", "unconstrained"),
  beta = NULL,
  prob = 0.95,
  ...
)

Arguments

Value

An object of class crr_benchmark: a list with table (one row per method) and fits (the underlying crr_fit objects, named by method).

Examples

if (requireNamespace("lpSolve", quietly = TRUE)) {
  con <- crr_constraints(rbind(c(1, 1, 0), c(0, 1, 1)), b = c(1, 1))
  sim <- simulate_crr(n = 50, p = 2, constraints = con, seed = 1)
  bm <- crr_benchmark(sim$Y, sim$X, con, beta = sim$beta,
                      n_iter = 200, warmup = 100, seed = 1)
  bm
}

Constraint system for combinatorial responses

Description

Defines the feasible set \{y \in \{0,1\}^d : Ay \le b\} of a combinatorial response. The matrix A must be integer-valued and, for the exactness of the dual-certificate augmentation used by crr(), totally unimodular.

Usage

crr_constraints(A, b, direction = "<=", check_tum = TRUE, dedup = TRUE)

Arguments

Value

An object of class crr_constraints with elements A, b (in <= form), m, d.

Examples

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

Sampler control parameters

Description

Tuning parameters for the MH-within-Gibbs sampler in crr().

Usage

crr_control(
  n_iter_hit_and_run = 50,
  rho = 1,
  zeta_block = "adaptive",
  zeta_target_accept = 0.6,
  max_dir_tries = 1000,
  bound_truncation = 1e+05,
  n_threads = 1
)

Arguments

Value

An object of class crr_control.

MCMC and regression diagnostics report

Description

Assembles a structured diagnostics report for a fitted model: per-parameter posterior summaries with effective sample sizes, sampling efficiency (ESS per second), split-Rhat convergence diagnostics, MH acceptance rates, in-sample regression fit (exact-match and per-coordinate accuracy of the fitted responses), and posterior predictive goodness-of-fit checks that calibrate the fit statistics against what the model itself predicts (see crr_ppc()). When the true coefficients are known (simulation), supplying beta adds estimation error (RMSE) and interval coverage. Potential problems (low ESS, high Rhat, low acceptance, extreme predictive p-values) are collected as warnings and flagged by print().

Usage

crr_diagnostics(object, beta = NULL, prob = 0.95, n_rep = 50)

Arguments

Value

An object of class crr_diagnostics: a list with elements table (per-parameter data frame with columns mean, sd, interval bounds, ess, ess_per_sec, rhat), fit_stats (exact-match rate and per-coordinate accuracy, or NULL if fitted responses are unavailable), ppc (a crr_ppc object, or NULL if skipped/unavailable), truth (RMSE, coverage, and the error matrix coef(object) - beta, or NULL if beta was not supplied), warnings (character vector), and the sampler configuration.

Examples

if (requireNamespace("lpSolve", quietly = TRUE)) {
  con <- crr_constraints(rbind(c(1, 1, 0), c(0, 1, 1)), b = c(1, 1))
  sim <- simulate_crr(n = 50, p = 2, constraints = con, seed = 1)
  fit <- crr(sim$Y, sim$X, con, n_iter = 200, warmup = 100, seed = 1)
  crr_diagnostics(fit, beta = sim$beta, n_rep = 20)
}

Effective sample sizes

Description

Effective sample size per coefficient, computed from post-warmup thinned draws and summed across chains (via coda::effectiveSize()).

Usage

crr_ess(object)

Arguments

Value

Named numeric vector of length p * d.

Examples

if (requireNamespace("lpSolve", quietly = TRUE)) {
  con <- crr_constraints(rbind(c(1, 1, 0), c(0, 1, 1)), b = c(1, 1))
  sim <- simulate_crr(n = 50, p = 2, constraints = con, seed = 1)
  fit <- crr(sim$Y, sim$X, con, n_iter = 200, warmup = 100, seed = 1)
  crr_ess(fit)   # effective sample size per coefficient; larger means better
                 # mixing (a few hundred is usually plenty for point estimates)
}

Posterior predictive goodness-of-fit checks

Description

Simulates replicated response matrices from the posterior predictive distribution (draw \beta from the retained draws, draw \zeta = X\beta + \varepsilon, map each row to its feasible maximizer) and compares observed test statistics against their predictive distribution. Statistics: the exact-match rate and per-coordinate accuracy of the point predictor fitted(object) applied to each replicate, and the per-coordinate marginal response frequencies.

Usage

crr_ppc(object, n_rep = 50, seed = NULL)

Arguments

Details

Point-prediction accuracy is bounded by the intrinsic utility noise, so a low raw exact-match rate does not by itself indicate misfit: the model fits well when the observed statistics are typical of the predictive distribution (two-sided p-values away from 0). Constrained fits require the 'lpSolve' package.

Value

An object of class crr_ppc: a list with observed (named statistic vector), replicated (n_rep x n_stats matrix), p_value (two-sided predictive p-values per statistic), and n_rep.

Examples

if (requireNamespace("lpSolve", quietly = TRUE)) {
  con <- crr_constraints(rbind(c(1, 1, 0), c(0, 1, 1)), b = c(1, 1))
  sim <- simulate_crr(n = 50, p = 2, constraints = con, seed = 1)
  fit <- crr(sim$Y, sim$X, con, n_iter = 200, warmup = 100, seed = 1)
  crr_ppc(fit, n_rep = 20, seed = 1)
}

Prior specification for crr models

Description

Currently a single family: independent Gaussian priors \beta_{kj} \sim N(0, sd^2) on all regression coefficients.

Usage

crr_prior(family = "gaussian", sd = 1)

Arguments

Value

An object of class crr_prior.

Examples

crr_prior()               # N(0, 1)
crr_prior(sd = 10)        # weakly informative

Split-Rhat convergence diagnostics

Description

Potential scale reduction factor per coefficient, computed on split half-chains so it is defined even for chains = 1.

Usage

crr_rhat(object)

Arguments

Value

Named numeric vector of length p * d.

Examples

if (requireNamespace("lpSolve", quietly = TRUE)) {
  con <- crr_constraints(rbind(c(1, 1, 0), c(0, 1, 1)), b = c(1, 1))
  sim <- simulate_crr(n = 50, p = 2, constraints = con, seed = 1)
  fit <- crr(sim$Y, sim$X, con, n_iter = 200, warmup = 100, seed = 1)
  crr_rhat(fit)  # values near 1 indicate convergence; above ~1.05 is a
                 # sign the chain needs to run longer
}

Draw latent utilities from their truncated-normal full conditional

Description

For each coordinate j in subset, draws \zeta_{\cdot j} \sim TN(\mu_{\cdot j}, 1) truncated below at (UA)_{\cdot j} where y = 1 and above where y = 0.

Usage

draw_utility(zeta, Mu, Y, UA, subset = seq_len(ncol(zeta)))

Arguments

Value

Matrix (n x d).

Check dual-certificate feasibility

Description

Verifies, per observation, that a dual matrix lies in \mathcal{D}(y_i, \zeta_i) \cap \{u \ge 0\}.

Usage

dual_feasible(constraints, zeta, Y, U, control = crr_control())

Arguments

Value

Integer 0/1 vector of length n.

Fitted responses

Description

Posterior-point-estimate responses on the training covariates: predict(object, type = "response").

Usage

## S3 method for class 'crr_fit'
fitted(object, ...)

Arguments

Value

Binary matrix (n x d).

Initialize dual certificates

Description

Computes the active-constraint mask active (1\{A y_i = b\}) and an initial dual matrix U with independent Exp(1) entries on the active support.

Usage

init_dual(constraints, Y)

Arguments

Value

List with U (n x m) and active (n x m binary mask).

Check feasibility of responses

Description

Check feasibility of responses

Usage

is_feasible(constraints, Y)

Arguments

Value

Logical vector of length n: does each row satisfy ⁠A y <= b⁠ and y in {0,1}?

Check total unimodularity

Description

Exhaustively checks whether every square submatrix of A has determinant in {-1, 0, 1}. Cost grows exponentially in min(nrow(A), ncol(A)); use only for small-to-moderate matrices.

Usage

is_tum(A)

Arguments

Value

TRUE or FALSE, with attribute "witness" (a list with rows, cols, determinant) when FALSE.

Examples

is_tum(rbind(c(1, 1, 0), c(0, 1, 1)))       # TRUE
is_tum(rbind(c(1, 1, 0), c(1, 0, 1), c(0, 1, 1)))  # FALSE

Diagnostic plots for a crr fit

Description

All plots use post-warmup thinned draws.

Usage

## S3 method for class 'crr_fit'
plot(
  x,
  pars = NULL,
  type = c("trace", "acf", "violin", "ess", "ess_time", "residual", "beta_diff"),
  beta = NULL,
  ...
)

Arguments

Details

• "trace": per-parameter trace plots, chains overlaid.

• "acf": per-parameter autocorrelation functions (chains pooled).

• "violin": violin plot of the marginal posterior of each coefficient, with median and central 95% interval.

• "ess": effective sample size per coefficient, with a dashed reference line at 100.

• "ess_time": effective sample size per second of sampling time — the efficiency measure used for method comparison in the paper.

• "residual": heat map of the training residuals Y - fitted(x) (entries in {-1, 0, 1}); rows are observations, columns response coordinates.

• "beta_diff": heat map of the coefficient estimation error coef(x) - beta against a known true coefficient matrix (supplied via beta, e.g. from simulate_crr()), with the error printed in each cell.

Value

x, invisibly.

Posterior predictions

Description

Posterior predictions

Usage

## S3 method for class 'crr_fit'
predict(object, newdata = NULL, type = c("utility", "response"), ...)

Arguments

Value

Matrix (n_new x d) of utilities or feasible responses.

Examples

if (requireNamespace("lpSolve", quietly = TRUE)) {
  con <- crr_constraints(rbind(c(1, 1, 0), c(0, 1, 1)), b = c(1, 1))
  sim <- simulate_crr(n = 50, p = 2, constraints = con, seed = 1)
  fit <- crr(sim$Y, sim$X, con, n_iter = 200, warmup = 100, seed = 1)
  head(predict(fit, type = "response"))  # feasible 0/1 predictions; each row
                                         # satisfies the constraints A y <= b
}

Generate a random totally unimodular constraint system

Description

For small problems (d * m <= 50 and m >= 2 fails), draws random {-1, 0, 1} matrices until one is totally unimodular, with b = 1. For larger problems, generates a random network (incidence-type) matrix with one +1 and one -1 per row — totally unimodular by construction — and Bernoulli right-hand sides. This mirrors the simulation design of the paper.

Usage

random_constraints(d, m)

Arguments

Value

A crr_constraints object.

Response residuals

Description

Training residuals Y - fitted(object), entries in {-1, 0, 1}.

Usage

## S3 method for class 'crr_fit'
residuals(object, ...)

Arguments

Value

Integer matrix (n x d).

Update dual certificates given latent utilities

Description

One Gibbs update of U \mid \zeta, y: for each observation, runs hit-and-run over the dual-certificate polytope \mathcal{D}(y_i, \zeta_i) \cap \{u \ge 0\} targeting the chosen kernel, and returns the endpoint.

Usage

sample_dual(
  constraints,
  zeta,
  Y,
  U,
  active,
  kernel = "exponential",
  control = crr_control()
)

Arguments

Value

Updated dual matrix (n x m).

Metropolis-Hastings update of the latent utilities

Description

One MH step of \zeta \mid U, \beta, y: proposes new utilities for a random block of coordinates from the truncated-normal proposal, forms the coordinatewise envelope \tilde\zeta (max where y=1, min where y=0), runs hit-and-run under \tilde\zeta to obtain a dual certificate U^\star, and accepts each observation's proposal iff its u_i^\star is feasible for the current \zeta_i.

Usage

sample_utility(
  constraints,
  zeta,
  Y,
  U,
  active,
  Mu,
  kernel = "exponential",
  control = crr_control(),
  block = NULL
)

Arguments

Value

List with zeta (updated matrix) and accept (0/1 vector, per-observation acceptance).

Simulate combinatorial response regression data

Description

Draws X_{ik} \sim N(0,1), coefficients \beta_{kj} \sim N(0,1) (unless supplied), latent utilities \zeta_i = \beta^\top x_i + \varepsilon_i with standard normal noise, and responses y_i = \arg\max_{y \in \mathcal{Y}} \zeta_i^\top y solved by integer programming.

Usage

simulate_crr(n, p, constraints, beta = NULL, seed = NULL)

Arguments

Value

List with Y (n x d), X (n x p), beta (p x d, the truth), zeta (n x d), and constraints.

Examples

if (requireNamespace("lpSolve", quietly = TRUE)) {
  con <- crr_constraints(rbind(c(1, 1, 0), c(0, 1, 1)), b = c(1, 1))
  sim <- simulate_crr(n = 20, p = 2, constraints = con, seed = 1)
}

Posterior summary of a crr fit

Description

Posterior summary of a crr fit

Usage

## S3 method for class 'crr_fit'
summary(object, ...)

Arguments

Value

A data frame with posterior mean, sd, and 2.5/50/97.5 percent quantiles per coefficient, computed from post-warmup thinned draws pooled across chains.

Conjugate Gaussian update of the regression coefficients

Description

Draws \beta_{\cdot j} \sim N(V X^\top \zeta_{\cdot j}, V) jointly for all response coordinates.

Usage

update_coef(X, zeta, precomp)

Arguments

Value

Coefficient matrix (p x d).