The Prior Influence Function in Variational Bayes
Abstract
In Bayesian analysis, the posterior follows from the data and a choice of a prior and a likelihood. One hopes that the posterior is robust to reasonable variation in the choice of prior, since this choice is made by the modeler and is often somewhat subjective. A different, equally subjectively plausible choice of prior may result in a substantially different posterior, and so different conclusions drawn from the data. Were this to be the case, our conclusions would not be robust to the choice of prior. To determine whether our model is robust, we must quantify how sensitive our posterior is to perturbations of our prior.
Variational Bayes (VB) methods are fast, approximate methods for posterior inference. As with any Bayesian method, it is useful to evaluate the robustness of a VB approximate posterior to changes in the prior. In this paper, we derive VB versions of classical nonparametric local robustness measures. In particular, we show that the influence function of Gustafson (2000) has a simple, easytocalculate closed form expression for VB approximations. We then demonstrate how local robustness measures can be inadequate for nonlocal prior changes, such as replacing one prior entirely with another. We propose a simple approximate nonlocal robustness measure and demonstrate its effectiveness on a simulated data set.
subsecref \newrefsubsecname = \RSsectxt \RS@ifundefinedthmref \newrefthmname = theorem \RS@ifundefinedlemref \newreflemname = lemma
1 Local robustness and the influence function
Bayesian robustness studies how changes to the model (i.e., the prior and likelihood) and to the data affect the posterior. If important aspects of the posterior are meaningfully sensitive to subjectively reasonable perturbations of the inputs, then the posterior is “nonrobust” to these perturbations. In this paper, we focus on quantifying the sensitivity of posterior means to perturbations of the prior – either infinitesimally mixing or completely replacing the original prior with another “contaminating prior”. Our methods allow fast estimation of sensitivity to any contaminating prior without refitting the model. We follow and extend the work of Gustafson (1996) and Gustafson (2000) to variational Bayes and to approximate nonlocal measures of sensitivity. For a more general review of Bayesian robustness, see Berger et al. (2000).
We will now define some terminology. Denote our data points by with . Denote our parameter by the vector . We will suppose that we are interested in the robustness of our prior to a scalar parameter where our prior can be written as . Let denote the posterior distribution of with prior given by and conditional on , as given by Bayes’ Theorem: .
A typical end product of a Bayesian analysis might be a posterior expectation of some function, , which is a functional of and . Local robustness considers how much changes locally in response to small perturbations in the value of (Gustafson, 2000). In the present work, we consider mixing our original prior, , with some known alternative functional form, :
(1) 
This is known as epsilon contamination (the subscript stands for “contamination”), and its construction guarantees that the perturbed prior is properly normalized. The contaminating prior, need not be in the same parametric family as , so as ranges over all possible priors, (LABEL:epsilon_contamination) represents an expressive class of perturbations. Under mild assumptions (given in A), the local sensitivity measure at the prior given by a particular is
(2) 
The definition in (LABEL:local_robustness) depends on a choice of , which we denote with a superscript on . At , we recover the local sensitivity around , which we denote .
Rather than choose some finite set of and calculate their corresponding , one can work with a single function that summarizes the effect of any , called the “influence function” (Gustafson, 2000). Observing that (LABEL:local_robustness) is a linear functional of when , , and are fixed, the influence function (when it exists) is defined as the linear operator that characterizes the dependence of on :
(3) 
At , we recover the local sensitivity around , which we denote . When perturbing a lowdimensional marginal of the prior, is an easytovisualize summary of the effect of sensitivity to an arbitrary using quantities calculated only under (see the example in 4 and the extended discussion in Gustafson (2000)). Additionally, the worst case prior in a suitably defined metric ball around is a functional of the influence function, as shown in Gustafson (2000).
2 Variational approximation and linear response
We now derive a version of (LABEL:local_robustness) for Variational Bayes (VB) approximations to the posterior. Recall that an variational approximate posterior is a distribution selected to minimize the KullbackLiebler (KL) divergence to across distributions in some class . Let denote the variational approximation to posterior . We assume that distributions in are smoothly parameterized by a finitedimensional parameter whose optimum lies in the interior of some feasible set .
We would like to calculate the local robustness measures of 1 for the variational approximation , but a direct evaluation of the covariance in (LABEL:local_robustness) can be misleading. For example, a common choice of the approximating family is the class of distributions that factorize across . This is known as the “mean field approximation” (Wainwright and Jordan, 2008). By construction, a mean field approximation does not model covariances between independent components of , so a naive estimate of the covariance in (LABEL:local_robustness) may erroneously suggest that the prior on one component of cannot affect the posterior on another.
However, for VB approximations, we can evaluate the derivative on the left hand side of (LABEL:local_robustness) directly. Using linear response variational Bayes (LRVB) (Giordano et al., 2016, 2015), we have
(4)  
It follows immediately from the definition in (LABEL:influence_definition) that we can define the variational influence function
(5) 
that captures the sensitivity of just as captures the sensitivity of .
The VB versions of epsilon sensitivity measures have some advantages and disadvantages relative to using Markov Chain Monte Carlo (MCMC) to evaluate the exact sensitivities in 1. Using MCMC samples from , one can form a Monte Carlo estimate of the covariance in (LABEL:local_robustness), though the sample variance may be infinite when has heavier tails than . The extent to which this is a real problem in practice will vary. Similarly, one must take care in numerically evaluating (LABEL:lrvb_epsilon_sensitivity), since naively sampling from may also result in infinite variance due to the term in the denominator. Since we have a closed form for , we can instead evaluate (LABEL:lrvb_epsilon_sensitivity) as an integral over using importance sampling, as described in D. Still, providing efficient estimates of (LABEL:lrvb_epsilon_sensitivity) for highdimensional, nonconjugate, heavytailed remains a challenge. Finally, in contrast to (LABEL:influence_definition), where we do not generally have a closed form expression for , every term in (LABEL:lrvb_influence_function) is known. This means it is easier to evaluate the influence function for VB approximations than from MCMC draws, especially far from the posterior.
3 Nonlocal approximation
Equation (3) quantifies the effect of adding an infinitesimal amount of the contaminating prior. In practice, we may also want to evaluate intermediate values of , particularly , which represents completely replacing with . Since may be quite different from , this is a nonlocal robustness measure. For MCMC samples, one can use importance sampling, which is essentially equivalent to evaluating the covariance in (LABEL:local_robustness) (with the same problem of infinite variance – see C). For VB, however, we either need to refit the model for each new prior (which may be time consuming) or somehow use the local information at . In this paper, we investigate the latter. For the remainder of this section, since our results are general, we will discuss using local information in . However, the reader should keep in mind that the ultimate goal is to apply the insights gained to the variational approximation .
One might hope to linearly extrapolate from to using the slope at . That is, we might hope that . However, as we will now show, this is not realistic when one of the two priors is more consistent with the data than the other. Inspection of (LABEL:influence_definition) shows that posterior expectations are highly sensitive to perturbations of priors which are inconsistent with the data: if is small in an area of the space where is not small, then the influence function will be quite large. The model will have high sensitivity to any contaminating prior, , that is more consistent with the model than at . In particular, this is true at where if is inconsistent with the data. In fact, as we show in A,
(6) 
When the model evidence is very different for and , e.g. when as in 4, the extrapolated slope can be quite different from the effect of replacing completely replacing with .
However, as grows away from zero and the new prior is taken into account, the influence function will shrink. Observe that as a function of (LABEL:local_robustness), one can show (see A) that
(7) 
For or , this bound is vacuous, since the ratio can be arbitrarily large in areas assigned positive probability by . However, for intermediate values of , such as , the bound is quite strong. In other words, contamination with can have great influence on when is close to the boundaries of , but once is taken into account with intermediate , its influence is tightly bounded by (LABEL:sensitivity_bound). In this sense, the value of is most atypical of its value across the interval at its endpoints. A reallife example of exactly this phenomenon is shown in 4 in Fig. (1).
This suggests replacing the derivative at with an average of the derivative over the interval . To do this, note that the difficulty of the integral in (LABEL:global_local_sens_evidence) is the complicated dependence of on in . However, we can approximate the integral by keeping fixed at so that only depends on through . Under this approximation, the integral can be evaluated analytically (see B), giving the contaminating pseudodensity, , which represents the approximate effect of integrating over from to :
(8) 
In the notation , “mv” stands for “mean value”, by analogy with the mean value theorem for functions of real numbers. As shown in 4, using with rather than can represent a significant improvement over (LABEL:influence_definition) in practice when extrapolating to . We will focus on using with the variational approximations described in 2.
4 Experiments
We demonstrate our methods using simulated data from a hierarchical model described in E. Here, we will demonstrate that our sensitivity measures accurately predict the changes in VB solutions. We discuss the close match between the VB and MCMC results in E. The results below are for the sensitivity of the expectation to the prior , though similar results are easily calculated for any other lowdimensional prior marginal or posterior expectation.
We generate data using true parameters that are far from so that our model is not robust to perturbations, as can be seen by the large values of the influence function, which is pictured in the top left panel of Fig. (1). The posterior mean is shown with a black dot, indicating that, though large, the influence function is very highly concentrated around the posterior mean. The top right panel of Fig. (1) indicates how depends on near . The slope is very steep at , reflecting the fact that takes values much larger than near the posterior where the influence function is very high. However, it very quickly levels off for only slightly above zero.
The top right panel of Fig. (1) indicates that extrapolating from the slope at will radically overestimate the effect of replacing with . This is confirmed in the bottom left panel, which has the actual change in on the xaxis and on the yaxis. Clearly, the extrapolation is unrealistic. However, the right panel of Fig. (1) demonstrates that accurately matches the effect of replacing with . Note the different ranges in the yaxis (which prohibit plotting the two graphs on the same scale). The error bars represent importance sampling error.
Acknowledgements
The paper was greatly improved through the insightful comments of our anonymous reviewers. Ryan Giordano was supported by the Director, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program of the U.S. Department of Energy under Contract No. DEAC0205CH1123.
References
 Barnard et al. [2000] John Barnard, Robert McCulloch, and XiaoLi Meng. Modeling covariance matrices in terms of standard deviations and correlations, with application to shrinkage. Statistica Sinica, 10(4):1281–1312, 2000.
 Berger et al. [2000] James O. Berger, David Ríos Insua, and Fabrizio Ruggeri. Robust Bayesian analysis. In David Ríos Insua and Fabrizio Ruggeri, editors, Robust Bayesian Analysis, volume 152. Springer Science & Business Media, 2000.
 Durrett [2010] Rick Durrett. Probability: theory and examples. Cambridge University Press, 2010.
 Giordano et al. [2015] Ryan Giordano, Tamara Broderick, and Michael Jordan. Linear response methods for accurate covariance estimates from mean field variational Bayes. arXiv preprint arXiv:1506.04088, 2015.
 Giordano et al. [2016] Ryan Giordano, Tamara Broderick, Rachael Meager, Jonathan Huggins, and Michael Jordan. Fast robustness quantification with variational Bayes. arXiv preprint arXiv:1606.07153, 2016.
 Gustafson [1996] Paul Gustafson. Local sensitivity of posterior expectations. The Annals of Statistics, 24(1):174–195, 1996.
 Gustafson [2000] Paul Gustafson. Local robustness in Bayesian analysis. In David Ríos Insua and Fabrizio Ruggeri, editors, Robust Bayesian Analysis, volume 152. Springer Science & Business Media, 2000.
 Lewandowski et al. [2009] Daniel Lewandowski, Dorota Kurowicka, and Harry Joe. Generating random correlation matrices based on vines and extended onion method. Journal of multivariate analysis, 100(9):1989–2001, 2009.
 Meager [2015] Rachael Meager. Understanding the impact of microcredit expansions: A Bayesian hierarchical analysis of 7 randomised experiments. arXiv preprint arXiv:1506.06669, 2015.
 Rubin [1981] Donald B Rubin. Estimation in parallel randomized experiments. Journal of Educational and Behavioral Statistics, 6(4):377–401, 1981.
 Stan Team [2015] Stan Team. Stan Modeling Language Users Guide and Reference Manual, Version 2.8.0, 2015. URL http://mcstan.org/.
 Wainwright and Jordan [2008] M. J. Wainwright and M. I. Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends® in Machine Learning, 1(12):1–305, 2008.
