Comparing Misspecified Models with Big Data: A Variational Bayesian Perspective

Comparing Misspecified Models with Big Data: A Variational Bayesian Perspective

\def\spacingset#1{ {#1}} \spacingset{1}

\if11 \fi

\if01 \fi

{\it Keywords:} Kullback–Leibler divergence; Information criterion; Model misspecification; Variational Bayes; Massive Data.

\spacingset{1.9}

Introduction

In numerous empirical studies, parametric models are commonly employed. However, parametric models inherently carry the risk of model misspecification. As George Box famously stated, “All models are wrong, but some are useful.” When a model is misspecified, it can result in inefficient or, in some cases, inconsistent estimation of key parameters. Furthermore, likelihood-based statistical inferences, such as hypothesis testing and goodness-of-fit assessments, are significantly affected. Therefore, developing robust methods to address model misspecification is of critical importance.

Model comparison is one of the most critical issues in statistical inference. For a partial list of studies, see granger1995comments, phillips1994posterior, phillips1995bayesian,phillips1996econometric, hansen2005test, and burnham2008model. There are essentially two strands of literature on model selection VehtariandOjanen,anderson2004model. The first strand aims to answer the question which model best explains the observed data. The Bayes factor (BF, kass1995bayes) and its variations belong to this strand. They compare models by examining “posterior probabilities” given the observed data and search for the “true” model. Bayes Information Criterion (BIC, schwarz1978estimating) is a large sample approximation to BF, although it is based on the maximum likelihood estimator (MLE). The second strand comes from a predictive perspective, answering the question which model gives the best predictions of future observations, which are generated by the same mechanism that gives the observed data. From the predictive perspective, many penalty-based information criteria have been proposed for model comparison. In the frequentist framework, the two most popular information criteria are the Akaike Information Criterion (AIC) proposed by Akaike1973 and the Takeuchi Information Criterion (TIC) introduced by takeuchi1976distribution. Both are asymptotically unbiased estimators of the expected Kullback-Leibler (KL) divergence between the data generating process (DGP) and the plug-in predictive distribution when the MLE is used. The plug-in predictive distribution is obtained by substituting parameter values with their optimal estimates to produce the plug-in estimated sampling distribution. The AIC assumes that all candidate models either nest the true model or are good approximations of the DGP, whereas the TIC allows for model misspecification, with its penalty term involving the inverse of the Hessian matrix. Under the Bayesian framework, Deviance Information Criterion (DIC), proposed by spiegelhalter2002bayesian, is one of the most popular penalty-based predictive information criteria. In a recent study, li2020deviance developed a variant of DIC for comparing misspecified models, while li2024deviance proposed a decision-theoretic interpretation of DIC, demonstrating that DIC is the Bayesian version of AIC.

In recent years, several model selection approaches utilizing the Variational Bayes (VB) method have been introduced. A common strategy in VB-based model selection is to use the evidence lower bound (ELBO) as a proxy for the logarithm of the marginal likelihood function, \(\log p(\mathbf{y})\), to perform Bayes factor (BF) comparisons. corduneanu2001variational investigated VB model selection in the context of mixture models, and used the ELBO as a proxy to determine the optimal number of components. you2014variational explored the application of VB to classical Bayesian linear models. They established that, under mild regularity conditions, VB-based estimators possess desirable frequentist properties, such as consistency. Additionally, they proposed two VB-specific information criteria: the Variational AIC (VAIC), which substitutes the VB posterior mean into the DIC, and the Variational Bayesian Information Criterion (VBIC), which uses the ELBO as a proxy for the marginal likelihood. They further showed that VAIC is asymptotically equivalent to the frequentist AIC, while VBIC is first-order equivalent to the BIC in linear regression. Zhang_2024 proposed using the ELBO as an alternative criterion for model selection and demonstrated its asymptotic equivalence to the BIC. However, in the context of misspecified models and the era of massive data, there has been relatively little research on Bayesian model selection from a predictive perspective. This gap highlights the need for further investigation into model selection methodologies that prioritize predictive performance in such settings.

In this paper, we propose two new penalty-based predictive information criteria for model comparison in the context of misspecified models with massive data. First, based on the variational posterior distribution, we demonstrate that, from a predictive perspective, two types of predictive distributions can be derived: the variational plug-in predictive distribution and the variational posterior predictive distribution. Second, we examine the risk functions associated with these two variational predictive distributions, defined as the expectations of the KL divergence between the DGP and the predictive distributions. Third, under certain regularity conditions, we establish that the proposed information criteria are asymptotically unbiased estimators of their corresponding risk functions. Finally, through simulations and real-world case studies, we illustrate the application of the proposed information criteria.

The paper is organized as follows. Section (ref) briefly reviews the literature on how to make statistical inferences about misspecified models and VB technique for misspecified models with massive data. Section (ref) investigates the risk functions of variational predictive distributions. Section (ref) introduces the statistical decision theory and proposes the new penalized-based information criterion to compare misspecified models with massive data. Section (ref) illustrates the new methods using two simulated big data and two real big data. Section (ref) concludes the paper. The Appendix collects the proof of the theoretical results and VB analytical expression of parametric models used in the paper.

Statistical Inference for Misspecified Models: A Review

MLE-based Inference under Model Misspecification

Let the observed data be $\mathbf{y}=\left( y_1,\cdots ,y_{n}\right)$, with an i.i.d. data generating process (DGP) denoted by $g(\mathbf{y})$. Consider a parametric model, denoted by $p(\mathbf{y}|\mbox{\boldmath${ \theta}$}) $ used to fit the data, where $\mbox{\boldmath${ \theta}$}$ is a $P$-dimensional parameter, and $\mbox{\boldmath${ \theta}$}\in {\mbox{\boldmath${\Theta}$}}\subseteq R^{P}$. The Kullback-Leibler (KL) divergence is used to measure the “distance” between $g(\mathbf{y})$ and $p(\mathbf{y}|{\mbox{\boldmath${\theta}$}})$, that is,

where $E_{g(\mathbf{y})}$ is with respect to the DGP $g(\mathbf{y})$. Let ${\mbox{\boldmath${\theta}$}}^{*}$ $\in \Theta \subset R^{p}$ the pseudo true value that minimizes the KL divergence

and ${\mbox{\boldmath${\hat{\theta}}$}}$ denoted as the quasi maximum likelihood (QML) estimator of ${\mbox{\boldmath${\theta}$}}$, which maximizes the log-likelihood function of the parametric model,

For simplicity, let $l_t\left(\mathbf{y}_t,\boldsymbol{\theta}\right)=\ln p\left(\mathbf{y}_t | \boldsymbol{\theta}\right)$ represent the conditional log-likelihood for the $t^{t h}$ observation for any $1 \leq t \leq n$. We suppress $l_t\left(\mathbf{y}_t, \boldsymbol{\theta}\right)$ as $l_t(\boldsymbol{\theta})$, so that the log-likelihood function $\ln p(\mathbf{y} | \boldsymbol{\theta})$ is expressed as $\sum_{t=1}^n l_t(\boldsymbol{\theta})$. Define $\nabla^j l_t(\boldsymbol{\theta})$ as the $j^{t h}$ order derivative of $l_t(\boldsymbol{\theta})$ and $\nabla^j l_t(\boldsymbol{\theta})=l_t(\boldsymbol{\theta})$ when $j=0$. Let $\mathbf{\hat{J}}({\mbox{\boldmath${\theta}$}})=\frac{1}{n}\sum_{t=1}^{n}\nabla l_t(\boldsymbol{\theta})\nabla l_t(\boldsymbol{\theta})^{\prime }-\frac{1}{n}\sum_{t=1}^{n}\nabla l_t(\boldsymbol{\theta})\sum_{t=1}^{n}\nabla l_t(\boldsymbol{\theta})^{\prime}$, $\mathbf{\hat{I}}({\mbox{\boldmath${\theta}$}})=-\frac{1}{n} \sum_{t=1}^{n}\nabla^2 l_t(\boldsymbol{\theta})$. White (1982) established the maximum likelihood (ML) theory for misspecified models, that is,

