如何在R代码中使用异常值测试

如果你担心异常值,而不是把它们扔掉,那么使用一个健壮的方法。例如,使用rlm而不是lm。

我同意德克的观点,这很难。我建议先看看你为什么会有异常值。离群值只是一个有人认为可疑的数字,它不是一个具体的“坏”值,除非你能找到离群值的原因,否则你可能不得不忍受不确定性。

有一件事你没有提到,那就是你在看什么样的异类。您的数据是围绕平均值聚集的,它们是否具有特定的分布,或者您的数据之间是否存在某种关系。

下面是一些例子

首先,我们将创建一些数据,然后用离群值来污染它;

> testout<-data.frame(X1=rnorm(50,mean=50,sd=10),X2=rnorm(50,mean=5,sd=1.5),Y=rnorm(50,mean=200,sd=25))
> #Taint the Data
> testout$X1[10]<-5
> testout$X2[10]<-5
> testout$Y[10]<-530

> testout
         X1         X2        Y
1  44.20043  1.5259458 169.3296
2  40.46721  5.8437076 200.9038
3  48.20571  3.8243373 189.4652
4  60.09808  4.6609190 177.5159
5  50.23627  2.6193455 210.4360
6  43.50972  5.8212863 203.8361
7  44.95626  7.8368405 236.5821
8  66.14391  3.6828843 171.9624
9  45.53040  4.8311616 187.0553
10  5.00000  5.0000000 530.0000
11 64.71719  6.4007245 164.8052
12 54.43665  7.8695891 192.8824
13 45.78278  4.9921489 182.2957
14 49.59998  4.7716099 146.3090
<snip>
48 26.55487  5.8082497 189.7901
49 45.28317  5.0219647 208.1318
50 44.84145  3.6252663 251.5620

以图形方式检查数据通常最有用(你的大脑比数学更擅长发现异常值)。

> #Use Boxplot to Review the Data
> boxplot(testout$X1, ylab="X1")
> boxplot(testout$X2, ylab="X2")
> boxplot(testout$Y, ylab="Y")

然后您可以使用测试。如果测试返回截止值,或者可能是异常值,则可以使用ifelse删除它。

> #Use Outlier test to remove individual values
> testout$newX1<-ifelse(testout$X1==outlier(testout$X1),NA,testout$X1)
> testout
         X1         X2        Y    newX1
1  44.20043  1.5259458 169.3296 44.20043
2  40.46721  5.8437076 200.9038 40.46721
3  48.20571  3.8243373 189.4652 48.20571
4  60.09808  4.6609190 177.5159 60.09808
5  50.23627  2.6193455 210.4360 50.23627
6  43.50972  5.8212863 203.8361 43.50972
7  44.95626  7.8368405 236.5821 44.95626 
8  66.14391  3.6828843 171.9624 66.14391 
9  45.53040  4.8311616 187.0553 45.53040
10  5.00000  5.0000000 530.0000       NA 
11 64.71719  6.4007245 164.8052 64.71719 
12 54.43665  7.8695891 192.8824 54.43665 
13 45.78278  4.9921489 182.2957 45.78278 
14 49.59998  4.7716099 146.3090 49.59998 
15 45.07720  4.2355525 192.9041 45.07720 
16 62.27717  7.1518606 186.6482 62.27717 
17 48.50446  3.0712422 228.3253 48.50446 
18 65.49983  5.4609713 184.8983 65.49983 
19 44.38387  4.9305222 213.9378 44.38387 
20 43.52883  8.3777627 203.5657 43.52883 
<snip>
49 45.28317  5.0219647 208.1318 45.28317 
50 44.84145  3.6252663 251.5620 44.84145

或者对于更复杂的例子,您可以使用stats来计算临界截止值,这里使用lund测试(参见lund,R.E.1975,“线性模型中离群值的近似测试表”,Technometrics,第17卷,第4期,第473-476页。和Prescott,P.1975,“线性模型中异常值的近似测试”,《技术度量》,第17卷,第1期,第129-132页。)

> #Alternative approach using Lund Test
> lundcrit<-function(a, n, q) {
+ # Calculates a Critical value for Outlier Test according to Lund
+ # See Lund, R. E. 1975, "Tables for An Approximate Test for Outliers in Linear Models", Technometrics, vol. 17, no. 4, pp. 473-476.
+ # and Prescott, P. 1975, "An Approximate Test for Outliers in Linear Models", Technometrics, vol. 17, no. 1, pp. 129-132.
+ # a = alpha
+ # n = Number of data elements
+ # q = Number of independent Variables (including intercept)
+ F<-qf(c(1-(a/n)),df1=1,df2=n-q-1,lower.tail=TRUE)
+ crit<-((n-q)*F/(n-q-1+F))^0.5
+ crit
+ }

