# An Empirical Bayes Optimal Discovery Procedure Based on Semiparametric Hierarchical Mixture Models

Authors: Hisashi Noma; Shigeyuki Matsui
Journal: Computational and Mathematical Methods in Medicine
 (2013)

---

## Abstract

Multiple testing has been widely adopted for genome-wide studies such as microarray experiments. For effective gene selection in these genome-wide studies, the optimal discovery procedure (ODP), which maximizes the number of expected true positives for each fixed number of expected false positives, was developed as a multiple testing extension of the most powerful test for a single hypothesis by Storey (Journal of the Royal Statistical Society, Series B, vol. 69, no. 3, pp. 347–368, 2007). In this paper, we develop an empirical Bayes method for implementing the ODP based on a semiparametric hierarchical mixture model using the “smoothing-by-roughening" approach. Under the semiparametric hierarchical mixture model, (i) the prior distribution can be modeled flexibly, (ii) the ODP test statistic and the posterior distribution are analytically tractable, and (iii) computations are easy to implement. In addition, we provide a significance rule based on the false discovery rate (FDR) in the empirical Bayes framework. Applications to two clinical studies are presented.


---

## Body

## 1. Introduction

The comprehensive, gene expression microarrays are a powerful tool for screening the differentially expressed genes among different phenotypes such as clinical subtypes and prognostic classes of disease from a large pool of candidate genes. Gene screening studies have the potential to be useful for elucidating disease biology and aggressiveness, identifying new therapeutic targets, and developing new molecular diagnostics for optimized medicine for individual patients [1–4]. The high dimensionality of microarray data, however, has posed special challenges in extracting a small number of relevant genes from a large quantity of noise variables by gene screening analysis. In addition to the control of false positives [5–7], improvement of the efficiency of gene screening analysis is important.For efficient screening of differentially expressed genes, Storey [8] developed the optimal discovery procedure (ODP), which can be interpreted as an extension of the most powerful test for a single hypothesis testing [9] to multiple testing. Storey defined an optimality criterion for multiple testing that maximizes the expected number of true positives (ETP) for each fixed level of expected false positives (EFP) [8]. The ODP was developed as a testing procedure that achieves this optimality; it improves the overall performance of multiple tests by “borrowing strength” across tests. However, in applying the ODP, the following two components must be estimated: (a) the true status of each significance test (null or alternative) and (b) the true probability distribution corresponding to each test [8]. To address these estimation problems, Storey et al. [10] constructed the generalized likelihood ratio statistic [11], using maximum likelihood estimates. Cao et al. [12] also proposed a Bayesian approach, and Woo et al. [13] developed a computationally efficient method called “modular ODP.”More recently, Noma and Matsui [14] developed an empirical Bayes approach for the ODP that can effectively circumvent the estimation problems of (a) and (b). Under a hierarchical, random effects model, they showed that the ODP was derived as a testing rule based on the marginal likelihood ratio statistic. In their numerical evaluations based on simulations, the empirical Bayes method nearly achieved the theoretical bound of ETP (i.e., average power) and performed well, compared with the existing methods under a broad range of conditions [14].A critical aspect of this empirical Bayes approach is that a fully parametric form in the hierarchical modelling must be specified. Although most of the empirical Bayes methods previously discussed in microarray studies [15–21] assume fully parametric natural conjugate models for the prior distribution, this assumption might be restrictive because of possible substantial diversity of effect size among informative genes caused by complex molecular pathways, whose distribution is generally unknown in microarray studies. In this article, we develop an empirical Bayes method for the ODP based on a semiparametric hierarchical mixture model, in which the prior distribution for effect sizes is not specified parametrically. We estimate the nonparametric component of the prior distribution by applying the “smoothing-by-roughening” approach [22, 23] and provide an effective computational method for the ODP statistic.This article is organized as follows. In Section2, we provide a brief overview of the theoretical results regarding the ODP methods by Storey [8] and Noma and Matsui [14]. In Section 3, we describe the semiparametric hierarchical mixture model and how to implement the ODP method. We present applications to prostate cancer and lymphoma clinical studies in Section 4. Finally, we provide a discussion in Section 5.


## 2. The Optimal Discovery Procedure in Multiple Significance Testing

### 2.1. Storey’s Optimal Discovery Procedure

