Dirichlet-Laplace priors for optimal shrinkage

Anirban Bhattacharya, Debdeep Pati, Natesh S. Pillai, David B. Dunson

Introduction

High-dimensional data have become commonplace in broad application areas, and there is an exponentially increasing literature on statistical and computational methods for big data. In such settings, it is well known that classical methods such as maximum likelihood estimation break down, motivating a rich variety of alternatives based on penalization and thresholding. There is a rich theoretical literature justifying the optimality properties of such penalization approaches , with fast algorithms and compelling applied results leading to routine use of L1L_{1} regularization in particular.

The overwhelming emphasis in this literature has been on rapidly producing a point estimate with good empirical and theoretical properties. However, in many applications, it is crucial to obtain a realistic characterization of uncertainty in the estimates of parameters and functions of the parameters, and in predictions of future outcomes. Usual frequentist approaches to characterize uncertainty, such as constructing asymptotic confidence regions or using the bootstrap, can break down in high-dimensional settings. For example, in regression when the number of subjects nn is much less than the number of predictors pp, one cannot naively appeal to asymptotic normality and resampling from the data may not provide an adequate characterization of uncertainty.

Most penalization approaches have a Bayesian interpretation as corresponding to the mode of a posterior distribution obtained under a shrinkage prior. For example, the wildly popular Lasso/L1L_{1} regularization approach to regression is equivalent to maximum a posteriori (MAP) estimation under a Gaussian linear regression model having a double exponential (Laplace) prior on the coefficients. Given this connection, it is natural to ask whether we can use the entire posterior distribution to provide a probabilistic measure of uncertainty. In addition to providing a characterization of uncertainty, a Bayesian perspective has distinct advantages in terms of tuning parameter choice, allowing key penalty parameters to be marginalized over the posterior distribution instead of relying on cross-validation. In addition, by inducing penalties through shrinkage priors, important new classes of penalties can be discovered that may outperform usual LqL_{q}-type choices.

From a frequentist perspective, we would like to be able to choose a default shrinkage prior that leads to similar optimality properties to those shown for L1L_{1} penalization and other approaches. However, instead of showing that a particular penalty leads to a point estimator having a minimax optimal rate of convergence under sparsity assumptions, we would rather like to show that the entire posterior distribution concentrates at the optimal rate, i.e., the posterior probability assigned to a shrinking neighborhood of the true parameter value converges to one, with the neighborhood size proportional to the frequentist minimax rate.

An amazing variety of shrinkage priors have been proposed in the Bayesian literature; however with essentially no theoretical justification for the performance of these priors in the high-dimensional settings for which they were designed. and provided conditions on the prior for asymptotic normality of linear regression coefficients allowing the number of predictors pp to increase with sample size nn, with requiring a very slow rate of growth and assuming p≤np\leq n. These results required the prior to be sufficiently flat in a neighborhood of the true parameter value, essentially ruling out shrinkage priors. considered shrinkage priors in providing simple sufficient conditions for posterior consistency in linear regression where the number of variables grows slower than the sample size, though no rate of contraction was provided.

In studying posterior contraction in high-dimensional settings, it becomes clear that it is critical to understand several aspects of the prior distribution on the high-dimensional space, including (but not limited to) the prior concentration around sparse vectors and the implied dimensionality of the prior. Specifically, studying the reduction in dimension induced by shrinkage priors is challenging due to the lack of exact zeros, with the prior draws being sparse in only an approximate sense. This substantial technical hurdle has prevented any previous results (to our knowledge) on posterior concentration in high-dimensional settings for shrinkage priors. In fact, investigating these properties is critically important not just in studying frequentist optimality properties of Bayesian procedures but for Bayesians in obtaining a better understanding of the behavior of their priors and choosing associated hyperparameters. Without such technical handle, it becomes an art to use intuition and practical experience to indirectly induce a shrinkage prior, while focusing on Gaussian scale families for computational tractability. Some beautiful classes of priors have been proposed by among others, with showing that essentially all existing shrinkage priors fall within the Gaussian global-local scale mixture family. One of our primary goals is to obtain theory that can allow evaluation of existing priors and design of novel priors, which are appealing from a Bayesian perspective in allowing incorporation of prior knowledge and from a frequentist perspective in leading to minimax optimality under weak sparsity assumptions.

A new class of shrinkage priors

For concreteness, we focus on the widely studied normal means problem (see, for example, and references therein); although most of the ideas developed in this paper generalize directly to high-dimensional linear and generalized linear models. In the normal means setting, one aims to estimate a nn-dimensional meanFollowing standard practice in this literature, we use nn to denote the dimensionality and it should not be confused with the sample size. based on a single observation corrupted with i.i.d. standard normal noise:

