In the univariate case for one sample problem, let X 1 < X 2 < … < X n be the order statistics of a random sample from a location scale population with the symmetric distribution function F ( ( x − θ ) / σ ) . For an integer, g < n denotes δ − trimmed mean, δ = g / n , by the following equation: X ˜ = ( X g + 1 + … + X n − g ) / ( n − 2 g ) , (1) where g is the number of observations trimmed from each end. Tukey and McLaughlin (1963) were the first to propose a robust analog of Student’s t -test by studentizing given by: t ˜ = ( X ˜ − θ ) / s ˜ T M , (2) where s ˜ T M 2 = [ ( g + 1 ) ( X g + 1 − X ˜ ) 2 + ( X g + 2 − X ˜ ) 2 + ⋯ + ( g + 1 ) ( X n − g − X ˜ ) 2 ] / h ( h − 1 )   , (3) is the Winsorized variance, and h = ( n − 2 g ) represents the “effective number of observations” obtained by trimming g observations from each end of the ordered observations. They proposed to approximate the null distribution of t ˜ by Student’s t distribution with ( h − 1 ) degrees of freedom ( d f ) .

Huber (1970) justified the studentization in light of the asymptotic normal distribution of the trimmed means. Specifically, he showed that when the underlying distribution F is symmetric, continuous with mean θ , and variance σ 2 and strictly increasing at points ± ξ, then asymptotically n → ∞ : n ( X ˜ − θ ) → N ( 0 , σ 2 ( δ ) ) , (4) where δ = F ( − ξ δ ) is the limit of the fraction g / n , and σ 2 ( δ ) = [ ∫ − ξ δ ξ δ x 2 d F + 2 δ ξ δ 2 ] / ( 1 − 2 δ ) 2 = b 2 ( δ ) σ 2 , (5)

It may be noted that σ 2 ( δ ) is nothing but a function in δ multiplied by σ 2 that can be explicitly evaluated for a given F . Huber also showed that as n → ∞ , then: n − 1 ( s ˜ 2 − b 2 ( δ ) σ 2 ) → N ( 0 , R 2 ( δ ) σ 4 ) , (6) where s ˜ 2 = [ ( g + 1 ) ( X g + 1 − X ˜ ) 2 + ( X g + 2 − X ˜ ) 2 + ⋯ + ( g + 1 ) ( X n − g − X ˜ ) 2 ] / n ( 1 − 2 δ ) 2 , and R 2 ( δ ) can be written as follows: R 2 ( δ ) ( 1 − 2 δ ) 4 σ 4 = ∫ − ξ δ ξ δ x 4 d F + 2 δ ( ξ δ 2 + 2 δ ξ δ f ( ξ δ ) ) 2 − [ ∫ − ξ δ ξ δ x 2 d F + 2 δ ( ξ δ 2 + 2 δ ξ δ f ( ξ δ ) ) ] 2 .

Now, if we fix the distribution F = Φ , where Φ represents the normal cumulative distribution function ( c d f ), then Eqs. 4, 6 can be written as follows: n ( X ˜ − θ ) → N ( 0 , b Φ 2 ( δ ) σ 2 ) , (7) n − 1 ( s ˜ 2 − b Φ 2 ( δ ) σ 2 ) → N ( 0 , R Φ 2 ( δ ) σ 4 ) , (8) where b Φ 2 ( δ ) and R Φ 2 ( δ ) 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 b Φ 2 ( δ ) and w Φ ∗ ( δ ) = b Φ 4 ( δ ) / R Φ 2 ( δ ) , since R Φ 2 ( δ ) appears only indirectly in calculation through w Φ ∗ ( δ ) . Hence, we have: b Φ 2 ( δ ) ≈ 1 + 0.48 δ + 1.21 δ 2 , (9) w Φ ( δ ) = ( n − 1 ) w Φ ∗ ( δ ) = ( n − 1 ) b Φ 4 ( δ ) / R Φ 2 ( δ ) ≈ ( n − 1 ) ( 0.5 − 1.62 δ + 1.91 δ 2 − 1.85 δ 3 ) . (10)