In this section, we briefly review the theoretical results regarding the ODP methods. We denote the datasets fromm genes as x1,x2,…,xm, where xk=(xk1,xk2,…,xkn)T(k=1,…,m). In microarray studies, xk corresponds to a set of log-transformed, normalized gene expressions from n samples for the kth gene (see Section 4). Also, suppose that the datasets are random vectors defined in a common probability space. The optimal goal when identifying differentially expressed genes is to find as many true positives as possible, without incurring too many false positives [10]. Storey [8] defined an optimal criterion of multiple testing as a rule of statistical significance that maximizes ETP for a certain fixed EFP level.He derived the following lemma that gives the multiple testing procedure that achieves the optimal criterion.Lemma 1 (Storey [8]).
Suppose thekth test has a null density fk(x) and alternative density gk(x)(k=1,…,m). Also, assume that the null hypothesis is true for the kth test (k=1,…,m0) and the alternative is true for k=m0+1,…,m, without loss of generality. The following significance threshold statistic achieves the ODP criterion:
(1)SODP(x)=gm0+1(x)+gm0+2(x)+⋯+gm(x)f1(x)+f2(x)+⋯+fm0(x).
Given a fixed cut-off λ(0≤λ<∞) chosen to attain an acceptable EFP level, the null hypothesis for the kth test is rejected if and only if SODP(xk)≥λ.It should be noted that the ODP statistic is composed of all the density functions onm tests. In other words, through borrowing strengths across m tests, the overall power is improved more than the conventional most powerful test for single hypothesis tests [9, 11]. Also, the optimality result is held regardless of correlation structure among the datasets x1,x2,…,xm.As mentioned in Section1, for applying SODP(x), there are two practical issues: (a) estimation of the true status of each hypothesis (null or alternative), corresponding to the numerator and denominator of SODP(x), and (b) estimation of the true probability density functions fk(x)’s and gk(x)’s at which xk’s are evaluated. Storey et al. [10] provided an practical estimating method, motivated by the generalized likelihood ratio test for single significance tests [11]. They noticed that the significance rule based on
(2)SODP(x)+1=f1(x)+⋯+fm0(x)+gm0+1(x)+⋯+gm(x)f1(x)+f2(x)+⋯+fm0(x)
is equivalent to that of SODP(x) in Lemma 1. For each test k=1,…,m, let f^k(x) be an estimate of the null density function with all the unknown parameters replaced by their maximum likelihood estimates under the constraints of the null hypothesis, and let g^k(x) be the analogous estimate for the alternative hypothesis using the unconstrained maximum likelihood estimates. Then, Storey et al. [10] proposed to use the canonical plug-in estimate:
(3)S^ODP(x)=g^1(x)+g^2(x)+⋯+g^m(x)w^1f^1(x)+w^2f^2(x)+⋯+w^mf^m(x),
where w^k(k=1,2,…,m) is an estimate of the status of the kth hypothesis. For estimating w^k, Storey et al. [10] ranked the tests based on a univariate statistic (e.g., the t-statistic for two-class comparison) and estimated it as w^k=1 when the kth test fell below the estimated proportion of true nulls by Storey [6], and w^k=0 otherwise. Other estimating procedures for SODP(x) were also proposed by Cao et al. [12] and Woo et al. [13].


### 2.2. Empirical Bayes Approach by Noma and Matsui [14]

As Storey [8] pointed out, the ODP is analogous to the well-known shrinkage estimation for multiple-point estimation [24, 25]. Because the shrinkage estimation is interpreted as the empirical Bayes estimation [26–28], hierarchical modelling would be a natural formulation for information sharing across tests, both in deriving the ODP and in developing estimation methods for implementing it [14].Suppose that the probability density functions forx1,x2,…,xm have the same parametric form f(x∣θk,ψk), where θk is the parameter of interest and ψk is the nuisance parameter for the kth test. We assume the following prior distribution for the parameters (θk,ψk):
(4)null:(θk,ψk)~G0(θ,ψ∣ξ0),alternative:(θk,ψk)~G1(θ,ψ∣ξ1),
where ξ0 and ξ1 are the hyperparameters of the prior distribution.Lemma 2 (Noma and Matsui [14]).
Under the empirical Bayes framework with the prior distribution (4), the following significance threshold function achieves the ODP criterion:
(5)RODP(x)=EG1[f(x∣θ,ψ)]EG0[f(x∣θ,ψ)]=∫f(x∣θ,ψ)dG1(θ,ψ∣ξ1)∫f(x∣θ,ψ)dG0(θ,ψ∣ξ0).
For a fixed cut-off λ(0≤λ<∞) chosen to reach an acceptable EFP level, the null hypothesis for the kth test is rejected if and only if RODP(xk)≥λ.See Noma and Matsui [14] for interpretations of RODP(x).For the specification of the prior distributionsG0(θ,ψ∣ξ0) and G1(θ,ψ∣ξ1), Noma and Matsui [14] adopted the empirical Bayes method [26, 27]. The empirical Bayes implementation is used to obtain estimates for ξ0 and ξ1 from the data, which can then be substituted into the ODP statistic, RODP(x),
(6)R^ODP(x)=∫f(x∣θ,ψ)dG1(θ,ψ∣ξ^1)∫f(x∣θ,ψ)dG0(θ,ψ∣ξ^0),
where ξ^0 and ξ^1 are the estimators of ξ0 and ξ1, respectively.In comparison with the existing ODP estimation methods based onSODP(x), the empirical Bayes approach has two advantages. First, it estimates only the hyperparameters ξ0 and ξ1 for the random effects models, and not the parameters in the density function xk(k=1,…,m), for which a large number of parameters (thousands or more in microarray studies) must be estimated as in the procedure of Storey et al. [10]. Second, in the empirical Bayes approach, the estimation of the true status, null or alternative, for each test is not needed.


## 2.1. Storey’s Optimal Discovery Procedure

