Introduction to betaregscale

Overview

The betaregscale package provides maximum-likelihood estimation of beta regression models for responses derived from bounded rating scales. Common examples include pain intensity scales (NRS-11, NRS-21, NRS-101), Likert-type scales, product quality ratings, and any instrument whose response can be mapped to the open interval \((0,1)\).

The key idea is that a discrete score recorded on a bounded scale carries measurement uncertainty inherent to the instrument. For instance, a pain score of \(y=6\) on a 0–10 NRS is not an exact value but rather represents a range: after rescaling to \((0,1)\), the observation is treated as interval-censored in \([0.55,0.65]\). The package uses the beta distribution to model such data, building a complete likelihood that supports mixed censoring types within the same dataset (Lopes, 2023).

Installation

# Development version from GitHub:
# install.packages("remotes")
remotes::install_github("evandeilton/betaregscale")
library(betaregscale)

The score scale and ncuts

ncuts is \(K\), the maximum score. Scores run over \(0, 1, \ldots, K\), so the scale has \(K + 1\) categories: an NRS-11 pain scale (0 to 10) has ncuts = 10, a 0 to 100 scale ncuts = 100. The lowest score is the left border and the highest the right border of the scale.

A scale that starts at 1 must be shifted to start at 0. A Likert item 1 to 5 becomes 0 to 4, with ncuts = 4:

likert <- c(1, 2, 3, 5, 4, 1)
kbl10(brs_check(likert - 1, ncuts = 4))
left right yt y delta
0.000 0.125 0.00 0 1
0.125 0.375 0.25 1 3
0.375 0.625 0.50 2 3
0.875 1.000 1.00 4 2
0.625 0.875 0.75 3 3
0.000 0.125 0.00 0 1

Without the shift, category 1 would be read as an interior score and the lower border of the scale would be a category that nobody can choose:

brs_check(likert, ncuts = 5)[, "delta"]
#> [1] 3 3 3 2 3 3

Likelihood and censoring types

The complete likelihood supports four censoring types, automatically classified by brs_check():

\(\delta\) Type Likelihood contribution
0 Exact (uncensored) \(f(y_i;a_i,b_i)\)
1 Left-censored (\(y=0\)) \(F(u_i;a_i,b_i)\)
2 Right-censored (\(y=K\)) \(1-F(l_i;a_i,b_i)\)
3 Interval-censored \(F(u_i;a_i,b_i)-F(l_i;a_i,b_i)\)

where \(f(\cdot)\) and \(F(\cdot)\) are the beta density and CDF, \([l_i,u_i]\) are the interval endpoints, and \((a_i,b_i)\) are the beta shape parameters derived from \(\mu_i\) and \(\phi_i\) via the chosen reparameterization. The likelihood is the product

\[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)}.\]

This is the complete likelihood of Lopes (2023). The table of censoring types in the dissertation swaps the labels \(\delta=1\) and \(\delta=2\); the package follows its equation. The cells are fixed by the scale, not by the observed value, so the censoring is non-informative.

Interval construction

A score \(s\) is mapped to a cell of \((0,1)\). The direction of the uncertainty interval is set by interval (the cases \(m\), \(r\) and \(l\) of Lopes, 2023):

interval cell of score \(s\) latent score of \(y^*\in(0,1)\)
"mid" (default) \([s-h,\,s+h]/K\) \(Ky^*\)
"right" \([s,\,s+1]/(K+1)\) \((K+1)y^*\)
"left" \([s,\,s+1]/(K+1)\) \((K+1)y^*-1\)

\(K\) is the maximum score (ncuts) and \(h\) the half-width (lim, default 0.5, used by "mid" only). The \(K+1\) cells of "right" and "left" are equal and partition \([0,1]\); this normalization is a package choice (the dissertation divides by \(K\)). "right" and "left" give the same fit and differ only in predict(type = "score"). The three directions are different coarsening models of the same scores. With lim = 0.5 each one gives a probability mass function over \(0,\dots,K\), so log-likelihoods, AIC and BIC of fits under different directions are comparable; the directions are not nested, so anova() refuses a likelihood-ratio test between them.

The borders need no special transformation: score 0 is left-censored and score \(K\) right-censored, and endpoints are only clamped to \([10^{-5},1-10^{-5}]\) as a numerical guard. The fit therefore does not depend on the sample size, unlike the \(\{y(n-1)/R+1/2\}/n\) transformation of Smithson and Verkuilen (2006).

