We consider a regression model with p independent variables (including the constant if there is an intercept)

where Y is an n × 1 vector of response, β is a p × 1 vector of unknown vector, X is an n × p matrix predictors, and ɛ is an vector of random errors with zero mean and variance matrix is σ²I_n. Linear regression method has been commonly used to analyze data because of its simplicity and applicability. However, this methodology is susceptible to disturbing effect caused by some outliers contained in the data. Therefore, the detection and examination of outliers or influential cases are important parts of data analysis. In detecting multiple outliers, masking phenomenon or swamping phenomenon usually occurs and makes it difficult to find true outliers.

Many methods resistant to masking effect and swamping effect are suggested. Marasinghe (1985) suggested a multistage procedure to improve the computational complexity required to obtain the most likely outlier subset. A generalized extreme studentized residual (GESR) procedure was proposed by Paul and Fung (1991). Kianifard and Swallow (1989, 1990) proposed to use test statistics based on recursive residuals for identifying outliers. Atkinson (1994) suggested a fast method for the detection of multiple outliers by using forward search. Pena and Yohai (1999) proposed a method which is fast and computationally feasible for a very large value of p. Jajo (2005) reviewed and compared outlier detection methods based on residuals.

One of general approaches for outlier detection is separating the data into a clean subset that is presumably free of outliers and a subset that contains all potential outliers. This approach is usually composed of several stages such as construction of a basic clean subset, calculation of residuals, outliers test and construction of a new clean subset. Hadi and Simonoff’s procedures (Hadi and Simonoff, 1993) adopted this approach and are summarized as follows.

A clean subset M of size s is constructed. The size of basic clean subset is s = int[(n + p − 1)/2]. Residuals d_i’s are calculated after fitting the model to the subset M. Residuals d_i’s correspond to the internally studentized residuals for the cases included in M and the scaled prediction errors for others:

Let |d|_(j) be the j^th order statistic of the |d_i| in ascending order. If |d|_(s+1) ≥ t_{(α/2(s+1),s−k)}, then the last (n − s) ordered cases are declared as outliers and the procedure stops. Otherwise, the procedure is repeated with a new clean subset consisting of the first (s+1) ordered observations. If n = s+1 no outlier is declared. Two methods for constructing the basic clean subset were suggested. The first method initially fits a model to the full data by using least squares (LS) estimate, orders the observations by their standardized residuals and takes the first p ordered observations as the early basic subset. After fitting a model to the early basic subset a new basic subset is formed by taking the first (p + 1) observations ordered by the residuals corresponding to (1.1). The basic subset is increased systematically until the size of the subset equals to int[(n + p − 1)/2]. The second method takes the backward-stepping approach (Rosner, 1975; Simonoff, 1988) based on the single linkage clustering (Hartigan 1981). A simulation study showed that the first method was more effective. Hadi and Simonoff (1993) demonstrated that their forward searching procedures are very effective and less affected by masking or swamping effects in comparison with other methods such as Marasinghe’s multistage procedure (1985), Paul and Fung’s (1991), and Kianifard and Swallow’s (1989, 1990) GESR procedure, OLS, and some high-breakdown methods like least median of squares (LMS) (Rousseeuw, 1984) and MM (Yohai, 1987).

The main purpose of this paper is to develop outliers detection methods which is invulnerable to masking and swamping effects. Newly suggested methods use least quantile squares (LQS) estimates for dividing the data into two groups, a clean set and an outlier set. LQS estimates are a generalization of LMS estimates. They have not been used as much as LMS because their breakdown points become small as the quartile value increases. But if the size of outliers is assumed to be fixed LQS estimates yield a good fit to the majority of data and residuals calculated from LQS estimates can be a reliable tool to detect outliers. In Section 2 the proposed methods are described with some examples. In Section 3 summary and concluding comments are presented.

For detecting outliers it is important to estimate a model immune from the influence of the outliers. The widely used estimator of parameters in the regression model is the LS estimator which involves the minimization of the sum of squared residuals. The drawback of the LS estimator is that it is sensitive to outliers in the data. A useful measure of the robustness of an estimator is the breakdown point of the estimator. Roughly speaking, it is the smallest fraction of arbitrary contamination that can cause the estimator to be infinite. For least squares, its breakdown points is 1/n which approaches zero as the sample size n increases.

Many robust estimators are suggested in the context of breakdown point. The most used high breakdown estimate is the LMS estimator, given by

where r_i(θ) is the the residual of the i^th observation. A generalization of the least median of squares estimator is the LQS estimate, which minimizes the k^th smallest squared residual.

We define k-LQS as

where $r_{(i)}^2(\theta )$ is the i^th order statistic of $r_{(i)}^2(\theta )$’s in descending order. Geometrically it corresponds to finding the narrowest strip covering k observations. The breakdown points of k-LQS decreases as k increases.