In this section, we briefly review the theoretical results regarding the ODP methods. We denote the datasets fromm genes as x1,x2,…,xm, where xk=(xk1,xk2,…,xkn)T(k=1,…,m). In microarray studies, xk corresponds to a set of log-transformed, normalized gene expressions from n samples for the kth gene (see Section 4). Also, suppose that the datasets are random vectors defined in a common probability space. The optimal goal when identifying differentially expressed genes is to find as many true positives as possible, without incurring too many false positives [10]. Storey [8] defined an optimal criterion of multiple testing as a rule of statistical significance that maximizes ETP for a certain fixed EFP level.He derived the following lemma that gives the multiple testing procedure that achieves the optimal criterion.Lemma 1 (Storey [8]).
Suppose thekth test has a null density fk(x) and alternative density gk(x)(k=1,…,m). Also, assume that the null hypothesis is true for the kth test (k=1,…,m0) and the alternative is true for k=m0+1,…,m, without loss of generality. The following significance threshold statistic achieves the ODP criterion:
(1)SODP(x)=gm0+1(x)+gm0+2(x)+⋯+gm(x)f1(x)+f2(x)+⋯+fm0(x).
Given a fixed cut-off λ(0≤λ<∞) chosen to attain an acceptable EFP level, the null hypothesis for the kth test is rejected if and only if SODP(xk)≥λ.It should be noted that the ODP statistic is composed of all the density functions onm tests. In other words, through borrowing strengths across m tests, the overall power is improved more than the conventional most powerful test for single hypothesis tests [9, 11]. Also, the optimality result is held regardless of correlation structure among the datasets x1,x2,…,xm.As mentioned in Section1, for applying SODP(x), there are two practical issues: (a) estimation of the true status of each hypothesis (null or alternative), corresponding to the numerator and denominator of SODP(x), and (b) estimation of the true probability density functions fk(x)’s and gk(x)’s at which xk’s are evaluated. Storey et al. [10] provided an practical estimating method, motivated by the generalized likelihood ratio test for single significance tests [11]. They noticed that the significance rule based on
(2)SODP(x)+1=f1(x)+⋯+fm0(x)+gm0+1(x)+⋯+gm(x)f1(x)+f2(x)+⋯+fm0(x)
is equivalent to that of SODP(x) in Lemma 1. For each test k=1,…,m, let f^k(x) be an estimate of the null density function with all the unknown parameters replaced by their maximum likelihood estimates under the constraints of the null hypothesis, and let g^k(x) be the analogous estimate for the alternative hypothesis using the unconstrained maximum likelihood estimates. Then, Storey et al. [10] proposed to use the canonical plug-in estimate:
(3)S^ODP(x)=g^1(x)+g^2(x)+⋯+g^m(x)w^1f^1(x)+w^2f^2(x)+⋯+w^mf^m(x),
where w^k(k=1,2,…,m) is an estimate of the status of the kth hypothesis. For estimating w^k, Storey et al. [10] ranked the tests based on a univariate statistic (e.g., the t-statistic for two-class comparison) and estimated it as w^k=1 when the kth test fell below the estimated proportion of true nulls by Storey [6], and w^k=0 otherwise. Other estimating procedures for SODP(x) were also proposed by Cao et al. [12] and Woo et al. [13].


## 2.2. Empirical Bayes Approach by Noma and Matsui [14]

As Storey [8] pointed out, the ODP is analogous to the well-known shrinkage estimation for multiple-point estimation [24, 25]. Because the shrinkage estimation is interpreted as the empirical Bayes estimation [26–28], hierarchical modelling would be a natural formulation for information sharing across tests, both in deriving the ODP and in developing estimation methods for implementing it [14].Suppose that the probability density functions forx1,x2,…,xm have the same parametric form f(x∣θk,ψk), where θk is the parameter of interest and ψk is the nuisance parameter for the kth test. We assume the following prior distribution for the parameters (θk,ψk):
(4)null:(θk,ψk)~G0(θ,ψ∣ξ0),alternative:(θk,ψk)~G1(θ,ψ∣ξ1),
where ξ0 and ξ1 are the hyperparameters of the prior distribution.Lemma 2 (Noma and Matsui [14]).
Under the empirical Bayes framework with the prior distribution (4), the following significance threshold function achieves the ODP criterion:
(5)RODP(x)=EG1[f(x∣θ,ψ)]EG0[f(x∣θ,ψ)]=∫f(x∣θ,ψ)dG1(θ,ψ∣ξ1)∫f(x∣θ,ψ)dG0(θ,ψ∣ξ0).
For a fixed cut-off λ(0≤λ<∞) chosen to reach an acceptable EFP level, the null hypothesis for the kth test is rejected if and only if RODP(xk)≥λ.See Noma and Matsui [14] for interpretations of RODP(x).For the specification of the prior distributionsG0(θ,ψ∣ξ0) and G1(θ,ψ∣ξ1), Noma and Matsui [14] adopted the empirical Bayes method [26, 27]. The empirical Bayes implementation is used to obtain estimates for ξ0 and ξ1 from the data, which can then be substituted into the ODP statistic, RODP(x),
(6)R^ODP(x)=∫f(x∣θ,ψ)dG1(θ,ψ∣ξ^1)∫f(x∣θ,ψ)dG0(θ,ψ∣ξ^0),
where ξ^0 and ξ^1 are the estimators of ξ0 and ξ1, respectively.In comparison with the existing ODP estimation methods based onSODP(x), the empirical Bayes approach has two advantages. First, it estimates only the hyperparameters ξ0 and ξ1 for the random effects models, and not the parameters in the density function xk(k=1,…,m), for which a large number of parameters (thousands or more in microarray studies) must be estimated as in the procedure of Storey et al. [10]. Second, in the empirical Bayes approach, the estimation of the true status, null or alternative, for each test is not needed.


## 3. Semiparametric Hierarchical Mixture Model and the ODP

### 3.1. Semiparametric Hierarchical Mixture Model

