On Familywise Error Rate Cutoffs under Pairwise Exchangeability

Abstract

In a pairwise exchangeable dependence setting for test statistics, the cutoffs of Sarkar et al. (2016) may be viewed as a first iteration improvement of Holm (1979)’s classical cutoffs under a convexity condition on the copula. The cutoffs of Seneta and Chen (1997) which improve Holm’s in the present exchangeability setting, are shown, after an analogous first iteration step, to lead to a refinement of Sarkar et al. (2016). Further, we show that the convexity condition can be circumvented in practice, computationally. Improvement by iteration limit of cutoffs is considered for both procedures. Comparisons between the effects of the several cutoff sets are made by way of plots of the familywise error rate against correlation $\rho$ in the classic setting of the multivariate Normal; and the distributional setting of the multivariate Generalized Hyperbolic for the important Variance Gamma type subfamily, for which a convexity condition cannot be analytically verified.

1 Introduction

1.1 Background

Seneta and Chen (2005) improved Holm’s step-down procedure by incorporating bivariate distribution information of the test statistics. This procedure (the CS procedure) is more powerful than Holm’s while still strongly controlling the familywise error rate at prespecified level $\alpha$ under arbitrary dependence structure. The cutoffs (critical values) for this procedure are of form

$$\begin{aligned} \min \, \left( \frac{\alpha }{n-i}, \, \frac{\alpha + \gamma _i}{n -i+1} \right) , \quad i=1, 2, \ldots ,n, \end{aligned}$$

for a set of size n of test statistics, where $\gamma _i$ incorporates the bivariate distributional information.

Numerical and simulation comparisons with other procedures by later authors (Sarkar et al. 2016; Lu 2014) have been implemented under the assumption of pairwise exchangeability, which considerably simplifies the expression for $\gamma _i$, while retaining the minimization constraint. We first indicate that under pairwise exchangeability the larger cutoffs $\frac{\alpha + \gamma _i}{n -i+1}$, $i = 1,2, \ldots ,n$, are valid (Seneta and Chen 1997; Lu 2014, Procedure III). The one-step iterative improvement of these larger cutoffs underpins the development of this paper.

The paper was motivated primarily by Sarkar et al. (2016), and its Supplementary Materials. These authors (see also Fu 2015) obtain a substantial improvement in the pairwise exchangeable setting, under a convexity condition, of Holm (1979)’s classical cutoffs, incorporating the bivariate distributional information as in Seneta and Chen (2005).

The stream of work of Sarkar and his collaborators on the Simes method for controlling the probability of Type I error is summarized in Sarkar (2008). Sarkar and Chang (1997) consider equicorrelated examples for a class of positively dependent multivariate distributions, and these are encompassed in our exchangeability setting.

We now formalize some of the above.

1.2 The Step-Down Procedure

We consider the multiple test problem when there are n hypotheses $H_1, H_2, \ldots , H_n$, and corresponding P-values $R_1, R_2, \ldots R_n,$ assuming the test statistics $X_1, X_2, \ldots , X_n$ are from a continuous distribution. Suppose that in such a multiple test procedure the property:

$$\begin{aligned} Pr(H_s, s \in I_{m},\text {are accepted}\mid H_s, s \in I_{m} \, \text {are true}) \ge 1- \alpha \end{aligned}$$

(1)

holds for prespecified size of test $\alpha$ (familywise error rate: FWE) for every $I_{m}$, where $I_{m}$ is any non-null subset of $\{1,2, \ldots , n \}$, for every $m = 1,2, \ldots , n.$ Then the FWE is said to be strongly controlled.

