应用稳健马氏距离评价实验室间比对能力考核
2017-06-05师耀龙吴晓凤柴文轩楚宝临
师耀龙，滕 曼，李 成，吴晓凤，柴文轩，杨 婧，楚宝临,付 强
1.中国环境监测总站，国家环境保护环境监测质量控制重点实验室，北京 1000122.河北农业大学科学技术研究院,河北 保定 071001
应用稳健马氏距离评价实验室间比对能力考核
师耀龙1，滕 曼1，李 成2，吴晓凤1，柴文轩1，杨 婧1，楚宝临1,付 强1
1.中国环境监测总站，国家环境保护环境监测质量控制重点实验室，北京 1000122.河北农业大学科学技术研究院,河北 保定 071001
基于稳健马氏距离的多元统计算法，综合评价了全国环境监测系统96个实验室4种有机氯农药监测能力考核的结果，从中筛选出25个考核结果存在异常的实验室。将多元统计结果与稳健Z比分数结果进行比较，发现基于稳健马氏距离筛选出的异常实验室其稳健Z比分数结果也存在一定的问题。同时简单介绍了稳健马氏距离及其相关的异常值筛选的算法原理，并详细介绍了其基于R语言的实现过程。
多元稳健统计；稳健马氏距离；实验室间比对；能力考核
环境监测实验室能力考核是中国环境监测总站(以下简称总站)组织的利用实验室间比对的方式考核环境监测实验室技术能力和质量管理水平的活动。通过实验室间比对，有助于总站了解领域内实验室相关监测项目的能力状况，及时发现和纠正影响检测水平的影响因素，促进领域内实验室检测能力的提高[1-2]。
伴随着能力考核工作的开展，其数据评价工作中出现了2个主要问题：①一些实验室在周期内参加了多轮能力考核，或在一轮能力考核中进行了多个样品的测试，传统的评价方法主要是对各实验室逐个项目或逐个轮次的评价[1-2]，如何用一项综合指标整体评价实验室能力尚待研究；②由于能力考核结果不服从或近似服从正态分布，需要对传统统计量进行稳健化，如果希望整体评价，则需要对其参与的各个项目构成的多元数据进行稳健统计，而稳健Z比分数和迭代法都是对能力考核的一元数据进行稳健统计的方法[3]，对能力考核的多元数据进行稳健统计的方法尚待研究。
鉴于此，研究详细讨论了稳健马氏距离在有机氯监测能力考核结果评价中的应用，并与稳健Z比分数评价结果进行了比较[1-4]。此外，由于马氏距离与多元数据的稳健统计方法在中国环境监测领域应用较少，研究对其一般原理与基于开源软件(R)的实现方法进行了详细介绍。
1 一般马氏距离
在能力考核多元数据处理中，马氏距离具有：①不受不同项目之间量纲不同的影响，所得距离为标准化距离；②计算过程考虑到了各项目之间相关性的存在，是建立在总体协方差矩阵上的标准化距离。马氏距离已应用于包括环境监测领域在内的各个领域内多元数据异常值筛选工作中[6]。但是，传统的统计方法中，异常值的存在会显著影响中心值和协方差矩阵的估计，使一般马氏距离不能正确反映各个观测的偏离程度。对于这类数据，需要通过稳健统计的方法，构建稳定的均值和协方差矩阵统计量。
2 稳健马氏距离及其相关R程序包介绍
开源R软件中的“robustbase”程序包[7-8]中的covMcd程序是基于Fast-MCD算法和PISON等人在2002年的改进完成的[7, 9-11]，广泛用于估计多元样本中的稳健统计量。其算法原理较为复杂，可简单的概括为从样本中选择h个观测，通过不断的迭代计算h个观测最小的协方差行列式，该协方差矩阵通过加权即可估计出稳健的统计量，详细的计算方法参见文献[7,10-12]。稳健估计结束后，以稳健的中心值和协方差，通过马氏距离计算各观测向量偏离中心值的稳健马氏距离，并可根据DM(x)2符合卡方分布的特点，筛选出数据集中的异常值，这一过程可由R软件中的“mvoutlier”程序包实现[13]。
pn(δ)=sup(G(u)-Gn(u))+u≥δ
同时，计算pcrit(δ,n,p)
计算完成后，比较pn(δ)和pcrit(δ,n,p)，若pn(δ)
“mvoutlier”程序包同时支持多种计算和画图功能，并能从单个变量观测整个数据集，分析哪些变量更容易导致异常值的出现，可用于能力考核中不合格站点的筛选与考核项目难易的判断。
研究将通过稳健马氏距离在有机氯监测能力考核中的应用实例来介绍稳健马氏距离在实验室能力考核中的应用。
3 稳健马氏距离在环境监测实验室能力考核中的应用
3.1 数据来源与R软件计算过程
总站于2013年开展了针对全国各地市环境监测站的有机氯监测能力考核，数据选取了96个实验室测定的α-六氯环己烷(简写为α-666)样品A、α-666样品B、γ-六氯环己烷(简写为γ-666)样品A、γ-666样品B、p, p′-双对氯苯基三氯乙烷(简写为p, p′-DDT)样品A、p, p′-DDT样品B、o, p′-双对氯苯基三氯乙烷(简写为o, p′-DDT)样品A和o, p′-DDT样品B的含量，构成一个96×8的矩阵，R计算过程如下：
>x<-read.table(file.choose(), header=TRUE, row.names=1)#导入txt格式矩阵
>library(robustbase)#调用robustbase
>library(mvoutlier)#调用mvoutlier>covMcd(x, alpha=0.75) #估计稳健统计量(h=0.75)
>cov(x) #估计一般协方差
>apply(x, 2, mean) #估计一般中心值
>res1<-dd.plot(x, quan=0.75) #计算一般马氏距离和稳健马氏距离
>mdc<-res1$md.cla
>write.table(mdc, file="classical_md.txt")#保存一般马氏距离
>mdr<-res1$md.rob
>write.table(mdr, file="robust_md.txt")#保存稳健马氏距离
>res2<-aq.plot(x, quan=0.75, alpha=0.05)#筛选不合格机构
>outliers<-which(res2$outliers==T)
>write.table(outliers,file="outliers.txt")#保存不合格机构
>uni.plot(x, quan=0.75, alpha=0.05)#从各个项目观测整体数据偏离情况
>res1<-covMcd(x, alpha=0.75)
>res2<-arw(x,res1$center,res1$cov,0.05)
>sqrt(res2$cn) #计算判断点位数据异常的临界值，若某机构马氏距离≥该值，则表明该机构能力考核结果偏离中心值较远，结果存在异常。
>q()
计算完成后，各监测机构测定结果偏离中心值的稳健马氏距离(robust_md.txt)、传统马氏距离(classical_md.txt)和基于稳健马氏距离检测得到不合格站点名称(outliers.txt)的分别储存在R软件默认的工作目录下，将25个不合格机构信息汇总进入原始数据集中进一步汇总分析。
3.2 结果与讨论
96家机构8种样品均值、协方差矩阵、稳健中心值和稳健协方差矩阵见表1～表3。比较稳健与非稳健均值和协方差矩阵后发现，由于94、95、96 3个极值的存在(表4)，严重影响了该数据集的均值和协方差。说明通过稳健统计估计出的稳健中心值和稳健协方差矩阵可以排除极值对正确估计数据集统计量的影响，更能体现有机氯能力考核数据的正常分布情况。
表2 96家监测机构稳健协方差矩阵Table 2 The covariance of 96 labs evaluated by robust multi-statistics
表3 96家监测机构非稳健协方差矩阵Table 3 The covariance of 96 labs evaluated by non-robust multi-statistics
表4 数据异常的25家监测机构稳健马氏距离、传统马氏距离、原始浓度和稳健Z比分数判定结果汇总Table 4 The robust mahalanobis distance, classical mehalanobis distance, concentrations and Z-scores of 25 outlier labs
由表4可见，当置信水平a=0.05时，通过“mvoutlier”程序包计算出该次能力考核判定合格与否的临界值为4.3，共筛选出25个稳健马氏距离≥4.3的能力考核结果存在异常的机构(机构名称用序号表示)。94、95、96站点其浓度值显著高于其他站点，为明显的异常值，由于这些极值对中心值和协方差矩阵存在影响，其他机构非稳健马氏距离明显小于稳健马氏距离，表明稳健估计能够排除极端值对于马氏距离计算的影响。
稳健马氏距离越高表明该机构所得结果偏离此次能力考核的中心值越远，可认为其在此次能力考核中的表现越差(表4)。当根据上文介绍的基于稳健马氏距离的算法判断某监测机构能力考核结果为异常值时，说明由于某些分析环节中出现错误，或质量体系中存在的一些问题，导致其与其他合格站点数据分布不一致。通过与Z比分数计算出的α-666、γ-666、p, p′-DDT和o, p′-DDT能力考核结果(满意、有问题或不满意)进行比较(表4)，发现基于稳健马氏距离判定出来能力考核结果异常的监测机构其Z比分数结果都存在着不满意或有问题的项目，说明基于稳健马氏距离的异常值筛选方法能够较好地从多元统计的角度发现能力考核中的不合格单位，其结果与Z比分数方法结果较为一致。同时发现，稳健马氏距离较高的站点，其Z比分数结果存在较多的不满意或有问题，可以考虑在今后的能力考核工作中根据马氏距离这一综合指标来对各个监测机构的整体分析能力做出评估。
4 结论
综上所述，稳健马氏距离方法在兼具了马氏距离的优点的同时，又能较好地排除异常值对其的影响。与传统的针对单轮次、单项目的稳健统计方法(如Z比分数方法、迭代法)相比，该方法可通过马氏距离这个单一指标根据多轮次、多项目的能力考核结果从稳健多元统计的角度对实验室能力或数据质量进行综合定量评价，并能有效筛选出结果存在问题的监测机构，可以在针对多项目、多轮次能力考核结果的综合评价工作中加以试用。该方法在实际工作中也存在着一定的局限性，如某些机构在年内只参加了个别轮次的能力考核，或在某轮次的能力考核中只进行部分样品的分析测试，稳健马氏距离无法对这些机构该年或该轮次的能力考核结果进行综合评估。
[1] 滕曼, 付强, 杨婧, 等. 2011年全国环境监测实验室地表水挥发性有机物检测能力分析[J]. 环境与健康杂志,2013,30(12):1 108-1 109.
TENG M, FU Q, YANG J, et al. Results analysis of proficiency assessment of VOCs monitoring in water [J]. Journal of Environment and Health,2013,30(12):1 108-1 109.
[2] 滕曼, 付强, 吴晓凤, 等. 环境监测实验室水中砷、汞监测能力考核结果评价[J]. 中国环境监测,2014,30(4):183-187.
TENG M, FU Q, WU X F, et al. Results analysis of proficiency assessment of As and Hg monitoring in ground water [J]. Environmental Monitoring of China,2014,30(4):183-187.
[3] 刑小茹, 马小爽, 田文,等. 实验室间比对能力验证中的两种稳健统计技术探讨[J]. 中国环境监测,2011,27(4):4-8.
XING X R, MA X S, TIAN W, et al. Two robust statistic techniques in proficiency testing by interlaboratory comparisons [J]. Environmental Monitoring of China,2011,27(4):4-8.
[4] 吴忠祥. 实验室能力验证中的分割水平检测样品与稳健统计技术[J]. 中国环境监测,2003,19(4):8-10.
WU Z X. Split-level test sample and robust statistics techniques in laboratory proficiency testing [J]. Environmental Monitoring of China, 2003, 19(4): 8-10.
[5] MAHALANOBIS P C. On the generalized distance in statistics [J]. Proceedings of the National Institute of Sciences (Calcutta),1936(2):49-55.
[6] NORSHAHIDA S, ABDUL A J, MOHDT L, et al.
Anomaly detection and sssessment of PM10functional data at several locations in the Klang Valley, Malaysia [J]. Atmospheric Pollution Research,2015(6):365-375.
[7] ROUSSEEUW P J, VAN D K. A fast algorithm for the minimum covariance determinant estimator [J]. Technometrics,1999(41):212-223.
[8] ROUSSEEUW P J, CROUX C, TODOROV V, et al. Robustbase: basic robust statistics [M].2015.
[9] TODOROV V, FILZMOSER P. An object-oriented framework for robust multivariate analysis[J]. Journal of Statistical Software,2009,32(3):1-47.
[10] PISON G, VAN A S, WILLEMS G. Small sample corrections for LTS and MCD [J]. Metrika, 2002(55):111-123.
[11] HUBERT M, ROUSSEEUW P J, VERDONCK T. A deterministic algorithm for robust location and scatter [J].Journal of Computational and Graphical Statistics,2012(21):618-637.
[12] 王斌会, 陈一非. 基于稳健马氏距离的多元异常值检测[J].统计与决策,2005(3):4-6.
WANG B H, CHEN Y F. Multivariate outlier detection based on the robust Mahalanobis distance [J]. Statistics & Decision,2005(3):4-6.
[13] FILZMOSER P, GARRETT R G, REIMANN C. Multivariate outlier detection in exploration geochemistry [J]. Computers & Geosciences,2005(31):579-587.
[14] GERVINI D. A robust and efficient adaptive reweighted estimator of multivariate location and scatter [J]. Journal of Multivariate Analysis,2003(84):116-144.
The Application of Robust Mahalanobis Distance in Proficiency Testing by Interlaboratory Comparisions
SHI Yaolong1,TENG Man1,LI Cheng2,WU Xiaofeng1,CHAI Wenxuan1,YANG Jing1,CHU Baolin1,FU Qiang1
1.State Environmental Protection Key Laboratory of Quality Control in Environmental Monitoring,China National Environmental Monitoring Centre,Bejing 100012,China2.Institute of Science and Technology,Agricultural University of Hebei,Baoding 071001,China
The result of proficiency testing of 4 organo-chlorine pesticide species among 96 enviromental monitoring labs was evaluated by the robust mahalanobis, and 25 labs were identified as outliers by this multi-statistics method. By the comparision between robust mahalanobis distance and Z-score, the outliers identified by robust mahalanobis distance were also indentified as outliers by Z-score. In addition, the fundamental and R implementation of outliers detection by robust mahalanobis distance were described in this article.
robust multi-statistics;robust Mahalanobis distance;interlaboratory comparisions;proficiency testing
2016-03-22;
2016-05-05
国家环保公益性行业科研专项“国家环境监测网环境空气自动监测(PM2.5、O3) 质量保证与质量控制技术体系研究与示范” (201409011)
师耀龙(1988-)，男，河北保定人，博士，工程师。
楚宝临
X830.3
A
1002-6002(2017)02- 0127- 05
10.19316/j.issn.1002-6002.2017.02.20
