# Bayesian linear regression: How to enforce constraint on the sum of coefficients?

I have a linear regression problem in which my $$X$$ matrix is not full rank. Here is a small example:

$$X = \left[\begin{array}{rrrr} -1 & 0 & 0 & 1 \\ 1 & 0 & -1 & 0 \\ 0 & 1 & -1 & 0 \\ 0 & 0 & 1 & -1 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & -1 & 1 \\ 1 & 0 & 0 & -1 \end{array}\right]$$

As you can see, each row sum is equal to 0, thus $$\mathbb{1}$$ is in the nullspace of $$X$$, indicating the matrix is not full rank. I would like to "fix" this problem by enforcing the constraint that the sum of all 4 coefficients is equal to 0, thus guaranteeing a unique solution. I can do this in a frequentist OLS pretty easily, by using the technique described here to get a new matrix $$X'$$, with one fewer dimension, and then proceed with the normal equations.

$$X' = \left[\begin{array}{rrr} -2 & -1 & -1 \\ 1 & 0 & -1 \\ 0 & 1 & -1 \\ 1 & 1 & 2 \\ -1 & 0 & 1 \\ -1 & -1 & -2 \\ 2 & 1 & 1 \end{array}\right]$$

However, I want to use a Bayesian approach by placing weakly informative priors around each of the coefficients. These priors are normal distributions that are not zero-centered. I should note I have priors that all 4 coefficients, including the coefficient for the dropped column.

My problem is that when I do Bayesian regression with these priors on the $$X'$$ dataset, the resulting posterior distribution for the 4th "hidden" coefficient is very wide and becomes almost nonsensical. I suspect this is because most default Bayesian software (I am using brms in R, which calls to Stan) assumes the priors are independent, when that is obviously not the case in my scenario (there has to be some negative correlation). Thus when I sample from the posteriors for the three coefficients, and use those to derive samples for the "hidden" coefficient (by taking the sum and multiplying by -1, as the coefficients must sum to 0), the posterior distribution is very wide, and many, many data points are required to get a tighter distribution that is about as wide as that for the other coefficients.

This is a problem because there is nothing "special" about the column that is dropped, yet the posterior for its coefficient will be a lot wider than that for the columns that are directly modeled.

How to resolve this?

One idea that I had is that I should somehow specify the prior as the entire joint distribution, but I'm not sure how to do this in brms or other software. I did some quick research on using covariance matrices with Bayesian approach, but I am not interested in estimating the matrix itself (placing a prior around it, etc.), moreso just relaxing the starting assumption that the priors are independent.

I should also note that in my full setting (~5000 data points, ~700 columns), the joint posterior clearly does have the required correlation between the coefficients, and the derived samples for the "hidden" coefficient are sensible. This issue only occurs when I have relatively few data points.

I found this question which pretty much describes the same setting that I have, but it does not address my specific question.

• You can adapt the procedure you linked for "frequentist" OLS: do MCMC with the exact same transformed matrix, and then apply the back-transformation to each posterior draw. Jan 14 at 20:43
• That's exactly what I did, my problem is that the posterior for the back-transformed coefficient is far more spread out than it "should" be (i.e., compared to the posteriors for the other coefficients, when really the posterior variances should be roughly the same) Jan 14 at 20:49
• OK, I figured out how to pass a multivariate normal as a full joint prior, but a) the time to sample from the posterior exploded; and b) it didn't fully rectify my issue (posterior was still wider than expected, but not as wide as it had been before) Jan 15 at 14:27
• oh yeah a) is cause you've got 700 cols and using the straightforward method would involve decomposing a 700x700 matrix on each iter. For b), I think we'd have to know more about the exact mvn you used. If you're using stan, here's another thing you can try: sample $\beta$ independently, per usual. Then define a new vector $\tilde{\beta}$ whose entry $\tilde{\beta}_i = \beta_i - \texttt{mean}(\beta)$. Now track $\tilde{\beta}$ instead of $\beta$ and use $\tilde{\beta}$ in the likelihood. Jan 15 at 14:48
• Alternatively, the stan page on intrinsic CARs (mc-stan.org/users/documentation/case-studies/icar_stan.html) talks about sum zero constraints in a slightly different setting (the prior is a correlated spatial prior). They seem to "solve" the problem by specifying: // soft sum-to-zero constraint on phi sum(phi) ~ normal(0, 0.001 * N); // equivalent to mean(phi) ~ normal(0,0.001) [phi is their beta] Jan 15 at 14:50

As I hinted at in the question, and as @John Madden explicitly laid out, it is possible to create a covariance matrix and define the entire joint prior distribution as a multivariate normal. In this setting, the correlation matrix has 1 on the diagonal, and $$\frac{-1}{P-1}$$ for every other element, where $$P$$ is the original number of features ($$P = 4$$, in the example from the question).

If all prior variances are equal, it can be shown that with this correlation matrix, the variance of the sum of all elements is equal to the variance of each individual element.

The other part of my question required implementation, ideally using brms in R. Here is a possible implementation, where p is as defined above, priorMean is a vector of length $$P-1$$ containing the prior means, and priorVar is a vector of length $$P-1$$ containing the prior variances (in my case, they were all roughly equal, but slightly different):

priorCor = matrix(-1/(p-1), nrow = p - 1, ncol = p - 1)
diag(priorCor) = 1
priorVcov = diag(priorVar ^ .5) %*% priorCor %*% diag(priorVar ^ .5)

model = brm(...,
prior = set_prior("multi_normal(priorMean, priorVcov)"),
stanvars = stanvar(priorMean) + stanvar(priorVcov)
...
)


where ... represents other arguments to brm (formula, data, chains, cores, etc.)

Unfortunately, this solution can drastically increase sampling time. Originally, with independent priors and ~700 columns, the model took about 10 seconds to generate 1750 samples (500 warmup, 1250 saved), on a subset of 250 data points. On the full sample of ~5000 data points, it took about 45 minutes. After defining the prior as a (correlated) multivariate normal distribution, to generate the same number of samples, the small dataset took about 45 minutes, and the full dataset took close to 3 hours.

• shame it takes so long in this software. because your covar is diagonal plus low rank, a lightning-fast custom implementation could be written. Jan 20 at 3:50
• What is interesting is that if a specify a covariance matrix that is diagonal (all off-diagonal is zero), the sampling still takes a very long time compared to specifying the same prior as a bunch of separate independent normals. But that is discussion for a different question. Thanks for your help! Jan 23 at 2:13
• presumably it is simply blindly performing "dense" matrix decompositions Jan 23 at 3:41
• oh now that I think of it there's a way to speed it up by taking the matrix decomposition ourselves and doing independent priors in the eigenbasis of the matrix if speed is important. Jan 23 at 15:11
• Thanks. I see the same behavior when I set sample_prior="only" in the brm call. Meaning just sampling from the multivariate normal itself is hard. But what is also weird is that I have another dataset with similar construction but only 350 columns (so half the size), and it only takes Stan a couple minutes to get that many samples just from the prior (compared to like an hour). I will post on Stan forums about this Jan 24 at 15:38