Bayesian Hierarchical Beta Regression

• Bayesian hierarchical beta regression is a modeling framework for bounded proportions that integrates a beta likelihood with multilevel regression structures.
• It links the conditional mean and precision to fixed and random effects using logit and log links, enhancing flexibility and interpretability.
• The method employs robust prior designs and computational strategies like MCMC and INLA to accommodate complex data and mitigate boundary issues.

Bayesian hierarchical beta regression is a class of models for continuous bounded responses, especially proportions measured on a continuous scale, in which the observation model is beta and one or more beta parameters are embedded in a multilevel regression structure. In the standard mean–precision parameterization, $$Y \mid \mu,\phi \sim \mathrm{Beta}(\mu\phi,(1-\mu)\phi)$$, with $\mathbb{E}(Y)=\mu$ and $\operatorname{Var}(Y)=\mu(1-\mu)/(1+\phi)$; hierarchical formulations then link $\mu$ and, when desired, $\phi$, to fixed and random effects, typically through logit and log links (Figueroa-Zuñiga et al., 2012). This framework was developed precisely because standard Gaussian, Poisson, and Binomial GLMMs do not address the case of continuous but bounded outcomes such as rates, percentages, indexes, and proportions, and because bounds are ignored if Gaussian mixed models are used instead (Bonat et al., 2014).

1. Distributional basis and inferential target

The canonical response domain for standard beta regression is the open unit interval $(0,1)$. Under the mean–precision parameterization used in mixed beta regression, a generic observation satisfies

$$Y \mid \mu,\phi \sim \mathrm{Beta}(\mu\phi,(1-\mu)\phi),$$

with density

$$f(y\mid \mu,\phi)=\frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1-\mu)\phi)}y^{\mu\phi-1}(1-y)^{(1-\mu)\phi-1},\qquad 0<y<1.$$

The parameter $\mu$ is the conditional mean, while $\phi>0$ is a precision parameter, so larger $\mathbb{E}(Y)=\mu$0 implies smaller variance for fixed $\mathbb{E}(Y)=\mu$1 (Figueroa-Zuñiga et al., 2012). When the original response lies in $\mathbb{E}(Y)=\mu$2, the standard construction first applies the linear rescaling $\mathbb{E}(Y)=\mu$3 and then models $\mathbb{E}(Y)=\mu$4 (Figueroa-Zuñiga et al., 2012).

This formulation is especially suited to grouped or clustered bounded data. It preserves the support of the outcome, gives a mean–variance relation driven by $\mathbb{E}(Y)=\mu$5 and $\mathbb{E}(Y)=\mu$6, and supports direct regression modeling of both expected value and dispersion. A recurrent misconception is that beta regression is only a mean model. In the hierarchical literature, both the mean and the precision can be endowed with their own regression structures, including random effects (Figueroa-Zuñiga et al., 2012).

A second boundary of the classical framework is equally central: standard beta likelihoods are defined on $\mathbb{E}(Y)=\mu$7, not on $\mathbb{E}(Y)=\mu$8. If any $\mathbb{E}(Y)=\mu$9, the log-likelihood contains $\operatorname{Var}(Y)=\mu(1-\mu)/(1+\phi)$0, so standard likelihood-based and Bayesian beta regression cannot be used directly (Kim et al., 16 Sep 2025). That restriction shaped much of the subsequent development of hierarchical beta-type models.

2. Mixed-effects formulation

In the hierarchical setting, observations are indexed by group. With $\operatorname{Var}(Y)=\mu(1-\mu)/(1+\phi)$1 groups and $\operatorname{Var}(Y)=\mu(1-\mu)/(1+\phi)$2 observations in group $\operatorname{Var}(Y)=\mu(1-\mu)/(1+\phi)$3, the mean submodel is written as

$\operatorname{Var}(Y)=\mu(1-\mu)/(1+\phi)$4

where $\operatorname{Var}(Y)=\mu(1-\mu)/(1+\phi)$5 and $\operatorname{Var}(Y)=\mu(1-\mu)/(1+\phi)$6 are fixed-effect and random-effect design vectors, respectively, and $\operatorname{Var}(Y)=\mu(1-\mu)/(1+\phi)$7 is a group-specific random-effect vector (Figueroa-Zuñiga et al., 2012). Given $\operatorname{Var}(Y)=\mu(1-\mu)/(1+\phi)$8 and $\operatorname{Var}(Y)=\mu(1-\mu)/(1+\phi)$9,

$\mu$0

