Easy way to do undertake simple regularised regression in rstanarm/brms?

Interfaces rstanarm
1

Hello,

I might be missing something in the documentation, but I’m trying to understand whether there is an easy way to undertake a regression (lm/glm) where all of my coefficients are subject to a common normal or student-t prior that is adapted as a free parameter.

Things I have looked at

From https://mc-stan.org/rstanarm/reference/stan_glm.html:

  • One could use eg prior=normal() or prior=student_t(), but these use fixed-width priors that are not fitted. (Even if the data is re-scaled so that the standard deviation of each column in the design matrix has a standard deviation of 1, this is still true).
  • One can use the hierarchical shrinkage families prior=hs() or prior=hs_plus() but these have geometries that are difficult to sample from, and thus can be very slow

Basically what I’m looking for

Taking the example from 1.1 Linear regression | Stan User’s Guide and modifying it, I’m really looking for something that does the following, where “THIS LINE NEW” has been added as annotation.

data {
  int<lower=0> N;   // number of data items
  int<lower=0> K;   // number of predictors
  matrix[N, K] x;   // predictor matrix
  vector[N] y;      // outcome vector
}
parameters {
  real alpha;           // intercept
  vector[K] beta;       // coefficients for predictors
  real<lower=0> sigma_beta; // THIS LINE NEW
  real<lower=0> nu_beta; // THIS LINE NEW
  real<lower=0> sigma;  // error scale
}
model {
  beta ~ student_t(nu_beta, 0, sigma_beta); // THIS LINE NEW 
  y ~ normal(x * beta + alpha, sigma);  // likelihood
}

Is there a way of easily doing this within rstanarm or brms? Or does this require a custom model like the one shown above?

Thank you!

1 Like
2

Hi,
I don’t think that’s currently straightforwardly possible. However you should be able to hack it together in brms by using stanvars (which let you insert specific pieces of Stan code into a brms model). A somewhat simpler example hopefully showing enought to let you achieve what you need:

dd <- data.frame(y = rnorm(10), x1 = rnorm(10), x2 = rnorm(10))
f <- y ~ x1 + x2
svars <- stanvar(scode = "real<lower=0> my_sigma;\n", block = "parameters") + 
                stanvar(scode = "target += student_t_lpdf(my_sigma | 3, 0, 1);\n", block = "model")

priors <- prior_string("normal(0, my_sigma)", class = "b")

# Check that Stan code is what we would expect it to be:
make_stancode(f, data = dd, prior = priors, stanvars = svars)

brm(f, data = dd, prior = priors, 
              stanvars = svars)

Best of luck with your model!