Stratified sampling in highly polluted data as an effective and reliable alternative to high breakdown point estimators

1. Introduction

Studies involving several real-life scenarios contains interest measurements that vary from the rest of the results, which are often called measurement errors. Some of these inconsistent measurements, however, are genuinely calculated observations that exhibit characters that differ from the rest of the data; we find those measurements to belong to another population. The presence of incoherent observations called outliers affects many statistical procedures and the validity of their results, especially when the distribution underlying them is normal. Nonetheless, there are many outlier concepts (Barnett & Lewis, 1995; Todorov et al., 2011; Ramsy & Ramsy, 2007); outliers are pure measurements (not necessary errors) that are incompatible with most of the data and do not behave in similar mode.

Outliers can occur for reasons such as; errors in data transmission or transcription, data coding and entry, ineffective survey design and processes, defective equipment or irregular spikes or, time-series data droppings. They can also occur because of natural causes as certain individuals can demonstrate characters well above or below the average. For example, some individuals may die earlier or live longer in life than most of their contemporaries in a given community (Fisher & Waclawski, 2009). Nevertheless, extreme natural events are uncommon, but in fact, they are representative of nature since they are correctly measured. Non-representative outliers; must be detected and rectified throughout the data editing process. Outliers with multiple measurements and/or variables are more likely to occur in datasets. A quick visual inspection still does not make them easy to detect (Hubert et al., 2008). Outliers introduce effects of masking which obscures real effects. Even a small cluster of observations may inflate the empirical covariance matrix to introduce the non-existent effect called swamping which makes some observations look like outliers when the true outliers pull the empirical matrix away from the non-outliers (van der Linde & Houle, 2006; Wilcox, 2011; Elliot & Stettler, 2006).

When assumptions of normality are broken due to the existence of outliers, the classical multivariate parameter estimation methods do not work well, because they are focused on; the empirical means, covariance and correlation matrices, and the least square fitting; all of which are highly influenced by only a few outliers that make them unstable (Croux et al., 1994). When addressing outliers’ problems, it is sometimes tempting to replace extreme values with plausible ones assumed by the researcher or remove them outright from the data set; however, extreme values are only removed if there are sufficient technical reasons for their aberrant behavior (Burke, 1988). Statistical criteria alone are not adequate to exclude extreme values from a data set, says Béguin and Hulliger (2008). For example, in fraud detection, network intrusion detection, and clinical disease diagnosis, outliers send a signal in many organizations and processes such as financial institutions (Murugavel & Punithavalli, 2011). The rest of the paper is organized as follows; robust estimation was addressed in Section 2, accurate estimation on the high breakdown point was addressed in Section 3. Partitioning algorithm discussed in Section 4. We discuss the simulation procedures used in the work in Section 5, we discuss parameter estimation in Section 6, results in Section 7, discussion of findings in Section 8, and the paper ended with the concluding remarks in Section 9.

2. Robust estimation

The classical parameter estimator is highly based on assumptions that in reality are difficult to satisfy. The classical mean and covariance measurements are in a given multivariate data set X n × p , for the i th observation denoted by x i = ( x i ⁢ 1 , … , x i ⁢ p ) ′ ; where x ¯ = 1 n ⁢ ∑ i = 1 n x i and S 2 = 1 ( n - 1 ) ⁢ ∑ i = 1 n ( x i - x ¯ ) ⁢ ( x i - x ¯ ) ′ are the mean and variance respectively.

