University of Sydney, School of Mathematics and Statistics, Research Report 99-9, March,1999.

A STEPWISE REJECTIVE TEST PROCEDURE WITH STRONG CONTROL OF FAMILYWISE ERROR RATE By T. CHEN^{1} and E. SENETA

Department of Statistics, University of Pittsburgh, Pittsburgh, PA 15260, U.S.A. and School of Mathematics and Statistics, University of Sydney, NSW 2006, Australia

SUMMARY

In Holm’s stepwise rejective procedure the critical values are replaced adaptively by larger ones using a degree 2 inequality in place of Boole’s. The refinement retains strong control of familywise error rate. There is a cost in calculational simplicity; but a substantial improvement in actual error rate, according to simulations. Some key words: Holm’s procedure; Degree 2 inequalities; Dependent test statistics; Error rate; Power; Multivariate t-distribution; Strong control; Simes’ test.

1. INTRODUCTION

We consider the multiple test problem where there are n hypotheses H_{1}, H_{2}, ⋯, H_{n}, and corresponding p-values R_{1}, ⋯, R_{n} , assuming the test statistics X_{1}, ⋯, X_{n} are from a continuous distribution. Suppose that in a multiple test procedure the property

P (H_{s}, sϵI , are accepted|H_{s}, sϵI, true) ≥ 1 − α  (1)

holds, for prespecified size of test (familywise error rate) α , where I is any non-null subset of {1, 2, ⋯ n} , and thus contains m items, 1 ≤ m ≤ n . Then the procedure is said to control strongly the familywise error rate (Hochberg and Tamhane, 1987, pp. 3, 7).

Let R_{(1)}, R_{(2)}, ⋯, R_{(n)} be the ordered p-values, and H_{(1)}, H_{(2)}, ⋯, H_{(n)} the corresponding hypotheses. The ‘Bonferroni’ multiple test procedure rejects the composite hypothesis {H_{(1)}, H_{(2)}, ⋯, H_{(n)}} if R_{(1)} ≤ α/n , and accepts it otherwise. This procedure was refined by Holm (1979) as follows. Examine whether R_{(1)} ≤ α/n ; if not, accept H_{(i)} , i = 1, ⋯ n as with Bonferroni; if so, reject H_{(1)} and examine whether R_{(2)} ≤ α/(n − 1). If the inequality is not satisfied accept H_{(2)}, ⋯, H_{(n)} ; otherwise reject H_{(2)} . Continue in this way. To summarize: if R_{(i)} ≤ α/(n − i + 1), i ≤ j − 1 , then at step j the remaining hypotheses are H_{(j)}, ⋯, H_{(n)} and the inequality next to check is R_{(j)} ≤ α/(n − j + 1) . The process may run at most until a decision is made on the basis of whether R_{(n)} ≤ α or not. Holm showed that his procedure strongly controls the familywise error rate. Inasmuch as it essentially depends on Boole’s (first Bonferroni) inequality, which is a degree 1 bound (e.g. Seneta, 1988), Holm’s procedure retains an elegant simplicity.

1Work of this author done in part at the University of Sydney

There have been a number of improvements on the Bonferroni-Holm degree 1 procedures, all of which are aimed at increasing power while retaining a simple structure of critical points (such as α/(n − j + 1) above). Hochberg’s (1988) and Hommel’s (1988) (1989) procedures both strongly control the familywise error rate providing each Simes (1986) test of H_{I} = {H_{i}, iϵI} is a level α-test. Recently Sarkar (1998) has shown, that the Simes property holds if the random vector (X_{1}, ⋯, X_{n}) is from an MT P_{2} distribution. This class includes the central multivariate t with the associated correlation matrix having a common and non-negative correlation. Since they are based on the Simes procedure, the step-up procedures of Hochberg/Hommel, now validated in the MT P_{2} situation, are more powerful than Holm’s (1979) procedure.

