Bayesian model selection and credible interval

Question

I have a dataset with three variables, where all variables are quantitatives. Let call it $y$, $x_1$ and $x_2$. I'm fitting a regression model in a Bayesian perspective via MCMC with rjags

I done a exploratory analysis and the scatterplot of $y\times x_2$ suggest that a quadratic term should be used. Then I fitted two models

(1) $y=\beta_0+\beta_1*x_1+\beta_2*x_2$

(2) $y=\beta_0+\beta_1*x1+\beta_2*x_2+\beta_3*x_1x_2+\beta_4*x_1^2+\beta_5*x_2^2$

In model 1 the effect size of each parameter is not small and the 95% credible interval not contains the value $0$.

In model 2 the effect size of parameters $\beta_3$ and $\beta_4$ are small and each of credible intervals for all parameters contains $0$.

The fact that a credible interval contains $0$ is enough to say that the parameter is not significant?

Then I adjusted the following model

(3)$y=\beta_0+\beta_1*x_1+\beta_2*x_2+\beta_3*x^2_2$

The effect size of each parameter is not small, but with exception of $\beta_1$ all credible intervals contains $0$.

Which is the right way to do variable selection in Bayesian statistics?

EDIT: I can use Lasso in any regression model, like Beta model? I'm using a model with variable dispersion where $$log(\sigma)=-\pmb{\delta}X$$ where $\pmb{\delta}$ is a vector. I should use Laplace prior in $\pmb{\delta}$ too?

EDIT2: I fitted two models, one with Gaussian priori for $\beta_j$, $\delta_j$ and one with Laplace(double-exponential).

The estimatives for the Gaussian model are

            Mean      SD  Naive SE Time-series SE
B[1]     -1.17767 0.07112 0.0007497      0.0007498
B[2]     -0.15624 0.03916 0.0004128      0.0004249
B[3]      0.15600 0.05500 0.0005797      0.0005889
B[4]      0.07682 0.04720 0.0004975      0.0005209
delta[1] -3.42286 0.32934 0.0034715      0.0034712
delta[2]  0.06329 0.27480 0.0028966      0.0028969
delta[3]  1.06856 0.34547 0.0036416      0.0036202
delta[4] -0.32392 0.26944 0.0028401      0.0028138

The estimatives for the Lasso model are

              Mean      SD  Naive SE Time-series SE
B[1]     -1.143644 0.07040 0.0007421      0.0007422
B[2]     -0.160541 0.05341 0.0005630      0.0005631
B[3]      0.137026 0.05642 0.0005947      0.0005897
B[4]      0.046538 0.04770 0.0005028      0.0005134
delta[1] -3.569151 0.27840 0.0029346      0.0029575
delta[2] -0.004544 0.15920 0.0016781      0.0016786
delta[3]  0.411220 0.33422 0.0035230      0.0035629
delta[4] -0.034870 0.16225 0.0017103      0.0017103
lambda    7.269359 5.45714 0.0575233      0.0592808

The estimatives for $\delta_2$ and $\delta_4$ reduced a lot in Lasso model, it means that I should remove this variables from the model?

EDIT3: The model with double exponential prior (Lasso) gaves me bigger Deviance, BIC and DIC values than the model with Gaussian priors and I even get a smaller values after removing the dispersion coefficient $\delta_2$ in the Gaussian model.

Answer 1

It is well known that building a model based on what is significant (or some other criterion such as AIC, whether a credible interval contains 0 etc.) is pretty problematic, particularly if you then do inference as if you had not done model building. Doing a Bayesian analysis does not change that (see also https://stats.stackexchange.com/a/201931/86652). I.e. you should not do variable selection, but rather model averaging (or something that could get you some zero coefficients, but reflects the whole modelling process, such as LASSO or elastic net).

Bayesian model choice is more typically framed as Bayesian model averaging. You have different models, each with a different prior probability. If the posterior model probability for a model becomes low enough, you are essentially entirely discarding the model. For equal prior weights for each model and flat priors, model averaging with weights proportional to $\exp(-\text{BIC}/2)$ for each model approximates this.

You can alternatively express the model averaging as a prior that is a mixture between a point mass (the weight of the point mass is the prior probability of the effect being exactly zero = the effect is not in the model) and a continuous distribution (e.g. spike-and-slab priors). MCMC sampling can be quite difficult for such a prior.

Carvalho et al. motivate the horseshoe shrinkage prior by suggesting that it works like a continuous approximation to a spike-and-slab prior. It is also a case of embedding the problem in a hierarchical model, where to some extent the size and presence of effects on some variables relax the required evidence for others a bit (through the global shrinkage parameter, this is a bit like false-discovery rate control) and on the other hand allow individual effects to stand on their own if the evidence is clear enough. There is a convenient implementation of it available from the brms R package that builds on Stan/rstan. There are a number of further similar priors such as the horseshoe+ prior and the whole topic is an area of ongoing research.

Answer 2

There are a number of formal methods for Bayesian variable selection. A slightly outdated review of Bayesian variable selection methods is presented in:

A review of Bayesian variable selection methods: what, how and which

A more recent review, which also includes a comparison of different methods and the performance of R packages where they are implemented is:

Methods and Tools for Bayesian Variable Selection and Model Averaging in Univariate Linear Regression

This reference is particularly useful in that it points you to specific R packages where you just need to plug in the response and the covariate values (and in some cases the hyperparameter values) in order to run the variable selection.

Another, quick and dirty and non-recommended, way of conducting "Bayesian" variable selection is to use stepwise selection (forward, backward, both) using BIC and the R command stepAIC(), which can be tweaked to perform selection in terms of BIC.

https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/stepAIC.html

Another quick and dirty way of testing $\beta_4=0$ is by using the Savage-Dickey density ratio and the posterior simulation you already got:

https://arxiv.org/pdf/0910.1452.pdf

Answer 3

The whole idea of Bayesian statistics is different from a frequentist approach. In this way I think to use the terms of significance is not accurate. I guess it is up to the reader to decide if the results (distribution) you get from your model for your $β$ 's are for him reliable or trustful. It always depends on the distribution itself. How skewed and wide is it and how much of the area is below zero?

You can also find a nice lecture about the topic here at 41:55 :

https://vimeo.com/14553953