Package 'betaregscale'

Akaike information criterion

Description

$-2\ell + k\,p$, with $p$ the number of estimated parameters.

Usage

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

Arguments

object	A fitted "brs" object.
...	Further fitted models (brs, brsmm or others with an AIC() method).
k	Penalty per parameter (default 2).

Value

Scalar AIC value; with several models, a data frame with columns df and AIC, one row per model, as stats::AIC().

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)
fit0 <- brs(y ~ 1, data = prep)
AIC(fit0, fit)

---

AIC for brsmm models

Description

$-2\ell + k\,p$, with $p$ the number of parameters (fixed effects, precision, packed random-effect and cutpoint parameters).

Usage

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

Arguments

object	A fitted "brsmm" object.
...	Further fitted models (brsmm, brs or others with an AIC() method).
k	Numeric penalty per parameter (default 2).

Value

Numeric scalar; with several models, a data frame with columns df and AIC, one row per model, as stats::AIC().

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)
AIC(brs(y ~ x1, data = prep), 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: fits under different interval maps or lim are not nested, and different ncuts describe different responses. Fits of the same scores under different maps (with lim = 0.5) can still be compared by AIC() and BIC(), since each map gives a probability mass function over the same scores (section 'Interval direction' of brs).

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.

A fit with cutpoints = "fixed" is nested in the same model with cutpoints = "boundary" at an interior point of the cutpoint parameters (section 'Boundary cutpoints' of brs), so their comparison is a regular $\chi^2_2$ test when both cutpoints are identified ($K \ge 4$ in practice; with $K = 3$ the fit warns that they are weakly identified); the heading names those rows.

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://acervodigital.ufpr.br/xmlui/handle/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)

# Interval maps are not nested: no LR test, but AIC is comparable
m2_right <- brs(y ~ time + group, data = nrs, ncuts = 10, interval = "right")
try(anova(m2, m2_right))
c(mid = AIC(m2), right = AIC(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://acervodigital.ufpr.br/xmlui/handle/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, ..., n = c("obs", "groups"))

Arguments

object	A fitted "brs" object.
...	Further fitted models (brs, brsmm or others with a BIC() method).
n	"obs" (default): the number of observations. "groups" applies to brsmm fits only (see BIC.brsmm) and is an error here.

Value

Scalar BIC value; with several models, a data frame with columns df and BIC, one row per model, as stats::BIC().

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 $k$ the number of parameters (fixed effects, precision, packed random-effect and cutpoint parameters) and $n$ the number of observations (n = "obs", the default, as lme4 does) or of groups (n = "groups"). There is no single sample size for a mixed model; counting groups penalizes less when groups have several observations each. Compare BIC values only between fits with the same convention.

Usage

## S3 method for class 'brsmm'
BIC(object, ..., n = c("obs", "groups"))

Arguments

object	A fitted "brsmm" object.
...	Further fitted models; with n = "groups" they must be "brsmm" fits, otherwise their own BIC() is used.
n	"obs" (default; nobs(object)) or "groups" (object$ngroups).

Value

Numeric scalar; with several models, a data frame with columns df and BIC, one row per model, as stats::BIC().

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)
BIC(fit, n = "groups")

---

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(),
  cutpoints = c("fixed", "boundary")
)

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).
cutpoints	"fixed" (default) keeps the cutpoints of the interval map; "boundary" also estimates the two boundary cutpoints (see the section 'Boundary cutpoints' of brs).

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 (maximized 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; with cutpoints = "boundary" the cells and cell centers of the estimated cutpoints), ncuts, lim, interval, repar, link, link_phi, cutpoints (the type) and cut (NULL under the fixed map; otherwise a list with par, the logit-scale estimates, value, the cutpoints on $(0, 1)$, jacobian, fixed, the map's cutpoints, and start, their logit-scale values).

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) maximizes 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 generalized 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 center 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 (or, when the Hessian is not negative definite, the largest absolute gradient exceeds 0.1). 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 normalization). "right" and "left" give the same fit and differ only in predict(type = "score"); "mid" and "right"/"left" are different coarsening models. For score data (with lim = 0.5) each map gives a probability mass function over the same scores $0, \ldots, K$, so logLik(), AIC and BIC of fits under different maps are comparable; the maps are not nested, so there is no likelihood-ratio test between them (anova.brs stops). Details: brs_check.

Boundary cutpoints