Two principal precision specifications appear in the mixed beta literature. In the constant-precision formulation, a single scalar $\mu$1 is shared across all observations. In the more general formulation,

$\mu$2

so the precision itself has a mixed-effects linear predictor with its own fixed effects $\mu$3 and random effects $\mu$4 (Figueroa-Zuñiga et al., 2012). The design matrices for mean and precision, $\mu$5 versus $\mu$6, may be chosen independently, which is a major source of modeling flexibility.

The random-effects law need not be Gaussian. A distinctive feature of the Bayesian mixed beta formulation is the use of multivariate Student-$\mu$7 random effects,

$\mu$8

with the Gaussian case recovered as a special or limiting case as the degrees of freedom increase (Figueroa-Zuñiga et al., 2012). This gives robustness to outliers in the random-effect space. In the latent Gaussian tradition, however, Gaussian random effects remain standard, particularly when INLA or spatial-temporal structures are used (Bonat et al., 2014).

The same basic template extends naturally to structured dependence. In the spatio-temporal obesity application, the mean predictor is

$\mu$9

with $\phi$0 a structured spatial effect and $\phi$1 a structured temporal effect; $\phi$2 follows an ICAR prior and $\phi$3 follows an AR(1) process (Rota et al., 7 Aug 2025). This places Bayesian hierarchical beta regression squarely within the broader class of latent Gaussian multilevel models.

3. Bayesian hierarchy, prior design, and regularization

The original mixed beta specification uses multivariate $\phi$4 priors for fixed effects,

$\phi$5

together with inverse-Wishart priors for random-effect covariance matrices and exponential priors for the degrees of freedom of $\phi$6-distributed random effects (Figueroa-Zuñiga et al., 2012). For a constant precision parameter $\phi$7, several prior families were studied: inverse-gamma, uniform-squared, a proposed beta-squared prior $\phi$8 with $\phi$9, and a $(0,1)$0-prior on $(0,1)$1 (Figueroa-Zuñiga et al., 2012). Sensitivity analysis in both simulated and real examples showed that the choice among these priors had limited impact on inferences about the regression structure, although DIC, EAIC, and EBIC could slightly favor particular specifications (Figueroa-Zuñiga et al., 2012).

INLA-based beta mixed models typically use diffuse Gaussian priors for fixed effects, Gamma priors for the beta precision parameter, and Wishart or Gamma constructions for Gaussian random-effect precision components (Bonat et al., 2014). In that line of work, prior elicitation for random-effect precisions is tied to the Wakefield–Fong construction, which encodes plausible ranges for latent effects through marginal $(0,1)$2-distributions (Bonat et al., 2014).

Later work broadened the prior architecture in two directions. One direction imposed linear inequality restrictions on beta-regression coefficients through truncated multivariate normal priors of the form $(0,1)$3, yielding a constrained Bayesian estimator that, in simulation, outperformed ordinary estimators and outperformed the ridge estimator in terms of the standard deviation and the mean squared error, even in the presence of multicollinearity (Seifollahi et al., 2024). The other direction targeted high-dimensional sparsity. In the spatio-temporal obesity model, Stochastic Search Variable Selection is implemented through a spike-and-slab hierarchy in which $(0,1)$4 has a small-variance spike or a larger-variance slab governed by a Bernoulli inclusion indicator $(0,1)$5 and a $(0,1)$6 prior on the inclusion probability $(0,1)$7 (Rota et al., 7 Aug 2025). In the high-dimensional sparse beta-regression framework, the coefficients satisfy

$(0,1)$8

with $(0,1)$9 and $$Y \mid \mu,\phi \sim \mathrm{Beta}(\mu\phi,(1-\mu)\phi),$$0, implemented through an inverse-Gamma expansion (Mai, 28 May 2025). There, the common precision parameter $$Y \mid \mu,\phi \sim \mathrm{Beta}(\mu\phi,(1-\mu)\phi),$$1 is treated as known in the main development (Mai, 28 May 2025).

These developments show that “hierarchical” in beta regression has come to denote more than random intercepts and slopes. It also includes hyperprior design, structural restrictions, and shrinkage hierarchies that control regularization and model size.

4. Posterior computation and model assessment