Let $R_{(1)}, R_{(2)}, \ldots , R_{(n)}$ be the ordered P-values, and $\Delta (i)$, $1 \le i \le n$ be a strictly increasing sequence of constants, with $0< \Delta (i)< 1,$ for each i. A step-down procedure begins by testing if $R_{(1)} < \Delta (1).$ If so, reject $H_{(1)}$ and go on to test if $R_{(2)} < \Delta (2).$ If not, accept all hypotheses. In general, if $R_{(i)} < \Delta (i), 1 \le i \le j-1$, then at step j the remaining hypotheses are $H_{(j)}, H_{(j+1)}, \ldots , H_{(n)},$ and the inequality next to be checked is $R_{(j)}< \Delta (j)$. If it holds, reject $H_{(j)}$ and continue. Otherwise accept $H_{(j)}, H_{(j+1)}, \ldots , H_{(n)}.$ The process may continue until a decision is made on the basis of whether $R_{(n)}< \Delta (n).$ The step-down procedure of Holm (1979) uses the set of constants

$$\begin{aligned} {\bar{\Delta }}(i) =\frac{\alpha }{n-i+1}, \, \, 1 \le i \le n, \end{aligned}$$

(2)

and Holm proved that with these constants the FWE is strongly controlled.

1.3 The Exchangeability Context

In the sequel we assume that $(X_i, X_j), i \ne j, i,j \in I_{m}$ are exchangeable (so that $(R_i, R_j)$ are also), when $I_{m}, m \ge 2$ is the index set of assumed true hypotheses, and that their joint bivariate distribution is the same for each such $I_{m}.$ The marginal distribution of each $R_i, i \in I_{m}, m \ge 1$ is U(0, 1). We introduce the notation $H(u) = Pr(R_1 \le u, R_2 \le u) , u \in (0,1)$, and $n_0 = n - i +1.$ Seneta and Chen (1997), Sections 4 and 5 showed that then the critical cutoff constants:

$$\begin{aligned} {\tilde{\Delta }}(i) = \frac{\alpha + (n_0 -1) H\left( \frac{\alpha }{n_0}\right) }{n_0}, \end{aligned}$$

(3)

are monotonic increasing with i (that is: decreasing with increasing $n_0$), and the step-down procedure based on them strongly controls the FWE. Moreover the step-down procedure with these constants provides tighter control on the FWE than Holm’s since

$$\begin{aligned} {\tilde{\Delta }}(i) > {\bar{\Delta }}(i) = \frac{\alpha }{n_0}, \quad 1 \le i \le n-1; \quad {\tilde{\Delta }}(n) = {\bar{\Delta }}(n) = \alpha . \end{aligned}$$

Note that $\frac{\alpha }{n_0}$ are Holm’s constants, and that $H(u) = Pr( \max (R_1, R_2) \le u )$, so that H(u), $0< u < 1$, is a cdf. It is also the copula $C(u_1, u_2) = Pr(R_1 \le u_1, R_2 \le u_2)$, $0< u_1, u_2 <1$, of the bivariate joint distribution of the exchangeable pair $(X_1, X_2)$, evaluated on the diagonal $u_1 = u_2 = u$.

Sarkar et al. (2016), motivated by Seneta and Chen (2005), make the overarching Assumption: H(u) is convex in $u \in (0,1)$, and consequently show that the monotonic increasing cut-offs:

$$\begin{aligned} {\hat{\Delta }}(i) = \frac{\alpha ^2/ n_0}{G_{n_0}(\alpha /n_0)}, \end{aligned}$$

(4)

where

$$\begin{aligned} G_{n_0}(u) = n_0 u - (n_0 -1) H(u), \quad 0< u < 1, \end{aligned}$$

(5)

provide tighter control on the FWE than

$$\begin{aligned} \Delta _{CS}(i)= \min \left\{ \frac{\alpha }{n_0 - 1},\, \frac{\alpha + (n_0 -1) H(\frac{\alpha }{n_0})}{n_0}\right\} . \end{aligned}$$

(6)

(which is a specialization to the exchangeable setting from a general context in Seneta and Chen (1997)) since

$$\begin{aligned} {\hat{\Delta }}(i) > {\Delta _{CS}(i)}, \quad 1 \le i \le n-1; \quad {\hat{\Delta }}(n) = \Delta _{CS}(n) = \alpha . \end{aligned}$$