Under the fixed map score $s$ has the cell $(c_{s-1}, c_s)$ with cutpoints fixed by interval, so the model is an ordinal threshold model whose thresholds are the scale labels. cutpoints = "boundary" estimates the two border cutpoints: $c_{lo}$ between scores 0 and 1 and $c_{hi}$ between $K - 1$ and $K$; all other cutpoints keep the values of the map. Score 0 then has the cell $(0, c_{lo})$, score 1 $(c_{lo}, u_1)$, score $K - 1$ $(l_{K-1}, c_{hi})$ and score $K$ $(c_{hi}, 1)$, where $u_1$ and $l_{K-1}$ are the ends given by the map. With $\epsilon = 10^{-5}$, the border clamp of the likelihood, the parameters are unconstrained and keep the order:

$c_{lo} = \epsilon + (u_1 - \epsilon)\, \mathrm{logit}^{-1}(g_{lo}), \qquad c_{hi} = l_{K-1} + (1 - \epsilon - l_{K-1})\, \mathrm{logit}^{-1}(g_{hi}),$

so the free cutpoints never enter the clamp and the likelihood is smooth in $(g_{lo}, g_{hi})$. The map is the interior point $g_{lo} = \mathrm{logit}\{(u_0 - \epsilon)/(u_1 - \epsilon)\}$, $g_{hi} = \mathrm{logit}\{(l_K - l_{K-1})/(1 - \epsilon - l_{K-1})\}$ (about $\mathrm{logit}(1/3)$ and $\mathrm{logit}(2/3)$ under "mid", about 0 under "right"/"left"), so the fixed model is nested in the boundary model. When both cutpoints are identified ($K \ge 4$ in practice, see below) anova.brs gives a regular likelihood-ratio test on 2 degrees of freedom. The fit starts at the map.

g_lo and g_hi are the last two coefficients, "(cut_lo)" and "(cut_hi)" (coef(fit, model = "cutpoints")); summary() reports $c_{lo}$ and $c_{hi}$ on $(0, 1)$ and on the latent score scale of interval (the scale of predict(type = "score")), with delta-method standard errors. Their gradient is analytic: $\pm f(c)/\{F(u_i) - F(l_i)\}$ for each row whose cell ends at $c$ (plus when $c$ is the upper end), times the Jacobian of the map above. The compiled Hessian covers the regression coefficients; the rows and columns of the cutpoint parameters come from numDeriv::jacobian() of the gradient (hessian_method = "cpp"). The fitted cells (fit$Y), score probabilities, expected scores, randomized quantile residuals, brs_bootstrap and brs_cv use the estimated cutpoints; brs_sim simulates from given ones.

The extension needs score data: every row the cell of a score $0, \ldots, K$ (no exact values, no analyst-supplied censoring or bounds from brs_prep), lim = 0.5, $K \ge 3$, and the scores next to each free cutpoint observed (0, 1, $K - 1$, $K$); otherwise brs() stops, since a cutpoint with an unobserved neighbor is not identified. With $K = 2$ the two cutpoints are not identified and brs() stops; with $K = 3$ they rest mainly on the shape of the beta distribution and a warning says they are weakly identified, so the $\chi^2_2$ reference is not trustworthy there. A border cell holding less than $10^{-4}$ of the range its cutpoint can take ($c_{lo} - \epsilon$ or $u_1 - c_{lo}$ relative to $u_1 - \epsilon$, and the same for $c_{hi}$) warns: that cell is collapsing and its cutpoint is not identified. The cell centers in fit$Y[, "yt"] (midpoint residuals, pseudo $R^2$) are those of the estimated cells: $c_{lo}/2$, $(c_{lo} + u_1)/2$, $(l_{K-1} + c_{hi})/2$ and $(c_{hi} + 1)/2$. The border centers differ by construction from those of the fixed map ($10^{-5}$ and $1 - 10^{-5}$ under "mid"), so pseudo $R^2$ and midpoint residuals are not directly comparable between fixed-map and boundary fits.

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://acervodigital.ufpr.br/xmlui/handle/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)

# Estimated boundary cutpoints: a 2-df likelihood-ratio test against the map
fit_b <- brs(y ~ time + group, data = nrs, ncuts = 10, cutpoints = "boundary")
summary(fit_b)$cutpoints
anova(fit, fit_b)

