knitr::opts_chunk$set( echo = TRUE, comment = "#>", out.width = "100%", dpi = 300, warning = FALSE ) library(performance) library(see) library(datawizard)
Beyond the challenge of keeping up-to-date with current best practices regarding the diagnosis and treatment of outliers, an additional difficulty arises concerning the mathematical implementation of the recommended methods. Here, we provide an overview of current recommendations and best practices and demonstrate how they can easily and conveniently be implemented in the R statistical computing software, using the {performance} package of the easystats ecosystem. We cover univariate, multivariate, and model-based statistical outlier detection methods, their recommended threshold, standard output, and plotting methods. We conclude by reviewing the different theoretical types of outliers, whether to exclude or winsorize them, and the importance of transparency.
Real-life data often contain observations that can be considered abnormal when compared to the main population. The cause of it can be hard to assess and the boundaries of "abnormal", difficult to define---they may belong to a different distribution (originating from a different generative process) or simply be extreme cases, statistically rare but not impossible.
Nonetheless, the improper handling of these outliers can substantially affect statistical model estimations, biasing effect estimations and weakening the models' predictive performance. It is thus essential to address this problem in a thoughtful manner. Yet, despite the existence of established recommendations and guidelines, many researchers still do not treat outliers in a consistent manner, or do so using inappropriate strategies [@simmons2011false; @leys2013outliers].
One possible reason is that researchers are not aware of the existing recommendations, or do not know how to implement them using their analysis software. In this paper, we show how to follow current best practices for automatic and reproducible statistical outlier detection (SOD) using R and the {performance} package [@ludecke2021performance], which is part of the easystats ecosystem of packages that build an R framework for easy statistical modeling, visualization, and reporting [@easystatspackage]. Installation instructions can be found on GitHub or its website, and its list of dependencies on CRAN.
The instructional materials that follow are aimed at an audience of researchers who want to follow good practices, and are appropriate for advanced undergraduate students, graduate students, professors, or professionals having to deal with the nuances of outlier treatment.
Although many researchers attempt to identify outliers with measures based on the mean (e.g., z scores), those methods are problematic because the mean and standard deviation themselves are not robust to the influence of outliers and those methods also assume normally distributed data (i.e., a Gaussian distribution). Therefore, current guidelines recommend using robust methods to identify outliers, such as those relying on the median as opposed to the mean [@leys2019outliers; @leys2013outliers; @leys2018outliers].
Nonetheless, which exact outlier method to use depends on many factors. In some cases, eye-gauging odd observations can be an appropriate solution, though many researchers will favour algorithmic solutions to detect potential outliers, for example, based on a continuous value expressing the observation stands out from the others.
One of the factors to consider when selecting an algorithmic outlier detection method is the statistical test of interest. Identifying observations the regression model does not fit well can help find information relevant to our specific research context. This approach, known as model-based outliers detection (as outliers are extracted after the statistical model has been fit), can be contrasted with distribution-based outliers detection, which is based on the distance between an observation and the "center" of its population. Various quantification strategies of this distance exist for the latter, both univariate (involving only one variable at a time) or multivariate (involving multiple variables).
When no method is readily available to detect model-based outliers, such as for structural equation modelling (SEM), looking for multivariate outliers may be of relevance. For simple tests (t tests or correlations) that compare values of the same variable, it can be appropriate to check for univariate outliers. However, univariate methods can give false positives since t tests and correlations, ultimately, are also models/multivariable statistics. They are in this sense more limited, but we show them nonetheless for educational purposes.
Importantly, whatever approach researchers choose remains a subjective decision, which usage (and rationale) must be transparently documented and reproducible [@leys2019outliers]. Researchers should commit (ideally in a preregistration) to an outlier treatment method before collecting the data. They should report in the paper their decisions and details of their methods, as well as any deviation from their original plan. These transparency practices can help reduce false positives due to excessive researchers' degrees of freedom (i.e., choice flexibility throughout the analysis). In the following section, we will go through each of the mentioned methods and provide examples on how to implement them with R.
Researchers frequently attempt to identify outliers using measures of deviation from the center of a variable's distribution. One of the most popular such procedure is the z score transformation, which computes the distance in standard deviation (SD) from the mean. However, as mentioned earlier, this popular method is not robust. Therefore, for univariate outliers, it is recommended to use the median along with the Median Absolute Deviation (MAD), which are more robust than the interquartile range or the mean and its standard deviation [@leys2019outliers; @leys2013outliers].
Researchers can identify outliers based on robust (i.e., MAD-based) z scores using the check_outliers()
function of the {performance} package, by specifying method = "zscore_robust"
.^[Note that check_outliers()
only checks numeric variables.] Although @leys2013outliers suggest a default threshold of 2.5 and @leys2019outliers a threshold of 3, {performance} uses by default a less conservative threshold of ~3.29.^[3.29 is an approximation of the two-tailed critical value for p < .001, obtained through qnorm(p = 1 - 0.001 / 2)
. We chose this threshold for consistency with the thresholds of all our other methods.] That is, data points will be flagged as outliers if they go beyond +/- ~3.29 MAD. Users can adjust this threshold using the threshold
argument.
Below we provide example code using the mtcars
dataset, which was extracted from the 1974 Motor Trend US magazine. The dataset contains fuel consumption and 10 characteristics of automobile design and performance for 32 different car models (see ?mtcars
for details). We chose this dataset because it is accessible from base R and familiar to many R users. We might want to conduct specific statistical analyses on this data set, say, t tests or structural equation modelling, but first, we want to check for outliers that may influence those test results.
Because the automobile names are stored as column names in mtcars
, we first have to convert them to an ID column to benefit from the check_outliers()
ID argument. Furthermore, we only really need a couple columns for this demonstration, so we choose the first four (mpg
= Miles/(US) gallon; cyl
= Number of cylinders; disp
= Displacement; hp
= Gross horsepower). Finally, because there are no outliers in this dataset, we add two artificial outliers before running our function.
library(performance) # Create some artificial outliers and an ID column data <- rbind(mtcars[1:4], 42, 55) data <- cbind(car = row.names(data), data) outliers <- check_outliers(data, method = "zscore_robust", ID = "car") outliers
What we see is that check_outliers()
with the robust z score method detected two outliers: cases 33 and 34, which were the observations we added ourselves. They were flagged for two variables specifically: mpg
(Miles/(US) gallon) and cyl
(Number of cylinders), and the output provides their exact z score for those variables.
We describe how to deal with those cases in more details later in the paper, but should we want to exclude these detected outliers from the main dataset, we can extract row numbers using which()
on the output object, which can then be used for indexing:
which(outliers) data_clean <- data[-which(outliers), ]
Other univariate methods are available, such as using the interquartile range (IQR), or based on different intervals, such as the Highest Density Interval (HDI) or the Bias Corrected and Accelerated Interval (BCI). These methods are documented and described in the function's help page.
Univariate outliers can be useful when the focus is on a particular variable, for instance the reaction time, as extreme values might be indicative of inattention or non-task-related behavior^[ Note that they might not be the optimal way of treating reaction time outliers [@ratcliff1993methods; @van1995statistical]].
However, in many scenarios, variables of a data set are not independent, and an abnormal observation will impact multiple dimensions. For instance, a participant giving random answers to a questionnaire. In this case, computing the z score for each of the questions might not lead to satisfactory results. Instead, one might want to look at these variables together.
One common approach for this is to compute multivariate distance metrics such as the Mahalanobis distance. Although the Mahalanobis distance is very popular, just like the regular z scores method, it is not robust and is heavily influenced by the outliers themselves. Therefore, for multivariate outliers, it is recommended to use the Minimum Covariance Determinant, a robust version of the Mahalanobis distance [MCD, @leys2018outliers; @leys2019outliers].
In {performance}'s check_outliers()
, one can use this approach with method = "mcd"
.^[Our default threshold for the MCD method is defined by stats::qchisq(p = 1 - 0.001, df = ncol(x))
, which again is an approximation of the critical value for p < .001 consistent with the thresholds of our other methods.]
outliers <- check_outliers(data, method = "mcd") outliers
Here, we detected 9 multivariate outliers (i.e,. when looking at all variables of our dataset together).
Other multivariate methods are available, such as another type of robust Mahalanobis distance that in this case relies on an orthogonalized Gnanadesikan-Kettenring pairwise estimator [@gnanadesikan1972robust]. These methods are documented and described in the function's help page.
Working with regression models creates the possibility of using model-based SOD methods. These methods rely on the concept of leverage, that is, how much influence a given observation can have on the model estimates. If few observations have a relatively strong leverage/influence on the model, one can suspect that the model's estimates are biased by these observations, in which case flagging them as outliers could prove helpful (see next section, "Handling Outliers").
In {performance}, two such model-based SOD methods are currently available: Cook's distance, for regular regression models, and Pareto, for Bayesian models. As such, check_outliers()
can be applied directly on regression model objects, by simply specifying method = "cook"
(or method = "pareto"
for Bayesian models).^[Our default threshold for the Cook method is defined by stats::qf(0.5, ncol(x), nrow(x) - ncol(x))
, which again is an approximation of the critical value for p < .001 consistent with the thresholds of our other methods.]
Currently, most lm models are supported (with the exception of glmmTMB
, lmrob
, and glmrob
models), as long as they are supported by the underlying functions stats::cooks.distance()
(or loo::pareto_k_values()
) and insight::get_data()
(for a full list of the 225 models currently supported by the insight
package, see https://easystats.github.io/insight/#list-of-supported-models-by-class). Also note that although check_outliers()
supports the pipe operators (|>
or %>%
), it does not support tidymodels
at this time. We show a demo below.
model <- lm(disp ~ mpg * disp, data = data) outliers <- check_outliers(model, method = "cook") outliers
Using the model-based outlier detection method, we identified a single outlier.
Table 1 below summarizes which methods to use in which cases, and with what threshold. The recommended thresholds are the default thresholds.
df <- data.frame( `Statistical Test` = c( "Supported regression model", "Structural Equation Modeling (or other unsupported model)", "Simple test with few variables (*t* test, correlation, etc.)"), `Diagnosis Method` = c( "**Model-based**: Cook (or Pareto for Bayesian models)", "**Multivariate**: Minimum Covariance Determinant (MCD)", "**Univariate**: robust *z* scores (MAD)"), `Recommended Threshold` = c( "_qf(0.5, ncol(x), nrow(x) - ncol(x))_ (or 0.7 for Pareto)", "_qchisq(p = 1 - 0.001, df = ncol(x))_", "_qnorm(p = 1 - 0.001 / 2)_, ~ 3.29"), `Function Usage` = c( '_check_outliers(model, method = "cook")_', '_check_outliers(data, method = "mcd")_', '_check_outliers(data, method = "zscore_robust")_'), check.names = FALSE )
Summary of Statistical Outlier Detection Methods Recommendations
x <- flextable::flextable(df, cwidth = 1.25) x <- flextable::theme_apa(x) x <- flextable::font(x, fontname = "Latin Modern Roman", part = "all") x <- flextable::fontsize(x, size = 10, part = "all") ftExtra::colformat_md(x)
All check_outliers()
output objects possess a plot()
method, meaning it is also possible to visualize the outliers using the generic plot()
function on the resulting outlier object after loading the {see} package (Figure 1).
plot(outliers)
@leys2018outliers report a preference for the MCD method over Cook's distance. This is because Cook's distance removes one observation at a time and checks its corresponding influence on the model each time [@cook1977detection], and flags any observation that has a large influence. In the view of these authors, when there are several outliers, the process of removing a single outlier at a time is problematic as the model remains "contaminated" or influenced by other possible outliers in the model, rendering this method suboptimal in the presence of multiple outliers.
However, distribution-based approaches are not a silver bullet either, and there are cases where the usage of methods agnostic to theoretical and statistical models of interest might be problematic. For example, a very tall person would be expected to also be much heavier than average, but that would still fit with the expected association between height and weight (i.e., it would be in line with a model such as weight ~ height
). In contrast, using multivariate outlier detection methods there may flag this person as being an outlier---being unusual on two variables, height and weight---even though the pattern fits perfectly with our predictions.
Finally, unusual observations happen naturally: extreme observations are expected even when taken from a normal distribution. While statistical models can integrate this "expectation", multivariate outlier methods might be too conservative, flagging too many observations despite belonging to the right generative process. For these reasons, we believe that model-based methods are still preferable to the MCD when using supported regression models. Additionally, if the presence of multiple outliers is a significant concern, regression methods that are more robust to outliers should be considered---like t regression or quantile regression---as they render their precise identification less critical [@mcelreath2020statistical].
The {performance} package also offers an alternative, consensus-based approach that combines several methods, based on the assumption that different methods provide different angles of looking at a given problem. By applying a variety of methods, one can hope to "triangulate" the true outliers (those consistently flagged by multiple methods) and thus attempt to minimize false positives.
In practice, this approach computes a composite outlier score, formed of the average of the binary (0 or 1) classification results of each method. It represents the probability that each observation is classified as an outlier by at least one method. The default decision rule classifies rows with composite outlier scores superior or equal to 0.5 as outlier observations (i.e., that were classified as outliers by at least half of the methods). In {performance}'s check_outliers()
, one can use this approach by including all desired methods in the corresponding argument.
outliers <- check_outliers(model, method = c("zscore_robust", "mcd", "cook")) which(outliers)
Outliers (counts or per variables) for individual methods can then be obtained through attributes. For example:
attributes(outliers)$outlier_var$zscore_robust
An example sentence for reporting the usage of the composite method could be:
Based on a composite outlier score [see the 'check_outliers()' function in the 'performance' R package, @ludecke2021performance] obtained via the joint application of multiple outliers detection algorithms [(a) median absolute deviation (MAD)-based robust z scores, @leys2013outliers; (b) Mahalanobis minimum covariance determinant (MCD), @leys2019outliers; and (c) Cook's distance, @cook1977detection], we excluded two participants that were classified as outliers by at least half of the methods used.
The above section demonstrated how to identify outliers using the check_outliers()
function in the {performance} package. But what should we do with these outliers once identified? Although it is common to automatically discard any observation that has been marked as "an outlier" as if it might infect the rest of the data with its statistical ailment, we believe that the use of SOD methods is but one step in the get-to-know-your-data pipeline; a researcher or analyst's domain knowledge must be involved in the decision of how to deal with observations marked as outliers by means of SOD. Indeed, automatic tools can help detect outliers, but they are nowhere near perfect. Although they can be useful to flag suspect data, they can have misses and false alarms, and they cannot replace human eyes and proper vigilance from the researcher. If you do end up manually inspecting your data for outliers, it can be helpful to think of outliers as belonging to different types of outliers, or categories, which can help decide what to do with a given outlier.
@leys2019outliers distinguish between error outliers, interesting outliers, and random outliers. Error outliers are likely due to human error and should be corrected before data analysis or outright removed since they are invalid observations. Interesting outliers are not due to technical error and may be of theoretical interest; it might thus be relevant to investigate them further even though they should be removed from the current analysis of interest. Random outliers are assumed to be due to chance alone and to belong to the correct distribution and, therefore, should be retained.
It is recommended to keep observations which are expected to be part of the distribution of interest, even if they are outliers [@leys2019outliers]. However, if it is suspected that the outliers belong to an alternative distribution, then those observations could have a large impact on the results and call into question their robustness, especially if significance is conditional on their inclusion, so should be removed.
We should also keep in mind that there might be error outliers that are not detected by statistical tools, but should nonetheless be found and removed. For example, if we are studying the effects of X on Y among teenagers and we have one observation from a 20-year-old, this observation might not be a statistical outlier, but it is an outlier in the context of our research, and should be discarded. We could call these observations undetected error outliers, in the sense that although they do not statistically stand out, they do not belong to the theoretical or empirical distribution of interest (e.g., teenagers). In this way, we should not blindly rely on statistical outlier detection methods; doing our due diligence to investigate undetected error outliers relative to our specific research question is also essential for valid inferences.
Removing outliers can in this case be a valid strategy, and ideally one would report results with and without outliers to see the extent of their impact on results. This approach however can reduce statistical power. Therefore, some propose a recoding approach, namely, winsorization: bringing outliers back within acceptable limits [e.g., 3 MADs, @tukey1963less]. However, if possible, it is recommended to collect enough data so that even after removing outliers, there is still sufficient statistical power without having to resort to winsorization [@leys2019outliers].
The easystats ecosystem makes it easy to incorporate this step into your workflow through the winsorize()
function of {datawizard}, a lightweight R package to facilitate data wrangling and statistical transformations [@patil2022datawizard]. This procedure will bring back univariate outliers within the limits of 'acceptable' values, based either on the percentile, the z score, or its robust alternative based on the MAD.
Finally, it is a critical part of a sound outlier treatment that regardless of which SOD method used, it should be reported in a reproducible manner. Ideally, the handling of outliers should be specified a priori with as much detail as possible, and preregistered, to limit researchers' degrees of freedom and therefore risks of false positives [@leys2019outliers]. This is especially true given that interesting outliers and random outliers are often times hard to distinguish in practice. Thus, researchers should always prioritize transparency and report all of the following information: (a) how many outliers were identified (including percentage); (b) according to which method and criteria, (c) using which function of which R package (if applicable), and (d) how they were handled (excluded or winsorized, if the latter, using what threshold). If at all possible, (e) the corresponding code script along with the data should be shared on a public repository like the Open Science Framework (OSF), so that the exclusion criteria can be reproduced precisely.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.