In this section, we consider applying the ODP to concrete microarray analyses of two-class comparisons (with classes denoted as “0” and “1”), for example, good prognosis and poor prognosis, on the basis of the expression levels ofm genes from n samples. The gene expression data considered here comprise normalized log ratios from two-color complementary DNA arrays or normalized log signals from oligonucleotide arrays (e.g., Affymetrix GeneChip). Let n0 and n1 be the sample sizes of classes 0 and 1, respectively. Let μ0k and μ1k(k=1,2,…,m) be the means of the gene expression levels for classes 0 and 1, respectively. We consider the detection of differentially expressed genes to be when μ0k≠μ1k. That is, testing problem for
(7)H0k:μ0k=μ1kversusH1k:μ0k≠μ1kaaaaaaaaaaaaaaaaa(k=1,2,…,m).We assume that gene expressions of thekth gene are normally distributed within each class,
(8)xk1,…,xkn0~N(μ0k,σk2),yk1,…,ykn1~N(μ1k,σk2),
where σk2 is a common within-class variance (k=1,2,…,m), as in Storey et al. [10]. Let dk be the difference of the mean of the two classes for the kth gene (k=1,2,…,m), μ1k-μ0k, which is the parameter of interest. Let x-k, y-k denote the sample means of the gene expression levels for classes 0 and 1, respectively. We consider the following translation:
(9)wki=n0n0-1(xki-x-k),zki=n0n0+1(yki-x-k).
Using this translation, the nuisance parameters μ0k’s are eliminated in the models of wki and zki. Specifically,
(10)wk1,…,wkn0~N(0,σk2),zk1,…,zkn1~N(θk,σk2),
where θk=n0/(n0+1)dk. Note that the testing problems in (7) are equivalent to testing
(11)H0k’:θk=0versusH1k’:θk≠0aaaaaaaaaaa(k=1,2,…,m).
Let wk=(wk1,…,wkn0) and zk=(zk1,…,zkn1).For applying the ODP to the testing problemsH0k’:θk=0(k=1,2,…,m) under model (10), we assume the following random effects distributions for the model parameters (θk,σk2):
(12)θk~π0δ(θ)+π1G(θ),σk2~Gamma-1(α,β),[p(σk2)=βαΓ(α)σk-2(α+1)exp(-βσk2)],
where δ(θ) in the first component of θk is the Dirac delta function, representing nondifferential expression between two classes. The second component of θk represents a component of the differential expression, and π0 and π1 are the mixing proportions for the components (π0+π1=1). The random effects distribution of the gene-specific variances σk2 is specified as the parametric conjugate model for analytical tractability. For modeling G(θ), instead of adopting a fully parametric conjugate normal distribution model as done by Noma and Matsui [14], we apply the “smoothing-by-roughening” approach of Laird and Louis [22] and Shen and Louis [23]. This approach was proposed to obtain the nonparametric maximum likelihood estimate of the prior distribution in empirical Bayes analysis [29, 30]. This approach starts with a smooth estimate of the prior distribution and uses the EM algorithm to “roughen” it towards the nonparametric maximum likelihood estimate. For modeling G(θ), we place supports at equally spaced discrete mass points at a series of nonzero points a=(a1,a2,…,aL), for sufficiently large values of L. We assume that
(13)Pr(θ=aj)=pj(j=1,2,…,L),
where p1+p2+⋯+pL=1. We denote the hyperparameters of the distribution function for the differential component as p=(p1,p2,…,pL). Shen and Louis [23] provided a detailed guideline for implementing the smoothing-by-roughening approach based on their numerical evaluations. Shen and Louis [23] recommended that the number of grid points should be at least 200. For simplicity, they suggested the grid points can be equally spaced on the support of the distribution. Increasing the number of grid points and spacing them according to the prior density will provide a better approximation of G(θ) and therefore will improve the estimation, particularly for the tail areas. However, too large numbers might not necessarily be beneficial. For more details involving analytical and numerical evaluations, see the guideline in Shen and Louis [23]. In this article, we refer to this hierarchical model as the “semiparametric hierarchical mixture model.” The important property of this model is that the posterior distributions of θk and σk2(k=1,2,…,m) can be obtained analytically. We provide the analytical representation of the posterior distributions in Section 3.4.The maximum likelihood estimates of the hyperparametersη=(π0,π1,p,α,β) are obtained by the EM algorithm [31]. Because the posterior distributions of θk and σk2 are obtained analytically, the E-step of the EM algorithm can be implemented easily. The details of the EM algorithm are presented in the appendix.


### 3.2. The Optimal Discovery Procedure

The ODP statisticRODP(w,z) for the testing problems (11) is obtained as the marginal likelihood ratio. The marginal likelihood of the null hypothesis corresponds to the marginal likelihood of the nondifferential component:
(14)h0(wk,zk)=∫0∞p(wk,zk∣θk=0,σk2)p(σk2∣α,β)dσk2=Γ(n/2+α)Γ(α)βα(2π)-n/2×(∑i=1n0wki2+∑i=1n1zki22+β)-(n/2+α).
Similarly, the marginal likelihood of the alternative hypothesis in (11) corresponds to that of the differential component:
(15)h1(wk,zk)=∑j=1L∫0∞p(wk,zk∣θk=aj,σk2)×p(σk2∣α,β)dσk2Pr(θk=aj)=Γ(n/2+α)Γ(α)βα(2π)-n/2×∑j=1L(∑i=1n0wki2+∑i=1n1(zki-aj)22+β)-(n/2+α)pj
under the semiparametric hierarchical mixture model.Thus, the ODP statistic (the marginal likelihood ratio) is given as(16)RODP(wk,zk)=(∑j=1L∫0∞p(wk,zk∣θk=aj,σk2)×p(σk2∣α,β)dσk2Pr(θk=aj)∑j=1L∫0∞p(wk,zk∣θk=aj,σk2))×(∫0∞p(wk,zk∣θk=0,σk2)p(σk2∣α,β)dσk2)-1=(∑j=1L(∑i=1n0wki2+∑i=1n1(zki-aj)22+β)-(n/2+α)pj)×((∑i=1n0wki2+∑i=1n1zki22+β)-(n/2+α))-1.
By plugging-in the maximum likelihood estimate of η, the empirical Bayes testing procedure is implemented based on
(17)R^ODP(wk,zk)=∑j=1L((∑i=1n0wki2+∑i=1n1(zki-aj)2)/2+β^)-(n/2+α^)p^j((∑i=1n0wki2+∑i=1n1zki2)/2+β^)-(n/2+α^).
Note that the ODP statistic is obtained in closed form using the maximum likelihood estimate of η.


