Confusion about Cross-Validation for Hierarchical Bayesian Regression Models I had two questions regarding model selection for a Hierarchical Bayesian (HB) Regression Model and the purpose of Cross-Validation.  
1).  I understand cross-validation as one way to perform model selection for a Hierarchical Bayesian Regression Model.  Is this understanding correct?
2).  If cross-validation can be used for model selection, would one go about using it in the following way?  Say you have 5 possible variables (x1, x2, x3, x4, x5), which I believe means that you have 2^5 possible models.  Do you run your HB regression assuming a particular model (e.g., x1 + x2), find the regression parameters obtained, cross-validate to find the groupwise or overall predictive error, and then choose the model that gives the minimum predictive error?
 A: I'm not disagreeing with the answer provided by jlimahaverford, but I feel a bit wider explanation would be useful.
The context is variable selection in (hierarchical) regression models estimated by Bayesian inference. In this context I would say that:


*

*No, it is actually rather uncommon to do a direct variable selection via CV. I guess you could do it without great harm, but the issue is a) it is not very Bayesian and therefore potentially difficult to justify b) CV on the posterior or the MAP? If MAP, why work Bayesian at all? Plus, may be difficult to get stable MAP estimates for hierarchical models c) most importantly: it is probably rather slow.

*IF you can really run the analysis on all possible models, the fractional Bayes factor, which is in some sense very similar to CV, seems preferable to me, because it has a straighforward Bayesian interpretation. See O'Hagan, A. (1995) Fractional Bayes Factors for Model Comparison. J. Roy. Stat. Soc. B Met., 57, 99-138.

*If this is not the case, you can use regularization such as the Bayesian ridge or lasso, which is what jlimahaverford refers to. In this case, you don't select variables, you use the full model, but with a penalty on parameter estimates different from 0 encoded in the prior. The Penalty may or may not be optimized by CV or some other method, see, e.g. Park, T. & Casella, G. (2008) The Bayesian Lasso. J. Am. Stat. Assoc., 103, 681-686.

*A related regularization approach that is often discussed (also here on CV) are spike-and-slab priors. See Ishwaran, H. & Rao, J. S. (2005) Spike and Slab Variable Selection: Frequentist and Bayesian Strategies. The Annals of Statistics, 33, pp. 730-773. I think Fabian Scheipl http://www.statistik.lmu.de/~scheipl/research.html has also done some interesting work in this area, maybe also regarding the comparison to the Bayesian lasso.
A: 1) Yes
2) No (or Maybe, depending on what you mean).
Hierarchical Bayesian Regression Models are parametric, generative models.  Regression is the technique of determining the parameters.  The fact that the models are Bayesian implies that (in the model) the parameters of the lowest level data generating distributions are themselves drawn from higher-level (prior) distributions.  Some of the parameters of these distributions are typically referred to as hyper-parameters.
In the regression procedure these hyper-parameters are held fixed, and cross validation is used to determine optimal values for these hyper-parameters.  The most common form of estimation used in these contexts is Maximum a Posteriori (MAP) regression.  In these instances, many commonly used priors lead to relatively simple terms in the objective function.  Adding these terms to the optimization function is called regularization.  In the following examples $w$ is the vector of parameters in the model.


*

*Normal Priors with fixed Variance -> L2 regularization (add $\lambda * \|w \|_2$) to the objective function.

*Laplacian Priors with fixed Variance -> L1 regularization (add $\lambda * \|w \|_1$) to the objective function.


Cross validation is used to determine $\lambda$.  Now, back to the "maybe."  L1 regularization is well known for the fact that it tends to lead to sparse parameter vectors ($w_i = 0$ for many $i$).  In theory, one would love to try regressing with every possible subset of features, testing the results on a validation set.  Unfortunately, this is computationally very challenging (impossible in contexts with thousands of variables).
I hope this answers your question and gives you the terminology necessary to find any other information needed.
EDIT
To be clear, there are two issues here.  Setting up the (parametric) model, and estimating the parameters.  Consider this very simple model
$$
w_i \sim \mathcal{N}(0, \sigma^2), i \in \{1, \dots, n \},
$$
$$
y_{i,j} \sim \mathcal{N}(w_i, \tau^2), j \in \{1, \dots, n_i \}.
$$
Our data are the "y"s and we know their corresponding "i"s.  So we've drawn many samples from normal distributions, whose means were themselves drawn from one normal distribution.  I will refer to our data as $Y$ and to some vector of possible "w"s as $w$.  Our goal is to find probable values for $w$, given $Y$.
$$
P(w | Y) = \frac{P(Y|w) \cdot P(w)}{P(Y)}.
$$
The numerator is easy to compute from our model, as we will see.  The denominator is more difficult to compute (in general) because we would have to marginalize $w$ out of the joint distribution (for which we have an equation).  But since the denominator is always positive we can maximize the whole thing by maximizing the numerator, or the log of the numerator (since log is monotone increasing).  We will actually minimize the negative log (convention).  So,
\begin{eqnarray}
P(w | Y) &\propto& P(Y|w) \cdot P(w) \\
&\propto& \left[\prod_i\prod_j exp\left(-\frac{(y_{i,j} - w_i)^2}{2\tau^2}\right) \right] \cdot \prod_i exp\left(-\frac{w_i^2}{2\sigma^2}\right).
\end{eqnarray}
Taking the negative log, multiplying by $2\tau^2$, and thinking of it as a function of $w$, we have
\begin{eqnarray}
f(w) &=& \left[ \Sigma_i \Sigma_j (y_{i,j} - w_i)^2 \right] + \frac{\tau^2}{\sigma^2} \cdot \left[ \Sigma_i w_i^2 \right] \\
&=&  \left[ \Sigma_i \Sigma_j (y_{i,j} - w_i)^2 \right] + \lambda \cdot \left[ \Sigma_i w_i^2 \right].
\end{eqnarray}
This is what we minimize.  The first term is the log-loss and the second term is the regularization.  If we minimized only the first term out estimates would be $w_i = \Sigma_j y_{i,j} / n_{i}$.  The effect of the regularization will be to pull all of these $w_i$ toward 0.  It turned out that variances of our distributions only mattered relative to one another, so I replaced their quotient with $\lambda$.  The larger $\lambda$ is, the more $w_i$ will be pulled toward 0.  In order to determine $\lambda$ we can train our model on a subset of our data, by minimizing the below function, then evaluate our model on the remainder of our data by computing only the logloss term.
