Denote the paired differences of the sample vectors X→1,X→2,…,X→m,Y→1,Y→2,…,Y→m described above, by D→i=X→i−Y→i, for i = 1, 2, …, m. Let us assume that the vectors D→1,D→2,…,D→m are independent and identically distributed (iid) with multivariate normal distribution, i.e., D→i~N(μ→,Σ), for i = 1, 2, …, m, where μ→ is the n-dimensional vector (μ₁, μ₂, …, μ_n) corresponding to the mean difference for each of the n features, whereas Σ is a n by n correlation matrix associated with the dependence between the n tested features.

Kim et al. [27] use block matrices to model the correlation matrix in its multivariate Gaussian model, which represents gene concentration levels. This structure assumes that groups of genes within the same regulatory pathway are correlated, while genes in different pathways are uncorrelated. Each block within the matrix corresponds to a distinct gene regulatory pathway, with specific variance and correlation coefficients defining the relationships among genes within that pathway. Motivated by the correlation form in [27], let us assume that Σ is a block matrix which can be expressed as

i.e., Σ has k blocks A₁, A₂, …, A_k. Note that in order to obtain such structure of the dependence matrix Σ, the tested features should be ordered (labeled) in a correct manner based on some prior information, e.g., in genetics one can apply knowledge from previous studies on the gene regulatory pathways. Dang et al. [28] use an interactive adjacency matrix to visualize binary relationships between proteins in biological pathways. To simplify complex networks, proteins with similar interaction patterns are grouped together, and these groups can then be “collapsed” into single nodes, effectively forming block matrices that represent relationships between groups rather than individuals. In case the number k and sizes of the blocks in Σ are unknown, Perreault et al. [29] propose an algorithm based on Kendall's rank correlation for identifying different disjoint groups in Σ when the features in each block are equicorrelated. This block structure with exchangeable variables, i.e., with constant correlation for any two features within a block, is also considered by Hartung [23] in the context of combining dependent statistics. In our study, we propose the following two generalizations of such dependence matrices, which allow the correlation coefficient in a block to decrease in an exponential or a linear manner.

For simplicity and clarity of the results in our simulation analysis, in Section 4 we assume that Σ has k blocks with the same size and the same structure. Thus, we consider the case A₁ = A₂ = ⋯ = A_k = A with A being an s by s correlation matrix, where s = n/k. We further assume that the common block matrix A can be parameterized by a single parameter θ∈(0, 1). We suggest the following two forms for A(θ):

where 1{·} is the indicator function and q is a real constant q≥1. It is important to highlight that A^*(θ) is positive definite for all values of θ∈(0, 1), whereas A′(θ) may not be positive definite for some combinations of the values of q≥1, the parameter θ∈(0, 1) and the matrix size s. Hence, A^*(θ), and consequently Σ(θ) based on A^*(θ), is a correlation matrix for all θ∈(0, 1), while A′(θ) may not be a correlation matrix for some θ∈(0, 1). However, A′(θ) is the equicorrelated matrix in the limit case q = 1, i.e., A′(θ) is positive definite for all θ∈(0, 1) when q = 1. For both structures A^*(θ) and A′(θ), the parameter θ is related to the dependence strength between the features in the block, whereas θ = 0 means that the features are independent. The case θ = 1 for A^*(θ) corresponds to a deterministic linear relation between the features that belong to the same block.

The form A^*(θ) is a Toeplitz type matrix, which is used often in stationary time series, [30, p. 311, Section 8.8.4], while A′(θ) is a linear generalization of the equicorrelated matrix, considered in [23] and [29]. Under the dependence structure A^*(θ), the correlation coefficients between the features decay exponentially, from θ for adjacent features to θ^s−1 for the furthest features, while the correlations in a block modeled by A′(θ) are decreasing in a linear manner from θ to θ/q. The equicorrelated dependence can be included in the form A′(θ) as the limit case q = 1. For simplicity, in this paper we consider only the case q = 2, i.e., correlation changing from θ to θ/2. It should be noted that when q = 2 the matrix A′(θ) is not always positive definite for arbitrary θ∈(0, 1) and block size s, however in our simulation study the values of θ and s are chosen such that A′(θ) is a correlation matrix, unless otherwise noted. In this work we further assume that all tests for the features within a given block are associated with either true null hypotheses or true alternatives. This corresponds to assuming that all genes in a pathway have similar behavior and consequently their expression levels are either significant or not. A general algorithm for estimating the value of θ by the observed data is presented in Section 3, while simulation results under the two forms A^*(θ) and A′(θ) are given in Section 4. Next, we describe the t-test procedures for comparing the two samples X and Y.