Appendices
Appendix A Epsilon sensitivity
Throughout the paper, we make the following assumptions:

Assumption 1. For all , is strictly positive where has positive measure.

Assumption 2. Both and are bounded as a function of .
Under these assumptions, (LABEL:local_robustness) follows directly from Gustafson [1996] result 8. For completeness, we reproduce a slightly simpler proof under the slightly less clearlyarticulated assumption that we can exchange integration and differentiation as described, for example, in Durrett [2010, Appendix A.5]. Denote by any distribution parameterized by the scalar (not necessarily a prior). Then by direct differentiation,
(9) 
By applying (LABEL:differentiate_under_integral) to , we see that , so we can subtract to get
(10) 
To derive (LABEL:local_robustness), we simply observe that
Note that the assumptions also suffice to assure that the covariance is bounded.
Next, we observe the simple relationship between epsilon sensitivity at and the effect of replacing one prior with another. First, defining the normalizing constants
by straightforward manipulation we have
Finally, we derive (LABEL:sensitivity_bound) using (LABEL:local_robustness):
where
Appendix B Mean Value Contaminating Prior
Under the assumption that ,
Where we have used
Consequently, applying (LABEL:influence_definition) with the pseudodensity
represents an approximation to the quantity , which is the effect of completely replacing with .
Appendix C Comparison with MCMC importance sampling
In this section, we show that using importance sampling with MCMC samples to calculate the local sensitivity (LABEL:local_robustness) is precisely equivalent to using the same MCMC samples to estimate the covariance in (LABEL:covariance_sensitivity) directly. Suppose, without loss of generality, we have samples drawn from
If we could calculate the normalizing constants, the importance sampling estimate for would be
Differentiating the weights,
It follows that
which is precisely the MCMC estimate of the covariance given by (LABEL:local_robustness). In particular, when , we have
and the importance sampling estimate for replacing with is
which confirms that the importance sampling estimate is exactly the Monte Carlo analogue of the relation (LABEL:global_local_sens_evidence).
In general, we do not know and and must use instead
Then
which simply replaces the (possibly intractable) expectation with its MCMC estimate.
Appendix D Variational Bayes importance sampling
To evaluate (LABEL:lrvb_epsilon_sensitivity) requires approximate integration for which we using importance sampling:
Note that the influence function can be evaluated once for a large number of draws from , and then the weights and prior density can be quickly calculated for any perturbation , allowing for fast computation of sensitivity to any with little additional overhead.
Since the influence function is mostly concentrated around , we set to be but with quadrupled variance (so that standard deviations are doubled). Note that this choice of is a poor approximation of , which is nominally the target distribution for importance sampling. However, since is very small far from , it is an adequate approximation of the integral (LABEL:lrvb_epsilon_sensitivity). Formally, suppose that is concentrated on a set in the sense that for some small . Then the absolute error in evaluating (LABEL:lrvb_epsilon_sensitivity) on the set only is also bounded by :
As long as is chosen to have lighter tails than (which is determined by ), will decay quickly away from , and we can choose centered on . Consequently, we can think of as approximating rather than .
Appendix E Microcredit model
We simulate data using a variant of the analysis performed in [Meager, 2015], though with somewhat different prior choices. In Meager [2015], randomized controlled trials were run in seven different sites to try to measure the effect of access to microcredit on various measures of business success. Each trial was found to lack power individually for various reasons, so there could be some benefit to pooling the results in a simple hierarchical model. For the purposes of demonstrating robust Bayes techniques with VB, we will focus on the simpler of the two models in [Meager, 2015] and ignore covariate information.
We will index sites with (here, ) and business within a site by . The total number of observations was . In site and business we observe whether the business was randomly selected for increased access to microcredit, denoted , and the profit after intervention, . We follow [Rubin, 1981] and assume that each site has an idiosyncratic average profit, and average improvement in profit, , due to the intervention. Given , , and , the _observed profit is assumed to be generated according to
The site effects, , are assumed to come from an overall pool of effects and may be correlated:
The effects and the covariance matrix are unknown parameters that require priors. For we simply use a bivariate normal prior. However, choosing an appropriate prior for a covariance matrix can be conceptually difficult [Barnard et al., 2000]. Following the recommended practice of the software package STAN [Stan Team, 2015], we derive a variational model to accommodate the nonconjugate LKJ prior [Lewandowski et al., 2009], allowing the user to model the covariance and marginal variances separately. Specifically, we use
Diagonal matrix  
Covariance matrix  
We can then put independent priors on the scale of the variances, , and on the covariance matrix, . We model the inverse of with a Wishart variational distribution, and use the following priors:
The necessary expectations have closed forms with the Wishart variational approximation, as derived in Giordano et al. [2016]. In addition, we put a normal prior on and an inverse gamma prior on :
(11)  
The prior parameters used were:
As seen in Fig. (2) , the means in VB and MCMC match closely.
In our simulation, and are in , so the domain of the prior is twodimensional and can be easily visualized. We consider the problem of estimating the effect of replacing the prior on with a product of independent centered Student t priors. In the notation of 1, we take
We leave all other priors the same, i.e. and . In our case, we used and . We will present sensitivity of , the first component of the first toplevel effect. In the notation of 1, we are taking . Most of the computation is in generating draws and values for the importance sampling of the influence function, which can be done once and then reused for any choice of and .