as n goes to infinity where the asymptotic variance takes the sandwich form. If the model is correctly specified, then

as n goes to infinity.

Bayesian Inference under Model Misspecification

Consider a statistical model indexed by a set of $P$ parameters, $\boldsymbol{\theta} \in \boldsymbol{\Theta} \subseteq \mathbb{R}^P$, with a prior distribution $p(\boldsymbol{\theta})$ defined over $\boldsymbol{\theta}$. By applying Bayes' theorem, the posterior distribution can be expressed as:

where $p(\mathbf{y}) = \int p(\mathbf{y}|\boldsymbol{\theta})p(\boldsymbol{\theta})d\boldsymbol{\theta}$ represents the marginal likelihood.

In most cases, the posterior distribution $p(\boldsymbol{\theta}|\mathbf{y})$ does not have a closed-form solution. Consequently, posterior sampling is typically conducted using Markov Chain Monte Carlo (MCMC) techniques gelman2003bayesian. Based on the random samples generated from posterior simulations, Bayesian statistical inference can be performed using the corresponding sample means and covariance matrices. For example, let $\{\boldsymbol{\theta}^{(j)} : j = 1, 2, \cdots, J\}$ denote the effective random samples generated from the posterior distribution after discarding burn-in samples. Bayesian estimates of $\boldsymbol{\theta}$ and the associated standard error can then be calculated as: $ \boldsymbol{\bar\theta}=\frac {1}{J}\sum_{j=1}^J\boldsymbol{\theta}^{(j)}, \widehat{Var(\boldsymbol{\theta}|\mathbf{y})} =\frac{1}{J-1}\sum_{j=1}^J(\boldsymbol{\theta}^{(j)}-\boldsymbol{\bar{\theta}}) (\boldsymbol{\theta}^{(j)}-\boldsymbol{\bar{\theta}})^{\prime}. $

These Bayesian estimates are consistent estimators of the posterior mean and covariance matrix. It is well documented in the literature that MCMC techniques are powerful and efficient for posterior simulation. Due to advances in MCMC, Bayesian methods have gained significant popularity for statistical inference and are now widely applied to a variety of complex models.

It is worth noting that the Bayesian large-sample theory exhibits a key difference from the QML large-sample theory, particularly for misspecified models. Unlike QML theory, Bayesian asymptotic results do not differ between correctly specified and misspecified models. In both cases, the Bayesian large-sample theory is given by:

in probability as $n \to \infty$ kleijn2012bernstein.

Variational Bayes for Misspecified Models with Massive Data

To compute \( p(\boldsymbol{\theta}|\mathbf{y}) \), the dominant paradigm in Bayesian statistics is MCMC, including the Metropolis-Hastings algorithm metropolis1953equation, hastings1970monte and the Gibbs sampler geman1984stochastic, among others. While MCMC provides a flexible and widely applicable method to sample from the posterior distribution of \(\boldsymbol{\theta}\), it faces significant challenges, particularly when applied to massive datasets where the sample size \(n\) is extremely large.

One notable scenario in which the log-likelihood becomes computationally intractable is when dealing with massive data bardenet2017markov, quiroz2019speeding. In such cases, the log-likelihood function is represented by the summation of numerous terms, making it prohibitively expensive to evaluate. Due to the high computational cost associated with likelihood evaluations for massive datasets, MCMC methods can require hours or even days to converge to a stationary posterior distribution.

Recently, to address the limitations of Bayesian inference based on MCMC for massive datasets, Variational Bayes (VB) methods jordan1999introduction, have garnered significant attention in the research community. VB offers an alternative to MCMC by solving the following optimization problem:

where \(p^{VB}(\boldsymbol{\theta}|\mathbf{y})\) denotes the Variational Bayesian posterior, and the goal is to approximate the posterior \(p(\boldsymbol{\theta}|\mathbf{y})\) using a tractable variational family \(\Gamma\). A commonly used variational family is the mean-field (MF) family, which assumes the factorized form: $q(\boldsymbol{\theta}) = \prod_{i=1}^P q_{\theta_i}(\theta_i).$ This simplification facilitates efficient optimization by reducing computational complexity.

Since VB formulates posterior inference as an optimization problem, it provides a computationally efficient alternative to MCMC, particularly in the context of massive datasets under Bayesian modeling attias2013inferring, bishop2006pattern. Empirical studies have shown that VB-based algorithms can be orders of magnitude faster than MCMC blei2017variational, gunawan2017fast. Beyond the classical mean-field VB, advances such as stochastic variational inference (SVI) hoffman2013stochastic have further enabled scalable Bayesian analysis for large-scale datasets.

The asymptotic properties of the VB posterior have been a topic of significant interest in the literature. Define the second-order derivative of the log-likelihood as $\mathbf{\bar{H}}_n(\boldsymbol{\theta}) := \frac{1}{n} \frac{\partial^2 \ln p(\mathbf{y}, \boldsymbol{\theta})}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^\prime}$, and take the expectation to obtain $\mathbf{H}_n(\boldsymbol{\theta}) := E[\mathbf{\bar{H}}_n(\boldsymbol{\theta})] $, then the normal approximation to the VB posterior can be expressed as:

where \(\mathbf{H}_n^d\) is a diagonal matrix whose diagonal elements match those of \(\mathbf{H}_n\). As established by han2019statistical and Zhang_2024, the KL divergence between the VB posterior \(p^{VB}(\boldsymbol{\theta}|\mathbf{y})\) and the normal approximation \(p^{VBN}(\boldsymbol{\theta}|\mathbf{y})\) converges to \(0\) in probability as \(n \to \infty\). wang2019variational proved that the total variation between the VB posterior and $p^{VBN}(\boldsymbol{\theta}|\mathbf{y})$ converges to \(0\) in probability as \(n \to \infty\).

Risk of Predictive Distributions on Misspecified Models based on Variational Bayes

In the literature, assessing the utility of a misspecified statistical model is typically achieved by examining its predictive performance bernardo1979expected. Given a set of future observations \(\mathbf{y}_f\), the predictive distribution is denoted by \(p_f(\mathbf{y}_f|\mathbf{y})\). A commonly used approach for quantifying the predictive performance of a misspecified model is to compute the KL divergence between the true data-generating process \(g(\mathbf{y}_f)\) and the predictive distribution \(p_f(\mathbf{y}_f|\mathbf{y})\), scaled by a factor of 2. This measure is expressed as: \[ 2 \times KL\left[ g\left( \mathbf{y}_f \right), p_f\left( \mathbf{y}_f|\mathbf{y} \right) \right] = 2 E_{\mathbf{y}_f}\left[ \ln \frac{g\left( \mathbf{y}_f \right)}{p\left( \mathbf{y}_f|\mathbf{y} \right)} \right], \] which can be rewritten as $2 \int \left[ \ln \frac{g\left( \mathbf{y}_f \right)}{p\left( \mathbf{y}_f|\mathbf{y} \right)} \right] g\left( \mathbf{y}_f \right) \, d\mathbf{y}_f.$ Building on this KL divergence, statistical decision theory allows the specification of a loss function associated with a decision \(d\) as: \[ \mathcal{L}(\mathbf{y}, d) = 2 \times KL\left[ g\left( \mathbf{y}_f \right), p\left( \mathbf{y}_f|\mathbf{y}, d \right) \right], \] where \(p(\mathbf{y}_f|\mathbf{y}, d)\) represents the predictive density based on decision \(d\). The corresponding risk function is then defined as good1952rational: \[ \text{Risk}(d) = E_{\mathbf{y}}\left[ \mathcal{L}(\mathbf{y}, d) \right] = \int \mathcal{L}(\mathbf{y}, d) g(\mathbf{y}) \, d\mathbf{y}. \] In the context of VB, two types of predictive distributions can be derived for prediction: the variational plug-in predictive distribution and the variational posterior predictive distribution. These two distributions correspond to different statistical decisions, resulting in two distinct risk functions. In the subsequent subsection, we evaluate these two risk functions and derive estimators for them. To facilitate this analysis, we first establish the necessary notations and outline mild regularity conditions.

