Using stacking to average Bayesian predictive distributions

TL;DR

This paper extends stacking methods to combine Bayesian predictive distributions, improving model averaging in complex settings using proper scoring rules, Pareto smoothed importance sampling, and regularization, with empirical validation showing its advantages.

Contribution

It generalizes stacking to predictive distributions with proper scoring rules, offering a more robust alternative to traditional Bayesian model averaging in the M-open setting.

Findings

Stacking of predictive distributions outperforms BMA and AIC-based methods.

Pareto smoothed importance sampling enables efficient leave-one-out computations.

BB-pseudo-BMA offers a computationally cheaper alternative with comparable performance.

Abstract

The widely recommended procedure of Bayesian model averaging is flawed in the M-open setting in which the true data-generating process is not one of the candidate models being fit. We take the idea of stacking from the point estimation literature and generalize to the combination of predictive distributions, extending the utility function to any proper scoring rule, using Pareto smoothed importance sampling to efficiently compute the required leave-one-out posterior distributions and regularization to get more stability. We compare stacking of predictive distributions to several alternatives: stacking of means, Bayesian model averaging (BMA), pseudo-BMA using AIC-type weighting, and a variant of pseudo-BMA that is stabilized using the Bayesian bootstrap. Based on simulations and real-data applications, we recommend stacking of predictive distributions, with BB-pseudo-BMA as an approximate alternative when computation cost is an issue.

Figures

Equations

$$p(\Delta |y)=\sum _{k=1}^{K}p(\Delta |M_k,y)p(M_k|y).$$

$$p(M_k|y)=\frac{p(y|M_k)p(M_k)}{{\sum }_{k=1}^{K}p(y|M_k)p(M_k)},$$

$$p(y|M_k)=\int p(y|{\theta }_k,M_k)p({\theta }_k|M_k)d{\theta }_k.$$

$$\hat{w}=\arg \underset{w}{\min }\sum _{i=1}^n{\left(y_i-\sum _k{w}_k{\hat{f}}_k^{(-i)}(x_i)\right)}^2.$$

$$\hat{\tilde{y}}=\sum _{k=1}^K{\hat{w}}_k{f}_k(\tilde{x}|{\hat{\theta }}_{k,y_{1:n}}).$$

$${\mathrm{E}}_{\tilde{y}}(u(\tilde{y},a(y,w))|y)=-\int ||\tilde{y}-\sum _{k=1}^K{w}_k{\hat{\tilde{y}}}_k|{|}^{2}p(\tilde{y}|y,M_t)d\tilde{y},$$

$$\hat{w}={\arg \max }_w\frac{1}{n}\sum _{i=1}^{n}u(y_i,a(y_{-i},w)),$$

$${\mathrm{E}}_k(\tilde{y}|{\hat{\theta }}_k,y)=\int \tilde{y}p(\tilde{y}|{\hat{\theta }}_{k,y},M_k)d\tilde{y},$$

$${\mathrm{E}}_k(\tilde{y}|y)=\int \tilde{y}p(\tilde{y}|{\theta }_k,M_k)p({\theta }_k|y,M_k)d\tilde{y}d{\theta }_k.$$

$$\int l(\tilde{y},a(y,w))p(\tilde{y}|y)d\tilde{y}-\frac{1}{n}\sum _{i=1}^{n}l(y_i,a(y_{-i},w))\text{xlongrightarrow}\text{L2}0,$$

$$U_n(w_k,{\theta }_k)=\sum _{i=1}^n\sum _{k=1}^K{w}_k{f}_k(y_i|{\theta }_k).$$

$$w_{k}=\frac{\exp\bigl{(}-n\mathrm{KL}(f^{*},f_{k})\bigr{)}}{\sum_{k=1}^{K}\exp\bigl{(}-n\mathrm{KL}(f^{*},f_{k})\bigr{)}}.$$

$$\mathrm{\widetilde{KL}_{1}}(f_{0},f_{k})=\int\!\!\int\mathrm{KL}\Bigl{(}f_{0}(\cdot|\theta_{0}),f_{k}(\cdot|\theta_{k})\Bigr{)}p(\theta_{k}|y,M_{k})p(\theta_{0}|y,M_{0})d\theta_{k}d\theta_{0},$$

$$\mathrm{\widetilde{KL}_{2}}(f_{0},f_{k})=\mathrm{KL}\Bigl{(}\int\!f_{0}(\cdot|\theta_{0})p(\theta_{0}|y,M_{0})d\theta_{0},\int\!f_{k}(\cdot|\theta_{k}))p(\theta_{k}|y,M_{k})d\theta_{k}\Bigr{)}.$$

$${\widetilde{\mathrm{KL}}}_2(f_0,f_k)=-\int \log p(\tilde{y}|y,M_k)p(\tilde{y}|y,M_0)d\tilde{y}+\int \log p(\tilde{y}|y,M_0)p(\tilde{y}|y,M_0)d\tilde{y}$$

$$w_{k}=\frac{\exp\bigl{(}n\int\!\log p(\tilde{y}|y,M_{k})p(\tilde{y}|y,M_{0})d\tilde{y}\bigr{)}}{\sum_{k=1}^{K}\exp\bigl{(}n\int\!\log p(\tilde{y}|y,M_{k})p(\tilde{y}|y,M_{0})d\tilde{y}\bigr{)}}$$

$$\min_{w}d\Bigl{(}\sum_{k=1}^{K}w_{k}p(\cdot|y,M_{k}),p_{t}(\cdot|y)\Bigr{)}.$$

$$\max_{w}S\Bigl{(}\sum_{k=1}^{K}w_{k}p(\cdot|y,M_{k}),p_{t}(\cdot|y)\Bigr{)},$$

$$p(\tilde{y}|y,M_k)=\int p(\tilde{y}|{\theta }_k,y,M_k)p({\theta }_k|y,M_k)d{\theta }_k.$$

$${\hat{p}}_{k,-i}(y_i)=\int p(y_i|{\theta }_k,M_k)p(\theta |y_{-i},M_k)d{\theta }_k.$$

$$\max_{w}\frac{1}{n}\sum_{i=1}^{n}S\Bigl{(}\sum_{k=1}^{K}w_{k}\hat{p}_{k,-i},y_{i}\Bigr{)},\quad s.t.\quad w_{k}\geq 0,\quad\sum_{k=1}^{K}w_{k}=1.$$

$$\hat{p}(\tilde{y}|y)=\sum _{k=1}^K{\hat{w}}_{k}p(\tilde{y}|y,M_k).$$

$$\underset{w}{\max }\frac{1}{n}\sum _{i=1}^n\log \sum _{k=1}^K{w}_{k}p(y_i|y_{-i},M_k),s.t.w_k\ge 0,\sum _{k=1}^K{w}_k=1.$$

$$\frac{1}{n}\sum_{i=1}^{n}S\Bigl{(}\sum_{k=1}^{K}w_{k}\hat{p}_{k,-i},y_{i}\Bigr{)}-\mathrm{E}_{\tilde{y}|y}S\Bigl{(}\sum_{k=1}^{K}w_{k}p(\tilde{y}|y,M_{k}),\tilde{y}\Bigr{)}\xlongrightarrow{\text{$L_{2}$}}0.$$

$$p(y_i|y_{-i},M_k)=\int p(y_i|{\theta }_k,M_k)p({\theta }_k|y_{-i},M_k)d{\theta }_k.$$

$$r_{i,k}^s=\frac{1}{p(y_i|{\theta }^s,M_k)}\propto \frac{p({\theta }^s|y_{-i},M_k)}{p({\theta }^s|y,M_k)}.$$

