KeepNotes blog

Stay hungry, Stay Foolish.

0%

识别离群点-ESD(Generalized ESD)

Posted on Edited on In ,
Symbols count in article: 2.8k Reading time ≈ 3 mins.

Generalized ESD Test (ESD)是Rosner教授基于Grubb's Test(或extreme studentized deviate (ESD) test)改进的识别离散值的方法

因为ESD的备择假设是数据集中有一个异常值，而现实情况下数据集中异常值不止一个；因此Rosner提出了GESD（泛化版ESD）

其适用于离散数未知的情况下，假设数据没有任何离散值，检验服从正态分布的单变量数据集，对数据集中的k个潜在离散值执行generalized ESD Test

在CLSI的EP09-A3 (Measurement Procedure Comparison and Bias Estimation Using Patient Samples)中对于离群点识别采用了generalized ESD方法，并对其计算步骤做了详细的说明

CLSI是美国【临床实验室标准化协会】的英文缩写，英文名为Clinical and Laboratory Standards Institute。目前在体外诊断试剂领域得到FDA认可的共识标准均属于CLSI制定的标准

Generalized ESD Test方法在R中已有包可调用，具体内部的源码还未查看，未能确定是否跟CLSI中算法一致，但从测试结果上来看是一致的，可用以下R包中的函数：

  • EnvStats包的rosnerTest函数
  • PMCMRplus包的gesdTest函数

Generalized ESD在EP中的计算过程如下：

  • 设置显著水平alpha，一般为0.05或者0.01
  • 指定离群点比例h，一般假定h=5%，如果有50个samples，那么潜在的离群点数目则为2，向下取整
  • 计算数据集均值和标准差，找出每个值与均值的绝对差值的最大值除以标准差，记作ESD_i，被认为是一个潜在的离群点；然后依次从数据集中剔除该点重复上述操作（计算ESD_i），这样得到h个ESD_i值（对应第i次计算的ESD值）
  • 计算critical vlaue: lambda_i
  • 最后统计每个i对应的EDS_i值是否大于lambda_i，取最大的i值作为离群点的数目；即上述1到h循环中，ESD_i大于lambda_i的条件下，最大的i为5，则说明数据集中有5个离群点

上述文字表达并不清楚，但考虑到CLSI是收费文档，在此就不补充原文了，具体算法简述可参考文献献：Comparison of methods for detecting outliers

根据上述，代码如下：

# Compute the critical value for ESD Test
esd.critical <- function(alpha, N, i) {
  v <- N - i - 1
  p <- 1 - alpha / (2 * (N - i + 1))
  t <- qt(p, v)
  lambda <- t * (N - i) / sqrt((N - i + 1) * (v + t^2))
  return(lambda)
}

# Calculate ESD(generalized ESD)
ESD_detect <- function(x, level = 0.05) {
  # Define values and vectors.
  x2 <- x
  n <- length(x)
  alpha <- 0.05
  h <- round(level * n) + 1
  
  # Removed from the sample base on ESDi
  res <- data.frame(stringsAsFactors = F)
  for (i in 1:h) {
    if (sd(x2) == 0){
      break
    } 
    temp <- abs(x2 - mean(x2)) / sd(x2)
    esdi <- max(temp)
    index <- which(temp == esdi)[1]
    
    lambda <- esd.critical(alpha, n, i)
    res <- rbind(res, data.frame(ID = x2[index], i = i, ESDi = esdi, lambda = lambda))
    
    ## iterated remove esdi
    x2 <- x2[temp != esdi]
  }
  
  # Keep the non-outliers
  # If ESDh and ESDh+1 are equal(a tie) then neither one should be seen as an outlier
  # The number of outliers is determined by finding the largest i such that ESDi > λi
  if (res$ESDi[h] == res$ESDi[h-1]) {
    res <- res[-((h-1):h),]
    if (length(which(res$ESDi > res$lambda)) > 0) {
      outlier <- res[1:which(res$ESDi > res$lambda),]
    }else{
      outlier <- NULL
    }
    
  }else{
    if (length(which(res$ESDi > res$lambda)) > 0) {
      outlier <- res[1:which(res$ESDi > res$lambda),]
    }else{
      outlier <- NULL
    }
  }
  
  out <- list(results = res, outliers = outlier)
  
  return(out)
}

若有理解错误的地方，请告知下；https://github.com/kaigu1990/CLSI-Statistic/blob/master/ESD.R

使用示例：

y = c(-0.25, 0.68, 0.94, 1.15, 1.20, 1.26, 1.26,
      1.34, 1.38, 1.43, 1.49, 1.49, 1.55, 1.56,
      1.58, 1.65, 1.69, 1.70, 1.76, 1.77, 1.81,
      1.91, 1.94, 1.96, 1.99, 2.06, 2.09, 2.10,
      2.14, 2.15, 2.23, 2.24, 2.26, 2.35, 2.37,
      2.40, 2.47, 2.54, 2.62, 2.64, 2.90, 2.92,
      2.92, 2.93, 3.21, 3.26, 3.30, 3.59, 3.68,
      4.30, 4.64, 5.34, 5.42, 6.01)

ESD_detect(y, level = 0.2)

参考资料

https://www.rdocumentation.org/packages/PMCMRplus/versions/1.4.4/topics/gesdTest
http://finzi.psych.upenn.edu/R/library/EnvStats/html/rosnerTest.html
https://zhuanlan.zhihu.com/p/34342025
https://www.statisticshowto.com/generalized-esd-test/

本文出自于http://www.bioinfo-scrounger.com转载请注明出处

Related Posts