基于Excel VBA的能力验证结果核密度图绘制及其应用

黄 欢,黄 宇,王舒乐

(华测检测认证集团股份有限公司 标准物质研究中心,深圳 518000)

在对能力验证结果进行统计分析时,根据GB/T 28043－2019《利用实验室间比对进行能力验证的统计方法》的要求,能力验证提供者应先启用技术和统计知识兼备的人员对参加者提交的结果数据进行直观检查,以确认数据的预期分布,识别离群值或特殊变异的来源。常采用直方图来直观检查结果数据是否呈单峰或对称分布,以及识别异常值。但是,直方图的组距一般对结果数量和切分点的变化反应敏感,难以确认出现簇状数据点时数据的真实分布。核密度图通常被简称为“密度图”,是一种对各点特定分布进行叠加形成的平滑分布曲线,常用于描述数据分布的一般形状。核密度图对于数据分布中可能存在的双峰或不对称分布的识别效果较好,是检查结果数据分布更有效的工具[1]。近年来,常见到尧敦图等能力验证结果统计分析中非必需图形绘制方法的报道,而关于能力验证结果核密度图绘制方法的报道较少[2-5]。核密度曲线计算和绘制过程较为复杂,不易通过常规Excel功能实现,而现有绘制方法均是采用R 语言完成的[1,6],对专业软件以及人员能力依赖性较高。鉴于此,本工作提出了基于Excel VBA 的能力验证结果核密度图绘制方法,以期为能力验证统计人员利用Excel绘制核密度图提供方法参考。

1 核密度图的绘制方法

每个点密度通常由以此点为中心的特定分布(一般为正态分布)来估计,把各点特定分布叠加并进行归一化,即得到“核密度估计”图像[1]。核密度图绘制过程一般包括6个步骤,分别为计算带宽、设置图形区间、设置绘图点数nk、计算各绘图点位置、计算各绘图点密度和生成核密度图。其中,分布标准差即为带宽(用σk表示),构建核密度图所用数据集X通常为参加者提交的系列数据结果(测定值x1,x2,…,xp,p为参加者结果数),也可能为由数据结果转化来的性能统计量值(如测定值的对数值)。

1.1 计算带宽

带宽主要反映核密度曲线的整体平滑程度,带宽越大,各点对核密度曲线形状的贡献越小,曲线越平滑。根据GB/T 28043－2019 规定,常采用Silverman提出的大拇指法则(或Scott法则)和根据国际纯粹与应用化学联合会(IUPAC)指南[7]要求计算带宽σk。依据大拇指法则计算时,先利用参加者提交的数据结果计算稳健标准差s*,进而利用σk＝0.9s*/p0.2(如果采用Scott法则计算,则将公式中系数0.9替换成1.06)计算带宽。依据IUPAC指南计算时,需要分两种情况:基于能力验证目标适用性确定的能力评定标准差σpt利用σk＝0.75σpt计算带宽时,后续采用比分数z或ζ进行能力评定;基于能力验证预先规定的最大允许误差δE利用σk＝0.2δE计算带宽时,后续利用参加者结果与指定值的差值D或参加者结果与指定值的相对差值百分比进行能力评定。

1.2 设置图形区间

利用参加者提供的数据结果选择图形区间,图形区间中最小绘图点为qmin,最大绘图点为qmax,需要分别满足公式(1)和公式(2)的要求。

1.3 设置绘图点数nk

通常,绘图点数nk＝200即可满足要求,若图形区间内存在极端离群值时,应适当增加绘图点数。

1.4 计算各绘图点位置

按照公式(3)计算图形区间内各绘图点位置,其中qi为第i个绘图点。

1.5 计算各绘图点密度

按照公式(4)计算图形区间各绘图点的密度。

式中:hi为第i个绘图点的密度;为标准正态分布密度函数;xj为数据集X 中第j个数据结果。

1.6 生成核密度图

以qi为横坐标,hi为纵坐标绘制平滑曲线,即为核密度图。

2 Excel VBA程序设计

2.1 编程思路