In Section 6 of Seneta and Chen (1997), a degree 2 step-down procedure is proposed which retains familywise control of error rate. This procedure is adaptive in that calculation at each step is determined by the joint outcome of all pairs of statistics in the experiment involved until the procedure stops. In view of the continuing interest in a general procedure with familywise control of error rate, we present here a substantial refinement of this procedure, in a form which resembles Holm’s. Specifically, the values α/(n−j+1), j ≥ 1 , are replaced by larger ones, thus increasing the power. We present, using simulation, a crude power comparison with the Bonferroni/Holm procedure and with the Hochberg procedure in the setting of multivariate t . A more extensive comparison with Hochberg/Hommel, is indicated, but outside our present intention of presenting briefly the streamlined procedure. By comparison with Bonferroni/Holm and Hochberg/Hommel, there is a loss of computational simplicity. The first step in particular requires the calculation of n(n − 1)/2 two dimensional probabilities (in specific situations such as Dunnett’s test, fewer two dimensional probabilities need to be computed) ; and our procedure may be useful only if n is small. It may, however, be possible to obtain approximations to two dimensional probabilities by bootstrapping in specific situations (these are being investigated), thus eliminating a need to know two dimensional distributions. PROCEDURE CS

Write for the moment R_{(i)} = R_{t_{i}} , so t_{i} is a random variable from the set {1, 2, ⋯, n}. Using the ordered p-values R_{t_{i}}, i = 1, ⋯, n observed, define the index sets K(·) by K(p) = {t_{p}, t_{p+1}, ⋯, t_{n}}, p = 1, 2, ⋯, n (these are random sets for p ≥ 2) and write

γ(p) = max(jϵK(p)) ∑_{iϵK(p)−\{j\}} P (R_{i} ≤ \frac{α}{n − p + 1}, R_{j} ≤ \frac{α}{n − p + 1}|H_{s}, sϵK(p), true)

(2)

for 1 ≤ p ≤ n − 1, with γ(n) = 0 . These may be calculated for successive p as far as required in what follows.

Step 1 :

R_{(1)} ≤ min(\frac{α}{n − 1}, \frac{α + γ(1)}{n}) ?

If yes, reject H_{(1)} and go to Step 2. If no, accept H_{(1)}, H_{(2)}, ⋯, H_{(n)} and stop. Continue in this way. If the i-th step is reached:

Step i :

R_{(i)} ≤ min(\frac{α}{n − i}, \frac{α + γ(i)}{n − i + 1}) ?

If yes, reject H_{(i)} and go to Step i + 1 . If no, accept H_{(i)}, H_{(i+1)}, ⋯, H_{(n)} and stop. If the n-th Step is reached:

Step n :

R_{(n)} ≤ α ?

If yes, reject H_{(n)} and stop . If no, accept H_{(n)} and stop. (Note that since γ(n) = 0, Step n is consistent with the others. Also, at the (n−1)^{th} step, min(α, \frac{α + γ(n − 1)}{2}) =

\frac{α + γ(n − 1)}{2} , since γ(n − 1) ≤ \frac{α}{2}). □

The procedure is therefore an adaptive one, and the γ(p)′s for 2 ≤ p ≤ n − 1 are random variables.

2. STRONG CONTROL OF FAMILYWISE ERROR RATE A key feature of the proof of the theorem is the use of the inequality (from which

(2) derives) of Kounias (1968)

P ( ⋃_{i=1}^{k} A_{i}) ≤ ∑_{i=1}^{k} P (A_{i}) − max(j = 1, ⋯, k) ∑_{i≠j}^{k} P (A_{i} ⋂ A_{j}),

a second-degree inequality.

LEMMA. Let I, of fixed size m , be as in Section 1. Define

γ = max(I) ∑_{jϵI−\{i\}} P (R_{i} ≤ \frac{α}{m}, R_{j} ≤ \frac{α}{m}|H_{s}, sϵI , true),

where γ = 0 when m = 1. Then

⋂_{iϵI} \{R_{i} > min (\frac{(α + γ)}{m}, \frac{α}{m − 1})\} ⊆ \{H_{s}, sϵI , are all accepted\}.

