On Familywise Error Rate Cutoffs under Pairwise Exchangeability

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 ρ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho$$\end{document} 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.


Background
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 under arbitrary dependence structure. The cutoffs (critical values) for this procedure are of form min n − i , + i n − i + 1 , i = 1, 2, … , n, 1 3 59 Page 2 of 13 for a set of size n of test statistics, where 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 i , while retaining the minimization constraint. We first indicate that under pairwise exchangeability the larger cutoffs + i n−i+1 , i = 1, 2, … , 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.

The Step-Down Procedure
We consider the multiple test problem when there are n hypotheses H 1 , H 2 , … , H n , and corresponding P-values R 1 , R 2 , … R n , assuming the test statistics X 1 , X 2 , … , X n are from a continuous distribution. Suppose that in such a multiple test procedure the property: holds for prespecified size of test (familywise error rate: FWE) for every I m , where I m is any non-null subset of {1, 2, … , n} , for every m = 1, 2, … , n. Then the FWE is said to be strongly controlled. Let R (1) , R (2) , … , R (n) be the ordered P-values, and Δ(i) , 1 ≤ i ≤ n be a strictly increasing sequence of constants, with 0 < Δ(i) < 1, for each i. A step-down procedure begins by testing if R (1) < Δ(1). If so, reject H (1) and go on to test if R (2) < Δ(2). If not, accept all hypotheses. In general, if R (i) < Δ(i), 1 ≤ i ≤ j − 1 , then at step j the remaining hypotheses are H (j) , H (j+1) , … , H (n) , and the inequality next to be checked is R (j) < Δ(j) . If it holds, reject H (j) and continue. Otherwise accept H (j) , H (j+1) , … , H (n) . The process may continue until a decision is made on the basis of whether R (n) < Δ(n). The stepdown procedure of Holm (1979) uses the set of constants and Holm proved that with these constants the FWE is strongly controlled.

The Exchangeability Context
In the sequel we assume that (X i , X j ), i ≠ j, i, j ∈ I m are exchangeable (so that (R i , R j ) are also), when I m , m ≥ 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 ∈ I m , m ≥ 1 is U(0, 1). We introduce the notation H(u) = Pr(R 1 ≤ u, R 2 ≤ u), u ∈ (0, 1) , and n 0 = n − i + 1. Seneta and Chen (1997), Sections 4 and 5 showed that then the critical cutoff constants: 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 Note that n 0 are Holm's constants, and that H(u) = Pr(max(R 1 , R 2 ) ≤ u) , so that H(u), 0 < u < 1 , is a cdf. It is also the copula C(u 1 , u 2 ) = Pr(R 1 ≤ u 1 , R 2 ≤ 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 ∈ (0, 1) , and consequently show that the monotonic increasing cut-offs: where provide tighter control on the FWE than (which is a specialization to the exchangeable setting from a general context in Seneta and Chen (1997)) since 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 1 , and hard-copy journal access may have been difficult. Sarkar et al. (2016)'s argument in any case implies that Δ (i) >Δ(i) , 1 ≤ i ≤ n − 1.
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 ↑ , u ↑ , 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 ↑ , u ↑ cannot be analytically verified.
Finally, under a condition such as H(u)∕u ↑ , u ↑ , there will be an iteration limit, . Both limits will satisfy the equation G n 0 (w(n 0 )) = , 0 < w(n 0 ) < . 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 as the means of comparison. This is in general, and here, defined as Thus for the step-down procedure 1)) is numerically calculated. Seneta and Chen (2005), Section 4, provide a tabulation, Table 4, with = 0.05 of (3) for n 0 = 1, 2, … , 8 against = 0.0 , 0.1, … , 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.

Testing Comparisons for the Multivariate Normal
is the cdf of standard normal, and X i , i = 1, 2, … , n is the multivariate normal sample with pairwise correlations . Then , and H(u) = Pr(R 1 ≤ u, R 2 ≤ 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 , 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 = Φ −1 X (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 , 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 .
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 , 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 = 0.8 also give the values of w 2 (n 0 ). Below, in Fig. 2, we compare graphically Δ (i) and w 2 (n 0 ) given by (4) and (8)
of , 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.
A plot of the FWER, defined for a given by (12), against for each of the cut-off sets, provides visual comparison. All simulations were done with the R software package (R Table 1 Comparing the cutoffs under multivariate Normal with n = 16 for = 0.5 & 0.8  Core Team 2021). In Fig. 1, we plotted the error rate in the multivariate normal setting at level = 0.05 with n = 16 for various values of common correlation and three cutoff constants Δ (i) , Δ (i) and Δ 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 Δ (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 for cutoff constants Δ (i) and the proposed w 2 (n 0 ) as in (8) under Fig. 1 Comparing the familywise error rate with cutoff constants Δ CS (i) , Δ (i) and Δ (i) given by (6), (3) and (4) respectively for various values of common correlation at level = 0.05 for the multivariate normal distribution with n = 16 Fig. 2 Comparing the familywise error rate with cutoff constants Δ (i) and w 2 (n 0 ) given by (4) and (8) respectively for various values of common correlation at level = 0.05 for the multivariate normal distribution with n = 16 the same setting as before. We can see that w 2 (n 0 ) provides tighter control of FWE over Δ (i) for high values of but the performance is indistinguishable otherwise.
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).

The Generalized Hyperbolic Distribution Tests
The multivariate skew GH distribution, GH (0, R n , , p, a, b), which we take as the joint distribution of our test-statistics X = (X 1 , X 2 , … , X n ) ⊤ to have is defined by its mean-variance mixing representation as 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 where Here K p ( ), > 0, is the modified Bessel function of the second kind (Erdélyi et al. 1954) with index p ∈ ℝ . 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, … , n , our ith null hypothesis as : H i ∶ i = i = 0, so that when all H i are true, and similarly for any subset of indices {1, 2, … , 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 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 ↑, u ↑ 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 = √ 2 , p = 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 = 0.05 with n = 16 for various values of common correlation and two cutoff constants Δ (i) and Δ 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 ( 0 ) provided tighter control of FWE for high values of but performance was indistinguishable otherwise.
In Fig. 4, we used the cut off constant of w 2 ( 0 ) computed under the multivariate normal and the symmetric GH setting for various values of common correlation . 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 Δ (i) for high values of .
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 ≤ u, R 2 ≤ u), u ∈ (0, 1) accurately for a variety of values of u, which first requires accurate calculation of bivariate probabilities P(X 1 ≤ x 1 , X 2 ≤ x 2 ) . Our solution was to convert to a univariate problem with Fig. 3 Comparing the familywise error rate with cutoff constants Δ (i) and w 2 (n 0 ) given by (4) and (8) respectively for various values of common correlation at level = 0.05 for the multivariate GH distribution with n = 16 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.