1 Introduction↩︎

In regression models with categorical predictors, parameter identifiability issues arise when all K levels of a factor are included. Consider a simple example where gender is included in a regression model: \[y_i = \beta_0 + \beta_{M} \text{Male}_i + \beta_{F}\text{Female}_i + ... + \epsilon_i.\] Without additional constraints, this model is overparameterized since knowing one gender category perfectly predicts the other, leading to multicollinearity. Traditional approaches like treatment coding (setting one level as reference) work but have interpretability disadvantages - the intercept represents the mean for the reference group rather than the population average, and coefficients represent deviations from the reference group rather than from the overall mean. This complicates interpretation, especially with complex models involving multiple categorical variables or when no natural reference category exists.

A more interpretable approach retains all K coefficients while imposing the constraint that \(\sum_{k=1}^K \beta_k = 0\) (sum-to-zero constraint), allowing \(\beta_0\) to represent the population average and each categorical coefficient to represent the deviation from this average. However, implementing such constraints in Bayesian settings is challenging. One method redefines the last coefficient as \(\beta_K = -\sum_{k=1}^{K-1} \beta_k\), which creates undesirable asymmetry among coefficients, complicates prior specification, and leads to difficult posterior geometry. Another approach implements "soft constraints" through highly informative priors (e.g., \(\sum_k \beta_k \sim \mathcal{N}(0, \varepsilon)\) with very small \(\varepsilon\)), but this approximation can lead to numerical instabilities. We propose that the optimal solution is to directly incorporate linear constraints into prior distributions, preserving both identifiability and interpretability while enabling proper uncertainty quantification. Our approach allows for multiple arbitrary linear constraints on parameter vectors and includes efficient sampling algorithms for implementation in probabilistic programming languages. For common zero-sum constrained shrinkage priors in statistical practice, we provide four specializations as special cases of our proposed method: Bayesian ridge, hierarchical Bayesian ridge, horseshoe prior, and horseshoe prior variants. These specializations provide practitioners with ready-to-use solutions.

2 Linear constraints via priors↩︎

In this section, we develop a principled approach for incorporating linear constraints into Bayesian regression models with multivariate normal priors.

Given coefficient vector \(\boldsymbol{\beta} = (\beta_1,\beta_2,...,\beta_K)^\top \in \mathbb{R}^K\), such that \[\beta_k \mid \lambda_k \sim \mathcal{N}\left( 0, \lambda_k^2 \right), ~ k=1,2,...,K\] where \(\lambda_k\) can be either a fixed hyperparameter or a random variable that forms a Gaussian scale mixture.

Consider a linear constraint \(A \boldsymbol{\beta} = b\), with constraint matrix \(A\in \mathbb{R}^{J \times K}\) (\(1\leq J \leq K-1\)) with full row rank and constant vector \(b \in \mathbb{R}^{K}\). We propose the following prior: \[\begin{align}\label{eq46multinormal95cond} \boldsymbol{\beta} \sim \mathcal{N}\left(D A^\top \left( A D A^\top \right)^{-1}b, \; D - D A^\top \left( A D A^\top \right)^{-1} A D \right) \end{align}\tag{1}\] where \(D = \text{diag}\left( \lambda_1^2, \lambda_2^2,...,\lambda_K^2 \right)\).

By adjusting the covariance structure in this manner, we effectively reduce the parameter space to values that almost surely satisfy the constraint. This approach implements hard parameter constraints solely through the prior distribution. Although this conditional method does modify the original prior distribution, we can largely preserve the desired marginal distributions by appropriately scaling the diagonal elements of the covariance matrix, minimizing distortion to the original model specification. The specific scaling factors and their effects on the resulting posterior distributions will be explored in detail in Section 4, along with practical guidelines for implementation in various modeling contexts.

3 Sampling algorithm↩︎

The the covariance matrix (Eq. 2 ) we derived is singular (see Lemma 2). Therefore, the multivariate normal implementations in most mathematical libraries cannot directly sample from it. This section gives a general sampling algorithm and Stan implementation.

We want to sample from \[\boldsymbol{\beta} \sim \mathcal{N}\Bigl(m^*,\,\Sigma^*\Bigr)\] with \[m^* = D A^\top (A D A^\top)^{-1} b,\quad \Sigma^* = D - D A^\top (A D A^\top)^{-1} A D,\] knowing that \(\Sigma^*\) is singular with rank \(K-J\), where \(J=\operatorname{rowrank}(A)\), \(1\leq J \leq K-1\) (see Lemma 3).