The paper Seneta and Chen (1997) in which the cutoffs (3) are justified in the exchangeable setting, was cited and used in Seneta and Chen (2005), but was not digitally available till now¹, and hard-copy journal access may have been difficult. Sarkar et al. (2016)’s argument in any case implies that $\hat{\Delta }(i) > \tilde{\Delta }(i)$, $1 \le i \le n-1$.

2 General Formulation of Cutoffs

In fact, any sequence $w_1(n_0)$, $n_0=1,2,\ldots ,n$ satisfying $w_1(1) = \alpha$ and for $n_0 \ge 2$

$$\begin{aligned} 0< w_1(n_0)< \alpha , \quad w_1(n_0)< w_1(n_0 -1), \quad G_{n_0}(w_1(n_0)) < \alpha , \end{aligned}$$

(7)

where $G_{n_0}(u)$ is defined by (5) and $H(u) = 0, u \in (0,1)$ or $H(u) = Pr(R_1 \le u, R_2 \le u)$, $u \in (0,1)$ provides a set of monotonic increasing cut-offs which assure strong control of the FWE in the step-down test procedure, by setting $\Delta (i)= w_1(n-i+1)$, $i = 1, 2, \ldots , n.$

This statement is implicit in Sarkar et al. (2016)’s reformulation of the justification in Seneta and Chen (2005), when specialized to the exchangeable setting.

The cutoffs (2), (3) satisfy the description of a sequence $w_1(n_0)$, $n_0=1,2,\ldots ,n$, when $H(0)=0$, $u \in (0,1)$. The cutoffs (4) satisfy the description of a sequence $w_2(n_0)$, $n_0=1,2, \ldots , n$, where $w_1(n_0)=\frac{\alpha }{n_0}$, that is: (2), according to the following theorem, whose proof is given in Appendix.

Theorem 1

Suppose $\alpha$, $0< \alpha <1$, is a fixed constant. Define

$$\begin{aligned} G_{n_0}(u)=n_0 u - (n_0-1)H(u), \quad 0 \le u \le \alpha , \end{aligned}$$

where $H(u), 0 \le u \le \alpha$ is any function satisfying $H(0)=0$, $0 \le H(u)/u<1$, $0<u \le \alpha$. Let $w_1(n_0)$, $n_0=1,2,\ldots , n$ satisfy $w_1(1) = \alpha$ and (7) for $n_0 \ge 2$. Define for $n_0 = 1,2,\ldots , n$:

$$\begin{aligned} w_2(n_0) {\mathop {=}\limits ^{def}} \frac{\alpha w_1(n_0)}{G_{n_0}(w_1(n_0))}. \end{aligned}$$

(8)

Then

$$\begin{aligned} w_2(1) = \alpha ; \quad \text { and for } \quad n_0 \ge 2, \quad 0< w_1(n_0)< w_2(n_0) < \alpha . \end{aligned}$$

(9)

If, additionally we assume $H(u)/u \uparrow$, $u \uparrow$, then for $n_0 \ge 2$, additionally:

$$\begin{aligned} w_2(n_0)< w_2(n_0 -1), \quad G_{n_0}(w_2(n_0)) < \alpha . \end{aligned}$$

(10)

If $w_1^{*}(n_0) > w_1(n_0)$, $n_{0} \ge 2$ and $w_1^{*}(n_0)$ satisfies the same conditions as $w_1(n_0)$, and if $H(u)/u \uparrow$, $u \uparrow$, then

$$\begin{aligned} w_2^{*}(n_0) > w_2(n_0), \quad n_0 \ge 2. \end{aligned}$$

(11)

In both Sarkar et al. (2016) and in its Supplementary Materials, much attention is focused on the problem of verifying the assumption of convexity of H(u), $0< u < 1$. According to the Theorem, the conditions for the copula H(u), to obtain the same conclusion, are less restrictive.

