Package 'betaregscale'

Akaike information criterion

Description

Akaike information criterion

Usage

## S3 method for class 'brs'
AIC(object, ..., k = 2)

Arguments

object	A fitted "brs" object.
...	Ignored.
k	Penalty per parameter (default 2).

Value

Scalar AIC value.

See Also

brs, logLik.brs, BIC.brs, brs_gof

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10)
)
prep <- brs_prep(dat, ncuts = 100)
fit <- brs(y ~ x1, data = prep)
AIC(fit)

---

AIC for brsmm models

Description

AIC for brsmm models

Usage

## S3 method for class 'brsmm'
AIC(object, ..., k = 2)

Arguments

object	A fitted "brsmm" object.
...	Currently ignored.
k	Numeric penalty per parameter.

Value

Numeric scalar.

See Also

brsmm, logLik.brsmm, BIC.brsmm, brs_gof

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10),
  id = factor(rep(1:4, each = 5))
)
prep <- brs_prep(dat, ncuts = 100)
fit <- brsmm(y ~ x1, random = ~ 1 | id, data = prep)
AIC(fit)

---

Likelihood-ratio comparison of nested beta interval models

Description

Compares fitted "brs" and "brsmm" models by log-likelihood, AIC, BIC and likelihood-ratio tests (Lopes, 2023, "Inferencia").

Usage

## S3 method for class 'brs'
anova(object, ..., test = c("Chisq", "none"))

Arguments

object	A fitted "brs" model.
...	Further fitted "brs" and/or "brsmm" models.
test	"Chisq" (default) or "none".

Details

The models are sorted by their number of parameters. For consecutive models the statistic is $LR = 2(\ell_1 - \ell_0)$, with Chi Df the difference in the number of parameters, and $p = P(\chi^2_{df} > LR)$. The models must be nested; this is not checked. They must also describe the same response: the same observations, interval, ncuts and (under "mid") lim, or the call stops, since different coarsenings are different likelihoods.

When the larger model adds one random-effect term (a "brs" model against a random-intercept "brsmm", or one more correlated random term), its variance lies on the boundary of the parameter space under $H_0$ and $LR$ follows the mixture $\frac12\chi^2_{df-1} + \frac12\chi^2_{df}$ (Self and Liang, 1987; Stram and Lee, 1994); for one variance component alone this is $\frac12\chi^2_0 + \frac12\chi^2_1$, i.e. half the naive p-value. The printed heading names the rows that use it. When more than one random term is added at once, the naive $\chi^2_{df}$ p-value is kept and flagged as conservative.

Value

An object of class "anova" (a data frame) with columns Df, logLik, AIC, BIC and, for test = "Chisq", Chisq, Chi Df and Pr(>Chisq); the attribute "heading" explains the p-values.

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

Self, S. G., and Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82(398), 605–610. doi:10.1080/01621459.1987.10478472

Stram, D. O., and Lee, J. W. (1994). Variance components testing in the longitudinal mixed effects model. Biometrics, 50(4), 1171–1177. doi:10.2307/2533455

See Also

anova.brsmm, summary.brs, logLik.brs

Examples

# Synthetic NRS-11 scores: 4 groups x 3 times. Simulated, not real data.
set.seed(2023)
nrs <- expand.grid(id = 1:80, time = c("6h", "12h", "24h"))
nrs$group <- factor(paste0("g", (nrs$id - 1) %% 4 + 1))
eta <- -1.3 + c(0, 0.75, 0.3)[nrs$time] + c(0, -0.1, 0.05, 0.1)[nrs$group]
shp <- brs_repar(mu = plogis(eta), phi = 0.3, repar = 2)
nrs$y <- round(10 * rbeta(nrow(nrs), shp$shape1, shp$shape2))

# Time only (m1) nested in time + group (m2): LR test on 3 df
m1 <- brs(y ~ time, data = nrs, ncuts = 10)
m2 <- brs(y ~ time + group, data = nrs, ncuts = 10)
anova(m1, m2)

# Fits under different interval directions are different response models
m2_right <- brs(y ~ time + group, data = nrs, ncuts = 10, interval = "right")
try(anova(m2, m2_right))

---

Likelihood-ratio comparison involving mixed models

Description

anova() for "brsmm" fits: the same table as anova.brs, with the chi-bar-square mixture $\frac12\chi^2_{df-1} + \frac12\chi^2_{df}$ for rows that add one random-effect term, whose variance is on the boundary under $H_0$. This is the test to use for a variance component: the Wald statistic of its log standard deviation is not meaningful (see summary.brsmm).

Usage

## S3 method for class 'brsmm'
anova(object, ..., test = c("Chisq", "none"))

Arguments

object	A fitted "brsmm" model.
...	Further fitted "brsmm" and/or "brs" models.
test	"Chisq" (default) or "none".

Value

An object of class "anova"; see anova.brs.

References

Self, S. G., and Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82(398), 605–610. doi:10.1080/01621459.1987.10478472

See Also

anova.brs, brsmm, summary.brsmm

Examples

set.seed(11)
g <- 20
d <- data.frame(id = factor(rep(1:g, each = 8)), x = runif(8 * g))
shp <- brs_repar(plogis(-0.4 + d$x + rnorm(g, sd = 0.6)[d$id]), phi = 0.25)
d$y <- round(10 * rbeta(nrow(d), shp$shape1, shp$shape2))
m0 <- brs(y ~ x, data = d, ncuts = 10)
m1 <- brsmm(y ~ x, random = ~ 1 | id, data = d, ncuts = 10)
anova(m0, m1)  # Pr(>Chisq) = half the chi2(1) tail

---

ggplot2 autoplot for brs models

Description

Produces ggplot2 diagnostics tailored to interval-censored scale models.

Usage

## S3 method for class 'brs'
autoplot(
  object,
  type = c("calibration", "score_dist", "cdf", "residuals_by_delta", "all"),
  bins = 10L,
  scores = NULL,
  newdata = NULL,
  n_grid = 200L,
  max_curves = 6L,
  residual_type = "rqr",
  title = NULL,
  xlab = NULL,
  ylab = NULL,
  ncol = 2L,
  theme = ggplot2::theme_minimal(),
  ...
)

Arguments

object	A fitted "brs" object.
type	Plot type: "calibration", "score_dist", "cdf", "residuals_by_delta", or "all" (produces all panels in a single grid).
bins	Number of bins for "calibration".
scores	Optional integer vector of scores for "score_dist". Defaults to all scores from 0 to ncuts.
newdata	Optional data frame of covariate scenarios for "cdf".
n_grid	Number of points on $(0,1)$ used to draw CDF curves.
max_curves	Maximum number of CDF curves shown when newdata is not provided.
residual_type	Residual type for "residuals_by_delta"; passed to residuals.brs.
title	Optional character: override the plot title via ggplot2::labs(title = ...). Ignored when type = "all".
xlab	Optional character: override the x-axis label. Ignored when type = "all".
ylab	Optional character: override the y-axis label. Ignored when type = "all".
ncol	Number of columns for the grid when type = "all". Defaults to 2.
theme	A ggplot2 theme object (e.g., ggplot2::theme_bw()) or a theme function. Applied to every panel. Defaults to ggplot2::theme_minimal().
...	Additional arguments forwarded to ggplot2::theme() and applied on top of theme. Use named theme element arguments, e.g. legend.position = "none".

Details

type = "calibration" bins predictions and compares mean observed vs mean predicted response in each bin.

type = "score_dist" compares observed score frequencies against expected frequencies implied by the fitted beta interval model.

Value

A ggplot2 object.

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

Hawker, G. A., Mian, S., Kendzerska, T., and French, M. (2011). Measures of adult pain: Visual Analog Scale for Pain (VAS Pain), Numeric Rating Scale for Pain (NRS Pain), McGill Pain Questionnaire (MPQ), Short-Form McGill Pain Questionnaire (SF-MPQ), Chronic Pain Grade Scale (CPGS), Short Form-36 Bodily Pain Scale (SF-36 BPS), and Measure of Intermittent and Constant Osteoarthritis Pain (ICOAP). Arthritis Care and Research, 63(S11), S240-S252. doi:10.1002/acr.20543

Hjermstad, M. J., Fayers, P. M., Haugen, D. F., et al. (2011). Studies comparing Numerical Rating Scales, Verbal Rating Scales, and Visual Analogue Scales for assessment of pain intensity in adults: a systematic literature review. Journal of Pain and Symptom Management, 41(6), 1073-1093. doi:10.1016/j.jpainsymman.2010.08.016

See Also