Let $ \mathbf{y}:=(y_{1},\ldots ,y_{n})$ and $l_{t}\left( \mathbf{y}_{t},{ \mbox{\boldmath${\theta}$}}\right) =\ln p(\mathbf{y}_{t}|{ \mbox{\boldmath${\theta}$}})$ be the conditional log-likelihood for the $ t^{th}$ observation for any $1\leq t\leq n$. For simplicity, we suppress $ l_{t}\left( \mathbf{y}_{t},{\mbox{\boldmath${\theta}$}}\right) $\ as $ l_{t}\left( {\mbox{\boldmath${\theta}$}}\right) $ so that the log-likelihood function $\ln p(\mathbf{y|}{\mbox{\boldmath${\theta}$})}$ is $ \sum_{t=1}^{n}l_{t}\left( {\mbox{\boldmath${\theta}$}}\right) $. And define $\bigtriangledown ^{j}l_{t}\left( {\mbox{\boldmath${\theta}$}} \right) $\ to be the $j^{th}$\ derivative of $l_{t}\left( { \mbox{\boldmath${\theta}$}}\right) $ and $\bigtriangledown ^{j}l_{t}\left( { \mbox{\boldmath${\theta}$}}\right) =l_{t}\left( {\mbox{\boldmath${\theta}$}} \right) $\ when $j=0$. We suppress the superscript when $j=1$, and

Then, the following regularity conditions can be imposed

Assumption 1: ${\mbox{\boldmath${\Theta}$}}\subset R^{P}$ is compact.

Assumption 2: The data $\mathbf{y}=(y_{1},\ldots ,y_{n})$ is independent and identically distributed.

Assumption 3: For all $t$, $l_{t}\left( {\mbox{\boldmath \ensuremath{{\theta}}}}\right) $\ is eight-times differentiable on $ \mbox{\boldmath\ensuremath{{\Theta}}}$ almost surely.

Assumption 4: For $j=0,1,2,3$, for any ${\mbox{\boldmath${\theta}$}, \mbox{\boldmath${\theta}$}^{\prime }\in \mbox{\boldmath${\Theta}$}}$, $ \left\Vert \bigtriangledown ^{j}l_{t}\left( {\mbox{\boldmath${\theta}$}} \right) -\bigtriangledown ^{j}l_{t}\left( {\mbox{\boldmath${\theta}$}} ^{\prime }\right) \right\Vert \leq c_{t}^{j}\left( \mathbf{y}_{t}\right) \left\Vert {\mbox{\boldmath${\theta}$}}-{\mbox{\boldmath${\theta}$}}^{\prime }\right\Vert $ in probability, where $c_{t}^{j}\left( \mathbf{y}_{t}\right) $ is a positive random variable with $\sup_{t}E\left\Vert c_{t}^{j}\left( \mathbf{y}_{t}\right) \right\Vert <\infty $ and $\frac{1}{n} \sum_{t=1}^{n}\left( c_{t}^{j}\left( \mathbf{y}_{t}\right) -E\left( c_{t}^{j}\left( \mathbf{y}_{t}\right) \right) \right) \overset{p}{ \rightarrow }0$.

Assumption 5: For $j=0,1,\ldots ,4$, there exists a function $M_{t}( \mathbf{y}_{t})$ such that for all ${\mbox{\boldmath${\theta}$}}\in {\ \mbox{\boldmath${\Theta}$}}$, $\bigtriangledown ^{j}l_{t}\left( {\ \mbox{\boldmath${\theta}$}}\right) $ exists, $\sup_{{\mbox{\boldmath${ \theta}$}\in \mbox{\boldmath${\Theta}$}}}\left\Vert \bigtriangledown ^{j}l_{t}\left( {\mbox{\boldmath${\theta}$}}\right) \right\Vert \leqslant M_{t}(\mathbf{y}_{t})$, and $\sup_{t}E\left\Vert M_{t}(\mathbf{y} _{t})\right\Vert ^{r+\delta }\leq M<\infty $ for some $\delta >0$ and $r>2$.

Assumption 6: Let $\mbox{\boldmath${\theta}$}_{n}^{p}$ be the pseudo-true value that minimizes the KL loss between the DGP and the candidate model

where $\left\{ \mbox{\boldmath${\theta}$}_{n}^{p}\right\} $ is the sequence of minimizers interior to ${\mbox{\boldmath${\Theta}$}}$ uniformly in $n$. For all $\varepsilon >0$,

where $N\left( \mbox{\boldmath${\theta}$}_{n}^{p},\varepsilon \right) $ is the open ball of radius $\varepsilon $ around $\mbox{\boldmath${\theta}$} _{n}^{p}$.

Assumption 7: The sequence $\left\{ \mathbf{H}_{n}\left( \mbox{\boldmath${\theta}$}_{n}^{p}\right) \right\} $ is negative definite and the sequence $\left\{ \mathbf{B}_{n}\left( \mbox{\boldmath${\theta}$} _{n}^{p}\right) \right\} $ is positive definite, both uniformly in $n$.

Assumption 8: The prior density $p\left( \mbox{\boldmath \ensuremath{{\theta}}}\right) $ is thrice continuously differentiable and $0<p\left( {\mbox{\boldmath\ensuremath{{\theta}}}}_{n}^{0}\right) <\infty $\ uniformly in $n$. Moreover, there exists an $n^{\ast }$ such that, for any $n>n^{\ast }$, the posterior distribution $p\left( {\ \mbox{\boldmath\ensuremath{{\theta}}}}|\mathbf{y}\right) $ is proper and $ \int \left\Vert {\mbox{\boldmath\ensuremath{{\theta}}}}\right\Vert ^{2}p\left( {\mbox{\boldmath\ensuremath{{\theta}}}}|\mathbf{y}\right) d{\ \mbox{\boldmath\ensuremath{{\theta}}}}<\infty $.

Assumptions 1-7 are well-known primitive conditions for developing the QML theory, namely consistency and asymptotic normality, for independent and identically distributed data; see, for example, gallant1988unified and WOOLDRIDGE19942639. Assumption 8 is the regular condition for prior density, see, for example, li2020deviance. Assumptions 1-8 are sufficient for the assumptions used by Zhang_2024 to develope the asymptotic properties of VB posterior distribution without latent variables.

Risk of VB Plug-in Predictive Distribution

Under VB inference, for a potentially misspecified model, let $\overline{\boldsymbol{\theta}}^{VB}$ denote the VB estimator of the parameter $\boldsymbol{\theta}$ which corresponds to the posterior mean of the variational posterior distribution $p^{VB}(\boldsymbol{\theta}|\mathbf{y})$. In cases where the posterior mean does not have a closed-form analytical solution, it can generally be approximated consistently using the sample mean $\boldsymbol{\bar\theta}^{VB}=\frac {1}{J}\sum_{j=1}^J\boldsymbol{\theta_{VB}}^{(j)} $, where $\boldsymbol{\theta_{VB}}^{(j)}$, $j=1,2,\cdots J$ are generated from $p^{VB}(\boldsymbol{\theta}|\mathbf{y})$.