The key idea is to sample in the full‐rank subspace where \(\Sigma^*\) is nondegenerate and then map the sample back to \(\mathbb{R}^K\). One common approach is to obtain an orthonormal basis for the null space of \(A\). We formulate the procedure below.

Therefore, the above algorithm first simplifies the sampling problem to a \(K-J\) subspace where the covariance is non-degenerate, and then lifts the samples to \(\mathbb{R}^K\). This achieves the construction of samples from the degenerate conditional Gaussian distribution while ensuring that the constraint \(A\boldsymbol{\beta} = b\) is satisfied.

4 Sum-to-zero specialization↩︎

In this section, we demonstrate how the covariance structure can be tailored to different modeling needs, especially sum-to-zero constraints.

Note that in GP regression, the constraint matrix \(A\) only required to have full row rank. Then it is possiable to set \(A = (1,1,...,1) \in\mathbb{R}^{1\times K}\) and \(b=0 \in \mathbb{R}\). Therefore, the constrain becomes \[\sum_{k=1}^K \beta_k = 0.\] i.e., sum-to-zero constraints.

We present four specialization of our multivariate normal structure for the commonly encountered sum-to-zero constraint in statistical practice: Bayesian ridge, herarchical Bayesian ridge, horseshoe prior, and variants of the horseshoe prior. For each prior type, we derive the explicit form of the covariance matrix, discuss the implications on mathematical details, and provide efficient computational approaches that can be directly implemented in probabilistic programming languages. These specializations provide practitioners with ready-to-use solutions for imposing sum-to-zero constraints while maintaining the desirable shrinkage properties of each prior family.

4.1 Bayesian ridge↩︎

Bayesian ridge regression uses a normal prior which imposes a global \(\ell_2\) regularization: \[\boldsymbol{\beta} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{I})\] where \(\boldsymbol{I}\) denotes the identity matrix. By choosing \(D=\boldsymbol{I}\) in Eq. 1 , the sum-to-zero constraint can be achieved by assuming the following prior \[\boldsymbol{\beta} \sim \mathcal{N}\left( \boldsymbol{0},~ \sigma^2\left( \boldsymbol{I}-\frac{1}{K}\boldsymbol{1}\boldsymbol{1}^\intercal \right) \right).\] \(\sigma^2 = \frac{K}{K-1}\) is chosen to ensure the diagonal element of covariance matrix to be \(\Sigma_{k,k} = \sigma^2 \left( 1 - \frac{1}{N} \right) = 1.\) So the marginal distribution are therefore preserved: \[\beta_k \sim \mathcal{N}(0, 1), k=1,2,...,K.\]

4.2 Hierarchical Bayesian ridge↩︎

It’s possiable to introducing an additional scale factor controling global shrinkage level. \[\boldsymbol{\beta} \sim \mathcal{N}\left( \boldsymbol{0}, \lambda^2 \boldsymbol{I} \right)\] The samller the \(\lambda\), the stronger the \(\ell_2\) penalty is. It is also a common choice to introduce partial pooling for variance parameters allowing adaptation to variability from the data \[\lambda \sim \text{Cauchy}^+(0,1).\] or any other distribution with support \(\mathbb{R}^+\).

By choosing \(D=\lambda^2\boldsymbol{I}\) in Eq. 1 , the sum-to-zero constraint can be achieved by assuming the following prior \[\boldsymbol{\beta} \sim \mathcal{N}\left( \boldsymbol{0},~ \lambda^{2}\sigma^2\left( \boldsymbol{I}-\frac{1}{K}\boldsymbol{1}\boldsymbol{1}^\intercal \right) \right).\] \(\sigma^2 = \frac{N}{N-1}\) is chosen so the marginal distribution are preserved: \[\beta_k \sim \mathcal{N}(0, \lambda^{2}), k=1,2,...,K.\]

When sampling, only need to carry out non-centered parameterization, \[\boldsymbol{\beta} = \lambda \boldsymbol{\beta^*}\] and then use Algorithm 2 to sample \(\boldsymbol{\beta^*}\).

4.3 Horseshoe prior↩︎

