## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 10,
  fig.height = 5,
  fig.align = "center",
  warning = FALSE,
  message = FALSE,
  digits = 4
)

kbl10 <- function(x, digits = 4, ...) {
  knitr::kable(utils::head(as.data.frame(x), 10), digits = digits, align = "c", ...)
}

## ----library------------------------------------------------------------------
library(betaregscale)

## ----simulate-fit-------------------------------------------------------------
set.seed(2026)
n <- 220
d <- data.frame(
  x1 = rnorm(n),
  x2 = rnorm(n),
  z1 = rnorm(n)
)

sim <- brs_sim(
  formula = ~ x1 + x2 | z1,
  data = d,
  beta = c(0.20, -0.45, 0.25),
  zeta = c(0.15, -0.30),
  ncuts = 100,
  repar = 2
)

fit_fixed <- brs(y ~ x1 + x2, data = sim, repar = 2)
fit_var <- brs(y ~ x1 + x2 | z1, data = sim, repar = 2)

## ----table--------------------------------------------------------------------
tab <- brs_table(
  fixed = fit_fixed,
  variable = fit_var,
  sort_by = "AIC"
)
kbl10(tab)

## ----me-----------------------------------------------------------------------
set.seed(2026) # For marginal effects simulation
me_mean <- brs_marginaleffects(
  fit_var,
  model = "mean",
  type = "response",
  interval = TRUE,
  n_sim = 120
)
kbl10(me_mean)

set.seed(2026) # Reset seed for reproducibility
me_precision <- brs_marginaleffects(
  fit_var,
  model = "precision",
  type = "link",
  interval = TRUE,
  n_sim = 120
)
kbl10(me_precision)

## ----scoreprob----------------------------------------------------------------
P <- brs_predict_scoreprob(fit_var, scores = 0:10)
dim(P)
kbl10(P[1:6, 1:5])
kbl10(
  data.frame(row_sum = rowSums(P)),
  digits = 4
)

## ----autoplot, eval = requireNamespace("ggplot2", quietly = TRUE)-------------
autoplot.brs(fit_var, type = "calibration")
autoplot.brs(fit_var, type = "score_dist", scores = 0:20)
autoplot.brs(fit_var, type = "cdf", max_curves = 4)
autoplot.brs(fit_var, type = "residuals_by_delta", residual_type = "rqr")

## ----cv-----------------------------------------------------------------------
set.seed(2026) # For cross-validation reproducibility
cv_res <- brs_cv(
  y ~ x1 + x2 | z1,
  data = sim,
  k = 3,
  repeats = 1,
  repar = 2
)

kbl10(cv_res)
kbl10(
  data.frame(
    metric = c("log_score", "rmse_yt", "mae_yt"),
    mean = colMeans(cv_res[, c("log_score", "rmse_yt", "mae_yt")], na.rm = TRUE)
  ),
  digits = 4
)