0
$\begingroup$

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:
            Median MAD_SD 2.5% 97.5%
(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?

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

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.