A DISCRETE MODIFICATION OF THE
BENJAMINI-YEKUTIELI PROCEDURE

arXiv:1612.04838v1 [] 14 Dec 2016

SEBASTIAN DÖHLER
DARMSTADT UNIVERSITY OF APPLIED SCIENCES

Abstract. The Benjamini-Yekutieli procedure is a multiple testing method
that controls the false discovery rate under arbitrary dependence of the pvalues. A modification of this and related procedures is proposed for the case
when the test statistics are discrete. It is shown that taking discreteness into
account can improve upon known procedures. The performance of this new
procedure is evaluated for pharmacovigilance data and in a simulation study.

1. Introduction
Consider the problem of testing n hypotheses H1 , . . . , Hn simultaneously. A classical approach to dealing with the multiplicity problem is to control the familywise
error rate (FWER), i.e. the probability of one or more false rejections. However,
when the number n of hypotheses is large, the ability to reject false hypotheses
is small. Therefore, alternative Type 1 error rates have been proposed that relax
control of FWER in order to reject more false hypothesis. One such errror rate is
the expected proportion of false rejections amongst all rejections, which is known
as the false discovery rate (FDR). Benjamini and Hochberg (1995) showed that the
so-called linear step-up procedure controls the FDR under independence, Benjamini
and Yekutieli (2001) proved that it also controls the FDR for certain types of positive dependence of the p-values. The latter publication also introduced a rescaled
version of the linear step up procedure - referred to as the Benjamini-Yekutieli (BY)
procedure hereafter - which guarantees FDR control under any type of dependence
among the p-values. In this paper we focus on the Benjamini-Yekutieli procedure.
Generally, the only available information on the distribution of p-values is that
it is stochastically larger than the uniform distribution. However, for discrete tests,
there are important practical situations in which the (conditional) p-value distributions under the null hypotheses are in fact known. Examples of such tests include
Fisher’s exact test (see e.g. Westfall and Wolfinger (1997) for an application to
clinical studies), the McNemar test (see e.g. Westfall et al. (2010) for an application to classification models) and the binomial test (see e.g. Chen and Doerge
(2015) for an application to next generation sequencing data). When the p-value
distributions are available it may be possible to improve classical multiple testing
procedures by incorporating discreteness of the distribution functions into the multiplicity adjustment while still controlling the Type 1 error rate. For the FWER,
such methods were obtained e.g. by Westfall and Wolfinger (1997) and Hommel and
Krummenauer (1998). For a review of such methods, see Gutman and Hochberg
(2007). For the FDR, Heller and Gur (2012), in the following abbreviated as [HG],
give an overview of existing approaches for taking discreteness into account. Moreover, they modify the Benjamini-Liu (BL) step-down procedure (Benjamini and
Liu, 1999) and prove FDR control for the case when the discrete test statistics are
independent as well as for a certain type of dependency. Recently, Chen and Doerge
(2015) have used grouping and weighting approaches for independent test statistics.
They point out that ’How to derive better FDR procedures in the discrete paradigm
remains an urgent but still unresolved problem.’ This paper is a contribution in
this direction, for the case where the p-values are arbitrarily dependent.
Date: March 12, 2018.
1

2

SEBASTIAN DÖHLER DARMSTADT UNIVERSITY OF APPLIED SCIENCES

In this paper we develop step-up procedures that control the FDR under arbitrary dependence of the p-values and improve upon known procedures when the
test statistics are discrete. Our approach is similar in spirit to the one described by
Westfall and Wolfinger (1997). The paper is organised as follows. First we present
some theoretical results, building on the work of Sarkar (2008a). Then we compare
the performance with some of the procedures considered by [HG] by analysing a
pharmacovigilance data set provided by [HG] in the R package discreteMTP and
in a simulation study. The paper concludes with a discussion.
2. Theoretical results
When testing hypotheses H1 , . . . , Hn we assume that corresponding p-values
P V1 , . . . , P Vn ∈ [0, 1] with distribution functions F1 , . . . , Fn under the null hypotheses are available. For x ∈ [0, 1] define the functions G(x) = F1 (x) + · · · + Fn (x) and
Gunif (x) = n · x. Following Lehmann and Romano (2005), the only distributional
assumption we make in this paper is the following.
Assumption 1. For any true hypothesis Hi we assume that the distribution of the
p-value P Vi is stochastically larger than a uniform random variable, i.e. Fi (u) ≤ u
for all u ∈ (0, 1).
Let C = {c ∈ [0, 1]n |c1 ≤ · · · ≤ cn } denote the set of non-decreasing critical values. For c ∈ C the corresponding step-up procedure rejects hypotheses
H(1) , . . . , H(k) , where k = max{i|P V(i) ≤ ci }. If no such i exists, no hypothesis is
rejected. For the corresponding step-down procedure, reject H(1) , . . . , H(k) , where
k = max{i|P V(j) ≤ cj , j = 1, . . . , i}. If P V(1) > c1 , no hypothesis is rejected.
Our main tool for proving FDR control for discrete step-up procedures is an
upper bound of the FDR in terms of the function G.
Proposition 1. Let P V1 , . . . , P Vn satisfy assumption 1. For any sequence c ∈ C
we have for the step-up procedure (with c0 = 0)
(1)

(2)

FDR(c) ≤

≤

n X
X
1
r=1 j∈I
n
X
r=1

r

(Fj (cr ) − Fj (cr−1 ))

1
(G(cr ) − G(cr−1 ))
r