$$p(y_i|y_{-i},M_k)\approx \frac{1}{\frac{1}{S}{\sum }_{s=1}^S{w}_{i,k}^s}.$$

$${{\widehat{\mathrm{elpd}}}_{\mathrm{loo}}}^k=\sum _{i=1}^n\log \left(\frac{{\sum }_{s=1}^S{w}_{i,k}^{s}p(y_i|{\theta }^s,M_k)}{{\sum }_{s=1}^S{w}_{i,k}^s}\right).$$

$${\mathrm{elpd}}^k=\sum _{i=1}^n\int p_t({\tilde{y}}_i)\log p({\tilde{y}}_i|y,M_k)d{\tilde{y}}_i,$$

$${\widehat{\mathrm{elpd}}}_{\mathrm{loo}}^k=\sum _{i=1}^n\log p(y_i|y_{-i},M_k).$$

Full text

Bayesian model averaging,

model combination,

proper scoring rule,

predictive distribution,

stacking,

Stan,

\startlocaldefs\endlocaldefs

,

and

1 Introduction

A general challenge in statistics is prediction in the presence of multiple candidate models or learning algorithms $\mathcal{M}=(M_{1},\dots,M_{K})$. Model selection—picking one model that can give optimal performance for future data—can be unstable and wasteful of information (see, e.g., Piironen and Vehtari, 2017). An alternative is model averaging, which tries to find an optimal model combination in the space spanned by all individual models. In Bayesian context, the natural target for prediction is to find a predictive distribution that is close to the true data generating distribution (Gneiting and Raftery, 2007; Vehtari and Ojanen, 2012).

Ideally, we prefer to attack the Bayesian model combination problem via continuous model expansion—forming a larger bridging model that includes the separate models $M_{k}$ as special cases (Gelman, 2004)—but in practice constructing such an expansion can require conceptual and computational effort, and so it makes sense to consider simpler tools that work with existing inferences from separately-fit models.

1.1 Bayesian model averaging

If the set of candidate models between them represents a full generative model, then the Bayesian solution is to simply average the separate models, weighing each by its marginal posterior probability. This is called Bayesian model averaging (BMA) and is optimal if the prior is correct, that is, if the method is evaluated based on its frequency properties evaluated over the joint prior distribution of the models and their internal parameters (Madigan et al., 1996; Hoeting et al., 1999). If $y=(y_{1},\dots,y_{n})$ represents the observed data, then the posterior distribution for any quantity of interest $\Delta$ is,

[TABLE]

In this expression, each model is weighted by its posterior probability,

[TABLE]

and this expression depends crucially on the marginal likelihood under each model,

[TABLE]

In Bayesian model comparison, the relationship between the true data generator and the model list $\mathcal{M}=(M_{1},\dots,M_{K})$ can be classified into three categories: $\mathcal{M}$-closed, $\mathcal{M}$-complete and $\mathcal{M}$-open. We adopt the following definition from Bernardo and Smith (1994); Key et al. (1999); Clyde and Iversen (2013):

• •

$\mathcal{M}$-closed means the true data generating model is one of $M_{k}\in\mathcal{M}$, although it is unknown to researchers.

• •

$\mathcal{M}$-complete refers to the situation where the true model exists and is out of model list $\mathcal{M}$. But we still wish to use the models in $\mathcal{M}$ because of tractability of computations or communication of results, compared with the actual belief model. Thus, one simply finds the member in $\mathcal{M}$ that maximizes the expected utility (with respect to the true model).

• •

$\mathcal{M}$-open refers to the situation in which we know the true model $M_{t}$ is not in $\mathcal{M}$, but we cannot specify the explicit form $p(\tilde{y}|y)=p(\tilde{y}|M_{t},y)$ because it is too difficult, we lack time to do so, or do not have the expertise, computational intractability, etc.

BMA is appropriate for $\mathcal{M}$-closed case. In the $\mathcal{M}$-open and $\mathcal{M}$-complete case, BMA will asymptotically select the one single model in the list that is closest in Kullback-Leibler (KL) divergence.

Furthermore, in BMA, the marginal likelihood depends sensitively on the specified prior $p(\theta_{k}|M_{k})$ for each model. For example, consider a problem where a parameter has been assigned a normal prior distribution with center 0 and scale 10, and where its estimate is likely to be in the range $(-1,1)$. The chosen prior is then essentially flat, as would also be the case if the scale were increased to 100 or 1000. But such a change would divide the posterior probability of the model by roughly a factor of 10 or 100.

1.2 Predictive accuracy

From another direction, one can simply aim to minimize out-of-sample predictive error, equivalently to maximize expected predictive accuracy. In this paper we propose a novel log score stacking method for combining Bayesian predictive distributions. As a side result we also propose a simple model weighting using Bayesian leave-one-out cross-validation.

1.3 Stacking

Stacking is a direct approach for averaging point estimates from multiple models. The idea originates with Wolpert (1992), and Breiman (1996) gives more details for stacking weights under different conditions. In supervised learning where the data are $((x_{i},y_{i}),i=1,\dots,n)$, and each model $M_{k}$ has a parametric form $\hat{y}_{k}=f_{k}(x|\theta_{k})$, stacking is done in two steps (Ting and Witten, 1999). In the first, baseline-level, step, each model is fitted separately and the leave-one-out (LOO) predictor $\hat{f}_{k}^{(-i)}(x_{i})=E(y_{i}|\hat{\theta}_{k,y_{-i}},M_{k})$ is obtained for each model $k$ and each data point $i$. Cross-validation or bootstrapping is used to avoid overfitting (LeBlanc and Tibshirani, 1996). In the second, meta-level, step, a weight for each model is obtained by minimizing the mean squared error, treating the leave-one-out predictors from the previous stage as covariates:

[TABLE]

Breiman (1996) notes that a positive constraint $w_{k}\geq 0,k=1,\dots K$, or a $K-1$ simplex constraint: $w_{k}\geq 0,\sum_{k}^{K}w_{k}=1$ enforces a solution. Better predictions may be attainable using regularization (Merz and Pazzani, 1999; Yang and Dunson, 2014). Finally, the point prediction for a new data point with feature vector $\tilde{x}$ is

[TABLE]

It is not surprising that stacking typically outperforms BMA when the criterion is mean squared predictive error (Clarke, 2003), because BMA is not optimized to this task. Wong and Clarke (2004) emphasize that the BMA weights reflect the fit to the data rather than evaluating the prediction accuracy. On the other hand, stacking is not widely used in Bayesian model combination because it only works with point estimates, not the entire posterior distribution (Hoeting et al., 1999).

Clyde and Iversen (2013) give a Bayesian interpretation for stacking by considering model combination as a decision problem when the true model $M_{t}$ is not in the model list. If the decision is of the form $a(y,w)=\sum_{k=1}^{K}w_{k}\hat{y}_{k}$, then the expected utility under quadratic loss is,

[TABLE]

where $\hat{\tilde{y}}_{k}$ is the predictor of new data $\tilde{y}$ in model k. The stacking weights are the solution to the LOO estimator:

[TABLE]

where $a(y_{-i},w)=\sum_{k=1}^{K}w_{k}\mathrm{E}(y_{i}|y_{-i},M_{k})$.

The stacking predictor for new data $\tilde{y}$ is $\sum_{k=1}^{K}\hat{w}_{k}\hat{\tilde{y}}_{k}$. The predictor in the $k$-th model $\hat{\tilde{y}}_{k}$ can be either the plug-in estimator ,

[TABLE]

