Posterior covariance information criterion for general loss functions

This subsection presents our methodologies. We work with the generalised posterior distribution given by

where $\pi(\theta)$ is a prior density on $\Theta$ and $s(x,\theta)$ is a score function that is the minus of a loss function. The score function may be different from the log-likelihood.

For the predictive evaluation of the generalised Bayesian approach, consider an evaluation function $\nu(x,\theta)$ for a future observation $x\in\mathcal{X}$ and a parameter vector $\theta\in\Theta$. Examples include the mean squared error $\nu(x,\theta)=\|x-\mathbb{E}[X\mid\theta]\|^{2}$, the $\ell_{1}$ error $\nu(x,\theta)=\|x-\mathbb{E}[X\mid\theta]\|$, the log likelihood $\nu(x,\theta)=\log p(x\mid\theta)$, the likelihood $\nu(x,\theta)=p(x\mid\theta)$, and the p-value $\nu(x,\theta)=\int_{x}^{\infty}p(x^{\prime}\mid\theta)dx^{\prime}$, where $\{p(x\mid\theta):\theta\in\Theta\}$ is a parametric model and $\mathbb{E}[\cdot\mid\theta]$ is the expectation with respect to $p(x\mid\theta)$.

On the basis of $\nu(x,\theta)$, we consider two types of predictive measures: the Gibbs generalisation error

and the plugin generalisation error

where $\mathbb{E}_{\mathrm{pos}}[\,\,\cdot\,\,\mid\,X^{n}]$ is the generalised posterior expectation, and $\mathbb{E}_{X_{n+1}}[\,\,\cdot\,\,]$ is the expectation with respect to $X_{n+1}$.

To estimate the Gibbs and plugin generalisation errors on the basis of current observations, we propose the Gibbs and the plugin posterior covariance information criteria $\mathrm{PCIC}_{\mathrm{G}}$ and $\mathrm{PCIC}_{\mathrm{P}}$:

where $\mathrm{Cov}_{\mathrm{pos}}[\,\,\cdot\,,\,\cdot\mid X^{n}]$ is the generalised posterior covariance. Table 1 summarises how we obtain estimates of the Gibbs and the plug-in generalisation errors. The first terms correspond to the empirical errors

Our proposed methods employ the posterior covariance as bias correction of empirical errors.

Before presenting theoretical results, we shall discuss the infinitesimal jackknife (IJ) interpretation of the proposed method. The IJ approach is a general methodology that approximates algorithms requiring the re-fitting of models, such as cross validation and the bootstrap methods. In the recent machine learning literature, this methodology has been rekindled as a linear approximation of LOOCV; see Beirami et al. (2017); Koh and Liang (2017); Giordano et al. (2019); Rad and Maleki (2020). We briefly explain about the IJ approximation of the leave-one-out $Z$-estimates given as

To describe the IJ approximation, we introduce the weighted $Z$-estimate

where the leave-one-out $Z$-estimate is denoted by

Then, by using the introduced weight $w$, the leave-one-out estimate $\widehat{\theta}^{(-i)}$ is linearly approximated as

which is known as the infinitesimal jackknife (IJ) approximation of leave-one-out estimates. By the implicit function theorem, the derivative respect to the weight is evaluated as

and then the IJ approximation becomes

What is the Bayesian version of this formula? Here we consider the Gibbs LOOCV

and the plug-in LOOCV

where $\mathbb{E}_{\mathrm{pos}}[\cdot\mid X^{-i}]$ is the expectation with respect to the leave-one-out generalised posterior distribution:

To investigate the IJ approximation of these cross validations, we define the weighted generalised posterior distribution

and denote by $\mathbb{E}_{\mathrm{pos}}^{w}[\,\,\cdot\,\,]$ the expectation with respect to $\pi_{w}(\theta\,;\,X^{n})$. Observe that we have

This, together with employing a structure analogous to equation (4), leads us to the following IJ approximations:

To evaluate the second terms in the IJ approximations, we employ the following variants of local sensitivity formula (Millar and Stewart, 2007; Giordano et al., 2018). For more details on Bayesian local sensitivity, refer to Thomas et al. (2018).

The proof is given in Appendix B. We should note that the final equation in the lemma is crucial for assessing the frequentist variance; see Giordano and Broderick (2023).

For the Gibbs LOOCV, Lemma 2.3 yields the form of the IJ approximation given by

