今天的文章介绍一篇论文[1]，涵盖了几种常见的识别异常值的方法，包括单变量方法，多变量方法，和基于模型的方法，以及在R中的实现。相信通过掌握上述的几种方法，可以助各位更加从容淡定的解决异常值的相关问题。
install.packages("performance")
install.packages("see")
install.packages("rempsyc")
install.packages("bigutilsr")
library(performance)
library(see)
library(rempsyc)
出于简化结果的目的，将提取前四个变量，再通过手动增加两行含异常值的样本，以及添加一个新的变量作为识别汽车的名字，代码如下： mydata <- rbind(mtcars[1:4], 42, 55)
mydata <- cbind(car_names = row.names(mydata), mydata)
mydata
mtcars本身一共包含32行，上述代码增加了两行，即33和34行。
另外，给4个汽车的变量添加了数值（最后两行）。最后两行是手动添加的异常值，目的在于使用后续的方法将它们检测出来。关于单变量方法，许多研究者会考虑使用z score来检测异常值，它代表样本偏离均数的标准差个数。如果偏离的太多，则将它们作为异常值处理。
虽然这种做法比较普遍，但是，这种方法存在缺陷。因为均数和标准差本身就容易受到异常值的影响，并且这个方法需假定正态分布。因此，目前处理异常值的相关指南往往推荐鲁棒方法。第一例子将使用z score的鲁棒版本，即robust z score （基于中位数绝对偏差），代码如下： outliers <- check_outliers(mydata, method = "zscore_robust", ID = "car_names")
outliers
从输出结果的第一行所知，这个方法检测到异常值所在行数，即33和34 （由我们自己创建）。第二行注明了使用的方法为zscore_robust以及阈值为3.291（超过这个数值则被判定为异常值）。
关于单变量方法，可以查看这个函数的帮助文档，如下图所示。关于为什么阈值设定为3.291，因为它是 p < .001 的双尾临界值的近似值。下方的代码可以帮助快速得到存在异常值的行数，有时也会比较实用，代码如下：
这对于快速剔除这个两行数据会非常帮助。假设现在希望剔除这两行，并且建立一个新的数据，取名叫mydata_clean，代码如下 : mydata_clean <- mydata[-which(outliers), ]
mydata_clean
与mydata相比，mydata_clean将33和34行剔除出去，可用于后续的分析。还可以将结果进行作图，可以高亮出异常值，代码如下： 关于多变量方法，可以使用函数的帮助文档进行查看，如下图：
outliers1 <- check_outliers(mydata, method = "mahalanobis_robust")
outliers1
从结果可知，这种方法识别出了3行包含异常值的样本。值得注意的是，不同方法得出的结果可能是不同的，这在统计里面是常态。因为与单变量方法不同的是，多变量方法考虑了数据的多变量距离。
关于基于模型的方法，也可以查看函数的帮助文档，如下： 这里以Cook's distance为例，先建立一个模型，再检测异常值，代码如下： mymodel <- lm(mpg ~ hp * cyl, data = mydata)
outliers2 <- check_outliers(mymodel, method = "cook")
outliers2
在这个模型中，此方法识别出了一个异常值，为34行样本。
关于如何选择不同的检测方法，原文作者制作了一个推荐的表格，如下：
![]()
参考文献的表1再举一个例子用于说明不同方法会产生不同的结果，以及相应的解决办法。创建一个数据集，最后两行为手动添加的异常值，代码如下：
mydata1 <- rbind(women[rep(seq_len(nrow(women)), each = 100), ],
c(100, 258),
c(100, 200))
summary(mydata1)
建立一个回归模型，再将结果作图（且包含回归直线），代码如下：
mymodel1 <- lm(weight ~ height, mydata1)
nice_scatter(mydata1, "height", "weight")
图片右侧两个点为手动添加的“异常值”，一个分布在回归直线上，而另外一个在回归直线下方。首先使用robust z socre法，代码如下：
outliers3 <- check_outliers(mymodel1, method = "zscore_robust")
which(outliers3)
与预期一致，因为这个方法为单变量方法，它只考虑数值在某一个样本的位置，它的算法中并没有考虑回归直线这个问题。
下一步，再使用Cook's distance法，代码如下：
outlier4 <- check_outliers(mymodel1, method = "cook")
which(outlier4)
此方法识别出了上述散点图中右侧下方的高杠杆点：1502。
那有没有方法可以考虑多种方法，给出一个“中庸”的结果呢？
答案是肯定的，这个函数可以给出复合异常值评分，为一种基于共识的方法。此方法比较民主，假设ABC三个人投票，原则是少数服从多数。AB两个都说让C来请客吃饭，C就必须要请客吃饭![]()
！这里以三种方法为例，即三个人投票 （cook, zscore_robust, iqr），R代码如下： outliers5 <- check_outliers(mymodel1,
method = c("cook", "zscore_robust", "iqr"))
which(outliers5)
从结果可知，应该至少有两种方法是支持两个样本存在异常值的。
check_outliers(mymodel1, method = c("cook"))
check_outliers(mymodel1, method = c("zscore_robust"))
check_outliers(mymodel1, method = c("iqr"))
最后，原文作者认为，处理异常值的时候最重要的是需要做到透明。最好是将所有关于异常值的处理的信息都进行汇报，比如一共检测到多少异常值；使用什么R包，什么函数，什么方法；它们是如何处理的，是直接剔除还是保留等等。
好啦，今天的内容就到这里。如果有帮助，记得分享给需要的人![]()
！1. Check your outliers! An introduction to identifying statistical outliers in R with easystats
【通过公众号菜单栏--线上课程--新手课程/回归课程】公众号核心成员担任SCI杂志Associate Editor!▌课程相关咨询可添加R师妹微信: kefu_rstats
