The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
131,337 characters
Comparing Misspecified Models with Big Data: A Variational Bayesian Perspective
\def\spacingset#1{\renewcommand{\baselinestretch}
{#1}\small\normalsize} \spacingset{1}
\if11
{
\title{\bf Comparing Misspecified Models with Big Data: A Variational Bayesian Perspective}
\author{Yong Li\thanks{
Li gratefully acknowledges \textit{the financial support of the Chinese Natural
Science fund (No.72273142, 72394392).} Yong Li, School of Economics, Renmin University
of China, Beijing, 1000872, China.}\hspace{.2cm}\\
Renmin University of China\\
[email removed]\\
\and
Sushanta K. Mallick \thanks{Sushanta K. Mallick, School of Business and Management, Queen Mary University of London}
\hspace{.2cm}\\
Queen Mary University of London\\
[email removed]\\
\and
Tao Zeng\thanks{Tao gratefully acknowledges \textit{the financial support of the National Natural
Science Foundation of China (No.72073121).} Tao Zeng, School of Economics, Academy of Financial Research and Beijing Research Center, Zhejiang
University, Zhejiang, 310058, China.}\\
Zhejiang University \\
[email removed]\\
\and
Junxing Zhang\thanks{Junxing Zhang, School of Economics,
Renmin University of China, Beijing, 1000872, China.}\\
Renmin University of China\\
[email removed]}
\maketitle
} \fi
\if01
{
\title{\bf Comparing Misspecified Models with Big Data, A Variational Bayesian Perspective}
\maketitle
\medskip
} \fi
\begin{abstract}
Optimal data detection in massive multiple-input multiple-output (MIMO) systems often requires prohibitively high computational complexity. A variety of detection algorithms have been proposed in the literature, offering different trade-offs between complexity and detection performance. In recent years, Variational Bayes (VB) has emerged as a widely used method for addressing statistical inference in the context of massive data. This study focuses on misspecified models and examines the risk functions associated with predictive distributions derived from variational posterior distributions. These risk functions, defined as the expectation of the Kullback-Leibler (KL) divergence between the true data-generating density and the variational predictive distributions, provide a framework for assessing predictive performance. We propose two novel information criteria for predictive model comparison based on these risk functions. Under certain regularity conditions, we demonstrate that the proposed information criteria are asymptotically unbiased estimators of their respective risk functions. Through comprehensive numerical simulations and empirical applications in economics and finance, we demonstrate the effectiveness of these information criteria in comparing misspecified models in the context of massive data.
\end{abstract}
\noindent
{\it Keywords:} Kullback–Leibler divergence; Information criterion; Model misspecification; Variational Bayes; Massive Data.
\spacingset{1.9}
\section{Introduction}
\label{sec:intro}
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.
\begin{comment}
One of the most important statistical inferential issues is model
comparison; see \citet{granger1995comments}, \citet{phillips1994posterior},\citet{phillips1995bayesian,phillips1996econometric},
\citet{hansen2005test}, \citet{burnham2008model} for a partial list of studies. From
a predictive point of view, many penalty-based information criteria have been
proposed for model comparison. In the frequentist framework, the
two most popular information criteria are \text{AIC} of \citet{Akaike1973} and
\text{TIC} of \citet{takeuchi1976distribution}. They are both asymptotically unbiased
estimators of the risk, which is the expectation of \text{KL} loss function
based on the maximum likelihood estimator(MLE). \text{AIC} generally assumes that all candidate models
nest the true model or are good approximations of the data generating process
(DGP). \text{TIC} allows candidate models to be misspecified and its penalty
term involves the inverse of the Hessian matrix. In the Bayesian framework,
in addition to Bayes factors (BFs) that compare models through their "posterior
probabilities" and try to find the "true" model, arguably the most
popular penalty-based information criterion is DIC proposed by \citet{spiegelhalter2002bayesian}. In a recent study, \citet{li2020deviance} develop a variant of DIC for comparing misspecified models. \citet{li2024deviance} propose a decision
theoretic explanation for DIC and show that \text{DIC} is the Bayesian
version of AIC.,As in literatures above, most information criterias under Baysian framework
are based on the output of Markov Chain Monte Carlo (MCMC) samples,
howerver, under the condition of massive data, one may face a dataset with
more that million observations, which is difficult in using MCMC to conduct
Baysian Modelling. As a alternative method, Variational Bayes (VB) can act
as a fast computational method for approximating posterior distributions by
optimize technique, especially under massive data. Despite VB has been
widely used to obtain approximated paramters in Bayesian statistics, it is
until recent years that some general theory have been proposed to explain
from a frequentist perspective why VB works so well \citet{Zhang_2024}.
\end{comment}
Model comparison is one of the most critical issues in statistical inference. For a partial list of studies, see \citet{granger1995comments}, \citet{phillips1994posterior}, \citet{phillips1995bayesian,phillips1996econometric}, \citet{hansen2005test}, and \citet{burnham2008model}. There are essentially two strands of literature on model selection \citep{VehtariandOjanen,anderson2004model}. The first strand aims to answer the question which model best explains the observed data. The Bayes factor (BF, \citealp{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, \citealp{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 \citet{Akaike1973} and the Takeuchi Information Criterion (TIC) introduced by \citet{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 \citet{spiegelhalter2002bayesian}, is one of the most popular penalty-based predictive information criteria. In a recent study, \citet{li2020deviance} developed a variant of DIC for comparing misspecified models, while \citet{li2024deviance} proposed a decision-theoretic interpretation of DIC, demonstrating that DIC is the Bayesian version of AIC.
\begin{comment}
Some model selection approaches using Variational Bayes method have been
introduced in recent years, commonly for Variational Bayes, one way is to
use the evidence lower bound (ELBO) as a proxy of logrithm likelihood
function $\log p\left(\mathbf{y}\right)$ to conduct BFs model comparision. \citet{corduneanu2001variational} studied on VB model selection on Mixture
models, and for the purpose of choosing the best number of components,they
use ELBO as a proxy of marginal likelihood to conduct model selection on
choosing the best number of component of mixture model, and they use a
Gaussian Mixture Model to show that by setting the mixing coefficients to
maximize the marginal log-likelihood, unwanted components can be suppressed,
and the appropriate number of components for the mixture can be determined
in a single training run without recourse to cross validation. \citet{you2014variational} studied Variational Bayes for a classical Bayesian
linear model, they proved that under mild regularity conditions, VB based
estimators enjoy some desirable frequentist properties such as consistency.
They also proposed two VB information criteria: the Variational AIC, which
substitue VB posterior mean in DIC, and the Variational BIC, which uses ELBO
as a proxy of marginal likelihood., then showed that VAIC is asymptotically
equivalent to the frequentist AIC and that the VBIC is first
order equivalent to the BIC in linear regression. \citet{Zhang_2024}
proposed ELBO as an alternative ctiterion for model selection and showed
that it is asymptotically equivalent to the Bayesian information criterion
(BIC). For misspecified models, however, there are few literatures talking
about Baysian model selection from predictive viewpoint in the current era
of massive data.
\end{comment}
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. \citet{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. \citet{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. \citet{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.
\begin{comment}
In this paper, for misspecified models with massive data, we propose two new
penalty-based predictive information criterions for model comparison. First,
based on variational posterior distribution, we show that from predictive
viewpoint, two kinds of predictive distributions, that is, variational
plug-in predictive distribution and variational posterior predictive
distribution can be obtained. Second, we investigate the risk functions for these two variational predictive distributions which are the expectation of the K-L divergence between the true density and the predictive distributions. Thrid, under some regularity conditions, we derive that these proposed two new information criterion are respectively
asymptotically unbiased estimator of the corresponding risk function.
At last, based on simulation and real studies, we illustrate the proposed information criterion using some popular examples in economics and finance.
\end{comment}
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{sec:review} briefly reviews the literature
on how to make statistical inferences about misspecified models and VB
technique for misspecified models with massive data. Section \ref{sec:risk} investigates
the risk functions of variational predictive distributions. Section \ref{sec:model selection}
introduces the statistical decision theory and proposes the new
penalized-based information criterion to compare misspecified models with
massive data. Section \ref{sec:simu and real} illustrates the new methods using two simulated big
data and two real big data. Section \ref{sec:conclude} concludes the paper. The Appendix
collects the proof of the theoretical results and VB analytical expression
of parametric models used in the paper.
\section{Statistical Inference for Misspecified Models: A Review} \label{sec:review}
\subsection{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,
\begin{eqnarray*}
& &KL[g(\mathbf{y}),p(\mathbf{y}|{\mbox{\boldmath${\theta}$}})]=\int g(
\mathbf{y})\ln \frac{g(\mathbf{y})}{p(\mathbf{y}|{\mbox{\boldmath${\theta}$}}
)}d\mathbf{y} \\
&=&E_{g(\mathbf{y})}\ln g(\mathbf{y})-E_{g(\mathbf{y})}\ln p(\mathbf{y}|{\mbox{\boldmath${
\theta}$}}),
\end{eqnarray*}
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
\begin{equation*}
{\mbox{\boldmath${\theta}$}}^{*}=\arg \min_{\theta }KL({\mbox{\boldmath${
\theta}$}})=\arg \max_{\theta }E_{g(\mathbf{y})}\ln p(\mathbf{y}|{\mbox{\boldmath${
\theta}$}}),
\end{equation*}
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,
\begin{equation*}
{\mbox{\boldmath${\hat{\theta}}$}}=\arg \max_{\theta }\ln p(\mathbf{y}|{
\mbox{\boldmath${\theta}$}}).
\end{equation*}
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,
\begin{equation}
\left(\mathbf{\hat{I}}^{-1}({\mbox{\boldmath${\hat{\theta}}$}})\mathbf{\hat{J}}({\mbox{\boldmath${\hat{\theta}}$}})\mathbf{\hat{I}}^{-1}({\mbox{\boldmath${\hat{\theta}}$}})\right)^{-1/2} \sqrt{n}({\mbox{\boldmath${\hat{\theta}}$}}-{\mbox{\boldmath${\theta}$}}^{*})\overset{d}{\rightarrow }N\left({\mbox{\boldmath${0}$}},\mathbf{I}\right) , \label{white}
\end{equation}
as n goes to infinity where the asymptotic variance takes the sandwich form. If the model is
correctly specified, then
\begin{equation}
\left(\mathbf{\hat{I}}^{-1}( {\mbox{\boldmath${\hat\theta}$}})\right)^{-1/2} \sqrt{n}({\mbox{\boldmath${\hat{\theta}}$}}-{\mbox{\boldmath${\theta}$}}^*)
\overset{d}{\rightarrow }N\left(\mathbf{0},\mathbf{I}\right) . \label{ml}
\end{equation}
as n goes to infinity.
\begin{comment}
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 \(p(\mathbf{y}|\boldsymbol{\theta})\) used to fit the data, where \(\boldsymbol{\theta}\) is a \(P\)-dimensional parameter vector, and \(\boldsymbol{\theta} \in \boldsymbol{\Theta} \subseteq \mathbb{R}^P\). The Kullback-Leibler (KL) divergence is commonly used to measure the "distance" between \(g(\mathbf{y})\) and \(p(\mathbf{y}|\boldsymbol{\theta})\), defined as:
\[
KL[g(\mathbf{y}), p(\mathbf{y}|\boldsymbol{\theta})] = \int g(\mathbf{y}) \ln \frac{g(\mathbf{y})}{p(\mathbf{y}|\boldsymbol{\theta})} d\mathbf{y} = E_{g(\mathbf{y})}\ln g(\mathbf{y}) - E_{g(\mathbf{y})}\ln p(\mathbf{y}|\boldsymbol{\theta}),
\]
where \(E_{g(\mathbf{y})}\) represents the expectation with respect to \(g(\mathbf{y})\). Let \(\boldsymbol{\theta}^* \in \boldsymbol{\Theta} \subset \mathbb{R}^P\) denote the pseudo-true value that minimizes the KL divergence:
\[
\boldsymbol{\theta}^* = \arg \min_{\theta} KL(\boldsymbol{\theta}) = \arg \max_{\theta} E_{g(\mathbf{y})} \ln p(\mathbf{y}|\boldsymbol{\theta}),
\]
and let \(\boldsymbol{\hat{\theta}}\) denote the quasi-maximum likelihood (QML) estimator of \(\boldsymbol{\theta}\), which maximizes the log-likelihood function of the parametric model:
\[
\boldsymbol{\hat{\theta}} = \arg \max_{\theta} \ln p(\mathbf{y}|\boldsymbol{\theta}).
\]
For simplicity, let \(l_t(\mathbf{y}_t, \boldsymbol{\theta}) = \ln p(\mathbf{y}_t | \boldsymbol{\theta})\) represent the conditional log-likelihood for the \(t\)-th observation (\(1 \leq t \leq n\)). We suppress this as \(l_t(\boldsymbol{\theta})\), so 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\)-th derivative of \(l_t(\boldsymbol{\theta})\), with \(\nabla^j l_t(\boldsymbol{\theta}) = l_t(\boldsymbol{\theta})\) when \(j=0\). Further, let:
\[
\mathbf{\hat{J}}(\boldsymbol{\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}}(\boldsymbol{\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, stating that:
\[
\left(\mathbf{\hat{I}}^{-1}(\boldsymbol{\hat{\theta}})\mathbf{\hat{J}}(\boldsymbol{\hat{\theta}})\mathbf{\hat{I}}^{-1}(\boldsymbol{\hat{\theta}})\right)^{-1/2} \sqrt{n}(\boldsymbol{\hat{\theta}} - \boldsymbol{\theta}^*) \overset{d}{\rightarrow} N\left(\boldsymbol{0}, \mathbf{I}\right),
\]
as \(n \to \infty\), where the asymptotic variance takes the sandwich form. If the model is correctly specified, then:
\[
\left(\mathbf{\hat{I}}^{-1}(\boldsymbol{\hat{\theta}})\right)^{-1/2} \sqrt{n}(\boldsymbol{\hat{\theta}} - \boldsymbol{\theta}^*) \overset{d}{\rightarrow} N\left(\boldsymbol{0}, \mathbf{I}\right),
\]
as \(n \to \infty\).
\end{comment}
\subsection{Bayesian Inference under Model Misspecification}
\begin{comment}
Consider a statistical model indexed by a set of $P$ parameters, ${\boldsymbol{\theta}} \in \mathbf{
\Theta }\subseteq R^{P}$, with a prior distribution $p({
\mbox{\boldmath${\theta}$}})$ defined over ${
\mbox{\boldmath${\theta}$}}$. By applying the Bayes' theorem, the posterior distribution can be obtained via
\begin{equation}
p({\mbox{\boldmath${\theta}$}} |\mathbf{y})=\frac{p(\mathbf{y}|{
\mbox{\boldmath${\theta}$}} )p({\mbox{\boldmath${\theta}$}} )}{p(\mathbf{y})}
\propto p({\mbox{\boldmath${\theta}$}} )p(\mathbf{y}|{\mbox{\boldmath${
\theta}$}} ), \label{post01}
\end{equation}
where $p(\mathbf{y})=\int p(\mathbf{y}|{\mbox{\boldmath${\theta}$}} )p({
\mbox{\boldmath${\theta}$}} )d{\mbox{\boldmath${\theta}$}} $ represents the marginal
likelihood.
Since the posterior distribution $p(\boldsymbol{\theta}|\mathbf{y})$ doesn't have analytical form in most cases, the posterior sampling generally can be realized via MCMC techniques \citep{gelman2003bayesian}. Based on the random observations obtained from the posterior sampling, the Bayesian statistical inference can be established via the corresponding sampling means and covariance matrix of the
generated observations. For example, let
$\{\boldsymbol{\theta}^{(j)},j=1,2,\cdots, J\}$ be the efficient random observations generated from posterior distribution after discarding the
burn-in samples of some length. After that, Bayesian estimates of
$\boldsymbol{\theta}$ and the corresponding standard error can be obtained via,
\begin{eqnarray*}
\boldsymbol{\bar\theta}=\frac {1}{J}\sum_{j=1}^J\boldsymbol{\theta}^{(j)},
\quad \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}.
\end{eqnarray*}
Clearly, these Bayesian estimates are consistent estimators of the
corresponding posterior means and variance matrix.In the literature, it is well-documented that MCMC techniques are very powerful and efficient approaches for posterior simulation. Just Due to the advance of MCMC techniques, Bayesin approaches won the popularity for statistical inference, hence, can be applied into a variety of complex models.
It is noted that for mssspecified models, Bayesian large sample theory has an important difference with ML large sample theory. Specifically,
unlike the ML theory, there is no difference in the Bayesian
asymptotic theory when the model is correctly specified and if the model is misspecified. In both cases, the Bayesian large sample theory is given by
\begin{equation*}
\left(\mathbf{\hat{I}}^{-1}(\hat{{\mbox{\boldmath${\theta}$}}})/n\right)^{-1/2} ({\mbox{\boldmath${\theta}$}}-\hat{{\mbox{\boldmath${\theta}$}}})|\mathbf{y}\overset{d}{\rightarrow }N\left(\mathbf{0},\mathbf{I}\right)
\end{equation*}
in probability as n goes to infinity \citep{kleijn2012bernstein}.
\end{comment}
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:
\begin{equation}
p(\boldsymbol{\theta}|\mathbf{y}) = \frac{p(\mathbf{y}|\boldsymbol{\theta})p(\boldsymbol{\theta})}{p(\mathbf{y})} \propto p(\boldsymbol{\theta})p(\mathbf{y}|\boldsymbol{\theta}), \label{post01}
\end{equation}
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 \citep{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:
\begin{equation*}
\left(\mathbf{\hat{I}}^{-1}(\hat{\boldsymbol{\theta}})/n\right)^{-1/2}(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}})|\mathbf{y} \overset{d}{\rightarrow} N\left(\mathbf{0}, \mathbf{I}\right),
\end{equation*}
in probability as $n \to \infty$ \citep{kleijn2012bernstein}.
\subsection{Variational Bayes for Misspecified Models with Massive Data}
\begin{comment}
To obtain the $p({\boldsymbol{\theta}}|\mathbf{y})$, the dominant paradigm in Bayes statistics is MCMC, including Metropolis-Hastings algorithm \citep{metropolis1953equation,hastings1970monte}, and the Gibbs sampler \citep{geman1984stochastic} et al. Although MCMC provides an easy and universal method to sample from Markov chain on $\mathbf{\theta}$, it still faces bottleneck, epsecially in massive data sets (the sample size n is very large). An important situation where the log-likelihood is computationally intractable is Massive Data \citep{bardenet2017markov,quiroz2019speeding}, where the log-likelihood function, is a sum of a very large number of terms and thus too expensive to compute. For the reasons of high computational cost of likelihood evaluation when dealing with massive datasets, MCMC may need thousands minutes or days to reach a stationary posterior distribution.
\end{comment}
To compute \( p(\boldsymbol{\theta}|\mathbf{y}) \), the dominant paradigm in Bayesian statistics is MCMC, including the Metropolis-Hastings algorithm \citep{metropolis1953equation, hastings1970monte} and the Gibbs sampler \citep{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 \citep{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.
\begin{comment}
In recent years, to overcome the disadvantages of MCMC-based Bayesian inference, as an alternative approach, Variational Bayes (VB) proposed by \citet{jordan1999introduction} gained a lot of attention for researchers in the literature. For convenience of presentation, here, we let $\mathbf{y}=\left(y_1,\dots,y_n\right)$ be the independently distributed observed data. As an alternative way, it is noted that VB solves the following optimization problem
\begin{equation*}
p^{VB}({\boldsymbol{\theta}}|\mathbf{y})=\arg \min_{q\left( {\boldsymbol{
\theta}}\right) \in \Gamma }KL\left[ q\left( {\boldsymbol{\theta}}\right) ,p({\boldsymbol{\theta}}|\mathbf{y})\right]
\end{equation*}
to look for an approximation for $p({\boldsymbol{\theta}}|\mathbf{y}
) $ from some friendly variational family $\Gamma $, a commonly used mean field (MF) family satisfies
$q\left( {\boldsymbol{\theta}}\right)=\prod_{i=1}^{P}q_{{\mathbf{\theta }}}\left( \theta _{i}\right)$,
where $p^{VB}({\boldsymbol{\theta}}|\mathbf{y})$ is denoted as Variational Bayesian posterior.
Because VB solves a problem of optimization,it is a computationally efficient alternative to MCMC, especially in the Massive Data context under Bayesian modelling \citep{attias2013inferring,bishop2006pattern} . Literatures have shown that, compared with MCMC, VB-based algorithm is orders of magnitude faster than MCMC \citep{blei2017variational,gunawan2017fast}. Beyond the classical mean-field VB, variational method based on stochastic optimization \citep{hoffman2013stochastic} has made scalable data analysis available.
The property of VB posterior is of interest in literature, we define the second order derivative of the logrithm of the likelihood as
$\mathbf{\bar{H}}_{n}({\mbox{\boldmath${\theta}$}}):=\frac{1}{n}\frac{
\partial ^{2}\ln p(\mathbf{y,}{\boldsymbol{\theta }})}{\partial {
\boldsymbol{\theta }}\partial {\boldsymbol{\theta }}^{\prime }}$, and
take expectation to get
$\mathbf{H}_{n}({\mbox{\boldmath${\theta}$}}):=E\left[ \mathbf{\bar{H}}_{n}({\mbox{\boldmath${\theta}$}})\right]$.
Let
\begin{equation*}
p^{VBN}({\boldsymbol{\theta}}|\mathbf{y})= \left( 2\pi \right)
^{-P/2}\left\vert -n\mathbf{H}_{n}^{d}\right\vert ^{1/2}\exp \left( -\frac{1}{2}\left( {\boldsymbol{\hat{\theta}}}_{n}-{\boldsymbol{\theta }}\right) ^{\prime}\left( -n\mathbf{H}_{n}^{d}\right) \left( {\boldsymbol{\hat{\theta}}}_{n}-{
\boldsymbol{\theta}}\right) \right) ,
\end{equation*}
where $\mathbf{H}_{n}^{d}$ is a diagonal matrix with the same diagonal elements as in $\mathbf{H}_{n}^{d}$. Following \citet{han2019statistical}, the KL
distance between the posterior distribution $p^{VB}({\boldsymbol{\theta
}}|\mathbf{y})$ and $p^{VBN}({\boldsymbol{\theta }}|\mathbf{y})$ goes to $0$ in probability as $n$ goes to infinity.
\end{comment}
Recently, to address the limitations of Bayesian inference based on MCMC for massive datasets, Variational Bayes (VB) methods \citep{jordan1999introduction}, have garnered significant attention in the research community. VB offers an alternative to MCMC by solving the following optimization problem:
\begin{equation*}
p^{VB}(\boldsymbol{\theta}|\mathbf{y}) = \arg \min_{q(\boldsymbol{\theta}) \in \Gamma} KL[q(\boldsymbol{\theta}), p(\boldsymbol{\theta}|\mathbf{y})],
\end{equation*}
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 \citep{attias2013inferring, bishop2006pattern}. Empirical studies have shown that VB-based algorithms can be orders of magnitude faster than MCMC \citep{blei2017variational, gunawan2017fast}. Beyond the classical mean-field VB, advances such as stochastic variational inference (SVI) \citep{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:
\begin{equation*}
p^{VBN}(\boldsymbol{\theta}|\mathbf{y}) = (2\pi)^{-P/2} \left| -n\mathbf{H}_n^d \right|^{1/2} \exp\left(-\frac{1}{2} (\boldsymbol{\hat{\theta}}_n - \boldsymbol{\theta})^\prime (-n\mathbf{H}_n^d)(\boldsymbol{\hat{\theta}}_n - \boldsymbol{\theta})\right),
\end{equation*}
where \(\mathbf{H}_n^d\) is a diagonal matrix whose diagonal elements match those of \(\mathbf{H}_n\). As established by \citet{han2019statistical} and \citet{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\). \citet{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\).
\section{Risk of Predictive Distributions on Misspecified Models based on Variational Bayes} \label{sec:risk}
\begin{comment}
In this section, we introduce how to evaluate the usefulness of a misspecified model from predictive viewpoint. In the literature, for the misspecified model, a typical method to assess the usefulness of a statistical model is to examine its predictive performance \citep{bernardo1979expected}. Given some future observations $\mathbf{y}_{f}$, it is assumed that some predictive distribution can be denoted as $p_{f}\left(\mathbf{y}_f|\mathbf{y}\right)$. A popular method to quantify the predictive ability with a misspecified model is to evaluate the KL divergence between $g\left(\mathbf{
y}_{f}\right) $ and a predictive distribution $p_f\left(\mathbf{y}_{f}|
\mathbf{y}\right) $ multiplied by $2$, that is,
\begin{eqnarray*}
&&2\times KL\left[ g\left( \mathbf{y}_{f}\right),p_f\left(\mathbf{y}_{f}|
\mathbf{y}\right)\right] =2E_{\mathbf{y}_{f}}\left[ \ln \frac{g\left(
\mathbf{y}_{rep}\right) }{p\left(\mathbf{y}_{f}|\mathbf{y}\right)}\right]
\\
&=&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) \mbox{d}
\mathbf{y}_{f}.
\end{eqnarray*}
Based on the above KL divergence, from the statistical decision theory, we can specify the loss function associated with some decision $d$ as
\begin{equation*}
\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] .
\end{equation*}
where $p\left(\mathbf{y}_{f}|\mathbf{y},d\right)$ is the predictive
density using some decision $d$. Given the loss function, the risk function can be defined as
\citep{good1952rational}
\begin{equation*}
Risk(d)=E_{\mathbf{y}}\left[ \mathcal{L}(\mathbf{y},d)\right] =\int
\mathcal{L}(\mathbf{y},d)g(\mathbf{y})\mbox{d}\mathbf{y},
\end{equation*}
As to Variational Bayesian inference, it is noted that two kinds of predictive distributions can be established for prediction,that is, Variational plug-in predictive distribution and Variational posterior predictive distribution. There two distributions can be described by different statistical decision, hence, lead to two different risk functions. In the following subsection, we need to evaluate these two risk functions and derive some estimators for them. For this purpose, here, we fix some notations and give some mildly regular conditions as follows.
\end{comment}
In the literature, assessing the utility of a misspecified statistical model is typically achieved by examining its predictive performance \citep{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 \citep{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
\begin{eqnarray*}
&&\mathbf{s}(\mathbf{y},{\mbox{\boldmath${\theta}$}}):=\frac{\partial \ln p(
\mathbf{y}|{\mbox{\boldmath${\theta}$}})}{\partial {\mbox{\boldmath${
\theta}$}}}=\sum_{t=1}^{n}\bigtriangledown l_{t}\left( {\mbox{\boldmath${
\theta}$}}\right) ,\;\mathbf{h}(\mathbf{y},{\mbox{\boldmath${\theta}$}}):=
\frac{\partial ^{2}\ln p(\mathbf{y}|{\mbox{\boldmath${\theta}$}})}{\partial {
\mbox{\boldmath${\theta}$}}\partial {\mbox{\boldmath${\theta}$}}^{\prime }}
=\sum_{t=1}^{n}\bigtriangledown ^{2}l_{t}\left( {\mbox{\boldmath${\theta}$}}
\right) , \\
&&\mathbf{s}_{t}({\mbox{\boldmath${\theta}$}}):=\bigtriangledown l_{t}\left(
{\mbox{\boldmath${\theta}$}}\right) ,\;\mathbf{h}_{t}({\mbox{\boldmath${
\theta}$}}):=\bigtriangledown ^{2}l_{t}\left( {\mbox{\boldmath${\theta}$}}
\right) , \\
&&\mathbf{B}_{n}\left( {\mbox{\boldmath${\theta}$}}\right) :=Var\left[ \frac{
1}{\sqrt{n}}\sum_{t=1}^{n}\bigtriangledown l_{t}\left( {\mbox{\boldmath${
\theta}$}}\right) \right] ,\mathbf{\bar{H}}_{n}({\mbox{\boldmath${\theta}$}}
):=\frac{1}{n}\sum_{t=1}^{n}\mathbf{h}_{t}({\mbox{\boldmath${\theta}$}}), \\
&&\mathbf{\bar{J}}_{n}({\mbox{\boldmath${\theta}$}}):=\frac{1}{n}
\sum_{t=1}^{n}\left[ \mathbf{s}_{t}({\mbox{\boldmath${\theta}$}})-\mathbf{
\bar{s}}_{t}({\mbox{\boldmath${\theta}$}})\right] \left[ \mathbf{s}_{t}({
\mbox{\boldmath${\theta}$}})-\mathbf{\bar{s}}_{t}({\mbox{\boldmath${\theta}$}
})\right] ^{\prime },\mathbf{\bar{s}}_{t}({\mbox{\boldmath${\theta}$}})=
\frac{1}{n}\sum_{t=1}^{n}\mathbf{s}_{t}({\mbox{\boldmath${\theta}$}}), \\
&&\mathcal{L}_{n}({\mbox{\boldmath${\theta}$}}):=\ln p({\mbox{
\boldmath${ \theta}$}}|\mathbf{y}),\mathcal{L}_{n}^{(j)}({
\mbox{\boldmath${\theta}$}}):=\partial ^{j}\ln p({\mbox{\boldmath${\theta}$}}
|\mathbf{y})/\partial {\mbox{\boldmath${\theta}$}}^{j}, \\
&&\mathbf{H}_{n}({\mbox{\boldmath${\theta}$}}):=\int \mathbf{\bar{H}}_{n}({
\mbox{\boldmath${\theta}$}})g\left( \mathbf{y}\right) d\mathbf{y},\;\mathbf{J
}_{n}({\mbox{\boldmath${\theta}$}})=\int \mathbf{\bar{J}}_{n}({
\mbox{\boldmath${\theta}$}})g\left( \mathbf{y}\right) d\mathbf{y}.
\end{eqnarray*}
Then, the following regularity conditions can be imposed
\textbf{Assumption 1}: ${\mbox{\boldmath${\Theta}$}}\subset R^{P}$ is
compact.
\textbf{Assumption 2}: The data $\mathbf{y}=(y_{1},\ldots ,y_{n})$ is
independent and identically distributed.
\textbf{Assumption 3: }For all $t$, $l_{t}\left( {\mbox{\boldmath
\ensuremath{{\theta}}}}\right) $\ is eight-times differentiable on $
\mbox{\boldmath\ensuremath{{\Theta}}}$ almost surely.
\textbf{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$.
\textbf{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$.
\textbf{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
\begin{equation*}
\mbox{\boldmath${\theta}$}_{n}^{p}=\arg \min_{{\mbox{\boldmath${\theta}$}}
\in {\mbox{\boldmath${\Theta}$}}}\frac{1}{n}\int \ln \frac{g(\mathbf{y})}{p(
\mathbf{y}|{\mbox{\boldmath${\theta}$}})}g(\mathbf{y})d\mathbf{y,}
\end{equation*}
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$,
\begin{equation}
\lim_{n\rightarrow \infty }\sup \sup_{{\Theta \backslash }N\left(
\mbox{\boldmath${\theta}$}_{n}^{p},\varepsilon \right) }\frac{1}{n}
\sum_{t=1}^{n}\left\{ E\left[ l_{t}\left( {\mbox{\boldmath${\theta}$}}
\right) \right] -E\left[ l_{t}\left( \mbox{\boldmath${\theta}$}
_{n}^{p}\right) \right] \right\} <0, \label{iden}
\end{equation}
where $N\left( \mbox{\boldmath${\theta}$}_{n}^{p},\varepsilon \right) $ is
the open ball of radius $\varepsilon $ around $\mbox{\boldmath${\theta}$}
_{n}^{p}$.
\textbf{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$.
\textbf{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, \citet{gallant1988unified} and
\citet{WOOLDRIDGE19942639}. Assumption 8 is the regular condition for prior density,
see, for example, \citet{li2020deviance}. Assumptions 1-8 are sufficient for the assumptions used by
\citet{Zhang_2024} to develope the asymptotic properties of VB posterior distribution without latent variables.
\subsection{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})$.
\begin{comment}
Following the literature of development of the popular information criterion AIC,TIC, DIC,etc, we assume that there are some future replicated data $\mathbf{y}_{rep}$ which has the same DGP with $\mathbf{y}$. More details about $\mathbf{y}_{rep}$, one can refer to the excellent textbook about model selection, that is, \citet{anderson2004model} and its reference therein. After that, for future data $\mathbf{y}_{rep}$, the VB plug-in predictive distribution can be specified as $p\left(\mathbf{y}_{rep}| \overline{\boldsymbol{\theta}}
^{VB}\right)$. Hence, for this predictive distribution, the loss function associated with the corresponding statistical decision denoted as $d_1$ can be given by
\end{comment}
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 \citet{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:
\begin{equation*}
\mathcal{L}(\mathbf{y},d_1)=2\times KL\left[ g\left( \mathbf{y}
_{rep}\right) ,p\left(\mathbf{y}_{rep}| \overline{\boldsymbol{\theta}}
^{VB}\right) \right].
\end{equation*}
In this context, the risk function can be expressed as:
\begin{eqnarray*}
&&Risk(d_1)=E_{\mathbf{y}}\left[ \mathcal{L}( \mathbf{y}
,d_1)\right]=2\times E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left[ \ln
\frac{g\left( \mathbf{y}_{rep}\right) }{p\left(\mathbf{y}_{rep}| \overline{
\boldsymbol{\theta}}^{VB}\right) }\right]\\
&&=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}| \overline{\boldsymbol{\theta}}^{VB}\right)
\right].
\end{eqnarray*}
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:
\begin{equation*}
Risk(d_1)=C+E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left[ -2\ln
p\left(\mathbf{y}_{rep}| \overline{\boldsymbol{\theta}}^{VB}\right) \right]
\end{equation*}
where $C=E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left[ 2\ln g\left(
\mathbf{y}_{rep}\right) \right]$.
\begin{comment}
From this risk function, we can see that the smaller the $Risk(d_1)$, the better the predictive distribution performs when using $p\left(\mathbf{y}_{rep}| \overline{\boldsymbol{\theta}}^{VB}\right) $ to predict $ \mathbf{y}_{rep}$. Unfortunately, generally, this risk function doesn't have analytical form. Hence, we still need to evaluate the risk function in order to examine the predictive behavior. In the following, we derive an asymptotic expansion of this risk function via the following theorem.
\end{comment}
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.
\begin{theorem}
\label{riskvtic} Under Assumptions 1-8, it can be shown that
\begin{eqnarray*}
&&E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left( -2\ln p\left(\mathbf{y}_{rep}|
\overline{\boldsymbol{\theta}}^{VB}\right) \right)
=E_{\mathbf{y}}\left( -2\ln p\left( \mathbf{y}|\widehat{{\mbox{\boldmath${\theta}$}}}_{n}\left( \mathbf{y}\right) \right) \right) -2
\mathbf{tr}\left[\mathbf{B}_{n}\mathbf{H}_{n}^{-1}\right] +o\left( 1\right).
\end{eqnarray*}
with
$
\mathbf{B}_{n} = \mathbf{B}_{n}\left(\boldsymbol{\theta}_n^p\right), \mathbf{H}_{n} = \mathbf{H}_{n}\left(\boldsymbol{\theta}_n^p\right)
$, where $\widehat{{\mbox{\boldmath${\theta}$}}}_{n}\left( \mathbf{y}\right)$
is the MLE estimator of $\boldsymbol{\theta}$.
\end{theorem}
\begin{remark}
Under Assumptions 1-8, it can be shown that when the model is correctly specified
\begin{eqnarray*}
&&E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left( -2\ln p\left(\mathbf{y}_{rep}|
\overline{\boldsymbol{\theta}}^{VB}\right) \right)
=E_{\mathbf{y}}\left( -2\ln p\left( \mathbf{y}|\widehat{{\mbox{\boldmath${\theta}$}}}_{n}\left( \mathbf{y}\right) \right) \right) -2
\mathbf{tr}\left[\mathbf{B}_{n}\mathbf{H}_{n}^{-1}\right] +o\left( 1\right)\\
&=&E_{\mathbf{y}}\left( -2\ln p\left( \mathbf{y}|\widehat{{\mbox{\boldmath${\theta}$}}}_{n}\left( \mathbf{y}\right) \right) \right) +2
\mathbf{tr}\left[\mathbf{H}_{n}\mathbf{H}_{n}^{-1}\right] +o\left( 1\right)\\
&=&E_{\mathbf{y}}\left( -2\ln p\left( \mathbf{y}|\widehat{{\mbox{\boldmath${\theta}$}}}_{n}\left( \mathbf{y}\right) \right) \right) +2\mathbf{P} +o\left( 1\right).
\end{eqnarray*}
\end{remark}
\subsection{Risk of VB Posterior Predictive Distribution}
\begin{comment}
Under Bayesian framework, for the replicate data $\mathbf{y}_{rep}$, the VB posterior predictive distribution corresponding to $p^{VB}(\boldsymbol{\theta}|\mathbf{y})$ is defined
as
\begin{equation}
p^{VB}(\mathbf{y}_{rep}|\mathbf{y})=\int p(\mathbf{y}_{rep}|{\mbox{\boldmath${
\theta}$}},\mathbf{y})p^{VB}({\mbox{\boldmath${\theta }$}}|\mathbf{y})d{
\mbox{\boldmath${\theta}$}}. \label{pred_VB}
\end{equation}
Like description in section 3.1, for a misspecified model, the KL
divergence between $g\left( \mathbf{y}_{rep}\right) $ and a predictive
distribution $p^{VB}\left( \mathbf{y}_{rep}|\mathbf{y}\right) $ multiplied
by $2$,
\begin{eqnarray*}
&&2\times KL\left[ g\left( \mathbf{y}_{rep}\right) ,p^{VB}\left( \mathbf{y}
_{rep}|\mathbf{y}\right) \right] =2E_{\mathbf{y}_{rep}}\left[ \ln \frac{
g\left( \mathbf{y}_{rep}\right) }{p^{VB}\left( \mathbf{y}_{rep}|\mathbf{y}
\right) }\right] \\
&=&2\int \left[ \ln \frac{g\left( \mathbf{y}_{rep}\right) }{p^{VB}\left(
\mathbf{y}_{rep}|\mathbf{y}\right) }\right] g\left( \mathbf{y}
_{rep}\right) \mbox{d}\mathbf{y}_{rep}
\end{eqnarray*}
is used to quantify the predictive ability. Similarly, the
loss function associated with the statistical decision $d_2$ to use the VB posterior predictive distributions for prediction can be denoted by
\begin{equation*}
\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] .
\end{equation*}
As a result, the corresponding risk function for decision $d_2$ can be written as,
\begin{eqnarray*}
&&Risk(d_2)=E_{\mathbf{y}}\left[ \mathcal{L}(\mathbf{y}
,d_2)\right] \\
&=&2\times E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left[ \ln
\frac{g\left( \mathbf{y}_{rep}\right) }{p^{VB}\left( \mathbf{y}_{rep}|
\mathbf{y}\right) }\right] \\
&=& 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^{VB}\left( \mathbf{y}_{rep}|\mathbf{y}\right) \right]
\end{eqnarray*}
Consequently, we also can have
\begin{equation*}
Risk(d_2)=C+E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left[ -2\ln
p^{VB}\left( \mathbf{y}_{rep}|\mathbf{y}\right) \right] .
\end{equation*}
From above, we can also see that the smaller the $Risk(d_2)$, the better the predictive distribution
performs when using $p^{VB}\left( \mathbf{y}_{rep}|\mathbf{y}\right) $ to
predict $g\left( \mathbf{y}_{rep}\right)$ among different models. In the following, we also derive an asymptotic expansion of this risk function via the following theorem.
\end{comment}
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:
\begin{equation}
p^{VB}(\mathbf{y}_{rep}|\mathbf{y})=\int p(\mathbf{y}_{rep}|{\mbox{\boldmath${
\theta}$}},\mathbf{y})p^{VB}({\mbox{\boldmath${\theta }$}}|\mathbf{y})d{
\mbox{\boldmath${\theta}$}}. \label{pred_VB}
\end{equation}
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:
\begin{eqnarray*}
&&2\times KL\left[ g\left( \mathbf{y}_{rep}\right) ,p^{VB}\left( \mathbf{y}
_{rep}|\mathbf{y}\right) \right] = 2E_{\mathbf{y}_{rep}}\left[ \ln \frac{
g\left( \mathbf{y}_{rep}\right) }{p^{VB}\left( \mathbf{y}_{rep}|\mathbf{y}
\right) }\right] \\
&=&2\int \left[ \ln \frac{g\left( \mathbf{y}_{rep}\right) }{p^{VB}\left(
\mathbf{y}_{rep}|\mathbf{y}\right) }\right] g\left( \mathbf{y}
_{rep}\right) \mbox{d}\mathbf{y}_{rep}
\end{eqnarray*}
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:
\begin{eqnarray*}
&&Risk(d_2)=E_{\mathbf{y}}\left[ \mathcal{L}(\mathbf{y}
,d_2)\right] =2\times E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left[ \ln
\frac{g\left( \mathbf{y}_{rep}\right) }{p^{VB}\left( \mathbf{y}_{rep}|
\mathbf{y}\right) }\right] \\
&=& 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^{VB}\left( \mathbf{y}_{rep}|\mathbf{y}\right) \right],
\end{eqnarray*}
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.
\begin{theorem}
\label{riskvpic} Under Assumptions 1-8, it can be shown
\begin{eqnarray*}
&&E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left( -2\ln p^{VB}\left( \mathbf{y}
_{rep}|\mathbf{y}\right) \right) \\
&=&E_{\mathbf{y}}\left( -2\ln p\left( \mathbf{y}|\widehat{{
\mbox{\boldmath${\theta}$}}}_{n}\left( \mathbf{y}\right) \right) \right)
+\ln \left( \left\vert -\mathbf{H}_{n}\left( -\mathbf{H}_{n}^{d}\right)
^{-1}+\mathbf{I}_{n}\right\vert \right) +\mathbf{tr}\left[ \mathbf{B}
_{n}\left( -\mathbf{H}_{n}\right) ^{-1}\right] \\
&&-\mathbf{tr}\left[ \left( -\mathbf{H}_{n}+\left( -\mathbf{H}
_{n}^{d}\right) \right) ^{-1}\left( \mathbf{B}_{n}+\left( -\mathbf{H}
_{n}^{d}\right) \mathbf{C}_{n}\left( -\mathbf{H}_{n}^{d}\right) \right)
\right] +\mathbf{tr}\left[ \left( -\mathbf{H}_{n}^{d}\right) \mathbf{C}_{n}
\right] +o\left( 1\right)
\end{eqnarray*}
where $\mathbf{C}_{n}=\mathbf{H}_{n}^{-1}\mathbf{B}_{n}\mathbf{H}_{n}^{-1}$,
$\mathbf{H}_{n}^{d}$ is a diagonal matrix with the same diagonal elements as
in $\mathbf{H}_{n}$.
\end{theorem}
\begin{remark}
If $\mathbf{H}_{n}$ is diagonal, that
is $\mathbf{H}_{n}^{d}=\mathbf{H}_{n}$, it can be shown that
\begin{equation}
\ln \left( \left\vert -\mathbf{H}_{n}\left( -\mathbf{H}_{n}^{d}\right) ^{-1}+
\mathbf{I}_{n}\right\vert \right) =\ln \left( \left\vert -\mathbf{H}
_{n}\left( -\mathbf{H}_{n}\right) ^{-1}+\mathbf{I}_{n}\right\vert \right)
=\ln \left( \left\vert 2\mathbf{I}_{n}\right\vert \right) =\mathbf{P}\ln 2,
\end{equation}
and
\begin{equation}
\begin{aligned}
&-\mathbf{tr}\left[ \left( -\mathbf{H}_{n}+\left( -\mathbf{H}
_{n}^{d}\right) \right) ^{-1}\left( \mathbf{B}_{n}+\left( -\mathbf{H}
_{n}^{d}\right) \mathbf{C}_{n}\left( -\mathbf{H}_{n}^{d}\right) \right)
\right] +\mathbf{tr}\left[ \left( -\mathbf{H}_{n}^{d}\right) \mathbf{C}_{n}
\right] \\
&=-\mathbf{tr}\left[ \left( -\mathbf{H}_{n}+\left( -\mathbf{H}_{n}\right)
\right) ^{-1}\left( \mathbf{B}_{n}+\left( -\mathbf{H}_{n}\right) \mathbf{C}
_{n}\left( -\mathbf{H}_{n}\right) \right) \right] +\mathbf{tr}\left[ \left( -
\mathbf{H}_{n}\right) \mathbf{C}_{n}\right] \\
&=-\mathbf{tr}\left[ \left( -2\mathbf{H}_{n}\right) ^{-1}\left( 2\mathbf{B}
_{n}\right) \right] +\mathbf{tr}\left[ \mathbf{B}_{n}\left( -\mathbf{H}
_{n}\right) ^{-1}\right] =0,
\end{aligned}
\end{equation}
then
\begin{eqnarray*}
&&E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left( -2\ln p^{VB}\left( \mathbf{y}
_{rep}|\mathbf{y}\right) \right) \\
&=&E_{\mathbf{y}}\left( -2\ln p\left( \mathbf{y}|\widehat{{
\mbox{\boldmath${\theta}$}}}_{n}\left( \mathbf{y}\right) \right) \right) +
\mathbf{P}\ln 2+\mathbf{tr}\left[ \mathbf{B}_{n}\left( -\mathbf{H}
_{n}\right) ^{-1}\right] +o\left( 1\right) .
\end{eqnarray*}
\end{remark}
\begin{corollary}
Under Assumptions 1-8, it can be shown that when the model is correctly specified
\begin{eqnarray*}
&&E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left( -2\ln p^{VB}\left( \mathbf{y}
_{rep}|\mathbf{y}\right) \right) \\
&=&E_{\mathbf{y}}\left( -2\ln p\left( \mathbf{y}|\widehat{{
\mbox{\boldmath${\theta}$}}}_{n}\left( \mathbf{y}\right) \right) \right)
+\ln \left( \left\vert -\mathbf{H}_{n}\left( -\mathbf{H}_{n}^{d}\right)
^{-1}+\mathbf{I}_{n}\right\vert \right) +\mathbf{P} \\
&&-\mathbf{tr}\left[ \left( -\mathbf{H}_{n}+\left( -\mathbf{H}
_{n}^{d}\right) \right) ^{-1}\left( -\mathbf{H}_{n}+\left( -\mathbf{H}
_{n}^{d}\right) \left( -\mathbf{H}_{n}\right) ^{-1}\left( -\mathbf{H}
_{n}^{d}\right) \right) \right] \\
&&+\mathbf{tr}\left[ \left( -\mathbf{H}_{n}^{d}\right) \left( -\mathbf{H}
_{n}\right) ^{-1}\right] +o\left( 1\right)
\end{eqnarray*}
where $\mathbf{H}_{n}^{d}$ is a diagonal matrix with the same diagonal
elements as in $\mathbf{H}_{n}$.
\end{corollary}
\begin{remark}
If $\mathbf{H}_{n}$ is
diagonal, that is $\mathbf{H}_{n}^{d}=\mathbf{H}_{n}$, then
\begin{eqnarray*}
&&E_{\mathbf{y}}E_{\mathbf{y}_{rep}}\left( -2\ln p^{VB}\left( \mathbf{y}
_{rep}|\mathbf{y}\right) \right) \\
&=&E_{\mathbf{y}}\left( -2\ln p\left( \mathbf{y}|\widehat{{
\mbox{\boldmath${\theta}$}}}_{n}\left( \mathbf{y}\right) \right) \right) +
\mathbf{P}\ln 2+\mathbf{P}+o\left( 1\right) .
\end{eqnarray*}
\end{remark}
\section{Predictive Information Criteria for Comparing
Misspecified Models with Massive Data based on VB}
\label{sec:model selection}
\begin{comment}
In this section, we show how to develop the newly predictive information criteria for model comparison under misspecification with massive data. Based on the risk functions investigated in Section \ref{sec:risk}, Section \ref{subsec:decision theory} introduces the statistical decision theory for model comparison. In Section \ref{subsec:plugin vbic}, we develop the information criterion named as VTIC on the basis of VB plug-in predictive distribution. In Section \ref{subsec:posterior vbic}, we develop the information criterion named as VPIC on the basis of VB posterior predictive distribution.
\end{comment}
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{sec:risk}, Section \ref{subsec:decision theory} introduces the framework of statistical decision theory for model comparison. In Section \ref{subsec:plugin vbic}, 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{subsec:posterior vbic}. At last, in Section \ref{subsec: bf and bic}, we then discuss BFs and BIC in the context of misspecified models.
\subsection{Statistical Decision Theory based on Risk Function for Model
Selection}
\label{subsec:decision theory}
\begin{comment}
In this section, from the predictive viewpoint, we extend the decisional
approach of Section \ref{sec:risk} to develop information criteria for model comparison.
Now suppose there are $K$ candidate models that are all misspecified and we
have to select a model. Denote these candidate models by $M_{k}$, $
k=1,2,\cdots ,K$. As argued in the last section, we do so by minimizing the
risk of the statistical decision. Following \citet{Akaike1973} and \citet{takeuchi1976distribution}, we assume that parameter ${\mbox{\boldmath${\theta}$}}$ is only
estimated from the sample $\mathbf{y}$. Unlike Akaike (1973) and \citet{takeuchi1976distribution}, we do not plug the MLE into the KL divergence function.
\end{comment}
In this section, from a predictive perspective, we extend the decisional framework introduced in Section \ref{sec:risk} 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.
\begin{comment}
Following the foundational work of \citet{Akaike1973} and \citet{takeuchi1976distribution}, we assume that the parameter \(\boldsymbol{\theta}\) is estimated solely from the observed sample \(\mathbf{y}\). However, unlike \citet{Akaike1973} and \citet{takeuchi1976distribution}, we do not substitute the maximum likelihood estimator (MLE) directly into the Kullback-Leibler (KL) divergence function.
\end{comment}
\begin{comment}
Assuming that the probabilistic behavior of observed data, $\mathbf{y}\in
\mathbf{Y,}$ is described by a set of probabilistic models such as $\left\{
M_{k}\right\} _{k=1}^{K}:=\left\{ p\left( \mathbf{y}|{\mbox{\boldmath${
\theta}$}}_{k},M_{k}\right) \right\} _{k=1}^{K}$ where parameter ${
\mbox{\boldmath${\theta}$}}_{k}$ is the set of parameters in model $M_{k}$.
Formally, the model selection problem can be taken as a decision problem to
select a model among $\left\{ M_{k}\right\} _{k=1}^{K}$ where the action
space has $K$ elements, namely, $\{d_{k}\}_{k=1}^{K}$, where $d_{k}$ means $
M_{k}$ is selected.
\end{comment}
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.
\begin{comment}
For the decision problem, similarly to Section \ref{sec:risk}, a loss function, $\mathcal{
L}(\mathbf{y},d_{k})$, which measures the loss of decision $d_{k}$ as a
function of $\mathbf{y}$, must be specified. Given the loss function, the
risk can be also defined as
\begin{equation*}
Risk(d_{k})=E_{\mathbf{y}}\left[ \mathcal{L}(\mathbf{y},d_{k})\right] =\int
\mathcal{L}(\mathbf{y},d_{k})g(\mathbf{y})\mbox{d}\mathbf{y,}
\end{equation*}
where $g(\mathbf{y})$ is the DGP of $\mathbf{y}$. Hence, the model selection
problem is equivalent to optimizing the statistical decision,
\begin{equation*}
k^{\ast }=\arg \min_{k}Risk(d_{k}).
\end{equation*}
Based on the set of candidate models $\left\{ M_{k}\right\} _{k=1}^{K}$, the
model $M_{k^{\ast }}$ with the decision $d_{k^{\ast }}$ is selected.
\end{comment}
For the decision-making process, as in Section \ref{sec:risk}, 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.
\begin{comment}
Again, the quantity that has been used to measure the quality of the
candidate model in terms of its ability to make predictions is the KL
divergence between $g\left( \mathbf{y}_{rep}\right) $ and some predictive
distribution $p\left( \mathbf{y}_{rep}|\mathbf{y},M_{k}\right) $
multiplied by $2$,
\begin{eqnarray*}
&&2\times KL\left[ g\left( \mathbf{y}_{rep}\right) ,p\left( \mathbf{y}
_{rep}|\mathbf{y},M_{k}\right) \right] =2E_{\mathbf{y}_{rep}}\left[ \ln
\frac{g\left( \mathbf{y}_{rep}\right) }{p\left( \mathbf{y}_{rep}|
\mathbf{y},M_{k}\right) }\right] \\
&=&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) \mbox{d}\mathbf{y}_{rep}.
\end{eqnarray*}
Similarly to Section 3, the loss function associated with decision $d_{k}$
is
\begin{equation*}
\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]
.
\end{equation*}
As a result, the model selection problem is,
\begin{eqnarray*}
k^{\ast } &=&\arg \min_{k}Risk(d_{k})=\arg \min_{k}E_{\mathbf{y}}\left[
\mathcal{L}(\mathbf{y},d_{k})\right] \\
&=&\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\} \\
&=&\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\} .
\end{eqnarray*}
Since $g\left( \mathbf{y}_{rep}\right) $ is the DGP, $E_{\mathbf{y}_{rep}}
\left[ 2\ln g\left( \mathbf{y}_{rep}\right) \right] $ is the same across all
candidate models, and, hence, is dropped from the above equation.
Consequently, we have
\begin{equation*}
k^{\ast }=\arg \min_{k}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]
.
\end{equation*}
The smaller the $Risk(d_{k})$, the better the candidate model performs when
using predictive distribution $p(\mathbf{y}_{rep}|\mathbf{y},M_{k})$ to
predict $g\left( \mathbf{y}_{rep}\right) $. The optimal decision makes it
necessary to evaluate the risk among candidate models.
\end{comment}
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.
\begin{comment}
It is noted that the action space is larger than before. Here, from predictive viewpoint, we not only need to select one model for prediction but also which predictive distribution is used. Denote the action
space by $\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 come 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}},
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\left(\mathbf{y}_{rep} \mid \mathbf{y}, M_k\right)$ is used. Let the KL
divergence function be represented as
\begin{equation*}
\left.\mathcal{L}\left(\mathbf{y}, d_{k^a}\right)=2 \times K L\left[g\left(
\mathbf{y}_{\text {rep }}\right), p\left(\mathbf{y}_{r e p}|\mathbf{y},
d_{k^a}\right)\right)\right]
\end{equation*}
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
\begin{equation*}
Risk\left(d_{k^a}\right)=E_{\mathbf{y}}\left(\mathcal{L}\left(\mathbf{
y}, d_{k^a}\right)\right)=\int \mathcal{L}\left(\mathbf{y}, d_{k^a}\right) g(
\mathbf{y}) \mathrm{d} \mathbf{y}
\end{equation*}
Hence, the model selection problem is equivalent to the following
statistical decision,
\begin{equation}
\min_{a \in\{1,2\}} \min_{k \in\{1, \cdots, K\}} Risk
\left(d_{k^a}\right).
\end{equation}
\end{comment}
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
\begin{equation*}
\left.\mathcal{L}\left(\mathbf{y}, d_{k^a}\right)=2 \times K L\left[g\left(
\mathbf{y}_{\text {rep}}\right), p\left(\mathbf{y}_{r e p}|\mathbf{y},
d_{k^a}\right)\right)\right]
\end{equation*}
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
\begin{equation*}
Risk\left(d_{k^a}\right)=E_{\mathbf{y}}\left(\mathcal{L}\left(\mathbf{
y}, d_{k^a}\right)\right)=\int \mathcal{L}\left(\mathbf{y}, d_{k^a}\right) g(
\mathbf{y}) \mathrm{d} \mathbf{y}.
\end{equation*}
Consequently, the model selection problem is equivalent to solving the following statistical decision problem:
\begin{equation}
\min_{a \in\{1,2\}} \min_{k \in\{1, \cdots, K\}} Risk
\left(d_{k^a}\right).
\end{equation}
\begin{comment}
Now we get the risk formula for model selection based on the statistical decision theory. Because we don not know the DGP process $g\left(\boldsymbol{y}\right)$ and $g\left(\boldsymbol{y}_{rep}\right)$, it will be unable to evaluate the risk of decision $d_{k^a}$. However, it will be helpful to evaluate the risk of decision $d_{k^a}$ by using an asymptotically unbiased estimator of $Risk\left(d_{k^a}\right)$. In literatures, it has been shown that different information criterias proposed is actually the asymptotically unbiased estimator of expectation of the loss function up to a constant under different statistical decision \citep{vrieze2012model}.
Traditionally, in the literature, model selection has been performed using
information criteria that measure the relative quality of a statistical
model, for a given set of data. The well-known penalized-based criteria can be divided into two frameworks. Under the frequentist framework, criteria includes Akaike information criterion (AIC) of \citet{Akaike1973}, Takeuchi
information criterion (TIC) of \citet{takeuchi1976distribution} and its variants, to name
just a few. Under the Bayesian framework, criteria includes deviance information
criterion (DIC) of \citet{spiegelhalter2002bayesian}, and its variants, including deviance information criterion under misspecified models ($\text{DIC}_M$) of \citet{li2020deviance}, to name
just a few. In fact, these popular information criterion have been established under the statistical decision theory discussed above. AIC, TIC, DIC and $\text{DIC}_M$ are all constructed by estimating the K-L divergence
between DGP and the corresponding predictive distributions. Hence, our developed approaches will follow the similar decision-theoretical framework. In the following, we first give two remarks to introduce these popular information criterion under this framework. After that, we propose our new information criterion in the following subsections.
\end{comment}
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 \citep{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 \citet{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.
\begin{remark}
Under some regularity conditions, under Bayesian framework, for misspecified models, \citet{li2020deviance} proposed the new version of DIC by \citet{spiegelhalter2002bayesian} named as so-called $
\text{DIC}_{M}^{k}$ for, that is, for model k,
\begin{equation}
\text{DIC}_{M}^{k}=-2\ln p(\mathbf{y}| \bar{\boldsymbol{\theta}} _{k},M_{k})+2P_{M}^{k},P_{M}^{k}=\mathbf{tr}\left\{ n\mathbf{\bar{\Omega}}
_{n}\left( \bar{\boldsymbol{\theta}}_{k} \right) V\left( \bar{\boldsymbol{\theta}}_{k} \right) \right\} , \label{dicm1}
\end{equation}
where $V\left( \bar{\boldsymbol{\theta}}_{k}\right) $ is the posterior
covariance matrix given by $V\left( \bar{\boldsymbol{\theta}}_{k}\right) =E\left[
\left( \boldsymbol{\theta}_{k}-\bar{\boldsymbol{\theta}}_{k}\right) \left( \boldsymbol{\theta}_{k}-\bar{\boldsymbol{\theta}}_{k}\right) ^{\prime }|\mathbf{y,}M_{k}
\right] $ and
$
\mathbf{\bar{\Omega}}_{n}\left( \widehat{\boldsymbol{\theta}}_{k}\right) =\frac{1}{n}
\sum_{t=1}^{n}\mathbf{s}_{t}\left( \widehat{\boldsymbol{\theta}}_{k}\right) \mathbf{
s}_{t}\left( \widehat{\boldsymbol{\theta}}_{k}\right) ^{\prime }$.
For this information criterion, \citet{li2020deviance} showed that the regular plug-in predictive distribution, $p\left( \mathbf{y}_{rep}|\mathbf{y},M_{k},d_{1}\right)
=p\left( \mathbf{y}_{rep}|\bar{\boldsymbol{\theta}}_{k},M_{k}\right) $ can be used for constructing the loss function and the corresponding risk function discussed above. Hence, from statistical decision viewpoint discussed above, when $a=1$,for misspecified models, it can be shown in \citet{li2020deviance} that
\begin{equation*}
Risk(d_{k^{1}})=E_{\mathbf{y}}\left( \mathcal{L}(\mathbf{y}
,d_{k^{1}})\right) =\int \mathcal{L}(\mathbf{y},d_{k^{1}})g(\mathbf{y})
\mbox{d}\mathbf{y}=E_{\mathbf{y}}\left[ \text{DIC}_{M}^{k}+2C\right] +o(1).
\end{equation*}
If the candidate models are restricted into correctly specified
models or good models which are good approximation to DGP, $\text{DIC}_{M}^{k}$ is reduced as a good approximation of $\text{DIC}^{k}$ of \citet{spiegelhalter2002bayesian} given by
\begin{equation}
\text{DIC}^{k}=-2\ln p(\mathbf{y}|\bar{\boldsymbol{\theta}}
_{k},M_{k})+2P_{D}^{k},P_{D}^{k}=\int 2\left[ \ln p(\mathbf{y}|\bar{\boldsymbol{\theta}}
_{k},M_{k})-\ln p(\mathbf{y}|\boldsymbol{\theta}
_{k},M_{k})\right]
\mbox{d}{\mbox{\boldmath${\theta}$}}. \label{dic2}
\end{equation}
It was shown in \citet{li2024deviance} that
\begin{equation*}
Risk(d_{k^{1}})=E_{\mathbf{y}}\left( \mathcal{L}(\mathbf{y}
,d_{k^{1}})\right) =\int \mathcal{L}(\mathbf{y},d_{k^{1}})g(\mathbf{y})
\mbox{d}\mathbf{y}=E_{\mathbf{y}}\left[ \text{DIC}^{k}+2C\right] +o(1),
\end{equation*}
More details about the theoretical development of $\text{DIC}^{k}$
and $\text{DIC}_{M}^{k}$, one can refer to \citet{spiegelhalter2002bayesian}, \citet{li2020deviance}, \citet{li2024deviance} and reference therein.
\end{remark}
\begin{remark}
For some misspecified model k, under frequentist framework, Takeuchi
information criterion (TIC) of \citet{takeuchi1976distribution} \footnote{TIC is originally developed by \citet{takeuchi1976distribution} for independent data and \citet{li2020deviance} relaxed this limitation to weakly dependent data} generally can be defined as
\begin{equation}
\text{TIC}^{k}=-2\ln p\left( \mathbf{y}|\widehat{\boldsymbol{\theta}}_{k}\right)
+2P_{T}^{k},P_{T}^{k}=-\mathbf{tr}\left\{ \mathbf{\bar{\Omega}}_{n}\left(
\widehat{\boldsymbol{\theta}}_{k}\right) \mathbf{\bar{H}}_{n}^{-1}\left(
\widehat{\boldsymbol{\theta}}_{k}\right) \right\} . \label{TICeq}
\end{equation}
From decision viewpoint, when $a=1$, the MLE, $\widehat{\boldsymbol{\theta}}_{k}$, replaced the
Bayesian estimator, ${\mbox{\boldmath${\bar{\theta}}$}}_{k}$ to formulate
the regular plug-in predictive distribution for constricting the risk function. Then, for misspecified models, it can be also shown in \citet{li2020deviance} that
\begin{equation*}
Risk(d_{k^{1}})=E_{\mathbf{y}}\left( \mathcal{L}(\mathbf{y}
,d_{k^{1}})\right) =\int \mathcal{L}(\mathbf{y},d_{k^{1}})g(\mathbf{y})
\mbox{d}\mathbf{y}=E_{\mathbf{y}}\left[ \text{TIC}^{k}+2C\right] +o(1).
\end{equation*}
Furthermore, when the candidate models are restricted into correctly
specified models or good models which are good approximation to DGP, TIC is reduced as the well-known AIC and it can
be shown in \citet{li2024deviance} that
\begin{equation*}
Risk(d_{k^{1}})=E_{\mathbf{y}}\left( \mathcal{L}(\mathbf{y}
,d_{k^{1}})\right) =\int \mathcal{L}(\mathbf{y},d_{k^{1}})g(\mathbf{y})
\mbox{d}\mathbf{y}=E_{\mathbf{y}}\left[ \text{AIC}^{k}+2C\right] +o(1),
\end{equation*}
where
\begin{equation}
\text{AIC}^{k}=-2\ln p(\mathbf{y}|\widehat{\boldsymbol{\theta}}_{k},M_{k})+2P^{k}
\label{aic2}
\end{equation}
\end{remark}
More details about the theoretical development of $\text{DIC}^{k}$
and $\text{DIC}_{M}^{k}$, one can refer to \citet{takeuchi1976distribution} , \citet{li2020deviance}, \citet{li2024deviance} and reference therein.
\subsection{Information Criterion for Comparing Misspecified Models based on
Variational Bayes Plug-in Predictive Distributions}
\label{subsec:plugin vbic}
Following the statistical decision theory shown in section \ref{subsec:decision theory}, 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 \citet{han2019statistical} and \citet{Zhang_2024}, we have
\begin{equation*}
{\mbox{\boldmath${\bar{\theta}}$}}_{k}^{VB}=\boldsymbol{\widehat{\theta}}
_{k}+O_{p}\left(n^{-3/4}\right),
\end{equation*}
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:
\begin{equation*}
\text{VDIC}_{M}^k=-2\ln p(\mathbf{y}|{\mbox{\boldmath${\bar{\theta}}$}}
_{k}^{VB},M_{k})+2 P^{k}_{VDIC_M},
\end{equation*}
where the penalty term $
P^k_{VDIC_M}$ for model $k$ is defined as
\begin{equation*}
\begin{aligned} P^{k}_{VDIC_M} = -\mathbf{tr}\left[
\mathbf{\bar{\Omega}}_{n}\left(
{\mbox{\boldmath${\bar{ \theta}}$}}_{k}^{VB}\right) \left(
\mathbf{\bar{H}}_{n}\left( {
\mbox{\boldmath${\bar{\theta}}$}}_{k}^{VB}\right) \right) ^{-1}\right].
\end{aligned}
\end{equation*}
\begin{theorem}
\label{theom vtic} Under Assumptions 1-8, we have,
\begin{equation*}
Risk(d_{k^1})=\int \text{VDIC}_M^{k}\times g(\mathbf{y})d\mathbf{y}+ 2C + o(1),
\text{ i.e., }E_{\mathbf{y}}(\text{VDIC}_M^{k})=Risk(d_{k^1}) -2C + o(1).
\end{equation*}
\end{theorem}
It can be proved that $\text{VDIC}_M^{k}$ is an asymptotically unbaised estimator of $Risk(d_{k^1})$ up to a constant.
\begin{remark}
For $\text{VDIC}_M$, $-2\ln p(\mathbf{y}|{\mbox{\boldmath${\bar{\theta}}$}}
_{k}^{VB},M_{k})$ can be understood as a Bayesian measure of fit, while $2
P^{k}_{VDIC_M}$ measures the model complexity. This feature of trade-off
between the goodness of fit of the model and the complexity of the model is
shared by other information criteria, such as \text{TIC} and $\text{DIC}_M$.
\end{remark}
\begin{remark}
Similar to \text{TIC} and $\text{DIC}_M$, $\text{VDIC}_M$ works for both correctly specified
and misspecified models.
\end{remark}
\subsection{Information Criterion for Comparing Misspecified Models based on
the VB Posterior Predictive Distribution}
\label{subsec:posterior vbic}
\begin{comment}
Again, following the statistical decision theory shown in section \ref{subsec:decision theory}, we may use $p^{VB}(\mathbf{y}_{rep}|\mathbf{y})$ for constructing the loss function and the corresponding risk function. Based on this posterior predictive distributions, a new information criteria can be developed to estimate $Risk(d_{k^2})$. Like other information criteria,
they can be used to select the optimal model. Let $\mathbf{\bar{\Omega}}
_{n}\left( {\mbox{\boldmath${\bar{\theta}}$}}_{k}^{VB}\right) $ , $\mathbf{
\bar{H}}_{n}\left( {\mbox{\boldmath${\bar{\theta}}$}}_{k}^{VB}\right) $ be
consistent estimators of $\mathbf{B}_{n}({\mbox{\boldmath${\theta}$}}
_{n}^{p})$ and $\mathbf{H}_{n}({\mbox{\boldmath${\theta}$}}_{n}^{p})$, and
the consistent estimator of $\mathbf{C}_{n}$ is $\mathbf{\hat{C}}_{n}\left( {
\mbox{\boldmath${\bar{\theta}}$}}_{k}^{VB}\right) =\left( \mathbf{\bar{H}}
_{n}\left( {\mbox{\boldmath${\bar{\theta}}$}}_{k}^{VB}\right) \right) ^{-1}
\mathbf{\bar{\Omega}}_{n}\left( {\mbox{\boldmath${\bar{\theta}}$}}
_{k}^{VB}\right) \left( \mathbf{\bar{H}}_{n}\left( {\mbox{\boldmath${\bar{
\theta}}$}}_{k}^{VB}\right) \right) ^{-1}$, where ${\mbox{\boldmath${\bar{
\theta}}$}}_{k}^{VB}$ is the mean of variational posterior density $
p^{VB}\left( \mbox{\boldmath${\theta}$}|\mathbf{y}\right) $.
When the misspecification is considered, we define the information criterion
using VB posterior predictive density as
\begin{equation*}
\text{VPIC}^{k}=-2\ln p(\mathbf{y}|{\mbox{\boldmath${\bar{\theta}}$}}
_{k}^{VB},M_{k})+2 P^{k}_{VPIC},
\end{equation*}
where it is named as Variational Predictive Information Criterion (VPIC) and the penalty term $
P^k_{VPIC}$ for model $k$ is defined as
\begin{equation*}
\begin{aligned} P^{k}_{VPIC} &= \frac{1}{2}\mathbf{tr}\left[
\mathbf{\bar{\Omega}}_{n}\left(
{\mbox{\boldmath${\bar{ \theta}}$}}_{k}^{VB}\right) \left(
-\mathbf{\bar{H}}_{n}\left( {
\mbox{\boldmath${\bar{\theta}}$}}_{k}^{VB}\right) \right) ^{-1}\right] +
\frac{1}{2}\ln \left( \left\vert \left( -\mathbf{\bar{H}}_{n}\left(
{\mbox{\boldmath${\bar{ \theta}}$}}_{k}^{VB}\right) \right) \left( -\mathbf{\bar{H}}_{n}^d\left(
{\mbox{\boldmath${\bar{ \theta}}$}}_{k}^{VB}\right) \right)^{-1} +\mathbf{I}_{n}\right\vert \right) \\
-&\frac{1}{2}\mathbf{tr}\left[\begin{array}{c}
\left(
-\mathbf{\bar{H}}_{n}\left(
{\mbox{\boldmath${ \bar{\theta}}$}}_{k}^{VB}\right) +\left( -\mathbf{\bar{H}}_{n}^d\left(
{\mbox{\boldmath${\bar{ \theta}}$}}_{k}^{VB}\right) \right) \right) ^{-1} \\
\times \left( \mathbf{\bar{\Omega}}_{n}\left(
{\mbox{\boldmath${\bar{\theta}}$} }_{k}^{VB}\right)+\left( -\mathbf{\bar{H}}_{n}^d\left(
{\mbox{\boldmath${\bar{ \theta}}$}}_{k}^{VB}\right) \right) \mathbf{\hat{C}}_{n}\left( {
\mbox{\boldmath${\bar{\theta}}$}}_{k}^{VB}\right) \left( -\mathbf{\bar{H}}_{n}^d\left(
{\mbox{\boldmath${\bar{ \theta}}$}}_{k}^{VB}\right) \right) \right)
\end{array}\right]\\
+&\frac{1}{2}\mathbf{tr}\left[ \left( -\mathbf{\bar{H}}_{n}^d\left(
{\mbox{\boldmath${\bar{ \theta}}$}}_{k}^{VB}\right) \right) \mathbf{\hat{C}}_{n}\left( {
\mbox{\boldmath${\bar{\theta}}$}}_{k}^{VB}\right) \right] . \end{aligned}
\end{equation*}
\end{comment}
Following the statistical decision theory outlined in Section \ref{subsec:decision theory}, 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}\)):
\begin{equation*}
\text{VPIC}^k = -2 \ln p(\mathbf{y}|{\boldsymbol{\bar{\theta}}}_k^{VB}, M_k) + 2 P_{VPIC}^k,
\end{equation*}
where the penalty term \(P_{VPIC}^k\) for model \(k\) is defined as:
\[
\begin{aligned}
P_{VPIC}^k = &\frac{1}{2} \mathbf{tr}\left[\mathbf{\bar{\Omega}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB}) \left(-\mathbf{\bar{H}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB})\right)^{-1}\right] \\
+ &\frac{1}{2} \ln \left\lvert \left(-\mathbf{\bar{H}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB})\right) \left(-\mathbf{\bar{H}}_{n}^d({\boldsymbol{\bar{\theta}}}_k^{VB})\right)^{-1} + \mathbf{I}_n \right\rvert \\
- &\frac{1}{2} \mathbf{tr}\Bigg[ \left(-\mathbf{\bar{H}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB}) + \left(-\mathbf{\bar{H}}_{n}^d({\boldsymbol{\bar{\theta}}}_k^{VB})\right)\right)^{-1} \\
&\hspace{1cm} \times \left( \mathbf{\bar{\Omega}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB}) + \left(-\mathbf{\bar{H}}_{n}^d({\boldsymbol{\bar{\theta}}}_k^{VB})\right) \mathbf{\hat{C}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB}) \left(-\mathbf{\bar{H}}_{n}^d({\boldsymbol{\bar{\theta}}}_k^{VB})\right) \right)\Bigg] \\
+ &\frac{1}{2} \mathbf{tr}\left[\left(-\mathbf{\bar{H}}_{n}^d({\boldsymbol{\bar{\theta}}}_k^{VB})\right) \mathbf{\hat{C}}_{n}({\boldsymbol{\bar{\theta}}}_k^{VB})\right].
\end{aligned}
\]
\begin{theorem}
\label{theom2} Under Assumptions 1-8, we have,
\begin{equation*}
Risk(d_{k^2})=\int \text{VPIC}^{k}\times g(\mathbf{y})d\mathbf{y} +2C +o(1),\text{
i.e., }E_{\mathbf{y}}(\text{VPIC}^{k})=Risk(d_{k^2})-2C +o(1)
\end{equation*}
\end{theorem}
It can be proved that $\text{VPIC}^{k}$ is an asymptotically unbaised estimator of $Risk(d_{k^2})$ up to a constant.
\begin{remark}
For $\text{VPIC}^{k}$, $-2\ln p(\mathbf{y}|{\mbox{\boldmath${\bar{\theta}}$}}
_{k}^{VB},M_{k})$ can be understood as a Bayesian measure of fit, while $2
P^{k}_{VPIC}$ measures the model complexity. $\text{VPIC}^{k}$ works for both correctly specified
and misspecified models.
\end{remark}
\subsection{BF and BIC type information criteria}
\label{subsec: bf and bic}
The BF and BIC belong to the first strand of model comparison in the section \ref{sec:intro}. 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, \citet{schwarz1978estimating} showed that the log-marginal likelihood can be approximated by
\begin{equation}
\label{bic_eq}
\begin{aligned}
\ln p\left(\mathbf{y} | M_k\right)=&\ln p\left(\mathbf{y} | \hat{\theta}_k, M_k\right)+\ln p\left(\hat{\boldsymbol{\theta}}_k | M_k\right)\\&+\frac{P_k \pi}{2}-\frac{P_k \ln n}{2}-\frac{\left|-\overline{\mathbf{H}}_n\left(\hat{\boldsymbol{\theta}}_k\right)\right|}{2}+O_p\left(\frac{1}{n}\right),
\end{aligned}
\end{equation}
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{bic_eq}) 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, \citet{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.
\begin{remark}
From the theoretical viewpoint, different criteria have different theoretical properties. BIC and BFs are consistent if the true model is one of the candidate models while $\mathrm{AIC}$, $\mathrm{TIC}$, $\mathrm{DIC}$, $\mathrm{DIC}_M$, $\mathrm{VDIC}_M$ and $\mathrm{VPIC}$ aim to provide the asymptotically unbiased estimator of the expected KL divergence between the DGP and a predictive distribution. When the true model is not included as a candidate model, which is often the case in practice, it is not clear what the best model selected by BIC and BFs can achieve. In this case, if one is concerned with the KL divergence between the DGP and a predictive distribution, it is expected that $\mathrm{TIC}$, $\mathrm{DIC}_M$ $\mathrm{VDIC}_M$ and $\mathrm{VPIC}$ perform better than BIC and BFs. Moreover, when the sample size is small, even when the true model is a candidate model, BIC and BFs may not select the true model. Again, if one is concerned with the KL divergence between the DGP and a predictive distribution, AIC and $\mathrm{DIC}$ can perform better than BIC and BFs. If one is considering model with massive data, in which MLE or MCMC methods can be intractable or costly, $\mathrm{VDIC}_M$ and $\mathrm{VPIC}$ will perform better.
\end{remark}
\section{Simulation and Empirical Studies}
\label{sec:simu and real}
\subsection{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 \citet{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
\begin{equation*}
\widehat{Risk\left(d_{k_*^1}\right)} = \frac{1}{1000} \sum_{i=1}^{1000} IC_{k^*}\left(
\mathbf{y}_i\right), \text{ for } \text{VDIC}_M^k,
\end{equation*}
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
\begin{equation*}
\widehat{Risk\left(d_{k_*^2}\right)} = \frac{1}{1000} \sum_{i=1}^{1000} IC_{k^*}\left(
\mathbf{y}_i\right), \text{ for } \text{VPIC}^k,
\end{equation*}
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 \citet{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
\begin{equation*}
\begin{aligned}
\widehat{Risk(d_{k^*})}_1 &= \frac{1}{1000}\sum_{i=1}^{1000} IC_{k^*}(\mathbf{y}_i), \text{ IC is } \text{VDIC}_M^k, \text{ and }\\
\widehat{Risk(d_{k^*})}_2 &= \frac{1}{1000}\sum_{i=1}^{1000} IC_{k^*}(\mathbf{y}_i), \text{ IC is } \text{VPIC}^k,
\end{aligned}
\end{equation*}
named as $\text{ELBO1}$, $\text{ELBO2}$, $\text{BIC1}$, and $\text{BIC2}$.
\subsubsection{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\).
\begin{figure}
\centering
\includegraphics[scale=0.4]{Figure_collection/linear_simu2.jpg}
\caption{The figure plots relative frequencies of the polynomial orders selected by different
criteria.}
\label{fig:poly_relative_k}
\end{figure}
\begin{table}[ht!]
\small
\centering
\caption{Averge $k^*$ selected under different criteria}
\begin{tabular}{cccccc}
\hline N & VPIC & $\text{VDIC}_M$ & ELBO & AIC & BIC \\
\hline 500 & 4.479 & 4.505 & 3.644 & 4.479 & 3.688 \\
5,000 & 5.792 & 5.798 & 5.546 & 5.799 & 5.166 \\
10,000 & 6.234 & 6.243 & 5.991 & 6.243 & 5.610 \\
100,000 & 7.906 & 7.883 & 8.191 & 7.887 & 6.989 \\
1,000,000 & 11.329 & 11.063 & 11.252 & 10.000 & 9.335 \\
\hline
\end{tabular}
\label{optimal k poly}
\end{table}
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{fig:poly_relative_k}. And the average values of $k^*$ is listed in the table \ref{optimal k poly}, all across 1,000 replications. Several interesting results can be found in Figure \ref{fig:poly_relative_k}. 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.
\begin{table}[ht!]
\small
\centering
\caption{Average risk of different criteria using polymomial regression (Scaled)}
\label{vbic for polymomial models}
\par
\begin{tabular}{cccccccc}
\hline & VPIC & $\text{VDIC}_M$ & ELBO1 & ELBO2 & AIC & BIC1 & BIC2 \\
\hline 500 & \textcolor{red}{1.42124} & 1.42177 & 1.42297 & 1.42355 & 1.42190 & 1.42268 & \textcolor{blue}{1.42326} \\
5000 & \textcolor{red}{14.19515} & 14.19565 & 14.19553 & 14.19605 & 14.19566 & 14.19710 & \textcolor{blue}{14.19764} \\
10000 & \textcolor{red}{28.38818} & 28.38866 & 28.38837 & 28.38886 & 28.38867 & 28.38987 & \textcolor{blue}{28.39038} \\
100000 & \textcolor{red}{283.82618} & 283.82660 & 283.82679 & 283.82715 & 283.82659 & 283.82735 & \textcolor{blue}{283.82779} \\
1000000 & \textcolor{red}{2837.99446} & 2838.00302 & 2838.01325 & \textcolor{blue}{2838.01344} & 2838.00639 & 2838.01196 & 2838.01227 \\
\hline
\end{tabular}
\end{table}
Table \ref{vbic for polymomial models} 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{vbic for polymomial models} 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.
\subsubsection{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
\begin{equation}
Y_i \mid X_i \overset{i . i . d .}{\sim} Bernoulli\left(\Phi
\left(X_i^{\prime} \beta \right)\right), \label{llk of probit}
\end{equation}
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,
\begin{table}[ht!]
\small
\centering
\caption{Candidate models of probit simulated data}
\label{candidate model of simulation2}
\begin{tabular}{cccc}
\hline Model & Numbers of variable & Model & Model specification\\
\hline $M_1$ & 1 & $Z_i=\left(1,x_{i1}\right)^{\prime}$ & Underfitting\\
$M_2$ & 1 & $Z_i=\left(1,x_{i2}\right)^{\prime}$ & Underfitting\\
$M_3$ & 1 & $Z_i=\left(1,x_{i3}\right)^{\prime}$ & Underfitting\\
$M_4$ & 2 & $Z_i=\left(1,x_{i1},x_{i2}\right)^{\prime}$ & Underfitting\\
$M_5$ & 2 & $Z_i=\left(1,x_{i1},x_{i3}\right)^{\prime}$ & Correctly specified\\
$M_6$ & 2 & $Z_i=\left(1,x_{i2},x_{i3}\right)^{\prime}$ & Underfitting\\
$M_7$ & 3 & $Z_i=\left(1,x_{i1},x_{i2},x_{i3}\right)^{\prime}$ & Overfitting\\
\hline
\end{tabular}
\end{table}
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.
\begin{table}[ht!]
\small
\centering
\caption{Average risk of different criteria using probit regression (Scaled)}
\label{vbic for probit simu}
\par
\begin{tabular}{cccccccc}
\hline & VPIC & $\text{VDIC}_m$ & ELBO1 & ELBO2 & AIC & BIC1 & BIC2 \\
\hline 500 & \textcolor{red}{0.59530} & 0.59624 & 0.59824 & \textcolor{blue}{0.59907} & 0.59626 & 0.59554 & 0.59643 \\
5,000 & \textcolor{red}{5.96466} & 5.96539 & 5.96474 & \textcolor{blue}{5.96562} & 5.96539 & 5.96473 & 5.96561 \\
10,000 & \textcolor{red}{11.92945} & 11.93039 & 11.92982 & \textcolor{blue}{11.93070} & 11.93040 & 11.92977 & 11.93065 \\
100,000 & \textcolor{red}{119.30140} & 119.30233 & 119.30170 & \textcolor{blue}{119.30258} & 119.30233 & 119.30169 & 119.30257 \\
1,000,000 & \textcolor{red}{1193.04493} & 1193.04586 & 1193.04522 & \textcolor{blue}{1193.04610} & 1193.04586 & 1193.04522 & \textcolor{blue}{1193.04610} \\
\hline
\end{tabular}
\end{table}
Table \ref{vbic for probit simu} 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.
\subsection{Empirical Studies}
\begin{comment}
In this subsection, we first use the linear model of different covariates to
investigate which model best predict the number of passengers transported in
the flight. In the second real study, we then investigate the credit risk
model, which is usually modelled by binary classification. In this study, we will use three new critera to choose which is the best model to explain the reason why the borrowers get loan from platform among six candidate probit models.
\end{comment}
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.
\subsubsection{US Domestic Flights Predictive Model}
\begin{comment}
In this section we use linear model with different covariates to investigate
which model best predict the number of passengers transported in the flight.
This dataset is related to the US domestic flights from 1990 to 2009, which
contains about $N=3.61$ Million data. Dataset we use is avalibale on Kaggle.
\citet{Chasiotis_Karlis_2024} use this data to fit
the linear regression model where $p=5$ measurements are selected. In this study, we can use the linear regression to model the relationships between
the dependent variable and covariates.
\end{comment}
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. \citet{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.
\begin{table}[ht!]
\small
\centering
\caption{Candidate model set for US domestic flights data}
\label{model descrpition in linear empirical}
\begin{tabular}{ccc}
\hline
Model & Description & Number of covariate\\ \hline
$M_1$ & $Y_i = \beta_0+\beta_1 x_1 + \beta_2 x_2 + \epsilon_i$ & 2\\
$M_2$ & $Y_i = \beta_0+\beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 +\epsilon_i$ & 3\\
$M_3$ & $Y_i = \beta_0+\beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4+ \beta_5 x_5 + \epsilon_i$ & 5\\
$M_4$ & $Y_i = \beta_0+\beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \dots + \beta_9 x_9 + \epsilon_i$ & 9\\ \hline
\end{tabular}
\end{table}
Table \ref{model descrpition in linear empirical} 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$.
\begin{comment}
Table \ref{Ic for air flight model} reports $\text{VPIC}^k$, $\text{VDIC}_M^k$
with corresponding penalty terms $P^k_{VPIC}$, $P^k_{VDIC_M}$,
and logrithm likelihood $\mathcal{L}(\boldsymbol{\theta},k)$ of candidate models
$\left\{M_k\right\}_{k=1}^4$. The difference of criteria is mainly caused by
the penalty term. First and most, both $\text{VPIC}^k$ and $\text{VDIC}_M^k$ choose the model $M_4$,
which means under the same criteria, we choose model $M_4$ as the better
model than other candidate models. What's more, among VPIC and $\text{VDIC}_M$, the $
\text{VPIC}^k$ of $M_4$ is the minimum. Thus we suggest choose the model $M_4$ and
use the VB posterior predictive distribution to make decision to reach the minimum predictive risk.
\end{comment}
Table \ref{Ic for air flight model} 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.
\begin{table}[ht!]
\small
\centering
\caption{Model selection results of 4 candidate models in US air flight data}
\label{Ic for air flight model}
\begin{tabular}{ccccc}
\hline & $\text{M}_1$ & $\text{M}_2$ & $\text{M}_3$ & $\text{M}_4$ \\
\hline $\text{VPIC}$ & 60255994.7307 & 60127455.0509 & 60123759.8081 & \textcolor{red}{60113021.1802} \\
$\text{VDIC}_M$ & 60256044.6951 & 60127505.9314 & 60123811.8674 & \textcolor{red}{60113073.3737} \\
$\text{ELBO}$ & -54214029.0097 & -54149777.3111 & -54147985.5373 & -54142673.4084 \\
$\text{AIC}$ & 60255940.7704 & 60127400.1087 & 60123703.9091 & 60112963.8209 \\
$\text{BIC}$ & 60255993.1500 & 60127465.5832 & 60123795.5734 & 60113107. 8649 \\
$\text{DIC}$ & 60255940.8067 & 60127400.1023 & 60123703.8001 & 60112963.8361 \\
$\text{DIC}_M$ & 60256042.8225 & 60127506.9269 & 60123810.6099 & 60113072.8852 \\
\hline
\end{tabular}
\end{table}
To show the difference among seven information criteria, we report a more detailed summary, shown in table \ref{detailed info of us air flight}. 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 \citet{li2024deviance} showed that $\text{DIC}$ is a Bayesian version of $\text{AIC}$.
\begin{table}[ht!]
\caption{Difference among fit term and penalty term}
\label{detailed info of us air flight}
\small
\centering
\begin{tabular}{ccccc}
\hline & $\text{M}_1$ & $\text{M}_2$ & $\text{M}_3$ & $\text{M}_4$ \\
\hline $\ell\left(\bar{\boldsymbol{\theta}}^{VB}\right)$ & -30127966.3852 & -30063695.0543 & -30061844.9545 & -30056470.9105 \\
$\ell\left(\bar{\boldsymbol{\theta}}\right)$ & -30127966.3856 & -30063695.0544 & -30061844.9552 & -30056470.9115 \\
$||\ell\left(\bar{\boldsymbol{\theta}}^{VB}\right)-\ell\left(\bar{\boldsymbol{\theta}}\right)||$ & 0.0003 & 0.0001 & 0.0007 & 0.0010 \\
$P^k_{\text{VPIC}}$ & 30.98 & 32.47 & 34.95 & 39.68 \\
$P^k_{\text{VDIC}_M}$ & 55.96 & 57.91 & 60.98 & 65.78 \\
$P^k_{\text{AIC}}$ & 4 & 5 & 7 & 11 \\
$P^k_{\text{BIC}}$ & 30.19 & 37.74 & 52.83 & 83.02 \\
$P^k_{\text{DIC}}$ & 4.02 & 5.00 & 6.94 & 11.01 \\
$P^k_{\text{DIC}_M}$ & 55.03 & 58.41 & 60.35 & 65.53 \\
\hline
\end{tabular}
\end{table}
\subsubsection{Credit Risk Analysis}
\begin{comment}
The credit risk model is the application of the Probit regression. It is
determined whether the loan should be given according to the borrower's
various information. When using binary classification, we will design $Y_i=1$
if borrower is allowed to give loans, and if $Y_i=0$, it means no loan. We
will use LendingClub dateset from Kaggle, ranging from 2007 to 2020 Q3.
\end{comment}
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{varaible descrpition} lists independent varaible and dependent variables that we are interested in.
\begin{table}[ht!]
\small
\centering
\caption{Variable description for credit risk model}
\label{varaible descrpition}
\begin{tabular}{ccc}
\hline
Variable & Symbol & Description \\ \hline
loan status & $Y_i$ & Current status of the loan, if loaned $Y_i=1$, else $
Y_i=0$ \\
annual inc & $AnuI_i$ & Annual income provided by the borrower during
registration \\
emp length & $Emp_i$ & Employment length in years. \\
dti & $DTI_i$ & Debt-to-Income Ratio, excluding mortgage and the requested
LC loan \\
loan amount & $Loanam_i$ & The amount of the loan applied for by the borrower
\\
term & $Term_i$ & The number of payments on the loan. \\ \hline
\end{tabular}
\end{table}
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{cred_risk_est}
\begin{table}[ht!]
\caption{Variational posterior mean and standard error of $\boldsymbol{\beta}$ in M1 and M2}
\label{cred_risk_est}
\centering
\begin{tabular}{cccccccc}
\hline & & $\beta_0$ & $\beta_1$ & $\beta_2$ & $\beta_3$ & $\beta_4$ & $\beta_5$ \\
\hline \multirow{2}{*}{M1} & $\mu_{VB}$ & 1.33 & 0.18 & 0.01 & -0.14 & -0.01 & -0.06 \\
& $\sigma^2_{VB}$ & 1.67E-02 & 1.62E-03 & 2.06E-04 & 1.29E-03 & 8.52E-05 & 1.74E-04 \\
\multirow{2}{*}{M2} & $\mu_{VB}$ & -1.18 & 0.49 & 0.01 & -0.06 & -0.02 & -0.11 \\
& $\sigma^2_{VB}$ & 7.69E-03 & 1.45E-03 & 1.92E-04 & 1.79E-03 & 2.24E-04 & 2.58E-04 \\
\hline
\end{tabular}
\end{table}
\begin{comment}
Table \ref{Ic for credit risk model} reports $\text{VPIC}^k$ and $\text{VDIC}_M^k$
with corresponding penalty terms $P^k_{VPIC}$, $P^k_{VDIC_M}$,
and logrithm likelihood $\mathcal{L}(\boldsymbol{\theta},k)$ of candidate models
$\left\{M_k\right\}_{k=1}^4$. The difference of criteria is mainly caused by
the penalty term. First and most, both $\text{VPIC}^k$ and $\text{VDIC}_M^k$ choose the model $M_4$,
which means under the same criteria, we choose model $M_4$ as the better
model than other candidate models. What's more, among VPIC and $\text{VDIC}_M$, the $
\text{VPIC}^k$ of $M_4$ is the minimum. Thus we suggest choose the model $M_4$ with the VB posterior predictive distribution to make decision to reach the minimum predictive risk.
\end{comment}
Table \ref{Ic for credit risk model} 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.
\begin{table}[ht!]
\small
\centering
\caption{Model selection results for the probit model and the logit model}
\label{Ic for credit risk model}
\begin{tabular}{cccccc}
\hline & $\text{VPIC}$ & $\text{VDIC}_M$ & $\text{ELBO}$ & $\text{AIC}$ & $\text{BIC}$ \\
\hline $\text{M}_1$ & \textcolor{red}{1572336.0802} & \textcolor{red}{1570922.7133} & -785536.9832 & 1570922.7620 & 1570996.9902 \\
$\text{M}_2$ & 1583798.1590 & 1582025.9001 & -825708.4686 & 1582024.7883 & 1582099.0165 \\
\hline
\end{tabular}
\end{table}
\section{Conclusion} \label{sec:conclude}
\begin{comment}
In the literature on machine learning research, VB technique has appealed a lot of attention for statistical inference of misspecified model with massive data. In this paper, firstly, for two possible predictive distributions based on variational posterior distribution, we investigate the risk function of these two predictive distributions which are
the expectation of the K-L divergence between the true density and the predictive distributions. Under some regularity conditions, we derive the corresponding asymptotically unbiased estimators for these two risk functions. After that, under the statistical decision theory, on the basis of risk functions, we propose two new information criterion for comparing misspecified models with massive data.
At last, we illustrate the proposed new information
criterions using some real studies in economics and finance under massive data.
\end{comment}
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.
\bibliographystyle{apalike}
\bibliography{ref2024}
\pagebreak