or the posterior mean,

[TABLE]

Le and Clarke (2017) prove that the stacking solution is asymptotically the Bayes solution. With some mild conditions on the distribution, the following asymptotic relation holds:

[TABLE]

where $l$ is the squared loss, $l(\tilde{y},a)=(\tilde{y}-a)^{2}.$ They also prove when the action is a predictive distribution $a(y_{-i},w)=\sum_{k=1}^{K}w_{k}p(y_{i}|y_{-i},M_{k})$, then the asymptotic relation still holds for negative logarithm scoring rules.

However, most early literature limited stacking to averaging point predictions, rather than predictive distributions. In this paper, we extend stacking from minimizing the squared error to maximizing scoring rules, hence make stacking applicable to combining a set of Bayesian posterior predictive distribution. We argue this is the appropriate version of Bayesian model averaging in the $\mathcal{M}$-open situation.

1.4 Fast leave-one-out cross-validation

Exact leave-one-out cross validation can be computationally costly. For example, in the econometric literature, Geweke and Amisano (2011, 2012) suggest averaging prediction models by maximizing predictive log score, only considering time series due to the computational cost of exact LOO for general data structures. In the present paper we demonstrate that Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO) (Vehtari et al., 2017) is a practically efficient way to calculate the needed leave-one-out predictive densities $p(y_{i}|y_{-i},M_{k})$ to compute log score stacking weights.

1.5 Akaike weights and pseudo Bayesian model averaging

Leave-one-out cross-validation is related to various information criteria (see, e.g. Vehtari and Ojanen, 2012). In case of maximum likelihood estimates, leave-one-out cross-validation is asymptotically equal to Akaike’s information criterion (AIC) Stone (1977). Given AIC =-2 log (maximum likelihood) + 2 (number of parameters), Akaike (1978) proposed to use $\exp(-\frac{1}{2}\mathrm{AIC})$ for model weighting (see also Burnham and Anderson, 2002; Wagenmakers and Farrell, 2004). More recently we have seen also Watanabe-Akaike information criterion (WAIC) (Watanabe, 2010) and leave-one-out cross-validation estimates used to compute model weights following the idea of AIC weights.

In Bayesian setting Geisser and Eddy (1979; see also, Gelfand 1996) proposed pseudo Bayes factors where marginal likelihoods $p(y|M_{k})$ are replaced with a product of Bayesian leave-one-out cross-validation predictive densities $\prod_{i=1}^{n}p(y_{i}|y_{-i},M_{k})$. Following the naming by Geisser and Eddy, we call AIC type weighting which uses Bayesian cross-validation predictive densities as pseudo Bayesian model averaging (Pseudo-BMA).

In this paper we show that the uncertainty in the future data distribution should be taken into account when computing such weights. We will propose a AIC type weighting using the Bayesian bootstrap and the expected log predictive density (elpd), which we call Pseudo-BMA+ weighting. We show that although Pseudo-BMA+ weighting gives better results than regular BMA or Pseudo-BMA weighting (in $\mathcal{M}$-open setting), it is still inferior to the log score stacking. Due to its simplicity we use Pseudo-BMA+ weights as initial guess for optimization procedure in the log score stacking.

1.6 Other model weighting approaches

Besides BMA, stacking and AIC type weighting, some other methods have been introduced to combine Bayesian models. Gutiérrez-Peña and Walker (2005) consider using a nonparametric prior in the decision problem stated above. Essentially they are fitting a mixture model with a Dirichlet process, yielding a posterior expected utility of,

[TABLE]

They then solve for the optimal weights $\hat{w}_{k}=\arg\max_{w_{k},\theta_{k}}U_{n}(w_{k},\theta_{k})$.

Li and Dunson (2016) propose model averaging using weights based on divergences from a reference model in $\mathcal{M}$-complete setting. If the true data generating density function is known to be $f^{*}$, then an AIC type (or Boltzmann-Gibbs type) weight can be defined as,

[TABLE]

The true model can be approximated with a reference model $M_{0}$ with density $f_{0}(.|\theta_{0})$ using nonparametric methods like Gaussian process or Dirichlet process, and $\mathrm{KL}(f^{*},f_{k})$ can be estimated by its posterior mean,

[TABLE]

or by the Kullback-Leibler divergence for posterior predictive distributions,

[TABLE]

Here, $\mathrm{\widetilde{KL}_{1}}$ corresponds to Gibbs utility, which can be criticized for not using the posterior predictive distributions (Vehtari and Ojanen, 2012), although asymptotically the two utilities are identical, and $\mathrm{\widetilde{KL}_{1}}$ is often computationally simpler than $\mathrm{\widetilde{KL}_{2}}$.

Let $p(\tilde{y}|y,M_{k})=\int\!f_{k}(\tilde{y}|\theta_{k})p(\theta_{k}|y,M_{k})d\theta_{k}$, $k=0,\dots,K$, then

[TABLE]

As the entropy of the reference model $\int\!\log p(\tilde{y}|y,M_{0})p(\tilde{y}|y,M_{0})d\tilde{y}$ is constant, the corresponding terms cancel out in the weight (1.2), leaving

[TABLE]

It is proportional to the exponential expected log predictive density, where the expectation is taken with respect to the reference model $M_{0}$. Comparing with definition 2.5 in Section 2.4, this method could be called Reference-Pseudo-BMA.

2 Theory and methods

We label classical stacking (1.1) as stacking of means because it combines the models by minimizing the mean squared error of the point estimate, which is the $L_{2}$ distance between the posterior mean and observed data. In general, we can use a proper scoring rule (or equivalently the underlying divergence) to compare distributions. After choosing that, stacking can be extended to combining the whole distribution.

2.1 Stacking using proper scoring rules

Adapting the notation of Gneiting and Raftery (2007), we label $Y$ as the random variable on the sample space $(\Omega,\mathcal{A})$ that can take values on $(-\infty,\infty)$. $\mathcal{P}$ is a convex class of probability measure on $\Omega$. Any member of $\mathcal{P}$ is called a probabilistic forecast. A function $S:\mathcal{P}\times\Omega\to\bar{\mathrm{R}}$ defines a scoring rule if $S(P,\cdot)$ is $\mathcal{P}$ quasi-integrable for all $P\in\mathcal{P}$. In the continuous case, distribution $P$ is identified with density function $p$.

For two probabilistic forecasts $P$ and $Q$, we write $S(P,Q)=\int S(P,\omega)dQ(\omega)$. A scoring rule $S$ is called proper if $S(Q,Q)\geq S(P,Q)$ and strictly proper if the equation holds only when $P=Q$ almost surely. A proper scoring rule defines the divergence $d:\mathcal{P}\times\mathcal{P}\to(0,\infty)$ as $d(P,Q)=S(Q,Q)-S(P,Q)$. For continuous variables, some popularly used scoring rules include:

Quadratic score: $\mathrm{QS}(p,y)=2p(y)-||p||_{2}^{2}$ with the divergence $d(p,q)=||p-q||_{2}^{2}$. 2. 2.

Logarithmic score: $\mathrm{LogS}(p,y)=\log(p(y))$ with $d(p,q)=\mathrm{KL}(q,p).$ The logarithmic score is the only proper local score assuming regularity conditions. 3. 3.

Continuous ranked score: $\mathrm{CRPS}(F,y)=-\int_{\mathrm{R}}(F(y^{\prime})-1(y^{\prime}\geq y))^{2}dy^{\prime}$ with $d(F,G)=\int_{\mathrm{R}}(F(y)-G(y))^{2}dy$, where $F$ and $G$ are the corresponding distribution function. 4. 4.