Using the aforementioned asymptotic theory, Mudholkar et al. (1991) refined the approximation for one-sample trimmed t -test and proposed a two-sample pooled trimmed t statistic as a robust analog of the two-sample pooled t -test. Now, we provide a brief account of the Behrens–Fisher 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: t ˜ Y , B F = ( X ˜ 1 − X ˜ 2 ) − ( θ 1 − θ 2 ) s ˜ 1 2 / h 1 + s ˜ 2 2 / h 2 , where X ˜ i for i = 1,2 are the δ i - trimmed means for the two samples obtained using (1), and their corresponding Winsorized variances ( s ˜ i 2 ) are obtained using (3) suggested to approximate it with a t -distribution with the d f obtained in a manner analogous to (12) with f i replaced by ( h i − 1 ) , for i = 1,2. 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 α = 0.01 , 0.05 ,  and  0.10. 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 α = 10 − 4 or 10 − 5 . 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 α = 10 − 3 and  10 − 4 . 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, F = Φ , the analogs of Eq. 7 and 8, after suppressing Φ , can be written as: n i ( X ˜ i − θ i )       →       N ( 0 , b 2 ( δ i ) σ i 2 ) ,   i = 1,2 , (13) n i − 1 ( s ˜ i 2 − b 2 ( δ i ) σ i 2 )       →       N ( 0 , R 2 ( δ i ) σ i 4 ) ,   i = 1,2. (14)

Then, using w i obtained from Eq. 10 and the asymptotic normality result from Eq. 14, one can approximate the distribution of 2 w i s ˜ i 2 / b i 2 ( δ i ) σ 2 ∼     χ 2 w i 2 , (15) which reduces to ( n i − 1 ) s i 2 / σ i 2 ∼ χ n i − 1 2 when there is no trimming, that is, δ i = 0 . Then, using the results of the asymptotic theory as stated in Eqs. 13, 14 and in a manner analogous to Yuen (1974) but replacing ( h i − 1 ) by ( 2 w i + 1 ) , we propose the refined two-sample robust Behrens–Fisher statistic based on trimmed means as: t ˜ B F = ( X ˜ 1 − X ˜ 2 ) − ( θ 1 − θ 2 ) s ˜ 1 2 ( 2 w 1 + 1 ) + s ˜ 2 2 ( 2 w 2 + 1 ) , (16) where X ˜ i , s ˜ i 2 , and w i , i = 1,2 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.

Remark: It may be noted that in the aforementioned derivation, the underlying distribution is fixed to normal, that is, F = Φ , to obtain the asymptotic distribution of the proposed test statistic, but the resulting test procedure is robust to the violations of the underlying assumptions of normality and homoscedasticity.

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 t approximation for the distribution of t ˜ B F given in (16). Again, we have assumed that the underlying populations are normally distributed.