For the j-th feature, j = 1, 2…, n, let us consider the Student's t-test for paired samples for testing H₀:μ_j = 0 against H₁:μ_j≠0. Let us denote the average of each column in the matrix of differences D={Xi(j)-Yi(j)}i=1,j=1m,n by

In the context of gene expression levels, this could be the differences between the logarithms of the expression levels, or equivalently, log(fold change). Then, the Student's t-statistic for each feature is

where σ^m(j)=∑i=1m(Di(j)-D¯j)2m-1 is the sample estimation of the standard deviation of the j-th difference. The joint distribution of the test statistics Tm(1),Tm(2),…,Tm(n) has the following asymptotic property, given in the next proposition.

Proposition 1. Under H₀ and the assumptions in Subsection 2.1 with an arbitrary correlation (dependence) matrix Σ, the joint distribution of Tm(1),Tm(2),…,Tm(n) is asymptotically normal such that

Proof. Let B=(bi,j)i,j=1m be an unitary matrix such that b1,i=m-1/2, for i = 1, 2, …, m. Then,

such that

and

with I_m×m being the identity matrix of size m.

Moreover, D¯jm is independent of Vm(j)=(m-1)∑i=2m(Zi(j))2 with

Thus, the t-test statistics are

Finally, by (3) and Slutsky's theorem it follows that

Let us consider the result of Proposition 1. In case the dependence correlation matrix Σ is known, then by (2) and the scaling property of the multivariate normal distribution:

where Z⁽¹⁾, Z⁽²⁾, …, Z⁽ⁿ⁾ are iid standard normal variables. Thus, instead of applying the B-H method to the p-values corresponding to the t-test statistics Tm(1),Tm(2),…,Tm(n), we can apply the B-H procedure to the p-values obtained by the decorrelated statistics

with Sm(j)d→N(0,1) as m → ∞, for j = 1, 2, …, n.

Since we assume that the correlation matrix Σ has a block structure, the statistic Sm(j) in (4) is a linear combination only of the statistics Tm(i) that belong to the same block as Sm(j), for j = 1, 2, …, n. Thus, the order of the blocks in transformation (4) is the same as in the raw t-test statistics. However, the statistic Sm(j) no longer corresponds to the j-th feature (gene) expression, but rather to a combination of features in the same block. Hence, if the B-H procedure detects Sm(j) as significant, this by itself does not necessarily mean that the j-th feature is significant. Because of our assumption that either the null hypothesis H₀ is true for all features in a specific block or H₁ holds for all features within a block, the FDR and power of a multiple testing procedure depend only on the number of rejected hypotheses within a given block instead of on the specific rejected tests. In Section 3.2, we propose a post-hoc procedure for identifying potentially missed significant features by the transformed B-H method.

In practice, the dependence matrix is unknown and needs to be estimated. Under our assumptions in Section 2.1, Σ can be parameterized by θ. Therefore, an estimation of Σ corresponds to an estimation θ^ of the parameter θ. For finding θ^, we propose a procedure based on minimizing a matrix norm of the difference between Σ(θ) and C^ for all values θ∈(0, 1), where C^ is the sample estimation of the correlation matrix Σ. For the matrix norm, based on our simulation results in Section 4.1, we suggest the Frobenius norm ||·||_F or the Max norm ||·||_M, i.e., the value of θ to be estimated by

or by

Alternative matrix norms that are additionally considered in Section 4.1 are the column-sum norm ||·||₁ and the spectral norm ||·||₂, defined for Σ(θ)-C^ as

where s_i(M) is the i-th singular value of a matrix M. Note that ||Σ(θ)-C^||1 is equivalent to the row-sum norm (infinity norm),