Building on the literature regarding the development of popular information criteria such as AIC, TIC, and DIC, we assume the existence of future replicated data \(\mathbf{y}_{rep}\), which shares the same DGP as the observed data \(\mathbf{y}\) and independent of \(\mathbf{y}\). For more details on the concept of \(\mathbf{y}_{rep}\), one may refer to the comprehensive discussion in the seminal textbook on model selection by anderson2004model and the references therein. For the future data \(\mathbf{y}_{rep}\), the VB plug-in predictive distribution can be expressed as $p\left(\mathbf{y}_{rep}| \overline{\boldsymbol{\theta}}^{VB}\right)$, where \(\overline{\boldsymbol{\theta}}^{VB}\) represents the VB estimator, typically the posterior mean of the variational posterior distribution. The predictive distribution provides a probabilistic framework for evaluating future observations based on the fitted model. Correspondingly, the loss function associated with the statistical decision, denoted as \(d_1\), can be specified as follows:

In this context, the risk function can be expressed as:

Since $E_{\mathbf{y}}E_{\mathbf{y}_{rep}} \left[ 2\ln g\left( \mathbf{y}_{rep}\right) \right] $ is the same across all statistical decisions, the risk function can be expressed as:

where $C=E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left[ 2\ln g\left( \mathbf{y}_{rep}\right) \right]$.

It is evident that a smaller value of \( \text{Risk}(d_1) \) indicates better performance of the predictive distribution \( p\left(\mathbf{y}_{rep}| \overline{\boldsymbol{\theta}}^{VB} \right) \) in predicting the replicate data \(\mathbf{y}_{rep}\). However, in general, this risk function does not have a closed-form analytical expression. Therefore, evaluating the risk function is essential for assessing the predictive behavior of the model.

To address this challenge, we derive an asymptotic expansion of the risk function, as presented in the following theorem. This derivation provides a practical approach to approximate the risk function in large-sample scenarios, offering insights into the predictive performance of the VB-based approach.

Risk of VB Posterior Predictive Distribution

Under the Bayesian framework, the VB posterior predictive distribution for the replicated data \(\mathbf{y}_{rep}\), corresponding to \(p^{VB}(\boldsymbol{\theta}|\mathbf{y})\), is defined as:

As described in Section 3.1, the KL divergence between the true data-generating process \(g\left( \mathbf{y}_{rep} \right)\) and the VB posterior predictive distribution \(p^{VB}\left( \mathbf{y}_{rep}|\mathbf{y}\right)\), multiplied by 2, is given by:

This divergence is used to quantify the predictive performance of the VB posterior predictive distribution. Accordingly, the loss function associated with the statistical decision \(d_2\), which involves using the VB posterior predictive distribution for prediction, is defined as: \[ \mathcal{L}(\mathbf{y}, d_2) = 2 \times KL\left[ g\left( \mathbf{y}_{rep} \right), p^{VB}\left( \mathbf{y}_{rep}|\mathbf{y} \right) \right]. \] The corresponding risk function for the decision \(d_2\) can be expressed as:

which can be further rewritten as: \[ \text{Risk}(d_2) = C + E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left[ -2 \ln p^{VB}\left( \mathbf{y}_{rep}|\mathbf{y} \right) \right], \] where \(C = E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left[ 2 \ln g\left( \mathbf{y}_{rep} \right) \right]\) is a constant that depends only on the DGP.

From this expression, it is evident that a smaller \(\text{Risk}(d_2)\) indicates better predictive performance of \(p^{VB}\left( \mathbf{y}_{rep}|\mathbf{y} \right)\) in approximating \(g\left( \mathbf{y}_{rep} \right)\). In the following, we derive an asymptotic expansion of this risk function via the following theorem.

Predictive Information Criteria for Comparing Misspecified Models with Massive Data based on VB

In this section, we outline the development of new predictive information criteria for model comparison in the context of misspecified models with massive data. Building on the risk functions analyzed in Section (ref), Section (ref) introduces the framework of statistical decision theory for model comparison. In Section (ref), we propose an information criterion, termed VDIC$_{M}$, based on the VB plug-in predictive distribution. We then present another information criterion, termed VPIC, which is constructed using the VB posterior predictive distribution in Section (ref). At last, in Section (ref), we then discuss BFs and BIC in the context of misspecified models.

Statistical Decision Theory based on Risk Function for Model Selection

In this section, from a predictive perspective, we extend the decisional framework introduced in Section (ref) to develop information criteria for model comparison. Suppose there are \(K\) candidate models, all of which may be misspecified, and the task is to select the most suitable model. These candidate models are denoted by \(M_{k}\), where \(k = 1, 2, \ldots, K\). As discussed in the previous section, this selection is achieved by minimizing the risk associated with the statistical decision.

Assume that the probabilistic behavior of the observed data \(\mathbf{y} \in \mathbf{Y}\) is described by a set of probabilistic models \(\{M_k\}_{k=1}^K := \{p(\mathbf{y}|\boldsymbol{\theta}_k, M_k)\}_{k=1}^K\), where \(\boldsymbol{\theta}_k\) represents the set of parameters associated with model \(M_k\). Formally, the model selection problem can be framed as a decision-making problem, where the goal is to select one model from \(\{M_k\}_{k=1}^K\). In this context, the action space comprises \(K\) elements, denoted by \(\{d_k\}_{k=1}^K\), where \(d_k\) indicates that model \(M_k\) is selected.

For the decision-making process, as in Section (ref), a loss function \(\mathcal{L}(\mathbf{y}, d_k)\) must be specified. This loss function quantifies the loss incurred by selecting decision \(d_k\). Given the loss function, the corresponding risk can be defined as: \[ \text{Risk}(d_k) = E_{\mathbf{y}}\left[ \mathcal{L}(\mathbf{y}, d_k) \right] = \int \mathcal{L}(\mathbf{y}, d_k) g(\mathbf{y}) d\mathbf{y}, \] where \(g(\mathbf{y})\) is the DGP. Consequently, the model selection problem is equivalent to optimizing the statistical decision by minimizing the risk: \[ k^\ast = \arg \min_k \text{Risk}(d_k). \] Based on the set of candidate models \(\{M_k\}_{k=1}^K\), the model \(M_{k^\ast}\), corresponding to the decision \(d_{k^\ast}\), is selected as the optimal model.

