IntroductionSometimes here at STATWORX, we have impromptu discussions about statistical methods. In one such discussion, one of my colleagues decided to declare (albeit jokingly) Bayesian statistics unnecessary. That made me ask myself: Why would I ever use Bayesian models in the context of a standard regression problem? Existing approaches such as Ridge Regression are just as good, if not better. However, the Bayesian approach has the advantage that it lets you regularize your model to prevent overfitting and meaningfully interpret the regularization parameters. Contrary to the usual way of looking at ridge regression, the regularization parameters are no longer abstract numbers but can be interpreted through the Bayesian paradigm as derived from prior beliefs. In this post, I’ll show you the formal similarity between a generalized ridge estimator and the Bayesian equivalent.
A (very brief) Primer on Bayesian StatsTo understand the Bayesian regression estimator, a minimal amount of knowledge about Bayesian statistics is necessary, so here’s what you need to know (if you don’t already): In Bayesian statistics, we think about model parameters (i.e., regression coefficients) probabilistically. In other words, the data given to us is fixed, and the parameters are considered random. That runs counter to the standard frequentist perspective in which the underlying model parameters are treated as fixed. At the same time, the data are considered random realizations of the stochastic process driven by those fixed model parameters. The end goal of Bayesian analysis is to find the posterior distribution, which you may remember from Bayes Rule:
While is our likelihood and is a normalizing constant, is our prior which does not depend on the data, . In classical statistics, is set to 1 (an improper reference prior) so that when the posterior ‘probability’ is maximized, really just the likelihood is maximized because it’s the only part that still depends on . However, in Bayesian statistics, we use an actual probability distribution in place of , a Normal distribution, for example. So let’s consider the case of a regression problem, and we’ll assume that our target, , and our prior follow normal distributions. That leads us to conjugate Bayesian analysis, in which we can neatly write down an equation for the posterior distribution. In many cases, this is not possible, and for this reason, Markov Chain Monte Carlo methods were invented to sample from the posterior – taking a frequentist approach, ironically. We’ll make the usual assumption about the data: is i.i.d. for all observations . This gives us our standard likelihood for the Normal distribution. Now we can specify the prior for the parameter we’re trying to estimate, . If we choose a Normal prior (conditional on the variance, ) for the vector or weights in , i.e. and an inverse-Gamma prior over the variance parameter it can be shown that the posterior distribution for is Normally distributed with mean
If you’re interested in a proof of this result check out Jackman (2009, p.526). Let’s look at it piece by piece:
- is our standard OLS estimator,
- is the mean vector of (multivariate normal) prior distribution, so it lets us specify what we think the average values of each of our model parameters are
- is the covariance matrix and contains our respective uncertainties about the model parameters. The inverse of the variance is called the precision
That would mean that we are infinitely uncertain about our prior beliefs that the mean vector of our prior distribution would vanish, contributing nothing to our posterior! Likewise, if our uncertainty decreases (and the precision thus increases), the prior mean, , would contribute more to the posterior mean. After this short primer on Bayesian statistics, we can now formally compare the Ridge estimator with the above Bayesian estimator. But first, we need to take a look at a more general version of the Ridge estimator.
Generalizing the Ridge estimatorA standard tool used in many regression problems, the standard Ridge estimator is derived by solving a least-squares problem from the following loss function:
While minimizing this gives us the standard Ridge estimator you have probably seen in textbooks on the subject, there’s a slightly more general version of this loss function:
Let’s derive the estimator by first re-writing the loss function in terms of matrices:
Differentiating with respect to the parameter vector, we end up with this expression for the gradient:
So, Minimizing over we get this expression for the generalized ridge estimator:
The standard Ridge estimator can be recovered by setting . Usually we regard as an abstract parameter that regulates the penalty size and as a vector of values (one for each predictor) that increases the loss the further these coefficients deviate from these values. When the coefficients are pulled towards zero. Let’s take a look at how the estimator behaves when the parameters, , and change. We’ll define a meaningful ‘prior’ for our example and then vary the penalty parameter. As an example, we’ll use the
diamondsdataset from the
ggplot2package and model the price as a linear function of the number of carats, in each diamond, the depth, table, x, y and z attributes As we can see from the plot, both with and without a prior, the coefficient estimates change rapidly for the first few increases in the penalty size. We also see that the ‘shrinkage’ effect holds from the upper plot: as the penalty increases, the coefficients tend towards zero, some faster than others. The plot on the right shows how the coefficients change when we set a sensible ‘prior’. The coefficients still change, but they now tend towards the ‘prior’ we specified. That’s because penalizes deviations from our , which means that larger values for the penalty pull the coefficients towards . You might be asking yourself how this compares to the Bayesian estimator. Let’s find out!
Comparing the Ridge and Bayesian EstimatorNow that we’ve seen both the Ridge and the Bayesian estimators, it’s time to compare them. We discovered, that the Bayesian estimator contains the OLS estimator. Since we know its form, let’s substitute it and see what happens:
This form makes the analogy much clearer:
- corresponds to , the matrix of precisions. In other words, since is the identity matrix, the ridge estimator assumes no covariances between the regression coefficients and a constant precision across all coefficients (recall that is a scalar)
- corresponds to , which makes sense, since the vector is the mean of our prior distribution, which essentially pulls the estimator towards it, just like ‘shrinks’ the coefficients towards its values. This ‘pull’ depends on the uncertainty captured by or in the ridge estimator.
Performance comparisonLastly, we’ll compare the predictive performance of the two models. Although we could treat the parameters in the model as hyperparameters, which we would need to tune, this would defy the purpose of using prior knowledge. Instead, let’s choose a previous specification for both models, and then compare the performance on a holdout set (30% of the data). While we can use the simple as our predictor for the Ridge model, the Bayesian model provides us with a full posterior predictive distribution, which we can sample from to get model predictions. To estimate the model I used the
|Bayesian Linear Model||1625.38||1091.36||44.15|
RecapIn this post, I’ve shown you how the ridge estimator compares to the Bayesian conjugate linear model. Now you understand the connection between the two models and how a Bayesian approach can provide a more readily interpretable way of regularizing your model. Normally would be considered a penalty size, but now it can be interpreted as a measure of prior uncertainty. Similarly, the parameter vector can be seen as a vector of prior means for our model parameters in the extended ridge model. As far as the Bayesian approach goes, we also can use prior distributions to implement expert knowledge in your estimation process. This regularizes your model and allows for incorporation of external information in your model. If you are interested in the code, check it out at our GitHub page!
- Jackman S. 2009. Bayesian Analysis for the Social Sciences. West Sussex: Wiley.