since Σ(θ)-C^ is a symmetric matrix. Thus, in Section 4.1 we present the estimation results based only on ||·||₁. More details about the listed matrix norms can be found in [30, Section 3.9].

Our simulation analysis in Section 4.2 shows that the FDR methods based on the p-values obtained from the transformed statistics Sm(1),Sm(2),…,Sm(n) in (4), with Σ=Σ(θ^), lead to a better control of the FDR in comparison to the FDR under the correlated t-statistics Tm(1),Tm(2),…,Tm(n). However, this performance improvement comes at the cost of reducing the power of the applied procedures. Therefore, next we propose an adjustment by combining those methods with the K-S tests for each block in Σ.

Let p1(l), p2(l), …, ps(l) be the raw p-values corresponding to the test statistics of the features in the l-th block of Σ, for l = 1, 2, …, k. Since their distribution is standard uniform under H₀ and the independence assumption, to test their uniformity we consider the one-sided K-S distances

where F^p(l)(x) is the empirical distribution function of p1(l), p2(l), …, ps(l), cf. [32, formula (4.3.3)]. Note that in (6) we consider only the one-sided differences F^p(l)(x)-x, since the test statistics obtained under H₁ are associated with smaller p-values resulting in a peak of F^p(l)(x) around the origin x = 0, i.e., F^p(l)(x)>x for small values of x. Hence, in addition to the rejected hypotheses by the used FDR correction method, we propose to reject all hypotheses corresponding to the blocks with significant K-S statistics Dl+, l = 1, 2, …, k. In such a manner the K-S tests serve as a post-hoc method for identifying potentially missed significant features. However, this might change the FDR and does not guarantee a strict FDR control. Note that the K-S tests can also be used as an ad-hoc tool, i.e., applied independently of the multiple comparison procedure. For the simulations in Section 4.2, we use the K-S tests as a post-hoc exploratory tool rather than a procedure with guaranteed FDR control. Our simulation study indicates that it is best to combine the K-S tests with the transformed B-H procedure based on the transformation (4) with Σ(θ^) and θ^ computed by (5). Thus, in Section 4.2 we present the K-S results only in combination with the transformed B-H algorithm.

Let us first analyze the estimators θ^ of the parameter θ based on different matrix norms, given in Section 3.1. To obtain various structures of the correlation matrix Σ(θ), we fix the sample size m and the number of features n to be 100, while varying the dependence strength θ∈{0.80, 0.90, 0.95} and the number of blocks k∈{5, 10, 20}. As a consequence, the size of each block s = n/k also takes values in the set {5, 10, 20}. Since Σ does not depend on the vector of mean differences, we set μ_j to be 0, for j = 1, 2, …, n. For solving the optimization problem in (5) under different matrix norms we use the function optimize() of the built-in Rstats package. The simulation results for the average bias (Avg. bias) and MSE of the estimators θ^ based on Frobenius ||·||_F, Max ||·||_M, Column-sum ||·||₁ and Spectral ||·||₂ norms are given in Table 1. The values obtained under the block form A^*(θ) are presented on the left-hand side, whereas those under A′(θ) on the right-hand side. The results under A′(θ) for n = 100, k = 20 and θ = 0.95 are missing, since in this case the A′(θ) is not positive definite, i.e., it is not a correlation matrix.

The average bias results suggest that under the Frobenius norm the estimations of θ have uniformly the smallest absolute bias. The estimators θ^ based on ||·||₁ and ||·||₂ tend to overestimate the value of θ, while those derived by ||·||_F and ||·||_M can have a positive or negative bias. However, in all cases the bias under ||·||_F is less than 0.13% and under ||·||_M is less than 5.35%, whereas for some simulation configurations the bias under ||·||₁ and ||·||₂ is up to 18.01% and 32.59%, respectively. In general, it can be concluded that the estimations derived by ||·||_F and ||·||_M have more stable bias performance compared to the ||·||₁ and ||·||₂ norms that have greater bias for larger number of blocks k. Furthermore, the ||·||₁ and ||·||₂ estimations are more sensitive to the dependence form, having smaller bias under A^*(θ) in comparison to A′(θ). Regarding the MSE, the estimations obtained by ||·||_F and ||·||_M clearly outperform those based on ||·||₁ and ||·||₂. Moreover, ||·||_F estimations seem to be more stable in terms of the MSE, indicating robustness in both the bias and variance of the estimator. In our simulations the evaluation of ||·||_F and ||·||_M required less time and resources to the other matrix norms. Therefore, we suggest using the norms ||·||_F and ||·||_M rather than ||·||₁ and ||·||₂ for estimating the value of θ.

