Several ideas could serve to justify model selection based on the ELPD instead of the model-posterior.
- Direct relation to frequentist-like MSE comparison for regression models with Gaussian likelihood and Gaussian prior. Consider the following regression models ($\sigma^2$ known):
Model $M_1$:
$$ \begin{align} y_i &\overset{iid}{\sim} N(\mu,\sigma^2), \quad i=1,\cdots,n \\
\mu &\sim N(\mu_0,\sigma^2_0).
\end{align} $$
Then, the log of the predictive distribution is:
$$
\log p(\tilde{y}\mid D,M_1) = -\frac{1}{2}\log\left[2\pi(\sigma^2+\sigma^2_1) \right] -\frac{1}{2}\frac{(\tilde{y}-\mu_1)^2}{\sigma^2+\sigma^2_1},
$$
where $\mu_1$ and $\sigma^2_1$ are the posterior mean and posterior variance of $\mu$.
Model $M_2$, with $x_i\in\mathbb{R}^d$:
$$ \begin{align} y_i\mid x_i &\overset{iid}{\sim} N(x_i^\top\beta,\sigma^2), \quad i=1,\cdots,n \\
\beta &\sim N_d(\beta_0,\Sigma_0). \end{align} $$
In this case, the log of the predictive distribution is:
$$
\log p(\tilde{y}\mid \tilde{x},D,M_2) = -\frac{1}{2}\log\left[2\pi(\sigma^2+\tilde{x}^\top\Sigma_2\tilde{x}) \right] -\frac{1}{2}\frac{(\tilde{y}-\tilde{x}^\top\beta_2)^2}{\sigma^2+\tilde{x}^\top\Sigma_2\tilde{x}},
$$
where $\beta_2$ and $\Sigma_2$ are the posterior mean and posterior covariance matrix of $\beta$.
When $n$ is large, $\sigma^2+\sigma^2_1\approx\sigma^2$ and $\sigma^2+\tilde{x}^\top\Sigma_2\tilde{x}\approx\sigma^2$. This is because, under no misspecification, the posterior variances $\sigma_1^2$ and $\Sigma_2$ shrink to zero with $n$. Thus, the ELPD of each model is approximately:
$$
\begin{aligned}
\mathbb{E}_{\tilde{y}} \log p(\tilde{y}\mid D,M_1) &\approx -\frac{1}{2}\log\left(2\pi\sigma^2 \right) -\frac{1}{2\sigma^2}\mathbb{E}_{\tilde{y}}(\tilde{y}-\mu_1)^2 = K_1-K_2^2\,\text{MSE}_1,\\
\mathbb{E}_{\tilde{y},\tilde{x}} \log p(\tilde{y}\mid\tilde{x}, D,M_2) &\approx -\frac{1}{2}\log\left(2\pi\sigma^2 \right) -\frac{1}{2\sigma^2}\mathbb{E}_{\tilde{y},\tilde{x}}(\tilde{y}-\tilde{x}^\top\beta_2)^2 = K_1-K_2^2\,\text{MSE}_2.
\end{aligned}
$$
So, in this context, and for a large sample size $n$, choosing the model with highest ELPD would mostly be equivalent to choosing the model with the lowest MSE. This is, selecting/comparing Bayesian models using the ELPD finds a direct analogue in the frequentist framework using the MSE. For other models beyond the Gaussian likelihood + Gaussian prior, the connection can be given in terms of different loss functions.
- The ELPD is directly related to several criteria from Information Theory and, in particular, to those based on Kullback-Leibler (KL) information. The KL divergence from the true data-generating distribution $p^*$ to the (marginal) predictive distribution in model $M$ is:
$$
D_{KL}(p^*\,|\!|\,p_{\text{pred}})=-\mathbb{E}_{\tilde{y}} \log p(\tilde{y}\mid D,M)-\mathbb{H}(p^*)\geq 0,
$$
where $\mathbb{H}(p^*)$ is the entropy of $p^*$. Most information criteria measures can be expressed as -2 times an estimated total ELDP, $\hat{L}=\sum_{i=1}^n\hat{L}_i$, in combination with a penalization term, $2\hat{\Lambda}(M)$. The latter is required to account for overfitting when no holdout data is available for evaluation.
For instance, for a finite-dimensional parametric model $M=M_\theta,\,\theta\in\mathbb{R}^k,\, k\in\mathbb{N}$, with $m$ samples from the posterior $p(\theta\mid D,M)$, one can compute the Watanabe-Akaike information criterion (Gelman's definition) as $\text{WAIC} =-2\cdot(\hat{L}_W-2\hat{\Lambda}_W)$, with:
$$
\hat{L}_W = \sum_{i=1}^n\log\left(\frac{1}{m}\sum_{j=1}^m p(y_i\mid \theta_j) \right),\quad \text{ and }\quad \hat{\Lambda}_W=\hat{L}_W-\sum_{i=1}^n\frac{1}{m}\sum_{j=1}^m\log p(y_i\mid \theta_j) .
$$
The WAIC employs a posterior-sample-based approximation of the total ELPD. For this reason, it is considered a more fully Bayesian approach than the Akaike IC (AIC) or the Deviance IC (DIC). Model selection using the WAIC is asymptotically equivalent to model selection using LOO-CV. Furthermore, if the true model is also contained in the family of candidate models $\mathcal{M}$, both are also asymptotically equivalent to model selection using the model-posteriors:
$$
p(M\mid D) = \frac{p(M)\,\int_\theta p(\theta\mid M)\,\prod_{i=1}^n p(y_i\mid\theta,M)\, \text{d}\theta}{{\sum_{\tilde{M}\in\mathcal{M}} p(\tilde{M})\,\int_\theta p(\theta\mid \tilde{M})\,\prod_{i=1}^n p(y_i\mid\theta,\tilde{M})\, \text{d}\theta}},\,\quad M\in\mathcal{M}.
$$
Yet, when the true model is not in the candidate model space $\mathcal{M}$, these approaches might not agree. This is because, in such a case, the candidate model with the lowest KL divergence from the true generating distribution might not be the same model with the highest posterior within $\mathcal{M}$.
So, the ELDP also helps bridging Bayesian analysis with model selection based on Information Theory and Information Geometry, and it helps defining "the best" candidate model as the one in $\mathcal{M}$ with the lowest KL divergence from the true data-generating distribution. For some purposes, this optimal projection is just enough.
- The issues with prior specification and high-dimensional model spaces.
Computing $p(M\mid D)$ for all $M\in\mathcal{M}$ has many practical and computational drawbacks. To start, it requires the specification of model-priors $p(M)$. If $\mathcal{M}$ is composed of a finite number of candidates, then an uninformative model-prior can be employed: $p(M)=|\mathcal{M}|^{-1},\ \forall M\in\mathcal{M}$. Other model-prior specifications can be more difficult to motivate.
More importantly, the model-posteriors depend sensibly on how the specification of parameter-priors $p(\theta\mid M_k)$ change across models $M_k\in\mathcal{M}$. If models are highly heterogeneous, then the specification of $p(\theta\mid M_k)$ should reflect model-adaptive prior-uncertainty. But, $\theta$ is presented as model-shared parameters, and so different specifications may lead to different model-posteriors $p(M\mid D)$.
Finally, as $p(M\mid D)$ is probably not analytically available, a Monte Carlo approach would be required. This, in itself, presents several computational challenges. It is particularly tedious when model space $\mathcal{M}$ is discrete, large, and with a high number of posterior modes. For instance, take the case of $M$ representing a directed acyclic graph (DAG) on a set of $d\in\mathbb{N}$ variables in the data, with no structural constrains (such as prior time-ordering). For a large $d$, computing the MAP model is NP-hard, as the number of DAGs is super-exponential, $O\left(2^{d^2}\right)$, and thus exploring the entire model space for sampling would be practically unfeasible.
Clearly, some of these troubles would be avoided by employing ELPD-based model comparison.
I share some references on the topic:
Burnham, K. P., & Anderson, D. R. (Eds.). (2002). Model selection and multimodel inference: a practical information-theoretic approach. New York, NY: Springer New York.
Watanabe, S., & Opper, M. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of machine learning research, 11(12).
Gelman, A., Hwang, J., & Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. Statistics and computing, 24, 997-1016.
Yao, Y., Vehtari, A., Simpson, D., & Gelman, A. (2018). Using Stacking to Average Bayesian Predictive Distributions (with Discussion). Bayesian Analysis, 13(3), 917-1007.
