Robust Behrens–Fisher Statistic Based on Trimmed Means and Its Usefulness in Analyzing High-Throughput Data

Abstract

In the context of high-throughput data, the differences in continuous markers between two groups are usually assessed by ordering the p-values obtained from the two-sample pooled t-test or Wilcoxon–Mann–Whitney test and choosing a stringent cutoff such as 10^–8 to control the family-wise error rate or false discovery rate All markers with p-values below the cutoff are declared to be significantly associated with the phenotype. This inherently assumes that the test procedure provides valid type I error estimates in extreme tails of the null distribution. The aforementioned tests assume homoscedasticity in the two groups, and the t-test further assumes underlying distributions to be normally distributed. Cao et al. (Biometrika, 2013, 100, 495–502) have shown that in the context of multiple hypotheses testing the approach based on may not be valid under non-normality and/or heteroscedasticity. Therefore, having a test statistic that is robust to these violations is needed. In this study, we propose a robust analog of Behrens–Fisher statistic based on trimmed means, conduct an extensive simulation study to compare its performance with other competing approaches, and demonstrate its usefulness by applying it to DNA methylation data used by Teschendorff et al. (Genome Res., 2010, 20, 440–446). An R program to implement the proposed method is provided in the Supplementary Material.

1 Introduction

In the context of genetic analysis, it is quite common to compare hundreds of thousands of genetic features such as gene expressions or DNA methylations between cases and controls. The well-known two-sample t-test or its robust analog Wilcoxon rank sum test has been commonly utilized to obtain p-values for comparing genetic features between the two groups at each locus. These p-values are then ordered and chosen to control the family-wise error rate or false discovery rate based on a cutoff; all p-values below that threshold are declared to have significantly different genetic features between two groups, for example, see Hochberg and Tamhane (1987), Storey (2002), and Benjamini and Yekutieli (2007). The validity of this approach is based on two underlying assumptions: 1) the p-values under the null hypothesis would be uniformly distributed, whereas the p-values under alternative hypothesis would tend to have values closer to 0; 2) the null distribution of the test statistic is well controlled even for stringent levels of However, Robins et al., (2000) have shown that the p-values will be uniformly distributed under (null model) only when the test statistic generating the p-values consists of a single distribution. On the other hand, if the null distribution of depends on the nuisance parameters, then the p-values will not be uniformly distributed.

When a study involves testing the differences between several thousands of genes between cases and controls, then it may be reasonable to assume that the sample size would be fixed for all comparisons. However, the p-values for comparing a gene expression profile between the two groups would still be affected by the effect size (standardized difference in the mean expression levels), underlying distributional assumptions (usually normality), and inequality of the variances for the two groups. In the context of a multiple-hypotheses setting, it is clear that where denotes the cutoff based on the test statistic, must be an increasing function in , and Cao et al. (2013) have shown that this monotonicity assumption could be violated and could lead to misleading results when the underlying distributions are not normal and/or have significantly different variances.

In the parametric setting, it is well known that the pooled two-sample t-test is optimal for comparing the two means when the underlying distributions are normal and have equal variances. When the variances are not equal, then one uses the Behrens–Fisher statistic with the Satterthwaite approximation; see Welch (1937) and Satterthwaite (1946).

However, it is also well known that these assumptions are rarely met in practice, and an alternative is to use non-parametric approaches that impose fewer conditions on the underlying distributions. For the two-sample location problem, one is often interested in comparing the medians of the two populations, and the widely used Wilcoxon–Mann–Whitney (WMW) test is distribution-free, if the two distributions are continuous and have the same shape. Pratt (1964) has shown that the test does not maintain type I error if the variances are different. Fligner and Policello (1981) proposed a modified WMW statistic for unequal variances but assumed the two populations to be symmetric. Brunner and Munzel (2000) and Neubert and Brunner (2007) further extended the non-parametric statistics to more general situations by further relaxing the underlying assumptions.

In general, it is well recognized, for example, Hampel et al. (1986) showed that the parametric tests would be optimal, that is, they would be valid and have the optimal power, when the underlying assumptions are satisfied. On the other extreme are the non-parametric tests that have minimal assumptions on the underlying distribution, but a price is paid in terms of loss of power. However, when there are modest departures from the target family (usually normal), robust methods serve as a viable alternative as they provide significant gain in power while maintaining the type I error control. In the context of the two-sample problem, assuming the underlying distributions to be in the neighborhood of the normal family with equal variances, Srivastava et al. (1992) and Mudholkar et al. (1991) have proposed robust test procedures based on L-statistics. In particular, the approach based on trimmed means, which is based on the concept of trimming the extreme observations, seems appealing and has shown better operating characteristics particularly for the distributions that are heavier tailed than normal.