where I ⊂ {1, . . . , n} denotes the set of true null hypotheses.
Comments. If P Vi ∼ U (0, 1) for i ∈ I, the bound (1) yields the statement of
Theorem 3.5 in Sarkar (2008a).
We have introduced assumption 1 mainly because it is commonly considered a
desirable property of p-values. However, as the following proof of the proposition
shows, we actually only need that Fj (0) = 0 and so this weaker assumption is
sufficient for the proposition to hold. For discrete tests this is relevant because there
are usually different ways of defining p-values, see Hirji (2006) for a comprehensive
review of such possibilities. For instance, [HG] consider in their analysis both exact
and mid-p-values, which may not satisfy assumption 1 but the weaker assumption
Fj (0) = 0. For mid-p-values this implies that the proposition remains true with
corresponding function G as long as P (mid − P Vi = 0) = 0. Thus the results
presented here carry over naturally and flexibly to mid-p-values (and other variants)
but since we are mainly interested in exact p-value methods we do not pursue this
further.
Proof. From Sarkar (2007), Lemma 3.3, we obtain (with k = 1 and C10 = I)

n X 
X
X
1{P Vj ≤cr }
(−j)
−
FDR ≤
P (P Vj ≤ c1 ) +
E P (Rn−1 ) ≥ r − 1|P Vj ) ·
r
r=2
j∈I

j∈I

1{P Vj ≤cr−1 }
r−1



DISCRETE FDR STEP-UP PROCEDURES

3

(−j)

where Rn−1 is the number of rejections for the step-up procedure based on {P Vi |i ∈
(−j)
I \ {j}} and the critical values c2 ≤ · · · ≤ cn . Since P (Rn−1 ) ≥ r − 1|P Vj ) ∈ [0, 1]
and

1{P Vj ≤cr }
r

−

1{P Vj ≤cr−1 }
r−1

≤

1
(1{P Vj ≤cr } − 1{P Vj ≤cr−1 } )
r

we obtain by taking expectations
FDR ≤

X

Fj (c1 ) +

n X
X
1
r=2 j∈I

j∈I

r

(Fj (cr ) − Fj (cr−1 )).

Since Fj (c0 ) = 0 by assumption 1, the result follows.



Under assumption 1, Benjamini and Yekutieli (2001) showed that the step-up
procedure (in the following abbreviated as BY) based on the critical constants
n

(3)

X1
i
ci =
α, where DBY =
n · DBY
i
i=1

controls the FDR at level α, while Sarkar (2008b) obtained a similar result for
i(i + 1)
α.
2n2
From the proposition we obtain a simple way of constructing step-up procedures
that control the FDR. The following result is partly a generalisation of Theorem
3.5 of Sarkar (2008a).
(4)

ci =

Corollary 1. Let P V1 , . . . , P Vn satisfy assumption 1. For any sequence 0 = ye0 ≤
ye1 ≤ · · · ≤ yen define
D=

n
X
1
i=1

i

(e
yi − yei−1 )

and for α ∈ (0, 1) let c ∈ C satisfy
(5)

α
yei
for i = 1, . . . , n.
D
(a) Then the step-up procedure with critical constants c controls the FDR at
level α under arbitrary dependence of the p-values.
(b) For the choice yei = i we obtain a modified Benjamini-Yekutieli procedure.
(c) For the choice yei = i(i + 1) we obtain a modification of Sarkar’s procedure.
G(ci ) ≤ yi :=

Proof. For statement (a) we have by (2)
FDR(c) ≤

n
X
1
i=1

≤

n−1
X
i=1

=

i

(G(ci ) − G(ci−1 )) =

n−1
X
i=1

1
1
yi + yn =
i(i + 1)
n

n
X
1
i=1

i

1
1
G(ci ) + G(cn )
i(i + 1)
n

(yi − yi−1 )

n
α X1
(e
yi − yei−1 ) = α
D i=1 i

where the second inequality follows from (5).
For statement (b) we have
D=

n
X
1
i=1

i

(i − (i − 1)) =

n
X
1
i=1

i

= DBY .

If G = Gunif we recover the Benjamini-Yekutieli given by (3).

4

SEBASTIAN DÖHLER DARMSTADT UNIVERSITY OF APPLIED SCIENCES

For statement (c) we have
n
X
1
D=
(i(i + 1) − (i − 1)i) = 2n
i
i=1
which is the rescaling constant from Sarkar’s procedure. If G = Gunif we recover
Sarkar’s constants from (4).

Comments. The corollary presents a generic method of constructing FDR controlling step-up procedures under general dependency. These procedures basically work
by choosing values y1 , . . . , yn and comparing these to the values of the function
G for the observed p-values and for arbitrary distribution of the p-values under
the null hypotheses. More specifically, the step-up procedure rejects hypotheses
H(1) , . . . , H(k) , where k = max{i|G(P V(i) ) ≤ yi }. Adjusted p-values for the procedure defined by (5) are obtained using e.g. Procedure 1.4 from Dudoit and van der
Laan (2007):



D
(6)
G(pv(j) ), 1
.
pf
v (i) = min
min
j=i,...,n
j
When P V1 , . . . , P Vn ∼ U (0, 1) we have G(x) = Gunif (x) = n · x and we obtain ci = n1 yi which yields the Benjamini-Yekutieli constants for yei = i and the
procedure of Sarkar (2008a) for yei = i(i + 1). When the p-values are distributed
discretely, the function G(x) = F1 (x) + · · · + Fn (x) may be considerably smaller
than Gunif , as explained in [HG]. Thus, procedures based on Gunif rather than on G
may be substantially (and needlessly) more conservative. We present an example
below. As Westfall and Wolfinger (1997) point out in the context of FWER control,
the degree of improvement ’depends upon the specific characteristics of the discrete
distributions. Larger gains are possible when the number of tests is large and where
many variables are sparse’.
For independent statistics, Heyse (2011) introduces a modification of the BH
procedure (DBH) by defining critical constants ci by G(ci ) ≤ i·α. This corresponds
to our modified BY procedure where D is set to 1. Since DBY = log(m) + γ with
γ ≈ 0.57721 the Euler-Mascheroni constant, this procedure is considerably more
powerful than the discrete BY procedure. However, as we discuss in the appendix,
the DBH procedure does not generally control the FDR at the desired level.
Example. Figure 1 presents the graphs of G and Gunif for the pharmacovigilance
data from [HG] which is described in more detail in the next section. In this case,
n = 2446. For x ∈ (0, 0.001) the shape of G is roughly linear and a numerical com1
parison yields GG
≈ 3.8
. This means that the classical step-up critical constants
unif
based on Gunif like the BY or Sarkar procedure can be multiplied by a factor of
roughly 3.8 and still maintain FDR control in this special discrete setting.
3. Empirical data
We revisit the pharmacovigilance data analysed by [HG], to which we also refer
for more details. This data is derived from a database for reporting, investigating and monitoring adverse drug reactions due to the Medicines and Healthcare
products Regulatory Agency in the United Kingdom. It contains the number of
reported cases of amnesia as well as the total number of adverse events reported
for each of the 2446 drugs in the database. Overall, there are 2051 cases of amnesia. Heller and Gur (2012) investigate the association between reports of amnesia
and suspected drugs by performing for each drug Fisher’s exact test (one-sided) for
testing for association between the drug and amnesia and adjusting for multiplicity
by using several (discrete) FDR procedures.
The primary goal of our analysis is to compare the (continuous) BenjaminiYekutieli and Sarkar procedures to the discrete modifications introduced in section 2. All these methods guarantee FDR control under arbitrary dependence
of the p-values. As a secondary analysis we also present results for the (continuous) Benjamini-Hochberg and Benjamini-Liu procedures as well as the discrete