> testoutlm<-lm(Y~X1+X2,data=testout)

> testout$fitted<-fitted(testoutlm)

> testout$residual<-residuals(testoutlm)

> testout$standardresid<-rstandard(testoutlm)

> n<-nrow(testout)

> q<-length(testoutlm$coefficients)

> crit<-lundcrit(0.1,n,q)

> testout$Ynew<-ifelse(abs(testout$standardresid)>crit,NA,testout$Y)

> testout
         X1         X2        Y    newX1   fitted    residual standardresid
1  44.20043  1.5259458 169.3296 44.20043 209.8467 -40.5171222  -1.009507695
2  40.46721  5.8437076 200.9038 40.46721 231.9221 -31.0183107  -0.747624895
3  48.20571  3.8243373 189.4652 48.20571 203.4786 -14.0134646  -0.335955648
4  60.09808  4.6609190 177.5159 60.09808 169.6108   7.9050960   0.190908291
5  50.23627  2.6193455 210.4360 50.23627 194.3285  16.1075799   0.391537883
6  43.50972  5.8212863 203.8361 43.50972 222.6667 -18.8306252  -0.452070155
7  44.95626  7.8368405 236.5821 44.95626 223.3287  13.2534226   0.326339981
8  66.14391  3.6828843 171.9624 66.14391 148.8870  23.0754677   0.568829360
9  45.53040  4.8311616 187.0553 45.53040 214.0832 -27.0279262  -0.646090667
10  5.00000  5.0000000 530.0000       NA 337.0535 192.9465135   5.714275585
11 64.71719  6.4007245 164.8052 64.71719 159.9911   4.8141018   0.118618011
12 54.43665  7.8695891 192.8824 54.43665 194.7454  -1.8630426  -0.046004311
13 45.78278  4.9921489 182.2957 45.78278 213.7223 -31.4266180  -0.751115595
14 49.59998  4.7716099 146.3090 49.59998 201.6296 -55.3205552  -1.321042392
15 45.07720  4.2355525 192.9041 45.07720 213.9655 -21.0613819  -0.504406009
16 62.27717  7.1518606 186.6482 62.27717 169.2455  17.4027250   0.430262983
17 48.50446  3.0712422 228.3253 48.50446 200.6938  27.6314695   0.667366651
18 65.49983  5.4609713 184.8983 65.49983 155.2768  29.6214506   0.726319931
19 44.38387  4.9305222 213.9378 44.38387 217.7981  -3.8603382  -0.092354925
20 43.52883  8.3777627 203.5657 43.52883 228.9961 -25.4303732  -0.634725264
<snip>
49 45.28317  5.0219647 208.1318 45.28317 215.3075  -7.1756966  -0.171560291
50 44.84145  3.6252663 251.5620 44.84145 213.1535  38.4084869   0.923804784
       Ynew
1  169.3296
2  200.9038
3  189.4652
4  177.5159
5  210.4360
6  203.8361
7  236.5821
8  171.9624
9  187.0553
10       NA
11 164.8052
12 192.8824
13 182.2957
14 146.3090
15 192.9041
16 186.6482
17 228.3253
18 184.8983
19 213.9378
20 203.5657
<snip>
49 208.1318
50 251.5620

编辑:我刚刚注意到我的代码中有一个问题。lund测试产生一个临界值,该临界值应与研究残余物的绝对值(即无符号)进行比较。

“这很难”。其中大部分依赖于上下文,您可能需要将其嵌入到应用程序中:

• 数据是否漂移、趋势或周期?
• 变异性是固定的还是本身是可变的?
• 是否有其他系列可用于“基准测试”?

除了离群值包,还有 qcc 包装,因为质量控制文献涵盖了这一领域。

您可以看到许多其他方面,例如 robust statistics Task View .

试试 outliers::score 功能。我不建议删除所谓的离群值,但知道你的极端观察是好的。

library(outliers)
set.seed(1234)
x = rnorm(10)
[1] -1.2070657  0.2774292  1.0844412 -2.3456977  0.4291247  0.5060559 -0.5747400 -0.5466319
[9] -0.5644520 -0.8900378
outs <- scores(x, type="chisq", prob=0.9)  # beyond 90th %ile based on chi-sq
#> [1] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
x[outs]  # most extreme
#> [1] -2.345698

你会找到更多的帮助 outlier detection here