# Illustrate brs_check with a 0-10 NRS scale
y_example <- c(0, 3, 5, 7, 10)
cr <- brs_check(y_example, ncuts = 10)
kbl10(cr)
left right yt y delta
0.00 0.05 0.0 0 1
0.25 0.35 0.3 3 3
0.45 0.55 0.5 5 3
0.65 0.75 0.7 7 3
0.95 1.00 1.0 10 2

The delta column shows that \(y=0\) is left-censored (\(\delta=1\)), \(y=10\) is right-censored (\(\delta=2\)), and all interior values are interval-censored (\(\delta=3\)). The same scores with right-direction cells:

kbl10(brs_check(y_example, ncuts = 10, interval = "right"))
left right yt y delta
0.0000 0.0909 0.0455 0 1
0.2727 0.3636 0.3182 3 3
0.4545 0.5455 0.5000 5 3
0.6364 0.7273 0.6818 7 3
0.9091 1.0000 0.9545 10 2

Values already in \((0,1)\) are exact observations (\(\delta=0\)), one observation at a time. Input that mixes such values with scores \(\geq 1\) is ambiguous and gives a warning:

kbl10(brs_check(c(0.3, 5, 10), ncuts = 10))
#> Warning: The response mixes values in (0, 1) with values >= 1: values in (0, 1)
#> are taken as exact (delta = 0), the others as scores on 0..ncuts. Rescale the
#> data if the (0, 1) values are scores (half-point scores: use y * 2 and ncuts *
#> 2).
left right yt y delta
0.30 0.30 0.3 0.3 0
0.45 0.55 0.5 5.0 3
0.95 1.00 1.0 10.0 2

Data preparation with brs_prep()

In practice, analysts may want to supply their own censoring indicators or interval endpoints rather than relying on the automatic classification of brs_check(). The brs_prep() function provides a flexible, validated bridge between raw analyst data and brs().

It supports four input modes:

Mode 1: Score only (automatic)

# Equivalent to brs_check - delta inferred from y
d1 <- data.frame(y = c(0, 3, 5, 7, 10), x1 = rnorm(5))
kbl10(brs_prep(d1, ncuts = 10))
left right yt y delta x1
0.00 0.05 0.0 0 1 0.5769
0.25 0.35 0.3 3 3 1.1427
0.45 0.55 0.5 5 3 -0.7074
0.65 0.75 0.7 7 3 -0.6186
0.95 1.00 1.0 10 2 -0.0083

Mode 2: Score + explicit censoring indicator

# Analyst specifies delta directly
d2 <- data.frame(
  y     = c(50, 0, 99, 50),
  delta = c(0, 1, 2, 3),
  x1    = rnorm(4)
)
kbl10(brs_prep(d2, ncuts = 100))
left right yt y delta x1
0.500 0.500 0.50 50 0 0.3389
0.000 0.005 0.00 0 1 1.4054
0.985 1.000 0.99 99 2 -0.8653
0.495 0.505 0.50 50 3 -0.8740

Mode 3: Interval endpoints with NA patterns

When the analyst provides left and/or right columns, censoring is inferred from the NA pattern. Bounds are on the latent score scale and must lie in \([-0.5, K+0.5]\) under "mid" (\([0,K+1]\) and \([-1,K]\) under "right" and "left"). An interval that reaches 0 or 1 on the unit scale becomes left- or right-censored; one that reaches both covers the whole scale, carries no information and gives a warning.

d3 <- data.frame(
  left  = c(NA, 20, 30, NA),
  right = c(5, NA, 45, NA),
  y     = c(NA, NA, NA, 50),
  x1    = rnorm(4)
)
kbl10(brs_prep(d3, ncuts = 100))
left right yt y delta x1
0.0 0.05 0.025 2.5 1 0.7743
0.2 1.00 0.600 60.0 2 -0.4019
0.3 0.45 0.375 37.5 3 -0.2151
0.5 0.50 0.500 50.0 0 0.6053

Mode 4: Analyst-supplied intervals

When the analyst provides y, left, and right simultaneously, their endpoints are used directly (rescaled by \(K\) under "mid"):

