Nonparametric relevanceshifted multiple testing procedures for the analysis of highdimensional multivariate data with small sample sizes
 Cornelia Frömke^{1}Email author,
 Ludwig A Hothorn^{2} and
 Siegfried Kropf^{3}Email author
DOI: 10.1186/14712105954
© Frömke et al; licensee BioMed Central Ltd. 2008
Received: 15 June 2007
Accepted: 27 January 2008
Published: 27 January 2008
Abstract
Background
In many research areas it is necessary to find differences between treatment groups with several variables. For example, studies of microarray data seek to find a significant difference in location parameters from zero or one for ratios thereof for each variable. However, in some studies a significant deviation of the difference in locations from zero (or 1 in terms of the ratio) is biologically meaningless. A relevant difference or ratio is sought in such cases.
Results
This article addresses the use of relevanceshifted tests on ratios for a multivariate parallel twosample group design. Two empirical procedures are proposed which embed the relevanceshifted test on ratios. As both procedures test a hypothesis for each variable, the resulting multiple testing problem has to be considered. Hence, the procedures include a multiplicity correction. Both procedures are extensions of available procedures for point null hypotheses achieving exact control of the familywise error rate. Whereas the shift of the null hypothesis alone would give straightforward solutions, the problems that are the reason for the empirical considerations discussed here arise by the fact that the shift is considered in both directions and the whole parameter space in between these two limits has to be accepted as null hypothesis.
Conclusion
The first algorithm to be discussed uses a permutation algorithm, and is appropriate for designs with a moderately large number of observations. However, many experiments have limited sample sizes. Then the second procedure might be more appropriate, where multiplicity is corrected according to a concept of datadriven order of hypotheses.
Background
Nowadays, highdimensional multivariate data are used in agriculture, biology and medicine. A recent example are microarray data, where two groups, for example normal and diseased tissue, are compared using tens of thousands of genes. The aim is to identify those genes with relevant over or underexpression. Therefore, only twosided tests are considered here. Nevertheless, directional onesided relevanceshifted versions are also available [1]. Distinguishing between formal statistical significance and biological relevance is a frequently discussed issue [2]. One reason is that the commonly used pointzero nullhypothesis H_{0} : μ_{2}  μ_{1} = 0 is often biologically inappropriate, because depending on sample size and variance, biologically irrelevant small differences can be marked as statistically different. Therefore, the relevanceshifted nullhypothesis H_{0} : μ_{2}  μ_{1}  δ = 0 should be used. Hereby, the problem of the choice of δ appears [3]. Instead of absolute relevance margins, the use of relative margins may be more appropriate in some applications. For example, a compound will be declared potentially mutagenic in the Ames mutagenicity assay if the number of revertants is at least double of those in the control; this is the socalled twofold rule [4]. Another example is the characterization of the antineoplastic activity of a new compound by its ratio of the mean tumor volume under treatment to that of the control [5]. For microarray experiments, a specific foldchange may be of interest as well. For example, Guo et al. [6] searched genes which were significant to the unadjusted α of 5% and have a fold change greater than 1.5. To analyze the relevanceshifted hypothesis in gene expression data, Li and Wong [7] propose using a confidence interval to estimate the fold change. However, this confidence interval requires the normality assumption and does not account for multiplicity. Both problems will be addressed in the subsequent proposals.
The problem occurs with the validation of normal assumptions for high dimensional data with small sample sizes. It seems to be hopeless to empirically characterize the distribution as multivariate normal. Hence nonparametric methods may be advantageous. In the current literature, examples for the analysis of multivariate studies using nonparametric methods can be found [8, 9]. The present article focusses on nonparametric approaches as well. Relevanceshifted Wilcoxon rank statistics are used as basic test statistics in both approaches considered in this paper. Parametric test procedures for relevanceshifted hypotheses can be found in Frömke [1].
Since a local decision is provided for each of the thousands of genes, the resulting multiplicity problem has to be considered, too. Otherwise, the probability to falsely reject a null hypothesis increases dramatically. To overcome the problem of multiplicity, several approaches are discussed in the literature. Aside from the classical familywise error rate (FWER) – the probability to reject at least one true null hypothesis – the false discovery rate (FDR) introduced by Benjamini and Hochberg [10] is often used [11], giving the expected proportion of falsely rejected hypotheses among the set of all rejected hypotheses. According to the definitions, the FDR criterion usually delivers more 'significant' genes because – in contrast to the FWER – a small rate of falsely positive results is tolerated. Therefore, the main arguments in favour of the FDR (or against the FWER) are the low detection rate of the FWER procedures for highdimensional data combined with a still sufficient type I error control for screening purposes. Nevertheless, many authors emphasize that the FWER criterion is necessary for confirmatory purposes [12–16]. Also the argument of the larger power of FDR procedures has to be qualified. As pointed out by Dudoit et al. [15] and Speed [16], the present FDR controlling procedures are usually based on independence assumptions between the single test statistics which are not acceptable in gene expression data (particularly, the BenjaminiHochberg procedure) or they are corrected for that problem and are then computer intensive and/or so much reduced in their power that the advantage with respect to the FWER procedures is more or less lost. Here, we focus on two FWER controlling procedures (in the strong sense, i.e., keeping the FWER under any constellation of true and false local hypotheses) and we will demonstrate that they may be well applied also in highdimensional data.
The simplest method to correct for multiplicity in both parametric and nonparametric settings is the αadjustment of Bonferroni. Here the unadjusted pvalues of the individual tests are compared with α/m, where m denotes the dimensionality of the data, that is, the number of observed variables. The modification by Holm [17] uses the threshold α/m only for the comparison with the smallest of the m individual pvalues. The next smallest are compared to α/(m  1), α/(m  2), ..., α/1. If one of the pvalues does not fall below the corresponding threshold, then the statistical procedure will stop and this null hypothesis as well as all succeeding ones will be accepted. However, even for lowdimensional data, the Bonferronitype methods are known to be conservative. Furthermore, the potential improvement using the Holm method is minimal in highdimensional data with only a small portion of variables with effects. One reason for the conservativeness is that the Bonferroni correction does not utilize the correlation structure of the variables. In rank tests with their discrete distributions, we have the additional problem that the procedures usually cannot fully exploit the prespecified error level but have to switch to the next possible pvalue less than or equal to the given threshold. Particularly for very small sample sizes and a high dimension m, it might thus even be impossible to reject the null hypothesis using Bonferroni type methods. Both procedures extended in this paper utilize the covariance structure as they are based on permutations of the whole multivariate observation vectors although the technical procedure does not show that explicitly.
In the following text the two original procedures, which will be extended in this article, are briefly presented. The first procedure is the wellknown permutation algorithm of Westfall and Young [18], as proposed by Dudoit et al. [19] for the analysis of microarray data. Just as the BonferroniHolm correction, this method is a stepdown procedure. However, it consists of a permutation algorithm to compute the null distribution of the pvalues. By permuting the variables, the algorithm takes their dependence structure into account. Given certain data conditions, above all not too small samples, this algorithm is probably the most powerful testing procedure for highdimensional data in the current literature.
 1.