Furthermore, by dividing the numerator and denominator of (16) by { b 1 2 σ 1 2 / ( 2 w 1 + 1 ) } + { b 2 2 σ 2 2 / ( 2 w 2 + 1 ) } , we get the numerator to be approximately N ( 0,1 ) , 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: [ { s ˜ 1 2 / ( 2 w 1 + 1 ) } + { s ˜ 2 2 / ( 2 w 2 + 1 ) } ] [ { b 1 2 σ 1 2 / ( 2 w 1 + 1 ) } + { b 2 2 σ 2 2 / ( 2 w 2 + 1 ) } ]   ≈   W ∗ d f W ∗ , (17) where W ∗ is the chi-square variate with degrees of freedom d f W ∗ such that E ( W ∗ / d f W ∗ ) = 1 and V a r ( W ∗ / d f W ∗ ) = 2 / ν ∗ . Then, following the logic of Satterthwaite approximation, the V a r ( W ∗ / d f W ∗ ) can be shown to be, after some algebraic simplification, the following: V a r ( W ∗ d f W ∗ ) = { R 1 2 σ 1 4 b 1 4 ( 2 w 1 + 1 ) 2 ( n 1 − 1 ) b 1 4 + R 2 2 σ 2 4 b 2 4 ( 2 w 2 + 1 ) 2 ( n 2 − 1 ) b 2 4 } ( b 1 2 σ 1 2 ( 2 w 1 + 1 ) + b 2 2 σ 2 2 ( 2 w 2 + 1 ) ) 2 , = ( b 1 2 σ 1 2 ( 2 w 1 + 1 ) ) 2 1 w 1 + ( b 2 2 σ 2 2 ( 2 w 2 + 1 ) ) 2 1 w 2 ( b 1 2 σ 1 2 ( 2 w 1 + 1 ) + b 2 2 σ 2 2 ( 2 w 2 + 1 ) ) 2 , = λ 2 w 1 + ( 1 − λ ) 2 w 2 , where λ = [ b 1 2 σ 1 2 / ( 2 w 1 + 1 ) ] / [ b 1 2 σ 1 2 / ( 2 w 1 + 1 )   +   b 2 2 σ 2 2 / ( 2 w 2 + 1 ) ] . Then, equating V a r ( W ∗ / d f W ∗ ) to 2 / ν ∗ and ν ∗ we obtain: 2 ν ⌢ ∗ = λ ^ 2 w 1 + ( 1 − λ ^ ) 2 w 2 , (18) where λ ^ = [ s ˜ 1 2 / ( 2 w 1 + 1 ) ] / [ s ˜ 1 2 / ( 2 w 1 + 1 )   +   s ˜ 2 2 / ( 2 w 2 + 1 ) ] is an estimate of λ . Then, for moderate to large samples, the test statistic in (16) can be approximated by a Student’s t distribution with ν ^ ∗ degrees of freedom obtained in (18). It may be noted that when there is no trimming, that is, δ 1 = δ 2 = 0 or g 1 = g 2 = 0 , the statistic reduces to the usual Behrens–Fisher statistic in (11) and the d f to f ^ in (12).

In order to render t ˜ B F usable in small samples, a finite sample approximation to its null distribution was obtained by approximating it by a scaled Student’s t distribution, that is, by A ∗ t ν ^ ∗ . This was done using an extensive simulation study in which two independent samples of sizes n 1 and n 2 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 δ 1 and δ 2 , ranging from 0–25%, one hundred thousand samples were generated to obtain the empirical estimate of the variance of A ∗ t ν ^ ∗ . The scaling factor was then obtained by equating the empirical variances of t ˜ B F with the variance of A ∗ t ν ^ ∗ , which is A ∗ 2 ν ^ ∗ / ( ν ^ ∗ − 2 ) . Regression methods were used to model the scaling factor A ∗ as a function of δ = ( δ 1 + δ 2 ) / 2 and ν ^ ∗ given in (18). The regression equation was obtained with the boundary condition that A ∗ → 1 as either δ → 0 or ν ^ ∗ → ∞ . Among various models considered, the following was found preferable on the grounds of accuracy and simplicity: A ∗ = 1 − 77 δ ν ^ ∗ + 1216 δ ν ^ ∗ 2 − 7186 δ ν ^ ∗ 3 + 17525 δ ν ^ ∗ 4 − 14881 δ ν ^ ∗ 5 . (19) That is, the test statistic (16) can be approximated by: t ˜ B F   ≈   A ∗ t ν ^ ∗ . (20)

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 ν m using the following equation instead of (18): 1 ν ⌢ m = λ ^ 2 w 1 + ( 1 − λ ^ ) 2 w 2 . (21) Then, obtain the scalar A m using the approach described before, as follows: A m = 1 − 29 δ ν ^ m + 215 δ 2 ν ^ m − 835 δ 3 ν ^ m − 273 δ ν ^ m ∗ 5 + 263 δ 2 ν ^ m ∗ 2 . (22) That is, the test statistic in (16) can be approximated by : t ˜ B F   ≈   A m t ν ⌢ m . (23)

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.