d4 <- data.frame(
  y     = c(50, 75),
  left  = c(48, 73),
  right = c(52, 77),
  x1    = rnorm(2)
)
kbl10(brs_prep(d4, ncuts = 100))
left right yt y delta x1
0.48 0.52 0.50 50 3 -0.6146
0.73 0.77 0.75 75 3 -0.7244

Using brs_prep with brs()

Data processed by brs_prep() is automatically detected by brs() - the internal brs_check() step is skipped:

set.seed(42)
n <- 200
dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n))
sim <- brs_sim(
  formula = ~ x1 + x2, data = dat,
  beta = c(0.2, -0.5, 0.3), phi = 0.3,
  link = "logit", link_phi = "logit",
  repar = 2
)
# Remove left, right, yt so brs_prep can rebuild them
prep <- brs_prep(sim[-c(1:3)], ncuts = 100)
fit_prep <- brs(y ~ x1 + x2,
  data = prep, repar = 2,
  link = "logit", link_phi = "logit"
)
summary(fit_prep, digits = 4)
#> 
#> Call:
#> brs(formula = y ~ x1 + x2, data = prep, link = "logit", link_phi = "logit", 
#>     repar = 2)
#> 
#> Quantile residuals:
#>     Min      1Q  Median      3Q     Max 
#> -3.3566 -0.7310 -0.0537  0.6635  2.9618 
#> 
#> Coefficients (mean model with logit link):
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)   0.1946     0.0996   1.954   0.0507 .  
#> x1           -0.5456     0.1079  -5.057 4.26e-07 ***
#> x2            0.2204     0.1067   2.066   0.0389 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Phi coefficients (precision model with logit link):
#>       Estimate Std. Error z value Pr(>|z|)    
#> (phi)   0.3326     0.0910   3.655 0.000257 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---
#> Log-likelihood: -804.9904 on 4 Df | AIC: 1617.9809 | BIC: 1631.1741 
#> Pseudo R-squared: 0.1372  (midpoint approx.; interpret with caution for heavily censored data) 
#> Number of iterations: 28 (BFGS) 
#> Censoring: 152 interval | 19 left | 29 right

Example 1: Fixed dispersion model

Simulating data

We simulate observations from a beta regression model with fixed dispersion, two covariates, and a logit link for the mean.

set.seed(4255)
n <- 250
dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n))

sim_fixed <- brs_sim(
  formula  = ~ x1 + x2,
  data     = dat,
  beta     = c(0.3, -0.6, 0.4),
  phi      = 1 / 10,
  link     = "logit",
  link_phi = "logit",
  ncuts    = 100,
  repar    = 2
)

kbl10(head(sim_fixed, 8))
left right yt y delta x1 x2
0.805 0.815 0.81 81 3 1.9510 -1.3482
0.585 0.595 0.59 59 3 0.7725 1.2966
0.245 0.255 0.25 25 3 0.7264 1.1585
0.965 0.975 0.97 97 3 0.0487 0.6334
0.285 0.295 0.29 29 3 -0.5445 0.1161
0.795 0.805 0.80 80 3 0.3600 -0.2545
0.355 0.365 0.36 36 3 0.7136 0.7694
0.745 0.755 0.75 75 3 -0.7274 1.0397

Each observation is centered in its interval. For example, a score of 67 on a 0–100 scale yields \(y_t=0.67\) with the interval \([0.665,0.675]\).

Fitting the model

fit_fixed <- brs(
  y ~ x1 + x2,
  data     = sim_fixed,
  link     = "logit",
  link_phi = "logit",
  repar    = 2
)
summary(fit_fixed)
#> 
#> Call:
#> brs(formula = y ~ x1 + x2, data = sim_fixed, link = "logit", 
#>     link_phi = "logit", repar = 2)
#> 
#> Quantile residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2.8953 -0.6209  0.0858  0.6729  3.4008 
#> 
#> Coefficients (mean model with logit link):
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.35460    0.08730   4.062 4.87e-05 ***
#> x1          -0.65902    0.09482  -6.950 3.64e-12 ***
#> x2           0.32880    0.09221   3.566 0.000363 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Phi coefficients (precision model with logit link):
#>       Estimate Std. Error z value Pr(>|z|)  
#> (phi)  0.15143    0.08202   1.846   0.0649 .
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---
#> Log-likelihood: -989.6966 on 4 Df | AIC: 1987.3932 | BIC: 2001.4790 
#> Pseudo R-squared: 0.2264  (midpoint approx.; interpret with caution for heavily censored data) 
#> Number of iterations: 29 (BFGS) 
#> Censoring: 198 interval | 14 left | 38 right