brs, plot.brs, autoplot.brs_bootstrap

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10),
  x2 = rep(c(0, 0, 1, 1), 5)
)
prep <- brs_prep(dat, ncuts = 100)
fit <- brs(y ~ x1 + x2, data = prep)
ggplot2::autoplot(fit, type = "calibration")
ggplot2::autoplot(fit,
  type = "calibration",
  title = "My calibration", ylab = "Observed"
)
ggplot2::autoplot(fit, type = "all")

---

ggplot2 autoplot for bootstrap results

Description

Produces visual summaries for objects returned by brs_bootstrap.

Usage

## S3 method for class 'brs_bootstrap'
autoplot(
  object,
  type = c("ci_forest", "dist", "qq", "stability"),
  parameter = NULL,
  title = NULL,
  caption = NULL,
  max_parameters = 12L,
  ci_level = NULL,
  theme = NULL,
  ...
)

Arguments

object	An object of class "brs_bootstrap".
type	Plot type: "ci_forest", "dist", "qq", or "stability".
parameter	Optional parameter name used by type = "dist", "qq", and "stability". If NULL, the first parameter is used.
title	Optional plot title override.
caption	Optional subtitles/titles for plot types. Accepts: a single string (used for the selected type); a character vector/list with up to four entries in the order ci_forest, dist, qq, stability.
max_parameters	Maximum number of parameters shown in type = "ci_forest".
ci_level	Confidence level used in type = "stability". Defaults to the level stored in object.
theme	Optional ggplot2 theme object (e.g., ggplot2::theme_bw()). If NULL, ggplot2::theme_minimal() is used.
...	Currently ignored.

Details

For type = "dist", "qq", and "stability", bootstrap draws must be present in attr(object, "boot_draws"), obtained by fitting with brs_bootstrap(..., keep_draws = TRUE).

Value

A ggplot2 object.

See Also

brs_bootstrap, brs, autoplot.brs

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10)
)
prep <- brs_prep(dat, ncuts = 100)
fit <- brs(y ~ x1, data = prep)
boot <- brs_bootstrap(fit, R = 50)
ggplot2::autoplot(boot, type = "ci_forest")

---

ggplot2 autoplot for marginal effects

Description

Produces visual summaries for objects returned by brs_marginaleffects.

Usage

## S3 method for class 'brs_marginaleffects'
autoplot(
  object,
  type = c("forest", "magnitude", "dist"),
  variable = NULL,
  top_n = 12L,
  title = NULL,
  caption = NULL,
  theme = NULL,
  ...
)

Arguments

object	An object of class "brs_marginaleffects".
type	Plot type: "forest", "magnitude", or "dist".
variable	Optional variable name for type = "dist".
top_n	Maximum number of variables shown in "magnitude" (ordered by |AME|).
title	Optional plot title override.
caption	Optional subtitle override.
theme	Optional ggplot2 theme object. If NULL, ggplot2::theme_minimal() is used.
...	Currently ignored.

Details

type = "dist" requires AME simulation draws stored in attr(object, "ame_draws"), which are available when marginal effects are computed with keep_draws = TRUE and interval = TRUE.

Value

A ggplot2 object.

See Also

brs_marginaleffects, brs, autoplot.brs

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10)
)
prep <- brs_prep(dat, ncuts = 100)
fit <- brs(y ~ x1, data = prep)
ame <- brs_marginaleffects(fit)
ggplot2::autoplot(ame, type = "forest")

---

ggplot2 autoplot for brsmm models

Description

Produces ggplot2 diagnostics tailored to mixed beta interval models.

Usage

## S3 method for class 'brsmm'
autoplot(
  object,
  type = c("calibration", "score_dist", "ranef_qq", "residuals_by_group",
    "ranef_caterpillar", "ranef_density", "ranef_pairs", "shrinkage", "all"),
  bins = 10L,
  scores = NULL,
  residual_type = c("response", "pearson"),
  max_groups = 25L,
  title = NULL,
  xlab = NULL,
  ylab = NULL,
  ncol = 2L,
  theme = ggplot2::theme_minimal(),
  ...
)

Arguments

object	A fitted "brsmm" object.
type	Plot type. One of "calibration", "score_dist", "ranef_qq", "residuals_by_group", "ranef_caterpillar", "ranef_density", "ranef_pairs", "shrinkage", or "all" (produces all panels in a single grid).
bins	Number of bins for "calibration".
scores	Optional integer vector of scores for "score_dist". Defaults to all scores from 0 to ncuts.
residual_type	Residual type for "residuals_by_group"; passed to residuals.brsmm.
max_groups	Maximum groups in "residuals_by_group".
title	Optional character: override the plot title via ggplot2::labs(title = ...). Ignored when type = "all".
xlab	Optional character: override the x-axis label. Ignored when type = "all".
ylab	Optional character: override the y-axis label. Ignored when type = "all".
ncol	Number of columns for the grid when type = "all". Defaults to 2.
theme	A ggplot2 theme object (e.g., ggplot2::theme_bw()) or a theme function. Applied to every panel. Defaults to ggplot2::theme_minimal().
...	Additional arguments forwarded to ggplot2::theme() and applied on top of theme. Use named theme element arguments, e.g. legend.position = "none".

Value

A ggplot2 object, or (when type = "all") invisibly returns the list of all panels after rendering the grid.

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10),
  id = factor(rep(1:4, each = 5))
)
prep <- brs_prep(dat, ncuts = 100)
fit_mm <- brsmm(y ~ x1, random = ~ 1 | id, data = prep)
ggplot2::autoplot(fit_mm, type = "calibration", bins = 4)
ggplot2::autoplot(fit_mm,
  type = "ranef_caterpillar",
  title = "My title", ylab = "Mode"
)
ggplot2::autoplot(fit_mm, type = "all")

---

Bayesian information criterion

Description

Bayesian information criterion

Usage

## S3 method for class 'brs'
BIC(object, ...)

Arguments

object	A fitted "brs" object.
...	Ignored.

Value

Scalar BIC value.

See Also

brs, logLik.brs, AIC.brs, brs_gof

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10)
)
prep <- brs_prep(dat, ncuts = 100)
fit <- brs(y ~ x1, data = prep)
BIC(fit)

---

BIC for brsmm models

Description

$-2\ell + \log(n)\,k$ with $n$ = the number of observations (nobs(object)), as lme4 does, and $k$ the number of parameters (fixed effects, precision and packed random-effect parameters). There is no single sample size for a mixed model: counting groups instead ($n =$ object$ngroups) penalises more and is a common alternative, -2 * logLik(object) + log(object$ngroups) * k. Compare BIC values only between fits with the same convention.

Usage

## S3 method for class 'brsmm'
BIC(object, ...)

Arguments

object	A fitted "brsmm" object.
...	Currently ignored.

Value

Numeric scalar.

See Also

brsmm, logLik.brsmm, AIC.brsmm, brs_gof

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10),
  id = factor(rep(1:4, each = 5))
)
prep <- brs_prep(dat, ncuts = 100)
fit <- brsmm(y ~ x1, random = ~ 1 | id, data = prep)
BIC(fit)

---

Fit a beta interval regression model

Description

Fits by maximum likelihood a beta regression for scores on a bounded scale, treated as coarsened (interval-censored) observations of a latent beta variable (Lopes, 2023). A one-part formula y ~ x1 + x2 keeps the second parameter constant (brs_fit_fixed); a two-part formula y ~ x1 + x2 | z1 + z2 also models it (brs_fit_var).

Usage

brs(
  formula,
  data,
  link = NULL,
  link_phi = NULL,
  ncuts = NULL,
  lim = NULL,
  repar = 2L,
  method = c("BFGS", "L-BFGS-B"),
  hessian_method = c("cpp", "numDeriv", "optim"),
  interval = NULL,
  start = NULL,
  control = list()
)

Arguments

