如何在R代碼中使用異常值測試

Question

作爲我的數據分析工作流程的一部分，我想測試異常值，然後在有和沒有異常值的情況下進行進一步計算。

我發現了異常值包，它有各種測試，但我不確定如何最好地將它們用於我的工作流。

Answer 1

我同意德克，這很難。我會建議先看看爲什麼你可能有異常值。異常值只是某人認爲可疑的一個數字，並不是一個具體的「壞」值，除非你能找到一個異常值的理由，否則你可能不得不忍受這種不確定性。

你沒有提到的一件事是你在看什麼樣的離羣值。你的數據是否集中在一個平均值上，它們是否有特定的分佈或者你的數據之間是否存在某種關係。

下面是一些例子

首先，我們將創建一些數據，然後用一個局外人玷污它;

> 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 

或者更復雜的例子，你可以用統計數據來計算臨界截止值，在這裏使用Lund測試（參見Lund，RE 1975，「用於線性模型中的離羣值的近似測試的表格」，Technometrics，第17卷，第4期，第473-476頁和Prescott， 。。。檢驗異常值的線性模型」，Technometrics，第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測試產生的臨界值應該與臨牀殘留的絕對值（即無符號）進行比較。

Answer 2

「很難」。這在很大程度上是依賴於上下文，你可能必須嵌入到你的應用程序是：

• 數據是否漂移，趨勢或循環？
• 變異性是固定的還是變量本身？
• 是否有其他系列可用於「基準測試」？

除了異常值包以外，還有qcc包，因爲質量控制文獻涵蓋此區域。

還有很多其他領域你可以看看，例如， robust statistics Task View。

Answer 3

如果您擔心異常值，而不是拋出異常值，請使用穩健的方法。例如，而不是lm，請使用rlm。

Answer 4

試試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