## ----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)

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

## -----------------------------------------------------------------------------
report <- crr_diagnostics(fit, beta = sim$beta)
report

## -----------------------------------------------------------------------------
crr_ess(fit)

## -----------------------------------------------------------------------------
crr_rhat(fit)

## -----------------------------------------------------------------------------
ppc <- crr_ppc(fit, n_rep = 100, seed = 1)
ppc

## -----------------------------------------------------------------------------
tab <- report$table
head(tab[order(tab$ess), c("parameter", "ess", "ess_per_sec", "rhat")])
head(tab[order(-tab$rhat), c("parameter", "rhat", "ess")])

## -----------------------------------------------------------------------------
plot(fit, type = "trace", pars = c("beta[1,1]", "beta[2,1]"))
plot(fit, type = "acf", pars = c("beta[1,1]", "beta[2,1]"))

## -----------------------------------------------------------------------------
plot(fit, type = "violin")

## -----------------------------------------------------------------------------
plot(fit, type = "ess")
plot(fit, type = "ess_time")

## ----fig.height = 5.5---------------------------------------------------------
plot(fit, type = "residual")

## -----------------------------------------------------------------------------
plot(fit, type = "beta_diff", beta = sim$beta)

## -----------------------------------------------------------------------------
bm <- crr_benchmark(sim$Y, sim$X, con, beta = sim$beta,
                    n_iter = 500, warmup = 250, seed = 42)
bm

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

## -----------------------------------------------------------------------------
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)

## -----------------------------------------------------------------------------
plot(fit_m, type = "ess_time")

## ----fig.height = 4-----------------------------------------------------------
plot(fit_m, type = "beta_diff", beta = sim_m$beta)

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