formula	A Formula-style formula with two parts: y ~ x1 + x2 | z1 + z2.
data	Data frame.
link	Link for the first parameter (the mean under repar = 1, 2; the shape $p$ under repar = 0). NULL (default) selects the link implied by repar: "logit" for repar = 1, 2, "log" for repar = 0. See the 'Reparameterizations and links' section of brs for the admissible values.
link_phi	Link for the second parameter. NULL (default) selects "logit" for repar = 2 (dispersion on $(0, 1)$) and "log" for repar = 0, 1 (positive shape/precision).
ncuts	Integer $K$: the maximum score, so that the scale is $0, 1, \ldots, K$ ($K + 1$ categories). NULL (default) uses the value stored by brs_prep or brs_sim, or 100 for raw data.
lim	Half-width of the score cell in $(0, 0.5]$ (interval = "mid" only). NULL (default) uses attr(data, "lim") from brs_prep, or 0.5; same rule as ncuts. Values below 0.5 warn (partial coarsening).
repar	Reparameterization scheme (default 2); see brs_repar.
method	Optimization method (default "BFGS").
hessian_method	How the Hessian for vcov() is computed: "cpp" (default; compiled chain rule on the linear predictors), "numDeriv" (numDeriv::hessian() of the log-likelihood) or "optim" (the finite-difference Hessian returned by optim).
interval	Direction of the uncertainty interval, "mid", "right" or "left" (see brs_check). NULL (default) uses attr(data, "interval") from brs_prep, or "mid"; same rule as ncuts.
start	Optional numeric vector of starting values (mean coefficients, then dispersion coefficients). NULL (default) uses compute_start(). Refits (bootstrap, jackknife) pass the parent estimate here.
control	Control list for optim; its entries are merged into the default list(maxit = 5000L).

Details

The scores run over $0, 1, \ldots, K$ with $K =$ ncuts: the scale has $K + 1$ categories and $K$ is its maximum. A scale that starts at 1, such as a Likert item 1–5, must be shifted to 0–4 and fitted with ncuts = 4; otherwise the lowest category is read as an interior score and the censoring of the lower border is lost.

Raw scores are converted to cells and censoring types by brs_check with ncuts, lim and interval; data from brs_prep or brs_sim are used as they are, with their stored ncuts, lim and interval. Values in $(0, 1)$ are exact observations. The model is

$Y_i \sim \mathrm{Beta}(a_i, b_i), \qquad g_1(\mu_i) = x_i^\top \beta, \qquad g_2(\phi_i) = z_i^\top \gamma,$

with $(a_i, b_i)$ obtained from $(\mu_i, \phi_i)$ by brs_repar.

Value

An object of class "brs": a list with, among others, par (estimates on the link scales), coefficients (mean and precision parts), value (maximised log-likelihood), hessian, convergence, diagnostics (see 'Fit diagnostics'), hatmu (first parameter per observation: the mean, or the shape $p$ under repar = 0), hatphi, Y (columns left, right, yt, y, delta), ncuts, lim, interval, repar, link and link_phi.

Likelihood

With cell $[l_i, u_i]$ and censoring type $\delta_i$,

$L(\beta, \gamma) = \prod_{i=1}^n f(y_i)^{I(\delta_i = 0)}\, F(u_i)^{I(\delta_i = 1)}\, \{1 - F(l_i)\}^{I(\delta_i = 2)}\, \{F(u_i) - F(l_i)\}^{I(\delta_i = 3)},$

where $f$ and $F$ are the beta density and distribution function with shapes $(a_i, b_i)$. The four factors are exact values ($\delta = 0$), the lowest score ($\delta = 1$, left-censored), the highest score ($\delta = 2$, right-censored) and the other scores ($\delta = 3$, interval-censored). This is eq. eqn_verossimilhanca_geral of Lopes (2023); the table of censoring types in the same section swaps the labels $\delta = 1$ and $\delta = 2$, and the package follows the equation.

Numerical details: endpoints are clamped to $[10^{-5}, 1 - 10^{-5}]$ and there is no probability floor; an interval probability uses lower-tail probabilities when the interval midpoint is at or below the mean $a/(a + b)$ and upper tails otherwise, which avoids cancellation; below $10^{-240}$ an endpoint Laplace approximation of the tail replaces it; a non-finite contribution becomes $-10^6$, as does a NaN parameter. The mean and the dispersion (repar = 2) are clamped to $[10^{-5}, 1 - 10^{-5}]$, the precision and the shape $p$ to $[10^{-5}, 10^8]$, and the beta shapes to $[10^{-12}, 10^8]$.

Estimation

optim (method = "BFGS", the default, or "L-BFGS-B"; maxit = 5000 unless control says otherwise) maximises the log-likelihood, which is evaluated in C++. The gradient uses the chain rule on the two linear predictors: the derivatives of each contribution with respect to $\eta_{1i}$ and $\eta_{2i}$ are Richardson central differences with step $10^{-4}\max(1, |\eta|)$, and the gradient is $X^\top d_1$ and $Z^\top d_2$. The Hessian behind vcov() is built the same way from the per-observation second derivatives (step $3 \times 10^{-4}\max(1, |\eta|)$; hessian_method = "cpp", the default):

$H = \left(\begin{array}{cc} X^\top W_{11} X & X^\top W_{12} Z \\ Z^\top W_{12} X & Z^\top W_{22} Z \end{array}\right),$

with $W_{jk}$ diagonal matrices of the second derivatives of each contribution in $(\eta_{1i}, \eta_{2i})$. "numDeriv" differentiates the log-likelihood with numDeriv::hessian() and "optim" takes the finite-difference Hessian of optim. Lopes (2023, "Estimacao") also uses BFGS with numerical derivatives; here they are taken on the linear predictors, not on the coefficients.

Starting values: start when given (bootstrap and jackknife refits pass the parent estimate); otherwise a quasi-binomial GLM of the cell midpoints for the mean, and for the second parameter the moment estimate on the midpoints as the intercept with zero slopes; under repar = 0 both shapes by the method of moments. The covariance matrix is $(-H)^{-1}$ at the estimate (vcov.brs); no generalised inverse is used.

Fit diagnostics

Before optim, the mean and precision model matrices are checked (pivoted QR, tolerance $10^{-7}$; condition number of the matrix with unit-length columns). After it, fit$diagnostics holds grad_norm (largest absolute gradient), grad_gain (the log-likelihood a Newton step would still gain, $\frac12 g^\top (-H)^{-1} g$), grad_step (that step in standard errors), min_eig, max_eig and min_eig_scaled (eigenvalues of $-H$, the last in correlation form), hessian_nd, n_clamped and clamped (observations on the clamps of the likelihood). Each problem gives one line:

“model matrix is rank deficient” (error)

Some columns are linear combinations of others (the message names them). Remove or recode them.

“model matrix is nearly collinear”

Condition number above $10^4$: estimates and standard errors are unstable. Drop or combine the named columns, or centre them.

“Optimizer did not converge”

optim stopped at its iteration limit or failed. Try method = "L-BFGS-B", rescale the covariates or simplify the model.

“Gradient not ~0 at the estimate”

A Newton step would still gain more than 0.01 in log-likelihood although optim reported convergence. Refit with the other method; the message asks to rescale the covariates when their scales differ by more than $10^3$ (e.g. a raw income next to a dummy).

“Hessian near-singular or not negative definite (SEs unreliable)”

The likelihood is nearly flat or curved the wrong way in some direction: a parameter that is not identified, or covariates on very different scales (the message then asks to rescale them). vcov() returns NA for variances it cannot estimate; use likelihood-ratio tests (anova.brs) or brs_bootstrap instead of Wald statements, rescale, or simplify the model.

“observations ... on the clamp boundary”

Fitted means, dispersions or shapes sit on the numerical limits above, typically because a group of observations lies entirely at one end of the scale (separation) or every observation is censored on the same side. The coefficients drift towards infinity and are not interpretable; merge sparse groups or remove the separating covariate.

brsmm adds a check of the variance components.

Reparameterizations and links

repar	parameters (Lopes, 2023)	link	link_phi
0	shapes $p, q > 0$; eq. eqn_beta_p1 (regression on the shapes is a package extension)	log, sqrt	log, sqrt
1	mean $\mu \in (0, 1)$, precision $\phi > 0$; "parametrizacao 1" (Ferrari and Cribari-Neto, 2004)	logit, probit, cauchit, cloglog	log, sqrt
2	mean $\mu \in (0, 1)$, dispersion $\phi \in (0, 1)$; "parametrizacao 2" (Bayer, 2011), written $\sigma$ there	logit, probit, cauchit, cloglog	logit, probit, cauchit, cloglog

The first link of each cell is the default (link = NULL, link_phi = NULL); any other combination is an error. "identity", "inverse" and "1/mu^2" are not accepted for positive parameters. With "sqrt" the inverse link is flat for $\eta \le 0$, and a warning is issued when a fitted linear predictor lies there. The mean and variance of $Y$ under each scheme are given in brs_repar. Under repar = 0 the object stores the shape $p$ in hatmu, while fitted(), predict(type = "response"), residuals and marginal effects use $E[Y] = p/(p + q)$.

Interval direction