Let us now focus on the effect of the sample size m and the number of features n on the performance of θ^ based on ||·||_F and ||·||_M. The mean, standard deviation (SD) and MSE of estimations with ||·||_F and ||·||_M norms, for θ = 0.80, k = 10, m∈{20, 50, 100, 150, 200} and under block structure A^*(θ), are illustrated in Figure 1 (n = 100) and Figure 2 (n = 1, 000). The corresponding results under the block matrix form A′(θ) are presented in Figure 3 (n = 100) and Figure 4 (n = 1, 000). Based on our simulations, when the blocks of Σ(θ) follow the form A′(θ), it seems that n has no effect on the estimation SD and MSE and there is a minor bias reduction of the ||·||_M estimator for larger values of n. In contrast, under the block structure A^*(θ), increasing the size n leads to a notable deterioration of the performance of θ^ derived by ||·||_M, while the mean bias, SD and MSE for the ||·||_F estimator are reduced when n is large. As expected, in all simulation settings the estimation accuracy is improved by increasing the sample size m, since the correlation sample estimation C^ is consistent and approximates Σ(θ) better for larger values of m. In general, the ||·||_F estimator outperforms the one based on ||·||_M in terms of the MSE measure, which summarizes both the estimation bias and variance. Thus, for transforming by (3) the test statistics obtained by the simulations in Section 4.2, we use the estimation θ^ in (5), derived by the Frobenius norm.

To study the effect on the estimation accuracy of θ^ when the block dependence is misspecified, we consider specific types of mismatches, namely estimations θ^ of θ based on the block structure A′(θ), when the true simulated form is A^*(θ), and estimations θ^ derived under A^*(θ), in the case of true underlying dependence simulated in the form A′(θ). The obtained average bias and MSE values under the same simulation setup as in Table 1 are given in Table 2. From the presented results we conclude that both the average bias and MSE of θ^ are increased under misspecification of the dependence form, especially in the case when θ^ is based on A′(θ), while the true simulated block structure is A^*(θ). The average bias of θ^ derived by the column-sum norm ||·||₁ is up to 51.51% (for θ = 0.80 and k = 5). When θ is estimated under A^*(θ), whereas the true underlying dependence is A′(θ), the estimations θ^ are closer to the true value of θ, with Frobenius and Max norms having average bias less than 6.58% when θ≥0.90. Note that for all matrix norms and in both misspecification scenarios, the absolute average bias and MSE decrease for larger values of the correlation coefficient θ and the number of blocks k in Σ(θ). In addition, the estimator θ^ based on the Frobenius norm ||·||_F seems more sensitive under misspecified dependence structure in comparison to the estimator θ^ derived by the Max norm ||·||_M. Similar to the case of correctly specified form, the norms ||·||_F and ||·||_M lead to the smallest MSE of θ^ in the most parameter settings. Nevertheless, the estimations based on ||·||_F outperform the others for large values of θ and k. Therefore, we recommend Frobenius norm to be used for the estimation θ^ of θ.

To compare the multiple testing procedures described in Section 3 under the block dependence structures A^*(θ) and A′(θ), we consider their mean FDR and power, obtained by averaging the associate results in each simulation. In Table 3 we present the average FDR and power of the B-H and B-Y procedures applied to the p-values corresponding to the raw and transformed t-statistics in Sections 2.2 and 3.1, respectively, for fixed α = 0.05, n = 100, m = 100, k = 10 and varying parameter values θ∈{0.80, 0.90, 0.95} and alternative differences μ_j∈{1/2, 1, 2}. The transformed statistics Sm(1),Sm(2),…,Sm(n) in (4) are computed by estimating Σ with Σ(θ^), where θ^ is calculated by (5). Since the asymptotic normality of Sm(j) can be used only for large values of the sample size m, the distribution of Sm(j) under H₀ can deviate from N(0,1) when m < 100. Our experience shows that when m = 20 or m = 50 the distribution of Sm(j) is better approximated by t(m−2) than by N(0,1) in terms of the critical values under H₀, although Sm(j) is a standardized linear transformation of t-statistics with m−2 degrees of freedom, which is not t(m−2) distributed. Since for m≥100 the Student's t distribution is practically identical to the standard normal distribution, in all simulations we have computed the p-values corresponding to Sm(j) by using the t(m−2) distribution. In addition, we present results for the combination of the B-H methods with the K-S tests, given in Section 3.2.

