在R中进行基于稳健马氏距离的异常检验

转自http://www.cnblogs.com/cloudtj/articles/5520230.html

对应书本《医学数据挖掘案例与实践》第一章1.1小节

前言

  我们研究的数据中经常包含着一些不同寻常的样本，这称之为异常值(Outlier)。这些异常值会极大的影响回归或分类的效果。异常值产生的原因有很多，其中可能是人为错误、数据测量误差，或者是实际确实存在这样的异常。为了使模型能够反映大部分数据的规律，所以在数据预处理阶段要进行异常值检测，为下一步分析奠定基础。还有一类情况是，当研究人员希望发现不平凡的事物时，异常值检测本身就是分析的首要目的。例如在信用卡欺诈、计算机入侵检测等问题中。此时由于样本的不平衡性，导致一般的分类方法无法使用，必须转而考虑异常检测方法。

  一种常用的异常检验思路是观察各样本点到样本中心的距离。如果某些样本点的距离太大，就可以判断是异常值。这里距离的度量一般使用马氏距离(Mahalanobis Distance)。因为马氏距离不受量纲的影响，而且在多元条件下，马氏距离还考虑了变量之间的相关性，这使得它优于欧氏距离。

  但是传统的马氏距离检测方法是不稳定的，因为个别异常值会把均值向量和协方差矩阵向自己方向吸引，这样算出来的样本马氏距离起不了检测异常值的所用。所以首先要利用迭代的思想构造一个稳健的均值和协方差矩阵估计量，然后计算稳健马氏距离(Robust Mahalanobis Distance)。这样使得异常值能够正确地被识别出来。

  在mvoutlier包中提供了基于稳健马氏距离的异常值检验方法。我们首先构造一个二维变量的人工数据，其中80个样本是标准正态分布，另一小撮别有用心的样本是均值为5，标准差为1的观测值。我们首先使用uni.plot函数在一维空间中观察这个数据。

library(mvoutlier)
    set.seed(1234)
    x <- cbind(rnorm(80), rnorm(80))
    y <- cbind(rnorm(10, 5, 1), rnorm(10, 5, 1))
    z <- rbind(x,y)
    # 一维数据的异常检验
    res1 <- uni.plot(z)
    # 返回异常值的编号
    which(res1$outliers==T)
    ################################
    > library(mvoutlier)
    > set.seed(1234)
    > x <- cbind(rnorm(80), rnorm(80))
    > y <- cbind(rnorm(10, 5, 1), rnorm(10, 5, 1))
    > z <- rbind(x,y)
    > # 一维数据的异常检验
    > res1 <- uni.plot(z)
    > # 返回异常值的编号
    > which(res1$outliers==T)
    [1] 20 62 81 82 83 84 85 86 87 88 89 90

上图中红色点表示疑似异常值，因为它偏离均值太远。更多时候我们会处理多元异常检测问题，此时用aq.plot函数来实行基于稳健马氏距离的异常值检验方法。下图中左上角图形为原始数据，右上角图形的X轴为各样本的稳健马氏距离排序，Y轴为距离的经验分布，红色曲线为卡方分布，蓝色垂线表示阀值，在阀值右侧的样本判断为异常值。左下和右下两张图均是用不同颜色来表示异常值，只是阀值略有不同。可以观察到那一小撮异常值被正确的判断出来，但也有两个正常值被误判为异常值，此时需要调整参数。

 # 基于稳健马氏距离的多元异常值检验
    res2 <-aq.plot(z)
    # 返回异常值的编号
    which(res2$outliers==T)
    ################################
    > res2 <-aq.plot(z)
    > which(res2$outliers==T)
     [1] 20 62 81 82 83 84 85 86 87 88 89 90

如果数据的维数过高，例如基因数据那样几千个变量，数据之间变得稀疏，从而使得距离不再有很大意义。此时可以融合主成分降维的思路来进行异常值检验。mvoutlier包中提供了pcout函数来进行高维空间异常检验。下面是以swiss数据集为例来判断异常值。

# 在高维空间中的异常值检验
    data(swiss)
    res3 <- pcout(swiss)
    # 返回异常值的编号
    which(res3$wfinal01==0)
    ################################
    > # 在高维空间中的异常值检验
    > data(swiss)
    > res3 <- pcout(swiss)
    > # 返回异常值的编号
    > which(res3$wfinal01==0)
        Delemont Franches-Mnt   Porrentruy        Broye        Glane 
               2            3            6            7            8 
         Gruyere       Sarine      Veveyse    La Vallee      Conthey 
               9           10           11           19           31 
       Entremont       Herens     Martigwy      Monthey   St Maurice 
              32           33           34           35           36 
          Sierre         Sion V. De Geneve 
              37           38           45