Such measures are seriously influenced by the inclusion in the data set of only a small fraction of outliers (Franklin et al., 2000), which makes the classical estimator unreliable in managing most data sets in real life. Due to the increasingly complex nature of several real-life situations, more particularly in surveys, efficient outlier detection and robust estimation methods are required to handle inconsistent measurements. Ignoring or improperly handling them can cause the forecasting model to neither represent the bulk of the data nor the outliers (Hubert et al., 2008; Davey & Savla, 2010; Rousseeuw & Hubert, 2018; Udea, 2009). An estimator or statistical technique is reliable because it offers valuable information even though it does not apply any of the premises used to validate the estimation process (Fox & Weisberg, 2010). The robust statistical methods are aimed at finding a fit close to the fit we would have without outliers (Verdonck et al., 2011), generating efficient estimates by reducing outliers’ effect while minimizing bias. Robust estimators are commonly used for estimating parameters in different distributions of polluted data; Aydina et al. (2018) looked at a robust estimation of scale and position parameters in shifted Gompertz distribution. The widely used robust estimators include S-estimates (Rousseeuw & Yohai, 1984), MM estimates (Yohai, 1987), the least median of squares (LMS) estimators and the least trimmed square (LTS) estimators (Rousseeuw, 1984) with a breakdown point of up to 50% but low asymptotic efficiency. Other estimators are the high breakdown point combined with Gaussian asymptotic distribution, generalized S-estimators (Rousseeuw & Yohai, 1984) achieving efficiencies of up to 33%. The MM-estimators introduced by Yohai (1987) and the Yohai and Zamar (1988) T-estimators that can achieve arbitrarily high efficiency without losing their 50% breakdown point, but pay for it with increased bias. Ahmed et al. (2014) compared the performance of several robust methods concerning their strengths and weaknesses in terms of efficiency, bias, and ability to accommodate high contamination levels. Many of the robust methods can accommodate outliers of up to 33% and retain efficiency while others can accommodate outliers of up to 50% but lose their efficiency. In practice, only a few estimators can accommodate up to 50% outliers while retaining their efficiency; however, they give highly biased estimates in most cases (Croux et al., 1994). Only when it is clear that the fitted model is accurate do we find high-breakdown point estimators, since these estimates do not allow the diagnosis of model misspecification (Cook et al., 1992).

3. Breakdown point

The breakdown point (Hampel, 1971) is a term used to measure the degree of robustness of an estimator and this is a very common quantitative characteristic. According to (Hubert & Debruyne, 2010; Rousseeuw & Hubert, 2018), the breakdown point calculates the smallest fraction of observations that must be replaced by arbitrary values to get the result to all limits, or similarly, it is the maximum fraction of outliers that can be added to a given sample without spoiling the result.

Donoho and Huber (1983) propose a finite sample breakdown point for a multivariate position estimator as follows; the cycle starts with a random sample X n , m that denotes the collection obtained by replacing m data points x i ⁢ 1 , … ⁢ x i ⁢ m including X n with arbitrary values, then a sample data points of size n , i.e. S n = ( x 1 , x 2 , … , x n ) is considered, and T n is used to denote regression estimator. The estimation of the coefficients of regression given by T n ⁢ ( S n ) and the estimator breakdown point at S n are defined by

ϵ n * ⁢ ( T n ; S n ) = 1 n ⁢ min ⁡ { m ∈ ( 1 , 2 , … , n ) : sup m ⁡ T n ⁢ ( X n ) - T n ⁢ ( X n , m ) = + ∞ }(1)

and for a multivariate estimator of scatter we have

ϵ n * ⁢ ( C n ; X n ) = 1 n ⁢ min ⁡ { m ∈ ( 1 , 2 , … , n ) : sup m ⁡ max i ⁡ { | log ⁡ ( λ i ⁢ ( C n ⁢ ( X n ) ) ) ⁢ log ⁡ ( λ i ⁢ ( C n ⁢ ( X n , m ) ) ) | } }(2)

For the eigenvalues of C n , that means the scatter estimator is broken if any of its own values can become arbitrary large or close to zero. High breakdown point estimators are robustness-based procedures, which are able to handle highly contaminated data, describing the ability of an estimator to withstand any degree of data contamination (Hampel, 1971). Two of the high breakdown point estimators considered for this analysis are the Least Median of Square Estimator (LMS) (Rousseeuw, 1984), which minimizes the median of square residuals in which Hampel (1975) expressed the definition of Least Median of Square (LMS) which minimizes square size.

med ⁡ { [ Y i - x i ′ ⁢ t ] 2 } , t ∈ ℝ p ⁢ and ⁢ 1 ⩽ i ⩽ n(3)

Rousseeuw (1984) showed that the estimator had a 50% breakdown point; it estimates β consistently with the β accuracy rate, which makes it highly efficient. In addition, Rousseeuw (1985) suggested the least trimmed square estimator (LTS), a robust alternative to the standard least square estimator and is given by,

r 2 = ∑ i = 1 h n ( Y i - x i ′ ⁢ t ) 2 := min , t ∈ ℝ p(4)