The quantity used to assess the predictive ability of a candidate model is the KL divergence between the DGP \(g\left( \mathbf{y}_{rep} \right)\) and a predictive distribution \(p\left( \mathbf{y}_{rep}|\mathbf{y},M_k \right)\), scaled by a factor of 2: \[ 2 \times KL\left[ g\left( \mathbf{y}_{rep} \right), p\left( \mathbf{y}_{rep}|\mathbf{y},M_k \right) \right] = 2 E_{\mathbf{y}_{rep}}\left[ \ln \frac{g\left( \mathbf{y}_{rep} \right)}{p\left( \mathbf{y}_{rep}|\mathbf{y},M_k \right)} \right], \] which can also be written as $ 2 \int \left[ \ln \frac{g\left( \mathbf{y}_{rep} \right)}{p\left( \mathbf{y}_{rep}|\mathbf{y},M_k \right)} \right] g\left( \mathbf{y}_{rep} \right) d\mathbf{y}_{rep}.$ Similar to the framework introduced in Section 3, the loss function associated with the decision \(d_k\) is defined as $ \mathcal{L}(\mathbf{y}, d_k) = 2 \times KL\left[ g\left( \mathbf{y}_{rep} \right), p\left( \mathbf{y}_{rep}|\mathbf{y},M_k \right) \right].$ Thus, the model selection problem is formulated as: \[ k^\ast = \arg \min_k \text{Risk}(d_k) = \arg \min_k E_{\mathbf{y}}\left[ \mathcal{L}(\mathbf{y}, d_k) \right], \] which can be further expanded as: \[ k^\ast = \arg \min_k \left\{ 2 \times E_{\mathbf{y}} E_{\mathbf{y}_{rep}}\left[ \ln \frac{g\left( \mathbf{y}_{rep} \right)}{p\left( \mathbf{y}_{rep}|\mathbf{y},M_k \right)} \right] \right\}. \] Rearranging terms gives: \[ k^\ast = \arg \min_k \left\{ E_{\mathbf{y}} E_{\mathbf{y}_{rep}}\left[ 2 \ln g\left( \mathbf{y}_{rep} \right) \right] + E_{\mathbf{y}} E_{\mathbf{y}_{rep}}\left[ -2 \ln p\left( \mathbf{y}_{rep}|\mathbf{y},M_k \right) \right] \right\}. \] Since \(g\left( \mathbf{y}_{rep} \right)\) is the DGP, the term \(E_{\mathbf{y}_{rep}}\left[ 2 \ln g\left( \mathbf{y}_{rep} \right) \right]\) is constant across all candidate models and can therefore be omitted from the equation. Consequently, the model selection problem simplifies to: \[ k^\ast = \arg \min_k \text{Risk}(d_k) = \arg \min_k E_{\mathbf{y}} E_{\mathbf{y}_{rep}}\left[ -2 \ln p\left( \mathbf{y}_{rep}|\mathbf{y},M_k \right) \right]. \] The smaller the value of \(\text{Risk}(d_k)\), the better the performance of the candidate model in using the predictive distribution \(p\left( \mathbf{y}_{rep}|\mathbf{y},M_k \right)\) to approximate \(g\left( \mathbf{y}_{rep} \right)\). Evaluating the risk among candidate models is therefore essential for making the optimal decision.

It is important to note that the action space in this context is larger than in previous cases. From a predictive perspective, we not only need to select a model for prediction but also determine which predictive distribution to use. The action space is denoted byby $\left\{d_{k^1}, d_{k^2}\right\}_{k=1}^K$ where $d_{k^a}(a \in(1,2)) $ means $M_k$ is selected, and the predictions are generated from $p\left( \mathbf{y}_{rep}|\mathbf{y}, M_k, d_a\right)$. If $a=1$, it means that the VB plug-in predictive distribution, $p\left(\mathbf{y}_{\text {rep}}|\mathbf{y}, M_k, d_1\right)=p\left(\mathbf{y}_{\text {rep}}|\overline{\boldsymbol{\theta}}^{VB}, M_k\right)$ is used; if $a=2$, it means that the VB posterior predictive distribution, $p\left(\mathbf{y}_{\text {rep }}|\mathbf{y}, M_k, d_2\right)=$ $p^{VB}\left(\mathbf{y}_{rep} \mid \mathbf{y}, M_k\right)$ is used. The KL divergence for this setup is defined as

where $p\left(\mathbf{y}_{r e p}|\mathbf{y}, d_{k^a}\right):=p\left(\mathbf{y }_{r e p} \mid \mathbf{y}, M_k, d_a\right)$. The risk associated with $ d_{k^a}$ is then given by

Consequently, the model selection problem is equivalent to solving the following statistical decision problem:

Since the DGPs \(g\left(\boldsymbol{y}\right)\) and \(g\left(\boldsymbol{y}_{rep}\right)\) are unknown, directly evaluating the risk associated with decision \(d_{k^a}\) is infeasible. However, it is possible to approximate the risk by using an asymptotically unbiased estimator of \(\text{Risk}\left(d_{k^a}\right)\). As noted in the literature, various information criteria proposed for model selection can be interpreted as asymptotically unbiased estimators of the expected loss function, up to a constant, under different statistical decision frameworks vrieze2012model.

Traditionally, model selection has been conducted using information criteria that assess the relative quality of statistical models for a given dataset. Under the frequentist framework, criteria such as AIC, TIC, and their variants have been widely applied. Under the Bayesian framework, criteria include DIC and its extensions, such as the deviance information criterion for misspecified models (\(\text{DIC}_M\)) proposed by li2020deviance. These information criteria have been shown to follow the principles of statistical decision theory discussed above. Specifically, AIC, TIC, DIC, and \(\text{DIC}_M\) are all constructed by estimating the KL divergence between the DGP and the corresponding predictive distributions. In this study, we develop new approaches that adhere to a similar decision-theoretical framework. To provide context, we first present two remarks that introduce these popular information criteria within this framework. Subsequently, we propose our new information criteria in the following subsections.

More details about the theoretical development of $\text{DIC}^{k}$ and $\text{DIC}_{M}^{k}$, one can refer to takeuchi1976distribution , li2020deviance, li2024deviance and reference therein.

Information Criterion for Comparing Misspecified Models based on Variational Bayes Plug-in Predictive Distributions

Following the statistical decision theory shown in section (ref), we utilize $\ln p(\mathbf{y}|{\mbox{\boldmath${\bar{\theta}}$}} _{k}^{VB},M_{k})$ to construct the loss function and the corresponding risk function. Subsequently, similar to existing information criteria such as AIC, TIC, DIC and $\text{DIC}_M$,we propose a new information criterion for model selection. Let $ \mathbf{\bar{\Omega}}_{n}\left( {\mbox{\boldmath${ \widehat{\theta}}$}} _{k}\right) $ , $\mathbf{\bar{H}}_{n}\left({\ \mbox{\boldmath${\widehat{ \theta}}$}}_{k}\right)$ be consistent estimators of $\mathbf{B}_{n}({ \mbox{\boldmath${\theta}$}}_{n}^{p})$ and $\mathbf{H}_{n}({ \mbox{\boldmath${\theta}$}}_{n}^{p})$ respectively. Based on the results of han2019statistical and Zhang_2024, we have

where ${\mbox{\boldmath${\bar{\theta}}$}}_{k}^{VB}$ is the mean of variational posterior density $p^{VB}\left( \mbox{\boldmath${\theta}$}| \mathbf{y}\right) $. Using this, we derive the consistent estimators of $\mathbf{B}_{n}({\mbox{\boldmath${\theta}$}}_{n}^{p})$ and $\mathbf{H}_{n}({ \mbox{\boldmath${\theta}$}}_{n}^{p})$ as $\mathbf{\bar{\Omega}}_{n}\left( { \mbox{\boldmath${\bar{\theta}}$}}_{k}^{VB}\right) $ and $\mathbf{\bar{H}} _{n}\left( {\ \mbox{\boldmath${\bar{\theta}}$}}_{k}^{VB}\right) $, respectively. To account for model misspecification, we define a new information criterion, termed the Variational Deviance Information Criterion under Model Misspecification ($\text{VDIC}_M$), using the variational plug-in predictive density:

where the penalty term $ P^k_{VDIC_M}$ for model $k$ is defined as

It can be proved that $\text{VDIC}_M^{k}$ is an asymptotically unbaised estimator of $Risk(d_{k^1})$ up to a constant.

Information Criterion for Comparing Misspecified Models based on the VB Posterior Predictive Distribution

Following the statistical decision theory outlined in Section (ref), we utilize \(p^{VB}(\mathbf{y}_{rep}|\mathbf{y})\) to construct the loss function and the corresponding risk function. Based on this posterior predictive distribution, a new information criterion can be developed to estimate \(\text{Risk}(d_{k^2})\).