It would be a daunting task to test the underlying distributional assumption and the homogeneity of variances at each locus and then use the appropriate test statistic based on the results of those tests. Even if one performed that, one will have to account for the increased number of tests being performed and the conditional nature of the p-values in the second stage. Pounds and Rai (2009) adopted the concept of an assumption adequacy averaging approach, which incorporates an assessment of the assumption of normality and weighs the results of the two alternative approaches based on whether the assumption of normality is satisfied or not. However, their approach does not extend to the situations where the assumptions of normality and homoscedasticity may be violated simultaneously.

The hallmarks of a “good” robust procedure should be that which is able to control the type I error rate and provide significant gain in power when the underlying assumptions are violated. Also, it should be able to maintain the type I error control with minimal loss in power when the underlying assumptions hold. The literature is filled with robust test procedures proposed for comparing two populations, and it is worthwhile to note that majority of them assess the type I error control at the traditional level of = 0.05, with few exceptions, for example, see Lee (1995). Fagerland and Sandvik (2009) compared the robustness of five two-sample location tests for skewed distributions and concluded that the tests conducted at were less robust than those conducted at However, it may be noted that in the context of testing multiple hypotheses or in the context of controlling , it is essential that the performance of a test procedure be evaluated at more stringent type I error rates such as 0.001, 0.0001, or even lower as we are looking for p-values in the tail of the distribution. It is also seen that in the context of designing genomic studies, often the sample size justification uses more stringent level of such as 10^–4 or 10^–5, as seen in Chow et al. (2008) and Kang et al. (2009). It is not difficult to visualize that, in this context, a test that is valid (could be somewhat conservative without significant loss in power) at stringent levels of is likely to give fewer false positives than the tests that are unable to control the type I error. Thus, in conducting genetic analysis, one must ensure that the test statistic used is not only robust to the violations of the underlying model assumptions but also exhibits good operating characteristics even at stringent levels of type I error.

In this study, we propose a robust analog of the Behrens–Fisher statistic based on trimmed means that is robust to the violations of underlying assumptions of normality and homogeneity. We use the asymptotic normality of trimmed means, derived by Huber (1970), to obtain the asymptotic distribution of the proposed test statistic. Then, using the first two moments of the proposed test statistics and regression methods, we obtain finite sample approximations to make the statistics useful in small sample sizes. The proposed statistic provides for a robust alternative for comparing two distributions that may not be normal and may be heteroscedastic. It is interesting to note that the proposed test statistic has the best performance for skewed distributions as well. However, in the context of high-throughput studies, the results based on the p-value approach will also be affected by the correlations among the test statistics by virtue of the correlation among genes. In this development, we focus our attention on developing a test statistic that will provide valid p-values when the assumptions of normality and homoscedasticity are violated with the understanding that one could use the approaches described in Benjamini and Yekutieli (2007) or Sun and Cai (2009) to conduct inference with correlated p-values. In Section 2, we provide the details of the motivating example. In Section 3, the background information regarding trimmed means and their asymptotic properties are discussed. In Section 4, a brief account of the Behrens–Fisher statistic for two samples is provided. In Section 5, we propose the new test statistic and obtain null distribution approximations for the proposed test statistic for finite samples. In Section 6, the details of the simulation study to evaluate the performance of the proposed test statistic in terms of type I error control and power properties and its comparison with the existing approaches are provided. In Section 7, the usefulness of the proposed test procedure is demonstrated by applying it to the DNA methylation data. Section 8 is dedicated to discussions and miscellaneous comments.

2 Motivating Example

Teschendorff et al., (2010) conducted a study to investigate the mechanism of diabetic nephropathy by comparing 27580 markers from a genome-wide methylation array between cases and controls. These data were submitted by the authors to the NCBI Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo) under accession nos. GSE20067 and can be easily downloaded. There were 97 cases who had type 1 diabetes (T1D) and nephropathy and 98 controls who had T1D but with no evidence of renal disease. The purpose of this study was to identify the markers that would be differentially expressed between the two groups. Cao et al. (2013) used the two-sample t-test and converted them to p-values to identify the proportions of DNA methylations that were different between the two groups. They further showed that the monotonicity assumption required for the validity of the -based approach was violated when the underlying assumption of normality and/or homoscedasticity underlying the two-sample t-test did not hold.

Cao et al. (2013) used the raw proportions of methylation which range between 0 and 100%. However, in practice, the logit transformation is often used before applying any test procedure. The benefit of using this transformation is that it transforms the proportions on a scale that ranges from −∞ to ∞ and possibly achieves variance stabilization; see Box (1953). We checked the assumption of normality in cases and controls and equality of variances between the two groups for all markers using Shapiro and Wilk (1965) and the F-test (Box, 1953) at different significance levels on the raw and logit-transformed methylation data. The results for the logit transformation (raw) are reported, and the results for the raw data were slightly worse.