DISCRETE FDR STEP-UP PROCEDURES

5

y

0.0

0.1

0.2

0.3

0.4

0.5

0.6

G and Gunif

0e+00

2e−04

4e−04

6e−04

8e−04

1e−03

x

Figure 1. Graphs of G (solid line) and Gunif (dotted line) for the
pharmacovigilance data from [HG]
Benjamini-Liu procedure developed by [HG] and the discrete BH procedure due to
Heyse (2011).
Table 1 presents the results of our analysis. For the BY and Sarkar procedures
the adjusted p-values are reduced roughly by a factor of 3–4 by the discrete modifications as expected. Except for the Sarkar procedure all discrete modifications
lead to some additional rejections. While the BH, BL and DBL procedures control
FDR only for certain types of dependence, the results in the first four columns
of the table are valid under any type of dependence. The discrete BY not only
rejects the same hypotheses as the discrete BL procedure, but, as a comparison of
columns 2 and 8 shows, its adjusted p-values are only half as large. For instance,
while the adjusted p-value for association between Fluoxetine and amnesia is 12.4%
for the discrete BL, it barely fails to be significant at the 5%-level for the discrete
BY. Thus, at least in this setting, the discrete BY procedure, while requiring less
restrictive dependence assumptions, may be more powerful than the BL procedure.
As expected, the DBH procedure rejects the most hypothesis, followed by the BH
procedure. While the focus of this paper is on FDR procedures for arbitrarily
dependent data, we would like to follow one of the reviewers who noted that for
safety data, treating the data as independent may be preferable. Thus, it is not our
intention of suggesting that the (discrete) BY procedure would be the ideal method
for safety analyses. Rather, the rationale for analysing this specific data set was to
facilitate a comparison to the work of HG.
4. Simulation study
We now investigate the power of the procedures from the previous section in
a simulation study similar to those described in Gilbert (2005) and Heller and
Gur (2012). Although our primary focus is on comparing the discrete BenjaminiYekutieli procedure with its continuous counterpart, we also evaluate the performance of the BH procedure as a classical benchmark as well as the discrete BH and
Benjamini-Liu procedures (we omit the continuous Benjamini-Liu procedure).
4.1. Simulated Scenarios. We simulate a two-sample problem in which a vector
of m independent binary responses (“adverse events”) is observed for each subject
in the two groups. In order to investigate the influence of discreteness, we consider
two cases:
a) High degree of discreteness: Each group consists of N = 25 subjects.

6

SEBASTIAN DÖHLER DARMSTADT UNIVERSITY OF APPLIED SCIENCES

Table 1. The 27 smallest adjusted p-values from the multiple
testing procedures for the pharmacovigilance data. The columns
from left to right are the adjusted p-values from the (continuous
and discrete) Benjamini-Yekutieli procedure, the (continuous and
discrete) Sarkar procedure, the (continuous and discrete) BH procedure and the (continuous and discrete) BL procedure.

BUPROPION
CITALOPRAM
DEXAMPHETAMINE
ETHANOL
FLUOXETINE
GABAPENTIN
INDOMETHACIN
LACOSAMIDE
LEVETIRACETAM
LITHIUM
LORAZEPAM
MEFLOQUINE
MIDAZOLAM
OXCARBAZEPINE
PAROXETINE
PREGABALIN
RIMONABANT
SERTRALINE
SIMVASTATIN
STRONTIUM RANELATE
TEMAZEPAM
TOPIRAMATE
TRIAZOLAM
VARENICLINE
VIGABATRIN
ZOLPIDEM
ZOPICLONE

BY

DBY

Sarkar

DSarkar

BH

DBH

BL

DBL

0.4161
0.0080
0.0606
1.0000
0.1956
0.0000
0.0000
0.2752
0.0446
0.0012
0.0000
0.0000
0.0189
1.0000
0.0000
0.0000
0.0000
0.5762
0.0000
0.0550
0.0007
0.0380
0.0000
0.0000
0.0418
0.0000
0.0000

0.1112
0.0018
0.0169
0.2709
0.0518
0.0000
0.0000
0.0710
0.0117
0.0003
0.0000
0.0000
0.0050
0.2912
0.0000
0.0000
0.0000
0.1547
0.0000
0.0156
0.0002
0.0104
0.0000
0.0000
0.0113
0.0000
0.0000

1.0000
0.2902
1.0000
1.0000
1.0000
0.0000
0.0000
1.0000
1.0000
0.0464
0.0000
0.0000
0.6480
1.0000
0.0000
0.0000
0.0000
1.0000
0.0000
1.0000
0.0312
1.0000
0.0000
0.0000
1.0000
0.0001
0.0000