Let \(\mathbf{\bar{\Omega}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB})\) and \(\mathbf{\bar{H}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB})\) be consistent estimators of \(\mathbf{B}_{n}({\boldsymbol{\theta}}_n^{p})\) and \(\mathbf{H}_{n}({\boldsymbol{\theta}}_n^{p})\), respectively. The consistent estimator of \(\mathbf{C}_{n}\) is given by: \[ \mathbf{\hat{C}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB}) = \left(\mathbf{\bar{H}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB})\right)^{-1} \mathbf{\bar{\Omega}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB}) \left(\mathbf{\bar{H}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB})\right)^{-1}, \] where \({\boldsymbol{\bar{\theta}}}_k^{VB}\) represents the mean of the variational posterior density \(p^{VB}(\boldsymbol{\theta}|\mathbf{y})\).

When accounting for model misspecification, we define a new information criterion based on the VB posterior predictive density, termed the Variational Predictive Information Criterion (\(\text{VPIC}\)):

where the penalty term \(P_{VPIC}^k\) for model \(k\) is defined as: \[

\]

It can be proved that $\text{VPIC}^{k}$ is an asymptotically unbaised estimator of $Risk(d_{k^2})$ up to a constant.

BF and BIC type information criteria

The BF and BIC belong to the first strand of model comparison in the section (ref). They compare competing models by examining model posterior probabilities and search for the“true" model. Both BFs and BIC enjoy the property of consistency, that is, when the true DGP is one of the candidate models, BFs and BIC select it with probability approaching 1 when the sample size goes to infinity.

Suppose there are two candidate models, $M_1$ and $M_2$. The BF of $M_1$ against $M_2$ is defined as $B_{12}=\frac{p\left(\mathbf{y} | M_1\right)}{p\left(\mathbf{y} | M_2\right)}$, where $p\left(y | M_k\right)$ is the marginal likelihood of model $M_k$ which is obtained by $$ p\left(\mathbf{y} | M_k\right)=\int_{\Theta_k} p\left(\mathbf{y} | \boldsymbol{\theta}_k, M_k\right) p\left(\boldsymbol{\theta}_k | M_k\right) \mathrm{d} \boldsymbol{\theta}_k, \quad \boldsymbol{\theta_k} \in \Theta_k, k=1,2, $$ where $\boldsymbol{\theta}_k$ is the set of parameters in $M_k, p\left(\mathbf{y} | \boldsymbol{\theta}_k, M_k\right)$ the likelihood function of $M_k, p\left(\boldsymbol{\theta}_k \mid M_k\right)$ the prior of $\boldsymbol{\theta}_{\boldsymbol{k}}$ in $M_k$. If $B_{12}>1$, $M_1$ is preferred to $M_2$ and vice versa.

Based on the Laplace approximation, schwarz1978estimating showed that the log-marginal likelihood can be approximated by

where $\hat{\boldsymbol{\theta}}_k$ is the MLE of $\boldsymbol{\theta}_k$ and $\overline{\mathbf{H}}_n\left(\hat{\boldsymbol{\theta}}_k\right)$, and $P_k$ is the dimension of $\boldsymbol{\theta}_k$. Ignoring all the $O_p(1)$ terms in ((ref)) and under noninformative priors such as $p\left(\boldsymbol{\theta}_k | \boldsymbol{M}_k\right) \propto 1$, Schwarz defined $\mathrm{BIC}_k$ as $\mathrm{BIC}_k =-2 \ln p\left(\mathbf{y}|\hat{\boldsymbol{\theta}}_k, M_k\right)+P_k \ln n$, where, as in AIC and TIC, $-2 \ln p\left(\mathbf{y} | \hat{\boldsymbol{\theta}}_k, M_k\right)$ is used to measure the model fit, but $P_k \ln n$ is the new penalty term. Obviously, $\mathrm{BIC}_k$ provides an approximation of $-2 \ln \left(\mathbf{y} | M_k\right)$.

Recently, Zhang_2024 showed that under regular conditions, the difference between the evidence lower bound, which is the by-product of VB algorithm, and $\mathrm{-BIC/2}$, is asymptotically constant as $n$ goes to infinity.

Simulation and Empirical Studies

Simulation Study

We begin by using two numerical simulation examples to evaluate the performance of our newly proposed criteria in the context of massive data. Both examples involve model misspecification. In the first study, we use polymomial regression to fit a nonlinear model, aming to select the model with the best predicitve among candidate models. Similarly, in the second study, we focus on identifying the "best" model among four candidate probit models. For each scenario, we conduct 1000 replications and apply our two newly developed information criteria to assess their effectiveness with the ELBO criterion proposed by Zhang_2024.

In every experiment, we simulate $\mathbf{y}_i$ and calculate $\text{VPIC}^k$, $\text{VDIC}_M^k$, $\text{ELBO}^k$, $\text{AIC}^k$ and $\text{BIC}^k$ of candidate model $M_k, k=1,\dots,K$. Each of the five criteria is used to selected a best model (call it $M_{k_i^*}$), we then record this model and the corresponding optimal criteria $IC(i)$. For $\text{VDIC}_M^k$, we use the VB plug-in predictive distribution $p^{VB}(\mathbf{y}_{rep}|\boldsymbol{\bar{\theta}}^{VB}_{k^{*}},M_{k^{*}})$ under the best model $M_{k*}$ to predict new data. Then we can estimate the risk by

where $Risk\left(d_{k_*^1}\right)=E_{\mathbf{y}}\left[\mathcal{L}\left(\mathbf{y}, d_{k^*}\right)\right]=E_{\mathbf{y}}\left[2\times KL \left[g(\mathbf{y}_{rep}),p(\mathbf{y}_{rep}|\mathbf{y},M_{k^*},d_1)\right]\right]$

For $\text{VPIC}^k$, we use the VB posterior predictive distribution $p^{VB}(\mathbf{y}_{rep}|\mathbf{y},M_{k*})$ under the best model $M_{k*}$ to predict new data. Then we can estimate the risk by

Same risk is calculate to estimate the risk of $\text{AIC}^k$.

For $\text{ELBO}^k$ and $\text{BIC}^k$, we will use two proxies to evaluate its risk. As Zhang_2024 noted, under some regular conditions, the difference between $-\text{BIC}^k/2$ and $\text{ELBO}^k$ is asymptotically to be constant as $n$ goes to infinity. For the reasons that $\text{BIC}$ is constructed as an approximation of the marginal likelihood $p(\mathbf{y})$, not from predictive perspective, averaging for $-2\times\text{ELBOs}$ and $\text{BIC}$ in all replication as the risk of both $\text{ELBO}$ and $\text{BIC}$ is not a proper way. We will use two proxies to see the relative risk of $\text{ELBO}$. In each experiment, when choosing the best model $M_{k_i^*}$ under $\text{ELBO}$ or $\text{BIC}$, we will use both $\text{VDIC}_M^{k_i^*}$ and $\text{VPIC}^{k_i^*}$ whose expectation is the KL loss as proxy. Then we can estimate the risk of $\text{ELBO}$ and $\text{BIC}$ by

named as $\text{ELBO1}$, $\text{ELBO2}$, $\text{BIC1}$, and $\text{BIC2}$.

Polymomial Regression