For the plug-in LOOCV, further calculi are required. If $\nu(X_{i},\theta)$ is 2-times continuous differentiable with respect to $\theta$ and its Hessian with respect to $\theta$ is bounded, then the Taylor expansion around $\widehat{\theta}=\mathbb{E}_{\mathrm{pos}}[\theta\mid X^{n}]$ yields

where $\theta=(\theta^{1},\ldots,\theta^{d})$, and $\mathrm{rem}$ is negligible under additional assumptions. Here Lemma 2.3 gives

Observe that the chain rule gives

implying that

This concludes that the IJ approximation for the plug-in LOOCV is

where the remaining term $\mathrm{rem}$ is negligible.

So, the IJ approximations (8) and (9) of LOOCV presents $\mathrm{PCIC}_{\mathrm{G}}$ and $\mathrm{PCIC}_{\mathrm{P}}$. In the next subsection, we will show that the residual between the IJ appximations and the expected generalisation error is negligible.

In the literature on information criteria, the IJ methodology has been used to show the asymptotic equivalence between information criteria and LOOCV (Stone, 1974; Watanabe, 2010a); see also Konishi and Kitagawa (1996). The IJ approximation of the LOOCV $Z$-estimate requires the second order differentiation and its inverse calculation (c.f., Beirami et al., 2017; Koh and Liang, 2017). The discussion here emphasizes that in Bayesian framework, these calculi are avoidable and there are user-friendly surrogates, that is, the posterior covariances. Recently, Giordano and Broderick (2023) analyse the Bayesian infinitesimal jackknife approximation for estimating frequentist’s covariance in details, showing the consistency in finite dimensional models and giving discussions on the behaviours in high-dimensional models.

On the basis of the IJ approximation discussed in the previous subsection, this section presents the theoretical support for the use of $\mathrm{PCIC}_{\mathrm{G}}$. The same result holds for $\mathrm{PCIC}_{\mathrm{P}}$, and it is omitted.

The following is the main theorem stating the asymptotic unbiasedness of the proposed criteria. The proof consists of two steps: (1) the analysis on the residual between the IJ appximation and the LOOCV; (2) the analysis on the expected values of the difference between the generalisation error and the LOOCV. Suppose that current $n$ observations $X^{n}=(X_{1},\ldots,X_{n})$ are independent and identically distributed from an unknown sampling distribution on a sample space $\mathcal{X}$ in $\mathbb{R}^{d}$, and our prediction target $X_{n+1}$ follows the same sampling distribution as that of each observation and is independent of $X^{n}$. The expectation $\mathbb{E}[\cdot]$ denotes the expectation with respect to $(X_{1,}\ldots,X_{n+1})$.

We remark that the first part of the proof is essentially the same as the cumulant expansion used in the proof of asymptotic unbiasedness of WAIC (Watanabe, 2018).

Differential privacy is a cryptographic framework designed to ensure the privacy of individual users’ data while enabling meaningful statistical analyses with a desired level of efficiency (Dwork, 2006). It provides a formal guarantee that the inclusion or exclusion of a single individual’s data does not significantly affect the output of the analysis, thereby protecting sensitive information. Differential privacy has many real-world applications (for example, see Desfontaines and Pejó, 2020), and many strategies have been developed (Abadi et al., 2016; Wang et al., 2015a; Minami et al., 2016; Mironov, 2017; Jewson et al., 2023).

The foundational concept of ensuring differential privacy is the $(\varepsilon,\delta)$-differential private learner. For $\varepsilon,\delta\geq 0$, an $(\varepsilon,\delta)$-differentially private learner is defined as a randomized estimator $\widehat{\theta}$ that satisfies the following property: for any pair of adjacent datasets $X^{n},\widetilde{X}^{n}\in\mathcal{X}^{n}$, that is, datasets with

and for any measurable set $A$, the inequality

holds. Here, $\varepsilon$ represents the privacy budget, quantifying the level of privacy , while $\delta$ accounts for a small probability of privacy violation.

A widely employed strategy for achieving $(\varepsilon,\delta)$-differential privacy is the one-posterior-sample (OPS) estimator (Wang et al., 2015b; Minami et al., 2016). This estimator is a single sample from a generalized posterior distribution given by

where $L(\cdot,\cdot)$ denotes a user-specified loss or score function. The hyperparameter $\beta$ is intricately controlled by the privacy parameters $(\varepsilon,\delta)$, the choice of the score function $L$, and the dataset size $n$. This approach leverages the flexibility of posterior sampling to balance the trade-off between privacy guarantees and statistical utility.