# 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)). A fit with cutpoints = "boundary" re-coarsens with its estimated boundary cutpoints and re-estimates them in every replicate; a replicate without one of the scores 0, 1, $K - 1$ or $K$ cannot be refitted and counts as failed. The reasons of failed replicates are kept in attr(, "fail_reasons") and printed.

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), "fail_reasons" (counts of the failure messages, or NULL), "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://acervodigital.ufpr.br/xmlui/handle/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 center (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 center: $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 center), 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. With lim = 0.5 each one gives a probability mass function over $0, \ldots, K$, so their log-likelihoods, AIC and BIC are comparable; the models are not nested, so anova.brs refuses a likelihood-ratio test between them. 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://acervodigital.ufpr.br/xmlui/handle/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://acervodigital.ufpr.br/xmlui/handle/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 (default 5).
repeats	Number of repeated k-fold runs (default 1); each repeat draws new folds with the current RNG state.
...	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).

With cutpoints = "boundary" in ..., held-out rows use the cells of the training fit's cutpoints. A training fold without one of the scores 0, 1, $K - 1$ or $K$ cannot be fitted: its row has converged = FALSE and the message in error.

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://acervodigital.ufpr.br/xmlui/handle/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

brsmm_cv for mixed models (folds by group).

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, p-values and Wald intervals. The cutpoint parameters of a cutpoints = "boundary" fit get NA z- and p-values (their zero is not a hypothesis of interest; see summary() and anova.brs).

References

Lopes, J. E. (2023). Modelos de regressao beta para dados de escala. Master's dissertation, Universidade Federal do Parana, Curitiba. URI: https://acervodigital.ufpr.br/xmlui/handle/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://acervodigital.ufpr.br/xmlui/handle/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://acervodigital.ufpr.br/xmlui/handle/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)

---

Log-likelihood contributions per observation or per group

Description

The terms whose sum is the log-likelihood: one per observation for a brs fit, and one per group (cluster) for a brsmm fit, the marginal log-likelihood of the group's observations with its random effects integrated out. With newdata, the contributions of new observations at the fitted parameters (held-out scoring).

Usage

brs_loglik_contrib(object, newdata = NULL, int_method = NULL, n_points = NULL)

Arguments

object	A fitted "brs" or "brsmm" object.
newdata	Optional data frame of new observations (see Details).
int_method, n_points	For brsmm fits: integration rule ("laplace", "aghq" or "qmc") and number of nodes (nodes per dimension for "aghq", points for "qmc"); NULL (default) uses those of the fit. Ignored for brs fits.

Details

For a brs fit the contribution of a score row is $\log P(S = s_i)$, the log probability of its cell (the log density for exact observations, $\delta = 0$), computed as in the compiled likelihood. For a brsmm fit the contribution of group $g$ is

$\log \int \prod_i P(S = s_{gi} \mid b)\, \varphi(b; 0, D)\, db,$

computed by the fit's int_method and n_points (or qmc_points) unless int_method / n_points are given. With the defaults and newdata = NULL, the contributions add up to logLik(object) (under "laplace", up to the tolerance of the inner mode search, about $10^{-8}$). A group without a positive-definite random-effect mode, which adds a penalty value to the fitted log-likelihood, gets NA with a warning.

With newdata, rows need the response, the covariates and (for brsmm) the grouping variable; rows with missing values in them are dropped. The groups need not be groups of the fit: their random effects are integrated out, so no random-effect estimate is used. Scores are coarsened with the cells of the fit (its estimated cutpoints under cutpoints = "boundary"). For score data the contributions are log probabilities of the observed scores (no Jacobian term), on the same scale as the log-likelihood of any model that gives a probability mass function on $0, \ldots, K$.

Value

Named numeric vector: by row name of the data (brs) or by group level (brsmm, with attribute "n", the number of observations of each group).

See Also

logLik.brsmm, brsmm_cv, brs_marginal

Examples

set.seed(21)
nrs <- expand.grid(id = 1:30, time = c("6h", "12h", "24h"))
b0 <- rnorm(30, sd = 0.8)
eta <- -1.3 + c(0, 0.75, 0.3)[nrs$time] + b0[nrs$id]
shp <- brs_repar(mu = plogis(eta), phi = 0.2, repar = 2)
nrs$y <- round(10 * rbeta(nrow(nrs), shp$shape1, shp$shape2))
m <- brsmm(y ~ time, random = ~ 1 | id, data = nrs, ncuts = 10)

ll_g <- brs_loglik_contrib(m)
head(ll_g)
c(sum = sum(ll_g), logLik = as.numeric(logLik(m)))
# Same groups, adaptive Gauss-Hermite with 31 nodes
head(brs_loglik_contrib(m, int_method = "aghq", n_points = 31))

---