Analyst Tools for betaregscale

Complete likelihood by censoring type

For each observation \(i\), let \(\delta_i\in\{0,1,2,3\}\) indicate exact, left-censored, right-censored, or interval-censored status. With beta CDF \(F(\cdot)\), beta density \(f(\cdot)\), and interval endpoints \([l_i,u_i]\), the contribution is:

\[ L_i(\theta)= \begin{cases} f(y_i; a_i, b_i), & \delta_i=0,\\ F(u_i; a_i, b_i), & \delta_i=1,\\ 1 - F(l_i; a_i, b_i), & \delta_i=2,\\ F(u_i; a_i, b_i)-F(l_i; a_i, b_i), & \delta_i=3. \end{cases} \]

This is the basis for fitting, prediction, and validation metrics.

Model-comparison metrics

brs_table() reports:

• \(\log L(\hat\theta)\),
• \(AIC=-2\log L(\hat\theta)+2k\),
• \(BIC=-2\log L(\hat\theta)+k\log n\),

where \(k\) is the number of estimated parameters and \(n\) is the sample size.

Average marginal effects (AME)

brs_marginaleffects() computes AME by central finite differences:

\[ \mathrm{AME}_j=\frac{1}{n}\sum_{i=1}^n\frac{\hat g_i(x_{ij}+h_j)-\hat g_i(x_{ij}-h_j)}{2h_j}, \qquad h_j=h\max\{\mathrm{sd}(x_j),1\}, \]

with \(\hat g_i\) on the requested prediction scale (response, the mean \(E[Y]\), or link) and \(h=10^{-5}\) by default. Only numeric covariates are used; code a factor as 0/1 dummies to get its contrast. For binary covariates \(x_j\in\{0,1\}\), it uses the discrete contrast \(\hat g(x_j=1)-\hat g(x_j=0)\).

Score-scale probabilities

For integer scores \(s\in\{0,\dots,K\}\) (\(K\) = ncuts, the maximum score), brs_predict_scoreprob() computes, under interval = "mid":

\[ P(Y=s)= \begin{cases} F(\mathrm{lim}/K), & s=0,\\ 1-F((K-\mathrm{lim})/K), & s=K,\\ F((s+\mathrm{lim})/K)-F((s-\mathrm{lim})/K), & 1\le s\le K-1. \end{cases} \]

Under "right" and "left" the cells are \([s,s+1]/(K+1)\), so \(P(Y=s)=F((s+1)/(K+1))-F(s/(K+1))\) with \(F(0)=0\) and \(F(1)=1\). The probabilities sum to 1 whenever the cells partition \([0,1]\) (always for "right"/"left", and for "mid" with lim = 0.5). They are directly aligned with interval geometry on the original instrument scale.

Cross-validation log score

In brs_cv(), fold-level predictive quality includes:

\[ \mathrm{log\_score}=\frac{1}{n_{test}}\sum_{i \in test}\log(p_i), \]

where \(p_i\) is the predictive contribution from the same censoring-rule piecewise definition shown above.

It also reports:

\[ \mathrm{RMSE}_{yt}=\sqrt{\frac{1}{n_{test}}\sum(y_{t,i}-\hat\mu_i)^2}, \qquad \mathrm{MAE}_{yt}=\frac{1}{n_{test}}\sum|y_{t,i}-\hat\mu_i|. \]

1) Simulate data and fit candidate models

set.seed(2026)
n <- 220
d <- data.frame(
  x1 = rnorm(n),
  x2 = rnorm(n),
  z1 = rnorm(n)
)

sim <- brs_sim(
  formula = ~ x1 + x2 | z1,
  data = d,
  beta = c(0.20, -0.45, 0.25),
  zeta = c(0.15, -0.30),
  ncuts = 100,
  repar = 2
)

fit_fixed <- brs(y ~ x1 + x2, data = sim, repar = 2)
fit_var <- brs(y ~ x1 + x2 | z1, data = sim, repar = 2)

2) Compare models in one table

tab <- brs_table(
  fixed = fit_fixed,
  variable = fit_var,
  sort_by = "AIC"
)
kbl10(tab)

3) Estimate average marginal effects

set.seed(2026) # For marginal effects simulation
me_mean <- brs_marginaleffects(
  fit_var,
  model = "mean",
  type = "response",
  interval = TRUE,
  n_sim = 120
)
kbl10(me_mean)

set.seed(2026) # Reset seed for reproducibility
me_precision <- brs_marginaleffects(
  fit_var,
  model = "precision",
  type = "link",
  interval = TRUE,
  n_sim = 120
)
kbl10(me_precision)

4) Predict score probabilities

P <- brs_predict_scoreprob(fit_var, scores = 0:10)
dim(P)
#> [1] 220  11
kbl10(P[1:6, 1:5])

kbl10(
  data.frame(row_sum = rowSums(P)),
  digits = 4
)

rowSums(P) should be close to 1 (up to numerical tolerance and score subset selection).

5) ggplot2 diagnostics

autoplot.brs(fit_var, type = "calibration")

autoplot.brs(fit_var, type = "score_dist", scores = 0:20)

autoplot.brs(fit_var, type = "cdf", max_curves = 4)

autoplot.brs(fit_var, type = "residuals_by_delta", residual_type = "rqr")

6) Repeated k-fold cross-validation

set.seed(2026) # For cross-validation reproducibility
cv_res <- brs_cv(
  y ~ x1 + x2 | z1,
  data = sim,
  k = 3,
  repeats = 1,
  repar = 2
)

kbl10(cv_res)

kbl10(
  data.frame(
    metric = c("log_score", "rmse_yt", "mae_yt"),
    mean = colMeans(cv_res[, c("log_score", "rmse_yt", "mae_yt")], na.rm = TRUE)
  ),
  digits = 4
)