Extensive simulation studies were performed to evaluate the performance of the proposed test statistics t ˜ B F in (16), denoted by T R I M , in terms of the type I error and power, and compared to the alternatives including the two-sample Behrens–Fisher test statistic denoted by t B F and t ˜ Y , B F 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 m M W W ; the non-parametric asymptotic statistics by Neubert and Brunner (2007) denoted by T N B A and its permutation version denoted by T N B P ; the traditional two-sample t-test denoted by t ; and the two-sample rank-sum test denoted by W . Various trimming proportions, such as 0.05, 0.10, 0.15, and 0.20, were considered, and the trimmed t-tests ( t ˜ B F ) corresponding to different trimming proportions were denoted by T R I M 0.05 ,   T R I M 0.10 ,   T R I M 0.15 , and T R I M 0.20 , respectively, for the approximation given in (20) and by m T R I M 0.05 ,   m T R I M 0.10 ,   m T R I M 0.15 , and m T R I M 0.20 , respectively, for the approximation given in (23). Similarly, t ˜ Y , B F for different trimming proportions were denoted by t ˜ Y , B F ( 0.05 ) , t ˜ Y , B F ( 0.10 ) , t ˜ Y , B F ( 0.15 ) ,   and  t ˜ Y , B F ( 0.20 ) , respectively.

It may be noted that the mWMW or its generalizations proposed by Neubert and Brunner (2007) test the general hypothesis of P ( X < Y ) = 0.5. 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 t ˜ Y , B F and m T R I M at levels of α = 0.05 ,   0.01 ,   0.001 ,   and  0.0001 by simulating two samples from normal, contaminated normal, combined normal and uniform, cauchy, and beta distributions of various sample sizes, n 1 and n 2 , varying from 20, 50, and 100, and the estimate of type I error estimates were obtained based on 10 6 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 α = 0.0001 , for various sample sizes and for 15% and 20% trimmings from the two samples the range of R for m T R I M is (0.01, 0.04), whereas for t ˜ Y , B F , 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 m T R I M nor t ˜ Y , B F can control the type I error, but it is also very clear that m T R I M performs significantly better compared to t ˜ Y , B F . 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 R is (0.35, 2.90) and (8.10, 17.0), corresponding to m T R I M and t ˜ Y , B F , respectively.

The results for the contaminated normal and combined normal and uniform suggest that both methods are conservative with m T R I M being somewhat more conservative than t ˜ Y , B F . 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 t ˜ Y , B F 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 10 − 4 or 10 − 5 . Since the performance of t ˜ Y , B F 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 n 1 (cases) and n 2 (controls) we simulated, n 1 and n 2 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.

For ease of readability, we have reported the ratio of empirical estimate of type I error/expected level of significance, that is, R = α ^ / α , 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, T R I M and m T R I M , corresponding to 15% trimming proportions ( T R I M 0.15 ,   m T R I M 0.15 ) 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 T R I M 0.15 , m M W W , t , and W 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 m T R I M 0.15 and its competitors t B F , T N B A , and T N B P .

Remark: In practice, the choice of trimming proportions would be critical. Our recommendation, supported by our simulation studies, is to use 15% trimming proportions because trimming less than that provides results that would be similar to the normal setting (type I error control often not maintained) and trimming more than that results in more conservative tests and loss of power. In general, higher proportion of trimming would be recommended for settings where the underlying distribution may be very heavy tailed. However, this would not be known a priori and would be a daunting task to check at each loci, particularly, in the context of high-throughput data, but a 15% trimming provides a balance between not trimming enough to trimming too much and generally provides reasonably good results for all underlying distributions (including asymmetric distribution) studied in our simulation studies.

It may also be noted that we compared the null distributions at two levels of α = 0.05 and α = 0.001. We wanted to compare the performance of the test procedures at more stringent levels of α in the range of 10 − 4 and 10 − 5 , but this was not feasible for the permutation test as that would have required us to generate millions of samples to get reasonable estimate of type I error. However, we did estimate type I error control at stringent levels of α for m T R I M 0.15 and found that the type I error control was strictly maintained for all the symmetric distributions and was somewhat conservative (data not shown).

Extensive simulations studies corresponding to all combination of parameters mentioned in Table 1 were conducted, and the results were very similar, so a summary of the simulation results for each distribution, corresponding to a specific parameter combination, is discussed below. It is worth noting that the performance of the test procedures is relatively good for α = 0.05 , but it gets progressively worse as the type I error becomes more stringent. Thus, we discussed the results corresponding to α = 0.001 only. The results corresponding to α = 0.05 for the same parameter combinations are available from the authors on request. Also, the results corresponding to contaminated normal and combined normal and uniform were similar, so we chose to report the results corresponding to the combined normal and uniform distribution.