At the conventional level of there were 88% (92%) of markers that failed the normality test in either cases or controls or had unequal variances between cases and controls, and 30% (46%) of markers failed the normality test for both cases and control and also failed the equal variance test between them. At a more stringent level of , 68% (78%) of the markers failed the normality test in either cases or controls or had unequal variances, and 4% (13%) of the markers failed the normality test for both cases and controls and had unequal variance.

Thus, it is obvious that the assumption of normality and/or equality of variance, in general, is questionable, and robust methods that are robust to such violations should be used. In the following section, we provide the background of the trimmed means and their asymptotic distribution.

3 Trimmed Means and Asymptotic Results

3.1 Trimmed Means

In the univariate case for one sample problem, let be the order statistics of a random sample from a location scale population with the symmetric distribution function For an integer, denotes trimmed mean, by the following equation:where is the number of observations trimmed from each end. Tukey and McLaughlin (1963) were the first to propose a robust analog of Student’s -test by studentizing given by:whereis the Winsorized variance, and represents the “effective number of observations” obtained by trimming observations from each end of the ordered observations. They proposed to approximate the null distribution of by Student’s distribution with degrees of freedom

3.2 Asymptotic Results

Huber (1970) justified the studentization in light of the asymptotic normal distribution of the trimmed means. Specifically, he showed that when the underlying distribution is symmetric, continuous with mean , and variance and strictly increasing at points ξ, then asymptotically :where is the limit of the fraction , and

It may be noted that is nothing but a function in multiplied by that can be explicitly evaluated for a given Huber also showed that as , then:where and can be written as follows:

Now, if we fix the distribution where represents the normal cumulative distribution function (), then Eqs. 4, 6 can be written as follows:where and are functions of alone and relatively complex expressions, but they can be very accurately approximated by cubic polynomials. We computed these expressions over a fine grid of from 0 to 0.25, with a step size of 0.01, and used regression methods to find the best polynomial fits. We approximated the functions and since appears only indirectly in calculation through . Hence, we have:

Using the aforementioned asymptotic theory, Mudholkar et al. (1991) refined the approximation for one-sample trimmed -test and proposed a two-sample pooled trimmed statistic as a robust analog of the two-sample pooled -test. Now, we provide a brief account of the Behrens–Fisher statistic.

4 Two-Sample Behrens–Fisher Trimmed Statistic

In the univariate setting for the two-sample problem, let and be the order statistics from two random samples from a location/scale population with a symmetric distribution function for respectively. Under the assumption of normality, , the well-known Behrens–Fisher statistic is:where and are the sample variances and estimates and respectively. Now, let and for then it is well known that the rejection region is a function of and , and many test procedures to conduct the test have been proposed, for example, see Welch (1937, 1947, 1949), Satterthwaite (1946), Lee and Gurland (1975), Cochran and Cox (1950), Wald (1955), and Pagurova (1968). However, the most commonly implemented test procedure is due to Satterthwaite, which approximates the distribution of the Behrens–Fisher statistic in (11) with a -distribution with degrees of freedom given by:where is obtained by substituting the sample variance in place of where . Now, utilizing the background information presented in Sections 3, 4, we present the derivation of the robust Behrens–Fisher statistic.

5 Trimmed Statistic and its Null Distribution

5.1 Trimmed t Statistic

For the two-sample problem discussed in Section 4, Yuen (1974) substituted trimmed means and Winsorized variances in place of means and variances in (11) and proposed a robust analog of the Behrens–Fisher statistic as:where for are the trimmed means for the two samples obtained using (1), and their corresponding Winsorized variances are obtained using (3) suggested to approximate it with a -distribution with the obtained in a manner analogous to (12) with replaced by , for However, in their simulation studies, they assumed equal amount of trimming for both samples, and the simulation studies were limited to smaller sample sizes; the performance of the null distribution was evaluated at nominal levels of However, as noted before, it is important, particularly in the context of analyzing genomic expression data, that a good robust test should be able to maintain good type I error control even at more stringent levels of such as or It is also critical in the context of developing robust procedures that the test performs optimally when the underlying assumptions are not violated. That is, the test should have well-controlled type I error when the normality assumption holds true. The performance of the proposed statistic to Yuen’s statistic was evaluated (Supplementary Table S1) and will be discussed later, but it was seen that when the underlying distribution is normal, Yuen’s test statistic cannot control type I error for stringent levels of Also, its performance was quite poor for skewed distributions. Thus, it is important to obtain a test that is valid at stringent levels of and valid for a wide variety of underlying distributions including skewed distributions. We carried out that by modifying the test statistic and obtaining a better null distribution approximation for the proposed test statistic.

Now, assuming the underlying distribution to be normal, that is, the analogs of Eq. 7 and 8, after suppressing , can be written as:

Then, using obtained from Eq. 10 and the asymptotic normality result from Eq. 14, one can approximate the distribution ofwhich reduces to when there is no trimming, that is, . Then, using the results of the asymptotic theory as stated in Eqs. 13, 14 and in a manner analogous to Yuen (1974) but replacing by we propose the refined two-sample robust Behrens–Fisher statistic based on trimmed means as:where and are obtained using Eqs. 1, 3, and 10, respectively. It may be noted that, unlike Yuen, in our development, equal amount of trimming for the two samples is not required.

5.2 Null Distribution Approximation

In this section, we have combined the large sample theory of Section 3 and the results of a Monte Carlo study to develop a scaled Student’s approximation for the distribution of given in (16). Again, we have assumed that the underlying populations are normally distributed.

Furthermore, by dividing the numerator and denominator of (16) by we get the numerator to be approximately and approximating the denominator within the square root sign by a chi-square variate divided by its degrees of freedom leads to the following approximation:where is the chi-square variate with degrees of freedom such that and Then, following the logic of Satterthwaite approximation, the can be shown to be, after some algebraic simplification, the following:where Then, equating to and we obtain:where is an estimate of . Then, for moderate to large samples, the test statistic in (16) can be approximated by a Student’s distribution with degrees of freedom obtained in (18). It may be noted that when there is no trimming, that is, or the statistic reduces to the usual Behrens–Fisher statistic in (11) and the to in (12).

In order to render usable in small samples, a finite sample approximation to its null distribution was obtained by approximating it by a scaled Student’s distribution, that is, by This was done using an extensive simulation study in which two independent samples of sizes and ranging from 10 to 100 in increments of 10 with same means but different variances in the ratio of 0.1, 0.25, 1, 4, and 10 from normal populations were generated. Then, for each combination of sample sizes, variances, and each combination of and ranging from 0–25%, one hundred thousand samples were generated to obtain the empirical estimate of the variance of The scaling factor was then obtained by equating the empirical variances of with the variance of which is Regression methods were used to model the scaling factor as a function of and given in (18). The regression equation was obtained with the boundary condition that as either or Among various models considered, the following was found preferable on the grounds of accuracy and simplicity:That is, the test statistic (16) can be approximated by:

However, based on the simulations, it was seen that with the finite sample correction in (20) the strict type I error control could not be achieved; see Table 2 (Column 8).

Then, using simulation studies and the logic that a slightly conservative test can be obtained by making the degrees of freedom smaller (the tails would become a little bit heavier), we obtained the modified degrees of freedom using the following equation instead of (18):Then, obtain the scalar using the approach described before, as follows:That is, the test statistic in (16) can be approximated by :

The performance of the test statistic (16) using the null distribution approximations given in (20) and (23) were evaluated using extensive simulation studies described in the next section.

6 Simulations

6.1 Simulation Setup

Extensive simulation studies were performed to evaluate the performance of the proposed test statistics in (16), denoted by in terms of the type I error and power, and compared to the alternatives including the two-sample Behrens–Fisher test statistic denoted by and corresponding to Yuen’s test statistic; the non-parametric analog of Behrens–Fisher statistic (modified Mann–Whitney–Wilcoxon test), proposed by Fligner and Policello (1981) and denoted by the non-parametric asymptotic statistics by Neubert and Brunner (2007) denoted by and its permutation version denoted by the traditional two-sample t-test denoted by and the two-sample rank-sum test denoted by . Various trimming proportions, such as 0.05, 0.10, 0.15, and 0.20, were considered, and the trimmed t-tests corresponding to different trimming proportions were denoted by and respectively, for the approximation given in (20) and by and respectively, for the approximation given in (23). Similarly, for different trimming proportions were denoted by respectively.

It may be noted that the mWMW or its generalizations proposed by Neubert and Brunner (2007) test the general hypothesis of However, the results from these tests can be interpreted as test of medians when the two underlying distributions are identical, except for the shift in location. Furthermore, it is not difficult to see, as noted by Neubert and Brunner’s (2007), that testing of the aforementioned hypothesis would be consistent with testing equality of two means when the underlying distributions are symmetric with possibly different variances.

Simulations were conducted where a single hypothesis was simulated and examined at increasingly stringent significance levels of α to mimic the situation of testing multiple hypotheses but assuming the underlying distribution to be same for the two groups. Five families of distributions were considered: normal, contaminated normal, combined normal and uniform (contaminated with normal/uniform distribution), cauchy (all symmetric continuous distributions), and transformed beta (skewed continuous distribution). Although, the theory and derivation of the test statistic in (16) assumes the underlying distributions to be symmetric, it may be argued that, after appropriate trimming, the “middle” of the skewed distribution may also resemble the “middle” of the normal distribution, and it may be reasonable to apply and evaluate the performance of the trimmed test statistic for skewed distributions as well. In addition, in situations where hundreds and thousands of gene expressions are compared, it is likely that some underlying distributions may be skewed in real-life setting. So we included beta distribution in our simulation studies to mimic such a situation.