Energy score: $\mathrm{ES}(P,y)=\frac{1}{2}\mathrm{E}_{P}||Y-Y^{\prime}||_{2}^{\beta}-\mathrm{E}_{P}||Y-y||_{2}^{\beta}$, where $Y$ and $Y^{\prime}$ are two independent random variables from distribution $P$. When $\beta=2$, this becomes $\mathrm{ES}(P,y)=-||\mathrm{E}_{P}(Y)-y||^{2}.$ The energy score is strictly proper when $\beta\in(0,2)$ but not when $\beta=2$. 5. 5.

Scoring rules depending on first and second moments: Examples include $S(P,y)=-\log\mathrm{det}(\Sigma_{P})-(y-\mu_{P})^{T}\Sigma_{p}^{-1}(y-\mu_{P})$, where $\mu_{P}$ and $\Sigma_{P}$ are the mean vector and covariance matrix of distribution $P$.

Now return to the problem of model combination after specifying the score rule $S$ and corresponding divergence $d$. The observed data are $y=(y_{1},\dots,y_{n})$. For simplicity, we remove all covariates $x$ in the notation. Suppose we have a set of probabilistic models $\mathcal{M}=(M_{1},\dots,M_{K})$; then the goal in stacking is to find an optimal super-model in the convex linear combination with the form $\mathcal{C}=(\sum_{k=1}^{K}w_{k}p(\cdot|M_{k})|\sum_{k}w_{k}=1,w_{k}\geq 0)$, such that its divergence to the true data generating model, denoted by $p_{t}(\cdot|y)$, is minimized:

[TABLE]

Or equivalently maximize the scoring rule of the predictive distribution,

[TABLE]

where $p(\tilde{y}|y,M_{k})$ is the predictive density of $\tilde{y}$ in model $k$:

[TABLE]

We label the leave-$i$-out predictive density in model $k$ as,

[TABLE]

where $y_{-i}=(y_{1},\dots,y_{i-1},y_{i+1},\dots,y_{n})$.

Then we define the stacking weights as the solution to the following optimization problem:

[TABLE]

Eventually, the combined estimation of the predictive density is

[TABLE]

When using logarithmic score (corresponding to Kullback-Leibler divergence), we call this stacking of predictive distributions:

[TABLE]

The choice of scoring rule can depend on the underlying application. Stacking of means (1.1) corresponds to the energy score with $\beta=2$. The reasons why we prefer stacking of predictive distributions (corresponding to the logarithmic score) to stacking of means are: (i) the energy score with $\beta=2$ is not a strictly proper scoring rule and can give rise to identification problems, and (ii) every proper local scoring rule is equivalent to the logarithmic score (Gneiting and Raftery, 2007).

2.2 Asymptotic behavior of stacking

The stacking estimate (2.1) finds the optimal predictive distribution within the convex set \mathcal{C}=\big{\{}\sum_{k=1}^{K}w_{k}p(\cdot|M_{k})\mid\sum_{k=1}^{K}w_{k}=1,w_{k}\geq 0\big{\}}, that is the closest to the data generating process with respect to the chosen scoring rule. This is different from Bayesian model averaging, which asymptotically with probability 1 will select a single model: the one that is closest in KL divergence to the true data generating process.

Solving for the stacking weights in (2.2) is an M-estimation problem. Under some mild conditions (Le and Clarke, 2017; Clyde and Iversen, 2013; Key et al., 1999), for either the logarithmic scoring rule or energy score (negative squared error) and a given set of weights $w_{1}\dots w_{k}$, as sample size $n\to\infty$, the following asymptotic limit holds:

[TABLE]

Thus the leave-one-out-score is a consistent estimator of the posterior score. In this sense, the stacking weights is the optimal combination weights asymptotically.

In terms of Vehtari and Ojanen (2012, Section 3.3), the proposed stacking with log score is $M_{*}$-optimal projection of the information in the actual belief model $M_{*}$ to $\hat{w}$, where explicit specification of $M_{*}$ is avoided by re-using data as a proxy for the predictive distribution of the actual belief model and $w_{k}$ are the free parameters.

2.3 Pareto smoothed importance sampling

One challenge in calculating the stacking weights proposed in (2.2) is that we need the leave-one-out (LOO) predictive density,

[TABLE]

Exact LOO requires refitting each model $n$ times. To avoid this onerous computation, we use the following approximate method. For the $k$-th model, we fit to all the data, obtaining simulation draws $\theta_{k}^{s}(s=1,\dots S)$ from the full posterior $p(\theta_{k}|y,M_{k})$ and calculate

[TABLE]

The ratio $r_{i,k}^{s}$ has a density function in its denominator and can be unstable, due to a potentially long right tail. This problem can be resolved using Pareto smoothed importance sampling (PSIS). For each fixed model $k$ and data $i$, we fit the generalized Pareto distribution to the 20% largest importance ratios $r_{i,k}^{s}$, and calculate the expected values of the order statistics of the fitted generalized Pareto distribution. We further truncate those values to get the smoothed importance weight $w_{i,k}^{s}$, which is used to replace $r_{i,k}^{s}$. For details of PSIS, see Vehtari et al. (2017). In the end, the LOO importance sampling is performed using

[TABLE]

When stacking using the logarithmic score, we are combining each model’s log predictive density. The PSIS estimate of the LOO expected log pointwise predictive density in the $k$-th model is,

[TABLE]

The reliability of the PSIS approximation can be determined by the estimated shape parameter $\hat{k}$ in the generalized Pareto distribution. For the left-out data points where $\hat{k}>0.7$, Vehtari et al. (2017) suggest replacing the PSIS approximation of those problematic cases by the exact LOO or $k$-fold validation.

One potential drawback of LOO is the large variance when the sample size is small. We see in the simulation that when the ratio of relative sample size to effective number of parameters is small, the weighting can be unstable. How to adjust this small sample behavior is left for the future research.

2.4 Pseudo-BMA

In our paper, we also consider an AIC type weighting using leave-one-out cross-validation estimated expected log predictive densities. As mentioned in Section 1.6, these weights estimate the same quantities as Li and Dunson (2016) that use the divergence from the reference model based inference.

To maintain comparability with the given dataset and to get easier interpretation of the differences in scale of effective number of parameters, we define the expected log pointwise predictive density (elpd) for a new dataset $\tilde{y}$ as a measure of predictive accuracy of a given model for the $n$ data points taken one at a time (Gelman et al., 2014). In model $M_{k}$,

[TABLE]

where $p_{t}(\tilde{y}_{i})$ denotes the true distribution of future data $\tilde{y}_{i}$.

Given observed data y, we estimate elpd using LOO:

[TABLE]

Then we define the Pseudo-BMA weighting for model $k$:

[TABLE]

However, this estimation doesn’t take into account the uncertainty resulted from having a finite number of proxy samples from the future data distribution. Taking into account the uncertainty would regularize the weights making them go further away from 0 and 1.

The computed estimate elpd is defined as the sum of $n$ independent components so it is trivial to compute their standard errors by computing the standard deviation of the $n$ pointwise values (Vehtari and Lampinen, 2002). Define

[TABLE]

and then we can calculate

[TABLE]

Simple modification of weights is to use the lognormal approximation:

[TABLE]