interval = "mid" (default) uses the cells $[s - \mathrm{lim}, s + \mathrm{lim}]/K$; "right" and "left" use $K + 1$ equal cells $[s, s + 1]/(K + 1)$ (the dissertation's $r$ and $l$, with a package normalisation). "right" and "left" give the same fit and differ only in predict(type = "score"); "mid" and "right"/"left" are different coarsening models, whose log-likelihoods and AIC are not comparable. Details: brs_check.

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

Ferrari, S. L. P., and Cribari-Neto, F. (2004). Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31(7), 799–815. doi:10.1080/0266476042000214501

Bayer, F. M. (2011). Modelagem e inferencia em regressao beta. PhD thesis, Universidade Federal de Pernambuco.

Hawker, G. A., Mian, S., Kendzerska, T., and French, M. (2011). Measures of adult pain: Visual Analog Scale for Pain (VAS Pain), Numeric Rating Scale for Pain (NRS Pain), McGill Pain Questionnaire (MPQ), Short-Form McGill Pain Questionnaire (SF-MPQ), Chronic Pain Grade Scale (CPGS), Short Form-36 Bodily Pain Scale (SF-36 BPS), and Measure of Intermittent and Constant Osteoarthritis Pain (ICOAP). Arthritis Care and Research, 63(S11), S240-S252. doi:10.1002/acr.20543

See Also

brs_check, brs_prep, brs_repar, summary.brs, predict.brs, residuals.brs, brs_predict_scoreprob, anova.brs, brsmm

Examples

# Synthetic NRS-11 pain scores (0-10, so ncuts = 10): 4 groups x 3
# post-operative times, patterned on the knee-surgery design of Lopes (2023).
# Simulated, not real data.
set.seed(2023)
nrs <- expand.grid(id = 1:80, time = c("6h", "12h", "24h"))
nrs$group <- factor(paste0("g", (nrs$id - 1) %% 4 + 1))
eta <- -1.3 + c(0, 0.75, 0.3)[nrs$time] + c(0, -0.1, 0.05, 0.1)[nrs$group]
shp <- brs_repar(mu = plogis(eta), phi = 0.3, repar = 2)  # dispersion 0.3
nrs$y <- round(10 * rbeta(nrow(nrs), shp$shape1, shp$shape2))

# Default cells [s - 0.5, s + 0.5] / 10; scores 0 and 10 are censored
fit <- brs(y ~ time + group, data = nrs, ncuts = 10)
summary(fit)
confint(fit)

# Post-fit checks: no log-likelihood left to gain, negative definite Hessian, no clamps
fit$diagnostics[c("grad_gain", "hessian_nd", "n_clamped")]

# Same scores read with right-direction cells [s, s + 1] / 11
fit_r <- brs(y ~ time + group, data = nrs, ncuts = 10, interval = "right")
cbind(mid = coef(fit), right = coef(fit_r))

# New patients: mean on (0, 1), latent score, expected score, P(S = s)
nd <- data.frame(time = c("6h", "12h", "24h"), group = "g1")
predict(fit, newdata = nd)
predict(fit, newdata = nd, type = "score")
predict(fit, newdata = nd, type = "expected_score")
round(brs_predict_scoreprob(fit, newdata = nd), 3)

# Randomized quantile residuals respect the censoring (N(0, 1) under the model)
set.seed(1)
summary(residuals(fit, type = "rqr"))
plot(fit, which = 1:2)

---

Parametric bootstrap confidence intervals for brs models

Description

Computes bootstrap-based confidence intervals for the parameters of a fitted "brs" model by repeatedly simulating data from the fitted model and re-estimating parameters. Only "brs" (fixed or variable-dispersion) objects are supported; "brsmm" is not supported.

Usage

brs_bootstrap(
  object,
  R = 199L,
  level = 0.95,
  ci_type = c("percentile", "basic", "normal", "bca"),
  max_tries = NULL,
  keep_draws = FALSE
)

## S3 method for class 'brs_bootstrap'
print(x, ...)

Arguments

object	A fitted "brs" object (fixed or variable dispersion).
R	Integer: number of bootstrap replicates (default 199).
level	Numeric: confidence level (default 0.95).
ci_type	Character: type of confidence interval. One of "percentile" (default), "basic", "normal", or "bca". See the section on the cost of "bca" below.
max_tries	Optional integer: maximum number of bootstrap attempts to obtain converged replicates. If NULL, uses max(3 * R, 50).
keep_draws	Logical: if TRUE, stores successful bootstrap parameter draws in attribute "boot_draws".
x	Object returned by brs_bootstrap.
...	Ignored.

Details

Each replicate draws a new response $y^*_i \sim \mathrm{Beta}(a_i, b_i)$ at the fitted shapes of the rows used by the fit and substitutes it into a copy of object$data; covariates, factor levels, transformations and the formula stay those of the original fit, which is then re-fitted with brs under the same ncuts, lim, interval, links, repar and optimizer. The observation mechanism of each row is reproduced: exact observations ($\delta = 0$) stay continuous; scores are re-coarsened on the fit's grid (the mapping of brs_sim). Rows whose bounds came from the analyst (brs_prep Modes 2–4, bounds that are not the cell of the row's score) keep their thresholds as fixed and non-informative (independent of $Y$), and $\delta$ is re-drawn by the cell of the partition they induce where $y^*$ falls: for an upper bound $c$, $y^* \le c$ gives $\delta = 1$ on $[\epsilon, c]$, otherwise $\delta = 2$ on $[c, 1 - \epsilon]$ (an interval $[l, u]$ gives three cells). When the original design had more thresholds than a row records (e.g. Mode 4 intervals cut from a finer instrument) the bootstrap is conservative, and for Mode 2 rows with a forced $\delta$ (a threshold built from the score itself, which is informative) it is only an approximation. The response must be a variable (not an expression such as I(y / 10)).

Each refit starts from the estimate of object and uses the compiled Hessian (hessian_method = "cpp"), so replicates are cheap.

Replicates that fail (refit error, non-convergence, non-finite estimates) are discarded and counted: attributes "n_failed" and "fail_rate", also printed. Intervals are computed from the bootstrap distribution of each parameter, with the method controlled by ci_type: "percentile" (default) uses the raw empirical quantiles; "basic" uses reflected empirical quantiles; "normal" uses a normal approximation from the bootstrap standard error (no quantiles); "bca" uses bias-corrected-and-accelerated adjusted quantiles. With parametric resampling and a nonparametric (leave-one-out) jackknife acceleration, "bca" is an approximation; a warning says so once per session, and failed jackknife refits are reported in "n_jack_failed".

Value

A data frame with columns parameter, estimate (original point estimate), se_boot (bootstrap standard error), ci_lower, ci_upper, mcse_lower, mcse_upper, wald_lower, wald_upper, and level. The attribute "n_success" gives the number of replicates that converged. Additional attributes include "R", "n_attempted", "n_failed", "fail_rate", "n_jack_failed" (BCa only), "ci_type", and optionally "boot_draws".

Methods (by generic)

• print(brs_bootstrap): Print method for bootstrap results

Cost of ci_type = "bca"

The bias-corrected and accelerated interval needs an acceleration constant, which is obtained here by a leave-one-out jackknife. That requires n additional model fits, one per observation, on top of the R bootstrap replicates: the total is R + n fits rather than R. The cost is therefore driven by the sample size, not by R, and grows quickly – for n = 1000 the jackknife alone dominates the run time by an order of magnitude. The other three interval types need only the R replicates. Prefer "percentile" or "basic" for exploratory work on large samples, and reserve "bca" for a final result.

See Also

confint.brs for Wald intervals; brs_sim for simulation; brs for fitting.

Examples

# Synthetic NRS-11 scores at 6h, 12h, 24h (time is a factor)
set.seed(3)
nrs <- data.frame(time = factor(rep(c("6h", "12h", "24h"), each = 40),
                                levels = c("6h", "12h", "24h")))
shp <- brs_repar(mu = plogis(-1.3 + c(0, 0.75, 0.3)[nrs$time]), phi = 0.3)
nrs$y <- round(10 * rbeta(nrow(nrs), shp$shape1, shp$shape2))
fit <- brs(y ~ time, data = nrs, ncuts = 10)

# Percentile intervals from 30 parametric replicates (use R >= 199 in practice)
set.seed(4)
bt <- brs_bootstrap(fit, R = 30)
bt
# Bootstrap next to Wald limits, and the bootstrap/Wald SE ratio
cols <- c("parameter", "ci_lower", "ci_upper", "wald_lower", "wald_upper")
as.data.frame(bt)[, cols]
round(bt$se_boot / sqrt(diag(vcov(fit))), 2)