1.0000
0.0661
0.4474
1.0000
1.0000
0.0000
0.0000
1.0000
0.3408
0.0105
0.0000
0.0000
0.1723
1.0000
0.0000
0.0000
0.0000
1.0000
0.0000
0.4340
0.0072
0.3373
0.0000
0.0000
0.3408
0.0000
0.0000

0.0497
0.0009
0.0072
0.1263
0.0233
0.0000
0.0000
0.0328
0.0053
0.0001
0.0000
0.0000
0.0023
0.1364
0.0000
0.0000
0.0000
0.0688
0.0000
0.0066
0.0001
0.0045
0.0000
0.0000
0.0050
0.0000
0.0000

0.0133
0.0002
0.0020
0.0323
0.0062
0.0000
0.0000
0.0085
0.0014
0.0000
0.0000
0.0000
0.0006
0.0348
0.0000
0.0000
0.0000
0.0185
0.0000
0.0019
0.0000
0.0012
0.0000
0.0000
0.0013
0.0000
0.0000

0.6864
0.0140
0.1386
0.9515
0.3956
0.0000
0.0000
0.5222
0.0948
0.0020
0.0000
0.0000
0.0350
0.9636
0.0000
0.0000
0.0000
0.8098
0.0000
0.1211
0.0011
0.0732
0.0000
0.0000
0.0846
0.0000
0.0000

0.2659
0.0032
0.0403
0.5541
0.1239
0.0000
0.0000
0.1721
0.0255
0.0004
0.0000
0.0000
0.0094
0.5931
0.0000
0.0000
0.0000
0.3599
0.0000
0.0357
0.0003
0.0205
0.0000
0.0000
0.0234
0.0000
0.0000

19

21

14

14

24

27

16

21

Number rejections

b) Moderate degree of discreteness: Each group consists of N = 100 subjects.
Then, the goal is to simultaneously test the m null hypotheses H0i : p1i = p2i , i =
1, . . . , m, where p1i and p2i are the success probabilities for the ith binary response
in group 1 and 2 respectively. We take m = 800, 2000 where m = m1 + m2 + m3
and data are generated so that the response is Bernoulli(0.01) at m1 positions for
both groups, Bernoulli(0.10) at m2 positions for both groups and Bernoulli(0.10)
at m3 positions for group 1 and Bernoulli(q) at m3 positions for group 2 where
q = 0.15, 0.25, 0.4 represents weak, moderate and strong effects respectively. The
null hypothesis is true for the m1 and m2 positions while the alternative hypothesis
is true for the m3 positions. We also take different configurations for the proportion
of false null hypotheses, m3 is set to be 10%, 30% and 60% of the value of m,
which represents the complete null hypothesis and small, intermediate and large
proportion of effects (the proportion of true nulls π0 is 1, 0.9, 0.7, 0.4, respectively).
Then, m1 is set to be 20%, 50% and 80% of the true nulls (m − m3 ) and m2 =
m − m1 − m3 .
For each of the 54 possible parameter configurations specified by m, m3 , m1 and
q3 , 10000 Monte Carlo trials are performed, that is, 10000 data sets are generated and for each data set, an unadjusted two-sided p-value from Fisher’s exact
test is computed for each of the m positions, and the multiple testing procedures
mentioned above are applied at level α = 0.05. The power of each procedure was
estimated as the fraction of the m3 false null hypotheses that were rejected, averaged over the 10000 simulations. For random number generation the R-function
rbinom was used. The two-sided p-values from Fisher’s exact test were computed
using the R-function fisher.test.
4.2. Simulation results.

DISCRETE FDR STEP-UP PROCEDURES

7

4.2.1. High degree of discreteness. Table 3 in the Appendix displays the (average)
power of the seven procedures under investigation. For weak and moderate effects,
i.e. q3 = 0.15 and q3 = 0.25, none of the procedure possesses relevant power. For
strong effects, the discrete procedures improve considerably on their standard counterparts. Both the discrete and the standard Sarkar procedure are outperformed
by the other procedures. For fixed m3 and q3 , the power of the discrete BenjaminiYekutieli procedure is slightly increasing in m1 . This may seem surprising, since
the m1 positions represent true null hypotheses. The reason for this behavior is
that because of the low success probability, these positions are accountable for the
‘most discrete’ part of the data. The more discrete the data are, the flatter the
function G (see section 2 and figure 1) will tend to be, leading to higher critical
values for the discrete procedure.
For strong effects, the results are also displayed in figure 2. In order to avoid
over-optimism we used, as explained above, for fixed m3 and q3 the configuration
with smallest m1 . Figure 2 shows that the discrete Benjamini-Yekutieli procedure
N=25

0.4
0.2
0.0

average power

0.6

m= 800

small

large prop.

0.2

0.4

0.6

m= 2000

0.0

average power

intermediate

small

intermediate

large prop.

Figure 2. Average power (high degree of discreteness) of the BYprocedure (orange), DBY-procedure (red), DBL-procedure (blue),
BH-procedure (green) and DBH (grey) for small, intermediate and
large proportions of alternatives with strong effects. The top panel
presents results for m = 800, the lower panel for m = 2000.
performs much better than its continuous counterpart. As expected, the discrete
Benjamini-Hochberg procedure performs best amongst all procedures. For small
and intermediate proportions of large effects, the discrete Benjamini-Yekutieli outperforms the continuous Benjamini-Hochberg procedure, which may be surprising
given that the standard Benjamini-Yekutieli procedure has a reputation for being
extremely conservative. This shows that by incorporating discreteness, even the
Benjamini-Yekutieli procedure may become a powerful tool in some settings. This
observation is similar to the findings of Westfall and Wolfinger (1997) with respect
to the Bonferroni procedure.
4.2.2. Moderate degree of discreteness. Table 4 in the Appendix displays the (average) power of the seven procedures under investigation, for N = 100. For weak
effects, again none of the procedure possesses any relevant power. For intermediate
and strong effects the situation is different from the highly discrete case. For strong
effects, all procedures exhibit roughly the same, extremely high, power. For moderate effects, figure 3 shows that while the discrete Benjamini-Yekutieli procedure
still improves on its continuous counterpart, the gains are now much smaller than
in the sparse case. A similar conclusion holds true for the DBH procedure.