Finally, Bayesian bootstrap (BB) can be used to compute uncertainties related to LOO estimation (Vehtari and Lampinen, 2002). Bayesian bootstrap (Rubin, 1981) makes simple non-parametric approximation to the distribution of random variable. Having samples of $z_{1},\dots,z_{n}$ from a random variable $Z$, it is assumed that posterior probabilities for all observed $z_{i}$ have the distribution $\mbox{Dirichlet}(1,\dots,1)$ and values of $Z$ that are not observed have zero posterior probability. Sampling from the uniform Dirichlet distribution gives BB samples from the distribution of $Z$ and thus samples of any parameter of this distribution can be obtained. In other words, each BB replication generates a set of posterior probability $\alpha_{1:n}$ for all observed $z_{1:n}$,

[TABLE]

This leads to one BB replication of any statistic $\phi(Z)$ that is of interest:

[TABLE]

The distribution over all replicated $\hat{\phi}(Z|\alpha)$ (i.e., generated by repeated sampling of $\alpha$) produces an estimation for $\phi(Z)$.

As the distribution of $\mathrm{\widehat{elpd}}^{k}_{\mathrm{loo},i}$ is often highly skewed, BB is likely to work better than the Gaussian approximation. In our model weighting, we can define

[TABLE]

We sample vectors $(\alpha_{1,b},\dots,\alpha_{n,b})_{b=1,\dots,B}$ from the Dirichlet $(\overbrace{1,\dots,1}^{n})$ distribution, and compute the weighted means,

[TABLE]

Then a Bayesian bootstrap sample of $w_{k}$ with size $B$ is,

[TABLE]

and the final adjusted weight of model $k$ is,

[TABLE]

which we call Pseudo-BMA+ weights.

3 Simulation examples

In this section, we first illustrate the advantage of stacking of predictive distributions with a Gaussian mixture model. Then we compare stacking, BMA, Pseudo-BMA, Pseudo-BMA+ and other averaging methods through a series of linear regression simulation, where stacking gives the best performance in most cases. Finally we apply stacking to two real datasets, averaging multiple models so as to better explain US Senate voting or well-switching pattern in Bangladesh.

3.1 Gaussian mixture model

This simple example helps us understand how BMA and stacking behave differently. It also illustrates the importance of the choice of the scoring rules in combining distributions. Suppose the observed data $y=(y_{i},i=1,\dots,n)$ come iid from a normal distribution N$(3.4,1)$, not known to the data analyst, and there are 8 candidate models, $\mbox{N}(1,1)$, $\mbox{N}(2,1)$,…, $\mbox{N}(8,1)$. This is an $\mathcal{M}$-open problem in that none of the candidates is the true model, and we have set the parameters so that the models are somewhat separate but not completely distinct in their predictive distributions.

For BMA with a uniform prior $\mbox{Pr}(M_{k})=\frac{1}{8},k=1,\dots,8$, we can write the posterior distribution explicitly:

[TABLE]

from which we see that \hat{w}_{3}^{\mathrm{BMA}}\xlongrightarrow{\text{P}}1 and \hat{w}_{k}^{\mathrm{BMA}}\xlongrightarrow{\text{P}}0 for $k\neq 3$ as sample size $n\to\infty$. Furthermore, for any given $n$ ,

[TABLE]

where $\mu=3.4$ and $\mu_{k}=k$ in this setting.

This example is simple in that there is no parameter to estimate within each of the models: $p(\tilde{y}|y,M_{k})=p(\tilde{y}|M_{k})$. Hence, in this case the weights from Pseudo-BMA and Pseudo-BMA+ are the same as the BMA weights, $\exp(\mathrm{\widehat{elpd}}_{\mathrm{loo}}^{k})/\sum_{k^{\prime}}\exp(\mathrm{\widehat{elpd}}_{\mathrm{{loo}}}^{k^{\prime}})$.

For stacking of means, we need to solve

[TABLE]

This is nonidentifiable because the solution contains any vector $\hat{w}$ satisfying

[TABLE]

For point prediction, the stacked prediction is always $\sum_{k=1}^{8}\hat{w}_{k}k=\frac{1}{n}\sum_{i=1}^{n}y_{i}$, but it can lead to different predictive distributions $\sum_{i=1}^{k}\hat{w}_{i}\mbox{N}(k,1)$. To get one reasonable result, we transform the least square optimization to the following normal model and assign a uniform prior to $w$:

[TABLE]

Then we could use the posterior means of $w$ as model weights.

For stacking of predictive distributions,we need to solve

[TABLE]

In fact, this example is a density estimation problem. Smyth and Wolpert (1998) first apply stacking to non-parametric density estimation, which they call stacked density estimation and now can be viewed as a special case of our stacking method.

We compare the posterior predictive distribution $\hat{p}(\tilde{y}|y)=\sum_{k}\hat{w}_{k}p(\tilde{y}|y,M_{k})$, for these three methods of model averaging. Figure 1 shows the predictive distributions in one simulation when the sample size $n$ varies from 3 to 200. We first notice that stacking of means (the middle row of graphs) gives an unappealing predictive distribution, even if its point estimate is reasonable. The broad and oddly spaced distribution here arises from nonidentification of $w$, and it demonstrates the general point that stacking of means does not even try to match the shape of the predictive distribution. The top and bottom row of graphs show how BMA picks up the single model that is closest in KL divergence, while stacking picks a combination; the benefits of stacking becomes clear for large $n$.

In this trivial non-parametric case, stacking of predictive distributions is almost the same as fitting a mixture model, except for the absence of the prior. The true model N$(3.4,1)$ is actually a convolution of single models rather than a mixture, hence no approach can recover the true one from the model list. From Figure 2 we can compare the mean squared error and the mean logarithmic score of this three combination methods. The average log scores and errors are calculated through 500 repeated simulation and 200 test data generating from the true distribution. The left panel shows logarithmic score (or equivalent, expected log predictive density) of the predictive distribution. Stacking of predictive distributions always gives a larger score except for extremely small $n$. In the middle panel, it shows the main squared error by considering the posterior mean of predictive distribution to be a point estimation, even if it is not our focus. In this case, it is not surprising to see that stacking of predictive distributions gives almost the same optimal mean squared error as the stacking of means, both of which are better than the BMA. Two distributions close in KL divergence are close in each moment, while the reverse doesn’t necessarily hold. This illustrates the necessary of matching the distribution, rather than matching the moments for the predictive distribution.

Finally it is worth pointing out that stacking depends only on the space expanded by all candidate models, while BMA or Pseudo-BMA weighting may by misled by such model expansion. If we add another N$(4,1)$ as the $9$th model in the model list above, stacking will not change at all in theory, even though it becomes non-strictly-convex and has infinite same-height mode. For BMA, it is equivalent to putting double prior mass on the original $4$th model, which doubles the final weights for it. The right panel of Figure 2 shows such phenomenon: we fix sample size $n$ to be 15 and add more and more N$(4,1)$ models. As a result, BMA (or Pseudo-BMA weighting) puts more weights on N$(4,1)$ and behaves worse, while stacking changes nowhere except for numerical fluctuation. It illustrates another benefit of stacking compared to BMA or Pseudo-BMA weights. We would not expect a combination method that would performs even worse as model candidate list expands, which may become a disaster when many similar weak models exist. We are not saying BMA can never work in this case. In fact some other methods are proposed to make BMA overcome such drawbacks. For example, George (2010) establishes dilution priors to compensate for model space redundancy for linear models, putting less weights on those models that are close to each other. Fokoue and Clarke (2011) introduce prequential model list selection to obtain an optimal model space. But we propose stacking as a more straightforward solution.

3.2 Linear subset regressions

The previous section demonstrates a simple example of combining several different nonparametric models. Now we turn to the parametric case. This example comes from Breiman (1996) who compares stacking to model selection. Here we work in a Bayesian framework.

