# Computationally verifying the equivalence of ridge regression estimates and Bayesian regression estimates

I'm trying to show that the numerical estimates of ridge regression's parameter estimates are the same as the MAP parameter estimates of a Bayesian regression model with normal prior distributions. So far I have not been able to achieve numerical agreement.

I'm just toying with the mtcars data.

Ridge estimates:

> library(MASS)
> fit_ridge <- lm.ridge(mpg ~ cyl + hp + disp + wt,
+ data = mtcars,
+ lambda = 1)
> fit_ridge
cyl           hp         disp           wt
39.133455180 -1.070747007 -0.019314307  0.003855685 -3.255577449


Bayesian model fit:

> library(rstanarm)
> fit_bayes <- stan_glm(mpg ~ cyl + hp + disp + wt,
+                       family = "gaussian",
+                       data = mtcars,
+                       prior = normal(0, 1),
+                       algorithm = "optimizing",
+                       prior_aux = NULL,
+                       prior_intercept = NULL)
+ > summary(fit_bayes)

Model Info:

function:     stan_glm
family:       gaussian [identity]
formula:      mpg ~ cyl + hp + disp + wt
algorithm:    optimizing
priors:       see help('prior_summary')
observations: 32
predictors:   5

Estimates:
(Intercept) 40.3    2.5   35.5 45.2
cyl         -1.2    0.6   -2.3 -0.1
hp           0.0    0.0    0.0  0.0
disp         0.0    0.0    0.0  0.0
wt          -3.7    0.9   -5.5 -2.0
sigma        2.3    0.3    1.8  3.0


Clearly there is some disagreement between the parameter estimates. Is my model specification incorrect? Or my understanding that these distinct estimators ought to be numerically the same?

• those bayes estimates are medians, not maxima.... – probabilityislogic May 12 at 23:45
• Are they? I specified algorithm = "optimizing" which I thought yields the posterior modes. mc-stan.org/rstanarm/reference/rstanarm-package.html – JTH May 13 at 0:14
• wondering if the ridge is also penalizing the intercept, but you set prior_intercept to null - this means flat prior for intercept right? – probabilityislogic May 13 at 9:48
• Correct, flat priors on the intercept and sigma in the Bayesian model. – JTH May 13 at 14:08
• is the ridge regression not ridging the intercept? – probabilityislogic May 14 at 8:37

I've been able to show numerical agreement using stan (instead of rstanarm) and my own matrix math (instead of lm.ridge or glmnet). The stan model is

$$y \sim N(X\beta, \sigma^2_e)$$ $$\beta \sim N(0, 1/\lambda)$$

which is implemented with the stan code:

data {
int<lower=0> n;
int<lower=0> p;
matrix[n, p] mat;
vector[n] y;
real<lower=0> lambda;
real<lower=0> sigma_e;
}

parameters {
vector[p] beta;
}

model {
// ridge prior on beta
beta ~ normal(0,  sqrt(1/(lambda)));
// likelihood model
y ~ normal(mat * beta, sigma_e);
}


The MAP is given by :

out <- optimizing(stan_fit,
data = list(mat = Xmat,
n = nrow(Xmat),
p = ncol(Xmat),
lambda = lambda,
sigma_e = 1,
y = y))


where Xmat is the fixed effects design matrix and y is a vector of responses.

The ridge estimator is given by

ridge_beta <- solve(t(Xmat) %*% Xmat + lambda * diag(ncol(Xmat))) %*% t(Xmat) %*% y


I'm assuming that $$\sigma_e$$ is a known constant.

In the examples I've tried, the resulting MAPs and ridge point estimates agree.

Edit: This shrinks the intercept coefficient too.