8

SEBASTIAN DÖHLER DARMSTADT UNIVERSITY OF APPLIED SCIENCES

N=100

0.4
0.2
0.0

average power

0.6

m= 800

small

intermediate

large prop.

0.4
0.2
0.0

average power

0.6

m= 2000

small

intermediate

large prop.

Figure 3. Average power (moderate degree of discreteness) of
the BY-procedure (orange), DBY-procedure (red), DBL-procedure
(blue), BH-procedure (green) and DBH (grey) for small, intermediate and large proportions of alternatives with intermediate effects.
The top panel presents results for m = 800, the lower panel for
m = 2000.
5. Discussion
In this paper we have introduced a discrete modification of the BenjaminiYekutieli procedure and investigated its performance for an empirical data set and
in a simulation study. We have shown that this procedure can perform remarkably
well when the number of tests is large and many variables are sparse. Thus we
recommend this procedure for discrete data with arbitrary dependency structure.
The proof of proposition 1 shows that due to (1) the FDR is actually controlled
by
FDR(c) ≤

n
X
1
r=1

r

(GI (cr ) − GI (cr−1 )),

P
where GI (x) = j∈I Fj (x) and I ⊂ {1, . . . , n} denotes the set of true hypotheses. Since n0 = |I| is unknown it might be possible to improve our procedures by
incorporating an appropriately chosen estimate of n0 . Such adaptive procedures
have been proposed for the BH procedure by Benjamini and Hochberg (2000), Benjamini et al. (2006) and Storey et al. (2004) and others under independence. Sarkar
(2008b) investigates such two-stage procedures also for the general dependence case
and obtains modifications of the BY procedure. It would be intersting to see if the
procedures introduced in this paper can be modified in a similar way and if this
leads to more powerful procedures.
The FDR is the expected value of the false discovery proportion FDP = V /R,
where R denotes the number of rejected hypotheses and V the number of falsely
rejected hypotheses. Guo et al. (2012) generalize the notion of FDP to the situation
where a small number k of false rejections is deemed acceptable. More specifically,
they introduce kFDP = V /R if V ≥ k and 0 otherwise, and the kFDR is defined
as E(kFDP). It can be shown that similar representations to our proposition 1 are
available for the kFDR. However, the upper bound in the proposition now involves
determining nk distribution functions and so this approach seems computationally
feasible only for special situations.
Finally, we mention that for independent p-values it still remains a challenge to
develop a discrete modification of the Benjamini-Hochberg procedure with guaranteed FDR control.

DISCRETE FDR STEP-UP PROCEDURES

9