Suppose the true model is

[TABLE]

In the model $\epsilon$ independently comes from N$(0,1)$. All the covariates $X_{j}$ are independently from N$(5,1)$. The number of predictors $J$ is 15. The coefficient $\beta$ is generated by

[TABLE]

where $\gamma$ is determined by fixing the signal-to-noise ratio such that

[TABLE]

The value $h$ determines the number of nonzero coefficients in the true model. For $h=1$, there are 3 “strong” coefficients. For $h=5$, there are 15 “weak” coefficients. In the following simulation, we fix $h=5$.

We consider the following two cases:

$\mathcal{M}$-open

Each subset contains only one single variable. Hence, the $k$-th model is a univariate linear regression with the $k$-th variable $X_{k}$. We have $K=J=15$ different models in total. One advantage of stacking and Pseudo-BMA weighting is that they are not sensitive to prior, hence even a flat prior will work, while BMA can be sensitive to the prior. For each single model $M_{k}:Y\sim\mbox{N}(\beta_{k}X_{k},\sigma^{2})$, we set prior $\beta_{k}\sim N(0,10)$, $\sigma\sim\mathrm{Gamma}(0.1,0.1)$. 2. 2.

$\mathcal{M}$-closed

Let model $k$ be the linear regression with subset $(X_{1},\dots,X_{k})$. Then there are still $K=15$ different models. Similarly, in model $M_{k}:Y\sim\mbox{N}(\sum_{j=1}^{k}\beta_{j}X_{j},\sigma^{2})$, we set prior $\beta_{j}\sim N(0,10),j=1,\dots,k$, $\sigma\sim\mathrm{Gamma}(0.1,0.1)$.

In both cases, we have seven methods for combine predictive densities: (1) stacking of predictive distributions, (2) stacking of means, (3) Pseudo-BMA, (4) Pseudo-BMA+, (5) best model selection by mean LOO value, (6) best model selection by marginal likelihood, and (7) BMA. A linear combination $p(\tilde{y}|y)=\sum_{k=1}^{K}\hat{w}_{k}p(\tilde{y}|M_{k})$ is what we have in the end as estimation for the posterior density of the new data $\tilde{y}$. We generate a test dataset $(\tilde{x}_{i},\tilde{y}_{i})$, $i=1,\dots,200$ from the underlying true distribution to calculate the out of sample score for the combined distribution under each method $k$:

[TABLE]

We loop the test simulation 100 times to get the expected predictive accuracy for each method.

Figure 3 shows the expected out-of-sample log predictive density for the seven methods, for a set of experiments with sample size $n$ ranging from 5 to 200. Stacking seems to outperform all other methods even for small $n$. Stacking of predictive distributions is asymptotically better than any other combination method. Pseudo-BMA+ weighting dominates naive Pseudo-BMA weighting. Finally, BMA performs similarly to Pseudo-BMA weighting, always better than any kind of model selection, but that advantage vanishes in the limit since BMA picks up one model. In this $\mathcal{M}$-open setting, we know model selection can never be optimal.

The results change when we move to the second case, in which the $k$-th model contains variables $X_{1},\dots,X_{k}$ so that we are comparing models of differing dimensionality. The problem is $\mathcal{M}$-closed because the largest subset contains all the variables, and we have simulated data from this model. Figure 4 shows the mean log predictive density of the seven combination methods in this case. For a large sample size $n$, almost all methods recover the true answer (putting weight 1 on the full model), except BMA and model selection based on marginal likelihood. The poor performance of BMA comes from the parameter priors: recall that the optimality of BMA arises when averaging over the priors and not necessarily conditional on any particular chosen set of parameter values. There is no general no way to obtain a “correct” prior that accounts for the complexity for BMA in an arbitrary model space. Model selection by LOO can recover the true model, while selection by marginal likelihood cannot due to the same prior problems. Once again, we see that BMA eventually become the same as model selection by marginal likelihood, which is much worse than any other methods asymptotically.

In this example, stacking is unstable for extremely small $n$. In fact, our computations for stacking of predictive distributions and Pseudo-BMA depend on the the PSIS approximation $\log p(y_{i}|y_{-i})$. If this approximation is crude, then the second step optimization cannot be accurate. It is known that the parameter $\hat{k}$ in the generalized Pareto distribution can be used to diagnose the accuracy of PSIS approximation. In our method, we replace PSIS approximation by running exact LOO for any data points with estimated $\hat{k}>0.7$ (see Vehtari et al., 2017).

Figure 5 shows the comparison of the mean elpd estimated by LOO,

[TABLE]

and the mean elpd calculated using 200 independent test data,

[TABLE]

for each model $k$ and each sample size in the simulation described above. The area of each dot in Figure 5 represents the relative complexity of the model as measured by effective number of parameters in the model divided by sample size. We evaluate the effective number of parameters using LOO (Vehtari et al., 2017). The sample size $n$ varies from 30 to 200 and variable size is fixed to be 20. Clearly, the relationship is far from the line $y=x$ for extremely small sample size, and the relative bias ratio ($\mathrm{elpd_{loo}}/\mathrm{elpd_{test}}$) depends on the complexity of the model. Empirically, we have found the approximation to be poor when the sample size is less than 5 times the number of parameters.

As a result, in the small sample case, LOO can lead to relatively large variance, which makes the stacking of predictive distributions and Pseudo-BMA/ Pseudo-BMA+ unstable, with performance improving quickly as $n$ grows.

3.3 Comparison with mixture models

Stacking is inherently a two-step procedure. In contrast, when fitting a mixture model, one estimates the model weights and the status within parameters in the same step. In a mixture model, given a model list $\mathcal{M}=\{M_{1},\dots,M_{k}\}$, each component in the mixture occurs with probability $w_{k}$. Marginalizing out the discrete assignments yields the joint likelihood

[TABLE]

The mixture model seems to be the most straightforward continuous model expansion. Nevertheless, there are several reasons why we may prefer stacking to the mixture model, though the latter one is a full Bayesian approach. First, the computation cost of mixture models can be relatively large. If the true model is a mixture model and the estimation of each model depends a lot on others, then it is worth paying the extra computation cost. However, it is not quite possible to combine several components in real application. It is more likely that the researcher is running combination among several mixture models with different pre-specified number of components. The model space can be always extended so it is infeasible to make such kind of full Bayesian inference.

Second, if every single model is close to one another and sample size is small, the mixture model can face non-identification or instability problem, unless a strong prior is added. Since the mixture model is relatively complex, this leads to a poor small sample behavior.

Figure 6 shows a comparison of mixture model and our model averaging methods in a numerical experiment, in which the true model is

[TABLE]

and there are 3 candidate models, each containing one covariate:

[TABLE]

In the simulation, we generate the design matrix by $\mathrm{Var}(X_{i})=1$ and $\mathrm{Cor}(X_{i},X_{j})=\rho$. $\rho$ determines how correlated these models are and it ranges from $-0.3$ to 0.9.

Figure 6 shows that both the performance of mixture models and single model selection are worse than any other model averaging method we suggest, even though the mixture model takes much longer time to run (about 30 more times) than stacking or Pseudo-BMA+. When the sample size is small, the mixture model is too complex to fit. On the other hand, stacking of predictive distributions and Pseudo-BMA+ outperform all other methods with a moderate sample size.

Clarke (2003) argues that the effect of (point estimation) stacking only depends on the space spanned by the model list, hence he suggests putting those “independent” models in the list. Figure 6 shows high correlations do not hurt stacking and Pseudo-BMA+ in this example.