The summary output follows the betareg package style, showing separate coefficient tables for the mean and precision submodels, with Wald z-tests and \(p\)-values based on the standard normal distribution.

Estimation and fit diagnostics

optim() (BFGS) maximizes the log-likelihood, evaluated in C++. The gradient and, by default, the Hessian use the chain rule on the two linear predictors, with per-observation derivatives by Richardson central differences (hessian_method = "cpp"; "numDeriv" is also available). Standard errors come from the inverse of the negative Hessian; no generalized inverse is used. brs() also accepts start and control (merged into list(maxit = 5000)).

After the fit, brs() checks the result and stores the numbers:

str(fit_fixed$diagnostics[c("grad_gain", "grad_step", "min_eig_scaled", "hessian_nd", "n_clamped")])
#> List of 5
#>  $ grad_gain     : num 3.06e-10
#>  $ grad_step     : num 2.15e-05
#>  $ min_eig_scaled: num 0.709
#>  $ hessian_nd    : logi TRUE
#>  $ n_clamped     : int 0

Each problem gives a one-line warning:

Warning Meaning What to do
model matrix is rank deficient (error) aliased columns remove the named columns
model matrix is nearly collinear condition number above \(10^4\) drop, combine or center the named columns
Gradient not ~0 at the estimate a Newton step would still gain more than 0.01 in log-likelihood refit with method = "L-BFGS-B"; rescale covariates whose scales differ widely (the message says so)
Hessian near-singular or not negative definite a direction is nearly flat (not identified, or covariates on very different scales) vcov() gives NA where it cannot estimate; rescale, use LR tests or the bootstrap, simplify
observations on the clamp boundary means or shapes at \(10^{-5}\) or \(1-10^{-5}\), e.g. a group entirely at one border merge sparse groups or drop the separating covariate

The classic case of the last two is a data set where every observation is censored on the same side: the likelihood has no finite maximum.

Goodness of fit

kbl10(brs_gof(fit_fixed))
logLik AIC BIC pseudo_r2
-989.6966 1987.393 2001.479 0.2264

Residual diagnostics

residuals() offers several types. Randomized quantile residuals (type = "rqr", Dunn and Smyth, 1996) draw \(u_i\) uniformly on \((F(l_i),F(u_i))\) (or \((0,F(u_i))\), \((F(l_i),1)\) at the borders) and return \(\Phi^{-1}(u_i)\): they are standard normal under the model whatever the censoring, and are the recommended type. Pearson, deviance and weighted residuals are evaluated at the cell center (the midpoint convention of Lopes, 2023) and are unreliable at the borders of the scale. summary() draws the quantile residuals without changing the RNG state.

The plot() method provides six diagnostic panels. By default, the first four are shown:

plot(fit_fixed, caption = NULL, gg = TRUE, title = NULL, theme = ggplot2::theme_bw())

For ggplot2 output (requires the ggplot2 package):

plot(fit_fixed, gg = TRUE)

Predictions

# Fitted means
kbl10(
  data.frame(mu_hat = head(predict(fit_fixed, type = "response"))),
  digits = 4
)
mu_hat
0.2019
0.5675
0.5638
0.6297
0.6795
0.5084 

# Conditional variance
kbl10(
  data.frame(var_hat = head(predict(fit_fixed, type = "variance"))),
  digits = 4
)
var_hat
0.0867
0.1320
0.1323
0.1254
0.1171
0.1344 

# Quantile predictions
kbl10(predict(fit_fixed, type = "quantile", at = c(0.10, 0.25, 0.5, 0.75, 0.90)))
q_0.1 q_0.25 q_0.5 q_0.75 q_0.9
0.0000 0.0006 0.0329 0.3130 0.7400
0.0335 0.2024 0.6326 0.9348 0.9943
0.0322 0.1974 0.6255 0.9323 0.9940
0.0624 0.2999 0.7476 0.9688 0.9982
0.0979 0.3959 0.8302 0.9855 0.9995
0.0168 0.1306 0.5167 0.8863 0.9865
0.0231 0.1596 0.5678 0.9096 0.9905
0.1981 0.5958 0.9385 0.9980 1.0000
0.1001 0.4013 0.8341 0.9862 0.9995
0.0088 0.0860 0.4218 0.8340 0.9755