Acknowledgments
The author would like to thank Irene Castro-Conde and Ruth Heller for helpful
discussions and Florian Junge for greatly improving the simulation code.
References
Benjamini, Y. and Y. Hochberg (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical
Society. Series B 57 (1), 289–300.
Benjamini, Y. and Y. Hochberg (2000). On the adaptive control of the false discovery rate in multiple testing with independent statistics. Journal of Educational
and Behavioral Statistics 25 (1), 60–83.
Benjamini, Y., A. M. Krieger, and D. Yekutieli (2006). Adaptive linear step-up
procedures that control the false discovery rate. Biometrika 93 (3), 491–507.
Benjamini, Y. and W. Liu (1999). A step-down multiple hypotheses testing procedure that controls the false discovery rate under independence. Journal of
Statistical Planning and Inference 82 (1-2), 163 – 170.
Benjamini, Y. and D. Yekutieli (2001). The control of the false discovery rate in
multiple testing under dependency. The Annals of Statistics 29 (4), 1165–1188.
Chen, X. and R. Doerge (2015). A weighted fdr procedure under discrete and
heterogeneous null distributions. arXiv:1502.00973 .
Dudoit, S. and M. J. van der Laan (2007). Multiple Testing Procedures and Applications to Genomics. Springer Series in Statistics. Springer. ISBN: 978-0-38749316-9.
Gilbert, P. (2005). A modified false discovery rate multiple-comparisons procedure
for discrete data, applied to human immunodeficiency virus genetics. Journal of
the Royal Statistical Society. Series C 54 (1), 143–158.
Guo, W., L. He, and S. Sarkar (2012). Further results on controlling the false
discovery proportion. Annals of Statistics. (under revision).
Gutman, R. and Y. Hochberg (2007). Improved multiple test procedures for discrete
distributions: New ideas and analytical review. Journal of Statistical Planning
and Inference 137 (7), 2380 – 2393.
Heller, R. and H. Gur (2012). False discovery rate controlling procedures for discrete
tests. arxiv:1112.4627v2 .
Heyse, J. F. (2011). A false discovery rate procedure for categorical data. pp. 43–
58. Recent Advances in Bio- statistics: False Discovery Rates, Survival Analysis,
and Related Topics.
Hirji, K. (2006). Exact Analysis of Discrete Data. Chapman & Hall/CRC.
Hommel, G. and F. Krummenauer (1998). Improvements and modifications of
tarone’s multiple test procedure for discrete data. Biometrics 54 (2), pp. 673–
681.
Lehmann, E. and J. P. Romano (2005). Generalizations of the familywise error
rate. Ann. Stat. 33 (3), 1138–1154.
Sarkar, S. K. (2007). Stepup procedures controlling generalized FWER and generalized FDR. Ann. Stat. 35 (6), 2405–2420.
Sarkar, S. K. (2008a). On Methods Controlling the False Discovery Rate.
Sankhya 70-A(2), 135–168.
Sarkar, S. K. (2008b). Two-stage stepup procedures controlling fdr. Journal of
Statistical Planning and Inference 138 (4), 1072 – 1084.
Storey, J. D., J. E. Taylor, and D. Siegmund (2004). Strong control, conservative
point estimation and simultaneous conservative consistency of false discovery
rates: a unified approach. Journal of the Royal Statistical Society: Series B
(Statistical Methodology) 66 (1), 187–205.
Westfall, P. and R. Wolfinger (1997). Multiple tests with discrete distributions.
The American Statistician 51 (1), 3–8.
Westfall, P. H., J. F. Troendle, and G. Pennello (2010). Multiple mcnemar tests.
Biometrics 66 (4), 1185–1191.

10

SEBASTIAN DÖHLER DARMSTADT UNIVERSITY OF APPLIED SCIENCES

Sebastian Döhler, University of Applied Sciences, CCSOR and Faculty of Mathematics and Science, D-64295 Darmstadt, Germany
E-mail address: sebastian.doehler@h-da.de

DISCRETE FDR STEP-UP PROCEDURES

11

Appendix 1
In this appendix we present a numerical example for which the procedure of
Heyse (2011) is anticonservative. [HG] already constructed a similar example, using
the definition (6) of adjusted p-values. However, as one referee pointed out, this
definition differs from the one originally given in Heyse (2011), which we recap here
(in our notation):
(7)

pf
v (n) = pv(n) ,

and

pf
v (i) = min{f
pv (i+1) , G(pv(i) )/i},

for i = 1, . . . , n − 1. Since G(pv(n) )/n ≤ pv(n) , definition (7) may yield larger
adjusted p-values than definition (6). In fact, as the referee pointed out, when the
example of [HG] is analysed using (7) instead of (6), the FDR level is no longer
larger than the nominal level.
Numerical example. In what follows we present a numerical example where even
the more conservative procedure given by definition (7) is anticonservative. Consider n = 4 p-value distributions (under the null hypotheses) defined by
P1 = 0.05 · δ{0.05} + 0.95 · δ{1} ,
P2 = 0.025 · δ{0.10} + 0.975 · δ{1} ,
P3 = 0.025 · δ{0.15} + 0.975 · δ{1}

and

P4 = δ{1} .
where δ{x} denotes the Dirac measure at x. Clearly, these distributions are stochastically larger than the uniform distribution and the function G = F1 + F2 + F3 + F4
can easily be determined. Now consider the configuration where all hypotheses
are true and use the discrete FDR procedure to control the FDR = FWER at
level α = 0.05 (since P V(4) = 1 definitions (6) and (7) coincide). This means
that all hypotheses with pf
v (j) ≤ α are declared significant. Table 2 illustrates
the situation in detail. All possible combinations of p-values are listed in the first
three columns (pv4 = 1 is omitted). The sorted values are shown in columns 3
to 6, the adjusted p-values are given in columns 7 to 9. Columns 10 to 12 list the
probability masses for the p-values in the first 3 columns. Since we assume independence, the joint probability of observing (pv1 , pv2 , pv3 ) is given by the product of
the probabilities (last column) and the joint probabilities for those constellations of
p-values for which at least one hypothtesis is rejected are printed in boldface. Thus
the probability of at least one rejection can be computed exactly and we obtain
FWER = FDR = 0.05059375.
In principle this type of example can be extended easily to larger n, however
since 2n−1 events have to be enumerated we soon run into computational difficulties. As a slight extension of the example above consider n = 10. Assume P1 as
above, P10 = δ{1} and P2 , . . . , P9 each with probability mass 0.00621 at locations
0.10, . . . , 0.45. In this case (the details are not presented) similar calculations as
above yield FWER = FDR = 0.05100062.
Mathematical considerations. The above example implies that the reasoning
provided in Heyse (2011) is inadequate for proving FDR control of this procedure.
We now try to identify which part of this reasoning may be problematic. Heyse
(2011) argued that this procedure controls the FDR by the same argument used
in Gilbert (2005). The latter paper appealed to equation (19) of Benjamini and
Yekutieli (2001), which we reproduce here (with our notation): For a given sequence
c1 ≤ · · · ≤ cn of critical constants the FDR can be expressed by
n
XX
1
(8)
Fj (cr ) · P (Cr(j) )
FDR =
r
r=1
j∈I

(j)

where the event Cr is defined on page 1178 of Benjamini and Yekutieli (2001). The
(j)
(j)
crucial point is that the events Cr and probabilities P (Cr ) are determined (in a

pv2

0.1000
0.1000
1.0000
1.0000
0.1000
0.1000
1.0000
1.0000

pv1

0.0500
1.0000
0.0500
1.0000
0.0500
1.0000
0.0500
1.0000

0.1500
0.1500
0.1500
0.1500
1.0000
1.0000
1.0000
1.0000

pv3
0.0500
0.1000
0.0500
0.1500
0.0500
0.1000
0.0500
1.0000

pv(1)
0.1000
0.1500
0.1500
1.0000
0.1000
1.0000
1.0000
1.0000

pv(2)
0.1500
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000

pv(3)
0.0333
0.0500
0.0500
0.1000
0.0375
0.0750
0.0500
1.0000

pf
v (1)
0.0333
0.0500
0.0500
1.0000
0.0375
1.0000
1.0000
1.0000

pf
v (2)
0.0333
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000

pf
v (3)
0.0500
0.9500
0.0500
0.9500
0.0500
0.9500
0.0500
0.9500

0.0250
0.0250
0.9750
0.9750
0.0250
0.0250
0.9750
0.9750

0.0250
0.0250
0.0250
0.0250
0.9750
0.9750
0.9750
0.9750
FWER

P (P V1 = pv1 ) P (P V2 = pv2 ) P (P V3 = pv3 )

Table 2. Details of the example

0.00003125
0.00059375
0.00121875
0.02315625
0.00121875
0.02315625
0.04753125
0.90309375
0.05059375

joint prob.

12
SEBASTIAN DÖHLER DARMSTADT UNIVERSITY OF APPLIED SCIENCES

DISCRETE FDR STEP-UP PROCEDURES

13

non-trivial way) by the given multiple testing procedure and the distribution of the
(j)
p-values. In Gilbert (2005), however, the values P (Cr ) are chosen (arbitrarily)
without further mathematical justification. The author seems to have been aware
that this reasoning does not constitute a valid mathematical proof of FDR control
since he states that this procedure ’usually’ (p. 151) controls the FDR. Thus, while
this procedure has some heuristic justification, its exact mathematical properties
remain unclear.
Discussion. We have shown that the DBH procedure of Heyse (2011) does not
generally guarantee FDR control at the desired level. We have also tried to identify
the mathematical reason for this. It might be argued that the magnitude of the
bias in the example(s) is rather small. However, it should be kept in mind that
these examples show the behaviour for specific constellations and it is not by any
means clear that these are ’worst case scenarios’ (which is what type 1 error control
is all about). With more mathematical and computational effort it may also be
possible to find more severe examples of FDR inflation (e.g. it is unclear what
may happen when m is large, some hypotheses are false, etc). Thus our example is
not meant to quantify the magnitude of bias, but rather to demonstrate that the
DBH procedure does not rigorously control FDR. To the best of our knowledge, no
(valid) upper FDR bound is available for the DBH procedure. Thus it still remains
an open problem to find a discrete modification of the BH procedure that rigorously
controls the FDR.

14

SEBASTIAN DÖHLER DARMSTADT UNIVERSITY OF APPLIED SCIENCES

Appendix 2

DISCRETE FDR STEP-UP PROCEDURES

15

Table 3. Average power of FDR procedures for N = 25
m

m3

m1

q3

DBH

BH

DBL

DBY

BY

DSarkar

Sarkar

800

80

144
144
144
360
360
360
576
576
576
112
112
112
280
280
280
448
448
448
64
64
64
160
160
160
256
256
256
360
360
360
900
900
900
1440
1440
1440
280
280
280
700
700
700
1120
1120
1120
160
160
160
400
400
400
640
640
640

0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4

0.0002
0.0157
0.4195
0.0005
0.0188
0.4511
0.0006
0.0233
0.4789
0.0002
0.0219
0.5095
0.0002
0.0242
0.5528
0.0003
0.0279
0.5875
0.0002
0.0279
0.5895
0.0002
0.0290
0.5903
0.0002
0.0306
0.6064
0.0002
0.0127
0.4293
0.0002
0.0149
0.4520
0.0002
0.0188
0.4700
0.0002
0.0186
0.5135
0.0002
0.0214
0.5539
0.0002
0.0265
0.5897
0.0001
0.0270
0.5907
0.0001
0.0288
0.5908
0.0001
0.0309
0.6059

0.0000
0.0003
0.0795
0.0000
0.0003
0.0795
0.0000
0.0003
0.0794
0.0000
0.0005
0.2142
0.0000
0.0005
0.2141
0.0000
0.0005
0.2140
0.0000
0.0007
0.3727
0.0000
0.0007
0.3725
0.0000
0.0007
0.3722
0.0000
0.0001
0.0726
0.0000
0.0001
0.0725
0.0000
0.0001
0.0725
0.0000
0.0001
0.2056
0.0000
0.0001
0.2055
0.0000
0.0001
0.2055
0.0000
0.0002
0.3749
0.0000
0.0002
0.3745
0.0000
0.0002
0.3742

0.0004
0.0084
0.1457
0.0007
0.0096
0.1466
0.0007
0.0152
0.1551
0.0002
0.0081
0.0942
0.0002
0.0082
0.0955
0.0005
0.0082
0.0961
0.0002
0.0031
0.0732
0.0002
0.0031
0.0732
0.0002
0.0032
0.0733
0.0002
0.0075
0.0887
0.0002
0.0081
0.0949
0.0003
0.0082
0.0975
0.0002
0.0030
0.0555
0.0002
0.0030
0.0556
0.0002
0.0030
0.0558
0.0001
0.0030
0.0451
0.0001
0.0030
0.0451
0.0002
0.0030
0.0451

0.0002
0.0080
0.2131
0.0003
0.0087
0.2232
0.0006
0.0097
0.2416
0.0002
0.0062
0.2350
0.0002
0.0064
0.2494
0.0002
0.0067
0.2650
0.0002
0.0053
0.2664
0.0002
0.0054
0.2687
0.0002
0.0055
0.2697
0.0001
0.0048
0.1981
0.0002
0.0052
0.2197
0.0002
0.0056
0.2291
0.0001
0.0042
0.2287
0.0001
0.0043
0.2302
0.0001
0.0045
0.2497
0.0001
0.0040
0.2449
0.0001
0.0041
0.2573
0.0001
0.0041
0.2662

0.0000
0.0000
0.0086
0.0000
0.0000
0.0086
0.0000
0.0000
0.0086
0.0000
0.0000
0.0241
0.0000
0.0000
0.0241
0.0000
0.0000
0.0241
0.0000
0.0000
0.0586
0.0000
0.0000
0.0586
0.0000
0.0000
0.0586
0.0000
0.0000
0.0046
0.0000
0.0000
0.0046
0.0000
0.0000
0.0046
0.0000
0.0000
0.0155
0.0000
0.0000
0.0155
0.0000
0.0000
0.0155
0.0000
0.0000
0.0453
0.0000
0.0000
0.0453
0.0000
0.0000
0.0453

0.0000
0.0001
0.0100
0.0000
0.0001
0.0101
0.0000
0.0001
0.0103
0.0000
0.0000
0.0212
0.0000
0.0000
0.0217
0.0000
0.0000
0.0223
0.0000
0.0000
0.1181
0.0000
0.0000
0.1221
0.0000
0.0000
0.1267
0.0000
0.0000
0.0036
0.0000
0.0000
0.0036
0.0000
0.0000
0.0036
0.0000
0.0000
0.0072
0.0000
0.0000
0.0073
0.0000
0.0000
0.0075
0.0000
0.0000
0.1026
0.0000
0.0000
0.1067
0.0000
0.0000
0.1118

0.0000
0.0000
0.0001
0.0000
0.0000
0.0001
0.0000
0.0000
0.0001
0.0000
0.0000
0.0001
0.0000
0.0000
0.0001
0.0000
0.0000
0.0001
0.0000
0.0000
0.0001
0.0000
0.0000
0.0001
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000

240

480

2000

200

600

1200

16

SEBASTIAN DÖHLER DARMSTADT UNIVERSITY OF APPLIED SCIENCES

Table 4. Average power of FDR procedures for N = 100
m

m3

m1

q3

DBH

BH

DBL

DBY

BY

DSarkar

Sarkar

800

80

144
144
144
360
360
360
576
576
576
112
112
112
280
280
280
448
448
448
64
64
64
160
160
160
256
256
256
360
360
360
900
900
900
1440
1440
1440
280
280
280
700
700
700
1120
1120
1120
160
160
160
400
400
400
640
640
640

0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4
0.15
0.25
0.4

0.0022
0.4429
0.9931
0.0027
0.5114
0.9950
0.0049
0.6095
0.9973
0.0026
0.6213
0.9975
0.0035
0.6593
0.9980
0.0057
0.7133
0.9988
0.0038
0.7163
0.9988
0.0044
0.7312
0.9989
0.0057
0.7539
0.9991
0.0011
0.4390
0.9930
0.0015
0.5116
0.9949
0.0030
0.6091
0.9973
0.0015
0.6212
0.9976
0.0020
0.6564
0.9981
0.0034
0.7125
0.9988
0.0024
0.7151
0.9988
0.0029
0.7294
0.9989
0.0037
0.7531
0.9990

0.0010
0.3013
0.9848
0.0010
0.3005
0.9847
0.0010
0.2996
0.9847
0.0011
0.4915
0.9945
0.0011
0.4907
0.9944
0.0011
0.4897
0.9944
0.0013
0.6227
0.9976
0.0013
0.6221
0.9976
0.0013
0.6215
0.9976
0.0003
0.2994
0.9848
0.0003
0.2986
0.9848
0.0003
0.2977
0.9847
0.0004
0.4928
0.9946
0.0004
0.4920
0.9945
0.0004
0.4910
0.9945
0.0006
0.6244
0.9976
0.0006
0.6238
0.9976
0.0006
0.6232
0.9976

0.0019
0.1396
0.8955
0.0028
0.1656
0.9204
0.0043
0.2086
0.9451
0.0018
0.1391
0.9203
0.0022
0.1458
0.9384
0.0029
0.1773
0.9621
0.0016
0.1308
0.9488
0.0018
0.1387
0.9643
0.0019
0.1448
0.9794
0.0009
0.0948
0.8520
0.0013
0.1088
0.8793
0.0020
0.1433
0.9208
0.0009
0.0895
0.8775
0.0010
0.0988
0.8954
0.0015
0.1189
0.9270
0.0008
0.0882
0.9204
0.0009
0.0914
0.9382
0.0009
0.0976
0.9615

0.0005
0.1806
0.9646
0.0006
0.2205
0.9748
0.0010
0.3085
0.9842
0.0004
0.3030
0.9840
0.0005
0.3429
0.9858
0.0007
0.3933
0.9901
0.0004
0.3832
0.9898
0.0005
0.4114
0.9909
0.0006
0.4285
0.9916
0.0002
0.1688
0.9646
0.0003
0.2066
0.9749
0.0005
0.2926
0.9830
0.0002
0.2902
0.9828
0.0002
0.3198
0.9855
0.0003
0.3688
0.9886
0.0002
0.3654
0.9885
0.0002
0.3860
0.9898
0.0003
0.4166
0.9910

0.0002
0.0829
0.9396
0.0002
0.0828
0.9396
0.0002
0.0828
0.9396
0.0002
0.1751
0.9690
0.0002
0.1750
0.9690
0.0002
0.1749
0.9690
0.0002
0.2672
0.9825
0.0002
0.2672
0.9825
0.0002
0.2671
0.9825
0.0000
0.0707
0.9336
0.0000
0.0706
0.9336
0.0000
0.0706
0.9336
0.0001
0.1604
0.9680
0.0001
0.1603
0.9680
0.0001
0.1603
0.9680
0.0001
0.2494
0.9803
0.0001
0.2494
0.9803
0.0001
0.2493
0.9803

0.0000
0.0144
0.9413
0.0000
0.0191
0.9460
0.0000
0.0344
0.9652
0.0000
0.0776
0.9847
0.0000
0.1131
0.9872
0.0000
0.1776
0.9908
0.0000
0.3458
0.9951
0.0000
0.3790
0.9957
0.0000
0.4179
0.9963
0.0000
0.0063
0.9409
0.0000
0.0083
0.9454
0.0000
0.0171
0.9648
0.0000
0.0601
0.9848
0.0000
0.0970
0.9872
0.0000
0.1672
0.9909
0.0000
0.3449
0.9951
0.0000
0.3776
0.9957
0.0000
0.4197
0.9963

0.0000
0.0033
0.8897
0.0000
0.0033
0.8896
0.0000
0.0033
0.8896
0.0000
0.0094
0.9712
0.0000
0.0094
0.9712
0.0000
0.0094
0.9712
0.0000
0.1111
0.9910
0.0000
0.1110
0.9910
0.0000
0.1110
0.9910
0.0000
0.0009
0.8886
0.0000
0.0009
0.8886
0.0000
0.0009
0.8886
0.0000
0.0028
0.9716
0.0000
0.0028
0.9716
0.0000
0.0028
0.9716
0.0000
0.0958
0.9910
0.0000
0.0957
0.9910
0.0000
0.0956
0.9910

240

480

2000

200

600

1200