The computation of k-LQS estimates is complex. Several algorithms for LQS estimates have been proposed (Stromberg, 1993; Carrizosa and Plastria, 1995; Watson 1998). An LQS estimate involves determining Chebyshev solutions on all set of n+1 components. Stromberg (1993) used least squares and developed an algorithm to compute exact LQS estimates by considering θ̂′s based on all _nC_p+1 possible subsets of size p + 1. Watson (1998) provided alternative algorithms by using Gaussian elimination and modified linear programming.

We now propose new procedures for the identification of outliers in the regression by using k-LQS estimates. Residuals calculated from LQS estimates are used to separate the data into a clean subset and an outliers subset. The processes considered in this paper including Hadi and Simonoff’s perform the procedure for identifying outliers twice. One is for organizing a clean subset, and the other is for detecting candidates of outliers to be tested. Hadi and Simonoff used LS residuals in both procedures, but the two procedures need not be the same.

We suggest to use LQS residuals for constructing a clean set of Mand then to use LS residuals d_i’s in (1.1) for detecting and testing outliers. The test is performed using a t-distribution in this procedure. If LQS residuals are used in computing d_i’s in (1.1), the distribution of d_i’s under the null hypothesis is not clear. Test procedure using LQS residuals needs p-values calculated by Monte Carlo simulations. The proposed procedures use k-LQS estimates using all observations in building a new clean subset, a new clean subset is formed independently from the previous clean subsets including the basic subset. However, Hadi and Simonoff’s procedure should organize a good basic subset because it uses the LS residuals calculated from the fitted model of the previous clean subset to build a new clean subset. The proposed methods are independent of the results of the previous stages and are expected to be more robust than Hadi and Simonoff’s procedure in building a clean set.

Three newly suggested procedures denoted as S1, S2, and S3 respectively are described as follows. Let “M1” denote Hadi and Simonoff’s procedure using the first method for building a basic clean subset.

• S1: This method differ from M1 only in finding the basic clean subset. The basic clean subset is the set of observations with the smallest squared residuals calculated from k-LQS regression where k = int[(n + p − 1)/2].

The second method uses LQS estimates for constructing both a basic clean subset and new clean subsets.

• S2: The basic clean subset is obtained by the same method used in S1. After outlyingness test, if needed, a new clean subset of size k is constructed by using k-LQS estimates.

As a forward searching method neither LS, used in M1 nor LQS is perfect to avoid masking and swamping. The third method is designed to be relatively resistant to both masking and swamping effects which may occur in the process of S2. Let M_s denote the clean set of size s and $M_s^c$ be the complement of M_s, the set of potential outliers of size (n−s). We calculate ${\gamma }_s=n(M_s^c\cap M_{s+1}^c)/n(M_{s+1}^c)$ where n(M_i) is the size of M_i. If γ_s is small enough, it is highly probable that masking or swamping occurred. Thus if γ_s is small, both two subsets M_s and M_s+1 are used as clean sets. The Step 3 for S3 is described as follows.

• S3: Find a basic clean subset by using LQS estimates, then compute and order residuals according to (1.1). If |d|_(s+1) ≥ t_{(α/2(s+1),s−k)}, or |d|_(s+1)< t_{(α/2(s+1),s−k)} and γ_s ≥ δ for some fixed value of δ (we use δ = 0.5), then form a new subset M_s+1 by taking the (s + 1) smallest residual observations calculated from (s+1)-LQS regression. Otherwise, continue S2 procedure with a new subset M_s+1 and also perform M1 procedure with M_s as a basic clean subset. The final outliers are determined as the union of outliers detected by two procedures.

We now show an example to illustrate the proposed procedures. Columns x₁ and y₁ in the Table 1 have 25 observations. The regular observations are generated from the model Y_i = X_i+ɛ_i, i = 8, . . . , 25 where X_i ~ U(0, 15) and ɛ_i ~ N(0, 0.5²). Seven outliers, case 1–7 are perturbed at x_i = −4, x_i = −5 + 0.2(i − 1) and X_i = 20 − 0.1(i − 1) to have Y_i = X_i + 4 or Y_i = X_i − 6 − 0.1(i − 1) respectively. Outlyingness of 7 cases are obvious as we see in Figure 1. Columns x₂ and y₂ are obtained by modifying 7 outliers closer to the regular trend.