An independent simulation study was undertaken to compare the type I error control of and at levels of by simulating two samples from normal, contaminated normal, combined normal and uniform, cauchy, and beta distributions of various sample sizes, and varying from 20, 50, and 100, and the estimate of type I error estimates were obtained based on replicates. A selection of the results in presented in Supplementary Table S1).

From Supplementary Table S1, it is clear that Yuen’s approach, based on the Welch-type approximation, cannot control type I error for normal distribution at stringent levels of For example, when for various sample sizes and for 15% and 20% trimmings from the two samples the range of for is (0.01, 0.04), whereas for it is (1.10, 2.20), resulting in twice as many false positives than expected. This finding is consistent with Lee’s (1995) observation that Welch’s approximation results in significantly higher percentage errors than that of Welch–Aspin (Welch (1947); and Aspin (1948)) or Lee–Gurland (Lee and Gurland (1975)) approximations for comparing means of two normal populations with unequal variances.

For beta distribution (skewed), neither nor can control the type I error, but it is also very clear that performs significantly better compared to For example, for the choice of scale and shift parameters as specified in the table, for 15% and 20% trimmings from the two samples, the range of is (0.35, 2.90) and (8.10, 17.0), corresponding to and respectively.

The results for the contaminated normal and combined normal and uniform suggest that both methods are conservative with being somewhat more conservative than This could be because of our fine tuning of the null distribution to obtain a better control of the null distribution when the underlying populations are normal. These simulations suggest that may be a reasonable alternative to the Behrens–Fisher statistic in the presence of heterogeneity when the tests are conducted at typical levels of such as 0.05 or 0.01. However, more refined approximation, such as the one proposed in this study, for the null distribution would be needed when the focus is on conducting the tests at more stringent levels of such as or Since the performance of was not satisfactory in controlling type I error at stringent levels of it was not included in further simulation comparisons.

In another independent simulation study, from each of the distributions mentioned before, random samples of sizes (cases) and (controls) we simulated, and varied from 20 to 100 and type I error estimates were estimated based on 10⁵ replicates. The significance levels considered for our evaluations were α = 0.05, 0.01, 0.005, and 0.001. To estimate the power, 10⁵ replicates were simulated for each case-control data. The empirical type I error rates and power estimates were calculated as the proportion of replicates with p-values less than α. Extensive simulation studies corresponding to various combinations of parameters listed in Table 1 were conducted.

To further evaluate the performance of the tests three possible scenarios arise for various combinations of σ₁ and σ₂ for normal, contaminated normal and contaminated with normal/uniform distributions as described later.

6.1.1 Normal

Scenario I:σ₁ = σ₂ represents the case of two normal distributions with equal variances.

Scenario II:σ₁ < σ₂ represents the case where variance of the second population is larger, and in evaluating power, it corresponds to the situation that the variance is larger for the second population with larger mean, that is,

Scenario III:σ₁ > σ₂ represents the case where variance of the first population is larger, and in evaluating power, it corresponds to the situation that the variance is smaller for the second population with larger mean (μ₂ > μ₁ = 0).

6.1.2 Contaminated Normal

Scenario I:σ1 = σ₂ represents the case of two normal distributions with unequal variances (no contamination).

Scenario II:σ₁ < σ₂ represents the situation that the variance of the contaminated part of the distribution is larger, that is, contamination with the normal distribution with “outliers”.

Scenario III:σ₁ > σ₂ represents the situation that the variance of the contamination part of the distribution is smaller, that is, contamination with the normal distribution with “inliers”.

6.1.3 Combined Normal and Uniform

Scenario I:σ1 = σ₂ represents the case of two contaminated normal distributions, contaminated with normal/uniform (N/U) distribution, with equal variances.

Scenario II:σ₁ < σ₂ represents that the variance of contamination part of the distribution with N/U is larger, that is, contamination is done with “outliers” coming from N/U distribution.

Scenario III:σ₁ > σ₂ represents that the variance of contamination part of the distribution with N/U is smaller, that is, contamination is done with “inliers” coming from N/U distribution.

For cauchy and transformed beta, the distributions were simulated for the parameters given in Table 1.

6.2 Simulation Results

For ease of readability, we have reported the ratio of empirical estimate of type I error/expected level of significance, that is, , for all tables and figures reporting type I error results so that for a well-controlled test the ratio should be close to 1.

Based on our extensive simulation studies, it is clear that for all the distributions under study, the empirical type I error rates were better controlled for both trimmed tests, and corresponding to 15% trimming proportions compared to other trimming proportions (data not shown). In addition, from Table 2 (α = 0.001) and for α = 0.05 (data not shown), it is clearly seen that the type I error control for normal distribution for and are not maintained, and they get progressively worse with increasing stringent levels of α. Therefore, for all figures evaluating type I error and power properties, for the five distributions, we only included and its competitors and

6.2.1 Normal Distribution