3.4 Variational inference with different initial values

In Bayesian inference, the posterior density of parameters $\theta=\{\theta_{1},\dots,\theta_{m}\}$ given observed data $y=\{y_{1},\dots,y_{n}\}$ can be difficult to compute. Variational inference can be used to give a fast approximation for $p(\theta|y)$ (for a recent review, see Blei et al., 2017). Among a family of distribution $\mathcal{Q}$, we try to find the member of that family such that the Kullback-Leibler divergence to the true distribution is minimized:

[TABLE]

One widely used variational family is mean-field family where parameters are assumed to be mutually independent $\mathcal{Q}=\{q(\theta):q(\theta_{1},\dots,\theta_{m})=\prod_{j=1}^{m}q_{j}(\theta_{j})\}$. Some recent progress is made to run variational inference algorithm in a black-box way. For example, Kucukelbir et al. (2017) implement Automatic Variational Inference in Stan. Assuming all parameters $\theta$ are continuous and model likelihood is differentiable, it transforms $\theta$ into real coordinate space ${\rm I\!R}^{m}$ through $\zeta=T(\theta)$ and uses normal approximation $p(\zeta|\mu,\sigma^{2})=\prod_{j=1}^{m}\hbox{N}(\zeta_{j}|\mu_{j},\sigma_{j}^{2})$. Plugging this into 3.1 leads to an optimization problem over $(\mu,\sigma^{2})$, which can be solved by stochastic gradient ascent. Under some mild condition, it eventually converges to a local optimum $q^{*}(\theta)$. However, $q^{*}(\theta)$ may depend on initialization since such optimization problem is in general non-convex, particularly when the true posterior density $p(\theta|y)$ is multi-modal.

Stacking of predictive distributions and Pseudo-BMA+ weighting can be used to average several sets of posterior draws coming from different approximation distribution. To do this, we repeat the variational inference $K$ times. At time $k$, we start from a random initial point and use stochastic gradient ascent to solve the optimization problem 3.1, ending up with an approximation distribution $q^{*}_{k}(\theta)$. Then we draw $S$ samples $\{\theta_{k}^{(1)},\dots,\theta_{k}^{(S)}\}$ from $q^{*}_{k}$ calculate the importance ratio $r_{i,k}^{s}$ defined in 2.4 as

[TABLE]

After this, the remaining steps follow as before. We obtain stacking or pseudo-BMA+ weights $w_{k}$ and average all approximation distribution as $\sum_{k=1}^{K}w_{k}q^{*}_{k}$.

Figure 7 gives a simple numerical example that the averaging strategy helps adjust the optimization uncertainty of initial values. Suppose the data is two-dimensional $y=(y^{(1)},y^{(2)})$ and the parameter is $(\mu_{1},\mu_{2})\in{\rm I\!R}^{2}$. The likelihood $p(y|\mu_{1},\mu_{2})$ is given by

[TABLE]

A prior is assigned on $\mu_{1}-\mu_{2}$:

[TABLE]

We generate two observations $(y^{(1)}_{1}=3,y^{(2)}_{1}=2)$ and $(y^{(1)}_{2}=-2,y^{(2)}_{2}=-2)$. The first panel shows the true posterior distribution of $\mathbf{\mu}=(\mu_{1},\mu_{2})$, which is bimodal. We run mean-field normal variational inference in Stan, with two initial values set to be $(\mu_{1},\mu_{2})=(5,5)$ and $(-5,-5)$ separately. This produces two distinct approximation distribution as shown in panel 2 and 3. We then draw 1000 samples each from these two approximation distribution and use stacking or Pseudo-BMA+ to combine them. The lower 2 panels show the averaged posterior distribution. Though neither can recover the true distribution, the averaged version is closer to it, measured by KL divergence to the true one.

3.5 Proximity and directional models of voting

Adams et al. (2004) use US Senate voting data from 1988 to 1992 to study voters’ preference on the candidates who propose policies that are similar to them. They introduce two similar variables that indicate the distance between voters and candidates. Proximity voting comparison represents the $i$-th voter’s comparison between the candidates’ ideological positions:

[TABLE]

where $x_{i}$ represents the $i$-th voter’s preferred ideological position, and $x_{D}$ and $x_{R}$ represent the ideological positions of the Democratic and Republican candidates, respectively. In contrast, the $i$-th voter’s directional comparison is defined by

[TABLE]

where $X_{N}$ is the neutral point of the ideology scale.

Finally, all these comparison is aggregated in the party level, leading to two party-level variable Democratic proximity advantage and Democratic directional advantage. The sample size is $n=94$.

For both of these two variables, there are two ways to measure candidates’ ideological positions $x_{D}$ and $x_{R}$, which leads to two different datasets. In the Mean candidate dataset, they are calculated by taking the average of all respondents’ answers in the relevant state and year. In the Voter-specific dataset, they are calculate by using respondents’ own placements of the two candidates. In both datasets, there are 4 other party-level variables.

The two variables Democratic proximity advantage and Democratic directional advantage are highly correlated. Montgomery and Nyhan (2010) point out that Bayesian model averaging is an approach to helping arbitrate between competing predictors in a linear regression model. They average over all $2^{6}$ linear subset models excluding those that contain both variables Democratic proximity advantage and Democratic directional advantage, (i.e., 48 models in total). Each subset regression is with the form

[TABLE]

Accounting for the different complexity, they used the hyper-$g$ prior (Liang et al., 2008). Let $\phi$ to be the inverse of the variance $\phi=\frac{1}{\sigma^{2}}$. The hyper-$g$ prior with a prespecified hyperparameter $\alpha$ is,

[TABLE]

The first two columns of Figure 8 show the linear regression coefficients as estimated using least squares. The remaining columns show the posterior mean and standard deviation of the regression coefficients using BMA, stacking of predictive distributions, and Pseudo-BMA+, respectively. Under all three averaging strategies, the coefficient of proximity advantage is no longer statistically significantly negative, and the coefficient of directional advantage is shrunk. As fit to these data, stacking puts near-zero weights on all subset models containing proximity advantage, whereas Pseudo-BMA+ weighting always gives some weight to each model. In this example, averaging subset models by stacking or Pseudo-BMA+ weighting gives a way to deal with competing variables, which should be more reliable than BMA according to our previous argument.

3.6 Predicting well-switching behavior in Bangladesh

Many wells in Bangladesh and other South Asian countries are contaminated with natural arsenic. People whose wells have arsenic levels that exceed a certain threshold are encouraged to switch to nearby safe wells (for background details, see Chapter 5 in Gelman and Hill, 2006). We are analyzing a dataset including $3020$ respondents to find factors predictive of the well switching among all people with unsafe wells. The outcome variable is

[TABLE]

And we consider the following inputs:

• •

The distance (in meters) to the closest known safe well,

• •

The arsenic level of the respondent’s well,

• •

Whether any member of the household is active in the community association,

• •

The education level of the head of the household.

We start with what we call Model 1, a simple logistic regression with all variables above as well as a constant term:

[TABLE]

Model 2 contains the interaction between distance and arsenic level.

[TABLE]

Furthermore, it makes sense to us a nonlinear model for logit switching probability as a function of distance and arsenic level. We can use spline to capture that. Our Model 3 contains the B-splines for distance and arsenic level with polynomial degree 2,

[TABLE]

where $B_{dis}$ is the B-spline basis of distance with the form \bigl{(}B_{dis,1}(dist),\dots,B_{dis,q}(dist)\bigr{)} and $\alpha_{dis},\alpha_{ars}$ are vectors. We also fix the number of knots to be 10 for both distance and arsenic level. Model 4 and 5 are the similar models with 3-degree and 5-degree B-splines, respectively.