---

Graphical and tabular censoring summary

Description

Produces a visual summary of the censoring structure in a fitted "brs" model or a response matrix produced by brs_check. The summary includes:

1. Bar chart of censoring-type counts, with percentage labels (censoring burden).

2. Boxplot (with jittered points) of the midpoint response $y_t$ by censoring type.

3. Boxplot of interval width ($u_i - l_i$) by censoring type.

4. Ordered interval map: $[l_i, u_i]$ segments for a sample of observations, colored by censoring type.

Usage

brs_cens(
  object,
  n_sample = 100L,
  which = 1:4,
  caption = NULL,
  gg = FALSE,
  title = "Censoring diagnostic overview",
  sub.caption = NULL,
  theme = NULL,
  palette = NULL,
  inform = FALSE,
  ...
)

Arguments

object	A fitted "brs" object, a matrix returned by brs_check, or a data frame returned by brs_prep (must contain columns left, right, yt, and delta).
n_sample	Integer: maximum number of observations to show in the interval plot (default 100). If the data has more observations, a random sample is drawn.
which	Integer vector selecting which panels to draw (default 1:4).
caption	Optional panel captions. Accepts a character vector (or list coercible to character) with up to 4 labels, in the order: burden, midpoint-by-type, width-by-type, ordered interval map.
gg	Logical: use ggplot2? (default FALSE).
title	Optional global title for the plotting page.
sub.caption	Optional subtitle/caption for the plotting page.
theme	Optional ggplot2 theme object (e.g., ggplot2::theme_bw()). If NULL, a minimal theme is used when gg = TRUE.
palette	Optional named character vector with colors for censoring types Exact, Left, Right, and Interval.
inform	Logical; if TRUE, prints brief interpretation messages about boundary and interval censoring intensity.
...	Further arguments (currently ignored).

Value

Invisibly returns a data frame with censoring counts and proportions, percentages, and interpretation flags.

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

Hawker, G. A., Mian, S., Kendzerska, T., and French, M. (2011). Measures of adult pain: Visual Analog Scale for Pain (VAS Pain), Numeric Rating Scale for Pain (NRS Pain), McGill Pain Questionnaire (MPQ), Short-Form McGill Pain Questionnaire (SF-MPQ), Chronic Pain Grade Scale (CPGS), Short Form-36 Bodily Pain Scale (SF-36 BPS), and Measure of Intermittent and Constant Osteoarthritis Pain (ICOAP). Arthritis Care and Research, 63(S11), S240-S252. doi:10.1002/acr.20543

Hjermstad, M. J., Fayers, P. M., Haugen, D. F., et al. (2011). Studies comparing Numerical Rating Scales, Verbal Rating Scales, and Visual Analogue Scales for assessment of pain intensity in adults: a systematic literature review. Journal of Pain and Symptom Management, 41(6), 1073-1093. doi:10.1016/j.jpainsymman.2010.08.016

Examples

y <- c(0, 3, 5, 7, 10)
Y <- brs_check(y, ncuts = 10)
brs_cens(Y)

prep <- brs_prep(data.frame(y = y), ncuts = 10)
brs_cens(prep)

---

Map scale scores to intervals on (0, 1) and censoring types

Description

Maps a score on $\{0, 1, \ldots, K\}$ ($K =$ ncuts, the maximum score; $K + 1$ categories) to a cell $[l_s, u_s]$ of $(0, 1)$ and to a censoring type $\delta$ of the complete likelihood (Lopes, 2023, eq. eqn_verossimilhanca_geral; see brs): $\delta = 0$ density $f(y)$, $\delta = 1$ $F(u)$, $\delta = 2$ $1 - F(l)$, $\delta = 3$ $F(u) - F(l)$. A value in $(0, 1)$ is exact ($\delta = 0$), observation by observation. A scale that starts at 1 (e.g. a Likert item 1–5) must be shifted to start at 0 (0–4, ncuts = 4): its lowest category is then the left-censored border.

Usage

brs_check(
  y,
  ncuts = 100L,
  lim = 0.5,
  delta = NULL,
  interval = c("mid", "right", "left")
)

Arguments

y	Numeric vector: scores on $\{0, 1, \ldots, K\}$ or values in $(0, 1)$.
ncuts	Integer $K$: the maximum score (default 100); the scale is $0, 1, \ldots, K$, with $K + 1$ categories. Must be $\geq \max(y)$.
lim	Numeric in $(0, 0.5]$: half-width of the cell under interval = "mid" (default 0.5, adjacent cells touch). Values below 0.5 give a partial coarsening (warning); ignored, with a warning, for "right"/"left".
delta	Integer vector or NULL. If NULL (default), censoring types are derived from the scores. If provided, it must have the length of y with elements in {0, 1, 2, 3}; it overrides the type per observation (see Details).
interval	Direction of the uncertainty interval: "mid" (default), "right" or "left"; see the section 'Interval direction'.

Details

With delta = NULL, each value in $(0, 1)$ is exact and each other value is a score with its cell and the type above; this is the rule of brs_prep too. Input that mixes values in $(0, 1)$ with values $\ge 1$ is ambiguous (proportions next to scores, or rescaled scores) and gives one warning; rescale it if the values in $(0, 1)$ are meant as scores (half-point scores 0, 0.5, 1, ...: use y * 2 and ncuts * 2). A user-supplied delta (the mechanism brs_sim uses in Monte Carlo studies) forces the type per observation and keeps the cell endpoints of the score:

$\delta$	$l_i$	$u_i$
0	cell centre (or $y$ when in $(0, 1)$)	same
1	$\epsilon$	$u_s$
2	$l_s$	$1 - \epsilon$
3	$l_s$	$u_s$

A forced $\delta = 1$ or 2 uses the cell of the observed score as the censoring threshold, so the threshold depends on $Y$: applied to every observation this is informative censoring and the likelihood has no finite maximum (brs_sim warns). Under "mid" with lim = 0.5: $u_0 = 0.5/K$, $l_K = (K - 0.5)/K$ and $[l_s, u_s] = [(s - 0.5)/K, (s + 0.5)/K]$. Scores outside $[0, K]$ are an error (as in brs_prep), and so is a $\delta = 3$ observation with $l_i \ge u_i$ (zero-probability interval). When the clamp squeezes a cell to zero width (a cell narrower than $10^{-5}$ at a border, i.e. ncuts too large or lim too small), the error says so.

yt is the cell centre: $s/K$ under "mid" (so the border scores sit at $\epsilon$ and $1 - \epsilon$) and $(s + 0.5)/(K + 1)$ otherwise; exact values are kept. It is the density argument for $\delta = 0$ and the point used by the midpoint residuals; censored contributions do not use it.

Data carrying the "is_prepared" attribute (from brs_prep or brs_sim) are used by brs as they are; brs_check() runs only for raw scores.

Value

A numeric matrix with one row per observation and columns left ($l_i$), right ($u_i$), yt (cell centre), y (the input) and delta.

Interval direction

interval is the direction of the uncertainty interval around the score (Lopes, 2023, "Mapeamento de intervalos para beta": $m = [s - 0.5, s + 0.5]$, $r = [s, s + 1]$, $l = [s - 1, s]$).

interval	cell of score $s$	latent score of $y^* \in (0, 1)$
"mid"	$[s - \mathrm{lim}, s + \mathrm{lim}]/K$	$K y^*$
"right"	$[s, s + 1]/(K + 1)$	$(K + 1) y^*$
"left"	$[s, s + 1]/(K + 1)$	$(K + 1) y^* - 1$

The latent score is the back-transformation used by predict(type = "score"): it lies in $[s - 0.5, s + 0.5]$, $[s, s + 1]$ or $[s - 1, s]$ when $y^*$ is in the cell of $s$. The $K + 1$ cells of "right" and "left" are equal and partition $[0, 1]$. This is a package choice that differs from the dissertation, which divides $r$ and $l$ by $K$; there the two directions differ by $1/K$ and chapter 4 reports opposite intercept biases for them, while here they give the same likelihood and the same coefficients. The three directions are different coarsening models of the same scores: their log-likelihoods are not comparable, and anova.brs refuses such comparisons. lim applies to "mid" only.

Scale change and borders

