# 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?

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.

• Thanks for the response but I think I'm a bit confused now. Could you clarify what you mean when you say in these contexts, MAP regression is most commonly used? Also, how does cross validation determine $\lambda$? – TSP Sep 1 '15 at 20:07
• I have edited my post, with some very concrete information. Let me know if you still don't understand. – jlimahaverford Sep 3 '15 at 19:47
• Thanks so much! I'm probably not understanding some fundamentals but here goes. So when running the HBM model on the training data, we will still obtain parameters, $w_train$. We then, using only the training data, and the training parameters solve for a value of $\lambda$? I'm still confused after this step, assuming this is right. Shouldn't we already have estimates of $\tau^2$ and $\sigma^2$ from running the model? What exactly are we plugging in to the log-loss term? I hope this question made sense. – TSP Sep 4 '15 at 0:01
• *For a fixed $\lambda$, we train by finding the $w$ that will minimize $f$. For different values of $\lambda$ we will get different values of $w$. How do we know which value $\lambda, w$ is best? By using $w$ to make predictions. The $w$ that makes the best predictions corresponds to the best $\lambda$. How do we measure the goodness of predictions. With the loss function (the term in $f$ that does not involve $\lambda$). The lower the loss function the better the predictions. – jlimahaverford Sep 6 '15 at 17:41

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.

• Great answer. The one comment I have is on the "If MAP, why work Bayesian at all?" While one of the major advantages of Bayesian Inference is obtaining a distribution for your parameters rather than a point estimate, having a reasonable way of introducing prior beliefs into statistical inference in itself is useful, and MAP is a highly utilized and general regression technique. – jlimahaverford Sep 28 '15 at 18:09
• Personally, I'm happy with the MAP - I guess the background of my comment is that in my field (ecological statistics), people typically only move to Bayesian analysis if the MLE becomes unstable, because a Bayesian analysis requires a justification for the frequentist reviewers who ask you why you didn't simply use the MLE. And in the end, with uninformative priors, MLE and MAP don't differ that much. – Florian Hartig Sep 29 '15 at 7:49