Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
36,264 characters · 10 sections · 49 citation commands
Robustifying Empirical Bayes
The Gaussian sequence model can be viewed as a compound decision problem with observed $X_i \sim \mathcal{N} (\theta_i, 1), \; i = 1, \dots , n$. The objective is to estimate the $\theta \in \mathbb{R}^n$ subject to quadratic loss. We will denote the standard Gaussian density and cumulative by $\varphi$ and $\Phi$ respectively. Observations are assumed to be exchangeable, so their marginal density is given by, \[ f_G (x) = \int \varphi (x| \theta) dG(\theta), \] for some mixing distribution $G$. Were $G$ known the optimal (Bayes) decision rule is given by Tweedie's formula, efron.11 \[ \hat \theta_i = \delta^B (x_i) = x_i + f_G'(x_i)/f_G (x_i). \]
When $G$ is unknown various shrinkage procedures have been proposed, initiated by the fundamental papers of stein56 and robbins56. The extensive literature on Stein shrinkage has offered a rich assortment of practical frequentist and Bayesian procedures for improving upon the naive maximum likelihood estimator, $\delta (x_i) = x_i$ in terms of quadratic loss. Among these procedures more recently the nonparametric maximum likelihood estimator (NPMLE) of kw, \[ \hat G = \mbox{argmax}_{G \in \mathcal{G}} \big \{ \sum_{i=1}^n \log (f_G (x_i) \big \}, \] has been proposed as a plug-in estimator for $G$, jz, km and soloff. This $G$-modeling strategy -- in the terminology of e19 -- performs well in simulations, e.g. km, gk16, kg26, relative to alternatives that attempt to estimate $f_G$ directly or that make a priori assumptions about the form of $G$. However, it is obviously subject to the criticism that the Gaussian assumption on the likelihood is quite strong, and priors are never terribly convincing.
In what follows we consider two basic strategies for robustifying empirical Bayes procedures. The first, following a proposal of hodgeslehmann, seeks protection from excessive confidence in our initial prior on $G$ by bounding pointwise risk thereby offering a compromise between minimax and Bayes decision rules. The second, acknowledges scepticism about the strictly Gaussian form of $\varphi$, the distribution of the model noise. We find that in accordance with familiar robustness lore that modest modifications of of our initial prior or the Gaussian noise assumption can significantly improve performance of empirical Bayes decision rules while sacrificing only modest performance in the event that the initial prior or the Gaussian noise assumptions are valid.
Bayes risk in the Gaussian sequence model is, \[ r(G, \delta) = \int R(\delta, \theta) dG(\theta), \] with \[ R(\delta, \theta) =\mathbb{E}_\theta[(\delta(X) - \theta)^2] = \int ( \delta(x)-\theta)^2 \varphi(x - \theta) dx. \] Plugging the optimal Bayes rule back into $r(G, \delta)$, we obtain Brown's identity, Brown71:
The second equality follows from Stein's lemma, the fourth from the fact that $\int f_G^{''} (x) dx = 0$, and the last equality from the definition of Fisher information for distributions with absolutely continuous densities.
It may seem curious that Bayes risk of the optimal empirical Bayes rule for the Gaussian sequence model reduces to Fisher information for a scalar location parameter of the convolution distribution $\Phi * G$. We will exploit the latter connection in the next section to consider least favorable alternative priors for an initial prior in which we lack complete confidence. Such modified priors offer some compromise between strictly Bayesian and minimax procedures. Choosing contamination models by minimizing Fisher information is a classical strategy for choosing alternatives for the univariate Gaussian location and regression problems following huber64. The convolution form of the Fisher information for Bayes risk leads us back to an alternative proposal of mallows78 as well.
A natural objection to many empirical Bayes procedures is that they place unjustified reliance on an initial prior. While such procedures may still perform well with respect to ensemble risk, they may also fail spectacularly for some subpopulations or individuals. This concern underlies the limited translation proposal of efronmorris. We will see that bounding minimax risk by modifying an initial prior can often serve to soften the impact of these failings.
In an effort to balance minimax and Bayes solutions, hodgeslehmann proposed solving, \[ \min_\delta r(G_0,\delta) \; s.t. \; \max_\theta R(\delta, \theta) \leq 1 + t, \] for an initial prior $G_0$ and some $t > 0$. Since $\max_\theta R(X,\theta)=1$ corresponding to the worst pointwise risk among all decision rules, achieved by the MLE estimator $\delta(X) = X$, the Hodges and Lehmann modified decision rule is thereby constrained to do uniformly well over the entire parameter space with risk bounded by $1+t$, while minimizing the Bayes risk under prior $G_0$, hence called the restricted Bayes rule. They show that this is equivalent to solving, \[ \min_\delta \; \max_{G \in \mathcal{G}_\epsilon(G_0)} r(G, \delta), \] with $\mathcal{G}_\epsilon(G_0) = \{ G = (1 - \epsilon) G_0 + \epsilon H \}$ for some $\epsilon \in (0,1)$ depending upon $t$, where $H$ is an arbitrary distribution for $\theta$. Due to the minimax theorem we can switch the minimization and the maximization, and the resulting optimal rule $\delta^*$ is the posterior mean of $\theta$ under the least favorable prior from the class $\mathcal{G}_{\epsilon}(G_0)$. berger85 comments that “it is very difficult to determine such $\delta^*$; furthermore, this 'optimal' $\delta^*$ is usually extremely messy and difficult to work with.” On the contrary, with the aid of modern convex optimization techniques we find them quite tractable and elegant.
bickel83 considers the case with $G_0$ having point mass one at zero. Then, by the Brown identity, the Hodges and Lehmann problem is equivalent to solving, \[ \max_{G \in \mathcal{G}_\epsilon(\delta_0)} (1 - I (\Phi * G)) = \min_{G \in \mathcal{G}_\epsilon(\delta_0)} I (\Phi * G) \] This is the problem posed by mallows78 motivated by robustness considerations for time-series problems with additive outliers. Mallows conjectured that the least favorable $G$ would be discrete, supported on the integers with mass declining exponentially. bickel83 reports a modified conjecture of Donoho that relaxes the spacing of the Mallows mass points, but is otherwise similar. Neither conjecture seems to be strictly correct, but numerical computations confirm the nearly exponential decay of the mass. bickel1983minimizing provide a detailed discussion of the discrete nature of the Mallows solutions based on the analyticity of the objective function. See also johnstone94.
As noted by bickel1983minimizing and marazzi the mallows78 problem is convex. marazzi suggests a gridding strategy that imposes an exponentially declining mass condition. This produces a remarkably accurate solution for an initial Gaussian prior employing generic optimization software. Using modern convex optimization software, we show that the problem can be efficiently solved numerically for any prior distribution. Our implementation embodied in the function HodgesLehmann in the our REBayes package for the R language employs the Mosek mosek optimizer and provides a general interface for computing either the Huber or Mallows solutions for an arbitrary initial prior, $G_0$.
We now describe our procedure for the simplest (Dirac) initial prior, $G_0 = \delta_0$. Our objective is to solve \[ \underset{f \in \mathcal{K} }{\min} \int \frac{f'(x)^2}{f(x)} dx \] with \[ \mathcal{K}_\epsilon =\Big \{f = (1 - \epsilon) \int \varphi(x- \theta)d \delta_0(\theta) + \epsilon \int \varphi(x- \theta) dH(\theta)\Big \} \] This can be solved, on a grid of $x$, $\{x_1< x_2 < \dots <x_M\}$ and a grid of $\theta$ as $\{ \theta_1 < \theta_2 < \dots < \theta_L\}$ as the rotated quadratic cone convex optimization problem: \[ \min \sum_{i = 1}^M w_i \] subject to
Provided that the grids are sufficiently finely spaced interior point optimization in Mosek is capable of producing very accurate solutions very efficiently.
In Figure (ref) we illustrate a Mallows marginal density, $f^M$, its corresponding mixing distribution, $P^M$, and plot the log mass of the discrete mass points of $P^M$ at their respective locations. At first glance, it seems that the mass points are approximately equally spaced and have mass that declines exponentially. However, on closer examination the spacing of the mass points in the right tail are estimated to be: $\{ 1.96, 1.80, 1.70, 1.61, 1.52, 1.39, 1.29, 1.37, 1.55, 1.70, 1.91\}$, which seems sufficiently non-uniform to call the uniform spacing conjecture into question. djm suggest an alternative computational strategy for the Mallows problem using a parametric model that assumes equal spacing of the mass points. The grid for evaluation of $f$ is equally spaced from -30 to 30 with 500 points of evaluation. The grid for evaluation of $P^M$ is also equally spaced from -20 to 20 with 4003 points of evaluation. For purposes of illustration the mass at $\theta = 0$ is taken to be 0.2.
In the previous example we have taken the initial prior, $G_0$ as Dirac, but there is no obstacle to starting from any other initial prior. To provide some additional intuition about the nature of the Mallows solution for general $G_0$, we illustrate in Figure (ref) a plot of the pointwise risk function $R(\delta^* , \theta) := \mathbb{E}_\theta[(\delta^*(X)-\theta)^2]$ of the decision rule $\delta^*(\cdot)$, constructed as the posterior mean of $\theta$ using the least favorable Mallows prior \[ G^* (\theta) = (1 - \epsilon ) G_0 (\theta) + \epsilon H^* (\theta) \] where $G_0 (\theta)$ is taken to be an equally weighted mixture of two point masses at -2 and 2. The tangencies in this plot with the horizontal dotted line marking out $\sup_\theta R(\delta^*, \theta) = 1 + t$ coincide with the location of the mass points of the solution of $H^*$ indicated in the plot by the vertical green lines. The bound, $1 + t$, is the pointwise risk bound chosen to constrain the Hodges-Lehmann restricted Bayes rule, which can be constructed using Mallows's least favorable prior $G^*$. Since the initial prior $G_0$ places all its mass on the two points $\{-2,2\}$ the Mallows modification hedges this bet by placing a considerable mass at zero and exponentially declining mass at a few points below -2 and above +2. This figure is strongly reminiscent of Figure 5.5 of lindsay illustrating the location of mass points of the NPMLE.
We now consider several special cases of the Hodges and Lehmann restricted Bayes approach. In each case we consider not only the Mallows equivalent form of the Hodges-Lehmann modification, but also a Huber alternative that relaxes the Mallows objective of minimizing the Fisher information over convolutions by minimizing over the entire class of contamination distributions for $X$. Taking the least favorable density and plug into the Tweedie formula gives rise the Huber procedure to estimate $\theta$ for each value of $x$.
The Mallows rule has the obvious advantage that it yields a Bayes decision rule while the corresponding Huber procedure does not. This is particularly evident in the third example of Casella and Strawderman where the unrestricted Bayes rule is minimax, so the Hodges-Lehman's restriction on point-wise risk is unbinding and the restricted Bayes rule coincides with the unrestricted, but the Huber procedure is inadmissible. On the other hand there is something attractive about the Huber rules that it can be shown that they are necessarily monotone. donohoreeves propose an alternative based on the huber74 spline that minimizes Fisher information over a Kolmogorov neighborhood of the marginal density of $X$ specified by a finite number of evaluations of its quantile function. They then apply the Tweedie formula with the resulting least favorable density. An implementation of this procedure is also included in the REBayes package with the function HuberSpline, although we do not pursue it further in this paper.
Having examined several examples of the Hodges and Lehmann restricted Bayes rules for some simple initial prior distributions, we now consider an empirical Bayes counterpart with $G_0$ estimated by maximum likelihood as proposed by kw and anticipated by r50. Given a sample from the compound decision problem posed in the introduction, we consider the nonparametric maximum likelihood estimator $\hat G$ as the initial $G_0$ and then proceed to construct a modified prior according to the principles laid out by Hodges and Lehmann.
We illustrate the consequences of this in Figure (ref). Data is generated from the standard Gaussian sequence model with $G_0 \sim U[0,3]$. Heavy black vertical lines indicate the original NPMLE $\hat G$ while the red vertical lines indicate the mass points of the modified prior. While some alteration of the mass in the center of the estimated mixing distribution can be seen, the main change is the new mass points in the tails which decline exponentially in accordance with the Mallows's conjecture. This feature resembles the efronmorris limited translation estimator that imposed linear shrinkage in the center of the distribution, but eschewed shrinkage in the tails. In baseball terms: a few extremely good hitters deserve their exalted averages. Note that the restricted and unrestricted prior decision rules illustrated in the right panel of the figure agree quite closely on the support of the true $\theta$'s, but diverge sharply beyond this support.
In the following theorem, we establish that the excess Bayes risk, the difference between an oracle Mallows rule $\delta^M$ with known $G_0$ and an empirical Mallows rule $\hat \delta^M$ with estimated (NPMLE) $\hat G$ vanishes asymptotically. The proof makes use of machinery from variational analysis and the important feature that the score function of the oracle and EB Mallows least favorable density is uniformly bounded. The proof appears in Appendix (ref).
To evaluate the cost of imposing restrictions on the prior of the Hodges-Lehmann type we consider three examples in this section:
For each of these settings we compute mean squared error (MSE) for each of the following decision rules:
The last four rules are evaluated for four distinct values of $\epsilon \in \{ 0.05, 0.1, 0.2, 0.4\}$. All the simulations are based on 500 replications, for each compound decision problem.
The linear rules work reasonably well for the Gaussian and Uniform settings, however they perform poorly in the two-point setting. The cost of the Hodges-Lehmann restricted priors is modest for small $\epsilon$, but not surprisingly grows substantially when $\epsilon$ is larger. With only $n = 100$ observations, the NPMLE rule, $\delta_{\hat G}^B$, is too variable, but for the larger sample sizes it is nearly competitive with the (oracle) Bayes rules.
For each $G_0$ we can evaluate for any $\epsilon$ the corresponding worst case pointwise risk of each rule. The minimax rule achieves a worst case pointwise risk of one, while the Bayes rule typically has unbounded pointwise risk. For the Mallows rule, this can be evaluated by $\mathbb{E}_{\theta^*}[(\delta^M(X)-\theta^*)^2]$ where $\theta^*$ is any mass point of $H^*$ as discussed in Section (ref). For the Huber rule, we can show that for all the $G_0$ we considered, $\sup_\theta R(\delta^H, \theta) = 1+k_\epsilon^2$ with $k_\epsilon = \sup_{x} | (\log f_\epsilon^H(X))'|$ in which $f_\epsilon^H$ is the least favorable Huber density for the given $G_0$ and $\epsilon$. Details along with some simulation results appear in Appendix (ref).
\input tabs1/sim1a \input tabs1/sim2a \input tabs1/sim3a
Rather than robustifying the prior an alternative strategy is to robustify the likelihood. We will consider two variants of this: the first following huber64 and the second following mallows78.
The classical procedure of Huber for estimating a location parameter is easily adapted to the Gaussian sequence compound decision problem. In place of the Gaussian likelihood in the NPMLE problem we simply insert the Huber log likelihood with density, \[ \varphi (u) =
\] where $\epsilon$ and $k$ are linked by $2 \varphi(k)/k - 2 \Phi(-k) = \epsilon/(1-\epsilon)$. This density is least favorable, that is has minimal Fisher information for location, in the contamination model, \[ \Psi_\epsilon= \{ \Psi = (1 - \epsilon) \Phi + \epsilon H \}, \] over all symmetric distributions $H$. When $\epsilon = 1/2$ the least favorable Huber distribution is Laplace, or double exponential, and can be viewed as least favorable against asymmetric noise as well as symmetric.
mallows78 proposes to consider minimizing $I(\Phi * G)$ over $\mathcal{G}$, the set of all distributions with mass $1 - \epsilon$ at zero, provides an alternative to the Huber contamination model. Rather than assuming iid innovations each arising from the Huber mixture model, Mallows considers an additive outlier model in which with probability $\epsilon$ innovations are standard Gaussian, but occasionally are generated by the convolution $\Phi * H$.
In Figure (ref) we contrast the Huber and Mallows decision rules with the traditional Gaussian rule. The true mixing distribution, $G$, is chosen to be $U[0,3]$ so the Gaussian rule is itself somewhat curved, not the linear rule we would expect were $G$ itself Gaussian. In contrast the Huber and Mallows rules impose a more aggressive form of shrinkage. With Gaussian $\varphi$ extreme observations can be confidently attributed to signal, while the heavier tailed $\varphi$ of the Huber and Mallows rules tend to attribute such observations to noise. As we have seen previously, the Mallows rule oscillates around the Huber rule in the tails, but otherwise their behavior is quite similar.
When the usual Gaussian $\varphi$ is replaced by either the Huber or Mallows alternative in the nonparametric maximum likelihood estimation of $G$ solutions are accordingly more concentrated with fewer extreme mass points. This effect accentuates the more aggressive shrinkage effect observed in Figure (ref).
Both the Mallows and Huber least favorable contamination models offer principled alternatives to the strictly Gaussian noise model. They preserve the convexity of the underlying NPMLE problem and therefore can be easily implemented in software. In the next section we compare performance of several variants of these procedures for a few simulated compound decision settings.
We consider the compound decision problem with observations generated from, \[ Y_i = \theta_i + U_i, \quad i = 1, \dots, n, \] with $\theta_i$ and $U_i$ independent and each generated iidly from $G$ and $\Psi$ respectively. There are two choices of $G$: Either $G \sim U[0,3]$ or $G \sim 0.9 \delta_0 + 0.1 \delta_3$. And three choices of $\Psi$: $\Psi \sim 0.8 \Phi + 0.2 \Phi(\cdot / 3)$, $\Psi \sim \text{Laplace}$ and $\Psi = \Phi$, which we label Tukey, Laplace and Gauss respectively.
In Table (ref) we compare mean squared error performance of ten options with an infeasible oracle procedure that “knows” both the $\Psi$ and $G$ distributions. The experiment has 500 replications each with sample size $n = 500$. The competing feasible decision rules are: GLmix, the Gaussian NPMLE; Laplace, the Laplacian NPMLE; HLmix$(\epsilon)$, the Huber NPMLE; MLmix$(\epsilon)$, the Mallows NPMLE, with $\epsilon \in \{ 0.20, 0.10, 0.05, 0.025 \}$.
\input tabs2/sima.tex
It is evident from the table that the Gaussian NPMLE, GLmix, performs best when the noise distribution is actually Gaussian; however, when $\Psi \neq \Phi$ it pays to consider one of the alternatives. The Mallows NPMLE procedures seem to perform slightly better than the corresponding Huber methods, while the Laplace NPMLE, LLmix, which can be regarded as a “median-type” estimator performs surprisingly well over all the experimental settings.
We have considered two distinct strategies for robustifying empirical Bayes decision rules for the Gaussian sequence model. In the first motivated by the seminal paper of Hodges and Lehmann we would like protection against deviations from an initial Bayes prior. In the second we seek protection against non-Gaussian behavior in the noise distribution. Both strategies rely on the classical robustness proposals of huber64 and mallows78. Some combination of the two strategies is obviously possible, but choice of tuning parameters remains a delicate issue.