Null Distribution: As seen from Figure 1A, it is clear that, in general, the type I error control at is well maintained with and being on the conservative side and being somewhat anti-conservative.

FIGURE 1

Power Properties: For Scenario I, with and (data not shown), it is seen that power for and are quite comparable with having slight advantage, as would be expected. does slightly worse followed by particularly when the sample size for one of the groups is smaller. For Scenario II, Figure 1B, taking and and once again the performance of and are quite comparable with doing slightly better. The performance of and are comparable with doing slightly better than For Scenario III, Figure 1C, taking and and the power for all tests is very low but clearly dominates all other tests, and the performance of is the worst, but not by much.

6.2.2 Combined Normal and Uniform

Null Distribution: From Figure 2A, it is clear that has the best type I error control as the percentage of times R > 1.5 among all cases for the tests. and are 0.89, 3.91, 20.89 and 6.40, respectively, when the true nominal level is Power Properties: For Scenario I, assuming and (data not shown), the performance of was comparable to with slight advantage for The power estimates are higher for , but it fails to control type I error. The performance of is the worst. For Scenario II, Figure 2B, and it is seen that, in general, has more power than and and is comparable to but does not control type I error well. For Scenario III, Figure 2C, and the inliers problem, the trimmed test does worse than and but the power, in general, is low. One should also keep in mind that the problem of “inliers” is less common in practice, and the trimmed test is designed to provide protection against “outliers” ; in the presence of outliers, the trimmed test performs well.

FIGURE 2

6.2.3 Cauchy Distribution

For cauchy distribution, from Figure 3A, it is very clear that and are conservative but could be anti-conservative particularly for small sample sizes. From Figure 3B for power estimates, it is clear that the performance of is the worst, and the performances of and are comparable with having a slight advantage. performs the best but one must keep in mind that often the type I error is not controlled, especially for smaller sample sizes.

FIGURE 3

6.2.4 Skewed Transformed Beta

For skewed transformed beta distribution, from Figure 4A, it is very clear that none of the four tests can control type I error well. The percentage of times the type I errors exceeds 1.5 threshold for and are 34.44, 50.00, 100.00, and 47.78, respectively. It may be noted that the percentage of times the type I errors exceeds the threshold of 2 for and are 14.44, 31.11, 93.33, and 33.33, respectively However, from Figure 4B, it is very clear that the performance of and are very comparable with performing slightly better than which is consistent with the observation noted in Fagerland and Sandvik (2009). It may be noted that both and have higher power estimates, but the type I error is poorly controlled.

FIGURE 4

7 Application to DNA Methylation Data

We downloaded the data from the NCBI Gene Expression Omnibus website and applied all the approaches discussed before to the example to evaluate their relative performances.

Figure 5 shows the histogram of p-values with the density estimates and the estimated FDR as a function of p-value cutoffs. The estimated non-null proportion are 0, 0.09, 0.28, 0.00005, and 0 corresponding to and, respectively. From Figure 5, it is clear that the estimated for and are not monotone functions of the p-value cutoff. Thus, conclusions drawn from such analysis would be misleading. The estimated of and are 1 for all p-value cutoffs, but and output 24 and 17 p-values of exactly 0. For the DNA methylation data, if we set the significance level at α = 0.0001, then the number of SNPs identified to be significantly associated with the phenotype were 6, 7, 1, 24, and 17 corresponding to and respectively. A histogram of the five most significant biomarkers for the two groups obtained based on is presented in Figure 6 to visually examine if the distributions of these markers can be perceived to be significantly different. For the marker cg15121304, it is clear that both distributions are skewed and probably no outliers, then based on our simulations results, we would expect and to perform similarly as seen with p-values of 3 × 10⁻⁵ and 1.49 × 10⁻⁵, respectively. Furthermore, we expect and to produce highly significant p-values as they are not able to control type I error, which is confirmed with p-values of 0 for both tests. The distributions for the 4th and 5th marker (cg00491404 and cg16098170) are highly skewed and possibly have outliers, and in such situations, as seen from simulations, the p-values based on and tests are essentially 0 as we expected as they cannot control type I error rate (false positives), whereas those based on and are not significant at level suggesting that the results based on these two methods are similar and probably more conservative and believable. Of course, realizing that, in general, the test would have significantly lower power than when the underlying assumptions are violated.

FIGURE 5

FIGURE 6

8 Discussion