In a fully Bayesian framework, it is common to place a beta prior on π\pi, leading to a beta-Bernoulli prior on the model size, which conveys an automatic multiplicity adjustment . In a beautiful recent paper, established that prior (3) with an appropriate beta prior on π\pi and suitable tail conditions on gθg_{\theta} leads to a minimax optimal rate of posterior contraction, i.e., the posterior concentrates most of its mass on a ball around θ0\theta_{0} of squared radius of the order of qnlog⁡(n/qn)q_{n}\log(n/q_{n}):

where M>0M>0 is a constant and sn2=qnlog⁡(n/qn)s_{n}^{2}=q_{n}\log(n/q_{n}). obtained consistency in model selection using point-mass mixture priors with appropriate data-driven hyperparameters.

2 Global-local shrinkage rules

Although point mass mixture priors are intuitively appealing and possess attractive theoretical properties, posterior sampling requires a stochastic search over an enormous space, leading to slow mixing and convergence . Computational issues and consideration that many of the θj\theta_{j}s may be small but not exactly zero has motivated a rich literature on continuous shrinkage priors; for some flavor of the vast literature refer to . noted that essentially all such shrinkage priors can be represented as global-local (GL) mixtures of Gaussians,

where τ\tau controls global shrinkage towards the origin while the local scales {ψj}\{\psi_{j}\} allow deviations in the degree of shrinkage. If gg puts sufficient mass near zero and ff is appropriately chosen, GL priors in (5) can intuitively approximate (3) but through a continuous density concentrated near zero with heavy tails.

GL priors potentially have substantial computational advantages over point mass priors, since the normal scale mixture representation allows for conjugate updating of θ\theta and ψ\psi in a block. Moreover, a number of frequentist regularization procedures such as ridge, lasso, bridge and elastic net correspond to posterior modes under GL priors with appropriate choices of ff and gg. For example, one obtains a double-exponential prior corresponding to the popular L1L_{1} or lasso penalty if ff is an exponential distribution. However, unlike point mass priors (3), many aspects of shrinkage priors are poorly understood, with the lack of exact zeroes compounding the difficulty in studying basic properties, such as prior expectation, tail bounds for the number of large signals, and prior concentration around sparse vectors. Hence, subjective Bayesians face difficulties in incorporating prior information regarding sparsity, and frequentists tend to be skeptical due to the lack of theoretical justification.

This skepticism is somewhat warranted, as it is clearly the case that reasonable seeming priors can have poor performance in high-dimensional settings. For example, choosing π=1/2\pi=1/2 in prior (3) leads to an exponentially small prior probability of 2−n2^{-n} assigned to the null model, so that it becomes literally impossible to override that prior informativeness with the information in the data to pick the null model. However, with a beta prior on π\pi, this problem can be avoided . In the same vein, if one places i.i.d. \mboxN(0,1)\mbox{N}(0,1) priors on the entries of θ\theta, then the induced prior on ∥θ∥\left\|\theta\right\| is highly concentrated around n\sqrt{n} leading to misleading inferences on θ\theta almost everywhere. Although these are simple examples, similar multiplicity problems can transpire more subtly in cases where complicated models/priors are involved and hence it is fundamentally important to understand properties of the prior and the posterior in the setting of (1).

3 Dirichlet-kernel priors

Let us revisit the global-local specification (5). Integrating out the local scales ψj\psi_{j}’s, (5) can be equivalently represented as a global scale mixture of a kernel K(⋅)\mathcal{K}(\cdot),

These traditional choices lead to a kernel which is bounded in a neighborhood of zero. However, if one instead uses a half Cauchy prior ψj1/2∼\mboxCa+(0,1)\psi_{j}^{1/2}\sim\mbox{Ca}_{+}(0,1), then the resulting horseshoe kernel is unbounded with a singularity at zero. This phenomenon coupled with tail robustness properties leads to excellent empirical performance of the horseshoe. However, the joint distribution of θ\theta under a horseshoe prior is understudied and further theoretical investigation is required to understand its operating characteristics. One can imagine that it concentrates more along sparse regions of the parameter space compared to common shrinkage priors since the singularity at zero potentially allows most of the entries to be concentrated around zero with the heavy tails ensuring concentration around the relatively small number of signals.

