Studies in Mathematical Sciences Vol. 4, No. 2, 2012, pp. 1-9 DOI: 10.3968/j.sms.1923845220120402.1623 ISSN 1923-8444 [Print] ISSN 1923-8452 [Online] www.cscanada.net www.cscanada.org A Procedure for Detecting a Pair of Outliers in Multivariate Dataset B.K. Nkansah[a],∗; B.K. Gordor[a] [a]Department of Mathematics and Statistics University of Cape Coast, Ghana ∗Corresponding author. Address: Department of Mathematics and Statistics University of Cape Coast, Ghana Received 22 January, 2012; accepted 19 May, 2012 Abstract The paper presents a procedure for detecting a pair of outliers in multivariate data. The procedure involves a reduction of the dimensionality of the dataset to only two dimensions along outlier displaying compo- nents, and then determines the orientation of a least squares ellipse that fits the scatter of points of the two dimensional dataset. Finally, the reduced data is projected unto a vector which is determined in terms of the orientation of the ellipse. The results show that if two observations constitute a pair of outliers in a data set, then the pair is extreme at either ends of the one-dimensional projection and separated clearly from the remaining observations. If the two outliers are not distinct on such a one-dimensional projection, three key rules are prescribed for successful determination of the right pair of outliers. Key words Multiple Outlier Detection; Outlier Displaying Component B.K. Nkansah, B.K. Gordor (2012). A Procedure for Detecting a Pair of Outliers in Multivariate Dataset. Studies in Mathematical Sciences, 4(2), 1-9. Available from: URL: http://www.cscanada.net/index. php/sms/article/view/j.sms.1923845220120402.1623 DOI: http://dx.doi.org/10.3968/j.sms.1923845220120402.1623 1. INTRODUCTION The detection of a single outlier in a multivariate dataset is one on which several methods converge. One main approach is the use of the Mahalanobis generalised squared distance U(xε ,Sε) = (xε − x̄)′S−1(xε − x̄) (1) On the detection of multiple outliers, several procedures have been developed. Most of these methods (e.g. Rousseeuw, 1985; Rousseeuw & Leroy, 1987) rely on the generalised distance in Equation (1) which is often used as a test statistic. Others (e.g. Caroni & Prescott, 1992) make use of the Wilk’s ratio statistic which determines the subset Tk of k observations among the sample of size n for which the ratio RT = | S(Tk) | | S | is minimum, where S is the matrix of sum of squares and cross-product and S(Tk) is the corresponding matrix with the observations in Tk removed from the sample. The results of some of these methods are divergent particularly on the problem of multiple outlier detection. A major cause of this divergence lies in the different approaches adopted at controlling the use of the general sample mean in the detection procedure. Another approach is the one by Gordor and Fieller (1999). Their approach, which is a graphical technique and referred to as the Outlier Displaying Component (ODC), is actually a displaying technique rather than a detecting one. This means that in the case of two outliers, for example, the outliers must be specified before the technique can be used to display them. However, very often these multiple outliers are not known. 1 B.K. Nkansah; B.K. Gordor/Studies in Mathematical Sciences Vol.4 No.2, 2012 We can identify some two basic difficulties that are common to these and many other attempts at multiple outlier detection. These are the amount of computations involved and the subjectivity that characterize these methods. In the next section we will discuss a procedure for detecting a pair of outliers from a one-dimensional plot of some univariate equivalent of the multivariate dataset. The process, which is a graphical approach, combines the technique of outlier displaying component and minimum ellipse fitting. The technique is illustrated by using some artificial data (see Appendix) which are generated from N(0, I) and consist of 50 four-dimensional observations. The reason for the choice of this artificial dataset will become apparent. In the end, we compare observations that constitute a pair of outliers in some well-known datasets by means of the proposed method and those obtained by the Wilk’s ratio statistic. 2. DESCRIPTION OF THE PROCEDURE The procedure is in three major phases. The first phase broadly deals with dimensionality reduction of the data set from p-dimensional (p ≥ 3) to 2-dimensional. The second deals with the determination of the orientation of a least squares ellipse that fits the scatter of points of the two dimensional dataset. The last phase then considers a projection vector in terms of the orientation of the ellipse obtained in the second phase. The three phases are expanded in this section. 2.1 Dimensionality Reduction of the Dataset Suppose that X = (x1, x2, · · · , xn) is a p-dimensional random sample with sample mean x̄. Let the single outlier in the sample be xε . The One-Outlier Displaying Component (ODC) vector [4] is given by β = S−1(xε − x̄) (2) where S is p × p sample sum of squares and cross-product matrix. Now, define the vector 1 = (1, 1, · · · , 1)′n×1. Then a projection of the mean-corrected data on β gives a univariate equivalent of the data given by t1 = (X − (1 × x̄′)) × β (3) Now, β is actually an eigenvector of the only non-zero eigenvalue of the matrix E = S−1(xε − x̄)(xε − x̄)′. This non-zero eigenvalue represents the squared generalised distance in Equation (1) of xε from the sample mean, x̄. Now, let Vp×p = (v1, v2, · · · , vp) be a matrix whose columns are eigenvectors corresponding to the eigenvalues of E, with the first column being β. We apply the Gram-Schmidt Orthog- onalization to convert the independent vectors of V into a set of orthogonal vectors, wk; (k = 1, 2, · · · , p) as follows: First, define w1 = v1. Then each wk is made orthogonal to the preceding w1,w2, · · · ,wk−1 by the relation wk = vk − k−1∑ j=1 v′kw j ‖w j‖2 w j, k = 2, 3, · · · , p (4) Since v1, v2, · · · , vp are linearly independent, wk , 0; k = 1, 2, · · · , p and v′kw j = 0, k , j Let uk = 1 ‖wk‖wk, k = 1, 2, · · · , p Thus, U = (u1,u2, · · · ,up) are orthonomal. Now the first vector u1 of U completely contains all information on the single outlier, xε . It will therefore not be useful to use u1 in the determination of a second component that contains information about the second outlier. Hence, we create a matrix U\u1 = (u2,u3, · · · ,up)p×p−1. After the orthogonalization, an (n − 1) × p transformed data set, Y(ε) = X \ xε × U \ u1 is obtained which excludes that of xε . Suppose ȳ is the mean vector and yη, the single outlier in Y(ε) obtained the usual way. Using ȳ and yη, the displaying component for Y is found as ϕ = S−1 Y (yη− ȳ), where SY is the sum of squares 2 B.K. Nkansah; B.K. Gordor/Studies in Mathematical Sciences Vol.4 No.2, 2012 and cross-product matrix of Y(ε). The component ϕ of dimension p − 1, is referred to as the Sub-Outlier Displaying Component (Sub-ODC) [4]. Now replace the ε th row of 1 with zero and then augment the n− 1× p data set Y(ε) to n× p by replacing the ε th row by a vector of zeros. The modified matrices 1 and Y(ε) are represented by 1aug and Y(ε)aug, respectively. A projection of the mean-corrected data Y(ε) on the Sub-ODC is taken to obtain a univariate data t2 = (Y(ε)aug − (1aug × ȳ′)) × ϕ (5) The data set Tn×2 = [t1 t2] gives a bivariate equivalent of X. Thus, Tn×2 provides a dimensionally reduced data set derived from the original data X. A scatter plot of Tn×2 may be generated in the ODC-Sub-ODC plane. Illustration In Figure 2.1 we have the scatter plot of the bivariate equivalent of the Artificial data set described in the previous section. It can be seen that observation 34 is extreme along the 1-ODC dimension (and is thus the single outlier). When observation 34 is deleted, the observation with the greatest generalized distance in the reduced data (Y(34)) is observation 45 which is extreme on the Sub-ODC. Figure 2.1 shows that observations 9 and 22 could also be candidates for an outlying pair. Figure 2.1 Scatter-plot of Artificial Data in ODC-Sub-ODC Plane After describing the next two phases of the method, it will then be possible to identify the two most extreme pair of observations in this dataset. 2.2 The Least Squares Ellipse and its Orientation This phase involves two main steps: (1) The determination of a least squares ellipse that fits the scatter of points, and (2) The determination of the orientation of the ellipse and its centre. 2.2.1 Determination of Least Squares Ellipse An ellipse is a special case of a general conic which can be described by an implicit second order polynomial E(x, y) = ( x y 1 )  a b 2 d 2 b 2 c e 2 d 2 e 2 f   x y 1  = 0 (6) with an ellipse-specific constraint b2 − 4ac < 0. (7) 3 B.K. Nkansah; B.K. Gordor/Studies in Mathematical Sciences Vol.4 No.2, 2012 The coefficients of the ellipse a, b, c, d, e, f are the parameters of ellipse and (x, y) are coordinates of points on it. The polynomial E(x, y) is the algebraic distance of the point (x, y) to the given conic. To fit a general ellipse to a set of data points (xi, yi), i = 1, 2, · · · , n, is to minimize the sum of squared algebraic distances of (xi, yi) to the conic subject to the constraint in Equation (7). Under a proper scaling, the inequality constraint, Equation (7), can be changed into an equality constraint 4ac − b2 = 1. Let a′ = (a, b, c, d, e, f ) andx′ = (x2, xy, y2, x, y, 1). The ellipse-specific fitting problem may be formulated as min a ‖Da‖2 subject to a ′ Ca = 1 (8) where Cn×6 is the constraint matrix and Dn×6 is the design matrix. Now let Se = D′D. By introducing the Lagrange Multiplier λ, we write the formulations in Equation (8) as ψ = a′Sea + λ(1 − a′Ca). Hence, differentiating ψ with respect to a and equating to zero, the conditions for optimal solution of a are obtained as Sea = λCa a ′ Ca = 1 (9) To determine the parameters, partitioning of the matrices in Equation (9) have been suggested in various ways (e.g. Halir & Flusser, 2000; Harker, O’leary & Zsombor-Murray, 2004). Going by the partitioning of [5], we write the design matrix, D, into two n × 3 matrices as follows: D1 =  x2 1 x1y1 y2 1 x2 2 x2y2 y2 2 ... ... ... x2 i xiyi y2 i ... ... ... x2 n xnyn y2 n  , D2 =  x1 y1 1 x2 y2 1 ... ... ... xi yi 1 ... ... ... xn yn 1  (10) The scatter matrix D′D is then given as a partitioned matrix as D ′ D = ( D′ 1D1 D′ 1D2 D′ 1D2 D′ 2D2 ) For easy representation of the element matrices of D′D, which is represented bySe, let D′ 1D1 = S11, D′ 1D2 = S12, and D′ 2D2 = S22. Similarly, the constraint matrix C (Equation (9)) can be expressed as C = ( C1 0 0 0 ) , where C1 =  0 0 2 0 −1 0 2 0 0  (11) and 0 is a 3 × 3 zero matrix. The vector of coefficients a is split accordingly as a = ( a1 a2 ) , where a1 =  a b c  and a2 =  d e f  (12) Using these partitions, the first of the conditions in Equation (9) becomes ( S11 S12 S′21 S22 ) ( a1 a2 ) = λ ( C1 0 0 0 ) ( a1 a2 ) 4 B.K. Nkansah; B.K. Gordor/Studies in Mathematical Sciences Vol.4 No.2, 2012 Hence, we obtain the equations S11a1 + S12a2 = λC1a1 S ′ 21a1 + S22a2 = 0 (13) Before proceeding to solve the system of equations, there is a condition to examine for the existence of real solution of the system. The matrix S22 =  Σx2 i Σxiyi Σxi Σxiyi Σy2 i Σyi Σxi Σyi n  is exactly a scatter matrix of fitting a line through a set of data points. It is known (Haralick & Shapiro, 1993) that this matrix is singular only if all the points lie on a line. In this study, the points (xi, yi) are generated by projections of p-dimensional dataset onto the ODC and the Sub-ODC. By the properties of the two dimensions already described, the scatter of points in the plane provided by these two dimensions cannot lie on a line. Thus, S22 is regular in this case. Consequently, there is a solution for fitting an ellipse through points in the ODC-Sub-ODC plane. Returning to the solution of a′ = (a1, a2), from the second of Equation (13), a2 = −S−1 22 S ′ 21a1 (14) Substituting this into the first of Equation (13) gives (S11 − S12S−1 22 S ′ 21)a1 = λC1a1 and since C1 is regular, C−1 1 (S11 − S12S−1 22 S ′ 21)a1 = λa1( C−1 1 (S11 − S12S−1 22 S ′ 21) − Iλ ) a1 = 0 It follows therefore that a1 is the eigenvector of the reduced scatter matrix M = C−1 1 (S11 − S12S−1 22 S ′ 21) (15) of size 3 × 3 and λ is the corresponding eigenvalue. The conditions of Equation (9) can now be restated as Ma1 = λa1 a ′ 1C1a1 = 1 a2 = −S−1 22 S ′ 21a1 a = ( a1 a2 ) (16) MATLAB codes have been outlined by [5] for computing the ellipse coefficients. 2.2.2 Derivation of the Orientation and Centre of the Ellipse The coefficient vector a in Equation (16) gives a tilted ellipse. Suppose that a′ = (a11, a12, a13, a21, a22, a23), and let a transformation of points (xi, yi) in the ODC-Sub-ODC plane be given by T : x ′ = cx + sy; y ′ = −sx + cy, where c = cosφ and s = sinφ. If this transformation should return a non-tilted ellipse, then the angle φ must be the required angle of orientation of the ellipse. Substituting this transformation in the equa- tion E = x·a,we obtain the expression En = x·a∗,where the components of a∗ = (a∗11, a ∗ 12, a ∗ 13, a ∗ 21, a ∗ 22, a ∗ 23) 5 B.K. Nkansah; B.K. Gordor/Studies in Mathematical Sciences Vol.4 No.2, 2012 are such that, for example, a∗12 = 2a11cs + (c2 − s2)a12 − 2a13cs. For a non-tilt ellipse, a∗12 = 0. After a few simplifications the orientation of the ellipse in Equation (16) is given as φ = 1 2 tan−1 ( a12 a13 − a11 ) (17) The coefficients of the non-tilt ellipse are then obtained by substituting sinφ and cosφ into the components of a∗. By this substitution, the non-tilt ellipse is given by the coefficients a∗ = (a∗11, a ∗ 13, a ∗ 21, a ∗ 22, a ∗ 23) ′ , where a∗23 = a23. It may be necessary to identify the centre of the main cloud of the data points. To determine this centre, (i.e. of the non-tilt ellipse) we express the equation a∗11x2 + a∗13y2 + a∗21x + a∗22y + a∗23 = 0 by square completion as ( x + a∗21 2a∗11 )2 a∗∗23 a∗11 + ( y + a∗22 2a∗13 )2 a∗∗23 a∗13 = 1 where a∗∗23 = −a23+ a∗221 4a∗11 + a∗222 4a∗13 . In this form, the centre (C) of the non-tilt ellipse is given as ( − a∗21 2a∗11 ,− a∗22 2a∗13 ) . Lastly, the centre of the tilt ellipse is obtained by T−1C. 2.3 Determination of a Projection Vector We first note that in the plane provided by the ODC and the Sub-ODC, the most extreme observation along the ODC dimension, xε , has coordinates xε ( U(xε ,S), 0 ) . Given the orientation φ of the minimum ellipse, a unit vector that is perpendicular to this orientation is v = (sinφ, −cosφ). We choose a point xB ε ( U(xε ,S), k ) , directly below xε , where k ≤ 0 is chosen such that k is equal to the smallest value along the Sub-ODC dimension. Now a vector that passes through xB ε and parallel to v is given by we =  ktanφ + U(xε ,S) −k − U(xε ,S) tanφ  (18) This vector, which is perpendicular to the orientation of the ellipse, will subsequently be referred to as the Perpendicular Elliptical Vector. This is to distinguish it from the Parallel Elliptical Vector given by ue =  k tanφ + U(xε ,S) −k − U(xε ,S)tanφ  (19) which passes through xB ε and parallel to the orientation of the ellipse. 2.4 Some Rules for Outlier Detection The use of the Elliptical vector in Equation (18) (or Equation (19)) involves a number of heuristics to suc- cessfully identify the pair of outliers. The following are three properties of the observations that constitute the pair of outliers after projection on the elliptical vector: 1. An observation that sticks out distinctly at one end of the (perpendicular) projection must necessarily be one of the outlying pair. 2. Two observations that stick out distinctly on either ends of the (perpendicular) projection constitute the pair of outliers. 6 B.K. Nkansah; B.K. Gordor/Studies in Mathematical Sciences Vol.4 No.2, 2012 3. If only one observation sticks out distinctly at one end of the (perpendicular) projection, but the other member of the pair does not stick out clearly on the other end, then the reduced two-dimensional data may be projected on the parallel elliptical vector. Each of the observations that stick out on the two separate projections then combine to form the pair of outliers. 3. ILLUSTRATION OF THE METHOD We illustrate the technique by completing the detection of a pair of outliers in the Artificial data in Section 2. Figure 3.1 is the one-dimensional plot of the projection of the reduced data displayed in Figure 2.1 on the perpendicular elliptical vector. Figure 3.1 Projection of Artificial Data on Elliptical Vector From Figure 3.1, it is observed that observations 9 and 45 are extreme on either ends of the projection. According to the second of our rules outlined above, 9 and 45 must be the pair of outliers. In fact, when we test for two outliers by calculating Wilk’s two-outlier statistic, r2 = minTRT , the pair (9, 45) has the smallest value of r2 (= 0.5656). The identification of these two observations as the pair of outliers shows that the most outlying pair in a dataset does not always include the single outlier. 4. A BACK-UP PROCEDURE FOR HIGHER DIMENSIONAL DATASET To successfully detect the right pair of outliers in data of dimension greater than 3, we further propose a ‘back-up’ procedure that could be implemented after the pair of outliers have been detected in a first round implementation of the proposed method. Now the incidence of an outlying observation is as a result of a marked deviation in the observation along one or more dimensions on which data were generated. Thus, a dimension along which there is smallest variation would not contribute significantly to outlier detection. Therefore, to assess the correctness of a detected pair (xε , xι) as outliers, we delete the dimension along which there is smallest variation, and then implement the procedure again. The pair of outliers that emerge in this second round of implementation are the right pair of outliers. Table 1 shows single outliers and a pair of outliers detected by means of the generalized distance (Gen. Dist.) from the sample mean, the proposed method and the Wilk’s ratio method in some widely studied data 7 B.K. Nkansah; B.K. Gordor/Studies in Mathematical Sciences Vol.4 No.2, 2012 sets of varying sizes, n and dimensionality, p. These data sets include the Milk Transportation Cost Data (Johnson & Wichern, 2002), the three types of Iris datasets (Johnson & Wichern, 2002; Anderson, 2003) and the US Food Price data (Sharma, 1996). Table 1 Performances of Proposed Method and Two Others Single Pair of outliers detected by Data Set Outlier Gen. Dist. Proposed Wilk’s Ratio Artificial data (p = 4, n = 50) 34 34, 45 9, 45 9, 45 Iris Setosa (p = 4, n = 50) 42 42, 44 42, 23 42, 23 Iris Versicolor (p = 4, n = 50) 19 19, 49 19, 49 19, 49 Iris Virginica (p = 4, n = 50) 19 19, 32 19, 18 19, 32 Transport-Cost (p = 3, n = 36) 9 9, 21 9, 21 9, 21 U.S Food Price (p = 5, n = 23) 10 10, 1 10, 16 10, 16 Table 1 shows that with exception of the Iris Virginica, the results of both the proposed method and the Wilk’s ratio method coincide in all the samples used. Also, with exception of the Artificial dataset, all others (which are real datasets) include the single outlier in the pair of outliers. 5. CONCLUSION The paper has considered a procedure for detecting a pair of outliers in multivariate data. We discover that if two observations constitute a pair of outliers, then there must exist a one-dimensional projection (given by the Perpendicular Elliptical Vector) on which the pair is extreme at either ends and separated clearly from the remaining observations. If the two outliers are not distinct on such a one-dimensional projection, a combination of two separate one-dimensional projections (the other given by the Parallel Elliptical Vector) may be obtained. Each of the observations that is extreme on the two separate projections combine to form the pair of outliers. The method proposed aims at reducing the dimensionality of the original data and then removes the influence of the generalised sample mean by eliminating all observations that are close to it. REFERENCES [1] Anderson, T. W. (2003). Introduction to Multivariate Statistical Analysis. New Jersey, Prentice Hall. [2] Caroni, C. & Prescott, P. (1992). Sequential Application of Wilk’s Multivariate Outlier Test. Applied Statistics, 41, 355-364. [3] Gordor, B. K. & Fieller, N. R. J. (1999). How to Display an Outlier in Multivariate Datasets. Journal of Applied Sciences & Technology, 4(2). [4] Halir R. & Flusser, J. (2000). Numerically Stable Direct Least Squares Fitting of Ellipses. Technical Report citeseer.nj.nec.com/350661.html, Department of Software Engineering, Charles University, Czech Republic. [5] Haralick R. M. & Shapiro, L., G. (1993). Computer and Robot Vision, 1. Addison Wesley. [6] Harker, M., O’leary P. & Zsombor-Murray, P. (2004). Direct and Specific Fitting of Conics to Scattered Data. Technical Report, Institute for Automation, University of Leoben. [7] Johnson, R. A. & Wichern, D. W. (2002). Applied Multivariate Statistical Analysis. New Jersey: Pren- tice Hall. 8 B.K. Nkansah; B.K. Gordor/Studies in Mathematical Sciences Vol.4 No.2, 2012 [8] Rousseeuw, P. J. (1985). Multivariate Estimation with High Breakdown Point. Mathematical Statistics and Application, B, 283-297. [9] Rousseeuw, P. J. & Leroy, A. (1987). Robust Regression and Outlier Detection. New York: Wiley. [10] Sharma, S. (1996). Applied Multivariate Techniques. New York: Wiley. Appendix: Artificial Data No. X1 X2 X3 X4 No. X1 X2 X3 X4 1 0.074 0.037 -0.479 0.296 26 -0.775 0.970 -1.420 0.169 2 -1.051 -0.095 0.098 -1.495 27 0.589 1.026 0.331 -1.977 3 1.767 -1.469 -0.592 2.392 28 -0.153 0.064 0.354 -0.306 4 -2.623 0.726 -1.317 0.106 29 1.808 -1.478 0.886 1.677 5 0.281 -0.650 -2.044 0.827 30 1.720 1.263 1.020 -0.417 6 1.717 -0.429 0.457 0.825 31 -1.050 0.147 1.306 2.716 7 -1.012 0.504 -0.514 0.244 32 -1.040 0.152 1.341 0.529 8 1.404 -0.549 0.558 1.473 33 0.058 1.614 0.208 -0.524 9 1.804 -1.849 -1.477 -0.921 34 1.179 1.914 -1.744 1.294 10 -1.423 -0.566 -0.342 0.120 35 -0.096 0.353 1.241 -0.584 11 0.560 -0.188 1.035 0.242 36 -0.564 -0.846 -0.892 0.070 12 0.605 1.324 -0.460 0.194 37 0.356 -0.952 -1.107 -0.354 13 0.480 1.954 -0.254 -0.213 38 -1.570 0.409 -0.141 -1.143 14 -1.166 1.151 -1.395 -0.438 39 0.470 -1.022 0.571 -0.772 15 1.203 1.511 1.219 -1.623 40 -0.480 0.972 -1.112 0.798 16 -1.135 -0.541 -0.392 1.028 41 0.222 -0.687 0.672 1.464 17 -0.027 0.742 -0.784 0.151 42 0.397 -0.003 0.444 0.160 18 0.580 -0.575 -0.234 0.118 43 -0.117 0.175 -0.758 -0.681 19 -1.111 0.952 -0.596 -0.733 44 -1.359 0.535 -2.054 0.335 20 -0.873 -1.631 -0.423 1.627 45 -1.640 1.721 1.726 1.181 21 1.315 -0.775 0.908 0.025 46 0.116 -0.400 0.975 3.444 22 0.957 0.191 -1.607 1.344 47 -1.391 -1.200 -0.305 0.088 23 -0.736 -0.494 -0.177 -0.738 48 0.726 -0.079 -1.076 0.284 24 -2.527 -0.661 -2.139 -0.391 49 1.149 0.126 0.193 -1.171 25 -0.352 -1.381 -1.146 0.801 50 -0.640 1.464 -1.224 -1.618 9