Confidence intervals

Wald confidence intervals based on the asymptotic normal approximation:

kbl10(confint(fit_fixed))
2.5 % 97.5 %
(Intercept) 0.1835 0.5257
x1 -0.8449 -0.4732
x2 0.1481 0.5095
(phi) -0.0093 0.3122 
kbl10(confint(fit_fixed, model = "mean"))
2.5 % 97.5 %
(Intercept) 0.1835 0.5257
x1 -0.8449 -0.4732
x2 0.1481 0.5095

Censoring structure

The brs_cens() function provides a visual and tabular overview of the censoring types in the fitted model:

brs_cens(fit_fixed, gg = TRUE, inform = TRUE)

Example 2: Variable dispersion model

In many applications, the dispersion parameter \(\phi\) may depend on covariates. The package supports variable-dispersion models using the Formula package notation: y ~ x1 + x2 | z1 + z2, where the terms after | define the linear predictor for \(\phi\). The same brs_sim() function is used for fixed and variable dispersion; the second formula part activates the precision submodel in simulation.

Simulating data

set.seed(2222)
n <- 250
dat_z <- data.frame(
  x1 = rnorm(n),
  x2 = rnorm(n),
  x3 = rbinom(n, size = 1, prob = 0.5),
  z1 = rnorm(n),
  z2 = rnorm(n)
)

sim_var <- brs_sim(
  formula = ~ x1 + x2 + x3 | z1 + z2,
  data = dat_z,
  beta = c(0.2, -0.6, 0.2, 0.2),
  zeta = c(0.2, -0.8, 0.6),
  link = "logit",
  link_phi = "logit",
  ncuts = 100,
  repar = 2
)

kbl10(head(sim_var, 10))
left right yt y delta x1 x2 x3 z1 z2
0.000 0.005 0.00 0 1 -0.3381 0.7980 1 -0.1859 1.7165
0.000 0.005 0.00 0 1 0.9392 -2.2067 0 -1.1473 0.8539
0.985 0.995 0.99 99 3 1.7377 0.0134 0 0.3364 0.9417
0.535 0.545 0.54 54 3 0.6963 0.2765 1 2.1178 0.1317
0.515 0.525 0.52 52 3 0.4623 -0.8707 1 0.5450 -0.3795
0.005 0.015 0.01 1 3 -0.3151 -1.0969 1 -0.4672 -0.8259
0.005 0.015 0.01 1 3 0.1927 1.8391 0 -1.2001 -0.1784
0.035 0.045 0.04 4 3 1.1307 -1.3548 0 -0.6740 -0.9092
0.315 0.325 0.32 32 3 1.9764 0.6912 1 0.9831 0.0286
0.825 0.835 0.83 83 3 1.2071 0.2571 0 0.6295 -0.2415

Fitting the model

fit_var <- brs(
  y ~ x1 + x2 | z1,
  data     = sim_var,
  link     = "logit",
  link_phi = "logit",
  repar    = 2
)
summary(fit_var)
#> 
#> Call:
#> brs(formula = y ~ x1 + x2 | z1, data = sim_var, link = "logit", 
#>     link_phi = "logit", repar = 2)
#> 
#> Quantile residuals:
#>     Min      1Q  Median      3Q     Max 
#> -3.2736 -0.5585  0.0259  0.6191  2.5197 
#> 
#> Coefficients (mean model with logit link):
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.31885    0.08541   3.733 0.000189 ***
#> x1          -0.48245    0.08961  -5.384 7.28e-08 ***
#> x2           0.23904    0.08270   2.890 0.003847 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Phi coefficients (precision model with logit link):
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.33338    0.08621   3.867  0.00011 ***
#> z1          -0.81258    0.09204  -8.829  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---
#> Log-likelihood: -936.0499 on 5 Df | AIC: 1882.0998 | BIC: 1899.7071 
#> Pseudo R-squared: 0.1164  (midpoint approx.; interpret with caution for heavily censored data) 
#> Number of iterations: 42 (BFGS) 
#> Censoring: 181 interval | 23 left | 46 right

Notice the (phi)_ prefix in the precision coefficient names, following the betareg convention.

Accessing coefficients by submodel