In (7), K\mathcal{K} is any symmetric (about zero) unimodal density with exponential or heavier tails; for computational purposes, we shall restrict attention to the class of kernels that can be represented as scale mixture of normals . While previous shrinkage priors in the literature obtain marginal behavior similar to the point mass mixture priors (3), our construction aims at resembling the joint distribution of θ\theta under a two-component mixture prior. Constraining ϕ\phi on Sn−1\mathcal{S}^{n-1} restrains the degrees of freedom of the ϕj\phi_{j}’s, offering better control on the number of dominant entries in θ\theta. In particular, letting ϕ∼\mboxDir(a,…,a)\phi\sim\mbox{Dir}(a,\ldots,a) for a suitably chosen aa allows (7) to behave like (3) jointly, forcing a large subset of (θ1,…,θn)(\theta_{1},\ldots,\theta_{n}) to be simultaneously close to zero with high probability.

We focus on the Laplace kernel from now on for concreteness, noting that all the results stated below can be generalized to other choices. The corresponding hierarchical prior given τ\tau,

is referred to as a Dirichlet–Laplace prior, denoted θ∣τ∼\mboxDLa(τ)\theta\mid\tau\sim\mbox{DL}_{a}(\tau).

To understand the role of ϕ\phi, we undertake a study of the marginal properties of θj\theta_{j} conditional on τ\tau, integrating out ϕj\phi_{j}. The results are summarized in Proposition 2.1 below.

Thus, marginalizing over ϕ\phi, we obtain an unbounded kernel K\mathcal{K}, so that the marginal density of θj∣τ\theta_{j}\mid\tau has a singularity at 0 while retaining exponential tails. A proof of Proposition 2.1 can be found in the appendix.

The parameter τ\tau plays a critical role in determining the tails of the marginal distribution of θj\theta_{j}’s. We consider a fully Bayesian framework where τ\tau is assigned a prior gg on the positive real line and learnt from the data through the posterior. Specifically, we assume a \mboxgamma(λ,1/2)\mbox{gamma}(\lambda,1/2) prior on τ\tau with λ=na\lambda=na. We continue to refer to the induced prior on θ\theta implied by the hierarchical structure,

as a Dirichlet–Laplace prior, denoted θ∼\mboxDLa\theta\sim\mbox{DL}_{a}.

There is a recent frequentist literature on including a local penalty specific to each coefficient. The adaptive Lasso relies on empirically estimated weights that are plugged in. instead propose to sample the penalty parameters from a posterior, with a sparse point estimate obtained for each draw. These approaches do not produce a full posterior distribution but focus on sparse point estimates.

4 Posterior computation

The proposed class of DL priors leads to straightforward posterior computation via an efficient data augmented Gibbs sampler. Note that the \mboxDLa\mbox{DL}_{a} prior (9) can be equivalently represented as

We detail the steps in the normal means setting noting that the algorithm is trivially modified to accommodate normal linear regression, robust regression with heavy tailed residuals, probit models, logistic regression, factor models and other hierarchical Gaussian cases. To reduce auto-correlation, we rely on marginalization and blocking as much as possible. Our sampler cycles through (i) θ∣ψ,ϕ,τ,y\theta\mid\psi,\phi,\tau,y, (ii) ψ∣ϕ,τ,θ\psi\mid\phi,\tau,\theta, (iii) τ∣ϕ,θ\tau\mid\phi,\theta and (iv) ϕ∣θ\phi\mid\theta. We use the fact that the joint posterior of (ψ,ϕ,τ)(\psi,\phi,\tau) is conditionally independent of yy given θ\theta. Steps (ii) - (iv) together give us a draw from the conditional distribution of (ψ,ϕ,τ)∣θ(\psi,\phi,\tau)\mid\theta, since

Steps (i) – (iii) are standard and hence not derived. Step (iv) is non-trivial and we develop an efficient sampling algorithm for jointly sampling ϕ\phi. Usual one at a time updates of a Dirichlet vector leads to tremendously slow mixing and convergence, and hence the joint update in Theorem 2.2 is an important feature of our proposed prior; a proof can be found in the Appendix. Consider the following parametrization for the three-parameter generalized inverse Gaussian (giG) distribution: Y∼\mboxgiG(λ,ρ,χ)Y\sim\mbox{giG}(\lambda,\rho,\chi) if f(y)∝yλ−1e−0.5(ρy+χ/y)f(y)\propto y^{\lambda-1}e^{-0.5(\rho y+\chi/y)} for y>0y>0.

The joint posterior of ϕ∣θ\phi\mid\theta has the same distribution as (T1/T,…,Tn/T)(T_{1}/T,\ldots,T_{n}/T), where TjT_{j} are independently distributed according to a \mboxgiG(a−1,1,2∣θj∣)\mbox{giG}(a-1,1,2|\theta_{j}|) distribution, and T=∑j=1nTjT=\sum_{j=1}^{n}T_{j}.