For each of the m variables separately, compute the interquartile range using the combined data of the two samples. Order the m variables for decreasing values of their empirical interquartile ranges. The interquartile ranges serve as selector statistics.
 2.
Calculate the twosample Wilcoxon test at the unadjusted level α in this order as long as significance is attained. Stop at the first nonsignificant result.
Assuming that all variables are measured on an equal scale and have similar variability within group, a large variability over the pooled samples for some of the variables could be a hint for large group differences. Therefore, the interquartile ranges of the pooled samples for each variable are used as selector statistics. The proof for the exact control of the familywise error rate utilizes, roughly speaking, the independence of rank and order statistics. If – contrary to the assumption – the variables have heterogeneous variability, then the procedure looses power. For example, Frömke [1] presents simulation studies, where the standard deviations vary by factor of 1.5 and the procedure looses approximately 10% power. The loss in power increases with increasing variability of standard deviations. However, the procedure controls the type I error in any case. A parametric counterpart of this procedure based on the theory of spherical distributions can be found in Kropf and Läuter [21].
where the operator '/' indicates a componentwise division of vectors. Thus, θ_{ j }denotes the true ratio of the treatment medians of variable j.
For each of the m variables, it shall be tested, whether the twosided null hypothesisH_{0,j}: θ_{ lower }≤ θ_{ j }≤ θ_{ upper }
can be rejected in favor of the alternativeH_{1,j}: θ_{ j }<θ_{ lower }or θ_{ j }> θ_{ upper },
where θ_{ lower }≤ 1 and θ_{ upper }≥ 1 denote the lower and the upper relevance threshold.
withδ= (δ_{1}, ..., δ_{ m })' and δ_{ j }= ln(θ_{ j }) (j = 1, ..., m)
where δ_{ lower }= ln(θ_{ lower }) ≤ 0 and δ_{ upper }= ln(θ_{ upper }) ≥ 0 denote the lower and the upper relevance threshold. In practice, the choice of δ_{ lower }and δ_{ upper }is dependent on the experimental question. For example, in microarray experiments the thresholds can be set to δ_{ lower }= δ_{ upper }= 0.4055, 0.6931 or 0.7885. This is equivalent to testing for a foldchange in gene expression of 1.5, 2 or 2.2 [6, 22, 23]. So an obvious approach would be to use the above mentioned two procedures after the logarithmic transformation and an additional shift by the relevance thresholds. However, this is associated with some problems. The shifted onesided tests control the familywise error rate on the threshold which was used for shifting. But here we have two onesided tests and two relevance thresholds, the lower and the upper one, and it is necessary to find some combination rule. It is likely that the second threshold (which is not used for the shift at that moment) is far enough from the first one so that a type I error caused by the onesided test at the opposite border of the null space is unlikely. But the whole parameter space between both thresholds belongs to the null hypothesis as well and there is no proof that the two basic procedures control the type I error or are conservative under these conditions though one would expect a monotone behaviour of the rejection probability for increasing deviations from the exact null point. Finally, a correction is necessary for the selector statistic in the second procedure with data dependent sequential testing. Otherwise, variables with no shift or a only small one (but within the tolerance region) could have larger expected values for the selector than variables under the alternative hypothesis. The procedure would then loose its power by stopping prematurely.
The modifications of the exact procedures for point null hypotheses described in the following section have been adapted to these problems in an empirical manner. No exact proof exists for the control of the familywise error rate. Therefore, results of simulation experiments are presented after the detailed description of the modified procedures and their demonstration in examples. An R package for the methods is available [24].
Results and Discussion
Algorithm
Relevanceshifted permutation algorithm
We will first introduce a relevanceshifted modification of the permutation algorithm for stepdown minP adjusted pvalues of Westfall and Young [18] for point null hypotheses. More strictly speaking, we are starting from a proposal from Ge et al. [25] which delivers the same results as that of Westfall and Young but is less time consuming. Whereas the original algorithm requires two permutation runs, one for the calculation of raw pvalues and a second one for multiplicity adjusting, Ge et al. [25] share the permutations of both runs.
In order to detect the deviation from the null hypothesis at both relevance thresholds, two passes are needed for each variable, one for relevant decrease and another one for relevant increase. Finally, out of the two onesided pvalues, a twosided one is computed for each variable. The passes themselves consist of two parts. The relevanceshifted permuted unadjusted pvalues from Wilcoxon's rank sum test are computed first. Then the unadjusted pvalues are corrected for multiplicity. As mentioned above, we will use the log transformed observations and relevance thresholds.
The proposed algorithm is given here in detail for the test on decrease:
Part 1: Permutation algorithm for raw pvalues

