## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6,
                      fig.height = 4)
library(bbqr)
set.seed(2026)

## ----fit----------------------------------------------------------------------
n <- 300
X <- matrix(rnorm(n * 5), n, 5, dimnames = list(NULL, paste0("x", 1:5)))
truth <- c(1.5, -1, 0.8, 0, 0)
y <- as.numeric(X %*% truth + rnorm(n) > 0)
dat <- data.frame(y = y, X)

fit <- bbqr(y ~ x1 + x2 + x3 + x4 + x5, data = dat,
            quantile = 0.5, penalty = "alasso",
            ndraw = 2000, burn = 500)
fit

## ----summary------------------------------------------------------------------
summary(fit)

## ----penalties----------------------------------------------------------------
pens <- c("none", "lasso", "alasso")
est <- sapply(pens, function(p)
  coef(bbqr(y ~ x1 + x2 + x3 + x4 + x5, data = dat, penalty = p,
            ndraw = 2000, burn = 500)))
round(est, 3)

## ----penalty-direction--------------------------------------------------------
cosine <- function(a, b) sum(a * b) / sqrt(sum(a^2) * sum(b^2))
round(apply(est[-1, ], 2, cosine, truth), 4)

## ----hierarchies--------------------------------------------------------------
prior("alasso")
prior("alasso", model = "v3")

## ----omega--------------------------------------------------------------------
om <- sapply(c(v3 = "v3", v5 = "v5"), function(m)
  bbqr(y ~ x1 + x2 + x3 + x4 + x5, dat, ndraw = 4000, burn = 1000,
       prior = prior("alasso", model = m))$omega)
signif(apply(om, 2, quantile, c(0, 0.01, 0.1, 0.5)), 3)

## ----anchors-table, echo = FALSE----------------------------------------------
knitr::kable(data.frame(
  anchor = c("sigma1", "free", "beta1", "norm1", "normslopes"),
  restriction = c("ALD inverse-scale held at 1", "none beyond the prior",
                  "first slope held at 1",
                  "||(beta0, beta)|| = 1", "||beta|| = 1, slopes only"),
  available = c("all", "lasso, alasso", "all", "all", "all")
), caption = "Anchors, and which penalties accept them.")

## ----anchor-demo--------------------------------------------------------------
f_sig  <- bbqr(y ~ x1 + x2 + x3 + x4 + x5, dat, anchor = "sigma1",
               ndraw = 2000, burn = 500)
f_norm <- bbqr(y ~ x1 + x2 + x3 + x4 + x5, dat, anchor = "norm1",
               ndraw = 2000, burn = 500)

c(sigma1_scale_fixed  = all(f_sig$sigma == 1),
  norm1_unit_norm     = all(abs(sqrt(f_norm$beta0^2 +
                                     rowSums(f_norm$beta^2)) - 1) < 1e-8),
  norm1_sigma_moves   = sd(f_norm$sigma) > 0,
  sigma1_derived      = f_sig$derived,
  norm1_derived       = f_norm$derived)

## ----easy---------------------------------------------------------------------
anch <- c("sigma1", "free", "beta1", "norm1", "normslopes")
easy <- sapply(anch, function(a)
  cosine(coef(suppressWarnings(
    bbqr(y ~ x1 + x2 + x3 + x4 + x5, dat, anchor = a,
         ndraw = 2000, burn = 500)))[-1], truth))
round(easy, 4)

## ----hard---------------------------------------------------------------------
lat  <- as.vector(X %*% truth) + rt(n, df = 3)
yh   <- as.numeric(lat + quantile(-lat, 0.85) > 0)
hard <- data.frame(y = yh, X)

res <- sapply(anch, function(a)
  cosine(coef(suppressWarnings(
    bbqr(y ~ x1 + x2 + x3 + x4 + x5, hard, quantile = 0.05,
         anchor = a, ndraw = 3000, burn = 750)))[-1], truth))
round(res, 4)

## ----other--------------------------------------------------------------------
p <- predict(fit, newdata = dat, type = "response")
round(quantile(p, c(0, .25, .5, .75, 1)), 3)

mean(predict(fit, newdata = dat, type = "class") == y)

## ----multi--------------------------------------------------------------------
mfit <- bbqr(y ~ x1 + x2 + x3 + x4 + x5, dat, quantile = c(0.25, 0.5, 0.75),
             ndraw = 1500, burn = 400)
round(coef(mfit), 3)

## ----cont---------------------------------------------------------------------
cd <- data.frame(x1 = rnorm(150), x2 = rnorm(150))
cd$y <- 2 + 1.5 * cd$x1 - cd$x2 + rnorm(150)

cfit <- cbqr(y ~ x1 + x2, cd, quantile = 0.5, penalty = "alasso",
             ndraw = 1500, burn = 400)
round(coef(cfit), 3)

## ----cont-sigma---------------------------------------------------------------
round(mean(cfit$sigma), 3)

## ----cont-pred----------------------------------------------------------------
round(head(predict(cfit, newdata = cd), 4), 3)