The summary of each step are finally provided below.

To sample θ∣ψ,ϕ,τ,y\theta\mid\psi,\phi,\tau,y, draw θj\theta_{j} independently from a \mboxN(μj,σj2)\mbox{N}(\mu_{j},\sigma_{j}^{2}) distribution with

The conditional posterior of ψ∣ϕ,τ,θ\psi\mid\phi,\tau,\theta can be sampled efficiently in a block by independently sampling ψj∣ϕ,θ\psi_{j}\mid\phi,\theta from an inverse-Gaussian distribution \mboxiG(μj,λ)\mbox{iG}(\mu_{j},\lambda) with μj=ϕjτ/∣θj∣,λ=1\mu_{j}=\phi_{j}\tau/|\theta_{j}|,\lambda=1.

Sample the conditional posterior of τ∣ϕ,θ\tau\mid\phi,\theta from a \mboxgiG(λ−n,1,2∑j=1n∣θj∣/ϕj)\mbox{giG}(\lambda-n,1,2\sum_{j=1}^{n}|\theta_{j}|/\phi_{j}) distribution.

To sample ϕ∣θ\phi\mid\theta, draw T1,…,TnT_{1},\ldots,T_{n} independently with Tj∼\mboxgiG(a−1,1,2∣θj∣)T_{j}\sim\mbox{giG}(a-1,1,2|\theta_{j}|) and set ϕj=Tj/T\phi_{j}=T_{j}/T with T=∑j=1nTjT=\sum_{j=1}^{n}T_{j}.

Concentration properties of Dirchlet–Laplace priors

The formulation (10) is analytically convenient since the joint distribution factors as a product of marginals and the marginal density can be obtained analytically in Proposition 3.1 below. The proof follows from standard properties of the modified Bessel function ; a proof is sketched in the Appendix.

The marginal density Π\Pi of θj\theta_{j} for any 1≤j≤n1\leq j\leq n is given by

is the modified Bessel function of the second kind.

Figure 1 plots the marginal density (11) to compare with other common shrinkage priors.

If an=1/na_{n}=1/n instead, then (12) holds when qn≳log⁡nq_{n}\gtrsim\log n.

A proof of Theorem 3.1 can be found in Section 6. To best of our knowledge, Theorem 3.1 is the first result obtaining posterior contraction rates for a continuous shrinkage prior in the normal means setting or the closely related high-dimensional regression problem. Theorem 3.1 posits that when the parameter aa in the Dirichlet–Laplace prior is chosen, depending on the sample size, to be n−(1+β)n^{-(1+\beta)} for any β>0\beta>0 small, the resulting posterior contracts at the minimax rate (2), provided ∥θ0∥22≤qnlog⁡4n\left\|\theta_{0}\right\|_{2}^{2}\leq q_{n}\log^{4}n. Using the Cauchy–Schwartz inequality, ∥θ0∥12≤qn∥θ0∥22\left\|\theta_{0}\right\|_{1}^{2}\leq q_{n}\left\|\theta_{0}\right\|_{2}^{2} for θ0∈l0[qn;n]\theta_{0}\in l_{0}[q_{n};n] and the bound on ∥θ0∥2\left\|\theta_{0}\right\|_{2} implies that ∥θ0∥1≤qn(log⁡n)2\left\|\theta_{0}\right\|_{1}\leq q_{n}(\log n)^{2}. Hence, the condition in Theorem 3.1 permits each non-zero signal to grow at a (log⁡n)2(\log n)^{2} rate, which is a fairly mild assumption. Moreover, in a recent technical report, the authors showed that a large subclass of global-local priors (5) including the Bayesian lasso lead to a sub-optimal rate of posterior convergence; i.e., the expression in (12) converges to whenever ∥θ0∥22/qn→∞\left\|\theta_{0}\right\|_{2}^{2}/q_{n}\to\infty. Therefore, Theorem 3.1 indeed provides a substantial improvement over a large class of GL priors.

The choice an=n−(1+β)a_{n}=n^{-(1+\beta)} will be evident from the various auxiliary results in Section 3.1, specifically Lemma 3.3 and Theorem 3.4. The conclusion of Theorem 3.1 continues to hold when an=1/na_{n}=1/n under an additional mild assumption on the sparsity qnq_{n}. In Table 1 of Section 4, detailed empirical results are provided with an=1/na_{n}=1/n as a default choice.

