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 normalisation 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, so anova() refuses to compare 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	-1.4000
0.25	0.35	0.3	3	3	0.2553
0.45	0.55	0.5	5	3	-2.4373
0.65	0.75	0.7	7	3	-0.0056
0.95	1.00	1.0	10	2	0.6216

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	1.1484
0.000	0.005	0.00	0	1	-1.8218
0.985	1.000	0.99	99	2	-0.2473
0.495	0.505	0.50	50	3	-0.2442

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.2827
0.2	1.00	0.600	60.0	2	-0.5537
0.3	0.45	0.375	37.5	3	0.6290
0.5	0.50	0.500	50.0	0	2.0650

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	-1.6310
0.73	0.77	0.75	75	3	0.5124

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) maximises 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 generalised 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 centre 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

Comparing link functions

The package supports several link functions for the mean submodel. We can compare them using information criteria:

links <- c("logit", "probit", "cauchit", "cloglog")
fits <- lapply(setNames(links, links), function(lnk) {
  brs(y ~ x1 + x2, data = sim_fixed, link = lnk, repar = 2)
})

# Estimates
est_table <- do.call(rbind, lapply(names(fits), function(lnk) {
  e <- brs_est(fits[[lnk]])
  e$link <- lnk
  e
}))
kbl10(est_table)

variable	estimate	se	z_value	p_value	ci_lower	ci_upper	link
(Intercept)	0.3546	0.0873	4.0619	0.0000	0.1835	0.5257	logit
x1	-0.6590	0.0948	-6.9505	0.0000	-0.8449	-0.4732	logit
x2	0.3288	0.0922	3.5657	0.0004	0.1481	0.5095	logit
(phi)	0.1514	0.0820	1.8463	0.0649	-0.0093	0.3122	logit
(Intercept)	0.2180	0.0532	4.0969	0.0000	0.1137	0.3223	probit
x1	-0.3964	0.0548	-7.2312	0.0000	-0.5038	-0.2890	probit
x2	0.1955	0.0547	3.5745	0.0004	0.0883	0.3026	probit
(phi)	0.1539	0.0819	1.8794	0.0602	-0.0066	0.3144	probit
(Intercept)	0.3062	0.0824	3.7164	0.0002	0.1447	0.4677	cauchit
x1	-0.6286	0.1100	-5.7153	0.0000	-0.8441	-0.4130	cauchit

# Goodness of fit
gof_table <- do.call(rbind, lapply(fits, brs_gof))
kbl10(gof_table)

	logLik	AIC	BIC	pseudo_r2
logit	-989.6966	1987.393	2001.479	0.2264
probit	-989.8231	1987.646	2001.732	0.2229
cauchit	-989.7611	1987.522	2001.608	0.1804
cloglog	-989.7198	1987.439	2001.525	0.1437

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 centre (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: 43 (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

Comparing link functions (variable dispersion)

links <- c("logit", "probit", "cauchit", "cloglog")
fits_var <- lapply(setNames(links, links), function(lnk) {
  brs(y ~ x1 + x2 | z1, data = sim_var, link = lnk, repar = 2)
})

# Estimates
est_var <- do.call(rbind, lapply(names(fits_var), function(lnk) {
  e <- brs_est(fits_var[[lnk]])
  e$link <- lnk
  e
}))
kbl10(est_var)

variable	estimate	se	z_value	p_value	ci_lower	ci_upper	link
(Intercept)	0.3188	0.0854	3.7333	0.0002	0.1515	0.4862	logit
x1	-0.4825	0.0896	-5.3841	0.0000	-0.6581	-0.3068	logit
x2	0.2390	0.0827	2.8904	0.0038	0.0770	0.4011	logit
(phi)_(Intercept)	0.3334	0.0862	3.8669	0.0001	0.1644	0.5024	logit
(phi)_z1	-0.8126	0.0920	-8.8287	0.0000	-0.9930	-0.6322	logit
(Intercept)	0.1972	0.0524	3.7601	0.0002	0.0944	0.3000	probit
x1	-0.2980	0.0541	-5.5038	0.0000	-0.4041	-0.1919	probit
x2	0.1477	0.0505	2.9237	0.0035	0.0487	0.2467	probit
(phi)_(Intercept)	0.3326	0.0862	3.8571	0.0001	0.1636	0.5015	probit
(phi)_z1	-0.8131	0.0920	-8.8379	0.0000	-0.9934	-0.6327	probit

# Goodness of fit
gof_var <- do.call(rbind, lapply(fits_var, brs_gof))
kbl10(gof_var)

	logLik	AIC	BIC	pseudo_r2
logit	-936.0499	1882.100	1899.707	0.1164
probit	-935.9539	1881.908	1899.515	0.1211
cauchit	-936.7522	1883.504	1901.112	0.0957
cloglog	-935.6783	1881.357	1898.964	0.0667

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.