The proposed trimmed analog of the Behrens–Fisher statistic is robust in the sense that it can strictly maintain the type I error rate compared to the alternatives currently available in the literature and at the same time can provide significant gain in power when the underlying distributions may be in the neighborhood of normal with possibly unequal variances. However, it is possible that for some other underlying distributions and for some parameter combinations, other test procedures may outperform However, based on the simulation studies, which include a broad range of symmetric heavy-tailed and skewed distributions, the example clearly outperforms its competitors in controlling the type I error rate, even at very stringent levels of and has shown comparable power properties for a broad range of distributions. Thus, it provides for a viable alternative for comparing two distributions even when the assumption of normality or homoscedasticity may not hold. In the context of multiple hypotheses, it may not be feasible to test the assumption of normality and homoscedasticity simultaneously for all the hypotheses and then appropriately incorporate the findings in testing the hypothesis of interest using the most appropriate test. Thus, a procedure that can be implemented in broad settings that has reasonable robustness properties is needed and, we feel that the proposed trimmed statistic meets that need. Although, our simulation studies have focused on testing single hypothesis at stringent levels of , but since our test is on the conservative side (without much loss in power), it is not hard to visualize that by using the proposed test statistic one should be able to minimize false discoveries in the context of multiple hypotheses setting. It may be noted that the implementation of the trimmed test is straightforward and quick since we can use the well-known t-distribution with modified degrees of freedom. The computing time for and (based on 10,000 permutations) were 0.001995087, 0.00199604, 0.009974957, and 9.892611 s, respectively, for one marker with 98 controls and 97 cases.

In a genome wide association study of a continuous outcome, often, we are interested in testing if the continuous outcomes corresponding to three genotypes are same or not. Furthermore, in a gene expression analysis of k samples, based on multiple dose levels, we could be interested in knowing if the gene expressions among k-sample are different or not. To address these issues, the commonly used parametric method will be ANOVA analysis if the data follow normal distribution; otherwise, the alternative of ANOVA will be the Kruskal and Wallis, 1952. However, similar to the two-sample comparison discussed in the study, it may be perceived that the two most popular methods may not be able to maintain type I error rate at a given significance level and may lose significant statistical power when the underlying distributions may not be normal and possibly heteroscedastic. We are currently in the process of investigating it and extending our approach to k-sample heteroscedastic case.

Often, the comparison in the two-sample case or k-sample case needs to be adjusted for covariates of interest that may be associated with the phenotype of interest. We are currently in the process of developing approaches that would provide robust comparisons after adjusting for the covariates.

Although the motivation and presentation of the method lies in identifying genetic features different between two groups, it is also readily applicable to any epidemiology studies of comparing continuous variables between two groups and any clinical trial of comparing continuous responses for two treatments. We have implemented the proposed method in R program (Supplementary Note S1). The method can be easily applied to compare the continuous variables between two groups from one to hundreds of thousands of tests.

References

• 1

  AspinA. A. (1948). An Examination and Further Development of a Formula Arising in the Problem of Comparing Two Mean Values. Biometrika35, 88–96. 10.1093/biomet/35.1-2.88

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  BenjaminiY.YekutieliD. (2007). The Control of False Discovery Rate in Multiple Testing under Dependence. Ann. Stat.29, 1165–1188.

  • Google Scholar
  • View reference in article

• 3

  BoxG. E. P. (1953). Non-Normality and Tests on Variances. Biometrika40 (3/4), 318–335. 10.2307/2333350

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  BrunnerE.MunzelU. (2000). The Nonparametric Behrens-Fisher Problem: Asymptotic Theory and a Small-Sample Approximation. Biom. J.42, 17–25. 10.1002/(sici)1521-4036(200001)42:1<17::aid-bimj17>3.0.co;2-u

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  CaoH.SunW.KosorokM. R. (2013). The Optimal Power Puzzle: Scrutiny of the Monotone Likelihood Ratio assumption in Multiple Testing. Biometrika100 (2), 495–502. 10.1093/biomet/ast001

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  ChowS. C.ShaoJ.WangH. (2008). Sample Size Calculations in Clinical Research. Chapman & Hall/CRC, Taylor and Francis Group.

  • Google Scholar
  • View reference in article

• 7

  CochranW. G.CoxG. M. (1950). Experimental Designs. New York: John Wiley.

  • Google Scholar
  • View reference in article

• 8

  FagerlandM. W.SandvikL. (2009). Performance of Five Two-Sample Location Tests for Skewed Distributions with Unequal Variances. Contemp. Clin. Trials30, 490–496. 10.1016/j.cct.2009.06.007

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  FlignerM. A.PolicelloG. E. (1981). Robust Rank Procedures for the Behrens-Fisher Problem. J. Am. Stat. Assoc.76, 162–168. 10.1080/01621459.1981.10477623

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  HampelF. R.RonchettiE. M.RousseeuwP. J.StahelW. A. (1986). Robust Statistics: The Approach Based on Influence Functions. New York: John Wiley and Sons.

  • Google Scholar
  • View reference in article

• 11

  HochbergY.TamhaneA. (1987). Multiple Comparison Procedures. New York: Wiley.

  • Google Scholar
  • View reference in article

• 12

  HuberP. J. (1970). “Studentizing Robust Estimates,” in Nonparametric Techniques in Statistical Inference. Editor PuriM. L. (Cambridge, England: Cambridge University Press), 453–463.

  • Google Scholar
  • View reference in article