The lower bound result alluded to in the previous paragraph precludes GL priors with polynomial tails, such as the horseshoe. We hope to address the polynomial tails case elsewhere, though based on strong empirical performance, we conjecture that the horseshoe leads to the optimal posterior contraction in a much broader domain compared to the Bayesian lasso and other common shrinkage priors.

In this section, we state a number of properties of the DL prior which provide a better understanding of the joint prior structure and also crucially help us in proving Theorem 3.1.

We first provide useful bounds on the joint density of the DL prior in Lemma 3.2 below; a proof can be found in the Appendix.

If min⁡1≤j≤∣S∣∣ηj∣>δ\min_{1\leq j\leq|S|}|\eta_{j}|>\delta for δ\delta small, then

If ∥η∥2≤m\|\eta\|_{2}\leq m for mm large, then

It is evident from Figure 1 that the univariate marginal density Π\Pi has an infinite spike near zero. We quantify the probability assigned to a small δ\delta-neighborhood of the origin in Lemma 3.3 below.

A proof of Lemma 3.3 can be found in the Appendix.

for some constant A>0A>0. If an=1/na_{n}=1/n instead, then (15) holds when qn≳log⁡nq_{n}\gtrsim\log n.

The posterior compressibility property in Theorem 3.4 ensures that the dimensionality of the posterior distribution of θ\theta (in an approximate sense) doesn’t substantially overshoot the true dimensionality of θ0\theta_{0}, which together with the bounds on the joint prior density near zero and infinity in Lemma 3.2 delivers the minimax rate in Theorem 3.1.

Simulation Study

The squared error loss corresponding to the posterior median averaged across simulation replicates is provided in Table 1. To offer further grounds for comparison, we have also tabulated the results for Lasso (LS), Empirical Bayes median (EBMed) as in The EBMed procedure was implemented using the package . , posterior median with a point mass prior (PM) as in and the posterior median corresponding to the horseshoe prior . For the fully Bayesian analysis using point mass mixture priors, we use a complexity prior on the subset-size, πn(s)∝exp⁡{−κslog⁡(2n/s)}\pi_{n}(s)\propto\exp\{-\kappa s\log(2n/s)\} with κ=0.1\kappa=0.1 and independent standard Laplace priors for the non-zero entries as in .Given a draw for ss, a subset SS of size ss is drawn uniformly. Set θj=0\theta_{j}=0 for all j∉Sj\notin S and draw θj,j∈S\theta_{j},j\in S i.i.d. from standard Laplace. The beta-bernoulli priors in (3) induce a similar prior on the subset size.

Even in this succinct summary of the results, a wide difference between the Bayesian Lasso and the proposed \mboxDL1/n\mbox{DL}_{1/n} is observed in Table 1, vindicating our theoretical results. The horseshoe performs similarly as the \mboxDL1/n\mbox{DL}_{1/n}. The superior performance of the \mboxDL1/n\mbox{DL}_{1/n} prior can be attributed to its strong concentration around the origin. However, in cases where there are several relatively small signals, the \mboxDL1/n\mbox{DL}_{1/n} prior can shrink all of them towards zero. In such settings, depending on the practitioner’s utility function, the singularity at zero can be softened using a \mboxDLa\mbox{DL}_{a} prior for a larger value of aa. In the next set of simulations, we report results for a=1/2a=1/2, whence computational gains arise as the distribution of TjT_{j} in (iv) turns out to be inverse-Gaussian (iG), for which exact samplers are available. In practice, one could as well use a discrete uniform prior on aa as in Section 5.

For visual illustration and comparison, we finally present the results from a single replicate in the first simulation setting with n=200n=200, qn=10q_{n}=10 and A=7A=7 in Figure 2 & 3. The blue circles indicate the entries of yy, while the red circles correspond to the posterior median of θ\theta. The shaded region corresponds to a 95%95\% point wise credible interval for θ\theta.

Prostate data application

We consider a popular dataset from a microarray experiment consisting of expression levels for 60336033 genes for 5050 normal control subjects and 5252 patients diagnosed with prostate cancer. The data takes the form of a 6033×1026033\times 102 matrix with the (i,j)(i,j)th entry corresponding to the expression level for gene ii on patient jj; the first 5050 columns correspond to the normal control subjects with the remaining 5252 for the cancer patients. The goal of the study is to discover genes whose expression levels differ between the prostate cancer patients (treatment) and normal subjects (control). A two sample tt-test with 100100 degrees of freedom was implemented for each gene and the resulting t-statistic tit_{i} was converted to a zz-statistic zi=Φ−1(T100(ti))z_{i}=\Phi^{-1}(T_{100}(t_{i})). Under the null hypothesis H0iH_{0i} of no difference in expression levels between the treatment and control group for the iith gene, the null distribution of ziz_{i} is \mboxN(0,1)\mbox{N}(0,1). Figure 4 shows a histogram of the zz-values, comparing it to a \mboxN(0,1)\mbox{N}(0,1) density with a multiplier chosen to make the curve integrate to the same area as the histogram. The shape of the histogram suggests the presence of certain interesting genes .