Lopes (2023, "Mudanca de escala intervalar") rescales scores by the range and handles the edge effect either with the transformation $y^* = \{y(n - 1)/R + 1/2\}/n$ (Smithson and Verkuilen, 2006), which depends on the sample size $n$ and the observed range $R$, or by moving 0 and 1 inwards by $\zeta = 10^{-4}$. The package does neither: each score is coarsened into its cell, using the scale maximum $K$ rather than the observed range; the border scores 0 and $K$ become left- and right-censored observations with contributions $F(u_0)$ and $1 - F(l_K)$; and the endpoints are clamped to $[\epsilon, 1 - \epsilon]$, $\epsilon = 10^{-5}$, only as a numerical guard. The border cells therefore keep their full probability, no data are moved towards 1/2, and the fit does not depend on $n$. The clamp does not enter the border contributions ($F(u_0)$ ignores $l_0$, $1 - F(l_K)$ ignores $u_K$); it matters only for exact values and for extreme user-supplied intervals.

The censoring type is derived from the score, before any clamping: $s = 0 \to \delta = 1$, $s = K \to \delta = 2$, otherwise $\delta = 3$.

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

Smithson, M., and Verkuilen, J. (2006). A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables. Psychological Methods, 11(1), 54–71. doi:10.1037/1082-989X.11.1.54

Hjermstad, M. J., Fayers, P. M., Haugen, D. F., et al. (2011). Studies comparing Numerical Rating Scales, Verbal Rating Scales, and Visual Analogue Scales for assessment of pain intensity in adults: a systematic literature review. Journal of Pain and Symptom Management, 41(6), 1073-1093. doi:10.1016/j.jpainsymman.2010.08.016

See Also

brs_prep for analyst-supplied censoring and intervals; brs; brs_sim.

Examples

# NRS-11 scores (K = 10) under the three interval directions
s <- 0:10
brs_check(s, ncuts = 10)                      # mid:   [s - 0.5, s + 0.5] / 10
brs_check(s, ncuts = 10, interval = "right")  # right: [s, s + 1] / 11
# "left" has the same cells as "right"; only the latent score read back differs
identical(brs_check(s, 10, interval = "left"), brs_check(s, 10, interval = "right"))

# Borders: score 0 -> delta 1 (F(u)), score 10 -> delta 2 (1 - F(l))
brs_check(c(0, 10), ncuts = 10)[, c("left", "right", "delta")]

# A Likert item 1-5: shift to 0-4 so that K = ncuts = 4 is its maximum
likert <- c(1, 2, 3, 5, 4)
brs_check(likert - 1, ncuts = 4)[, c("left", "right", "delta")]
# Unshifted, category 1 would be an interior cell and 0 a category never used
brs_check(likert, ncuts = 5)[, "delta"]

# Values already in (0, 1) are exact observations (delta 0)
brs_check(c(0.12, 0.5, 0.97))
# Per observation: 0.3 is exact, 5 and 10 are scores (with a warning)
brs_check(c(0.3, 5, 10), ncuts = 10)

# A forced delta keeps the cell of the score
brs_check(c(30, 60), ncuts = 100, delta = c(1L, 2L))

---

Internal coefficient table (deprecated, use brs_est() or summary())

Description

Deprecated convenience wrapper. Use brs_est for coefficient estimates or summary.brs for a full model summary.

Usage

brs_coef(fit, alpha = 0.05)

Arguments

fit	A fitted "brs" object.
alpha	Significance level.

Value

A list with components est (from brs_est) and gof (from brs_gof).

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

Ferrari, S. L. P., and Cribari-Neto, F. (2004). Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31(7), 799–815. doi:10.1080/0266476042000214501

Hawker, G. A., Mian, S., Kendzerska, T., and French, M. (2011). Measures of adult pain: Visual Analog Scale for Pain (VAS Pain), Numeric Rating Scale for Pain (NRS Pain), McGill Pain Questionnaire (MPQ), Short-Form McGill Pain Questionnaire (SF-MPQ), Chronic Pain Grade Scale (CPGS), Short Form-36 Bodily Pain Scale (SF-36 BPS), and Measure of Intermittent and Constant Osteoarthritis Pain (ICOAP). Arthritis Care and Research, 63(S11), S240-S252. doi:10.1002/acr.20543

Hjermstad, M. J., Fayers, P. M., Haugen, D. F., et al. (2011). Studies comparing Numerical Rating Scales, Verbal Rating Scales, and Visual Analogue Scales for assessment of pain intensity in adults: a systematic literature review. Journal of Pain and Symptom Management, 41(6), 1073-1093. doi:10.1016/j.jpainsymman.2010.08.016

See Also

brs_est, brs_gof, summary.brs

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10)
)
prep <- brs_prep(dat, ncuts = 100)
fit <- brs(y ~ x1, data = prep)
suppressWarnings(brs_coef(fit))  # deprecated; use brs_est()

---

K-fold cross-validation for brs models

Description

Performs repeated k-fold cross-validation for brs models.

Usage

brs_cv(formula, data, k = 5L, repeats = 1L, ...)

Arguments

formula	Model formula passed to brs.
data	Data frame.
k	Number of folds.
repeats	Number of repeated k-fold runs.
...	Additional arguments forwarded to brs (e.g., repar, link, interval, method).

Details

The log_score is the mean log predictive contribution under the complete likelihood contribution implied by each observation's censoring type (delta).

Value

A data frame with one row per fold and columns: repeat, fold, n_train, n_test, log_score, rmse_yt, mae_yt, converged, and error. The object has class "brs_cv".

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

Hawker, G. A., Mian, S., Kendzerska, T., and French, M. (2011). Measures of adult pain: Visual Analog Scale for Pain (VAS Pain), Numeric Rating Scale for Pain (NRS Pain), McGill Pain Questionnaire (MPQ), Short-Form McGill Pain Questionnaire (SF-MPQ), Chronic Pain Grade Scale (CPGS), Short Form-36 Bodily Pain Scale (SF-36 BPS), and Measure of Intermittent and Constant Osteoarthritis Pain (ICOAP). Arthritis Care and Research, 63(S11), S240-S252. doi:10.1002/acr.20543

Hjermstad, M. J., Fayers, P. M., Haugen, D. F., et al. (2011). Studies comparing Numerical Rating Scales, Verbal Rating Scales, and Visual Analogue Scales for assessment of pain intensity in adults: a systematic literature review. Journal of Pain and Symptom Management, 41(6), 1073-1093. doi:10.1016/j.jpainsymman.2010.08.016

Examples

# Synthetic NRS-11 scores: 4 groups x 3 times. Simulated, not real data.
set.seed(3)
nrs <- expand.grid(id = 1:40, time = c("6h", "12h", "24h"))
nrs$group <- factor(paste0("g", (nrs$id - 1) %% 4 + 1))
eta <- -1.3 + c(0, 0.75, 0.3)[nrs$time] + c(0, -0.1, 0.05, 0.1)[nrs$group]
shp <- brs_repar(mu = plogis(eta), phi = 0.3, repar = 2)
nrs$y <- round(10 * rbeta(nrow(nrs), shp$shape1, shp$shape2))

# 3-fold CV of two nested models, time only (m1) and time + group (m2);
# log_score is the mean held-out log-likelihood contribution (higher is better)
set.seed(5)
cv1 <- brs_cv(y ~ time, data = nrs, k = 3, ncuts = 10)
set.seed(5)
cv2 <- brs_cv(y ~ time + group, data = nrs, k = 3, ncuts = 10)
c(m1 = mean(cv1$log_score), m2 = mean(cv2$log_score))
cv2

---

Coefficient estimates with inference

Description

Coefficient estimates with inference

Usage

brs_est(object, alpha = 0.05)

Arguments

object	A fitted "brs" object.
alpha	Significance level (default 0.05).

Value

Data frame of estimates, standard errors, z-values, and p-values.

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

Ferrari, S. L. P., and Cribari-Neto, F. (2004). Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31(7), 799–815. doi:10.1080/0266476042000214501

Hawker, G. A., Mian, S., Kendzerska, T., and French, M. (2011). Measures of adult pain: Visual Analog Scale for Pain (VAS Pain), Numeric Rating Scale for Pain (NRS Pain), McGill Pain Questionnaire (MPQ), Short-Form McGill Pain Questionnaire (SF-MPQ), Chronic Pain Grade Scale (CPGS), Short Form-36 Bodily Pain Scale (SF-36 BPS), and Measure of Intermittent and Constant Osteoarthritis Pain (ICOAP). Arthritis Care and Research, 63(S11), S240-S252. doi:10.1002/acr.20543

Hjermstad, M. J., Fayers, P. M., Haugen, D. F., et al. (2011). Studies comparing Numerical Rating Scales, Verbal Rating Scales, and Visual Analogue Scales for assessment of pain intensity in adults: a systematic literature review. Journal of Pain and Symptom Management, 41(6), 1073-1093. doi:10.1016/j.jpainsymman.2010.08.016

