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

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

## ----likert-------------------------------------------------------------------
likert <- c(1, 2, 3, 5, 4, 1)
brs_check(likert - 1, ncuts = 4)

## ----repar--------------------------------------------------------------------
brs_repar(mu = 0.3, phi = 0.2, repar = 2)         # mean 0.3, dispersion 0.2
brs_repar(mu = 0.3, phi = 1 / 0.2 - 1, repar = 1) # same law: precision 4

## ----check--------------------------------------------------------------------
brs_check(c(0, 3, 5, 10), ncuts = 10)
brs_check(c(0, 3, 5, 10), ncuts = 10, interval = "right")

## ----fit----------------------------------------------------------------------
set.seed(4255)
d <- data.frame(x1 = rnorm(200), x2 = rbinom(200, 1, 0.5))
sim <- brs_sim(~ x1 + x2, data = d, beta = c(-0.4, 0.5, -0.6),
               phi = qlogis(0.15), ncuts = 10)
head(sim)

fit <- brs(y ~ x1 + x2, data = sim)
summary(fit)

## ----pmf----------------------------------------------------------------------
P <- brs_predict_scoreprob(fit)
round(P[1:3, ], 3)
range(rowSums(P))
c(sum = sum(brs_loglik_contrib(fit)), logLik = as.numeric(logLik(fit)))

## ----coef---------------------------------------------------------------------
coef(fit, model = "mean")
confint(fit)

## ----diagnostics--------------------------------------------------------------
str(fit$diagnostics[c("grad_gain", "min_eig_scaled", "n_clamped")])

## ----prep---------------------------------------------------------------------
d_prep <- data.frame(
  y     = c(5, 0, NA, NA, 5),
  delta = c(3, NA, NA, NA, 0),
  left  = c(NA, NA, 2, NA, NA),
  right = c(NA, NA, 5, 3, NA)
)
brs_prep(d_prep, ncuts = 10)

## ----maps---------------------------------------------------------------------
raw <- data.frame(y = sim$y, x1 = sim$x1, x2 = sim$x2)
fits <- lapply(c(mid = "mid", right = "right", left = "left"), function(m)
  brs(y ~ x1 + x2, data = raw, ncuts = 10, interval = m))
data.frame(logLik = sapply(fits, logLik), AIC = sapply(fits, AIC),
           BIC = sapply(fits, BIC))
tryCatch(anova(fits$mid, fits$right), error = conditionMessage)

## ----boundary-----------------------------------------------------------------
set.seed(1)
heap <- brs_sim(~ x1 + x2, data = d, beta = c(0, 0.5, -0.6),
                phi = qlogis(0.2), ncuts = 10, cutpoints = c(0.12, 0.88))
table(heap$y)

fit_fix <- brs(y ~ x1 + x2, data = heap)
fit_bnd <- brs(y ~ x1 + x2, data = heap, cutpoints = "boundary")
anova(fit_fix, fit_bnd)
coef(fit_bnd, model = "cutpoints")