The classical Bonferroni correction for multiple testing flags only 66 genes as significant, while the two-group empirical Bayes method of found 139139 significant genes, being much less conservative. The local Bayes false discovery rate (fdr) control method identified 5454 genes as non-null. For detailed analysis of this dataset using existing methods, refer to .

To apply our method, we set up a normal means model zi=θi+ϵi,i=1,…,6,033z_{i}=\theta_{i}+\epsilon_{i},i=1,\ldots,6,033 and assign θ\theta a \mboxDLa\mbox{DL}_{a} prior. Instead of fixing aa, we use a discrete uniform prior on aa supported on the interval [1/6,000,1/2][1/6,000,1/2], with the support points of the form 10(k+1)/6,000,k=0,1,…,K10(k+1)/6,000,k=0,1,\ldots,K. Such a fully Bayesian approach allows the data to dictate the choice of the tuning parameter aa which is only specified up to a constant by the theory and also avoids potential numerical issues arising from fixing a=1/na=1/n when nn is large. Updating aa is straightforward since the full conditional distribution of aa is again a discrete distribution on the chosen support points.

We implemented the Gibbs sampler in Section 2.4 for 10,000 draws discarding a burn-in of 5,000. Mixing and convergence of the Gibbs sampler was satisfactory based on examination of trace plots, with the 5,000 retained samples having an effective sample size of 2369.2 averaged across the θi\theta_{i}’s. The computational time per iteration scaled approximately linearly with the dimension. The posterior mode of aa was at 1/201/20.

In this application, we expect there to be two clusters of ∣θi∣|\theta_{i}|s, with one concentrated closely near zero corresponding to genes that are effectively not differentially expressed and another away from zero corresponding to interesting genes for further study. As a simple automated approach, we cluster ∣θi∣|\theta_{i}|s at each MCMC iteration using kmeans with 2 clusters. For each iteration, the number of non-zero signals is then estimated by the smaller cluster size out of the two clusters. A final estimate (MM) of the number of non-zero signals is obtained by taking the mode over all the MCMC iterations. The MM largest (in absolute magnitude) entries of the posterior median are identified as the non-zero signals.

Using the above selection scheme, our method declared 128 genes as non-null. Interestingly, out of the 128 genes, 100 are common with the ones selected by EBMed. Also all the 54 genes obtained using FDR control form a subset of the selected 128128 genes. Horseshoe is overly conservative; it selected only 11 gene (index: 610) using the same clustering procedure; the selected gene was the one with the largest effect size (refer to Table 11.2 in ).

Proof of Theorem 3.1

For a sequence of positive real numbers rnr_{n} to be chosen later, let δn=rn/n\delta_{n}=r_{n}/n. Define Dn=∫∏i=1nfθ0i(yi)/fθi(yi) dΠ(θ)\mathcal{D}_{n}=\int\prod_{i=1}^{n}f_{\theta_{0i}}(y_{i})/f_{\theta_{i}}(y_{i})\,d\Pi(\theta). Let

Let rn2=qnlog⁡nr_{n}^{2}=q_{n}\log n. The proof of Theorem 3.1 is completed by deriving an upper bound to βS,j,i\beta_{S,j,i} in the following Lemma 6.1 akin to Lemma 5.4 in .

log⁡βS,j,i≤∣S∣log⁡(2j)+C(∣S∣+∣S0∣)log⁡n+C′rn2\log\beta_{S,j,i}\leq|S|\log(2j)+C(|S|+|S_{0}|)\log n+C^{\prime}r_{n}^{2}.