# Full parameter vector
round(coef(fit_var), 4)
#>       (Intercept)                x1                x2 (phi)_(Intercept) 
#>            0.3188           -0.4825            0.2390            0.3334 
#>          (phi)_z1 
#>           -0.8126

# Mean submodel only
round(coef(fit_var, model = "mean"), 4)
#> (Intercept)          x1          x2 
#>      0.3188     -0.4825      0.2390

# Precision submodel only
round(coef(fit_var, model = "precision"), 4)
#> (phi)_(Intercept)          (phi)_z1 
#>            0.3334           -0.8126

# Variance-covariance matrix for the mean submodel
kbl10(vcov(fit_var, model = "mean"))
(Intercept) x1 x2
(Intercept) 0.0073 -3e-04 0.0002
x1 -0.0003 8e-03 0.0001
x2 0.0002 1e-04 0.0068

Diagnostics for variable dispersion

plot(fit_var)

Advanced analyst functions

The package includes analyst-facing helpers for uncertainty quantification, effect interpretation, score-scale communication, and predictive validation.

Parametric bootstrap confidence intervals

Each replicate draws a new response from the fitted model and refits the original formula on a copy of the data. Each row keeps its observation mechanism: exact values stay continuous, scores are re-coarsened on the fit’s grid, and analyst thresholds from brs_prep() are kept as fixed, non-informative thresholds (the censoring type is re-drawn by the side of them where the new value falls). Refits start from the fitted estimate. Failed replicates are counted and printed. Use R of at least 199 in practice; ci_type = "bca" costs one extra refit per observation.

set.seed(101)
boot_ci <- brs_bootstrap(
  fit_fixed,
  R = 30,
  level = 0.95,
  ci_type = "percentile",
  keep_draws = TRUE
)
kbl10(head(boot_ci, 10))
parameter estimate se_boot ci_lower ci_upper mcse_lower mcse_upper wald_lower wald_upper level
(Intercept) 0.3546 0.0982 0.1768 0.5276 0.0291 0.0298 0.1835 0.5257 0.95
x1 -0.6590 0.1044 -0.8213 -0.4645 0.0213 0.0344 -0.8449 -0.4732 0.95
x2 0.3288 0.1088 0.1437 0.5298 0.0244 0.0287 0.1481 0.5095 0.95
(phi) 0.1514 0.0913 0.0432 0.3201 0.0101 0.0193 -0.0093 0.3122 0.95 
autoplot.brs_bootstrap(
  boot_ci,
  type = "ci_forest",
  title = "Bootstrap (percentile) vs Wald intervals"
)

Average marginal effects (AME)

set.seed(202)
ame <- brs_marginaleffects(
  fit_fixed,
  model = "mean",
  type = "response",
  interval = TRUE,
  n_sim = 60,
  keep_draws = TRUE
)
kbl10(ame)
variable ame std.error ci.lower ci.upper model type n
x1 -0.1425 0.0175 -0.1764 -0.1123 mean response 250
x2 0.0711 0.0184 0.0378 0.1055 mean response 250 
autoplot.brs_marginaleffects(ame, type = "forest")

Score probabilities on the original scale

prob_scores <- brs_predict_scoreprob(fit_fixed, scores = 0:10)
kbl10(prob_scores[1:6, 1:6])
score_0 score_1 score_2 score_3 score_4 score_5
0.3600 0.0758 0.0406 0.0289 0.0228 0.0190
0.0393 0.0280 0.0192 0.0157 0.0136 0.0122
0.0403 0.0285 0.0195 0.0159 0.0137 0.0123
0.0252 0.0205 0.0147 0.0122 0.0108 0.0098
0.0171 0.0155 0.0115 0.0097 0.0087 0.0080
0.0587 0.0363 0.0240 0.0191 0.0163 0.0145 
kbl10(
  data.frame(row_sum = rowSums(prob_scores)[1:6]),
  digits = 4
)
row_sum
0.6136
0.1774
0.1802
0.1340
0.1044
0.2261

Repeated k-fold cross-validation

