Mathematical model
Assume observations \(i = 1, \dots,
n_j\) within groups \(j = 1, \dots,
G\), with group-specific random-effects vector \(\mathbf{b}_j \in \mathbb{R}^{q_b}\).
Linear predictors
\[
\eta_{\mu,ij} = x_{ij}^\top \beta + w_{ij}^\top \mathbf{b}_j,
\qquad
\eta_{\phi,ij} = z_{ij}^\top \gamma
\]
\[
\mu_{ij} = g^{-1}(\eta_{\mu,ij}),
\qquad
\phi_{ij} = h^{-1}(\eta_{\phi,ij})
\]
with \(g(\cdot)\) and \(h(\cdot)\) chosen by link and
link_phi. The random-effects design row \(w_{ij}\) is defined by
random = ~ terms | group.
Beta parameterization
For each \((\mu_{ij},\phi_{ij})\),
repar maps to beta shape parameters \((a_{ij},b_{ij})\) via
brs_repar().
Conditional contribution by censoring type
Each observation contributes:
\[
L_{ij}(b_j;\theta)=
\begin{cases}
f(y_{ij}; a_{ij}, b_{ij}), & \delta_{ij}=0\\
F(u_{ij}; a_{ij}, b_{ij}), & \delta_{ij}=1\\
1 - F(l_{ij}; a_{ij}, b_{ij}), & \delta_{ij}=2\\
F(u_{ij}; a_{ij}, b_{ij}) - F(l_{ij}; a_{ij}, b_{ij}), &
\delta_{ij}=3
\end{cases}
\]
where \(l_{ij},u_{ij}\) are interval
endpoints on \((0,1)\), \(f(\cdot)\) is beta density, and \(F(\cdot)\) is beta CDF.
Random-effects distribution
\[
\mathbf{b}_j \sim \mathcal{N}(\mathbf{0}, D),
\]
where \(D\) is a symmetric
positive-definite covariance matrix. Internally, brsmm()
optimizes a packed lower-Cholesky parameterization \(D = LL^\top\) (diagonal entries on
log-scale for positivity).
Group marginal likelihood
\[
L_j(\theta)=\int_{\mathbb{R}^{q_b}}
\prod_{i=1}^{n_j} L_{ij}(b_j;\theta)\,
\varphi_{q_b}(\mathbf{b}_j;\mathbf{0},D)\,d\mathbf{b}_j
\]
\[
\ell(\theta)=\sum_{j=1}^G \log L_j(\theta)
\]
Laplace approximation used by brsmm()
Define \[
Q_j(\mathbf{b})=
\sum_{i=1}^{n_j}\log L_{ij}(\mathbf{b};\theta)+
\log\varphi_{q_b}(\mathbf{b};\mathbf{0},D)
\] and \(\hat{\mathbf{b}}_j=\arg\max_{\mathbf{b}}
Q_j(\mathbf{b})\), with curvature \[
H_j = -\nabla^2 Q_j(\hat{\mathbf{b}}_j).
\] Then
\[
\log L_j(\theta) \approx
Q_j(\hat{\mathbf{b}}_j) +
\frac{q_b}{2}\log(2\pi) -
\frac{1}{2}\log|H_j|.
\]
brsmm() maximizes the approximated \(\ell(\theta)\) with
stats::optim(), and computes group-level posterior modes
\(\hat{\mathbf{b}}_j\). For \(q_b = 1\), this reduces to the scalar
random-intercept formula.
Simulating clustered scale data
The next helper simulates data from a known mixed model to illustrate
fitting, inference, and recovery checks.
sim_brsmm_data <- function(seed = 3501L, g = 24L, ni = 12L,
beta = c(0.20, 0.65),
gamma = c(-0.15),
sigma_b = 0.55) {
set.seed(seed)
id <- factor(rep(seq_len(g), each = ni))
n <- length(id)
x1 <- rnorm(n)
b_true <- rnorm(g, mean = 0, sd = sigma_b)
eta_mu <- beta[1] + beta[2] * x1 + b_true[as.integer(id)]
eta_phi <- rep(gamma[1], n)
mu <- plogis(eta_mu)
phi <- plogis(eta_phi)
shp <- brs_repar(mu = mu, phi = phi, repar = 2)
y <- round(stats::rbeta(n, shp$shape1, shp$shape2) * 100)
list(
data = data.frame(y = y, x1 = x1, id = id),
truth = list(beta = beta, gamma = gamma, sigma_b = sigma_b, b = b_true)
)
}
sim <- sim_brsmm_data(
g = 12,
ni = 20,
beta = c(0.20, 0.65),
gamma = c(-0.15),
sigma_b = 0.55
)
str(sim$data)
#> 'data.frame': 240 obs. of 3 variables:
#> $ y : num 0 83 99 56 65 10 8 1 65 98 ...
#> $ x1: num -0.3677 -2.0069 -0.0469 -0.2468 0.7634 ...
#> $ id: Factor w/ 12 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
Fitting brsmm()
fit_mm <- brsmm(
y ~ x1,
random = ~ 1 | id,
data = sim$data,
repar = 2,
int_method = "laplace",
method = "BFGS",
control = list(maxit = 1000)
)
summary(fit_mm)
#>
#> Call:
#> brsmm(formula = y ~ x1, random = ~1 | id, data = sim$data, repar = 2,
#> int_method = "laplace", method = "BFGS", control = list(maxit = 1000))
#>
#> Randomized Quantile Residuals:
#> Min 1Q Median 3Q Max
#> -2.7630 -0.6319 0.0210 0.7206 3.5456
#>
#> Coefficients (mean model with logit link):
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.36832 0.15437 2.386 0.017 *
#> x1 0.63301 0.09465 6.688 2.27e-11 ***
#> ---
#> 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.15937 0.08464 -1.883 0.0597 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Random-effects parameters (Cholesky scale):
#> Estimate Std. Error z value Pr(>|z|)
#> logSD.(Intercept)|id -0.7973 0.2925 -2.726 0.00641 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---
#> Mixed beta interval model (Laplace)
#> Observations: 240 | Groups: 12
#> Log-likelihood: -1008.2467 on 4 Df | AIC: 2024.4934 | BIC: 2038.4160
#> Pseudo R-squared: 0.1364
#> Number of iterations: 35 (BFGS)
#> Censoring: 212 interval | 8 left | 20 right
Random intercept + slope example
The model below includes a random intercept and random slope for
x1:
fit_mm_rs <- brsmm(
y ~ x1,
random = ~ 1 + x1 | id,
data = sim$data,
repar = 2,
int_method = "laplace",
method = "BFGS",
control = list(maxit = 1200)
)
summary(fit_mm_rs)
#>
#> Call:
#> brsmm(formula = y ~ x1, random = ~1 + x1 | id, data = sim$data,
#> repar = 2, int_method = "laplace", method = "BFGS", control = list(maxit = 1200))
#>
#> Randomized Quantile Residuals:
#> Min 1Q Median 3Q Max
#> -2.9374 -0.6419 0.0104 0.6845 2.3799
#>
#> Coefficients (mean model with logit link):
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.3533 0.1575 2.244 0.0249 *
#> x1 0.6292 0.1057 5.955 2.61e-09 ***
#> ---
#> 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.17469 0.08485 -2.059 0.0395 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Random-effects parameters (Cholesky scale):
#> Estimate Std. Error z value Pr(>|z|)
#> logSD.(Intercept)|id -0.7689 0.2859 -2.689 0.00716 **
#> cov.x1:(Intercept)|id -0.1576 0.1072 -1.470 0.14161
#> logSD.x1|id -4.8354 10.4522 -0.463 0.64364
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---
#> Mixed beta interval model (Laplace)
#> Observations: 240 | Groups: 12
#> Log-likelihood: -1007.0654 on 6 Df | AIC: 2026.1308 | BIC: 2047.0147
#> Pseudo R-squared: 0.1364
#> Number of iterations: 32 (BFGS)
#> Censoring: 212 interval | 8 left | 20 right
Covariance structure of random effects:
kbl10(fit_mm_rs$random$D)
| 0.2149 |
-0.0731 |
| -0.0731 |
0.0249 |
kbl10(
data.frame(term = names(fit_mm_rs$random$sd_b), sd = as.numeric(fit_mm_rs$random$sd_b)),
digits = 4
)
| (Intercept) |
0.4635 |
| x1 |
0.1578 |
kbl10(head(ranef(fit_mm_rs), 10))
| -0.0040 |
0.0013 |
| 0.2585 |
-0.0880 |
| 0.0430 |
-0.0146 |
| 0.1528 |
-0.0518 |
| -0.0354 |
0.0119 |
| 0.6182 |
-0.2102 |
| -0.4101 |
0.1396 |
| -0.1081 |
0.0366 |
| -0.4720 |
0.1606 |
| -0.6703 |
0.2280 |
Additional studies of random effects (numerical and visual)
Following practices from established mixed-models packages, the
package now allows for a dedicated study of the random effects focusing
on:
- \(D\) structure and
correlation;
- empirical distribution of modes by group;
- empirical shrinkage intensity;
- specific visual diagnostics for the random components.
re_study <- brsmm_re_study(fit_mm_rs)
print(re_study)
#>
#> Random-effects study
#> Groups: 12
#>
#> Random-effects (VarCorr):
#> Name Std.Dev. Corr
#> re1 0.4635
#> re2 0.1578 -0.9987
#>
#> ICC (latent logistic scale): 0.0613
#>
#> Summary by term (SD_model = model SD; shrinkage = Var(modes)/Var(model)):
#> term sd_model mean_mode sd_mode shrinkage_ratio shapiro_p
#> (Intercept) 0.4635 6e-04 0.4160 0.8053 0.8381
#> x1 0.1578 -2e-04 0.1415 0.8036 0.8369
kbl10(re_study$summary)
| (Intercept) |
0.4635 |
6e-04 |
0.4160 |
0.8053 |
0.8381 |
| x1 |
0.1578 |
-2e-04 |
0.1415 |
0.8036 |
0.8369 |
| 0.2149 |
-0.0731 |
| -0.0731 |
0.0249 |
| 1.0000 |
-0.9987 |
| -0.9987 |
1.0000 |
Suggested visualizations for random effects:
if (requireNamespace("ggplot2", quietly = TRUE)) {
autoplot.brsmm(fit_mm_rs, type = "ranef_caterpillar")
autoplot.brsmm(fit_mm_rs, type = "ranef_density")
autoplot.brsmm(fit_mm_rs, type = "ranef_pairs")
autoplot.brsmm(fit_mm_rs, type = "ranef_qq")
}
Core methods
Coefficients and random effects
coef(fit_mm, model = "random") returns packed
random-effect covariance parameters on the optimizer scale
(lower-Cholesky, with a log-diagonal). For random-intercept models, this
simplifies to \(\log \sigma_b\).
kbl10(
data.frame(
parameter = names(coef(fit_mm, model = "full")),
estimate = as.numeric(coef(fit_mm, model = "full"))
),
digits = 4
)
| (Intercept) |
0.3683 |
| x1 |
0.6330 |
| (phi)_(Intercept) |
-0.1594 |
| (re_chol_logsd)_(Intercept)|id |
-0.7973 |
kbl10(
data.frame(
log_sigma_b = as.numeric(coef(fit_mm, model = "random")),
sigma_b = as.numeric(exp(coef(fit_mm, model = "random")))
),
digits = 4
)
kbl10(head(ranef(fit_mm), 10))
| -0.0420 |
| 0.2152 |
| 0.0468 |
| 0.2022 |
| -0.0720 |
| 0.6098 |
| -0.3624 |
| -0.1838 |
| -0.4129 |
| -0.6105 |
For random intercept + slope models:
kbl10(
data.frame(
parameter = names(coef(fit_mm_rs, model = "random")),
estimate = as.numeric(coef(fit_mm_rs, model = "random"))
),
digits = 4
)
| (re_chol_logsd)_(Intercept)|id |
-0.7689 |
| (re_chol)_x1:(Intercept)|id |
-0.1576 |
| (re_chol_logsd)_x1|id |
-4.8354 |
kbl10(fit_mm_rs$random$D)
| 0.2149 |
-0.0731 |
| -0.0731 |
0.0249 |
Variance-covariance, summary and likelihood criteria
vc <- vcov(fit_mm)
dim(vc)
#> [1] 4 4
sm <- summary(fit_mm)
kbl10(sm$coefficients)
| (Intercept) |
0.3683 |
0.1544 |
2.3860 |
0.017 |
-0.1594 |
0.0846 |
-1.8829 |
0.0597 |
-0.7973 |
0.2925 |
-2.7259 |
0.0064 |
| x1 |
0.6330 |
0.0947 |
6.6877 |
0.000 |
-0.1594 |
0.0846 |
-1.8829 |
0.0597 |
-0.7973 |
0.2925 |
-2.7259 |
0.0064 |
kbl10(
data.frame(
logLik = as.numeric(logLik(fit_mm)),
AIC = AIC(fit_mm),
BIC = BIC(fit_mm),
nobs = nobs(fit_mm)
),
digits = 4
)
| -1008.247 |
2024.493 |
2038.416 |
240 |
Fitted values, prediction and residuals
kbl10(
data.frame(
mu_hat = head(fitted(fit_mm, type = "mu")),
phi_hat = head(fitted(fit_mm, type = "phi")),
pred_mu = head(predict(fit_mm, type = "response")),
pred_eta = head(predict(fit_mm, type = "link")),
pred_phi = head(predict(fit_mm, type = "precision")),
pred_var = head(predict(fit_mm, type = "variance"))
),
digits = 4
)
| 0.5234 |
0.4602 |
0.5234 |
0.0936 |
0.4602 |
0.1148 |
| 0.2801 |
0.4602 |
0.2801 |
-0.9441 |
0.4602 |
0.0928 |
| 0.5736 |
0.4602 |
0.5736 |
0.2967 |
0.4602 |
0.1126 |
| 0.5424 |
0.4602 |
0.5424 |
0.1701 |
0.4602 |
0.1142 |
| 0.6920 |
0.4602 |
0.6920 |
0.8096 |
0.4602 |
0.0981 |
| 0.4792 |
0.4602 |
0.4792 |
-0.0834 |
0.4602 |
0.1149 |
kbl10(
data.frame(
res_response = head(residuals(fit_mm, type = "response")),
res_pearson = head(residuals(fit_mm, type = "pearson"))
),
digits = 4
)
| -0.5234 |
-1.5446 |
| 0.5499 |
1.8052 |
| 0.4164 |
1.2410 |
| 0.0176 |
0.0520 |
| -0.0420 |
-0.1342 |
| -0.3792 |
-1.1188 |
Diagnostic plotting methods
plot.brsmm() supports base and ggplot2 backends:
plot(fit_mm, which = 1:4, type = "pearson")
if (requireNamespace("ggplot2", quietly = TRUE)) {
plot(fit_mm, which = 1:2, gg = TRUE)
}
autoplot.brsmm() provides focused ggplot
diagnostics:
if (requireNamespace("ggplot2", quietly = TRUE)) {
autoplot.brsmm(fit_mm, type = "calibration")
autoplot.brsmm(fit_mm, type = "score_dist")
autoplot.brsmm(fit_mm, type = "ranef_qq")
autoplot.brsmm(fit_mm, type = "residuals_by_group")
}
Prediction with newdata
If newdata contains unseen groups,
predict.brsmm() uses a random effect equal to zero for
those levels.
nd <- sim$data[1:8, c("x1", "id")]
kbl10(
data.frame(pred_seen = as.numeric(predict(fit_mm, newdata = nd, type = "response"))),
digits = 4
)
| 0.5234 |
| 0.2801 |
| 0.5736 |
| 0.5424 |
| 0.6920 |
| 0.4792 |
| 0.6047 |
| 0.4668 |
nd_unseen <- nd
nd_unseen$id <- factor(rep("new_cluster", nrow(nd_unseen)))
kbl10(
data.frame(pred_unseen = as.numeric(predict(fit_mm, newdata = nd_unseen, type = "response"))),
digits = 4
)
| 0.5338 |
| 0.2886 |
| 0.5839 |
| 0.5528 |
| 0.7009 |
| 0.4897 |
| 0.6147 |
| 0.4773 |
The same logic applies to random intercept + slope models:
kbl10(
data.frame(pred_rs_seen = as.numeric(predict(fit_mm_rs, newdata = nd, type = "response"))),
digits = 4
)
| 0.5294 |
| 0.2858 |
| 0.5793 |
| 0.5483 |
| 0.6965 |
| 0.4853 |
| 0.6101 |
| 0.4730 |
kbl10(
data.frame(pred_rs_unseen = as.numeric(predict(fit_mm_rs, newdata = nd_unseen, type = "response"))),
digits = 4
)
| 0.5305 |
| 0.2871 |
| 0.5803 |
| 0.5493 |
| 0.6971 |
| 0.4865 |
| 0.6110 |
| 0.4742 |
Statistical tests and validation workflow
Wald tests (from summary)
summary.brsmm() reports Wald \(z\)-tests for each parameter: \[
z_k = \hat\theta_k / \mathrm{SE}(\hat\theta_k).
\]
sm <- summary(fit_mm)
kbl10(sm$coefficients)
| (Intercept) |
0.3683 |
0.1544 |
2.3860 |
0.017 |
-0.1594 |
0.0846 |
-1.8829 |
0.0597 |
-0.7973 |
0.2925 |
-2.7259 |
0.0064 |
| x1 |
0.6330 |
0.0947 |
6.6877 |
0.000 |
-0.1594 |
0.0846 |
-1.8829 |
0.0597 |
-0.7973 |
0.2925 |
-2.7259 |
0.0064 |
Evolutionary scheme and Likelihood Ratio (LR) test selection
A practical workflow of increasing complexity:
brs(): no random effect (ignores clustering);brsmm(..., random = ~ 1 | id): random intercept;brsmm(..., random = ~ 1 + x1 | id): random intercept +
slope.
In the first jump (brs to brsmm with
intercept), the hypothesis \(\sigma_b^2 =
0\) lies on the boundary of the parameter space. Thus, the
classical asymptotic \(\chi^2\)
reference distribution should be interpreted with caution. In the second
jump (intercept to intercept + slope), the Likelihood Ratio (LR) test
with a \(\chi^2\) distribution is
commonly used as a practical diagnostic for goodness-of-fit gains.
# Base model without a random effect
fit_brs <- brs(
y ~ x1,
data = sim$data,
repar = 2
)
# Reuse the mixed models already fitted:
# fit_mm : random = ~ 1 | id
# fit_mm_rs : random = ~ 1 + x1 | id
tab_lr <- anova(fit_brs, fit_mm, fit_mm_rs, test = "Chisq")
kbl10(
data.frame(model = rownames(tab_lr), tab_lr, row.names = NULL),
digits = 4
)
| M1 (brs) |
3 |
-1014.913 |
2035.826 |
2046.268 |
NA |
NA |
NA |
| M2 (brsmm) |
4 |
-1008.247 |
2024.493 |
2038.416 |
13.3331 |
1 |
0.0003 |
| M3 (brsmm) |
6 |
-1007.065 |
2026.131 |
2047.015 |
2.3626 |
2 |
0.3069 |
Operational decision rule (analytical):
- If the second jump (RI to RI+RS) does not improve the fit (high
p-value), prefer the random-intercept model for parsimony.
- If there is a robust gain, adopt the RI+RS model and validate
parameter stability (especially
sd_b and the \(D\) matrix) via sensitivity and residual
diagnostics.
Residual diagnostics (quick checks)
r <- residuals(fit_mm, type = "pearson")
kbl10(
data.frame(
mean = mean(r),
sd = stats::sd(r),
q025 = as.numeric(stats::quantile(r, 0.025)),
q975 = as.numeric(stats::quantile(r, 0.975))
),
digits = 4
)
| 0.0305 |
0.9788 |
-1.8421 |
1.5132 |
Parameter recovery experiment
A single-fit recovery table can be produced directly from the
previous fit:
est <- c(
beta0 = unname(coef(fit_mm, model = "mean")[1]),
beta1 = unname(coef(fit_mm, model = "mean")[2]),
sigma_b = unname(exp(coef(fit_mm, model = "random")))
)
true <- c(
beta0 = sim$truth$beta[1],
beta1 = sim$truth$beta[2],
sigma_b = sim$truth$sigma_b
)
recovery_table <- data.frame(
parameter = names(true),
true = as.numeric(true),
estimate = as.numeric(est[names(true)]),
bias = as.numeric(est[names(true)] - true)
)
kbl10(recovery_table)
| beta0 |
0.20 |
0.3683 |
0.1683 |
| beta1 |
0.65 |
0.6330 |
-0.0170 |
| sigma_b |
0.55 |
0.4505 |
-0.0995 |
For a Monte Carlo recovery study, repeat simulation and fitting
across replicates:
mc_recovery <- function(R = 50L, seed = 7001L) {
set.seed(seed)
out <- vector("list", R)
for (r in seq_len(R)) {
sim_r <- sim_brsmm_data(seed = seed + r)
fit_r <- brsmm(
y ~ x1,
random = ~ 1 | id,
data = sim_r$data,
repar = 2,
int_method = "laplace",
method = "BFGS",
control = list(maxit = 1000)
)
out[[r]] <- c(
beta0 = unname(coef(fit_r, model = "mean")[1]),
beta1 = unname(coef(fit_r, model = "mean")[2]),
sigma_b = unname(exp(coef(fit_r, model = "random")))
)
}
est <- do.call(rbind, out)
truth <- c(beta0 = 0.20, beta1 = 0.65, sigma_b = 0.55)
data.frame(
parameter = colnames(est),
truth = as.numeric(truth[colnames(est)]),
mean_est = colMeans(est),
bias = colMeans(est) - truth[colnames(est)],
rmse = sqrt(colMeans((sweep(est, 2, truth[colnames(est)], "-"))^2))
)
}
kbl10(mc_recovery(R = 50))
References
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.
Pinheiro, J. C. and Bates, D. M. (2000). Mixed-Effects Models in
S and S-PLUS. Springer. DOI: 10.1007/b98882. Validated online via:
https://doi.org/10.1007/b98882.
Rue, H., Martino, S., and Chopin, N. (2009). Approximate Bayesian
inference for latent Gaussian models by using integrated nested Laplace
approximations. Journal of the Royal Statistical Society: Series
B, 71(2), 319-392. DOI: 10.1111/j.1467-9868.2008.00700.x. Validated
online via: https://doi.org/10.1111/j.1467-9868.2008.00700.x.