More to the point, the validity of (10) and (11) can be checked directly from the values of $w_2(n_0)$, and the values of {$w_i(n_0)$, $w_i^{*}(n_0)$, $i=1, 2$} respectively in a testing setting, without knowledge such as $H(u)/u \uparrow$, $u \uparrow$, which removes the dependence on such a condition in statistical application. We make use of this in Section 4, in an example when the condition $H(u)/u \uparrow$, $u \uparrow$ cannot be analytically verified.

Finally, under a condition such as $H(u)/u \uparrow$, $u \uparrow$, there will be an iteration limit, $w_L(n_0)$ to the iteration sequence $\{w_i(n_0), i \ge 1 \},$ where $w_{i+1} =\alpha w_i(n_0)/G_{n_0}(w_i(n_0))$ on account of the monotonicity $0< w_i(n_0)< w_{i+1}(n_0) < \alpha .$ Likewise $w^{*}_L(n_0)$ for the analogous sequence $\{w^{*}_i(n_0), i \ge 1 \},$ with $w^{*}_L(n_0) \ge w_L(n_0).$ Both limits will satisfy the equation $G_{n_0}(w(n_0))= \alpha , 0< w(n_0)< \alpha$. For our starting points (2), (3) of the two sequences the limit values appear to be coincident in the cases we have considered. In practice, a further iteration of the sequence with starting point (3) will produce slightly higher cutoffs than both (4) and (8).

For comparisons of the differential effects of various cutoffs, we use the familywise error rate (FWER) as a function of $\rho$ as the means of comparison. This is in general, and here, defined as

$$Pr(\mathrm{Reject \;at \;least \;one} \;H_{i}, \;{i}= 1,2, \ldots , {n}\; \mid {H}_{i}, {i}= 1,2, \ldots , \;n\; \;\mathrm{are \;all \;true}).$$

Thus for the step-down procedure

$$\begin{aligned} \begin{aligned} \text {FWER}&= Pr(R_{(1)}< \Delta (1) \mid H_i, i= 1,2, \ldots , n \ \text {are all true}) \\&= Pr(R_{(1)}< w_j(n) \mid H_i, i= 1,2, \ldots , n \ \text {are all true}) \\&= Pr\left( \min _{i=1,2, \ldots , n}X_i < F_X^{-1}(\Delta (1)) \mid H_i, i= 1,2, \ldots , n \, \text {are all true} \right) \end{aligned} \end{aligned}$$

(12)

for fixed $j=1,2$, in terms of sample values $X_i$, $i=1,2, \ldots , n$, where $F_X(x)$ is the marginal cdf of $X_i$. The expression (12) may be estimated from repeated simulation of the sample $X_i$, $i=1,2, \ldots , n$, once $F_X^{-1}(\Delta (1))$ is numerically calculated.

3 Testing Comparisons for the Multivariate Normal

Seneta and Chen (2005), Section 4, provide a tabulation, Table 4, with $\alpha = 0.05$ of (3) for $n_0 = 1,2, \ldots , 8$ against $\rho = 0.0$, 0.1, $\ldots$, 0.9, 0.95, 0.99, 1.00 for test statistics from the multivariate normal for a two-tailed test for zero mean, and investigate familywise error rate and power in this setting when $n=3$, by simulation. This was for a two-tailed test, capitalizing on the symmetry about zero of the marginal normal distribution under the null hypothesis. Here $R_i = 2\Phi _X(- |X_i |)$, where $\Phi _X(x)$ is the cdf of standard normal, and $X_i$, $i =1,2,\ldots , n$ is the multivariate normal sample with pairwise correlations $\rho$. Then

$$\begin{aligned} H\left( \frac{\alpha }{n_0}\right) = Pr( |X_1 |> t_0, |X_2 |> t_0), \end{aligned}$$

(13)

where $t_0= \Phi _X^{-1}\left( 1- \frac{\alpha }{2n_0}\right)$, and $H(u) = Pr(R_1 \le u, R_2 \le u)$ as before.