Next, we can add a bivariate spline to capture a nonlinear interaction,

[TABLE]

where $B_{dis,ars}$ is the bivariate spline basis with the degree to be $2\times 2,3\times 3,5\times 5$ in Model $6,7$ and $8$ respectively.

Figure 9 shows the inference results in all 8 models, which are summarized by the posterior mean, 50% confidence interval and 95% confidence interval of the probability of switching from an unsafe well as a function of distance or arsenic level. Any other variables such as assoc and edu are fixed at their means. It is not obvious to pick one best model. Spline models give a more flexible shape, but also introduce more variance for posterior estimation.

Finally, we run stacking of predictive distributions and Pseudo-BMA+ to combine these 8 models. The calculated model weights are listed above each panel in Figure 9. For both combination methods, Model 5 (univariate splines with 5th degree) accounts for the majority share. It is also worth pointing out that Model 8 is the most complicated one, but both stacking and Pseudo-BMA+ avoid overfitting by assigning a very small weight on it.

Figure 10 shows the posterior mean, 50% confidence interval, and 95% confidence interval of the switching probability in the stacking-combined model. Pseudo-BMA+ weighting gives a similar combination result for this example. At first glance, the combination looks quite similar to Model 5, while it may not seem necessary to put an extra 0.09 weight on Model 1 in stacking combination since Model 1 is completely contained in Model 5 if setting $\alpha_{dis}=\alpha_{ars}=0$. However, Model 5 is not perfect since it predicts that the posterior mean of switching probability will decrease as a function of distance to the nearest safe well, for very small distances. In fact, without further control, it is not surprising to find boundary fluctuation as a main drawback for higher order splines. Fortunately, we notice this decrease trend around the left boundary is a little bit flatter in the combined distribution since the combination contains part of straightforward logistic regression (in stacking weights) or lower order splines (in Pseudo-BMA+ weights). In this example the sample size $n=3020$ is large, hence we have reasons to believe stacking of predictive distributions gives the optimal combination.

4 Discussion

4.1 Sparse structure and high dimensions

Yang and Dunson (2014) propose to estimate a linear combination of point forecasts, $f=\sum_{k=1}^{K}w_{k}f_{k}$, using a Dirichlet aggregation prior, w\sim\mathrm{Dirichlet}\bigl{(}\frac{\alpha}{K^{\gamma}},\dots,\frac{\alpha}{K^{\gamma}}\bigr{)}, to pull toward structure, and estimating the weights $w_{k}$ using adaptive regression rather than cross-validation. They show that the combination under this setting can achieve the minimax squared risk among all convex combinations,

[TABLE]

where $F_{0}=(f:||f||_{\infty}\leq 1).$

Similar to our problem, when the dimension of model space is high, it can make sense to assign a strong prior can to the weights in estimation equation (2.2) to improve the regularization, using a hierarchical prior to pull toward sparsity if that is desired.

4.2 Constraints and regularity

In point estimation stacking, the simplex constraint is the most widely used regularization so as to overcome potential problems with multicollinearity. Clarke (2003) suggests relaxing the constraint to make it more flexible.

When combining distributions, there is no need to worry about multicollinearity except in degenerate cases. But in order to guarantee a meaningful posterior predictive density, the simplex constraint becomes natural, which is satisfied automatically in BMA and Pseudo-BMA weighting. As mentioned in the previous section, stronger priors can be added.

Another assumption is that the separate posterior distributions are combined linearly and with weights that are positive and sum to 1. There could be gains from going beyond convex linear combinations. For instance, in the subset regression example when each individual model is a univariate regression, the true model distribution is a convolution instead of a mixture of each possible models distribution. Both of them lead to the additive model in the point estimation, so stacking of the means is always valid, while stacking of predictive distributions is not possible to recover the true model in the convolution case.

Our explanation is that when the model list is large, the convex span should be large enough to approximate the true model. And this is the reason why we prefer adding stronger priors to make the estimation of weights stable in high dimensions.

4.3 General recommendations

The methods discussed in this paper are all based on the idea of fitting models separately and then combining the estimated predictive distributions. This approach is limited in that it does not pool information between the different model fits: as such, it is only ideal when the $K$ different models being fit have nothing in common. But in that case we would prefer to fit a larger super-model that includes the separate models as special cases, perhaps using an informative prior distribution to ensure stability in inferences.

That said, in practice it is common for different sorts of models to be set up without any easy way to combine them, and in such cases it is necessary from a Bayesian perspective to somehow aggregate their predictive distributions. The often-recommended approach of Bayesian model averaging can fail catastrophically in that the required Bayes factors can depend entirely on arbitrary specifications of noninformative prior distributions. Stacking is a more promising general method in that it is directly focused on performance of the combined predictive distribution. Based on our theory, simulations, and examples, we recommend stacking (of predictive distributions) for the task of combining separately-fit Bayesian posterior predictive distributions. As an alternative, Pseudo-BMA+ is computationally cheaper and can serve as an initial guess for stacking. The computations can be done in R and Stan, and the optimization required to compute the weights connects directly to the predictive task.

Appendix A. Implementation in Stan and R

The $n\times K$ matrix of cross-validated log likelihood values, $\{p(y_{i}|y_{-i},M_{k})\}_{i=1,\dots,n,k=1,\dots,K}$, can be computed from the generated quantities block in a Stan program, following the approach of Vehtari et al. (2017). For the example in Section 3.2, the $k$-th model is a linear regression with the $k$-th covariates. We put the corresponding Stan code in the file regression.stan:

data { int n; int p; vector[n] y; matrix[n, p] X; } parameters { vector[p] beta; real<lower=0> sigma; } transformed parameters { vector[n] theta; theta = X * beta; } model { y ~ normal(theta, sigma); beta ~ normal(0, 10); sigma ~ gamma(0.1, 0.1); } generated quantities { vector[n] log_lik; for (i in 1:n) log_lik[i] = normal_lpdf(y[i] | theta[i], sigma); }

In R we can simulate the likelihood matrices from all $K$ models and save them as a list:

library("rstan") log_lik_list <- list() for (k in 1:K){

Fit the k-th model with Stan

fit <- stan("regression.stan", data=list(y=y, X=X[,k], n=length(y), p=1)) log_lik_list[[k]] <- extract(fit)[["log_lik"]] }

The function model_weights() in the loo package111 The package can be downloaded from https://github.com/stan-dev/loo/tree/yuling-stacking in R can give model combination weights according to stacking of predictive distributions, Pseudo-BMA and Pseudo-BMA+ weighting. We can choose whether to use the Bayesian bootstrap to make further regularization for Pseudo-BMA, if computation time is not a concern.

model_weights_1 <- model_weights(log_lik_list, method="stacking") model_weights_2 <- model_weights(log_lik_list, method="pseudobma", BB=TRUE)

in one simulation with six models to combine, the output gives us the computed weights under each approach:

The stacking weights are: [1,] "Model 1" "Model 2" "Model 3" "Model 4" "Model 5" "Model 6" [2,] "0.25" "0.06" "0.09" "0.25" "0.35" "0.00"

The Pseudo-BMA+ weights using Bayesian Bootstrap are: [1,] "Model 1" "Model 2" "Model 3" "Model 4" "Model 5" "Model 6" [2,] "0.28" "0.05" "0.08" "0.30" "0.28" "0.00"

For reasons discussed in the paper, we generally recommend stacking for combining separate Bayesian predictive distributions.