set.seed(303) # For cross-validation reproducibility
cv_res <- brs_cv(
  y ~ x1 + x2,
  data = sim_fixed,
  k = 5,
  repeats = 2,
  repar = 2,
)
kbl10(cv_res)
repeat fold n_train n_test log_score rmse_yt mae_yt converged error
1 1 200 50 -4.0896 0.3338 0.2852 TRUE NA
1 2 200 50 -4.0065 0.3482 0.3080 TRUE NA
1 3 200 50 -3.7595 0.3069 0.2693 TRUE NA
1 4 200 50 -3.9121 0.3660 0.3281 TRUE NA
1 5 200 50 -4.0921 0.3834 0.3422 TRUE NA
2 1 200 50 -4.2336 0.3669 0.3191 TRUE NA
2 2 200 50 -3.8826 0.3231 0.2829 TRUE NA
2 3 200 50 -3.8884 0.3151 0.2787 TRUE NA
2 4 200 50 -4.0835 0.3745 0.3272 TRUE NA
2 5 200 50 -3.8069 0.3697 0.3366 TRUE NA 
kbl10(
  data.frame(
    metric = c("log_score", "rmse_yt", "mae_yt"),
    mean = c(
      mean(cv_res$log_score, na.rm = TRUE),
      mean(cv_res$rmse_yt, na.rm = TRUE),
      mean(cv_res$mae_yt, na.rm = TRUE)
    )
  ),
  digits = 4
)
metric mean
log_score -3.9755
rmse_yt 0.3488
mae_yt 0.3077

S3 methods reference

The following standard S3 methods are available for objects of class "brs":

Method Description
print() Compact display of call and coefficients
summary() Detailed output with Wald tests and goodness-of-fit
coef(model=) Extract coefficients (full, mean, or precision)
vcov(model=) Variance-covariance matrix (full, mean, or precision)
confint(model=) Wald confidence intervals
logLik() Log-likelihood value
AIC(), BIC() Information criteria
nobs() Number of observations
formula() Model formula
model.matrix(model=) Design matrix (mean or precision)
fitted() Fitted mean values
residuals(type=) Residuals: response, pearson, deviance, rqr, weighted, sweighted
predict(type=) Predictions: response, link, precision, variance, quantile, score, expected_score
plot(gg=) Diagnostic plots (base R or ggplot2)

Reparameterizations

The package supports three reparameterizations of the beta distribution, controlled by the repar argument. Each models different parameters, so the admissible links differ; the first link of each cell is the default.

repar parameters shapes \((a,b)\) link link_phi
0 shapes \(p,q>0\) \((p,q)\) log, sqrt log, sqrt
1 mean \(\mu\), precision \(\phi>0\) \((\mu\phi,(1-\mu)\phi)\) logit, probit, cauchit, cloglog log, sqrt
2 mean \(\mu\), dispersion \(\phi\in(0,1)\) \((\mu(1-\phi)/\phi,(1-\mu)(1-\phi)/\phi)\) logit, probit, cauchit, cloglog logit, probit, cauchit, cloglog

In every scheme \(E[Y]=a/(a+b)\) and \(\mathrm{Var}[Y]=E[Y](1-E[Y])/(1+a+b)\). Under repar = 1 this is \(\mu(1-\mu)/(1+\phi)\) (Ferrari and Cribari-Neto, 2004); under repar = 2 it is \(\phi\,\mu(1-\mu)\), so the dispersion is \(\phi=1/(1+a+b)\), the variance relative to its maximum \(\mu(1-\mu)\) (not a coefficient of variation). Under repar = 0 the first parameter is the shape \(p\) and the mean is \(p/(p+q)\).

# Precision parameterization: mu = 0.5, phi = 10 (high precision)
brs_repar(mu = 0.5, phi = 10, repar = 1)
#>   shape1 shape2
#> 1      5      5

# Mean-variance parameterization: mu = 0.5, phi = 0.1 (low dispersion)
brs_repar(mu = 0.5, phi = 0.1, repar = 2)
#>   shape1 shape2
#> 1    4.5    4.5

References

  • Lopes, J. E. (2023). Modelos de regressão beta para dados de escala. Master’s dissertation, Universidade Federal do Paraná, Curitiba. https://hdl.handle.net/1884/86624.

  • Dunn, P. K., and Smyth, G. K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5(3), 236–244.

  • 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. Validated online via: https://doi.org/10.1080/0266476042000214501.

  • 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. Validated online via: https://doi.org/10.1037/1082-989X.11.1.54.

  • 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. Validated online via: https://doi.org/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. Validated online via: https://doi.org/10.1016/j.jpainsymman.2010.08.016.