Data Depth Trimming Counterpart of the Classical 𝑡 (or 𝑇²) Procedure

Abstract

The classical 𝑡 (or 𝑇² in high dimensions) inference procedure for unknown mean 𝜇∶𝑋±𝑡_𝛼(𝑛−1𝑆_𝑛/√𝑛(or{𝜇∶𝑛(𝑥−𝜇)^′𝑆⁻¹(𝑥−𝜇)≤𝜒²_(1−𝛼)(𝑝)}) is so fundamental in statistics and so prevailing in practices; it is regarded as an optimal procedure in the mind of many practitioners. It this manuscript we present a new procedure based on data depth trimming and bootstrapping that can outperform the classical 𝑡 (or 𝑇² in high dimensions) confidence interval (or region) procedure.

1. Introduction

Let 𝑋^𝑛∶={𝑋₁,…,𝑋_𝑛} be a random sample from distribution 𝐹 with an unknown mean parameter 𝜇. The most prevailing procedure for estimating 𝜇 is the classical 𝑡-confidence interval. A 100(1−2𝛼)% confidence interval (CI) for 𝜇 and large 𝑛 is 𝑋±𝑡_𝛼𝑠(𝑛−1)√𝑛, (1.1) where 𝑋=(1/𝑛)∑^𝑛_𝑖=1𝑋_𝑖 is the standard sample mean, 𝑠=(1/(𝑛−1))∑^𝑛_𝑖=1(𝑋_𝑖−𝑋)² is the standard sample deviation, and 𝑡_𝑟(𝑁) is the 𝑟th upper quantile of a 𝑡 distribution with degrees of freedom 𝑁. The rule of thumb in most textbooks for the sample size 𝑛 is: 𝑛<15, do not use 𝑡 procedure, 15<𝑛≤40, do not use it if outliers present, use it if 𝑛>40. The procedure is based on the large sample property and the central limit theorem. So it is not exact but an approximation for large sample size 𝑛 and arbitrary population distribution.

In higher dimensions, the counterpart to procedure (1.1) is the celebrated Hotelling's 𝑇² procedure: A 100(1−𝛼)% confidence region for the unknown vector 𝜇 and large 𝑛 is the region: {𝜇∶𝑛(𝑥−𝜇)^′𝑆⁻¹(𝑥−𝜇)≤𝜒²_(𝛼)(𝑝)}, (1.2) where 𝑆 is the sample covariance matrix and 𝜒²_(𝛼)(𝑝) is the upper 𝛼th quantile of 𝜒² distribution with degrees of freedom 𝑝.

Procedure (1.1) and (1.2) are so prevailing in practices that in many practitioners, mind they are regarded as optimal and unbeatable procedure. Are they really unbeatable? In this manuscript we introduce a new procedure that can outperform these seemingly optimal procedures.

The rest of the paper is organized as follows. Section 2 introduces the new procedure and Section 3 conducts some simulation studies. The paper ends in Section 4 with some concluding remarks.

2. A New Procedure for the Unknown 𝜇

2.1. A Univariate Location Estimator

It is well known that the sample mean in the 𝑡 procedure is extreme sensitive to outliers, heavy tailed distributions, or contamination. The procedure therefore is not robust. So naturally, one would replace the sample mean with a robust counterpart. We will utilize a special univariate location estimator 𝜇 to replace the sample mean 𝑋 in the 𝑡 procedure.

Now we consider a special univariate “projection depth-trimmed mean” (PTM_𝛽) for 𝑋^𝑛 in ℝ¹, 𝛽>0(see Wu and Zuo [1] in ℝ¹, also see Zuo [2] for a multidimensional PTM_𝛽)

𝜇(𝑋^𝑛)∶=PTM_𝛽(𝑋^𝑛1)=𝑘^𝑘_𝑖=1𝑋_𝑗𝑖, (2.1) where 𝑗₁,…,𝑗_𝑘 are distinct numbers from {1,…,𝑛} such that PD(𝑋_𝑗𝑖,𝑋^𝑛)≥𝛽 for some 𝛽>0,𝑖=1,…,𝑘, and PD(𝑋_𝑖,𝑋^𝑛)∶=1/(1+|𝑋_𝑖−Med(𝑋^𝑛)|/MAD(𝑋^𝑛)); where Med(𝑋^𝑛) is standard sample median and MAD∶=med{|𝑋_𝑖−Med(𝑋^𝑛)|,𝑖=1,…,𝑛} is standard median of absolute deviations (see Zuo and Serfling [3] and Zuo [4] for the study of PD in high-dimensions).

Let 𝐹_𝑛 be the empirical distribution based on 𝑋^𝑛 which places mass 1/𝑛 at points 𝑋_𝑖, 𝑖=1,…,𝑛. We sometimes write 𝐹_𝑛 (for 𝑋^𝑛) for convenience. Let 𝐹 be the distribution of 𝑋_𝑖. Replacing 𝑋^𝑛 with 𝐹 in the above definition, we obtain the population version. For example, the popular version of PTM_𝛽 for 𝐹∈ℝ¹ is

PTM_𝛽∫(𝐹)=^{PD(𝑥,𝐹)>𝛽}𝑥𝑑𝐹(𝑥)∫_{PD(𝑥,𝐹)>𝛽}.𝑑𝐹(𝑥) (2.2)

2.2. The New Procedure

Let 𝑋^∗_𝑛={𝑋^∗₁,…,𝑋^∗_𝑛} be a random sample from the empirical distribution 𝐹_𝑛. It is often called a bootstrap sample. Let 𝑌^𝑚∶={𝑋^∗_𝑛1,…,𝑋^∗_𝑛𝑚} be 𝑚 bootstrap samples from 𝐹_𝑛.

We calculate 𝑦_𝑗∶=PTM_𝛽(𝑋^∗_𝑛𝑗) for 𝑗=1,…,𝑚. Now we calculate depth of 𝑦_𝑗 with respect to 𝑦^𝑚∶={𝑦₁,…,𝑦_𝑚}: PD(𝑦_𝑗,𝑦^𝑚) and then order 𝑦_𝑗 with respect to their depth from the smallest to the largest: 𝑦₍₁₎,…,𝑦_(𝑚) where PD(𝑦₍₁₎,𝑦^𝑚),≤,…,≤PD(𝑦_(𝑚),𝑦^𝑚).

Finally, we simply delete first 100⋅2𝛼⋅𝑚% points from 𝑦₍₁₎,…,𝑦_(𝑚). Then the interval (or closed convex hull in high-dimensions) formed by 𝑦_{(⌊100⋅2𝛼⋅𝑚%⌋+1)},…,𝑦_(𝑚) is our 100(1−2𝛼)% confidence interval for 𝜇, where ⌊⋅⌋ is the floor function.

3. Simulation Study

Now we conduct simulation study to examine the performance of the new and classical 𝑡 (or 𝑇²) procedure based on 2000 (replication) samples from various distribution 𝐹(including 𝑁(0,1),𝑡(3), and others). Set 𝛼=0.025 and 𝛽=0.078; we consider the combinations of 𝑛=100 with the bootstrap number 𝑚=300, 500, 1000, and 2000.

We will confine attention to the average length (or area) of the confidence interval (or region) from both procedures as well as their coverage frequency of true parameter 𝜇 (which is assumed to be the mean of the 𝐹), which ideally should be close to 95%. If both procedures can reach the nominal level 95%, then it is obviously better to have a shorter (or smaller) confidence interval (or region) or smaller average length (or area) of the intervals (or regions).

3.1. One Dimension

Table 1 lists the simulation results at the normal and 𝑡(3) distributions.

Table 1 

Average coverage (length) of 95% CI's by 𝑡 and PTM_𝛽.

Inspecting the table immediately reveals that the bootstrap number 𝑚 affects the average coverage of the new procedure, with the increase of 𝑚 it gets closer to the nominal level 95%, while the average length of intervals gets slightly larger. Of course, it does not affect the 𝑡 procedure which has nothing to do with bootstrap. Overall, both procedures are indeed (roughly) 95% procedure and the new one produces an interval on the average about 2%-3% shorter than that of the classical 𝑡 procedure even at 𝑁(0,1) case, and it becomes 12%-13% shorter in the 𝑡(3) case.

Figure 1 displays the typical single run results from two procedures based on 100 sample points from 𝑁(0,1) (a) and 𝑡(3) (b). We see even at 𝑁(0,1) case, the new procedure outperforms the classical 𝑡 procedure with a 95% confidence interval [−0.1422550,0.2463009]10% shorter than that of 𝑡[−0.1644447,0.2661397], both cover the target parameter 𝜇=0. At 𝑡(3) case, new procedure produces an interval [−0.0874566,0.357108]60% shorter than that of 𝑡∶[−0.1644486,0.9482213]. Both cover the target parameter 𝜇=0.

(a)

(b)

Figure 1 

95% confidence intervals by 𝑡 (red one) and new procedures (blue one) for the mean of 𝐹 based on 100 sample points from 𝐹: (a) 𝐹=𝑁(0,1) and (b) 𝐹=𝑡(3).

In our simulation studies, we also compare our new procedure with the existing bootstrap percentile confidence procedure (i.e., it orders means of 𝑚 bootstrap samples and then just to trim the upper and lower ⌈100⋅𝑚⋅𝛼%⌉ points, the left points form an interval which is called bootstrap percentile conference interval, where ⌈𝑥⌉ is the ceiling function of 𝑥), our new procedure also outperforms this one. But the later performs better than the classical 𝑡 procedure in term of the average length of intervals at the same confidence level.

Our experiments with 𝑛 also reveal that small 𝑛 (the real situation in practice) is in favor of our new procedure. Note that this is the exact case where it is difficult to determine if the data are close to normal and hence to decide if one is able to use the classical CI. This is what we expected since the classical CI is based on normal 𝐹 (or on the large sample property for large 𝑛). But this does not mean that the classical CI has an edge over the new procedure at really large sample size 𝑛 (say, 10,000) even for the perfect 𝑁(0,1) case.

In addition to the distributions we considered in Table 1, we also conduct simulation studies to compare the performance of the new and classical 𝑡 procedure at contaminated normal model: (1−𝜖)𝑁(0,1)+𝜖𝑁(𝜇,𝜎²) with different choices of 𝜖 and (𝜇,𝜎²) since we know in practice, there is never a pure (exact) 𝑁(0,1); we may have just a slight departure from the pure normal or some contamination. Our results reveal that the new procedure is overwhelmingly more robust than the classical 𝑡; this is what we would expect since the 𝑡 procedure depends on the sample mean which is notorious for its extreme sensitivity to outliers or contaminations. We also compare the performance of the two procedures at Cauchy distribution since we know that sample mean 𝑥 performs extremely well at symmetric light tailed distributions like 𝑁(0,1) but not so at heavy tailer ones like cauchy distribution.

We first display the typical single run results of 95% confidence intervals in Figure 2 to demonstrate the difference between the two procedures.

(a)

(b)

Figure 2 

95% CIs by 𝑡 (red one) and new procedures (blue one) for the mean of 𝐹 based on 100 sample points from 𝐹: (a) 𝐹=0.95𝑁(0,1)+0.5𝑁(1.5,0.12) and (b) 𝐹=𝑡(1).

Here in Figure 2, on the left-hand side are 95% CI's by 𝑡 (red one) and by our new procedure (blue one) at the model 0.95𝑁(0,1)+0.05𝑁(1.5,0.1²) with an interval from 𝑡: [−0.09723137,0.3047054] and from new procedure [−0.06083767,0.2763990] which is 16% longer than that of 𝑡. These intervals are supposed to estimating the mean parameter 𝜇 in this case is 𝜖⋅𝑢=0.075. So both intervals cover the unknown parameter 𝜇.

On the right-hand side are 95% CIs by 𝑡 (red one) and by new procedure (blue one) at the Cauchy distribution with an interval from 𝑡: [−19.02593,6.527279] and from new procedure: [−0.8909354,0.5884936] which is 94% shorter.

Of course, the single run results may not represent the overall performance of the two procedures. So we conduct a simulation over 2000 replications. The results are listed in Table 2.

Table 2 

Average coverage (length) of 95% CIs by 𝑡 and PTM_𝛽.

Inspecting the table immediately reveals that the classical 𝑡 procedure becomes useless in the heavy tailed Cauchy distribution case: exceeding the nominal level 95% and reaching 98% with an extremely wide confidence interval, no informative any more. At the same time, the new procedure can roughly reach the nominal level 95% (it is about 97%) and provide a meaningful estimation about the underlying unknown parameter. We list the results from the contaminated model with just 5% contamination to a pure 𝑁(0,1) model with the contamination also come from a normal distribution centered at 1.5 and with a small variance 0.01. Under such a potential real situation, the classical 𝑡   95% procedure becomes again useless since it can never reach the nominal level 95%, it is a roughly 96% procedure with an interval slightly longer than that of the new procedure, while the new procedure still is a reasonable 95% procedure with an interval on the average 2%-4% shorter than that of 𝑡 one.

3.2. Higher Dimensions

In higher dimensions, with the multivariate version of PTM and PD (see Zuo [4], Zuo [2]) it is straightforward to extend our new procedure described in Section 2. That is, with the 𝑚 bootstrap sample: 𝑌^𝑚={𝑋^∗_𝑛1,…,𝑋^∗_𝑛𝑚} we calculate 𝑦_𝑗∶=PTM_𝛽(𝑋^∗_𝑛𝑗) for 𝑗=1,…,𝑚. Then we calculate the projection depth of 𝑦_𝑗 with respect to 𝑦^𝑚∶={𝑦₁,…,𝑦_𝑚}: PD(𝑦_𝑗,𝑦^𝑚) and then we order 𝑦_𝑗's with respect to their depth from smallest to largest: 𝑦₍₁₎,…,𝑦_(𝑚) where PD(𝑦₍₁₎,𝑦^𝑚),≤,…,≤PD(𝑦_(𝑚),𝑦^𝑚). The final step is the same as before: trimming first 100⋅2𝛼⋅𝑚% points from 𝑦_𝑗's the left formed a convex hull, that is our 100(1−2𝛼)% confidence region for 𝜇. We will examine the performance of this one and the classical Hotelling's 𝑇² given in (1.2) in term of their average area of confidence regions as well as their coverage frequency of true parameter 𝜇 (which is assumed to be the mean of the 𝐹). The latter ideally should be close to 95%. If both procedures can reach the nominal level 95%, then it is obviously better to have a smaller confidence region or smaller average area of confidence regions.

We first display single run results of two procedures at bivariate standard normal distribution 𝑁₂(0,1) and bivariate 𝑡 distribution with 3 degrees of freedom 𝑡₂(3) in Figure 3.

(a)

(b)

Figure 3 

95% confidence regions by 𝑇2 (red one) and new procedures (blue one) for the mean of 𝐹 based on 100 sample points from 𝐹: (a) 𝐹=𝑁2(0,1), (b) 𝐹=𝑡2(3).

Of course, single run result may not represent the overall performance of the two procedures. To see if the single run results are repeatable now we list the average of coverage and the area of the confidence regions based on two procedures in 2000 replications in Table 3. Here we set 𝛽=0.1, 𝑛=100 and 𝛼=0.025.

Table 3 

Average coverage (length) of 95% confidence regions by 𝑇² and PTM_𝛽.

Inspecting the Table reveals that the two procedures are indeed (roughly) 95% confidence procedures. Therefore it make sense to compare their average area of confidence regions. The table entries show that the new procedure produce a confidence region on the average 11%-22% smaller than that of the classical Hotelling's 𝑇² procedure in term of area even at 𝑁₂(0,1). This becomes 32%-40% in 𝑡₂(3) case.

4. Concluding Remarks

From the last section we see that the new procedure has some advantages over the classical (seemingly optimal) procedures. But we know that we cannot get all the advantages of the new procedure for free. What kind of price we have to pay here? For all the advantage of the new procedures possess over the classical ones, the price it has to pay is the intensive computing in the implement of the procedure. In our simulation study, there are 4 million basic operations (the case 𝑛=100, replication 𝑅=2000, and bootstrap number 𝑚=2000). Computing the data depth in two or higher dimensions is very challenging. Fortunately, there is a R package (called ExPD2D) for the exact computation of projection depth of bivariate data already developed by Zuo and Ye [5] and is part of CRAN now. For high-dimensional computation, see Zuo [6]. In one dimension it is straightforward. One can compute the sample median in linear time (i.e., the worst case time complexity is 𝑂(𝑛)) by employing special technic (see any computer science Algorithm textbook), for further discussion about the property of related remedian, see H. Chen and Z. Chen [7]. Fortunately, in practices, only one replication is needed. Also with the everlasting advance in computing power, the computation burden should not be an excuse for not using a better procedure.

A natural question is Why the new procedure has advantage over the classical one? The procedure clearly depends on bootstrap and data depth. Is it due to bootstrap or data depth? Who is the main contributor? If one just uses bootstrap, can one have some advantages? The answer for the latter is positive, Indeed, in our simulation we compare the classical one with the bootstrap percentile procedure, it reveals that the bootstrap percentile one does have some mild advantage over the classical one but still is inferior to our new procedure. So both bootstrap and data depth make contributions to the advantages of the new procedure. But remember, it is data depth that allow the bootstrap percentile procedure (which originally was defined only in one dimension) implementable in high-dimensions: to order sample bootstrap mean vectors. Without the data depth, it is impossible to implement the procedure in high-dimensions. So overall, it is data depth that makes the major contribution towards the advantages of the new procedure.

We also like to point out at this point that there is different new procedure introduced and studied in Zuo [8], where depth-weighted mean used in the procedure instead of the depth-trimmed mean used in our current procedure. However, our simulation studies indicate that our current new procedure is superior to the one in Zuo [8] which confines attention mainly to one dimension.

Our empirical evidence for the new procedure in one and higher dimensions is very promising, but we still need some theoretical developments and justifications, which is beyond scope of this paper and will be pursued elsewhere. A heuristic argument is because the bootstrap percentile confidence interval has advantage over the classic confidence interval procedure in term of at the same nominal level it can produce an asymptotically shorter interval (see Hall [9], and Falk and Kaufmann [10]). But the classical bootstrap percentile interval procedure is limited to one dimension, here we use data depth to ordering high-dimensional estimators so that we can extend the procedure to high-dimensions. The advantage of bootstrap percentile confidence interval carries on to high-dimensions.

One question left about our new procedure in practices is how does one choose the 𝛽 value? Well, there are at least two ways to deal with this 𝛽 value problem. First, one can chose a fixed value, our empirical experience indicates a value between 0.01-0.1 will serve most of our purposes. Or (second), dynamically choose 𝛽 value by minimizing some objective function which could be your interval length in our simulation case or variance in the efficiency evaluation case. With such a data dependant 𝛽, one natural question raised is: Is the theory in Zuo [2] established based on the fixed constant 𝛽 still holds? Fortunately, all still hold if we employ a more powerful tool (empirical process theory) from Pollard [11] or van der Vaart and Wellner [12] to handle this situation with a data dependent 𝛽.

There are a number of depth functions and related depth estimators (see Tukey [13], Liu [14], Zuo and Serfling [3], and Bai and He [15]), but among them projection depth function used here is the most favorite one (see Zuo [4, 16]). Furthermore, the computation of depth functions all are very challenging but we have some algorithm at hand for the projection depth function, this is yet another motivation for us to pick the projection depth function in this paper.

Finally, we comment that findings in this paper are consistent with the results obtained in Bai and Saranadasa (BS) [17] which shows the Effect of high-dimension, that is, there are better procedures than the classical inference procedures like Hotrlling's 𝑇² one which is inferior compared to other procedures like Dempster's nonexact test (Dempster [18]) and BS proposed test even for moderately large dimension and sample sizes.

Acknowledgment

This research was partially supported by NSF Grants DMS-0234078 and DMS-0501174.