### 3.3. Assessing Significance

For assessing significance, we can obtain an empirical Bayes estimator of FDR [5–7] under the semiparametric hierarchical mixture model as follows. For any given threshold λ, we denote the significant region Λ(λ)={(w,z)|RODP(w,z)≥λ} for each test. Let Ξ denote the set of indices for the significant genes. In the framework of Bayesian selection rules, the FDR accords to the misclassification error rate of differential/non-differential classification [6, 7, 32], and one of the well-known Bayesian estimator of the FDR [33, 34] is
(18)FDR^(λ)=1|Ξ|∑k∈ΞPr(H0k∣wk,zk),
where
(19)Pr(H0k∣wk,zk)=π0h0(wk,zk)π0h0(wk,zk)+π1h1(wk,zk)=(π0(∑i=1n0wki2+∑i=1n1zki22+β)-(n/2+α))×((∑i=1n0wki2+∑i=1n1(zki-aj)22+β)-(n/2+α)π0(∑i=1n0wki2+∑i=1n1zki22+β)-(n/2+α)+π1∑j=1L(∑i=1n0wki2+∑i=1n1(zki-aj)22+β)-(n/2+α)pj)-1.
Accordingly, an empirical Bayes estimator of FDR^(λ) for the significant region Λ(λ)={(w,z)∣RODP(w,z)≥λ} can be obtained by plugging-in the maximum likelihood estimate of η.Note that we conducted small simulations for checking correctness of the proposed ODP procedure and these results were provided in the Supplementary Materials available online athttp://dx.doi.org/10.1155/2013/568480.


### 3.4. The Posterior Distribution

Under the semiparametric hierarchical mixture model, the posterior distributions ofθk and σk2(k=1,2,…,m) can be obtained analytically as noted above; this gives a computational advantage for implementing Bayesian inference. The posterior distribution will be used to derive estimates and confidence intervals of dk’s, in addition to the statistical significance measures such as the q-value.Letγ0k and γ1k denote indicator variables to which the kth gene component belongs (k=1,2,…,m). The joint posterior distribution of (θk,σk2) is expressed as a mixture of differential and non-differential components:
(20)p(θk,σk2∣wk,zk)=p(θk,σk2∣wk,zk,γ0k=1)Pr(γ0k=1∣wk,zk)+p(θk,σk2∣wk,zk,γ1k=1)Pr(γ1k=1∣wk,zk),
where p(θk,σk2∣wk,zk,γ0k=1)=δ(θ)p(σk2∣wk,zk). The weight Pr(γ0k=1∣wk,zk)=Pr(H0k∣wk,zk) is given in Section 3.3, and Pr(γ0k=1∣wk,zk)=1-Pr(γ1k=1∣wk,zk). The posterior distribution of the differential component is expressed as
(21)p(θk=aj,σk2∣wk,zk,γ1k=1)=Pr(θk=aj∣σk2,wk,zk)p(σk2∣wk,zk),
where
(22)Pr(θk=aj∣σk2,wk,zk)=p(wk,zk∣θk=aj,σk2)pj∑l=1Lp(wk,zk∣θk=al,σk2)pl,(23)p(σk2∣wk,zk)=p(wk,zk∣σk2)p(σk2)∫0∞p(wk,zk∣σk2)p(σk2)dσk2=∑j=1Lp(wk,zk∣θk=aj,σk2)p(σk2)pj∑j=1L∫0∞p(wk,zk∣θk=aj,σk2)p(σk2)dσk2pj,(24)∫0∞p(wk,zk∣θk=aj,σk2)p(σk2)dσk2=Γ(n/2+α)Γ(α)βα(2π)-n/2×(∑i=1n0wki2+∑i=1n1(zki2-aj)22+β)-(n/2+α).
The marginal posterior distribution of θk can be obtained as follows:
(25)Pr(θk=aj∣wk,zk)=p(wk,zk∣θk=aj)Pr(θk=aj)∑l=1Lp(wk,zk∣θk=al)Pr(θk=al)=((∑i=1n0wki2+∑i=1n1(zki2-aj)2)/2+β)-(n/2+α)pj∑l=1L{((∑i=1n0wki2+∑i=1n1(zki2-al)2)/2+β)-(n/2+α)pl}.


## 3.1. Semiparametric Hierarchical Mixture Model