Here, we demonstrate the application of $\mathrm{PCIC}_{\mathrm{G}}$ to understanding the predictive behaviour of OPS estimators with different values of $\beta$. We use two sets of classification data from UCI Machine Learning Repository (Dua and Graff, 2017), namely, the banknote authentication data set and the adult data set. The banknote authentication data set classifies genuine and forged banknote-like specimens based on four image features. The adult data set predicts whether income exceeds $50$K/yr based on 14 features from census data. We work with the generalised posterior distribution based on the logistic regression

where $\sigma(\cdot)$ is the sigmoid function: $\sigma(x):=1/(1+\exp(-x))$. We consider predictive evaluation of OPS estimators using three major classification losses, the Brier loss, the misclassification loss, and the spherical loss:

where for the $i$ observation, $p$ is given by $\sigma(X_{i}\theta)$ using $\theta$. The misclassification loss is non-differentiable with respect to $\theta$ and our theoretical results may not imply accurate generalisation error estimation; so, we confirm the applicability by numerical experiments. First, we randomly split the whole data into a training dataset of a sample size 50, 20 times. The test dataset is set to be the remaining. We then calculate the empirical errors using the training dataset, and the average of generalisation errors using the test dataset. We used 3980 MCMC samples after thinning out by 5 and a burnin period of length 100. The MCMC algorithm here employs the Gibbs sampler based on Polya-Gamma augmentation (Polson et al., 2013).

Figures 1 and 2 display the average values of generalisation error estimates relative to the average of the generalisation errors. The exact LOOCV is very close to the average of the generalisation error irrelevant to the loss and the dataset, but its computational cost is almost prohibitive. GDIC is not close to the average of the generalisation error. BPSIC works only for $\beta=1$ and a differentiable evaluation function. The proposed method $\mathrm{PCIC}_{\mathrm{G}}$ successfully estimates the generalisation errors of OPS estimators for all three evaluation functions, including the non-differentiable evaluation function. IS-CV and PSIS-CV are comparable to the proposed method; therefore, we cannot reach a conclusion regarding the performance difference of this setting.

We then apply the proposed methods to predictive evaluation of Bayesian hierarchical modeling. Here we shall work with a simple two-level normal model:

where $v_{j}^{2}$s are assumed to be given, $n_{j}$ is the sample size of the $j$-th group ($j=1,\ldots,J$), and $J$ is the number of groups. Here we assume that $n_{j}$ is equal to 1 and denote $X_{ij}$ by $X_{j}$, which is always possible by the sufficiency reduction. The prior distributions of $\mu$ and $\tau$ are specified by the scale mixture of the normal distribution with the half Cauchy distribution and by the half normal distribution, respectively:

where $\gamma$ and $\sigma_{\tau}$ are hyper parameters.

We measure the prediction accuracy of the information fusion in the Bayesian hierarchical structure across groups by (unscaled) $\ell_{2}^{2}$ loss

The working posterior in this application is rewritten as

where $\phi(x\,;\,a,b)$ is the normal density with mean $a$ and variance $b$. So, $\mathrm{PCIC}_{\mathrm{G}}$ becomes

To check the behaviour, we use two datasets commonly used in the demonstration of the Bayesian hierarchical modeling. One is the 8-schools data: the dataset consists of coaching effects with the standard deviations in $J=8$ different high schools in New Jersey; see Table 5.2 of Gelman et al. (2013). The other is the meta-analysis data from Yusuf et al. (1985): the dataset consists of $J=22$ clinical trials of $\beta$-blockers for reducing mortality after myocardial infarction; see Table 5.4 of Gelman et al. (2013). We apply empirical log-odds transformation so as to employ the normal model described above as in Section 5.8 of Gelman et al. (2013).

We design two types of experiments. One is changing the data scale by multiplying the data scaling factor $\delta_{X}\in\{1,2,\ldots,9,10\}$ to $X_{j}$s. The other is changing the prior scale by multiplying the prior scaling factor $\delta_{\gamma}\in\{0.1i\,:\,i=1,\ldots,10\}$ to $\gamma$. We fix $\sigma_{\tau}=10$ and set $\gamma=10$ as initial values. In the experiments, we first randomly split the original datasets of $J$ groups into the training data of the remaining $\lfloor J/2\rfloor$ groups and the test data of $J-\lfloor J/2\rfloor$ groups 20 times. We then obtain 5000 posterior samples using PyMC (Abril-Pla et al., 2023), calculate the empirical errors and the generalisation error estimates using the training data set, and calculate the average of generalisation errors using the test data set.