Fix the thresholds δ_{ lower }= ln(θ_{ lower }) and δ_{ upper }= ln(θ_{ upper }).

Create the pseudosample vectors ${y}_{ik}^{\ast}=\left({y}_{ik1}^{\ast},\mathrm{...},{y}_{ikm}^{\ast}\right)$' with ${y}_{1kj}^{\ast}$ = y_{1kj}+ δ_{ lower }and ${y}_{2kj}^{\ast}$ = y_{2kj}(i = 1, 2; k = 1, ..., n_{ i }; j = 1, ..., m).

In the b^{ th }permutation step, b = 0, ..., B (b = 0 corresponds to unpermuted data) do:

For each variable, compute the onesided Wilcoxon rank sum statistic W_{1b}, ..., W_{ mb }for the pseudosamples:${W}_{jb}={\displaystyle \sum _{k=1}^{{n}_{2}}{r}_{2jkb}}$

Permute the n_{1} + n_{2} = N pseudosample vectors ${y}_{ik}^{\ast}$ (i = 1, 2; k = 1, ..., n_{ i }).

Calculate the onesided raw pvalues for hypothesis H_{0,j}: δ_{ j }≥ δ_{ lower }as$\begin{array}{c}{p}_{j,lower}^{\ast}=\frac{\#\{b:b>0\text{and}{W}_{jb}\le {W}_{j0}\}}{B}\\ \begin{array}{cc}\text{for}& j=1,\mathrm{...},m\end{array}.\end{array}$(10)
Part 2: Permutation algorithm for stepdown minP adjusted pvalues

Renumber the m variables such that ${p}_{1,lower}^{\ast}\le \mathrm{...}\le {p}_{m,lower}^{\ast}.$

Prepare three matrices for further computation:
Two empty matrices P = (p_{ jb }) of size m × B and Q = (q_{ jb }) of size (m + 1) × B are filled successively from the bottom to the top in the course of the following algorithm.

Set q_{m+1,b}= 1 for b = 1, ..., B.

For j = m, m  1, ..., 1 do:
which is in row j of matrix W for each W_{ jb }the proportion of test statistics W_{ jb' }equal to or smaller than W_{ jb }.

Determine the j th row of matrix Q as the successive minimaq_{ jb }= min(q_{j+1,b}, p_{ jb }), b = 1, ..., B.

Enforce monotonicity of ${\stackrel{~}{p}}_{j,lower}^{\ast}$:$\begin{array}{l}{\stackrel{~}{p}}_{1,lower}^{\ast}:={\stackrel{~}{p}}_{1,lower}^{\ast},\hfill \\ {\stackrel{~}{p}}_{j,lower}^{\ast}:=\mathrm{max}({\stackrel{~}{p}}_{j1,lower}^{\ast},{\stackrel{~}{p}}_{j,lower}^{\ast})\hfill \\ \begin{array}{cc}\text{for}& j=2,\mathrm{...},m.\end{array}\hfill \end{array}$(15)

Revoke the renumbering of the variables in the beginning of Part 2.
For the test on increase, repeat the entire procedure with the pseudosample vectors ${y}_{ik}^{\ast}=\left({y}_{ik1}^{\ast},\mathrm{...},{y}_{ikm}^{\ast}\right)$', where ${y}_{1kj}^{\ast}$ = y_{1kj}+ δ_{ upper }and ${y}_{2kj}^{\ast}$ = y_{2kj}(j = 1, ..., m) and the rank sum test on increase to achieve the onesided multiplicity adjusted pvalues on increase ${\stackrel{~}{p}}_{j,upper}^{\ast}$. Finally, the twosided adjusted pvalues are given by ${\stackrel{~}{p}}_{j}^{\ast}=\mathrm{min}(2\cdot {\stackrel{~}{p}}_{j,lower}^{\ast},2\cdot {\stackrel{~}{p}}_{j,upper}^{\ast})$.
Procedure with a datadriven order of relevanceshifted hypotheses
An empirical extension for the nonparametric procedure of Kropf et al. [20] for relevanceshifted hypotheses will now be proposed:

Select the two relevance thresholds δ_{ lower }= ln(θ_{ lower }) and δ_{ upper }= ln(θ_{ upper }).

Determine the pseudosample vectors ${y}_{ik}^{\ast}=\left({y}_{ik1}^{\ast},\mathrm{...},{y}_{ikm}^{\ast}\right)$ with ${y}_{1kj}^{\ast}$ = y_{1kj}+ δ_{ lower }and ${y}_{2kj}^{\ast}$ = y_{2kj}(i = 1, 2; k = 1, ..., n_{ i }; j = 1, ..., m) and calculate for each variable the onesided Wilcoxon rank sum statistic W_{j,lower}= $\sum}_{k=1}^{{n}_{2}}{r}_{2jk$ (j = 1, ..., m), again using the ranks determined over the two combined pseudosamples,

Replace δ_{ lower }by δ_{ upper }and repeat exactly the former step to compute W_{j,upper}(j = 1, ..., m).

Use the permutation algorithm described in the previous procedure or suitable tables to derive the corresponding unadjusted onesided pvalues p_{j,lower}and p_{j,upper}, respectively, for the variablewise Wilcoxon statistics.

Compute the unadjusted twosided pvalues p_{ j }as p_{ j }= min(2·p_{j,lower}, 2·p_{j,upper}) (j = 1, ..., m) for each variable.

In order to prepare the determination of selector statistics, calculate the sample medians for the j th (logarithmic but not shifted) variable in sample 1 and 2, ${\tilde{y}}_{1j}$ and ${\tilde{y}}_{2j}$, respectively, and, once again, derive pseudosample values by$\begin{array}{l}{y}_{1kj}^{\ast \ast}=\{\begin{array}{lll}{y}_{1kj}+{\delta}_{lower}\hfill & \text{if}\hfill & {\stackrel{~}{y}}_{2j}{\stackrel{~}{y}}_{1j}<0\hfill \\ {y}_{1kj}+{\delta}_{upper}\hfill & \text{if}\hfill & {\stackrel{~}{y}}_{2j}{\stackrel{~}{y}}_{1j}\ge 0\hfill \end{array}\hfill \\ {y}_{2kj}^{\ast \ast}=\{\begin{array}{lll}{y}_{2kj}{\stackrel{~}{y}}_{2j}+{\stackrel{~}{y}}_{1j}+{\delta}_{upper}\hfill & \text{if}\hfill & 0\le {\stackrel{~}{y}}_{2j}{\stackrel{~}{y}}_{1j}<{\delta}_{upper}\hfill \\ {y}_{2kj}{\stackrel{~}{y}}_{2j}+{\stackrel{~}{y}}_{1j}+{\delta}_{lower}\hfill & \text{if}\hfill & {\delta}_{lower}<{\stackrel{~}{y}}_{2j}{\stackrel{~}{y}}_{1j}<0\hfill \\ {y}_{2kj}\hfill & \text{else},\hfill \end{array}\hfill \end{array}$
(k = 1, ..., n_{1} or k = 1, ..., n_{2}, respectively; j = 1, ..., m).

Compute a selector statistic for each variable as the interquartile range IRQ_{ j }(difference of percentiles 75 and 25) from the combined sample values ${y}_{ikj}^{\ast \ast}$, (i = 1, 2; k = 1, ..., m).

Sort the m pvalues p_{ j }for decreasing values of the corresponding selector statistics IQR_{ j }.

In this order, compare the corresponding pvalue with the unadjusted α for each variable j. The original variable has a significantly relevant ratio of medians if p_{ j }<α and all previously tested null hypotheses have been rejected, too.

Stop at the first nonsignificance and accept the null hypothesis for all further variables.
The different formulae for the selector statistic depending on the difference of the two group medians (positive or negative, within or without the tolerance region) ensure an appropriate sorting of the variables giving the procedure a high power.
In the following sections, the BonferroniHolm procedure will be used for comparison. The unadjusted pvalues will also be taken from twosample Wilcoxon tests with the pseudosample values as in the above two methods. The onesided pvalues p_{j,lower}and p_{j,upper}will be determined separately for each of the m variables, each with the corresponding shift in the pseudosamples. These unadjusted pvalues can be either taken from the first pass of the minP algorithm or from the second procedure. Then – as above – twosided pvalues will be calculated using p_{ j }= min(2·p_{j,lower}, 2·p_{j,upper}) (j = 1, ..., m) and will be used as the basis for the BonferroniHolm adjustment.
Testing
Performance on simulated data
To confirm the control of the FWER, extended simulation studies were applied to the new procedures. A small part of the results for twosided testing is shown in the following two tables. All scenarios were tested with three levels of relevance thresholds. For comparison with Kropf et al. [20], in one type of setting the thresholds were set to θ_{ lower }= θ_{ upper }= 1 (δ_{ lower }= δ_{ upper }= 0). In this case, the procedure with a datadriven order of relevanceshifted hypotheses reduces to the exact nonparametric procedure with a datadriven order of pointzero hypotheses for two unpaired samples applied to the logarithmized data. In the remaining two types of settings, the thresholds were set to ${\theta}_{lower}^{1}$ = θ_{ upper }= 1.5
(δ_{ lower }= δ_{ upper }= 0.4055) or to the extreme case of ${\theta}_{lower}^{1}$ = θ_{ upper }= 5 (δ_{ lower }= δ_{ upper }= 1.6094). Fifty observed variables were tested in all simulated situations.
To test if the procedures control the FWER in the weak sense, which is in the special case where all null hypotheses are true and the simulated FWER is less or equal to the selected α, 25 variables were generated with expected values of μ_{1j}= μ_{2j}/θ_{ upper }= 100 and true standard deviations of σ_{1j}= σ_{2j}/θ_{ upper }. The remaining 25 had μ_{1j}/θ_{ upper }= μ_{2j}= 100 and σ_{1j}/θ_{ upper }= σ_{2j}.
Furthermore, the control in the strong sense is important. In this case, the FWER is protected if some null hypotheses may be true and others false but the probability to reject any true null hypothesis is less or equal to α. For the assessment of the control, 45 variables were simulated under the null hypothesis and had the same true mean values as for the weak control; 22 were set to a nonrelevant decrease and 23 to an increase. From the remaining five variables, two had a relevant ratio of treatment means with μ_{1j}= 100·θ_{ upper }+ 50 and μ_{2j}= 100 with σ_{1j}= 10·θ_{ upper }+ 5, σ_{2j}= 10 and the other three had μ_{1j}= 100 and μ_{2j}= 100·θ_{ upper }+ 50 with σ_{1j}= 10, σ_{2j}= 10·θ_{ upper }+ 5.
All variables had equal pairwise correlations ρ and equal variances 'on a logarithmic scale'. Together with the sample size per group, these parameters differed between the individual simulation settings and are noted in the table. If not stated otherwise, the random numbers were generated from a normal distribution, the nominal FWER was set α = 5% and the empirical FWER was computed with 10,000 simulation runs each. The modified WestfallYoung permutation algorithm is shown as 'permutation' in the following tables and figures and the procedure with a datadriven order of hypotheses is shown as 'selector'.
Simulation results of the FWER for normal distributed data. This table shows simulation results of the relevanceshifted WestfallYoung permutation algorithm using 500,000 permutation runs ('permutation') and the procedure with a datadriven order of relevanceshifted hypotheses ('selector') for different levels of sample sizes, variances and correlation coefficients using normal distributed data.
selector  permutation  

${\theta}_{lower}^{1}$ = θ_{ upper }  n _{ i }  σ _{ ij }  ρ  weak  strong  weak  strong 
1  5  10  0.1  3.42  3.05  0  0 
1.5  5  10  0.1  2.53  1.98  0  0 
5  5  10  0.1  2.53  1.87  0  0 
1  10  10  0.1  4.39  4.22  4.65  4.12 
1.5  10  10  0.1  2.70  2.81  3.30  3.14 
5  10  10  0.1  2.70  2.48  3.30  3.01 
1  10  10  0.5  4.26  4.22  4.74  4.61 
1.5  10  10  0.5  3.08  3.07  3.36  3.71 
5  10  10  0.5  3.08  2.89  3.36  3.59 
1  5  10  0.9  3.12  3.13  1.42  0.91 
1.5  5  10  0.9  2.83  2.61  1.37  1.45 
5  5  10  0.9  2.83  2.80  1.37  1.44 
1  5  15  0.1  3.40  2.52  0  0 
1.5  5  15  0.1  2.48  1.73  0  0 
5  5  15  0.1  2.48  1.98  0  0 
1  10  15  0.1  4.27  4.21  4.65  4.10 
1.5  10  15  0.1  2.67  2.78  3.30  3.11 
5  10  15  0.1  2.67  2.51  3.30  2.97 
1  20  10  0.9  5.01  4.96  4.90  4.58 
1.5  20  10  0.9  4.04  3.97  3.29  3.64 
5  20  10  0.9  3.90  3.97  3.29  3.64 
Simulation results of the FWER for nonnormal distributed data. This table shows the results for similar settings as in Table 2. Again, different levels of sample sizes, variances and correlation coefficients were tested, however nonnormal distributed data was generated.
selector  permutation  

${\theta}_{lower}^{1}$ = θ_{ upper }  n _{ i }  σ  ρ  weak  strong  weak  strong 
1  5  10  0.1  3.31  2.87  0  0 
1.5  5  10  0.1  2.40  1.57  0  0 
5  5  10  0.1  2.40  1.78  0  0 
1  10  10  0.1  4.15  4.32  4.88  4.18 
1.5  10  10  0.1  2.48  2.64  3.46  3.42 
5  10  10  0.1  2.48  2.37  3.46  3.15 
1  5  10  0.9  3.19  2.96  0.80  0.39 
1.5  5  10  0.9  2.78  2.45  0.83  1.12 
5  5  10  0.9  2.78  2.64  0.84  1.07 
1  10  10  0.9  4.33  4.48  4.65  4.25 
1.5  10  10  0.9  3.51  3.59  3.26  0.33 
5  10  10  0.9  3.51  3.49  3.26  0.33 
1  5  5  0.9  3.16  3.13  0.80  0.46 
1.5  5  5  0.9  2.77  2.8  0.84  1.12 
5  5  5  0.9  2.77  2.44  0.84  1.12 
1  5  15  0.9  3.17  2.73  0.80  0.31 
1.5  5  15  0.9  2.76  2.44  0.80  1.09 
5  5  15  0.9  2.76  2.56  0.84  1.03 
The results in the tables show that the new procedures control the FWER empirically. Likewise, the FWER is protected for twosided testing in further simulated situations, including other settings of the true mean values, skewed data, variances and correlations among the variables.
Extended results for onesided testing using the procedure with a datadriven order of relevanceshifted hypotheses are also given [1]. Small increases of the FWER occurred in that case. The largest increase for the nominal α level of 5% was 6.3%. Error rates for the permutation algorithm corresponding to the onesided case have not yet been been analyzed.
The control of the FWER is a premise of a statistical test. However, the aim of the experiments discussed here is to find variables which discriminate two kinds of treatments with a high probability. Hence, graphical representations of the simulation results in terms of the power of the new procedures compared to a standard technique will now be shown.
The simulation setting is nearly the same as for the control of the FWER in the strong sense. However, the setting of the expected values of variables under H_{0} was changed. For the control of the FWER, these variables had a ratio of means set to one of the margins of the null hypothesis because this choice resulted in the largest empirical FWER compared to variables with ratios closer to 1. A more realistic setting was selected for the simulation of the power, where a variable under H_{0} received a random ratio of means. This random value was a number θ_{ lower }≤ τ ≤ θ_{ upper }. Two sets of values were created to generate τ. One set included all values between 1 and θ_{ upper }in steps of 0.05. To receive an equal amount of ratios in the second set, all values between 1 and ${\theta}_{lower}^{1}$ in steps of 0.05 were computed and the second set took the inverse of these values. The sets were combined and τ was chosen separately for each variable. If τ < 1 then expectation values were set to μ_{1j}= 100/τ and μ_{2j}= 100 and the standard deviations were set to σ_{1j}= σ_{2j}/τ. Otherwise the true mean values were μ_{1j}= 100 and μ_{2j}= τ·100 with σ_{1j}= σ_{2j}/τ.
As for the simulations of the FWER, α = 5% and each result consisted of 10,000 simulation runs. In the tested scenarios, the thresholds were set to ${\theta}_{lower}^{1}$ = θ_{ upper }= 2 (δ_{ lower }= δ_{ upper }= 0.6931) and σ_{ ij }= 10.
All further settings of the parameters are given in the captions of the figures. The figures show the ratio of detected false null hypotheses as an estimation of the proportional power, which is defined as the average probability of rejecting the false null hypotheses [19]. The power of the exact relevanceshifted Wilcoxon rank sum test on ratio with the multiplicity correction of BonferroniHolm ('BonferroniHolm') is plotted together with simulation results of the two new procedures.
As in Figure 3, the procedure with a datadriven order of hypotheses is more powerful than the permutation algorithm if the sample sizes are small. However, using a larger sample size the permutation algorithm is preferable. The BonferroniHolm correction achieves no power, because the procedure is too discrete. If an experiment consists of 5,000 variables, a sample size of 12 per group is required to achieve a twosided pvalue of 3.7%. For example, using a sample size of 11 per group, the smallest achievable twosided pvalue is 14.2%. Irrespective of the effect size, this pvalue cannot be less than α = 5%.
The choice of the procedure with the best power does not only depend on the sample size. In particular with an increasing α, the permutation algorithm and the BonferroniHolm correction are more powerful than the procedure with a datadriven order of hypotheses with sample sizes as low as 7 or 10. The choice is also dependent on the correlation among the variables as shown in Figure 1. Furthermore, the power of the procedure with a datadriven order of hypotheses reduces drastically in the case of variance heterogeneity among the variables. To be powerful, the procedure requires approximately homogeneous variances after the logarithmic transformation. Corresponding simulation results to these influences are given in Frömke [1]. Although the impact of a varying number of variables was not examined, it can be assumed to have significant effects as well.
Implementation
Method comparison using a publicly available dataset
This section illustrates the application of the two procedures using a subset of the microarray study published by Khan et al. [27]. The entire data set consists of four subgroups of small, round blue cell tumors (SRBCTs) of childhood. Cell lines are available for all four subgroups, and biopsy material is available for two subgroups. The subset of the original study used here incorporates the biopsy material, which consists of 13 samples of the Ewing family of tumors (EWS) and 10 samples of the rhabdomyosarcoma (RMS). Furthermore, all 2,308 genes of the original data set will be analyzed. For the following analysis, significant twofold under or over expressions to an α = 5% are sought. Hence, the thresholds are set to ${\theta}_{lower}^{1}$ = θ_{ upper }= 2 corresponding to δ_{ lower }= δ_{ upper }= 0.6931.
Results from the publicly available dataset. Results from the publicly available dataset: This table shows the results of the relevanceshifted WestfallYoung permutation algorithm using 500,000 permutation runs, the procedure with a datadriven order of relevanceshifted hypotheses and the BonferroniHolm correction. Note that the fourth and the sixth column are necessary for the the procedure with a datadriven order of relevanceshifted hypotheses only.
procedure  image id.  ratio  selector statistic  pvalue  test decision  ranking 

permutation  770394  0.051    0.00278  (reject H_{0})  6 
814260  0.032    0.00278  (reject H_{0})  75  
244618  24.918    0.00348  (reject H_{0})  7  
207274  4259.257    0.01983  (reject H_{0})  2  
43733  0.040    0.04832  (reject H_{0})  9  
selector  207274  4259.257  6.881  0.00003  reject H_{0}  2 
122159  169.102  6.204  0.02136  reject H_{0}  40  
296448  1445.051  6.041  0.00145  reject H_{0}  1  
34849  0.728  5.239  1.00000  accept H_{0}, stop procedure    
BonferroniHolm  770394  0.051    0.00403  (reject H_{0})  6 
244618  24.918    0.00403  (reject H_{0})  7  
814260  0.032    0.00403  (reject H_{0})  75 
On top the table shows the results for the five significant genes found by the relevanceshifted WestfallYoung permutation algorithm after 500,000 permutation runs. In contrast, the procedure with a datadriven order of relevanceshifted hypotheses detects three significant genes, where one of them was also found using the above method. Three genes are also found with the BonferroniHolm adjustment. They are completely different genes compared to the former procedure, but they were also found using the modified WestfallYoung method.
In this analysis, the permutation algorithm detects more significant variables than both the procedure with a datadriven order of relevanceshifted hypotheses and the αadjustment of BonferroniHolm. As shown in the former section, this can be explained with the general performance of these three methods for the present case of moderately large sample sizes in both groups. However, the procedure with a datadriven order of sequential testing is the only one that found the gene 296448, which according to Khan et al. [27] is the most important one.
Conclusion
The comparison of two groups of individuals with many variables is a common problem in biological studies. In the current literature, procedures are proposed which perform local tests for each variable and correct for multiplicity. Most of these procedures test the pointzero or pointone null hypotheses of a difference or ratio in treatment effects of 0 or 1, respectively. A parametric procedure is available for relevanceshifted hypotheses [7]. In this article, two nonparametric procedures are proposed which perform a local relevanceshifted test on ratio on the two samples for each variable and include a multiplicity correction as well. They are extensions of the WestfallYoung permutation algorithm [18] and of a sequential procedure with datadriven order of hypotheses [20], which consider pointnull hypotheses in their original form.
Both new procedures utilize the correlation structure. In the proofs of the original versions, this can be seen in the fact that they consider permutations of the whole observation vectors and not separate permutations for single variables. In the technical procedures, the influence of the correlation among the variables is not seen explicitly because univariate test statistics and selector functions are calculated. But it is present in the ordering of variables, which is part of both procedures in some way. When the variables are highly correlated, then the relation of their Wilcoxon test statistics or interquartile ranges effectively reflects the relation in the degree of violation of the corresponding null hypothesis. The less these correlations are, the more this relation is overlaid with random influences.
As not all modifications, applied to the pointnull versions, could be covered by the theoretical considerations, simulation experiments were carried out for the control of the FWER and for the assessment of the power. In these experiments, the FWER was always controlled for the twosided test versions discussed in this paper.
The power of the two new proposals and of the BonferroniHolm method was similar to the original procedures for pointnull hypotheses (cf. Kropf et al. [20]). The procedure with datadriven sequential hypothesis testing uses a nonparametric measure of variability in the pooled samples as an additional source of information. This provides an advantage in small samples if the variances of the different variables are more or less homogeneous in the data. This advantage is lost and even reversed with increasing sample sizes. As discussed in Kropf et al. [20], this is due to the fact, that the probability to detect a difference in the unadjusted tests (which is the major input in the other test procedures) increases faster than the probability of the correct ordering of variables with and without deviations from the null hypothesis. Therefore, this ordering becomes the critical part in the sequential procedure for at least moderately large samples. However, data from microarray and proteomics experiments are commonly characterized by a very large number of variables and small sample sizes. The analysis of such experiments using standard multivariate approaches is inappropriate. The proposed procedures can be used instead, particularly if relevance shifted hypotheses are of interest.
Declarations
Acknowledgements
The authors are very grateful to Dr. Frank Bretz and to the anonymous reviewers for their helpful comments and constructive suggestions.
Authors’ Affiliations
References
 Frömke C: Relevanceshifted tests for high dimensional data with small sample sizes. PhD thesis. University of Hannover, Institute of Biostatistics; 2006. [http://www.biostat.unihannover.de/research/thesis/]Google Scholar
 Hauschke D, Schall R, Luus HG: Statistical Significance. In Encyclopedia of Biopharmaceutical Statistics. 1st edition. Edited by: Chow SC. New York: Marcel Dekker; 2000:493–497.Google Scholar
 Lange S, Freitag G: Choice of delta: requirements and reality – Results of a systematic review. Biom J 2005, 47: 12–27. 10.1002/bimj.200410085View ArticlePubMedGoogle Scholar
 Cariello NF, Piegorsch WW: The Ames test: The twofold rule revisited. Mutat Res 1996, 369: 23–31. 10.1016/S01651218(96)900440View ArticlePubMedGoogle Scholar
 Hothorn LA: Statistical analysis of in vivo anti cancer experiments: Tumor growth inhibition. Drug Inf J 2006, 40: 229–238.Google Scholar
 Guo L, Fang H, Collins J, Fan X, Dial S, Wong A, Mehta K, Blann E, Shi L, W T, Dragan YP: Differential gene expression in mouse primary hepatocytes exposed to the peroxisome proliferatoractivated receptor a agonists. BMC Bioinformatics 2006, 7(Suppl 2):S18. 10.1186/147121057S2S18PubMed CentralView ArticlePubMedGoogle Scholar
 Li C, Wong WH: Modelbased analysis of oligonucleotide arrays: model validation, design issue and standard error application. Genome Biol 2001, 2: research0032.1–0032.11. 10.1186/gb200128research0032Google Scholar
 Schaid DJ, McDonnell SK, Hebbring SJ, Cunningham JM, Thibodeau SN: Nonparametric tests of association of multiple genes with human disease. Am J Hum Genet 2005, 76: 780–793. 10.1086/429838PubMed CentralView ArticlePubMedGoogle Scholar
 Troyanskaya OG, Garber ME, Brown PO, Botstein D, Altman RB: Nonparametric methods for identifying differentially expressed genes in microarray data. Bioinformatics 2002, 18: 1454–1461. 10.1093/bioinformatics/18.11.1454View ArticlePubMedGoogle Scholar
 Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Statist Soc B 1995, 57: 289–300.Google Scholar
 Polacek DC, Passerini AG, Shi C, Francesco NM, Manduchi E, Grant GR, Powell S, Bischof H, Winkler H, Stoeckert CJ, Davies PF: Fidelity and enhanced sensitivity of differential transcription profiles following linear amplification of nanogram amounts of endothelial mRNA. Physiol Genomics 2003, 13: 147–156.View ArticlePubMedGoogle Scholar
 Westfall PH, Kropf S, Finos L: Weighted FWEcontrolling methods in highdimensional situations. In Recent Developments in Multiple Comparison Procedures. Volume 47. Edited by: Benjamini Y, Bretz F, Sarkar S. Institute of Mathematical Statistics Lecture NotesMonograph Series; 2004:143–154.View ArticleGoogle Scholar
 Chich JF, David O, Villers F, Schaeffer B, Lutomski D, Huet S: Statistics for proteomics: Experimental design and 2DE differential analysis. Journal of Chromatography B 2007, 849: 261–272. 10.1016/j.jchromb.2006.09.033View ArticleGoogle Scholar
 Witt E, McClure J: Statistics for Microarrays. Design, Analysis and Inference. Chichester: John Wiley & Sons; 2004.Google Scholar
 Dudoit S, H YY, Callow MJ, Speed TP: Statistical methods for identifying differentially expressed genes in replicated cDNA microarray experiments. 2000.Google Scholar
 Speed T: Statistical analysis of gene expression microarray data. Boca Raton: Chapman & Hall/CRC; 2003.View ArticleGoogle Scholar
 Holm S: A simple sequentially rejective multiple test procedure. Scand J Statist 1979, 6: 65–70.Google Scholar
 Westfall PH, Young SS: Resamplingbased multiple testing: Examples and methods for pvalue adjustment. New York: John Wiley & Sons; 1993.Google Scholar
 Dudoit S, Shaffer JP, Boldrick JC: Multiple hypothesis testing in microarray experiments. Statist Sci 2003, 18: 71–103. 10.1214/ss/1056397487View ArticleGoogle Scholar
 Kropf S, Läuter J, Eszlinger M, Krohn K, Paschke R: Nonparametric multiple test procedures with datadriven order of hypotheses and with weighted hypotheses. J Statist Plann Inference 2004, 125: 31–48. 10.1016/j.jspi.2003.07.021View ArticleGoogle Scholar
 Kropf S, Läuter J: Multiple tests for different sets of variables using a datadriven ordering of hypotheses, with an application to gene expression data. Biom J 2002, 44: 789–800. Publisher Full Text 10.1002/15214036(200210)44:7<789::AIDBIMJ789>3.0.CO;2#View ArticleGoogle Scholar
 Zimmermann N, King NE, Laporte J, Yang M, Mishra A, Pope SM, Muntel EE, Witte DP, Pegg AA, Foster PS, Hamid Q, Rothenberg M: Dissection of experimental asthma with DNA microarray analysis identifies arginase in asthma pathogenesis. The Journal of Clinical Investigations 2003, 111: 1863–1874. 10.1172/JCI200317912View ArticleGoogle Scholar
 Iyer VR, Eisen MB, Ross DT, Schuler G, Moore T, Lee JCF, Trent JM, Staudt LM, Hudson J, Boguski MS, Lashkari D, Shalon D, Botstein D, Brown PO: The transcriptional program in the response of human fibroblasts to serum. Science 1999, 283: 83–87. 10.1126/science.283.5398.83View ArticlePubMedGoogle Scholar
 Frömke C.: See R program npmrest.[http://www.biostat.unihannover.de/software/index_en.html]
 Ge Y, Dudoit S, Speed TP: Resamplingbased multiple testing for microarray data hypothesis. Test 2003, 12: 1–44. 10.1007/BF02595811View ArticleGoogle Scholar
 Fleishman AI: A method for simulating nonnormal distributions. Psychometrika 1978, 43: 521–532. 10.1007/BF02293811View ArticleGoogle Scholar
 Khan J, Wei JS, Ringnér M, Saal LH, Ladanyi M, Westermann F, Berthold F, Schwab M, Antonescu CR, Peterson C, Meltzer PS: Classification and diagnostic prediction of cancers using gene expression profiling and artificial neural networks. Nat Med 2001, 7: 673–679. 10.1038/89044PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.