Let vq(r)v_{q}(r) denote the qq-dimensional Euclidean ball of radius rr centered at zero and ∣vq(r)∣|v_{q}(r)| denote its volume. For the sake of brevity, denote vq=∣vq(1)∣v_{q}=|v_{q}(1)|, so that ∣vq(r)∣=rqvq|v_{q}(r)|=r^{q}v_{q}. The numerator of (18) can be clearly bounded above by ∣v∣S∣(2jrn)∣sup⁡∣θj∣>δn ∀ j∈SΠS(θS)|v_{|S|}(2jr_{n})|\sup_{|\theta_{j}|>\delta_{n}\,\forall\,j\in S}\Pi_{S}(\theta_{S}). Since the set {∥θS0−θ0S0∥2<rn}\{\left\|\theta_{S_{0}}-\theta_{0S_{0}}\right\|_{2}<r_{n}\} is contained in the ball v∣S0∣(∥θ0S0∥2+rn)={∥θS0∥2≤∥θ0S0∥2+rn}v_{|S_{0}|}(\|\theta_{0S_{0}}\|_{2}+r_{n})=\{\|\theta_{S_{0}}\|_{2}\leq\|\theta_{0S_{0}}\|_{2}+r_{n}\} and ∥θ0S0∥2=∥θ0∥2\|\theta_{0S_{0}}\|_{2}=\|\theta_{0}\|_{2}, the denominator of (18) can be bounded below by ∣v∣S0∣(rn)∣inf⁡v∣S0∣(tn)ΠS0(θS0)|v_{|S_{0}|}(r_{n})|\inf_{v_{|S_{0}|}(t_{n})}\Pi_{S_{0}}(\theta_{S_{0}}), where tn=∥θ0∥2+rnt_{n}=\|\theta_{0}\|_{2}+r_{n}. Putting together these inequalities and invoking Lemma 3.2, we have

Using vq≍(2πe)q/2q−q/2−1/2v_{q}\asymp(2\pi e)^{q/2}q^{-q/2-1/2} (see Lemma 5.3 in ) and rn2≥qn=∣S0∣r_{n}^{2}\geq q_{n}=|S_{0}|, we can bound log⁡{rn∣S∣v∣S∣/(rn∣S0∣v∣S0∣)}\log\{r_{n}^{|S|}v_{|S|}/(r_{n}^{|S_{0}|}v_{|S_{0}|})\} from above by C(∣S∣log⁡n+rn2)C(|S|\log n+r_{n}^{2}). Therefore, we have

Now, since ∥θ0∥22≤qnlog⁡4n\left\|\theta_{0}\right\|_{2}^{2}\leq q_{n}\log^{4}n and rn2=qnlog⁡nr_{n}^{2}=q_{n}\log n, we have tn≲qn1/2log⁡2nt_{n}\lesssim q_{n}^{1/2}\log^{2}n and hence ∣S0∣3/4tn1/2≲qnlog⁡n=rn2|S_{0}|^{3/4}t_{n}^{1/2}\lesssim q_{n}\log n=r_{n}^{2}. Substituting in (20), we have

Substituting the upper bound for βS,j,i\beta_{S,j,i} obtained in Lemma 6.1, and noting that ∣S∣≤Aqn\left|S\right|\leq Aq_{n} and ∣NS,j∣≤eC∣S∣\left|N_{S,j}\right|\leq e^{C\left|S\right|}, the expression in the left hand side of (16) can be bounded above by

When an=n−(1+β)a_{n}=n^{-(1+\beta)}, the conclusion of Lemma 6.1 remains unchanged and the proof of Theorem 3.4 does not require qn≳log⁡nq_{n}\gtrsim\log n. The rest of the proof remains exactly the same.

Appendix

When a=1/na=1/n, ϕj∼\mboxBeta(1/n,1−1/n)\phi_{j}\sim\mbox{Beta}(1/n,1-1/n) marginally. Hence, the marginal distribution of θj\theta_{j} given τ\tau is proportional to

Substituting z=ϕj/(1−ϕj)z=\phi_{j}/(1-\phi_{j}) so that ϕj=z/(1+z)\phi_{j}=z/(1+z), the above integral reduces to

In the general case, ϕj∼\mboxBeta(a,(n−1)a)\phi_{j}\sim\mbox{Beta}(a,(n-1)a) marginally. Substituting z=ϕj/(1−ϕj)z=\phi_{j}/(1-\phi_{j}) as before, the marginal density of θj\theta_{j} is proportional to

The above integral can clearly be bounded below by a constant multiple of

The above expression clearly diverges to infinity as ∣θj∣→0|\theta_{j}|\to 0 by the monotone convergence theorem.

Proof of Theorem 2.2

Integrating out τ\tau, the joint posterior of ϕ∣θ\phi\mid\theta has the form

We now state a result from the theory of normalized random measures (see, for example, (36) in ). Suppose T1,…,TnT_{1},\ldots,T_{n} are independent random variables with TjT_{j} having a density fjf_{j} on (0,∞)(0,\infty). Let ϕj=Tj/T\phi_{j}=T_{j}/T with T=∑j=1nTjT=\sum_{j=1}^{n}T_{j}. Then, the joint density ff of (ϕ1,…,ϕn−1)(\phi_{1},\ldots,\phi_{n-1}) supported on the simplex Sn−1\mathcal{S}^{n-1} has the form

