R語言：異常數據處理
前言
在數據處理中，尤其在作函數擬合時，異常點的出現不僅會很大程度的改變函數擬合的效果，而且有時還會使得函數的梯度出現奇異梯度，這就導致算法的終止，從而影響研究變量之間的函數關系。為了有效的避免這些異常點造成的損失，我們需要采取一定的方法對其進行處理，而處理的第一步便是找到異常點在數據中的位置。
什么是異常值？如何檢測異常值？
目錄
1. 單變量異常值檢測
2. 使用LOF（local outlier factor，局部異常因子）進行異常檢測
3. 通過聚類的方法檢驗異常值
4. 檢驗時間序列數據里面的異常值
5. 討論
主要程序包
install.packages(c("DMwR","dprep"))
library(DMwR)
library(dprep）
1. 單變量異常值檢測
這節主要講單變量異常值檢測，並演示如何將它應用到多元（多個自變量）數據中。使用函數boxplot.stats()實現單變量檢測，該函數根據返回的統計數據生成箱線圖。在上述函數的返回結果中，有一個參數out，它是由異常值組成的列表。更明確的說就是里面列出了箱線圖中箱須線外面的數據點。其中參數coef可以控制箱須線從箱線盒上延伸出來的長度，關於該函數的更多細節可以通過輸入‘？boxplot.ststs’查看。
畫箱線圖：
set.seed（3147）
#產生100個服從正態分布的數據
x <- rnorm(100)
summary(x)
#輸出異常值
boxplot.stats(x)$out
#繪制箱線圖
boxplot(x)
如上的單變量異常檢測可以用來發現多元數據中的異常值，通過簡單搭配的方式。在下例中，我們首先產生一個數據框df，它有兩列x和y。之后，異常值分別從x和y檢測出來。然后，我們獲取兩列都是異常值的數據作為異常數據。
x <- rnorm(100)
y <- rnorm(100)
# 生成一個包含列名分別為x與y的數據框df
df <- data.frame(x, y)
rm(x,y)
head(df)
# 連接數據框df
attach(df)
# 輸出x中的異常值
(a <- which(x %in% boxplot.stats(x)$out))
# 輸出y中的異常值
(b <- which(y %in% boxplot.stats(y)$out))
# 斷開與數據框的連接
detach(df)
# 輸出x,y相同的異常值
(outlier.list1 <- intersect(a,b))
plot(df)
# 標注異常點
points(df[outlier.list1,], col="red", pch="+", cex=2.5)
# x或y中的異常值
(outlier.list2 <- union(a, b))
plot(df)
points(df[outlier.list2,], col="blue", pch="x", cex=2)
當有三個以上的變量時，最終的異常值需要考慮單變量異常檢測結果的多數表決。當選擇最佳方式在真實應用中進行搭配時，需要涉及領域知識。
2. 使用LOF（local outlier factor，局部異常因子）進行異常檢測
LOF（局部異常因子）是一種基於密度識別異常值的算法。算法實現是：將一個點的局部密度與分布在它周圍的點的密度相比較，如果前者明顯的比后者小，那么這個點相對於周圍的點來說就處於一個相對比較稀疏的區域，這就表明該點事一個異常值。（使用LOF，一個點的局部密度會與它的鄰居進行比較。如果前者明顯低於后者（有一個大於1 的LOF值），該點位於一個稀疏區域，對於它的鄰居而言，這就表明，該點是一個異常值。）LOF算法的缺點是它只對數值型數據有效。
lofactor()函數使用LOF算法計算局部異常因子，並且它在DMwR和dprep包中是可用的。下面將介紹一個使用LOF進行異常檢測的例子，k是用於計算局部異常因子的鄰居數量。下圖呈現了一個異常值得分的密度圖。
> library(DMwR)
> # 移除“Species”這個鳶尾花類別列數據
> iris2 <- iris[,1:4]
> # k是計算局部異常因子所需要判斷異常點周圍的點的個數
> outlier.scores <- lofactor(iris2, k=5)
> # 繪制異常值得分的密度分布圖
> plot(density(outlier.scores))
> # 挑出得分排前五的數據作為異常值
> outliers <- order(outlier.scores, decreasing = T)[1:5]
> # 輸出異常值
> print(outliers)
[1] 42 107 23 110 63
接下來對鳶尾花數據進行主成分分析，並利用產生的前兩個主成分繪制成雙標圖來顯示異常值。
> n <- nrow(iris2)
> labels <- 1:n
> # 除了異常值以外所有的數據用"."標注
> labels[-outliers] <- "."
> biplot(prcomp(iris2), cex=.8, xlabs=labels)
上面的代碼中，prcomp()實現對數據集iris2的主成分分析，biplot()取主成分分析結果的前兩列數據也就是前兩個主成分繪制雙標圖。上圖中，x軸和y軸分別代表第一、二主成分，箭頭指向了原始變量名，其中5個異常值分別用對應的行號標注。
我們也可以通過pairs()函數繪制散點圖矩陣來顯示異常值，其中異常值用紅色的'+'標注：
> # 使用rep()生成n個"."
> pch <- rep(".", n)
> pch[outliers] <- "+"
> col <- rep("black", n)
> col[outliers] <- "red"
> pairs(iris2, pch=pch, col=col)
Rlof包，對LOF算法的並行實現。它的用法與lofactor()相似，但是lof()有兩個附加的特性，即支持k的多元值和距離度量的幾種選擇。如下是lof()的一個例子。在計算異常值得分后，異常值可以通過選擇前幾個檢測出來。注意，目前包Rlof的版本在MacOS X和Linux環境下工作，但並不在windows環境下工作，因為它要依賴multicore包用於並行計算。
library(Rlof)
outlier.scores <- lof(iris2, k=5)
# 嘗試使用不同的k值
# try with different number of neighbors (k=5,6,7,8,9 and 10)
outlier.scores <- lof(iris2, k=c(5:10))
3. 通過聚類的方法檢驗異常值
另外一種異常檢測的方法是聚類。通過把數據聚成類，將那些不屬於任務一類的數據作為異常值。比如，使用基於密度的聚類DBSCAN，如果對象在稠密區域緊密相連，它們將被分組到一類。因此，那些不會被分到任何一類的對象就是異常值。
我們也可以使用k-means算法來檢測異常。使用k-means算法，數據被分成k組，通過把它們分配到最近的聚類中心。然后，我們能夠計算每個對象到聚類中心的距離（或相似性），並且選擇最大的距離作為異常值。
#####################################
iris2 <- iris[,1:4]
kmeans.result <- kmeans(iris2, centers=3)
#class(kmeans.result)
# 輸出簇中心
kmeans.result$centers
#length(kmeans.result$centers)
# 分類結果
kmeans.result$cluster
#mode(kmeans.result$cluster)
# 計算數據對象與簇中心的距離
centers <- kmeans.result$centers[kmeans.result$cluster, ]
#class(centers)
distances <- sqrt(rowSums((iris2 - centers)^2))
# 挑選出前5個最大距離
outliers <- order(distances, decreasing = T)[1:5]
# 輸出異常值
print(outliers)
print(iris2[outliers,])
# 畫出聚類結果
plot(iris2[,c("Sepal.Length", "Sepal.Width")], pch="o", col=kmeans.result$cluster, cex=0.3)
# 繪制類（簇）中心，用"*"標記
points(kmeans.result$centers[,c("Sepal.Length", "Sepal.Width")], col=1:3, pch=8, cex=1.5)
# 畫出異常值，用"+"標記
points(iris2[outliers,c("Sepal.Length", "Sepal.Width")], pch="+", col=4, cex=1.5)
#####################################
#####################################
> iris2 <- iris[,1:4]
> kmeans.result <- kmeans(iris2, centers=3)
> #class(kmeans.result)
> # 輸出簇中心
> kmeans.result$centers
Sepal.Length Sepal.Width Petal.Length Petal.Width
1 6.314583 2.895833 4.973958 1.7031250
2 5.175758 3.624242 1.472727 0.2727273
3 4.738095 2.904762 1.790476 0.3523810
> #length(kmeans.result$centers)
> # 分類結果
> kmeans.result$cluster
[1] 2 3 3 3 2 2 2 2 3 3 2 2 3 3 2 2 2 2 2 2 2 2 2 2 3 3 2 2 2 3 3 2
[33] 2 2 3 2 2 2 3 2 2 3 3 2 2 3 2 3 2 2 1 1 1 1 1 1 1 3 1 1 3 1 1 1
[65] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1
[97] 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[129] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
> #mode(kmeans.result$cluster)
> # 計算數據對象與簇中心的距離
> centers <- kmeans.result$centers[kmeans.result$cluster, ]
> #class(centers)
> distances <- sqrt(rowSums((iris2 - centers)^2))
> # 挑選出前5個最大距離
> outliers <- order(distances, decreasing = T)[1:5]
> # 輸出異常值
> print(outliers)
[1] 119 118 132 123 106
> print(iris2[outliers,])
Sepal.Length Sepal.Width Petal.Length Petal.Width
119 7.7 2.6 6.9 2.3
118 7.7 3.8 6.7 2.2
132 7.9 3.8 6.4 2.0
123 7.7 2.8 6.7 2.0
106 7.6 3.0 6.6 2.1
> # 畫出聚類結果
> plot(iris2[,c("Sepal.Length", "Sepal.Width")], pch="o", col=kmeans.result$cluster, cex=0.3)
> # 繪制類（簇）中心，用"*"標記
> points(kmeans.result$centers[,c("Sepal.Length", "Sepal.Width")], col=1:3, pch=8, cex=1.5)
> # 畫出異常值，用"+"標記
> points(iris2[outliers,c("Sepal.Length", "Sepal.Width")], pch="+", col=4, cex=1.5)
#####################################
在上圖中，聚類中心被標記為星號，異常值標記為'+'
4. 檢驗時間序列數據里面的異常值
本部分講述一個對時間序列數據進行異常檢測的例子。在本例中，時間序列數據首次使用stl()進行穩健回歸分解，然后識別異常值。
########################################
# 使用穩健回歸擬合
f <- stl(AirPassengers, "periodic", robust=TRUE)
(outliers <- which(f$weights < 1e-8))
# 繪圖布局
op <- par(mar=c(0, 4, 0, 3), oma=c(5, 0, 4, 0), mfcol=c(4, 1))
plot(f, set.pars=NULL)
sts <- f$time.series
# 畫出異常值，用紅色"x"標記
points(time(sts)[outliers], 0.8*sts[,"remainder"][outliers], pch="x", col="red")
par(op)
########################################
> # 使用穩健回歸擬合
> f <- stl(AirPassengers, "periodic", robust=TRUE)
> (outliers <- which(f$weights < 1e-8))
[1] 79 91 92 102 103 104 114 115 116 126 127 128 138 139 140
> # 繪圖布局
> op <- par(mar=c(0, 4, 0, 3), oma=c(5, 0, 4, 0), mfcol=c(4, 1))
> plot(f, set.pars=NULL)
> sts <- f$time.series
> # 畫出異常值，用紅色"x"標記
> points(time(sts)[outliers], 0.8*sts[,"remainder"][outliers], pch="x", col="red")
> par(op)
########################################
在上圖中，異常值用紅色標記為'x'
5. 討論
LOF算法擅長檢測局部異常值，但是它只對數值數據有效。Rlof包依賴multicore包，在Windows環境下失效。對於分類數據的一個快速穩定的異常檢測的策略是AVF(Attribute Value Frequency)算法。
一些用於異常檢測的R包包括：
extremevalues包：單變量異常檢測
mvoutlier包：基於穩定方法的多元變量異常檢測
outliers包：對異常值進行測驗