In this section, we consider applying the ODP to concrete microarray analyses of two-class comparisons (with classes denoted as “0” and “1”), for example, good prognosis and poor prognosis, on the basis of the expression levels ofm genes from n samples. The gene expression data considered here comprise normalized log ratios from two-color complementary DNA arrays or normalized log signals from oligonucleotide arrays (e.g., Affymetrix GeneChip). Let n0 and n1 be the sample sizes of classes 0 and 1, respectively. Let μ0k and μ1k(k=1,2,…,m) be the means of the gene expression levels for classes 0 and 1, respectively. We consider the detection of differentially expressed genes to be when μ0k≠μ1k. That is, testing problem for
(7)H0k:μ0k=μ1kversusH1k:μ0k≠μ1kaaaaaaaaaaaaaaaaa(k=1,2,…,m).We assume that gene expressions of thekth gene are normally distributed within each class,
(8)xk1,…,xkn0~N(μ0k,σk2),yk1,…,ykn1~N(μ1k,σk2),
where σk2 is a common within-class variance (k=1,2,…,m), as in Storey et al. [10]. Let dk be the difference of the mean of the two classes for the kth gene (k=1,2,…,m), μ1k-μ0k, which is the parameter of interest. Let x-k, y-k denote the sample means of the gene expression levels for classes 0 and 1, respectively. We consider the following translation:
(9)wki=n0n0-1(xki-x-k),zki=n0n0+1(yki-x-k).
Using this translation, the nuisance parameters μ0k’s are eliminated in the models of wki and zki. Specifically,
(10)wk1,…,wkn0~N(0,σk2),zk1,…,zkn1~N(θk,σk2),
where θk=n0/(n0+1)dk. Note that the testing problems in (7) are equivalent to testing
(11)H0k’:θk=0versusH1k’:θk≠0aaaaaaaaaaa(k=1,2,…,m).
Let wk=(wk1,…,wkn0) and zk=(zk1,…,zkn1).For applying the ODP to the testing problemsH0k’:θk=0(k=1,2,…,m) under model (10), we assume the following random effects distributions for the model parameters (θk,σk2):
(12)θk~π0δ(θ)+π1G(θ),σk2~Gamma-1(α,β),[p(σk2)=βαΓ(α)σk-2(α+1)exp(-βσk2)],
where δ(θ) in the first component of θk is the Dirac delta function, representing nondifferential expression between two classes. The second component of θk represents a component of the differential expression, and π0 and π1 are the mixing proportions for the components (π0+π1=1). The random effects distribution of the gene-specific variances σk2 is specified as the parametric conjugate model for analytical tractability. For modeling G(θ), instead of adopting a fully parametric conjugate normal distribution model as done by Noma and Matsui [14], we apply the “smoothing-by-roughening” approach of Laird and Louis [22] and Shen and Louis [23]. This approach was proposed to obtain the nonparametric maximum likelihood estimate of the prior distribution in empirical Bayes analysis [29, 30]. This approach starts with a smooth estimate of the prior distribution and uses the EM algorithm to “roughen” it towards the nonparametric maximum likelihood estimate. For modeling G(θ), we place supports at equally spaced discrete mass points at a series of nonzero points a=(a1,a2,…,aL), for sufficiently large values of L. We assume that
(13)Pr(θ=aj)=pj(j=1,2,…,L),
where p1+p2+⋯+pL=1. We denote the hyperparameters of the distribution function for the differential component as p=(p1,p2,…,pL). Shen and Louis [23] provided a detailed guideline for implementing the smoothing-by-roughening approach based on their numerical evaluations. Shen and Louis [23] recommended that the number of grid points should be at least 200. For simplicity, they suggested the grid points can be equally spaced on the support of the distribution. Increasing the number of grid points and spacing them according to the prior density will provide a better approximation of G(θ) and therefore will improve the estimation, particularly for the tail areas. However, too large numbers might not necessarily be beneficial. For more details involving analytical and numerical evaluations, see the guideline in Shen and Louis [23]. In this article, we refer to this hierarchical model as the “semiparametric hierarchical mixture model.” The important property of this model is that the posterior distributions of θk and σk2(k=1,2,…,m) can be obtained analytically. We provide the analytical representation of the posterior distributions in Section 3.4.The maximum likelihood estimates of the hyperparametersη=(π0,π1,p,α,β) are obtained by the EM algorithm [31]. Because the posterior distributions of θk and σk2 are obtained analytically, the E-step of the EM algorithm can be implemented easily. The details of the EM algorithm are presented in the appendix.


## 3.2. The Optimal Discovery Procedure

The ODP statisticRODP(w,z) for the testing problems (11) is obtained as the marginal likelihood ratio. The marginal likelihood of the null hypothesis corresponds to the marginal likelihood of the nondifferential component:
(14)h0(wk,zk)=∫0∞p(wk,zk∣θk=0,σk2)p(σk2∣α,β)dσk2=Γ(n/2+α)Γ(α)βα(2π)-n/2×(∑i=1n0wki2+∑i=1n1zki22+β)-(n/2+α).
Similarly, the marginal likelihood of the alternative hypothesis in (11) corresponds to that of the differential component:
(15)h1(wk,zk)=∑j=1L∫0∞p(wk,zk∣θk=aj,σk2)×p(σk2∣α,β)dσk2Pr(θk=aj)=Γ(n/2+α)Γ(α)βα(2π)-n/2×∑j=1L(∑i=1n0wki2+∑i=1n1(zki-aj)22+β)-(n/2+α)pj
under the semiparametric hierarchical mixture model.Thus, the ODP statistic (the marginal likelihood ratio) is given as(16)RODP(wk,zk)=(∑j=1L∫0∞p(wk,zk∣θk=aj,σk2)×p(σk2∣α,β)dσk2Pr(θk=aj)∑j=1L∫0∞p(wk,zk∣θk=aj,σk2))×(∫0∞p(wk,zk∣θk=0,σk2)p(σk2∣α,β)dσk2)-1=(∑j=1L(∑i=1n0wki2+∑i=1n1(zki-aj)22+β)-(n/2+α)pj)×((∑i=1n0wki2+∑i=1n1zki22+β)-(n/2+α))-1.
By plugging-in the maximum likelihood estimate of η, the empirical Bayes testing procedure is implemented based on
(17)R^ODP(wk,zk)=∑j=1L((∑i=1n0wki2+∑i=1n1(zki-aj)2)/2+β^)-(n/2+α^)p^j((∑i=1n0wki2+∑i=1n1zki2)/2+β^)-(n/2+α^).
Note that the ODP statistic is obtained in closed form using the maximum likelihood estimate of η.


## 3.3. Assessing Significance