The Horseshoe prior [1] adopts a local-global hierarchical structure for adaptive shrinkage in regression analysis. For each regression coefficient, \[\begin{align} \beta_k &= \mathcal{N}(0, \tau^2\lambda_k^2), k=1,2,...,K \\ \lambda_j &\sim C^+(0, 1) \\ \tau &\sim C^+(0, 1) \end{align}\] The global scale parameter \(\tau\) controls the overall shrinkage level, while the local scale parameter \(\lambda_j\) adjusts the individual shrinkage for each coefficient.

By choosing \(D=\tau^2 \text{diag}\left( \lambda_1^2, \lambda_2^2,...,\lambda_K^2 \right)\) in Eq. 1 , the sum-to-zero constraint can be achieved by assuming the following prior \[\begin{align} \pmb{\beta} \mid D \sim \mathcal{N}\left(0, \Sigma\right), ~ \Sigma=D - D \mathbf{1} \left( \mathbf{1}^\top D \mathbf{1} \right)^{-1} \mathbf{1}^\top D \end{align}\]

The constraint of sum to zero introduces negative correlations among \(\beta_k\), which are reflected in the off-diagonal elements of the adjusted covariance matrix \(\tilde{\Sigma}\). The covariance between \(\beta_k\) and \(\beta_j\) is given by: \[\operatorname{Cov}(\beta_k, \beta_j) = -\tau^2 \frac{\lambda_k^2 \lambda_j^2}{\sum_{l=1}^K \lambda_l^2}, \quad j \neq k\]

The marginal variance of each \(\beta_k\) is: \[\operatorname{Var}(\beta_k) = \tau^2 \left( \lambda_k^2 - \frac{\lambda_k^4}{\sum_{l=1}^K \lambda_l^2} \right)\]

Compared to the desired \(\tau^2 \lambda_k^2\), this variance is slightly reduced due to the presence of the constraint. Therefore, in the constrained covariance matrix, the marginal variance could be compensated by multiplying \(D\) by the factor \(K / (K - 1)\), which assumes that all \(\lambda_k\) are equal.

It can also be used to first separate the global shrinkage parameter \(\tau\) by non-central parameterization before running the algorithm.

4.4 Horseshoe-like prior↩︎

The proposed method is also applicable to variants of the horseshoe prior. The regularized horseshoe prior (RHS) [2] is a class of hierarchical priors designed to address sparsity in high dimension coefficients. In contrast to the traditional horseshoe prior, RHS improves model robustness by imposing moderate regularization on large coefficients through adjustments to the slab width. RHS is formulated as \[\begin{align} \beta_k \mid \zeta_k, c, \tau &\sim \mathcal{N}\left( 0, \tau^2 \tilde{\zeta}_k^2 \right), \quad k=1,2,...,K \\ \zeta_k &\sim \text{Student-}t^+_{\nu_1} (0,1) \\ c^2 &\sim \text{Inv-Gamma}(\nu_2/2, \nu_2s^2/2) \\ \tau &\sim \text{Student-}t^+_{\nu_3} (0, \tau_0) \end{align}\] where \(\tilde{\zeta}_k^2=\frac{c^2 \zeta_k^2}{c^2 + \tau^2 \zeta_k^2}\), and \(\tau_0 = \frac{p_0}{K-p_0} \frac{\tilde{\sigma}}{\sqrt{n}}\). \(p_0\in \left\{ 1,...,K-1 \right\}\) is a hyperparameter describing the prior belief on effective number of non-zero coefficients. \(\tilde{\sigma}^2\) is the pseudo variance define by the likelihood and link function.

By choosing \(D=\tau^2 \text{diag}\left( \tilde{\zeta}_1^2, \tilde{\zeta}_2^2,...,\tilde{\zeta}_K^2 \right)\) in Eq. 1 , the sum-to-zero constraint can be achieved by assuming the following prior \[\begin{align} \pmb{\beta} \mid D \sim \mathcal{N}\left(0, \Sigma\right), ~ \Sigma=D - D \mathbf{1} \left( \mathbf{1}^\top D \mathbf{1} \right)^{-1} \mathbf{1}^\top D. \end{align}\]

When sampling, only need to carry out non-centered parameterization, \[\boldsymbol{\beta} = \tau \boldsymbol{\beta^*}\] and then use Algorithm 3 to sample \(\boldsymbol{\beta^*}\).

5 Supplementary materials↩︎

5.1 Supplementary proofs↩︎

5.2 Supplementary algorithms↩︎

References↩︎