The original mixed beta formulation was implemented in BUGS/WinBUGS with MCMC, conceptually as a Gibbs sampler cycling through covariance matrices, degrees of freedom, fixed effects, precision parameters, and latent group-specific quantities (Figueroa-Zuñiga et al., 2012). Because the beta likelihood and $$Y \mid \mu,\phi \sim \mathrm{Beta}(\mu\phi,(1-\mu)\phi),$$2-random effects are generally non-conjugate, WinBUGS handles many updates internally through adaptive Metropolis, slice sampling, or other built-in mechanisms (Figueroa-Zuñiga et al., 2012). In the simulation example, 100,000 iterations were run, the first 10,000 were discarded, and the retained draws showed multivariate Gelman–Rubin $$Y \mid \mu,\phi \sim \mathrm{Beta}(\mu\phi,(1-\mu)\phi),$$3; in the Prater gasoline example, 200,000 iterations with 10,000 discarded produced $$Y \mid \mu,\phi \sim \mathrm{Beta}(\mu\phi,(1-\mu)\phi),$$4, together with stationary and unimodal posterior diagnostics (Figueroa-Zuñiga et al., 2012).

INLA provides a deterministic alternative for beta mixed models with latent Gaussian structure. In the industrial workers’ life-quality application, INLA produced posterior marginals very similar to MCMC and likelihood analysis, while being “suitable and faster” and easier to use for fitting alternative models and priors; DIC, LML, and CPO were used for model comparison (Bonat et al., 2014). In that study, a random-intercept model was preferred to a random-intercept-plus-random-slope alternative (Bonat et al., 2014).

More recent computation has exploited auxiliary-variable schemes. In the high-dimensional sparse beta model, Pólya–Gamma augmentation is used so that conditional on latent $$Y \mid \mu,\phi \sim \mathrm{Beta}(\mu\phi,(1-\mu)\phi),$$5, the likelihood becomes Gaussian in the linear predictor, yielding a multivariate normal full conditional for $$Y \mid \mu,\phi \sim \mathrm{Beta}(\mu\phi,(1-\mu)\phi),$$6 and inverse-Gamma full conditionals for the Horseshoe scale parameters (Mai, 28 May 2025). This supports a fully Gibbs sampler and underpins the paper’s posterior consistency and convergence-rate results (Mai, 28 May 2025). Custom MCMC code also appears in hierarchical SLTB regression, where Metropolis–Hastings and Gibbs steps are combined for hierarchical linear and nonlinear models on $$Y \mid \mu,\phi \sim \mathrm{Beta}(\mu\phi,(1-\mu)\phi),$$7 (Kim et al., 16 Sep 2025).

A computational direction rather than an implemented beta model appears in the neural posterior estimation framework “metabeta.” That model is trained and evaluated for linear mixed-effects regression with Gaussian observation noise, not for a beta likelihood, but its authors note that the hierarchical structure is essentially the same and that adaptation to Bayesian hierarchical beta regression would amount to replacing the Gaussian likelihood with a beta likelihood in the simulator and posterior calculations (Kipnis et al., 8 Oct 2025). This suggests an amortized route for mixed beta models when repeated inference over many related datasets is required.

5. Empirical domains and substantive interpretation

The main applied use-case is grouped or correlated proportion data: repeated measures of proportions, multilevel data with individuals nested within groups or sites, and bounded continuous outcomes for which intra-group correlation must be modeled (Figueroa-Zuñiga et al., 2012). The central interpretive advantage is the simultaneous modeling of the conditional mean and the conditional precision. In the Prater gasoline conversion data, the mean submodel used a random intercept by batch,

$$Y \mid \mu,\phi \sim \mathrm{Beta}(\mu\phi,(1-\mu)\phi),$$8

and the best precision model was

$$Y \mid \mu,\phi \sim \mathrm{Beta}(\mu\phi,(1-\mu)\phi),$$9

so higher end point increased both the expected proportion converted and the precision, meaning the response was both higher and less variable at higher temperatures (Figueroa-Zuñiga et al., 2012).

A second representative application is the life-quality index of industry workers, with companies nested in states. There the final selected beta mixed model used a random intercept for state, showing that hierarchical beta regression can capture cluster-level heterogeneity in bounded socioeconomic indices without abandoning the beta likelihood (Bonat et al., 2014). In the obesity analysis for Italian regions from 2010 to 2022, the model incorporated structured spatial and temporal random effects, gender, and exogenous predictors; the analysis found regional heterogeneity and dependence over time, and the structured spatial and temporal random effects, along with gender, emerged as the primary determinants of obesity prevalence across Italian regions, while the role of exogenous covariates was minimal at the regional level (Rota et al., 7 Aug 2025).