See Also

brs, brs_gof, brs_hessian, summary.brs

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10)
)
prep <- brs_prep(dat, ncuts = 100)
fit <- brs(y ~ x1, data = prep)
brs_est(fit)

---

Goodness-of-fit measures

Description

Goodness-of-fit measures

Usage

brs_gof(object)

Arguments

object

A fitted "brs" or "brsmm" object.

Value

Data frame with logLik, AIC, BIC, and pseudo-R-squared.

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

Hawker, G. A., Mian, S., Kendzerska, T., and French, M. (2011). Measures of adult pain: Visual Analog Scale for Pain (VAS Pain), Numeric Rating Scale for Pain (NRS Pain), McGill Pain Questionnaire (MPQ), Short-Form McGill Pain Questionnaire (SF-MPQ), Chronic Pain Grade Scale (CPGS), Short Form-36 Bodily Pain Scale (SF-36 BPS), and Measure of Intermittent and Constant Osteoarthritis Pain (ICOAP). Arthritis Care and Research, 63(S11), S240-S252. doi:10.1002/acr.20543

Hjermstad, M. J., Fayers, P. M., Haugen, D. F., et al. (2011). Studies comparing Numerical Rating Scales, Verbal Rating Scales, and Visual Analogue Scales for assessment of pain intensity in adults: a systematic literature review. Journal of Pain and Symptom Management, 41(6), 1073-1093. doi:10.1016/j.jpainsymman.2010.08.016

See Also

brs, brs_est, brs_hessian

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10)
)
prep <- brs_prep(dat, ncuts = 100)
fit <- brs(y ~ x1, data = prep)
brs_gof(fit)

---

Extract the Hessian matrix

Description

Extract the Hessian matrix

Usage

brs_hessian(object)

Arguments

object

A fitted "brs" or "brsmm" object.

Value

Numeric Hessian matrix.

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

Hawker, G. A., Mian, S., Kendzerska, T., and French, M. (2011). Measures of adult pain: Visual Analog Scale for Pain (VAS Pain), Numeric Rating Scale for Pain (NRS Pain), McGill Pain Questionnaire (MPQ), Short-Form McGill Pain Questionnaire (SF-MPQ), Chronic Pain Grade Scale (CPGS), Short Form-36 Bodily Pain Scale (SF-36 BPS), and Measure of Intermittent and Constant Osteoarthritis Pain (ICOAP). Arthritis Care and Research, 63(S11), S240-S252. doi:10.1002/acr.20543

Hjermstad, M. J., Fayers, P. M., Haugen, D. F., et al. (2011). Studies comparing Numerical Rating Scales, Verbal Rating Scales, and Visual Analogue Scales for assessment of pain intensity in adults: a systematic literature review. Journal of Pain and Symptom Management, 41(6), 1073-1093. doi:10.1016/j.jpainsymman.2010.08.016

See Also

brs, vcov.brs, brs_est

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10)
)
prep <- brs_prep(dat, ncuts = 100)
fit <- brs(y ~ x1, data = prep)
brs_hessian(fit)

---

Marginal effects for brs models

Description

Computes average marginal effects (AME) for numeric covariates in the mean or precision submodel of a fitted "brs" object.

Usage

brs_marginaleffects(
  object,
  newdata = NULL,
  model = c("mean", "precision"),
  type = c("response", "link"),
  variables = NULL,
  h = 1e-05,
  interval = TRUE,
  level = 0.95,
  n_sim = 400L,
  keep_draws = FALSE
)

Arguments

object	A fitted "brs" object.
newdata	Optional data frame for evaluation; defaults to the data used in fitting.
model	Character; "mean" (default) or "precision".
type	Character prediction scale: "response" (default) or "link".
variables	Optional character vector of covariate names. Defaults to all numeric covariates in the selected submodel.
h	Finite-difference step for non-binary numeric covariates.
interval	Logical; compute interval estimates via simulation.
level	Confidence level for interval estimates.
n_sim	Number of parameter draws when interval = TRUE.
keep_draws	Logical; if TRUE and interval = TRUE, stores AME simulation draws in attribute "ame_draws".

Details

AMEs for a numeric covariate are computed by a central difference on predictions, with the step size scaled by the covariate's standard deviation:

$\mathrm{AME}_j = \frac{1}{n}\sum_{i=1}^{n} \frac{\hat{g}_i(x_{ij} + h_j) - \hat{g}_i(x_{ij} - h_j)}{2 h_j}, \qquad h_j = h \cdot \max(\mathrm{sd}(x_{\cdot j}), 1),$

where $\hat{g}_i$ is the selected prediction scale and h is the (small) base step supplied by the caller.

For binary covariates coded as 0/1, the effect is computed as the average discrete difference $\hat{g}(x_j=1)-\hat{g}(x_j=0)$.

If interval = TRUE, uncertainty is approximated by asymptotic parameter simulation from $\mathcal{N}(\hat{\theta}, \hat{V})$.

Value

A data frame with one row per variable and columns: variable, ame, std.error, ci.lower, ci.upper, model, type, and n. The returned object has class "brs_marginaleffects" and attributes with analysis metadata.

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

Hawker, G. A., Mian, S., Kendzerska, T., and French, M. (2011). Measures of adult pain: Visual Analog Scale for Pain (VAS Pain), Numeric Rating Scale for Pain (NRS Pain), McGill Pain Questionnaire (MPQ), Short-Form McGill Pain Questionnaire (SF-MPQ), Chronic Pain Grade Scale (CPGS), Short Form-36 Bodily Pain Scale (SF-36 BPS), and Measure of Intermittent and Constant Osteoarthritis Pain (ICOAP). Arthritis Care and Research, 63(S11), S240-S252. doi:10.1002/acr.20543

Hjermstad, M. J., Fayers, P. M., Haugen, D. F., et al. (2011). Studies comparing Numerical Rating Scales, Verbal Rating Scales, and Visual Analogue Scales for assessment of pain intensity in adults: a systematic literature review. Journal of Pain and Symptom Management, 41(6), 1073-1093. doi:10.1016/j.jpainsymman.2010.08.016

Examples

# Synthetic NRS-11 scores with a 0/1 time dummy and a numeric covariate (age)
set.seed(3)
nrs <- data.frame(t12 = rep(0:1, each = 60), age = round(rnorm(120, 40, 12)))
shp <- brs_repar(mu = plogis(-1.3 + 0.75 * nrs$t12 + 0.03 * (nrs$age - 40)),
                 phi = 0.3, repar = 2)
nrs$y <- round(10 * rbeta(nrow(nrs), shp$shape1, shp$shape2))
fit <- brs(y ~ t12 + age, data = nrs, ncuts = 10)

# Average marginal effects on E[Y] (0/1 dummy: discrete change; age: slope),
# with intervals from 50 draws of the estimates (use n_sim >= 400 in practice)
set.seed(6)
brs_marginaleffects(fit, n_sim = 50)

# On the 0-10 score scale (mid cells): multiply by K = 10
set.seed(6)
ame <- brs_marginaleffects(fit, n_sim = 50)
transform(ame[, c("variable", "ame", "ci.lower", "ci.upper")],
          ame = 10 * ame, ci.lower = 10 * ci.lower, ci.upper = 10 * ci.upper)

---

Predict score probabilities from a fitted brs model

Description

Computes predicted probabilities for integer scores on the original scale $\{0, 1, \ldots, K\}$ implied by the fitted beta interval model.

Usage

brs_predict_scoreprob(
  object,
  newdata = NULL,
  scores = NULL,
  format = c("matrix", "long"),
  id_col = "id"
)

Arguments

object	A fitted "brs" object.
newdata	Optional data frame for prediction.
scores	Optional integer vector of scores to evaluate. Defaults to all scores in 0:ncuts.
format	Output format: "matrix" (default) or "long".
id_col	Column name for observation id when format = "long".

Details

For a score $s$ with cell $[l_s, u_s]$ of the fit's interval (see brs_check), $P(S = s) = F(u_s) - F(l_s)$ with $l_0 = 0$ and $u_K = 1$, where $F$ is the beta CDF under the fitted $(\mu_i, \phi_i)$. Under "mid" this is $F(\mathrm{lim}/K)$ for $s = 0$, $1 - F((K - \mathrm{lim})/K)$ for $s = K$ and $F((s + \mathrm{lim})/K) - F((s - \mathrm{lim})/K)$ otherwise; under "right"/"left" the cells are $[s, s + 1]/(K + 1)$. The rows sum to 1 whenever the cells partition $[0, 1]$ (always for "right"/"left"; lim = 0.5 for "mid").