where ( r 2 ) 1 : n ⩽ ( r 2 ) 2 : n ⩽ … ⩽ ( r 2 ) n : n are the ordered squared residuals (the residuals are first squared and then ordered), this criteria does not count the largest squared residuals thereby allowing the LTS fit to steer clear of outliers. When h = [ n / 2 ] , LTS which is a β -consistent estimator locates that half of the observations which has the smallest estimated variance, and in that case it has a breakdown point of 50%. The LTS is insensitive to corruption due to outliers if the outliers make up less than 50% of the set (Mount et al., 2014). The LMS estimator on the other hand minimizes the mean squared residual; the advantage LTS has over LMS is that it is statistically more effective than LMS.

Our proposed approach assumes that the high concentration of pollutants divides the population into k strata or classes that are mutually exclusive. We consider three levels of contamination; 30%, 40%, and 50%; therefore, the partitions may have the following structure ( γ % ⩽ [ 50 ⁢ to ⁢ 70 ] % ⩽ δ % ) in which the proportion of outliers is below or above the minimum and maximum observations of the original data set. The percentage of the contamination γ and δ are such that; γ + δ ⩽ 50 percent. In this work, we assume that the contaminants divide the data into two groups which do not overlap ( k = 2 ) . We limit our method to one-directional contamination, which breaks the data into two mutual subpopulations of size N 1 and N 2 belonging to specific communities where N 1 + N 2 = N is of course stratified by the population.

5. Simulation

Simulations and parameter estimation was performed in R (R Core Team, 2018). We simulate a population of size 10,000 units of multivariate data set (one dependent and four explanatory variables), and draw random samples of sizes n 1 = 50 , n 2 = 200 and n 3 = 500 from the population and the response variable was contaminated at 30%, 40% and 50% levels yielding 27 scenarios. A random seed was set for all the simulation and sampling scenarios to ensure reproducibility of results, and we flag the outliers for identification and comparison purposes. The contaminated samples were analysed using two high breakdown point estimators, the least median of squares (LMS) and the least trimmed squares (LTS). For the proposed method (Nmeth), we first contaminate the 10,000 simulated units at 30%, 40%, and 50% levels. We believe that the high level of contamination splits the population into two distinct sub-populations of size N 1 and N 2 , with either 70/30%, 60/40% or 50/50% stratum population sizes respectively, where the first stratum contains the original measurements that are less or equal to the limit of the pre-contamination simulated data. The second stratum contains all the contaminated units, which are now larger than the maximum set in the original. Stratification is used subject to certain precision limitations, to achieve cost minimization or variability in sampling (Barcaroli 2014). Khan and Wesolowski (2019) considered the question of allocating samples in two and three stages with a view to optimizing multi-domain and population-efficient allocation. Their approach minimizes relative variances in all domains (controlled by given priority weights) as well as the overall relative variance under total constraints (expect) such an approach allows for solutions in multi-stage stratified simple random sampling without replacement (srswor) schemes which are direct generalization of the Neyman-type allocation in multi-stage stratified srswor scheme.

We use simple random sampling without replacement (srswor) within each stratum under Neyman allocation to obtain samples of size n 1 and n 2 where n 1 + n 2 = n . Our assumption is to optimize allocation by minimizing the relative variances in each stratum while taking in to account the sampling weights. The proposed method is easy to implement and can be used with any sample size that can be stratified into two or more non-overlapping groups containing at least two measurements. Table 1 shows the distribution of sample sizes drawn from the stratified population. Having selected the samples under the various contamination levels, we use the simple ordinary least squares (OLS) to fit the regression model to the data. Our assumption here is that the sample is optimally constructed to ensure the OLS is working effectively.

For the estimation of the breakdown point, we use linear regression with arbitrary error distribution to estimate parameters under LMS and LTS of Rousseeuw and Yohai (1984). The various models for regression are;

y i = x i ′ ⁢ β + ε i , for ⁢ i = 1 , … , n(5)

If y i is an observed response, x i is a vector of explanatory variables in dimensions p × 1 and β is a p × 1 vector of unknown parameters. Classically ε i , i = 1 , … , n is assumed to be independent and distributed in the same way as N ⁢ ( 0 , σ 2 ) , for any ${\sigma }^2>0$ , the regression class is defined as;

β ^ = arg ⁢ min β ∑ i = 1 k w i ⁢ ρ ⁢ ( | r ⁢ ( β ) | ( i ) )(6)