We begin with a simple experiment to compare alternative model selection criteria when the true DGP is not included in the set of candidate models. In other words, all candidate models are misspecified. Following Ding et al. (2019), we generate data from the following model $$ y_i=\ln \left(1+46 x_i\right)+e_i, e_i \sim N(0,1), i=1, \ldots, N, $$ where $x_i=0.7(i-1) / n$ which is fixed under repeated sampling by design. In practice, researchers do not know the functional form. Suppose the following set of polynomial regressions is considered, $$ M_k: y_i=\sum_{j=0}^{k-1} \beta_{k, j+1} x_i^j+u_i $$ where $k=1, \ldots,\left\lfloor \ln\left(N\right)\right\rfloor$ and $u_i$ is assumed to be i.i.d. $N\left(0, \sigma^2\right)$. When $k \rightarrow \infty$ as $N \rightarrow \infty$, the polynomial regression is related to the sieve estimator which uses progressively more complex models to estimate an unknown function as more data becomes available. In our experiment, we estimate and compare all the candidate models $\left\{M_k, k=1, \ldots,\left\lfloor \ln\left(n^{3/4}\right) \right\rfloor\right\}$. Let $\mathrm{x}^j=\left(x_1^j, x_2^j, \ldots, x_N^j\right)^{\prime}, \mathbf{X}_k=\left(\mathrm{x}^0, \mathrm{x}^1, \ldots, \mathrm{x}^{k-1}\right)$, and $\mathbf{X}=\left(\mathrm{x}^0, \mathrm{x}^1, \ldots, \mathrm{x}^{\left[\ln(N)\right]-1}\right)$. In $M_k$, function $f\left(\boldsymbol{\beta}_k,\mathbf{X}_k\right) = \sum_{j=0}^{k-1} \beta_{k, j+1} x_i^j$ is used to approximate $\ln \left(1+46 x_i\right)$. Let $\boldsymbol{\beta}_k=\left(\beta_1, \ldots, \beta_k\right)^{\prime}$ so that $\boldsymbol{\theta}_k=\left(\boldsymbol{\beta}_k^{\prime}, \sigma^2\right)$, and the number of parameters is $k+1$.

For Bayesian analysis, we assign priors to $\boldsymbol{\beta}_k$ and $\sigma^2 = h^{-1}$ as follows: \[ \boldsymbol{\beta}_k \sim N(\tilde{\mu}, h^{-1} \tilde{V}), \quad h \sim \text{Gamma}(a, b), \] the hyperparameters of the priors are set as \(a = 1\), \(b = 1\), \(\tilde{\mu} = \mathbf{0}\), and \(\tilde{V} = 10^5 \times \mathbf{I}_k\). The optimal VB posterior of $\beta$ and $h$, which is $q(\beta, h)=q(\beta)q(h)$, approximats the true posterior $p\left(\beta,h|y\right)$, see more details in appendix B.1.

In the simulation study, the sample size varies from \(N = 500\) to \(N = 1,000,000\). For each sample size, we simulate the DGP 1000 times. In the \(i\)-th replication, a dataset of size \(N\) is simulated, and the values of \(\text{VPIC}^k\), \(\text{VDIC}_M^k\), \(\text{ELBO}^k\), \(\text{AIC}^k\) and \(\text{BIC}^k\) are computed for the candidate models \(M_k, k = 1, \dots, \left\lfloor \ln\left(N\right) \right\rfloor\).

The relative frequencies of the selected models by each of three criteria (namely VPIC, $\text{VDIC}_M$, ELBO, AIC and BIC) are reported in Figure (ref). And the average values of $k^*$ is listed in the table (ref), all across 1,000 replications. Several interesting results can be found in Figure (ref). First, the models selected by the ELBO and $\text{BIC}$ tend to be parsimonious than those selected by $\text{VPIC}$, $\text{VDIC}_M$ and $\text{AIC}$, this result is not surprising as BIC has a larger penalty term than AIC. Second, as $N$ increases, the average $k^*$ s selected by VPIC and $\text{VDIC}_M$ tends to be similar, suggested that they tend to select the same model. Though under regular conditions, the difference between $\text{BIC}$ and $\text{ELBO}$ are constant as $N$ goes infinity, in our simulation, we find that the averge $k^*$ selected by of $\text{BIC}$ and $\text{ELBO}$ tends to be different as $N$ increases. Third, as the sample size increases, the average $k^*$ s selected by all criteria tend to increase. This is not surprising as the true DGP is not a candidate model.

Table (ref) reports the results of risks. We report $(\widehat{\operatorname{Risk}}-1-\ln (2 \pi))$ scaled by $10^3$ instead of $\widehat{\operatorname{Risk}}$ to better highlight differences in the risks under different criteria. We focus on the risk of $\mathrm{VPIC}$, $\mathrm{VDIC}_M$, $\mathrm{AIC}$ and two proxies of risk $\mathrm{ELBO}$ and $\mathrm{BIC}$. In our simulation experiment, $\mathrm{\text{VDIC}_M}$ and $\mathrm{\text{VPIC}}$ have smaller risks than $\mathrm{ELBO}$, $\mathrm{AIC}$ and $\mathrm{BIC}$. The most important result from Table (ref) is that $\text{VPIC}$ leads to a much smaller value of the expected KL divergence than the other criteria. Results obtained from this Monte Carlo study indicate that if one's objective is to get a best prediction for the future data, we should not only consider how to choose a "best" model and estimator the parametric in this model. We should take predictive distribution into consideration, that means we should use $\mathrm{\text{VDIC}_M}$ and $\mathrm{VPIC}$, not only one criterion, compare these criteria and get a minimum. Then we choose this "optimal" model and use the corresponding predictive distribution to predict the future data.

Probit Regression

In this subsection, we report a generalized linear model (GLM) example, using probit regression. We have linear predictor $Z_i=X_i^{\prime} \beta$ based on vector $X_i$, and we choose the probit link $g\left(E\left[Y_i|X_i\right]\right)=g\left(p_i\right)=Z_i$ as link function, the inverse of the link function $g^{-1}\left(\cdot\right)=\Phi\left(\cdot\right)$ is the cumulative distribution function (cdf) of standard normal distribution, it is shown that

where $\beta$ is $p\times1$ vector. For the Bayesian analysis, we assume a normal prior \(\beta \sim N(\tilde{\mu}, \tilde{V})\), where $\tilde{\mu} = \mathbf{0}$ and $\tilde{V} = 10^5 \times \mathbf{I}_p$, then employ the mean-field VB method to derive the optimal VB posterior distribution \(q(\beta)\). For further details, refer to Appendix B.2

In this simulation study, we define the DGP as $p=4$, $\beta = \left(\beta_0, \beta_1, \beta_2, \beta_3\right)^{\prime}$ with $\beta_0=-0.2$, $\beta_1=0.3$, $ \beta_2=0$, $\beta_3 = 0.7$, $X_i = \left(1, x_{i1},x_{i2}, x_{i3}\right)^{\prime}$, and $ N$ ranging from $N=500$ to $N=1,000,000$. We define such a model to simulate the scenario of under-fitting and over-fitting. Similarly to the first simulation study, we consider seven candidate models, as detailed below,

We replicated DGP for 1000 times, in the $i^{th}$ replication, we generate the data with sample size $N$, and calculate $\text{VPIC}^k$, $\text{VDIC}_M^k$, $\text{ELBO}^k$, $\text{AIC}^k$ and $\text{BIC}^k$ with $M_k=M_1,\dots,M_7$. Then we compare the performance of these criteria.

Table (ref) presents the average risk associated with two different information criteria for seven candidate models under \(N\) ranging from 500 to 1,000,000. Each column in the table reports the risk when choosing the optimal candidate model \(M_{k^*}\). The risk of \(\text{VPIC}^k\) is consistently lower than that of \(\text{VDIC}_M^k\). Also, same like the results in the first simulation, our Monte Carlo experiment has shown that predictive risk under choosing the optimal candidate model from Variational predictive distribution is lower than that from variational lower bound.

Empirical Studies

In this subsection, we first analyze a linear model with different covariates to identify the model that best predicts the number of passengers transported by flight. In the second study, we examine a credit risk model, typically formulated as a binary classification problem. These real data studies aim to show the performance of our two proposed new criteria, and to present that these VB based information criteria can well behave under big data analysis.