\underline{Proof.} For any fixed realization (sample point) of the experiment within

⋂_{iϵI} \{R_{i} > min (\frac{(α + γ)}{m}, \frac{α}{m − 1})\} ,  (3)

if the R_{i}, iϵI, are in fact the largest m values R_{(i)}, i = n − m + 1, ⋯, n , then I = K(n−m+1), and γ = γ(n−m+1), so H_{i}, iϵI , will be accepted by the procedure, by at most the (n − m + 1)^{th} Step.

If the R_{i}, iϵI, are not the largest m values, then the smallest value of R_{i}, iϵI , is R_{(n−m+1−k)} for some k, 1 ≤ k ≤ n − m, and on account of (3) either

(a) R_{(n−m+1−k)} > (α + γ)/m; or (b) R_{(n−m+1−k)} > α/(m − 1) .

If (a), then R_{(n−m+1−k)} > α/m ≥ α/(m + k − 1) = α/(n − (n − m − k + 1)), so H_{(i)}, n − m − k + 1 ≤ i ≤ n , and so H_{i}, iϵI, are all accepted by at most the (n − m − k + 1)^{th} Step. If (b), then R_{(n−m+1−k)} > α/(m − 1) > α/m ≥ α/(m + k − 1), with the same conclusion. □ THEOREM: If the set {H_{i}, iϵI} is the set of true hypotheses (where I is any non-null subset of {1, ⋯, n}), then (1) holds for PROCEDURE CS.

\underline{Proof.} For notational convenience put A = \{H_{s}, sϵI, true\} . P (H_{s}, sϵI, are accepted |A)

≥ P (⋂_{iϵI} \{R_{i} > min (\frac{(α + γ)}{m}, \frac{α}{m − 1})\} |A)

by the Lemma;

= 1 − P (⋃_{iϵI} \{R_{i} ≤ min (\frac{(α + γ)}{m}, \frac{α}{m − 1})\} |A)

≥ 1 − ∑_{iϵI} P (R_{i} ≤ min (\frac{(α + γ)}{m}, \frac{α}{m − 1}) |A)

+ max(iϵI) ∑_{jϵI−\{i\}} P (R_{i} ≤ min (\frac{(α + γ)}{m}, \frac{α}{m − 1}) , R_{j} ≤ min (\frac{(α + γ)}{m}, \frac{α}{m − 1}) |A)

by Kounias’ inequality;

≥ 1 − min((α + γ), mα/(m − 1)) + max(iϵI) ∑_{jϵI−\{i\}} P (R_{i} ≤ \frac{α}{m}, R_{j} ≤ \frac{α}{m}|A)

by the uniformity of distribution of each R_{i}ϵI, since min((α + γ)/m, α/(m − 1)) < 1; = 1 − min((α + γ), mα/(m − 1)) + γ ≥ 1 − (α + γ) + γ = 1 − α .

COROLLARY. For any given I , if γ ≥ α/(m − 1), then P (H_{s}, sϵI, are accepted |A) ≥ 1 − mα/(m − 1) + γ . □

It will be seen that the test procedure could have been constructed, with increase of complexity at (2), round the sharper inequality of Hunter (1976), of which there are non-graph theoretic formulations (Margaritescu, 1986; Seneta, 1988). Possible use of degree 2 inequalities such as ours is mentioned in this context by Simes (1986, Section 5) and Shaffer (1986), and use of them for confidence intervals is made by Stoline (1983) (see Seneta, 1993).

The quantity γ defined in the LEMMA above is required for the proof of the THEOREM following it, which establishes strong control of the adaptive test procedure, but is not needed in the adaptive test procedure itself. Note that γ ≤ (m − 1)α/m ≤ α.

3. EXAMPLE AND SIMULATIONS

We shall measure power by

P (Reject at least one H_{i}, i = 1, ⋯, n).