When columns x₁ and y₁ are used S1, S2, and S3 correctly detected cases 1 to 7 as outliers but M1 identified only case 1 as an outlier. M1 fails because the basic clean subset contains none of true outliers. If columns x₁ and y₂ are used M1 and S1 detected only case 1 as an outlier, whereas S2 and S3 succeeded to detect 7 cases as outliers correctly. Table 2 and Table 3 contain subsets of outliers detected from the data of columns x₂ and y₂. The first column of the tables represents the size of subsets of outliers, decreasing from half of the data size down the table and the second column is the size of the clean subsets. The last column shows outcomes of the outlier test in which “Y” means significant. M1 and S1 show the same results because they form the same basic clean subset of size 12. They failed to find the true outliers because the basic subset contains none of true outliers. S2 starts with a wrongly identified basic clean subset and chooses subsets of potential outliers successfully containing all of the true outliers in the early stages. But it eventually fails to detect the true outliers because swamping occurred at s = 18. By applying S3 to the data, S2 is executed and M1 is also executed with the basic clean subset M₁₄ where M₁₄ = {1 2 3 4 5 6 7 12 15 23 25} because $M_{14}^c$ is quite different from $M_{13}^c=\left\{\hspace{0.17em}1\hspace{0.17em}8\hspace{0.17em}9\hspace{0.17em}10\hspace{0.17em}13\hspace{0.17em}15\hspace{0.17em}16\hspace{0.17em}17\hspace{0.17em}21\hspace{0.17em}22\hspace{0.17em}24\hspace{0.17em}25\right\}$. As seen in Table 3 M1 with a basic clean subset of M₁₄ detects {1 2 3 4 5 6 7} as outliers; consequently, the final detection of outliers using S3 is {1 2 3 4 5 6 7}, the union of {1 2 3 4 5 6 7} and {1}.

A Monte Carlo experiment for comparing the power of M1, S1, S2, and S3 are performed. We use the same outlier patterns and simulation scheme as used in Hadi and Simonoff (1993). The simulation results are limited because the pattern of the predictor, the position of outliers and the sample size are fixed. The data sets are generated from the model Y_i = X_i + ɛ_i where X_i ~ U(0, 15) and ɛ_i ~ N(0, 1²) except two models H3(s) and H3h2(s) which use ɛ_i ~ N(0, 0.7²). Outliers are planted at X_i = 20 − 0.05(i − 1) for high leverage cases and at X_i = −7.5 + 0.05(i − 1) for low leverage cases. The y values are given by Y_i = X_i + 4. Table 4 shows the simulation results.

The first column of Table 4 indicates models. H and L refer to high and low leverage outlier respectively, a number indicates the number of outliers and “s” indicates a model using ɛ_i ~ N(0, 0.7²). Notations used in the second column are defined as: p₁ is the probability of identifying true outliers exactly. p₂ is the probability of identifying at least one true outlier. p₃ is the probability of identifying a clean observation as an outlier. The measures related to two negative effects of outliers are calculated as Pr(masking) = 1 − p₂ and Pr(swamping) = 1 − p₃.

The first row of Table 4 shows that all methods control the size of Type I error at 0.05 approximately. S1, S2, and S3 provided better overall performance than M1.

S3 gave the best results among the examined methods. S2 and S3 had a similar performance. The comparison of H3 and H3(s) shows that the difference in performance among M1, S1, S2 and S3 becomes larger if the ordinary observations are generated from the model with a small variance. The model H3h2(s) used a data set containing a group of 3 high leverage outliers and an additional group of 2 outliers positioned at lower right corner. The y values of two outliers are given by Y_i = X_i − 4 for x_i = −5 + 0.05(i − 1). The results of the model H3h2(s) show that the general patterns do not change when there is more than one group of outliers.

Hadi and Simonoff’s procedure and the suggested methods have low breakdowns as n and p increase. But as the number of outliers is bigger and n is smaller, we can expect that the effect of LQS-residuals increases compared to LS-residuals in selecting a clean set. The proposed methods are applied to data sets having multiple predictors. We summarize the results.

Modified wood gravity data set (Rousseeuw and Leroy, 1987) contains 20 observations of five predictors including four outliers. M1, S1, S2 and S3 succeed in identifying true outliers. Hawkins, Bradu and Kass data set (Hawkins et al., 1984) consists of 75 observations in three predictors. The first ten observations are high leverage outliers. Hadi and Simonoff reported that their method succeeded in detecting the cases 1–10 as outliers. But if the LMS approach is used in constructing a basic clean set, M1 identified cases 11, 12, 13, 14 as outliers. S1 also gave the same results whereas S2 and S3 correctly detected the outliers.

In this article we have proposed three procedures for the identification of the outliers using k-LQS estimates. All three methods are related to the method proposed by Hadi and Simonoff (1993). The first method finds the basic clean subset using LQS estimates. The second method adopts the k-LQS approach for forming a clean subset of size k for the detection of outliers of size (n − k). The third method is designed to be more resistant to swamping or masking problems by combining the first method and the second method. The procedures search a fixed number of outliers using k-LQS estimates; however, they do not require presetting the size of outliers. The size of outliers are determined by t-tests for the no-outlier hypothesis.

As compared with other methods the LQS approach needs more computation. But in large regression problems a subsampling algorithm (Rousseeuw and Leroy, 1987) for approximating LQS estimates can be used to avoid excessive computation.