For assessing significance, we can obtain an empirical Bayes estimator of FDR [5–7] under the semiparametric hierarchical mixture model as follows. For any given threshold λ, we denote the significant region Λ(λ)={(w,z)|RODP(w,z)≥λ} for each test. Let Ξ denote the set of indices for the significant genes. In the framework of Bayesian selection rules, the FDR accords to the misclassification error rate of differential/non-differential classification [6, 7, 32], and one of the well-known Bayesian estimator of the FDR [33, 34] is
(18)FDR^(λ)=1|Ξ|∑k∈ΞPr(H0k∣wk,zk),
where
(19)Pr(H0k∣wk,zk)=π0h0(wk,zk)π0h0(wk,zk)+π1h1(wk,zk)=(π0(∑i=1n0wki2+∑i=1n1zki22+β)-(n/2+α))×((∑i=1n0wki2+∑i=1n1(zki-aj)22+β)-(n/2+α)π0(∑i=1n0wki2+∑i=1n1zki22+β)-(n/2+α)+π1∑j=1L(∑i=1n0wki2+∑i=1n1(zki-aj)22+β)-(n/2+α)pj)-1.
Accordingly, an empirical Bayes estimator of FDR^(λ) for the significant region Λ(λ)={(w,z)∣RODP(w,z)≥λ} can be obtained by plugging-in the maximum likelihood estimate of η.Note that we conducted small simulations for checking correctness of the proposed ODP procedure and these results were provided in the Supplementary Materials available online athttp://dx.doi.org/10.1155/2013/568480.


## 3.4. The Posterior Distribution

Under the semiparametric hierarchical mixture model, the posterior distributions ofθk and σk2(k=1,2,…,m) can be obtained analytically as noted above; this gives a computational advantage for implementing Bayesian inference. The posterior distribution will be used to derive estimates and confidence intervals of dk’s, in addition to the statistical significance measures such as the q-value.Letγ0k and γ1k denote indicator variables to which the kth gene component belongs (k=1,2,…,m). The joint posterior distribution of (θk,σk2) is expressed as a mixture of differential and non-differential components:
(20)p(θk,σk2∣wk,zk)=p(θk,σk2∣wk,zk,γ0k=1)Pr(γ0k=1∣wk,zk)+p(θk,σk2∣wk,zk,γ1k=1)Pr(γ1k=1∣wk,zk),
where p(θk,σk2∣wk,zk,γ0k=1)=δ(θ)p(σk2∣wk,zk). The weight Pr(γ0k=1∣wk,zk)=Pr(H0k∣wk,zk) is given in Section 3.3, and Pr(γ0k=1∣wk,zk)=1-Pr(γ1k=1∣wk,zk). The posterior distribution of the differential component is expressed as
(21)p(θk=aj,σk2∣wk,zk,γ1k=1)=Pr(θk=aj∣σk2,wk,zk)p(σk2∣wk,zk),
where
(22)Pr(θk=aj∣σk2,wk,zk)=p(wk,zk∣θk=aj,σk2)pj∑l=1Lp(wk,zk∣θk=al,σk2)pl,(23)p(σk2∣wk,zk)=p(wk,zk∣σk2)p(σk2)∫0∞p(wk,zk∣σk2)p(σk2)dσk2=∑j=1Lp(wk,zk∣θk=aj,σk2)p(σk2)pj∑j=1L∫0∞p(wk,zk∣θk=aj,σk2)p(σk2)dσk2pj,(24)∫0∞p(wk,zk∣θk=aj,σk2)p(σk2)dσk2=Γ(n/2+α)Γ(α)βα(2π)-n/2×(∑i=1n0wki2+∑i=1n1(zki2-aj)22+β)-(n/2+α).
The marginal posterior distribution of θk can be obtained as follows:
(25)Pr(θk=aj∣wk,zk)=p(wk,zk∣θk=aj)Pr(θk=aj)∑l=1Lp(wk,zk∣θk=al)Pr(θk=al)=((∑i=1n0wki2+∑i=1n1(zki2-aj)2)/2+β)-(n/2+α)pj∑l=1L{((∑i=1n0wki2+∑i=1n1(zki2-al)2)/2+β)-(n/2+α)pl}.


## 4. Applications

We illustrate our method using real data sets from two clinical studies with microarrays. Because previous simulation studies have indicated that the Storey’s method [10] is more powerful than many other multiple testing methods [10, 12, 35], we chose to use this method as a reference.


### 4.1. Prostate Cancer Example

Efron [36, 37] analyzed a modified gene expression data set from Singh et al. [38] with m=6033 genes for n=102 samples; n1=52 prostate cancer patients and n0=50 healthy controls. For the proposed empirical Bayes approach, we placed the grid points for G(θ) equally spaced on [-1,1] by 0.01 (except for 0; L=200). The hyperparameters were estimated as π^0=0.716,π^1=0.284,α^=8.624,andβ^=7.373, and the estimated prior distribution for G(θ) is given in Figure 1. The estimated prior distribution was skewed and multimodal. This indicates that analytically tractable parametric models (for example, the normal distribution) might be inadequate in this case.Figure 1
The nonparametric estimate for the differential component of the prior distribution of the prostate cancer data [38].Figure2 presents the results of multiple testing by the method of Storey et al. [10] and the empirical Bayes method developed in this article. This plot summarizes the correspondence of the numbers of significant genes and q-values [6, 7] for varying values of the cut-off λ. For a fair comparison, the same FDR-controlling method presented in Section 3.3 was applied to these two methods. For any levels of the q-values, the number of significant genes identified by the proposed method was greater than that of the Storey et al. [10] method.Figure 2
Results of the multiple testing for the prostate cancer data [38].


### 4.2. Lymphoma Example