Hierarchical beta-type models also appear in settings where zeros are common and spatial structure matters. In ecological percent-cover data, zero-inflated beta regression distinguishes structural zeros due to unsuitability from zeros due to left censoring of a latent beta regression, and spatial random effects improved predictive performance (Tang et al., 2021). In delay discounting and alcohol-use applications, hierarchical SLTB models were used to retain a beta-like interpretation while allowing responses on the closed unit interval (Kim et al., 16 Sep 2025).

Across these domains, the same substantive pattern recurs: fixed effects govern systematic shifts in the central tendency of a bounded outcome, random effects absorb group-level or structured dependence, and precision modeling determines whether variability itself changes across covariate profiles or clusters.

6. Boundary values, robustness, and related generalizations

Classical Bayesian hierarchical beta regression is restricted to $$f(y\mid \mu,\phi)=\frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1-\mu)\phi)}y^{\mu\phi-1}(1-y)^{(1-\mu)\phi-1},\qquad 0<y<1.$$0. That restriction has led to several distinct extensions. One route is the Scale-Location-Truncated Beta model, which applies a scale–location transform of a beta variable, then truncates to $$f(y\mid \mu,\phi)=\frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1-\mu)\phi)}y^{\mu\phi-1}(1-y)^{(1-\mu)\phi-1},\qquad 0<y<1.$$1. The resulting SLTB density is nearly identical to the beta density on $$f(y\mid \mu,\phi)=\frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1-\mu)\phi)}y^{\mu\phi-1}(1-y)^{(1-\mu)\phi-1},\qquad 0<y<1.$$2 but has finite positive density at 0 and 1; the paper implements this in both hierarchical linear and hierarchical nonlinear Bayesian models (Kim et al., 16 Sep 2025). A second route is the endpoint-heterogeneous beta family, where interior observations follow a standard beta law and boundary observations use limiting beta distributions summarized through an outcome-dependent $$f(y\mid \mu,\phi)=\frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1-\mu)\phi)}y^{\mu\phi-1}(1-y)^{(1-\mu)\phi-1},\qquad 0<y<1.$$3; its hierarchical panel formulation uses random intercepts and, when needed, bounded likelihood contributions for endpoint-dominated groups (Hahn, 2023). A third route is zero-inflated beta regression, which combines a Bernoulli point mass at zero with a left-censored extended beta mechanism and can include spatial Gaussian-process random effects (Tang et al., 2021).

Robustness concerns have also generated alternatives that depart from beta regression in a stricter sense. The TNBbeta distribution is introduced as a median-based and concentration-based model on $$f(y\mid \mu,\phi)=\frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1-\mu)\phi)}y^{\mu\phi-1}(1-y)^{(1-\mu)\phi-1},\qquad 0<y<1.$$4 with a boundary parameter $$f(y\mid \mu,\phi)=\frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1-\mu)\phi)}y^{\mu\phi-1}(1-y)^{(1-\mu)\phi-1},\qquad 0<y<1.$$5; with $$f(y\mid \mu,\phi)=\frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1-\mu)\phi)}y^{\mu\phi-1}(1-y)^{(1-\mu)\phi-1},\qquad 0<y<1.$$6, it has positive finite density at 0 and 1, and Pólya–Gamma augmentation yields simple auxiliary-variable Gibbs samplers for regression models that can incorporate spatial structure and exact zero observations (Lederman et al., 10 Jun 2026). Likewise, cobin and micobin regression are presented as scalable and robust alternatives for continuous proportional data; micobin supports responses exactly at the boundary and uses Kolmogorov–Gamma data augmentation for Gibbs sampling, including in hierarchical cases involving nested, longitudinal, or spatial data (Lee et al., 21 Apr 2025).

These developments clarify two points that standard presentations can obscure. First, Bayesian hierarchical beta regression is a mature mixed-model framework for bounded continuous data, not merely a beta analogue of logistic regression. Second, the inability of the classical beta likelihood to admit exact $$f(y\mid \mu,\phi)=\frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1-\mu)\phi)}y^{\mu\phi-1}(1-y)^{(1-\mu)\phi-1},\qquad 0<y<1.$$7 or $$f(y\mid \mu,\phi)=\frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1-\mu)\phi)}y^{\mu\phi-1}(1-y)^{(1-\mu)\phi-1},\qquad 0<y<1.$$8 values does not force a single remedy. It has instead produced a family of hierarchical Bayesian models—truncated, endpoint-heterogeneous, zero-inflated, and robust alternatives—that preserve the central multilevel logic while modifying the observation law (Kim et al., 16 Sep 2025).