where | r ⁢ ( β ) | ( 1 ) ⩽ | r ⁢ ( β ) | ( 2 ) ⩽ … ⩽ | r ⁢ ( β ) | ( n ) are the ordered absolute values of the regression residuals r i = y i - x i ′ ⁢ β , w i is the weights, and ρ is a strictly increasing continuous function, such that ρ ⁢ ( 0 ) = 0 . The breakdown of the regression estimators Eq. (5) is equal to ( n - k ) / n if w i ⩾ 0 for i = 1 , … , n , $w_k>0$ for $k=\max \{i:w_i>0\}$ , the index k is within the boundaries ( n + p + 1 ) 2 ⩽ k ⩽ n - p - 1 , n ⩾ 3 ⁢ ( p + 1 ) and the general position of the data points x i ∈ R p .

In this section, we present calculated results from all the simulation scenarios used to test LMS, LTS and NewMeth’s efficiency according to the proposed design in Table 2. We fit regression lines to the data using all three techniques and, the Akaike Knowledge Criterion (AIC) calculated for each simulation scenario in order to see which model works best. We used a linear regression model to determine the intensity of each model that is highly resistant to outliers.

8. Discussion

Table 2 results show the efficiency of the two high breakdown point estimators (LMS and LTS) and our proposed method (Nmeth) at different contamination rates and sample sizes. Based on the various scenarios, we compare the AIC values for each estimator; with sample size 50 units, the output of both LMS and LTS for 30%, 40% and 50% contamination rates are not substantially different, they remain the same with minimal variations.

For sample size of 50 units, our proposed method with 30% contamination did not perform well since AIC Nmeth = 704.7168 $>{\text{AIC}}_{\text{LMS}}=$ 701.7 and AIC LTS = 702.712 , but with 40% contamination level AIC Nmeth = 714.1782 717.3 and AIC LTS = 715, and with 50% contamination level; AIC Nmeth = 702.4988 733.8 and AIC LTS = 730.9. Therefore, our proposed approach is better suited to the data than both LMS and LTS, meaning our method performed better with contamination levels of 40% and 50%. LMS and LTS performed better than out-of-the-box approach under 30% contamination for the sample size of 200 units, the performance of LMS and LTS however does not vary significantly from each other. With higher levels of contamination for the sample size 200 units, that is, 40% contamination level; AIC Nmeth = 2803.784 2852.2 and AIC LTS = 2852.9 and 50% contamination level; AIC Nmeth = 2815.43 2911.1 and AIC LTS = 2897.8, showing that our method performed better than the two robust high breakdown point estimators did, but both performed equally the same. Similarly, the output trend for the sample size of 500 units remains the same with our method performing worst compared to the two methods with 30% contamination where AIC Nmeth = 7070.05 $>{\text{AIC}}_{\text{LMS}}=$ 6968.8 and AIC LTS = 6963.8. However, at 40% contamination level; AIC Nmeth = 7064.944 7151.8 and AIC LTS = 7138.5. Similarly, at 50% contamination level; AIC Nmeth = 6947.975 7296.5 and AIC LTS = 7284.7 indicating a better performance by the proposed method than the two high breakdown point estimators.

9. Concluding remarks

In this paper, we present alternative method to high breakdown point estimators for handling highly polluted data set by using stratified sampling with Neyman allocation, taking into account the variance within the stratum and the homogeneity of the measurements for optimum estimation. This study aimed to explore new, easy-to-use, accurate, and robust methods for addressing highly polluted data set to high breakdown point estimators that can effectively accommodate outliers of up to 50 per cent. The AIC values for sample sizes 50, 200, and 500 for LMS, LTS, and Nmeth indicate that LMS and LTS performed better than Nmeth when fitting a robust regression line to the 30% contaminated data, but LMS and LTS performed similarly. For 40% and 50% contamination rates, Nmeth performed better than LMS and LTS generating smaller AIC values under the sample sizes of 50, 200 and 500 units. It is necessary to remember that the LMS and LTS yield the same results in all variations of the sample sizes and contamination levels. The analysis shows that our proposed approach (Nmeth) worked well when compared with LMS and LTS at contamination rates between 40% and 50%. This suggests that besides the two high breakdown point estimators (LMS and LTS) we may use Nmeth as an alternative in addressing highly contaminated data collection. We evaluate the performance of our proposed method using only simulated data but it is expected to perform effectively in real-life scenarios and produce comparable outcomes. The accuracy and effectiveness of the proposed approach may be of interest to the statistical analysis community at large.