Sarkar et al. (2016) compare graphically the effect of their cutoffs (4) in the two-tail setting but with $n=16$, and not with (3), but with the less effective (6). Thus a graphical comparison of the effects of all 3 sets of cutoff constants (3), (6), and (4), with increasing $\rho$, as regards error rate in the multivariate normal setting is one of our aims. We shall do so with $n=16$, and for a lower-tailed test setting only. In this setting more simply $R_i = \Phi _X^{-1}(X_i)$.

In the notation of the Theorem, we aim at the comparison of the effect of the Sarkar et al.’s cut offs $w_2(n_0)$ given by (4) with the effect of the cutoffs $w_1(n_0)$ given by (3), and of the corresponding $w_2(n_0)$ cutoffs given by (8). The cut offs are presented numerically, for two values of $\rho$, in the first two columns in the body of Table 1. Since clearly (2) < (3), the corresponding sequences of cutoffs have the same inequality relation, as the Theorem asserts. Consequently the (8) cutoff values deriving from $w_1(n_0)$ and given by (3) provide tighter control on the FWE for each given $\rho$.

The third column ($w_2'(n_0)$) in Table 1 is the iteration of the second ($w_2(n_0)$) column. For these two values of $\rho$, and for lower values, the values in the these two columns coincide to 4 decimal places. Table 2 reports $w_2'(n_0)$ to the accuracy shown in the table. The columns up to $\rho = 0.8$ also give the values of $w_2(n_0).$ Below, in Fig. 2, we compare graphically $\hat{\Delta }(i)$ and $w_2(n_0)$ given by (4) and (8) respectively. Apart from high values of $\rho$, the two plots of FWER are indistinguishable, which suggests that in practice it is advisable to use the values of the next iteration, $w_2'(n_0)$ in general. This analytical iteration step requires little more computational effort.

Table 1 Comparing the cutoffs under multivariate Normal with $n=16$ for $\rho = 0.5$ & 0.8

Table 2 Cutoffs $w_2'(n_0)$ with various $\rho$ and $n=16$ for the multivariate normal

A plot of the FWER, defined for a given $\rho$ by (12), against $\rho$ for each of the cut-off sets, provides visual comparison. All simulations were done with the R software package (R Core Team 2021). In Fig. 1, we plotted the error rate in the multivariate normal setting at level $\alpha = 0.05$ with $n = 16$ for various values of common correlation $\rho$ and three cutoff constants $\tilde{\Delta }(i)$, $\hat{\Delta }(i)$ and $\Delta _{CS}(i)$ in (3), (4) and (6) respectively. One hundred thousand (100,000) independent replications were used in all simulations. The smoothed line is provided by a generalised additive model (GAM) fit. Figure 1 parallels Fig. 1a of Sarkar et al. (2016) under our lower tail setting. We can see how $\hat{\Delta }(i)$ has tighter control of FWE than the other two cutoff constants.

In Fig. 2, we plotted the error rate in the multivariate normal setting for various values of common correlation $\rho$ for cutoff constants $\hat{\Delta }(i)$ and the proposed $w_2(n_0)$ as in (8) under the same setting as before. We can see that $w_2(n_0)$ provides tighter control of FWE over $\hat{\Delta }(i)$ for high values of $\rho$ but the performance is indistinguishable otherwise.

Fig. 1

Comparing the familywise error rate with cutoff constants $\Delta _{CS}(i)$, $\tilde{\Delta }(i)$ and $\hat{\Delta }(i)$ given by (6), (3) and (4) respectively for various values of common correlation $\rho$ at level $\alpha = 0.05$ for the multivariate normal distribution with $n=16$

Fig. 2

Comparing the familywise error rate with cutoff constants $\hat{\Delta }(i)$ and $w_2(n_0)$ given by (4) and (8) respectively for various values of common correlation $\rho$ at level $\alpha = 0.05$ for the multivariate normal distribution with $n=16$

While we computed the FWER in this section using simulation once cutoffs have been numerically calculated, similarly to Sarkar et al. (2016) and Lu (2014), we remark that, alternatively, the FWER can be calculated with numerical integration in this multivariate normal setting, by adopting Gupta et al. (1973).

4 The Generalized Hyperbolic Distribution Tests

The multivariate skew GH distribution, $GH(0,R_n,{\theta }, p,a,b),$ which we take as the joint distribution of our test-statistics $X = (X_1, X_2, \ldots ,X_n)^{\top }$ to have is defined by its mean-variance mixing representation as

$$\begin{aligned} X = {\mu }+ {\theta } W + \sqrt{W}Z \quad \end{aligned}$$

(14)

where ${\mu } = (\mu _1, \mu _2, \ldots ,\mu _n)^{\top }, {\theta } = (\theta _1,\theta _2, \ldots , \theta _n)^{\top }$, and $W \sim GIG(p,a,b)$ is independently distributed of ${Z} \sim N(0, R_n)$. Here $R_n =(1-\rho ) {I}_n + \rho 1_n 1_n^{\top }$, with $-1<\rho <1$, so that $R_2 = \left( {\begin{matrix} 1 &{} &{} \rho \\ \\ \rho &{} &{} 1\end{matrix}}\right)$. We assume that the random variable W has a (univariate) Generalised Inverse Gaussian (GIG) distribution, denoted by GIG(p, a, b), that is, it has density

$$\begin{aligned} \begin{aligned} f_{GIG}(w)&= \frac{1}{2\overline{K}_{p}(a,b)}w^{p-1}\exp (-\frac{1}{2}(a^2w^{-1}+b^2w)), \quad w>0; \\&= 0, \quad \text {otherwise;} \end{aligned} \end{aligned}$$

where

$$\begin{aligned} \overline{K}_{p}(a,b) = {\left\{ \begin{array}{ll} (\frac{a}{b})^{p}K_{p}(ab), &{}{p} \in \mathbb {R}, \;\mathrm{if} \;a,b>0;\\ b^{-2p}\Gamma (p)2^{p-1}, &{}{p ,b}>0, \;\mathrm{if}\; a=0;\\ a^{2p}\Gamma (-p)2^{-p-1}, &{}{a}>0 \;\mathrm{and}\; {p} <0, \mathrm{if} \;b=0. \end{array}\right. } \end{aligned}$$

(15)

Here $K_{p}(\omega ), \ \omega >0,$ is the modified Bessel function of the second kind (Erdélyi et al. 1954) with index $p \in \mathbb {R}$. The multivariate GH distribution was first introduced but only briefly discussed in Barndorff-Nielsen (1977). Its application was considered in detail by Blaesild (1981).

We take for, $i=1,2,\ldots , n$, our ith null hypothesis as : $H_i: \mu _i = \theta _i = 0,$ so that when all $H_i$ are true,

$$\begin{aligned} {X} = \sqrt{W}{Z}, \end{aligned}$$

(16)

and similarly for any subset of indices $\{1,2, \ldots , n \},$ so that exchangeability holds, and the other assumptions of Section 1.3, The Exchangeability Context, hold. Then the distribution functions of the $X_i$ are the same: $F_i(u) = F_1(u)$. The marginal density of each of $X_i$ is then given by

$$\begin{aligned} f_{X_1}(x) = \frac{\overline{K}_{p-1/2}((x^2 +a^2)^{1/2}, b)}{{\sqrt{2\pi }\, {\overline{K}}_p(a,b)}},\quad x \in \mathbb {R}. \end{aligned}$$

(17)

The case $b=0$ of (16) encompasses the symmetric multivariate t-distribution, which has been studied as a central example in Sarkar et al. (2016). The case $b>0$ describes the Variance Gamma type family. The bivariate case ($n=2$), which is central, as we have seen, to the role of the copula H(u) in our exchangeable setting, has been studied extensively, by Fung and Seneta (2021), from a related standpoint, but we have not been able to establish the general property $H(u)/u \uparrow$, $u \uparrow$ in the exchangeable setting, unlike for the various examples considered by Sarkar et al. (2016) and Sarkar and Chang (1997).

The Variance Gamma distribution, which is a special case, when $a=0$, $b = \sqrt{\frac{2}{\nu }}$, $p=\frac{1}{2}$, is important in financial mathematics modelling (e.g. Fung and Seneta 2010). It is therefore important for purposes of iterative improvement, to see on examples of the Variance Gamma type family with parameters numerically specified, whether (10) and (11) check directly.

In Fig. 3, we plotted the error rate in the multivariate GH setting with $a = 0, p = 4,\;b = 0.5$, at level $\alpha = 0.05$ with $n = 16$ for various values of common correlation $\rho$ and two cutoff constants $\hat{\Delta }(i)$ and $\Delta _{CS}(i)$ in (4) and (6) respectively. One hundred thousand (100,000) independent replications were used in all simulations. The smoothed line is once again provided by a GAM fit. Similar to the multivariate Normal setting, we can see that $w_2(\eta _0)$ provided tighter control of FWE for high values of $\rho$ but performance was indistinguishable otherwise.

In Fig. 4, we used the cut off constant of $w_2(\eta _0)$ computed under the multivariate normal and the symmetric GH setting for various values of common correlation $\rho$. The parameters used for the GH are the same as before. We can see that GH provided a tighter control when the data are highly correlated. $w_2(n_0)$ provided a tighter control of FWE over $\hat{\Delta }(i)$ for high values of $\rho$.

Here we provide some additional information regarding the simulations we conducted in this section. To calculate the cutoffs in general under pairwise exchangeability, we need to calculate $H(u) = Pr(R_1 \le u, R_2 \le u) , u \in (0,1)$ accurately for a variety of values of u, which first requires accurate calculation of bivariate probabilities $P(X_1 \le x_1, X_2 \le x_2)$. Our solution was to convert to a univariate problem with

$$\begin{aligned} P(X_1 \le x_1, X_2 \le x_2) = \int ^{x_2}_{-\infty } P(X_1 \le x_1 \mid X_2 = y) f_{X_2}(y)\,dy, \end{aligned}$$

which can be evaluated with high precision with numerical integration as the multivariate GH is closed under conditioning (Blaesild 1981). We had tried the ghyp package of Weibel et al. (2020) that allows one to evaluate the joint cdf of multivariate GH but it uses the Monte Carlo method. Given we often need to compute the joint tail probability with high threshold, the Monte Carlo method leads to low accuracy even when the number of simulations we used were quite substantial. While the algorithm that calculates the multivariate normal joint cdf also has a Monte Carlo component, the values fluctuate a lot less there. While one can always increase the number of simulations to achieve the required accuracy, it is not feasible when combined with the number of replications we need to run to obtain the FWER. All the code used in the paper can be obtained by contacting the authors.

Fig. 3

Comparing the familywise error rate with cutoff constants $\hat{\Delta }(i)$ and $w_2(n_0)$ given by (4) and (8) respectively for various values of common correlation $\rho$ at level $\alpha = 0.05$ for the multivariate GH distribution with $n=16$

Fig. 4

Comparing the familywise error rate with cutoff constants $w_2(n_0)$ given by (8) for various values of common correlation $\rho$ at level $\alpha = 0.05$ for the multivariate normal and GH distributions with $n = 16$

Appendix: Proof of Theorem 1

Proof

That $w_2(1) = \alpha$ follows from $w_1(1)=\alpha$ and (7), recalling the definition of $G_{n_0}(u)$. That $w_2(n_0)< \alpha , \ n_0 \ge 2,$ follows from $H(w_1(n_0))/w_1(n_0) <1,$ which implies $G_{n_0}(w_1(n_0))/w_1(n_0) >1.$ That $w_1(n_0) < w_2(n_0)$ follows from $G_{n_0}(w_1(n_0)) <\alpha$ in (7).

Next, notice that $G_{n_0-1}(u)< G_{n_0}(u),\, 0 <u \le \alpha .$ This is tantamount to:

$$\begin{aligned} (n_0-1) u - (n_0-2)H(u)< n_0 u - (n_0-1)H(u), \text {that is:}\,\, H(u)< u, \, 0 <u \le \alpha , \end{aligned}$$

which is true by assumption.

Then if, additionally, $H(u)/u \uparrow$ as $u \uparrow , 0 <u \le \alpha ,$ it follows that $G_{n_0}(u)/u \downarrow$ as $u \uparrow$, and since $G_{n_0-1}(u) < G_{n_0}(u),$

$$\begin{aligned} \begin{aligned} \frac{w_1(n_0-1)}{G_{n_0-1}(w_1(n_0-1))}&> \frac{w_1(n_0-1)}{G_{n_0}(w_1(n_0-1))} \\&> \frac{w_1(n_0)}{G_{n_0}(w_1(n_0))} \end{aligned} \end{aligned}$$

since by initial assumption (7) , $w_1(n_0) < w_1(n_0-1)$.

Thus from the definition (8):

$$\begin{aligned} \begin{aligned} w_2(n_0)=&\frac{\alpha w_1(n_0)}{G_{n_0}(w_1(n_0 ))}\\&<\frac{\alpha w_1(n_0-1)}{G_{n_0-1}(w_1(n_0-1))}= w_2(n_0-1). \end{aligned} \end{aligned}$$

Therefore, for fixed n,

$$\begin{aligned} w_2(n-(i-1)+1) < w_2(n-i+1), \end{aligned}$$

so $w_2(n-i+1)$ is increasing with i.

From (8),

$$\begin{aligned} G_{n_0}(w_2(n_0))= G_{n_0}(\frac{\alpha w_1(n_0)}{G_{n_0}(w_1(n_0))}). \end{aligned}$$

(18)

Since $G_{n_0}(u)/u\downarrow\mathrm{as}\,u\uparrow,0<u\leq\alpha$ if we continue to assume $H(u)/u \uparrow , u \uparrow$, it follows that

$$\begin{aligned} \begin{aligned} \frac{G_{n_0}(cu)/cu}{G_{n_0}(u)/u}< & {} 1, \, \text {for} \,\, c> 1, \text {providing}\, cu \le \alpha , \text {so that} \\ \frac{G_{n_0}(cu)}{G_{n_0}(u)}< & {} c, \text {for} \, c > 1, \text {providing}\, cu \le \alpha . \end{aligned} \end{aligned}$$

(19)

Now, taking

$$\begin{aligned} \begin{aligned} c= & {} \frac{\alpha }{G_{n_0}(w_1(n_0))}, u = w_1(n_0), \text {so that} \\ cu= & {} \frac{\alpha w_1(n_0)}{G_{n_0}(w_1(n_0))}= w_2(n_0), \end{aligned} \end{aligned}$$

it follows from (7) that $c>1$ and from (8) and the fact that $w_2(n_0) < \alpha$, that (19) amounts to: $G_{n_0}(w_2(n_0))< \alpha ,$ as asserted in (10).

Finally, equivalent to $w_2^{*}(n_0) > w_2(n_0)$ is

$$\begin{aligned} \begin{aligned} \frac{\alpha w_1^{*}(n_0)}{G_{n_0}(w_1^{*}(n_0))}> & {} \frac{\alpha w_1(n_0)}{G_{n_0}(w_1(n_0))}.\, \, \text { That is:}\\ \frac{G_{n_0}(w_1^{*}(n_0))}{w_1^{*}(n_0)}< & {} \frac{G_{n_0}(w_1(n_0))}{w_1(n_0)} \end{aligned} \end{aligned}$$

and this is so since $w_1^{*}(n_0) > w_1(n_0)$ by assumption, and $G_{n_0}(u)/u \downarrow , u \uparrow , 0 < u \le \alpha ,$ when $H(u)/u \uparrow , u \uparrow .$

This completes the proof of all parts of the theorem. $\Box$