This has the advantage that when all of the H_{i}, i = 1, ⋯, n hold, from (1) this value will be ≤ α, and its closeness to the nominal error α will measure the actual conservativeness of the error rate. According to the COROLLARY above, if we take I = {1, 2, ⋯, n} then γ > α/(n − 1) results in a bound < α . This suggests that the degree of conservativeness of the PROCEDURE CS is related to the strength of positive association between the R_{i}’s (and hence of X_{i}’s ) from the definition of γ . This is confirmed by Table 1 below.

We take the test statistics to be exchangeable under corresponding null hypotheses, so from (2)

γ(p) = (n − p)P (R_{1} ≤ \frac{α}{n − p + 1}, R_{2} ≤ \frac{α}{n − p + 1}),  (4)

which is thus non-random in this special setting. More specifically we consider uppertail tests where X_{i} = |T_{i}|, i = 1, 2, ⋯, n with T_{1}, T_{2}, ⋯, T_{n} defined by T_{i} = W_{i}/\sqrt{χ_{ν}^{2}/ν}, i = 1, 2, ⋯, n where the W_{i}’s are multivariate normal with EW_{i} = μ_{i}, Var(W_{i}) = 1, i = 1, 2, ⋯, n , Corr(W_{i}, W_{j}) = ρ, i ≠ j, and are independently distributed of χ_{ν}^{2} . Thus under H_{i} : μ_{i} = 0, i = 1, ⋯, n, the T_{i}, i = 1, ⋯, n have jointly a multivariate exchangeable t distribution with parameters n, ρ(ρ ≥ −1/(n − 1)), ν as in Dunnett’s tests (Miller, 1981). We take n = 3, ν = 16, α = 0.05 , and consider 0 ≤ ρ ≤ 1 . We can calculate from tables giving upper-tail values P (T_{1} ≤ a, T_{2} ≤ a) for various a and ρ = 0, ±0.1, ⋯, ±0.9 (Krishnaiah, Armitage and Breiter, 1969) our values of Δ(1) = (α + γ(1))/n. Some of these are shown in Table 1. Notice that in our setting Δ(1) > α/2 = 0.025 at ρ = 1 , but is < 0.025 for ρ ≤ 0.9 . Our measure of power when ρ ≤ 0.9 (in fact for ρ upto approximately 0.95) is thus P (R_{(1)} ≤ Δ(1)), whereas for the Bonferroni-Holm procedure the measure of power is P (R_{(1)} ≤ α/n) , which is smaller; and remains smaller than our measure of power for ρ very close to 1 viz. P (R_{(1)} ≤ α/(n − 1)). (A more sensitive measure of power would separate out Holm from Bonferroni.) Table 2 displays the power at ρ = 0.9 when μ_{1} = δ, μ_{2} = 2|δ|, μ_{3} = 3|δ| for the Bonferroni-Holm, Hochberg and CS procedures.

Table 1. Values of Δ(1) and Error Rate (E.R.) (n = 3, ν = 16, α = 0.05)

ρ  0  0.5  0.8  0.9  1
Δ(1)  0.0171 0.0183 0.0211 0.0228 0.0278
E.R. (α/n) 0.049 0.042 0.033 0.028
E.R. (Δ(1)) 0.050 0.046 0.042 0.040

Table 2. Power at ρ = 0.9 (n = 3, ν = 16, α = 0.05)

δ  -1  -0.5  0  0.5  1
α/n  0.685 0.184 0.028 0.159 0.633
Hochberg 0.689 0.187 0.034 0.164 0.634
Δ(1)  0.749 0.227 0.040 0.193 0.685

The error rate (E.R.) entries in Table 1 were produced from a simulation of 20,000 independent sets of values of the triple T_{i}, i = 1, 2, 3 at each ρ. These values for ρ = .9 are given again in Table 2 at δ = 0. The other values of Table 2 were also produced from 20,000 triples.

Overall, the simulations support a conclusion that our proposed procedure is most effective as regards power when test statistics are strongly positively dependent. The error rate is closer to the nominal value α irrespective of degree of dependence, and is not much affected by it. The indication is that Procedure CS controls error rate well, and has significantly better power than Hochberg.