where ϕn=1−∑j=1n−1ϕj\phi_{n}=1-\sum_{j=1}^{n-1}\phi_{j}. Setting fj(x)∝1xδe−∣θj∣/xe−x/2f_{j}(x)\propto\frac{1}{x^{\delta}}e^{-|\theta_{j}|/x}e^{-x/2} in (22), we get

We aim to equate the expression in (23) with the expression in (21). Comparing the exponent of ϕj\phi_{j} gives us δ=2−a\delta=2-a. The other requirement n−1−nδ=λ−n−1n-1-n\delta=\lambda-n-1 is also satisfied, since λ=na\lambda=na. The proof is completed by observing that fjf_{j} corresponds to a \mboxgiG(a−1,1,2∣θj∣)\mbox{giG}(a-1,1,2|\theta_{j}|) when δ=2−a\delta=2-a.

Proof of Proposition 3.1

Proof of Lemma 3.2

Letting h(x)=log⁡Π(x)h(x)=\log\Pi(x), we have log⁡ΠS(η)=∑1≤j≤∣S∣h(ηj)\log\Pi_{S}(\eta)=\sum_{1\leq j\leq|S|}h(\eta_{j}).

We first prove (13). Since Π(x)\Pi(x), and hence h(x)h(x), is monotonically decreasing in ∣x∣|x|, and ∣ηj∣>δ|\eta_{j}|>\delta for all jj, we have log⁡ΠS(η)≤∣S∣h(δ)\log\Pi_{S}(\eta)\leq|S|h(\delta). Using Kα(z)≍z−αK_{\alpha}(z)\asymp z^{-\alpha} for ∣z∣|z| small and Γ(a)≍a−1\Gamma(a)\asymp a^{-1} for aa small, we have from (11) that Π(δ)≍a−1∣δ∣(a−1)\Pi(\delta)\asymp a^{-1}|\delta|^{(a-1)} and hence h(δ)≍(1−a)log⁡(δ−1)−log⁡a−1+C≤Clog⁡(δ−1)h(\delta)\asymp(1-a)\log(\delta^{-1})-\log a^{-1}+C\leq C\log(\delta^{-1}).

We next prove (14). Noting that Kα(z)≳e−z/zK_{\alpha}(z)\gtrsim e^{-z}/z for ∣z∣|z| large (section 9.7 of ), we have from (11) that −h(x)≤log⁡a−1+3/2log⁡∣x∣+2∣x∣-h(x)\leq\log a^{-1}+3/2\log|x|+\sqrt{2}\sqrt{|x|} for ∣x∣|x| large. Using Cauchy–Schwartz inequality twice, we have (∑j=1∣S∣∣ηj∣)4≤∣S∣3∥η∥22(\sum_{j=1}^{|S|}\sqrt{|\eta_{j}|})^{4}\leq|S|^{3}\|\eta\|_{2}^{2}, which implies ∑j=1∣S∣∣ηj∣≤∣S∣3/4∥η∥21/2≤∣S∣3/4m1/2\sum_{j=1}^{|S|}\sqrt{|\eta_{j}|}\leq|S|^{3/4}\|\eta\|_{2}^{1/2}\leq|S|^{3/4}m^{1/2}.

Proof of Lemma 3.3

where C>0C>0 is a constant independent of δ\delta. Using a bound for the incomplete gamma function from Theorem 2 of ,

for δ\delta small. The proof is completed by noting that (1/2)a(1/2)^{a} is bounded above by a constant and C+log⁡(1/δ)≤2log⁡(1/δ)C+\log(1/\delta)\leq 2\log(1/\delta) for δ\delta small enough.

Proof of Theorem 3.4

where Nn′\mathcal{N}^{\prime}_{n} and Dn′\mathcal{D}^{\prime}_{n} respectively denote the numerator and denominator of the expression in (26). Observe that

where An′\mathcal{A}^{\prime}_{n} is a subset of σ(y(n))\sigma(y^{(n)}) as in Lemma 5.2 of (replacing θ\theta by θS0c\theta_{S_{0}^{c}} and θ0\theta_{0} by ) defined as

with Pθ0(Anc)≤e−rn2P_{\theta_{0}}(\mathcal{A}_{n}^{c})\leq e^{-r_{n}^{2}} for some sequence of positive real numbers rnr_{n}. We set rn2=qnr_{n}^{2}=q_{n} here. With this choice, from (27),

References