Based on our simulations, the B-Y algorithm is very conservative with average FDR much smaller than the significance level α = 0.05 and the FDR results of the other multiple testing procedures. This is in agreement with the way it is constructed, see [20]. Moreover, in all simulation setups the B-H method has the largest average power compared to the other multiple testing approaches. Thus, under our dependence assumptions, the B-Y is outperformed by the B-H algorithm. However, the average FDR of the B-H procedure is less than α = 0.05 as well. Since the dependence structure of Σ is part of the normalization (4), the transformed B-H and B-Y techniques lead to an increase of the FDR, which is also shown by the results in Table 3. Similar to the transformations applied to the raw p-values, the transformed B-H controls the FDR in a better manner compared to the transformed B-Y. Nevertheless, this improvement in the FDR comes at the expense of average power. More specifically, both methods (transformed B-H and B-Y) suffer from low power in case μ_j∈{1/2, 1}, i.e., when the alternative H₁ is closer to the null hypothesis H₀. In order to fix this disadvantage, we suggest to combine the transformed B-H procedure with the K-S tests, which allows a balance to be achieved between controlling the FDR and preserving power. Further note that the raw B-H and B-Y procedures have almost identical FDR and power results for μ_j = 1 and μ_j = 2, indicating that their performance is stable above a certain threshold. Regarding the influence of the dependence strength, associated with θ, from our simulations it can be concluded that increasing the value of θ, i.e., increasing the correlation between the features in a block, has a negative effect on the FDR control and reduces the average power for all multiple testing methods.

Let us concentrate now on the performance of the multiple testing procedures under various values of the sample size m. Our simulation results for n = 100, k = 10, θ = 0.80, μ_j = 2 and m∈{20, 50, 100, 150, 200} are illustrated in Figures 5, 6 under block structures A^*(θ) and A′(θ), respectively. Along with the average FDR and power, we also provide their standard deviation (SD), estimated by the simulation values. It seems that for all testing procedures the FDR stabilizes when m≥100. The transformed B-H technique shows the best control of the FDR for m≥100, with good approximation of the significance level α = 0.05 and smaller SD than the one based on the raw B-H algorithm. One can observe that the SD of FDR for the B-Y and transformed B-Y methods are very small since the corresponding FDR values are close to 0 and thus have a small room for variation. Since the transformed B-H and B-Y have an unsatisfactory average power when the small sample size m is small, it appears that for m < 100 the transformed schemes should be combined with the K-S tests. However, the raw B-H procedure outperforms the others in terms of average power when m < 100, whereas the transformed B-H method has the best balance of FDR control and power for larger samples (m≥100). Hence, in case the alternative hypothesis H₁ is closer to the null H₀, i.e., μ_j is closer to 0, we recommend to combine the transformed approaches with the K-S tests when the sample size m is not large enough. In addition, we should mention that simulation results, obtained for the number of features n∈{500, 1000} and which are not included in this paper, suggest similar conclusions as those under n = 100.

A comparison between the results shown in Figures 5, 6 suggests that the presented methods have larger average power and FDR that is closer to the desired level α under dependence A^*(θ) than under the block form A′(θ). This is also supported by the results for θ = 0.80 in Table 3 and is due the fact that when θ = 0.80 the matrix A^*(θ) is closer to the identity matrix, which corresponds to independence between all features. However, for larger values of θ the structure A^*(θ) presents stronger correlation than the form A′(θ), leading to a better FDR control and larger power of the multiple testing procedures under A′(θ) compared to those under A^*(θ), as shown by the results for θ∈{0.90, 0.95} in Table 3.