While our procedure may be useful only for small n (1) holds without any restriction on the continuous joint distribution of test statistics.

Finally, our computational results on E.R. are consistent with those of Sarkar and Chang (1997, Table 2), inasmuch as

P (R_{(i)} ≥ \frac{α}{n − i + 1}, 1 ≤ i ≤ n) ≥ P (R_{(i)} ≥ \frac{iα}{n}, 1 ≤ i ≤ n)

ACKNOWLEDGEMENT

E.S. wishes to thank Harry Khamis and Barbara Mann for helping to motivate this work. The simulations were mostly done by Greta Lai in the course of an MA project during 1997.

REFERENCES

HOCHBERG, Y. (1988). A sharper Bonferroni procedure for multiple tests of

significance. Biometrika 75, 800-802. HOCHBERG, Y. and ROM, D. (1995). Extensions of multiple testing

procedures based on Simes’ test. J. Statist. Plan. Inf. 48, 141-152. HOCHBERG, Y. and TAMHANE, A.C. (1987). Multiple Comparison Procedures.

New York: Wiley. HOLM, S. (1979). A simple sequentially rejective multiple test procedure.

Scandinavian Journal of Statistics 6, 65-70. HOMMEL, G. (1988). A stagewise rejective multiple test procedure based on

a modified Bonferroni test. Biometrika 75, 383-386. HOMMEL, G. (1989). A comparison of two modified Bonferroni procedures.

Biometrika 76, 624-625. HUNTER, D. (1976). An upper bound for the probability of a union.

Journal of Applied Probability 13, 597-603. KOUNIAS, E.G. (1968). Bounds for the probability of a union of events, with

applications. Ann. Math. Statist. 39, 2154-2158. KRISHNAIAH, P.R., ARMITAGE, J.V. and BREITER, M.C. (1969). Tables for the

Probability Integrals of the Bivariate t Distribution. ARL 69-0060.

Applied Mathematics Research laboratory, Aerospace Research Laboratories, Ohio: Wright Patterson AFB.

MARGARITESCU, E. (1986). A note on Bonferroni’s inequalities.

Biometrical Journal 28, 937-943. MILLER, R.G. (1981). Simultaneous Statistical Inference, 2nd Edition,

New York: Springer. ROM, D.M. (1990). A sequentially rejective test procedure based on a modified

Bonferroni inequality. Biometrika 77, 663-666.

SAMUEL-CAHN, E. (1996). Is Simes improved Bonferroni procedure conservative?

Biometrika 83, 928-933. SARKAR, S.K., (1998). Some probability inequalities for ordered MT P_{2}

random variables: a proof of the Simes conjecture. Ann. Statist., 26, 494-504. SARKAR, S.K. and CHANG,C.K. (1997) The Simes method for multiple hypothesis testing with positively dependent test statistics. J. Amer. Statist. Assoc., 92,1601-1608. SENETA, E. (1988). Degree, iteration and permutation in improving

Bonferroni-type bounds. Australian J. Statist. 30A, 27-38. SENETA, E. (1993). Probability Inequalities and Dunnett’s Test. In F.M. Hoppe,

Ed: Multiple Comparisons, Selection, and Applications in Biometry

Chapter 3, pp 29-45. New York: Marcel Dekker. SENETA, E. and CHEN, T.H. (1997). A sequentially rejective test procedure.

Theory of Stochastic Processes 4 (to appear). SHAFFER, J.P. (1986). Modified sequentially rejective multiple test procedure.

J. American Statistical Association 81, 826-831. SIMES, R.J. (1986). An improved Bonferroni procedure for multiple tests of

significance. Biometrika 73, 751-754. STOLINE, M.R. (1983). The Hunter method of simultaneous inference and its

recommended use for applications having large known correlation structures, Journal of the American Statistical Association 78, 366-370.

WRIGHT, S.P. (1992). Adjusted P-Values for simultaneous inference.

Biometrics 48, 1005-1013.