Figure 3 displays the generalisation error estimates relative to the average generalisation errors along with the change of the data scaling factor $\delta_{X}$. Figure 4 displays those along with the change of the prior scaling factor $\delta_{\gamma}$. As the data scaling factor becomes largers, all esitmates deteriorate (becomes far from the generalisation error) but the deteriorating rates of IS-CV and PCIC are slower than that of PSIS-CV. Fro any prior scaling factor, IS-CV and PCIC are closer to the generalisation error than PSIS-CV. GDIC is a relatively good metric for the meta-analysis data but not for the 8-school data.

We next apply the proposed method to Bayesian regression the presence of influential observations. The presence of influential observations impacts the variability of the case-deletion importance sampling weights, as discussed in Peruggia (1997), which results in the instability of IS-CV.

Here, we focus on the performance comparison between our method and IS-CV in Bayesian regression with influential observations. We work with the following Bayesian regression model in Peruggia (1997): Let $n$ be the sample size and let $X_{i}\in\mathbb{R}^{d}$ ($i=1,\ldots,n$) be the $d$ covariates. Then,

where $\mathrm{IG}$ denotes the inverse gamma distribution, $\mathrm{IW}$ denotes the inverse Wishart distribution, and $I_{d\times d}$ denotes the $d\times d$ identity matrix. Three quantities $\sigma_{a},\,\sigma_{b},\,\kappa$ are hyperparameters, and we fix $\sigma_{a}=\sigma_{b}=\kappa=1$. The score function $s(x,\theta)$ in this subsection is simply the log-likelihood. For the loss functions, we use $\ell_{2}$, $\ell_{1}$, scaled $\ell_{2}$ and scaled $\ell_{1}$ losses:

We begin with the following synthetic Bayesian regression model: for a given $R$, $\widetilde{X}_{1},\ldots,\widetilde{X}_{n}$ are fixed to

and then $Y_{1},\ldots,Y_{n}$ are given as

where $\varepsilon_{i}$s are i.i.d. from the standard Gaussian distribution, and $\beta_{0}$, $\beta_{1}$, and $\sigma$ are unknown parameters. We assign (0,1,1) to the true values of $(\beta_{0},\beta_{1},\sigma)$. We set 50 to the sample size $n$. For each $R$, we simulate the training data set 50 times and obtain 50 values of $\mathrm{PCIC}_{\mathrm{G}}$, PSIS-CV, and IS-CV. We then calculate the average of the generalisation errors using a test data set with a sample size of 10. By using the same Gibbs sampler as in Peruggia (1997), we obtained 2000 MCMC samples after thinning out by 10 and a burnin period of length 10000.

Figure 5 displays the comparison between $\mathrm{PCIC}_{\mathrm{G}}$ and IS-CV. Consider $\nu_{\ell_{2}}$. For $R\leq 2$, the performance of the two methods does not differ. For $R\geq 3$, the bias of IS-CV relative to the average of the generalisation errors, is larger than that of $\mathrm{PCIC}_{\mathrm{G}}$. Consider $\nu_{\mathrm{scaled}\,\ell_{1}}$. In this case, for all $R$, the bias of IS-CV, relative to the average of the generalisation errors, is larger than those of $\mathrm{PCIC}_{\mathrm{G}}$ and PSIS-CV. In particular, PSIS-CV performs the best for this loss.

We next employ two real datasets containing influential observations. One is the the stack loss data: The data consists of $n=21$ days of measurements with three condition records $x_{i}$. The other is the Gesell data: This consists of $n=21$ children’s records of the age $x_{i}$ at which they first spoke and their Gesell Adaptive Score $y_{i}$. These datasets are known to contain influential observations; the indices for influential observations are $i=21$ and $i=18$, respectively. In the experiments, we first randomly split the original datasets into the training data of $\lfloor n/2\rfloor$ samples and the test data of $n-\lfloor n/2\rfloor$ samples 20 times. We then obtain posterior samples with different sample sizes, calculate the empirical errors and the generalisation error estimates using the training data set, and calculate the average of generalisation errors using the test data set.