US Domestic Flights Predictive Model

In this section, we analyze a linear model with different covariates to identify the model that best predicts the number of passengers transported by flight. The data set used in this analysis pertains to US domestic flights from 1990 to 2009 and contains approximately \(N = 3.61\) million observations. This data set is publicly available on Kaggle. Chasiotis_Karlis_2024 employed this dataset to fit a linear regression model, selecting \(p = 5\) measurements as covariates. In this study, we utilize linear regression to explore the relationships between the dependent variable PASSENGERS (number of passengers, $y$) and the selected covariates, including SEATES (number of seats available on flight, $x_1$), FLIGHTS (number of flights between two locations, $x_2$), DISTANCE (distance flown between origin and destination, $x_3$), ORIGIN POP (origination city's population, $x_4$), DESTINATION POP (destination city's population, $x_5$), ORIGIN LONG (origination airport longitude, $x_6$), DESTINATION LONG (destination airport longitude, $x_7$), ORIGIN LAT (origination airport latitude, $x_8$), DESTINATION LAT (origination airport latitude, $x_9$). To conduct the model selection problem of this dataset, we consider four candidate models, and we list the candidate models and related considerations.

Table (ref) lists the variables we use in the linear regression model. For model comparison, we use mean-field VB to obtain the variational posterior estimators, and then compute the two new proposed information criteria $\text{VPIC}^{k}$ and $\text{VDIC}_M^{k}$ for all candidate models. In choosing the optimal model, to compare the performance of our new proposed method, with other commonly used criteria, we also report $\text{ELBO}^k$, $\text{AIC}^{k}$, $\text{BIC}^k$, $\text{DIC}^k$ and $\text{DIC}^k_M$.

Table (ref) presents the values of \(\text{VPIC}^k\) and \(\text{VDIC}_M^k\), along with $\text{ELBO}^k$ and conventional (or benchmark) information criteria like $\text{AIC}^k$, $\text{BIC}^k$, $\text{DIC}$ and $\text{DIC}_M$. For the candidate models \(\{M_k\}_{k=1}^4\). Importantly, both \(\text{VPIC}^k\) and \(\text{VDIC}_M^k\) select model \(M_4\), same as the benchmark information criteria, indicating that \(M_4\) is preferred over the other candidate models under the same criteria. Based on these results, we recommend selecting model \(M_4\) and using the VB posterior predictive distribution for decision-making to achieve the minimum predictive risk.

To show the difference among seven information criteria, we report a more detailed summary, shown in table (ref). As is known, common information criterion are constructed with two terms: one is the fit term $D(\boldsymbol{\theta})$ equals $-2*\ell\left(\boldsymbol{\theta}\right)$, where $\ell\left(\boldsymbol{\theta}\right)$ is the logarithm likelihood function, and penalty term $P_{IC}$ depend on different measures. If one conducts model selection under Bayes framework, one aims to use the true posterior mean of $\boldsymbol{\theta}$, which is $\bar{\boldsymbol{\theta}}$, or turn to use VB posterior mean $\bar{\boldsymbol{\theta}}^{VB}$ in fit term. Some results can be found in this table. First, as we report in the table, the difference between $\ell\left(\bar{\boldsymbol{\theta}}\right)$ and $\ell\left(\bar{\boldsymbol{\theta}}^{VB}\right)$ is very small, showing that the true posterior mean and VB posterior mean tend to converge to the same value as the size of observed data $N$ goes to infinity. It should also be noted that the inference time between $\bar{\boldsymbol{\theta}}$ and $\bar{\boldsymbol{\theta}}^{VB}$ differs in application, to obtain $\bar{\boldsymbol{\theta}}$ in this around 3 million data, we expend 7159.35 seconds using MCMC, however, 2.54 seconds is used to obtain $\bar{\boldsymbol{\theta}}^{VB}$ under VB as we recorded. As N becomes larger or a more complicated model incoming, one may have to turn to used VB based information criteria rather than using other criteria. Second, both the penalty term $P_{\text{VDIC}_M}$ and $P_{\text{DIC}_M}$ are similar, indicating that $\text{VDIC}_M$ behaves like $\text{DIC}_M$. In addition, $P_{DIC}$ is similar to $P_{AIC}$, as li2024deviance showed that $\text{DIC}$ is a Bayesian version of $\text{AIC}$.

Credit Risk Analysis

The credit risk analysis is an application of binary classification model, including probit regression and logistic regression, used to determine whether a loan should be granted based on various borrower-specific information. In the context of binary classification, we define \(Y_i = 1\) if a loan is approved for the borrower, and \(Y_i = 0\) if it is not. For this study, we utilize the LendingClub dataset, which is publicly available on Kaggle. This data downloaded from Kaggle has about 3 million, and covers the period from 2007 to the third quarter of 2020. By referring filtering process in \href{https://www.kaggle.com/code/kevinkollcaku/loan-classification}{Loan Classification}, we finally got 1.74 million data points. Tabel (ref) lists independent varaible and dependent variables that we are interested in.

We use probit regression and logistic regression to model the factors that affect personal loans, linear combination is $Z_i = \beta_0+\beta_1 \log AnuI_i + \beta_2 Emp_i + \beta_3 DTI_i + \beta_4 \log Loanam_i + \beta_5 Term_i$. Candidate models are $M_k, k=1,2$, which can be listed as $$Y_i|Z_i \overset{i.i.d.}{\sim} Bernoulli\left(\mu \left(Z_i\right)\right),$$ where $\mu\left(Z_i\right)$ differs in $M_1: \mu\left(Z_i\right) =\Phi\left(Z_i\right)$, and $M_2: \mu\left(Z_i\right) =logit\left(Z_i\right)$, where $\Phi\left(\cdot\right)$ is the cumulative density function (CDF) of standard normal distribution, and $logit\left(\cdot\right)$ is the logit link function. For choosing the best model, we use mean-field VB to obtain the variational posterior mean estimator and compute the $\text{VPIC}^{k}$, $\text{VDIC}_M^{k}$, $\text{ELBO}^k$. Benchmark criteria, including $\text{AIC}^k$ and $\text{BIC}^k$ are also calculated for all candidate models. The estimator of two models are reported in table (ref)

Table (ref) presents the values of \(\text{VPIC}^k\) and \(\text{VDIC}_M^k\), along with $\text{AIC}^k$ and $\text{BIC}^k$ for models \(\{M_k\}_{k=1}^2\). The primary differences between the criteria is mainly due to the logarithm likelihood function (or fit term), which is no surprising as the prior of $\boldsymbol{\beta}$ is vauge. Importantly, both \(\text{VPIC}^k\) and \(\text{VDIC}_M^k\) identify model \(M_1\) as the best among the candidate models, indicating its superiority under these criteria. These VB based criteria suggest that the probit model is better than the logit model. Based on these findings, we recommend selecting model \(M_1\) and employing the VB posterior predictive distribution for decision-making to minimize predictive risk.

Conclusion

In this paper, we propose two novel penalty-based predictive information criteria for model comparison in the context of misspecified models with massive data. First, leveraging the VB posterior distribution, we demonstrate that two types of predictive distributions can be derived from a predictive perspective: the variational plug-in predictive distribution and the variational posterior predictive distribution. Second, we investigate the risk functions associated with these two variational predictive distributions, which are defined as the expectations of the KL divergence between the DGP and the predictive distributions. Third, under specific regularity conditions, we prove that the proposed information criteria are asymptotically unbiased estimators of their respective risk functions. Finally, through comprehensive numerical simulations and empirical applications in the fields of economics and finance, we demonstrate the performance of the proposed information criteria for model comparison of misspecified models in the context of massive data.