Value

If format = "matrix", a numeric matrix with one row per observation and one column per requested score.

If format = "long", a long data frame with columns id_col, score, and prob.

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

Hawker, G. A., Mian, S., Kendzerska, T., and French, M. (2011). Measures of adult pain: Visual Analog Scale for Pain (VAS Pain), Numeric Rating Scale for Pain (NRS Pain), McGill Pain Questionnaire (MPQ), Short-Form McGill Pain Questionnaire (SF-MPQ), Chronic Pain Grade Scale (CPGS), Short Form-36 Bodily Pain Scale (SF-36 BPS), and Measure of Intermittent and Constant Osteoarthritis Pain (ICOAP). Arthritis Care and Research, 63(S11), S240-S252. doi:10.1002/acr.20543

Hjermstad, M. J., Fayers, P. M., Haugen, D. F., et al. (2011). Studies comparing Numerical Rating Scales, Verbal Rating Scales, and Visual Analogue Scales for assessment of pain intensity in adults: a systematic literature review. Journal of Pain and Symptom Management, 41(6), 1073-1093. doi:10.1016/j.jpainsymman.2010.08.016

Examples

dat <- data.frame(
  y = c(
    0, 5, 20, 50, 75, 90, 100, 30, 60, 45,
    10, 40, 55, 70, 85, 25, 35, 65, 80, 15
  ),
  x1 = rep(c(1, 2), 10),
  x2 = rep(c(0, 0, 1, 1), 5)
)
prep <- brs_prep(dat, ncuts = 100)
fit <- brs(y ~ x1, data = prep)
pmat <- brs_predict_scoreprob(fit)
head(pmat[, 1:5])
plong <- brs_predict_scoreprob(fit, scores = 0:10, format = "long")
head(plong)

---

Pre-process analyst data for beta interval regression

Description

Turns analyst data into the cells $[l_i, u_i]$ and censoring types $\delta_i$ used by brs and brsmm. Scores run over $0, 1, \ldots, K$ with $K =$ ncuts the maximum; shift a scale that starts at 1 (a Likert item 1–5 becomes 0–4, ncuts = 4). Four input modes are recognised per row:

1. Score only (y): the cell of the score and $\delta$ from the score, exactly as brs_check (0 $\to$ 1, $K$ $\to$ 2, otherwise 3; a value in $(0, 1)$ $\to$ 0, per observation).

2. Score and delta: the analyst's censoring type, with the cell of the score.

3. Bounds only (left and/or right, y missing): an interval, or a one-sided censoring when one bound is NA.

4. Score and both bounds: the analyst's interval.

Covariate columns are kept unchanged.

Usage

brs_prep(
  data,
  y = "y",
  delta = "delta",
  left = "left",
  right = "right",
  ncuts = 100L,
  lim = 0.5,
  interval = c("mid", "right", "left")
)

Arguments

data	Data frame with the response columns and covariates.
y	Character: name of the score column (default "y").
delta	Character: name of the censoring-type column (default "delta"); values in {0, 1, 2, 3} or NA.
left, right	Character: names of the lower and upper bound columns (defaults "left", "right"), on the latent score scale.
ncuts	Integer $K$: the maximum score (default 100); the scale is $0, 1, \ldots, K$ ($K + 1$ categories). Must be at least the largest y; bounds may reach the latent range given in Details.
lim	Numeric in $(0, 0.5]$: half-width of the score cell under interval = "mid" (default 0.5); see brs_check.
interval	Direction of the uncertainty interval, "mid" (default), "right" or "left"; see brs_check.

Details

A non-missing delta always wins. Otherwise the type comes from the pattern of left, right and y (here $K = 10$, interval = "mid"):

left	right	y	cell on $(0, 1)$	$\delta$
NA	3	NA	$[\epsilon, 0.3]$ (below 3)	1
7	NA	NA	$[0.7, 1 - \epsilon]$ (above 7)	2
2	5	NA	$[0.2, 0.5]$	3
4	6	5	$[0.4, 0.6]$ (analyst interval)	3
-0.5	4	NA	$[\epsilon, 0.4]$ (reaches 0)	1
3	10.5	NA	$[0.3, 1 - \epsilon]$ (reaches 1)	2
(no columns)		5	$[0.45, 0.55]$ (Mode 1)	3
NA	NA	5	0.5 (exact reading)	0
NA	NA	0	$[\epsilon, 0.05]$	1
NA	NA	10	$[0.95, 1 - \epsilon]$	2

When the data have left/right columns, a row that gives only an interior score (both bounds NA) is an exact value at the cell centre, with a density contribution; without those columns the same score is interval-censored (Mode 1). An analyst interval that reaches 0 (or 1) on the unit scale is left- (or right-) censored; one that reaches both borders covers the whole scale, stays $\delta = 3$ and gives a warning, since such a row carries no information about the parameters. A row with neither a score nor a bound is an error (a delta alone defines no interval).

Score-based rows (Modes 1 and 2) use the cells of interval, as in brs_check: $[s - \mathrm{lim}, s + \mathrm{lim}]/K$ for "mid", $[s, s + 1]/(K + 1)$ for "right"/"left"; a forced $\delta = 1$ keeps $u_s$ and a forced $\delta = 2$ keeps $l_s$. Analyst bounds $L$ (Modes 3 and 4) are values of the latent score of the chosen direction (the scale of predict(type = "score")) and map to $L/K$ ("mid"), $L/(K + 1)$ ("right") or $(L + 1)/(K + 1)$ ("left"), so that the dissertation's intervals $[s - 0.5, s + 0.5]$, $[s, s + 1]$ and $[s - 1, s]$ all give the cell of score $s$. They must lie in the latent range of the direction: $[-0.5, K + 0.5]$ ("mid"), $[0, K + 1]$ ("right") or $[-1, K]$ ("left"); a bound outside it is an error. Rows with no score get y = the latent score of their cell centre, so that model.frame() keeps them.

Endpoints are clamped to $[\epsilon, 1 - \epsilon]$, $\epsilon = 10^{-5}$ (why this replaces the edge-effect transformation of Lopes, 2023: section 'Scale change and borders' of brs_check). A $\delta = 3$ row with left >= right is an error (zero probability; use $\delta = 0$), and so is a cell the clamp squeezes to zero width (the message tells whether ncuts is too large for the cells or the analyst interval lies inside the clamp). As in brs_check, values in $(0, 1)$ mixed with values $\ge 1$ give one warning. Unusual combinations, e.g. $\delta = 1$ with $y \neq 0$, give a warning but are kept.

Value

A data frame with columns left, right (cell on $(0, 1)$), yt (cell centre), y (score, or the filled latent score for rows without one) and delta, followed by the covariates. Attributes "is_prepared" (TRUE), "ncuts", "lim" and "interval" are reused by brs and brsmm; an explicit different value there is ignored with a warning.

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://hdl.handle.net/1884/86624.

See Also

brs_check for the cells and the automatic classification; brs for fitting the model.

Examples

# Mode 1: score only; delta from the score (0 -> 1, K -> 2, else 3)
d1 <- data.frame(y = c(0, 3, 10), x = c(1.2, 0.4, -0.3))
brs_prep(d1, ncuts = 10)

# Mode 2: score + analyst delta (the same score read as exact and as a cell)
d2 <- data.frame(y = c(4, 4), delta = c(0, 3))
brs_prep(d2, ncuts = 10)

# Mode 3: only left and/or right bounds (NA pattern gives delta)
d3 <- data.frame(left = c(NA, 7, 2), right = c(3, NA, 5))
brs_prep(d3, ncuts = 10)

# Mode 4: score with analyst bounds (used as given, divided by K)
d4 <- data.frame(y = 5, left = 4, right = 6)
brs_prep(d4, ncuts = 10)

# An analyst interval reaching a border is one-sided censoring
brs_prep(data.frame(left = c(-0.5, 3), right = c(4, 10.5)), ncuts = 10)

# A Likert item 1-5: shift it to 0-4, so that ncuts = 4 is the maximum
lk <- data.frame(item = c(1, 2, 5, 3, 4), x = c(0.1, 0.5, 0.9, 0.3, 0.7))
brs_prep(data.frame(y = lk$item - 1, x = lk$x), ncuts = 4)

# Right-direction cells [s, s + 1] / 11; the choice is stored as attributes
p <- brs_prep(d1, ncuts = 10, interval = "right")
p[, c("left", "right", "delta")]
attributes(p)[c("ncuts", "lim", "interval")]

---