本工作主要采用Excel VBA 编程完成核密度图数据点的计算和图形的绘制。首先根据参加者提交的数据集(x1,x2,…,xp)计算获得p、qmin、qmax和σk,然后进行VBA 编程,基于已知数据qmin、qmax和nk利用公式(3)计算qi系列值,基于已知数据p、σk、nk、qi以及数据集(x1,x2,…,xp)利用公式(2)计算hi系列值,最后绘制核密度图。

2.2 绘制核密度图

本工作采用Microsoft 365 MSO(版本2202)绘制核密度图。首先,新建一个Excel工作簿,调出开发工具(依次点击文件-选项-自定义功能区,勾选开发工具),并开启所有宏(依次点击文件-选项-信任中心-信任中心设置-宏设置,勾选启用所有宏)。其次,把工作表1(sheet1)作为绘制核密度图的工作表,并将参加者提交的数值系列录入单元格C2－C36(以GB/T 28043－2019示例E.6“牛奶中大肠杆菌菌落数量”测试数据为例编写程序,p＝35)。然后,按ALT＋F11进入VBA 编辑器,插入模块并修改模块名称为核密度图,在编程窗口编写如下程序代码并运行。

表1 大米粉中镉含量结果Tab.1 Results of cadmium content in rice flour mg·kg-1

最后,在Excel中对获得的图表进行格式调整,包括增加、删减或修改图表元素(标题、图例和数据标签等)以及设置坐标轴格式等,核密度图绘制过程结束。

3 基于Excel VBA绘制的核密度图的验证

为验证基于Excel VBA 绘制的核密度图的可靠性和准确性,采用Excel VBA 编程绘制GB/T 28043－2019示例E.6大肠杆菌菌落数量结果的核密度图(采用σpt计算带宽),同时与采用R 语言程序绘制的核密度图[1]进行比对。其中,R 语言程序使用density()函数计算各点对应的核密度,带宽bw采用bw＝0.75σpt计算,通过标准正态分布密度函数建立核函数,绘图点数为512 个,说明除了绘图点数,两种方法的重要计算参数和公式基本一致。根据GB/T 28043－2019中10.3.2节规定,在无极端离群值存在下,绘图点数达到200个即可满足要求,因此认为两种方法的绘图点数差异不会对核密度图的绘制造成显著差异。两种方法绘制的核密度图见图1。

图1 基于Excel VBA 和R 语言绘制的核密度图Fig.1 Kernel density plots by Excel VBA and the R programming language

由图1可知,两种方法绘制的曲线的形态以及变化趋势几乎无差异,两条曲线基本重合,说明基于Excel VBA 编程绘制核密度图的方法有效且可靠。

4 核密度图在能力验证结果统计分析中的应用

合理合规的数值结果是确保能力验证工作顺利实施的关键,而核密度图有助于在直观检查环节中确认数据结果的预期分布,识别意外变异或异常值来源[1,8]。例如,从上述GB/T 28043－2019 示例E.6大肠杆菌菌落数量结果的核密度图结果来看,数据呈不对称分布,不符合预期的正态分布,这说明在进行后续能力验证统计分析计算指定值时,不宜采用计算数据算术均值和标准差的经典统计方法或以数据服从正态分布为前提的稳健统计方法,而应采用自助法等非参数估计方法[9-10]。以大米粉中镉含量能力验证数据结果(见表1)为基础,采用Excel VBA 编程绘制核密度图[见图2,带宽采用σpt(0.030 mg·kg－1)计算],所得核密度曲线呈现明显的双峰分布(由于不同方法、污染样本或描述不清楚的操作指令导致的数据结果的混合分布[1])。此时,统计人员应先明确造成数据结果混合分布的原因,再根据差异来源进行数据分组和统计分析。如果无法明确原因和完成分组,则建议使用参照值确定能力验证指定值[11]。

5 结语

本工作使用Excel中VBA 程序处理GB/T 28043－2019中示例E.6的结果数据,得到的核密度图与采用R 语言绘制的完全一致,说明了Excel VBA 绘制方法的可靠性和准确性。提出的核密度图绘制方法只需要Excel软件即可完成,对专业分析工具依赖性低,降低了获取数据结果核密度图的门槛,有利于推广核密度图在数据统计分析中的应用。