Automatic detection of discordant outliers via the Ueda's method

Marmolejo-Ramos et al. Journal of Statistical Distributions and Applications (2015) 2:8 DOI 10.1186/s40488-015-0031-y RESEARCH Open Access Automatic detection of discordant outliers via the Ueda’s method Fernando Marmolejo-Ramos1†*, Jorge I. Vélez2,3† and Xavier Romão4 *Correspondence: fernando.marmolejo.ramos@ pshychology.su.se †Equal contributors 1Gösta Ekman Laboratory, Department of Psychology, Stockholm University, Stockholm, Sweden Full list of author information is available at the end of the article Abstract The importance of identifying outliers in a data set is well known. Although various outlier detection methods have been proposed in order to enable reliable inferences regarding a data set, a simple but less known method has been proposed by Ueda (1996/2009). Since this new method, called Ueda’s method, has not been systematically analysed in previous research, a simulation study addressing its performance and robustness is presented. Although the method was derived assuming that the underlying data is normally distributed, its performance was analysed using data from various outlier-prone distributions commonly found in several research fields. The results obtained enable us to define the strengths and weaknesses of the method along with its limits of applicability. Furthermore, an unforeseen field of application of the method, which requires further studies was also identified. Keywords: Ueda’s method; Outliers; Skewed distributions; AIC; Normal distribution; Statistical simulation Classification codes: 97K80; 68U20 1 Introduction Identifying outliers is essential to the data analyst in order to make reliable inferences on the data at hand. Various methods have been proposed for the identification of outliers and the logic of those methods depends directly on how outliers are defined (Aggarwal 2013; Chandola et al. 2009). In cognitive science, specifically in experimental psychology and neuroscience, outliers are accommodated via data transformations (e.g., Box-Cox transformation) or eliminated via truncation (e.g. observations below and/or above the extremes of a pre-set range of values are removed). An approach frequently used is the z-value test for outliers (Songwon, 2006) in which observations are converted to z-scores in order to see which observations fall a pre-set number of standard deviations (SDs) away from the mean. Although there have been proposed SD values to be used given specific sample sizes (see Van Selst and Jolicoeur 1994), in practice researchers use ±2, ±2.5 and ±3SDs as benchmarks regardless (examples of the usage of these benchmarks can be found in Bertels et al. 2010; Havas et al. 2007; Otte et al. 2011; a comprehensive simulation study comparing these and other methods can be found in Marmolejo-Ramos et al.2015). The z-score approach implicitly assumes that the data comes from a normal distribution, thus forcing the actual distribution of the data set to adopt a bell shape. © 2015 Marmolejo-Ramos et al. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Marmolejo-Ramos et al. Journal of Statistical Distributions and Applications (2015) 2:8 Page 2 of 14 Many outlier detection methods are available for univariate data sets (e.g. see Hayes and Kinsella 2003; Thode 2002; Verma 1997; Verma et al. 2014). Among these, the method proposed by Ueda (Ueda 1996/2009) which still assumes the underlying data to be normally distributed and based on the Akaike’s Informatin Criterion (AIC) (Kitagawa 1979), presents interesting features. However, this method has not been systematically studied and, to be best of the authors’ knowledge, it has only been featured in just a handful of studies involving genetic data (Kadota et al. 2003a,b; Tsuyuzaki et al. 2013). The aim of this paper is to study the performance and robustness of the Ueda’s method to detect outliers via computer simulations, as well as determine its applicability to other types of data. It should be noted that, in the context of the proposed study, outliers are not just seen as errors in the data, as considered in other situations (see Hoaglin et al. 1983 for an example). Instead, outliers are generally viewed as data values which are numerically distant from the bulk of the sample, thus requiring particular consideration given their significance. In the first section, the Ueda’s method is described and a modification of it, aimed at automating the selection of the number of candidate outliers, is proposed. In addition, some examples are provided to illustrate this modification. In order to study the performance of the Ueda’s method, various outlier-prone distributions commonly found in actual research are used. For comparative purposes, a 95 % confidence level is used (i.e., the nominal type I error probability is 5 %) and the normal distribution is also included in the study. Finally, and based on the results obtained, potential applications of and research topics involving the Ueda’s method are discussed in the context of applied and statistical research. 1.1 Background Let Xbe a random variable with probability distribution fX(x),andlet x1,x2,...,xn,xn+1,...,xn+s−1,xNbe a random sample of size N=n+sfrom this distribution, with nand s={1, 2, 3, ...}being the number of regular and outlying observations, respectively. The Ueda’s method is a simplified version of the use of the Akaike Information Criterion (AIC) proposed by Kitagawa (Kitagawa 1979). The aim of the AIC, estimated as AIC =2k−2ln(L),(1) is to provide a measure for model selection. In the expression above, kis the number of independently adjusted parameters and Lis the maximum likelihood value for nnumber of regular observations. Thus, the best probabilistic model for the regular observations is that maximising L= n  i=1 fX(xi),(2) where nrepresents the regular observations in a given univariate data set and fX(xi)represents the value of the density function for the ith observation (i=1, 2, ...,n).Inthe normal distribution, the AIC becomes: AIC =2nln ˆσ−2lnn!+2s(3) Marmolejo-Ramos et al. Journal of Statistical Distributions and Applications (2015) 2:8 Page 3 of 14 where ˆσ=N−1/2    N  i=1 (xi−¯ x)2(4) is the standard deviation of the full sample, and ¯ xis the sample mean. Ueda (Ueda 1996/2009) noticed that n−1ln n!≈1forn∼5−9, ≈2forn∼10 −28, ≈3forn∼29 −82, and so forth (see Fig. 1a). By considering the term 2s−2lnn!inthe right side of Eq. (3) and as n=N−sand n−1ln n!≈1forn∼5−9, it follows that −2lnn!+2s≈2s−2n(5) =2s−2(N−s) =4s−2N =4s−constant This same procedure can be used for other values of nin order to obtain 6s−constant, 8s−constant, and so on (see Table two in Ueda 1996/2009). Note that the resulting values 4s,6s,8s,...are correction terms that depend on n,butthatareinaccuratesincethetotal number of samples Nis constant. To ameliorate the problem, Ueda used the correction factor √2sln n! nwhich still depends on n, but includes the number of outliers sin the full sample (see Fig. 1b). Based on this, the proposed test statistic to identify outliers using this method is given by (Ueda 1996/2009): Ut=1 2AIC (6) ≈nln ˆσ−√2sln n! n =(N−s)ln ˆσ−√2sln n! n· In Eq. (6), the term nln ˆσprovides an index of the predictability of true outliers from a set of candidate outliers (i.e., particularly from both ends of the tails of a distribution), whereas the term √2sln n! npenalises the model according to the number of parameters (i.e., the number of outliers s). Specifically, low values in Fig. 1 Graphical representation of the results reported originally by Ueda (Ueda 1996/2009) in Tables one (a) and three (b). Here, N=n+sis the sample size, nis the number of non outlying observations, and s represents the number of candidate outliers Marmolejo-Ramos et al. Journal of Statistical Distributions and Applications (2015) 2:8 Page 4 of 14 the first term indicate that true outliers have been found, whereas high values in the second term signal unreliability of the model due to having many parameters (Kadota et al. 2003a,b). Thus, the lowest AIC associated with a combination of candidate outliers signals the best model that detects a set of potential outliers. 1.2 Calculating Ut Based on Ueda (1996/2009), the following procedure is suggested to determine outliers in a full sample: 1. Order the full sample to obtain xord =x(1),x(2),...,x(N),with x(i)being the i th order statistic (i=1, 2, ...,N). 2. Calculate the z -score for each x(i)as z(i)=x(i)−¯ x ˆς(7) with ˆς=N i=1(x(i)−¯ x)2 N−1 and ¯ xthe sample mean (i=1, 2, ...,N). 3. Provided a number of steps s≥0, calculate the test statistic Utas in Eq. (6). Observe that s=0means that no outliers are present in the data, and any other value of s implies otherwise. Furthermore, as the original sample is ordered, s>0 also implies that the outliers are at either end of the sample (e.g. on the lower or upper tail of the distribution); which end of the ordered sample is tackled first, heavily depends on how the calculation of the test statistic is performed (see §1.3). 4. Suppose that s=1. It could be the case that the outlier is located either at the beginning of the ordered sample, i.e. x(1)istheoutlier,orattheendofit,i.e.x(N)is the outlier. If the former, then the test statistic Utis calculated as in Eq. (6) after removing x(1)from the ordered sample. If the latter, x(N)is removed and Ut subsequently calculated. 5. Step 4 continues for other values of s>1. As a result of this process, a collection of Utvalues is obtained. 1.3 On the number of steps The previous section described how the test statistic Utis calculated provided a value of s≥0. One question that remains to be answered is how many values of Utneed to be calculated. Let Cbe such a number and suppose that up to soutliers are to be detected in the sample. Then, from the collection of Utvalues obtained in step 5 of §1.2, the following matrix is constructed: Ut= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ U00 U01 U02 ···U0m U10 U11 U12 ···U1m U20 U21 U22 ···U2m . . .. . .. . ..... . . Um0Um1Um2···Umm ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦m×m (8) where m=s+1, U00 is the Utstatistic when no outliers are detected, and Uij is the Utstatistic when iand joutliers are detected in the lower and upper tails of xord, Marmolejo-Ramos et al. Journal of Statistical Distributions and Applications (2015) 2:8 Page 5 of 14 respectively (i,j=1, 2, ...,m). This is equivalent to removing x(i),x(i−1),...,x(1)and x(N−j),x(N−j+1),...,x(N)from xord. From Utit is clear that the number of calculations to detect up to soutliers in xord is, at the most, m2. However, this number can be reduced to m2−2ifU00 and Umm, included for convenience and comparison purposes, are not calculated. Hence, m2−2≤C≤m2. To illustrate this, consider the squared upper-left block of Utwhose entries are U00, U10,U01 and U11, and which correspond to the test statistic in Eq. (6) when no outliers are detected and when x(1),x(N),andbothx(1)and x(N)are removed from the full sample, respectively. 1.4 Automatic detection of the number of discordant outliers In the original version of the Ueda’s method, the analyst sets the number of expected outliers in the data set in advance, i.e. the value of sfollowing the notation considered herein. Carling (2000) showed that, for several probability distributions and various sample sizes, approximately 5 % to 50 % of the total number of observations Ncan be, in general, set as the minimum and maximum reasonably expected number of potential outliers, respectively. In the automation process for selecting the number of candidates outliers sin a sample of size N,sis constrained to be 1 ≤s≤smax,withsmax being the maximum number of outliers to be detected. In this inequality, the minimum value of sis set to 1 such that at least one outlier is detected by the Ueda’s method in a full sample. Because of how the Utstatistic is calculated (see §1.2), it is straightforward to show that smax =⎧ ⎪ ⎨ ⎪ ⎩ 1 2(N−1)if nis odd, N 2−1ifnis even. (9) Based on the value of smax,theUtmatrix is then constructed as in §1.3. It can be seen that the number of outliers in the sample scorresponds to the sum of the indexes of the entry of Utfor which the Utstatistic is minimum, and the location of these outliers in xord is given by the entry’s indexes. For instance, if the minimum value of Utis U12,then the number of outliers in the sample is s=1+2=3 and the values will be given by x(1), x(N−1)and x(N)in xord. 1.5 Examples This section presents three examples illustrating the use of the Ueda’s method and the proposed implementation to automatically detect candidate outliers. Here, one normal, one positivelyand one negatively-skewed distributions are shown. In order to have control on the outliers to be found via the Ueda’s method, some observations were added to each one of these distributions to play the role of outliers. In the three scenarios, the outlying observations were sampled from normal distributions and placed on the right tails of the normal and positively-skewed distributions, and on the left tail of the negatively-skewed distribution. Example 1. Normally distributed data. Consider a random sample x1,x2,...,x25 from aN(300, 102)distribution, and a second random sample y1,y2,...,y5of outlier observations from a N(400, 52)distribution. Then, the full sample of size N=n+s=30 will Marmolejo-Ramos et al. Journal of Statistical Distributions and Applications (2015) 2:8 Page 6 of 14 be given by x1,x2,...,x25,y1,y2,...,y5. Although the true number of outliers is s=5, smax issetto10(seeAppendix6.1)soupto10observationsareautomaticallytestedfor their potentiality to be outliers using the Ueda’s method. As shown in Fig. 2a, the number of outliers detected in the full sample is s=5. Example 2. Positively-skewed distributed data. One hundred observations were generated from an Ex-Gaussian distribution with parameters μ=100, σ=10 and ν=10 and 20 outlying observations were generated from a N(200, 202).UsingtheUeda’smethod with smax =25, a total of 22 outliers are detected. Fig. 2b depicts the Utstatistic as a function of smax. Note that the minimum value of Utis reached when smax =22, which suggests that the number of outliers in the full sample of size N=120 is s=22 and not s=20 as initially introduced. This example shows that the Ueda’s method detects those outlier observations initially introduced in addition to a few more that were not thought to be outliers in the full sample. Closer inspection reveals that these two additional observations are located in the upper tail of the distribution and may, as the Ueda’s method shows, be signalled as potential outliers (see Fig. 2b). Example 3. Negatively-skewed distributed data. Consider a random sample of size n=50 from a Beta distribution with parameters α=5andβ=2. In this case, three observations generated from a N(0.1, 0.022)distribution are added for the Ueda’s method to detect them in the combined sample of size N=53 using smax =5. As shown in Fig. 2c, the number of outliers detected is the one introduced. In the case of the normal distribution, the outlying observations were readily captured by the Ueda’s method. However, examples 2 and 3 represent cases of distributions prone to having outliers (in the sense they were defined in Section 1). That is, these types of Fig. 2 Results of Ueda’s method for examples a1, b2andc3. First row: Probability densities for the non contaminated (green), the contaminating (pink) and full (light blue) samples. Second row: Utstatistic as a function of smax. The value of smax for which Utis the lowest (red dot) indicates the number of candidate discordant outliers in the full sample Marmolejo-Ramos et al. Journal of Statistical Distributions and Applications (2015) 2:8 Page 7 of 14 distributions can have outliers as they have elongated tails (see Gleason 1993) accounting for some observations that fall away from where most of the data tend to cluster. Although in the case of the negatively-skewed distribution only the observations placed on the left tail were appropriately caught by the method, this did not happen in the case of the positively-skewed distribution. In that case, two observations that were not part of the added normally-distributed observations were signalled as outliers. This specific case demonstrates that skewed distributions can in fact have outlying observations on the tails and that such observations do occur due to the natural way in which such types of distributions can take place (e.g., positively-skewed data is the norm in human reaction time research). What is interesting is that the Ueda’s method seems to be sensitive to these types of naturally-occurring outlying observations. As skewed distributions are prone to having outliers and are typically found in actual research, it is thus essential to determine how the Ueda’s method performs in these cases. 2 Simulation study To study the performance of the Ueda’s method, a simulation study was carried out using a similar simulation scheme to that presented by Vélez and Correa (2014), and following the recommendations by Salazar and Baena (2009). In brief, nobservations of a particular probability distribution (see Table 1) are generated and a measure of interest is calculated. The following algorithm was implemented in R (R Core Team 2013) to evaluate the performance and robustness of the Ueda’s method when no outliers are introduced into the data: 1. Draw n random observations from a probability distribution fX(x|θ), where θis the parameter vector (see Table 1). 2. Apply the Ueda’s method and determine the number of outliers detected. Denote this number as ˆs. 3. Repeat steps 1 and 2, B times, and calculate R=N0/B,whereN0is the number of times in the B samples of size n that ˆs=0, i.e. no outliers are detected. After generating n={10, 20, 30, 50, 100, 200, 500}observations from each probability distribution in Table 1, Rwas calculated. These distributions and sample sizes were chosen because of what is often seen in real-world applications. Indeed, the Ex-Gaussian, Gamma, Weibull and Lognormal distributions are positively skewed distributions used Table 1 Probability distributions considered in this study. The probability density function f(·)is shown in the second column, and the parameter vector θdefining each distribution in the third column Distribution f(·)θ Normal 1 √2πσ e−(x−μ)2 2σ2μ,σ Student’s t( ν+1 2) (ν/2)√νπ 1+x2 ν−ν+1 2ν Tukey λSee text in §2 λ Beta (α+β) (α)(β) xα−1(1−x)β−1α,β Ex-Gaussian 1 ν√2πe σ2 2ν2−x−μ ν·[(x−μ)/σ]−σ/ν −∞ e−y2 2dy μ,σ,ν Gamma 1 (α)βαxα−1e−x βα,β Weibull β αβxβ−1e−(x α)βα,β Lognormal 1 x√2πσ e−(ln x−μ)2 2σ2μ,σ Marmolejo-Ramos et al. Journal of Statistical Distributions and Applications (2015) 2:8 Page 8 of 14 in cognitive science to model behavioural and neurological data (e.g., Leiva et al. 2015). The Beta distribution has been used to model soil data (Haskett et al. 1995) and rates and proportions (Ferrari and Cribari-Neto 2004), the Student’s thas been used to fit share price changes (Praetz 1972), and the Tukey-λhas been fitted to solar radiation data (Öztürk and Dale 1982). In all simulation scenarios, a total of B=10, 000 replicates were used. In what follows, it is described how the simulation scenarios were constructed for each probability distribution f. Every single probability distribution fin Table 1 is defined by a set of parameters θ. Provided f, a simulation scenario is defined as a combination of nand specific values of θ for that f. For instance, in the Normal distribution, θ=(μ,σ). Without loss of generality, μwas fixed at 0 and σ={0.5, 1, 1.5, ...,5}for each sample size, to obtain a total of 70 simulation scenarios. In the Student’s t-distribution, 350 scenarios were studied and each of them was defined by the combination of the sample size nand ν={1, 2, ...,50}. Let p∈[0, 1] and x=F−1(p),withFbeing the cumulative distribution function of a random variable X.IntheTukey-λdistribution, the quantile and density function are respectively given by (Chalabi et al. 2014) F−1(p|λ) =λ+pλ−(1−p)λ λ(10) and f(x|λ) =fF−1(p)=1 pλ−1+(1−p)λ−1·(11) Forspecificvaluesofλ,theTukey-λdistribution resembles the characteristics of some known probability distributions. For instance, the values λ=−1, λ=0.14 and λ=1 correspond to the Cauchy(0, π),N(0, 1.46362)and Uniform(−1, 1)distributions, respectively. In the present simulation strategy, λwas varied in the interval [-2, 2] in steps of 0.2, excluding λ=0 and including λ=0.14. Thus, a total of 147 scenarios were evaluated. For the Beta(α,β) distribution, the parameters αand βvaried in the rectangle [a,b]×[a,b]witha=0.5 and b=5 using increments of h=0.5 within each margin. The Ex-Gaussian distribution is defined by θ=(μ,σ,ν). To generate observations from this distribution, μ={200, 300, 400, 500, 600},σwas fixed at 20, and ν= {400, 300, 200, 100, 50}. In the Gamma(α,β) and Weibull(α,β) distributions, α,β={2, 4, 6, 8, 10}.Asimilar approach was used for the Lognormal distribution after replacing μby αand σby β. 3 Results As shown in Fig. 3a, the results for the normal distribution indicate that the proportion of no outliers detected in the data by Ueda’s method is below the expected 95 % when n<100, and that the standard deviation σseemstohavenoeffectonthis.Evenso,the proportion of times the method finds no outliers is >90 % across sample sizes. Conversely, the method tends to find more outliers than expected when the degrees of freedom of the Student’s t-distribution are small (ν<10, Fig. 3b). This is particularly clear when both the degrees of freedom and the sample size are small (in which case the proportion of not finding outliers decreases to ∼0.75). Graphical inspections of Student’s t-distributions corroborate this finding: when ν=1, the Student’s t-distributions presents some observations that fall very far away from the mean, whereas for ν=10 the observations Marmolejo-Ramos et al. Journal of Statistical Distributions and Applications (2015) 2:8 Page 9 of 14 Fig. 3 Results for the aNormal, bStudent’s t,andcTukey-λdistributions. Higher values of Rare presented in dark red, and lower values in dark blue. The parameter ncorresponds to the full sample size spread out more evenly around the mean. In the Tukey-λdistribution, more outliers tend to be detected when the λshape parameter decreases (λ<0, Fig. 3c). That the Ueda’smethodfindsmoreoutliersinthisdistributionisaccentuatedwhenbothnand λ are small. It is worth noting that when λis large and positive, an uniform distribution is obtained and the Ueda’s method does not find outliers. Although uniform distributions are not normal, they are symmetric and observations are spread evenly across the entire data range; these specific results thus suggest that the Ueda’s method focuses on the shape and symmetry of the distribution. Overall, compared to the normal distribution case, the Ueda’s method tends to find more outliers in the Student’s tand Tukey-λdistributions, especially when their parameters generate distributions that are too platykurtic or too leptokurtic. Fig. 4 Results for the Ex-Gaussian, Beta, Gamma, Weibull and Lognormal distributions for fixed sample sizes. Higher values of Rare presented in dark red, and lower values in dark blue/gray