Open-source · github.com/YuZh98/combreg

combreg

Bayesian regression for combinatorial response data: each response is a binary vector constrained to an integral polytope, $\mathcal{Y} = {,y \in {0,1}^d : A y \le b,}$. Ignoring the constraint and fitting an unconstrained model leads to biased estimation, and enumerating all feasible outcomes is computationally impractical in high-dimensional settings. In Zheng, Ghosh & Duan (2026), an augmented likelihood that respects the combinatorial constraints is proposed, along with a Metropolis–Hastings-within-Gibbs (MH-Within-Gibbs) sampler that scales to high-dimensional problems.

The package implements the MH-Within-Gibbs sampler of Zheng, Ghosh & Duan (2026), with hit-and-run dual updates in C++ (OpenMP-parallel across observations), an unconstrained probit baseline, one-call benchmarking, MCMC and regression diagnostics, and utilities for constraint validation and data simulation.

Installation

Install the released version from CRAN:

1

install.packages("combreg")

Or the development version from GitHub:

1
2

# install.packages("remotes")
remotes::install_github("YuZh98/combreg")

Building from source requires a C++17 toolchain (the package links to RcppArmadillo); OpenMP is used when available and falls back to a serial build otherwise.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

library(combreg)

# a small matching-type constraint: A y <= b on y in {0,1}^3
A <- rbind(c(1, 1, 0),
           c(0, 1, 1))
con <- crr_constraints(A, b = c(1, 1))

# simulate constrained responses, then fit
sim <- simulate_crr(n = 300, p = 2, constraints = con, seed = 1)
fit <- crr(sim$Y, sim$X, con, n_iter = 2000, warmup = 1000, chains = 2, seed = 1)
summary(fit)

# alternatively, predictors can also be given as a formula + data frame
df  <- data.frame(x1 = sim$X[, 1], x2 = sim$X[, 2])
fit <- crr(sim$Y, ~ 0 + x1 + x2, con, data = df,
           n_iter = 2000, warmup = 1000, chains = 2, seed = 1)

# diagnostics, plots, and interoperability
crr_diagnostics(fit)                 # structured MCMC + regression report
plot(fit, type = "trace")            # also: "acf", "violin", "ess", "ess_time", "residual"
coda::as.mcmc(fit)                   # -> coda ecosystem
posterior::as_draws(fit)             # -> posterior / bayesplot ecosystem
predict(fit, type = "response")

# benchmark the constrained sampler against the unconstrained probit baseline
crr_benchmark(sim$Y, sim$X, con, beta = sim$beta,
              n_iter = 2000, warmup = 1000, seed = 1)

Overview

Two vignettes and a R script cover the details:

• vignette("combreg"): the model, why constraint-ignoring baselines are biased, constraint setup, the adaptive block controller, and building custom models from the exported primitives.
• vignette("diagnostics"): the diagnostics and benchmarking workflow.
• inst/scripts/reproduce-table2.R: runs the simulation in Zheng, Ghosh & Duan (2026) and reproduces its Table 2 (coefficient RMSE for varying (d, m)).

Citation

1

citation("combreg")

Zheng, Y., Ghosh, M., & Duan, L. (2026). Statistical Modeling of Combinatorial Response Data. Journal of the American Statistical Association, advance online publication. https://doi.org/10.1080/01621459.2026.2714575

License

MIT © Hugh Zheng