• 13

  JanssenA. (1997). Studentized Permutation Tests for non-i.i.D. Hypotheses and the Generalized Behrens-Fisher Problem. Stat. Probab. Lett.36, 9–21. 10.1016/s0167-7152(97)00043-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  KangG.YeK.LiuN.AllisonD. B.GaoG. (2009). Weighted Multiple Hypothesis Testing Procedures. Stat. Appl. Genet. Mol. Biol.8, 1–22. 10.2202/1544-6115.1437

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  KruskalW. H.WallisW. A. (1952). Use of Ranks in One-Criterion Variance Analysis. J. Am. Stat. Assoc.47 (260), 583–621. 10.1080/01621459.1952.10483441

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  LeeA. F. S. (1995). Coefficients of lee-gurland Two-Sample Test on normal Means. Commun. Stat. - Theor. Methods24 (7), 1743–1768. 10.1080/03610929508831583

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  LeeA. F. S.GurlandJ. (1975). Size and Power of Tests for equality of Means of Two normal Populations with Unequal Variances. J. Am. Stat. Assoc.70, 933–941. 10.1080/01621459.1975.10480326

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  MudholkarA.MudholkarG. S.SrivastavaD. K. (1991). A Construction and Appraisal of Pooled Trimmed-Tstatistics. Commun. Stat. - Theor. Methods20, 1345–1359. 10.1080/03610929108830569

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  NeubertK.BrunnerE. (2007). A Studentized Permutation Test for the Non-parametric Behrens-Fisher Problem. Comput. Stat. Data Anal.51, 5192–5204. 10.1016/j.csda.2006.05.024

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  PagurovaV. I. (1968). On a Comparison of Means of Two normal Samples. Theor. Probab. Appl.13, 527–534. 10.1137/1113069

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  PoundsS.RaiS. N. (2009). Assumption Adequacy Averaging as a Concept for Developing More Robust Methods for Differential Gene Expression Analysis. Comput. Stat. Data Anal.53 (5), 1604–1612. 10.1016/j.csda.2008.05.010

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  RobinsJ. M.van DER VaartA.VenturaV. (2000). Asymptotic Distribution ofPValues in Composite Null Models. J. Am. Stat. Assoc.95, 1143–1156. 10.1080/01621459.2000.10474310

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  SatterthwaiteF. E. (1946). An Approximate Distribution of Estimates of Variance Components. Biometrics Bull.2, 110–114. 10.2307/3002019

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  ShapiroS. S.WilkM. B. (1965). An Analysis of Variance Test for Normality (Complete Samples). Biometrika52 (3–4), 591–611. 10.1093/biomet/52.3-4.591

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  SrivastavaD. K.MudholkarG. S.MudholkarA. (1992). Assessing the Significance of Difference between Two Quick Estimates of Location. J. Appl. Stat.19 (3), 405–416. 10.1080/02664769200000036

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  StoreyJ. D. (2002). A Direct Approach to False Discovery Rates. J. R. Stat. Soc. Ser. B Stat. Methodol.64, 479–498. 10.1111/1467-9868.00346

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  SunW.Tony CaiT. (2009). Large-scale Multiple Testing under Dependence. J. R. Stat. Soc. Ser. B Stat. Methodol.71, 393–424. 10.1111/j.1467-9868.2008.00694.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  TeschendorffA. E.MenonU.Gentry-MaharajA.RamusS. J.WeisenbergerD. J.ShenH.et al (2010). Age-dependent DNA Methylation of Genes that Are Suppressed in Stem Cells Is a Hallmark of Cancer. Genome Res.20 (4), 440–446. 10.1101/gr.103606.109

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  TukeyJ. W.McLaughlinD. H. (1963). Less Vulnerable Confidence and Significance Procedures for Location Based on a Single Sample: Trimming/Winsorization I. Sankhya, A25, 331–352.

  • Google Scholar
  • View reference in article

• 30

  WaldA. (1955). Selected Papers in Statistics and Probability. New York: McGraw-Hill.

  • Google Scholar
  • View reference in article

• 31

  WelchB. L. (1949). Further Note on Mrs. Aspin’s Tables and on Certain Approximations to the Tabled Function. Biometrika36, 293–296.

  • Google Scholar

• 32

  WelchB. L. (1947). The Generalization of “Student’s” Problem when Several Difference Population Variances Are Involved. Biometrika34 (1-2), 28–35. 10.2307/2332510

  • CrossRef
  • Google Scholar

• 33

  WelchB. L. (1937). The Significance of the Difference between Two Means when the Population Variances Are Unequal. Biometrika29, 350–362.

  • Google Scholar

• 34

  YuenK. K. (1974). The Two-Sample Trimmed T for Unequal Population Variances. Biometrika61 (1), 165–170. 10.2307/2334299

  • CrossRef
  • Google Scholar