Dave et al. [39] analyzed Affymetrix HG-U133A and U133B (Affymetrix, Santa Clara, CA) microarray data for 191 biopsy specimens from patients with untreated follicular lymphoma with m=44,928. The data are available from the BRB-ArrayTools Data Archive for Human Cancer Gene Expression [40]. We compared patients who died within 5 years (poor prognosis; n1=51) with patients who survived for more than 5 years (good prognosis; n0=109).For the maximum likelihood estimate of the semiparametric prior distribution, the hyperparameters were estimated asπ^0=0.831,π^1=0.169,α^=1.266,andβ^=0.311. The estimated prior distribution for G(θ) is given in Figure 3. (The grid points were placed on [-1,1] equally spaced by 0.01 except for 0; L=200.) Compared with the prostate cancer data, large masses of the estimated G(θ) are located in the regions of small θ, and thus the signal on differential expression contained in the data would be weaker than that in the prostate cancer data. In Figure 4, the relationships of the numbers of significant genes and their q-values are presented. The number of significant genes identified by the empirical Bayes method was greater than that of the method of Storey et al. [10].Figure 3
The nonparametric estimate for the differential component of the prior distribution of the lymphoma data [39].Figure 4
Results of the multiple testing for the lymphoma data [39].


## 4.1. Prostate Cancer Example

Efron [36, 37] analyzed a modified gene expression data set from Singh et al. [38] with m=6033 genes for n=102 samples; n1=52 prostate cancer patients and n0=50 healthy controls. For the proposed empirical Bayes approach, we placed the grid points for G(θ) equally spaced on [-1,1] by 0.01 (except for 0; L=200). The hyperparameters were estimated as π^0=0.716,π^1=0.284,α^=8.624,andβ^=7.373, and the estimated prior distribution for G(θ) is given in Figure 1. The estimated prior distribution was skewed and multimodal. This indicates that analytically tractable parametric models (for example, the normal distribution) might be inadequate in this case.Figure 1
The nonparametric estimate for the differential component of the prior distribution of the prostate cancer data [38].Figure2 presents the results of multiple testing by the method of Storey et al. [10] and the empirical Bayes method developed in this article. This plot summarizes the correspondence of the numbers of significant genes and q-values [6, 7] for varying values of the cut-off λ. For a fair comparison, the same FDR-controlling method presented in Section 3.3 was applied to these two methods. For any levels of the q-values, the number of significant genes identified by the proposed method was greater than that of the Storey et al. [10] method.Figure 2
Results of the multiple testing for the prostate cancer data [38].


## 4.2. Lymphoma Example

Dave et al. [39] analyzed Affymetrix HG-U133A and U133B (Affymetrix, Santa Clara, CA) microarray data for 191 biopsy specimens from patients with untreated follicular lymphoma with m=44,928. The data are available from the BRB-ArrayTools Data Archive for Human Cancer Gene Expression [40]. We compared patients who died within 5 years (poor prognosis; n1=51) with patients who survived for more than 5 years (good prognosis; n0=109).For the maximum likelihood estimate of the semiparametric prior distribution, the hyperparameters were estimated asπ^0=0.831,π^1=0.169,α^=1.266,andβ^=0.311. The estimated prior distribution for G(θ) is given in Figure 3. (The grid points were placed on [-1,1] equally spaced by 0.01 except for 0; L=200.) Compared with the prostate cancer data, large masses of the estimated G(θ) are located in the regions of small θ, and thus the signal on differential expression contained in the data would be weaker than that in the prostate cancer data. In Figure 4, the relationships of the numbers of significant genes and their q-values are presented. The number of significant genes identified by the empirical Bayes method was greater than that of the method of Storey et al. [10].Figure 3
The nonparametric estimate for the differential component of the prior distribution of the lymphoma data [39].Figure 4
Results of the multiple testing for the lymphoma data [39].


## 5. Discussion

For efficient gene screening using high-dimensional microarray data, multiple testing methodologies are effective tools for controlling false positives findings. Although the multiple testing can provide a relevant framework to ensure control of false positives for a set of significant genes (e.g., FDR [5–7]), researchers would like to find as many true positive genes as possible for further investigations. The optimal discovery procedure (ODP), introduced by Storey [8], would provide an efficient solution for this requirement. Storey et al. [10] provided desirable results in many numerical evaluations [10, 12, 35].In this article, we proposed an empirical Bayes ODP method based on a semiparametric hierarchical model, through relaxing the parametric prior assumption in the empirical Bayes ODP method by Noma and Matsui [14]. Hierarchical mixture models and empirical Bayes methods form a basis for information sharing across genes and provide a framework for efficient multiple testing [15–21]. Although empirical Bayes methods are more robust to prior misspecifications than eliciting a single prior, parametric empirical Bayes methods can suffer from a lack of robustness when the true prior is not from the assumed parametric family, as pointed out by Morris [28], such as a failure to incorporate the mixture structure or specify the form of the effect-size distribution (G). Actually, the estimated prior distribution in the prostate cancer example (Section 4.1) had a multimodal and skewed shape. In such examples, natural conjugate families might be inadequate for modeling an unknown random effects distribution. The “smoothing-by-roughening” approach [22, 23] is one of the flexible modeling methods in empirical Bayes inference [41–43] and could be a basis for effective empirical Bayes inference for a broad range of situations. An extension of our method is to specify a nonparametric distribution for the random effect distribution of σk2 through applying the smoothing-by-roughening approach, although it requires substantial computation for estimating a large number of parameters and can suffer from unstable parameter estimates.Previously, many efficient gene selection methods have been proposed for microarray studies, not only the multiple testing methodologies, for example, Bayesian ranking methods [20, 21, 43]. For evaluating practical values of these methods, comprehensive numerical investigations, including large-scale simulations and applications to many real datasets, would be a worthwhile subject.



---

*Source: 568480-2013-04-10.xml*