黃 歡,黃 宇,王舒樂(lè)

(華測(cè)檢測(cè)認(rèn)證集團(tuán)股份有限公司 標(biāo)準(zhǔn)物質(zhì)研究中心,深圳 518000)

在對(duì)能力驗(yàn)證結(jié)果進(jìn)行統(tǒng)計(jì)分析時(shí),根據(jù)GB/T 28043－2019《利用實(shí)驗(yàn)室間比對(duì)進(jìn)行能力驗(yàn)證的統(tǒng)計(jì)方法》的要求,能力驗(yàn)證提供者應(yīng)先啟用技術(shù)和統(tǒng)計(jì)知識(shí)兼?zhèn)涞娜藛T對(duì)參加者提交的結(jié)果數(shù)據(jù)進(jìn)行直觀檢查,以確認(rèn)數(shù)據(jù)的預(yù)期分布,識(shí)別離群值或特殊變異的來(lái)源。常采用直方圖來(lái)直觀檢查結(jié)果數(shù)據(jù)是否呈單峰或?qū)ΨQ分布,以及識(shí)別異常值。但是,直方圖的組距一般對(duì)結(jié)果數(shù)量和切分點(diǎn)的變化反應(yīng)敏感,難以確認(rèn)出現(xiàn)簇狀數(shù)據(jù)點(diǎn)時(shí)數(shù)據(jù)的真實(shí)分布。核密度圖通常被簡(jiǎn)稱為“密度圖”,是一種對(duì)各點(diǎn)特定分布進(jìn)行疊加形成的平滑分布曲線,常用于描述數(shù)據(jù)分布的一般形狀。核密度圖對(duì)于數(shù)據(jù)分布中可能存在的雙峰或不對(duì)稱分布的識(shí)別效果較好,是檢查結(jié)果數(shù)據(jù)分布更有效的工具[1]。近年來(lái),常見(jiàn)到堯敦圖等能力驗(yàn)證結(jié)果統(tǒng)計(jì)分析中非必需圖形繪制方法的報(bào)道,而關(guān)于能力驗(yàn)證結(jié)果核密度圖繪制方法的報(bào)道較少[2-5]。核密度曲線計(jì)算和繪制過(guò)程較為復(fù)雜,不易通過(guò)常規(guī)Excel功能實(shí)現(xiàn),而現(xiàn)有繪制方法均是采用R 語(yǔ)言完成的[1,6],對(duì)專業(yè)軟件以及人員能力依賴性較高。鑒于此,本工作提出了基于Excel VBA 的能力驗(yàn)證結(jié)果核密度圖繪制方法,以期為能力驗(yàn)證統(tǒng)計(jì)人員利用Excel繪制核密度圖提供方法參考。

1 核密度圖的繪制方法

每個(gè)點(diǎn)密度通常由以此點(diǎn)為中心的特定分布(一般為正態(tài)分布)來(lái)估計(jì),把各點(diǎn)特定分布疊加并進(jìn)行歸一化,即得到“核密度估計(jì)”圖像[1]。核密度圖繪制過(guò)程一般包括6個(gè)步驟,分別為計(jì)算帶寬、設(shè)置圖形區(qū)間、設(shè)置繪圖點(diǎn)數(shù)nk、計(jì)算各繪圖點(diǎn)位置、計(jì)算各繪圖點(diǎn)密度和生成核密度圖。其中,分布標(biāo)準(zhǔn)差即為帶寬(用σk表示),構(gòu)建核密度圖所用數(shù)據(jù)集X通常為參加者提交的系列數(shù)據(jù)結(jié)果(測(cè)定值x1,x2,…,xp,p為參加者結(jié)果數(shù)),也可能為由數(shù)據(jù)結(jié)果轉(zhuǎn)化來(lái)的性能統(tǒng)計(jì)量值(如測(cè)定值的對(duì)數(shù)值)。

1.1 計(jì)算帶寬

帶寬主要反映核密度曲線的整體平滑程度,帶寬越大,各點(diǎn)對(duì)核密度曲線形狀的貢獻(xiàn)越小,曲線越平滑。根據(jù)GB/T 28043－2019 規(guī)定,常采用Silverman提出的大拇指法則(或Scott法則)和根據(jù)國(guó)際純粹與應(yīng)用化學(xué)聯(lián)合會(huì)(IUPAC)指南[7]要求計(jì)算帶寬σk。依據(jù)大拇指法則計(jì)算時(shí),先利用參加者提交的數(shù)據(jù)結(jié)果計(jì)算穩(wěn)健標(biāo)準(zhǔn)差s*,進(jìn)而利用σk＝0.9s*/p0.2(如果采用Scott法則計(jì)算,則將公式中系數(shù)0.9替換成1.06)計(jì)算帶寬。依據(jù)IUPAC指南計(jì)算時(shí),需要分兩種情況:基于能力驗(yàn)證目標(biāo)適用性確定的能力評(píng)定標(biāo)準(zhǔn)差σpt利用σk＝0.75σpt計(jì)算帶寬時(shí),后續(xù)采用比分?jǐn)?shù)z或ζ進(jìn)行能力評(píng)定;基于能力驗(yàn)證預(yù)先規(guī)定的最大允許誤差δE利用σk＝0.2δE計(jì)算帶寬時(shí),后續(xù)利用參加者結(jié)果與指定值的差值D或參加者結(jié)果與指定值的相對(duì)差值百分比進(jìn)行能力評(píng)定。

1.2 設(shè)置圖形區(qū)間

利用參加者提供的數(shù)據(jù)結(jié)果選擇圖形區(qū)間,圖形區(qū)間中最小繪圖點(diǎn)為qmin,最大繪圖點(diǎn)為qmax,需要分別滿足公式(1)和公式(2)的要求。

1.3 設(shè)置繪圖點(diǎn)數(shù)nk

通常,繪圖點(diǎn)數(shù)nk＝200即可滿足要求,若圖形區(qū)間內(nèi)存在極端離群值時(shí),應(yīng)適當(dāng)增加繪圖點(diǎn)數(shù)。

1.4 計(jì)算各繪圖點(diǎn)位置

按照公式(3)計(jì)算圖形區(qū)間內(nèi)各繪圖點(diǎn)位置,其中qi為第i個(gè)繪圖點(diǎn)。

1.5 計(jì)算各繪圖點(diǎn)密度

按照公式(4)計(jì)算圖形區(qū)間各繪圖點(diǎn)的密度。

式中:hi為第i個(gè)繪圖點(diǎn)的密度;為標(biāo)準(zhǔn)正態(tài)分布密度函數(shù);xj為數(shù)據(jù)集X 中第j個(gè)數(shù)據(jù)結(jié)果。

1.6 生成核密度圖

以qi為橫坐標(biāo),hi為縱坐標(biāo)繪制平滑曲線,即為核密度圖。

2 Excel VBA程序設(shè)計(jì)

2.1 編程思路

本工作主要采用Excel VBA 編程完成核密度圖數(shù)據(jù)點(diǎn)的計(jì)算和圖形的繪制。首先根據(jù)參加者提交的數(shù)據(jù)集(x1,x2,…,xp)計(jì)算獲得p、qmin、qmax和σk,然后進(jìn)行VBA 編程,基于已知數(shù)據(jù)qmin、qmax和nk利用公式(3)計(jì)算qi系列值,基于已知數(shù)據(jù)p、σk、nk、qi以及數(shù)據(jù)集(x1,x2,…,xp)利用公式(2)計(jì)算hi系列值,最后繪制核密度圖。

2.2 繪制核密度圖

本工作采用Microsoft 365 MSO(版本2202)繪制核密度圖。首先,新建一個(gè)Excel工作簿,調(diào)出開(kāi)發(fā)工具(依次點(diǎn)擊文件-選項(xiàng)-自定義功能區(qū),勾選開(kāi)發(fā)工具),并開(kāi)啟所有宏(依次點(diǎn)擊文件-選項(xiàng)-信任中心-信任中心設(shè)置-宏設(shè)置,勾選啟用所有宏)。其次,把工作表1(sheet1)作為繪制核密度圖的工作表,并將參加者提交的數(shù)值系列錄入單元格C2－C36(以GB/T 28043－2019示例E.6“牛奶中大腸桿菌菌落數(shù)量”測(cè)試數(shù)據(jù)為例編寫(xiě)程序,p＝35)。然后,按ALT＋F11進(jìn)入VBA 編輯器,插入模塊并修改模塊名稱為核密度圖,在編程窗口編寫(xiě)如下程序代碼并運(yùn)行。

表1 大米粉中鎘含量結(jié)果Tab.1 Results of cadmium content in rice flour mg·kg-1

最后,在Excel中對(duì)獲得的圖表進(jìn)行格式調(diào)整,包括增加、刪減或修改圖表元素(標(biāo)題、圖例和數(shù)據(jù)標(biāo)簽等)以及設(shè)置坐標(biāo)軸格式等,核密度圖繪制過(guò)程結(jié)束。

3 基于Excel VBA繪制的核密度圖的驗(yàn)證

為驗(yàn)證基于Excel VBA 繪制的核密度圖的可靠性和準(zhǔn)確性,采用Excel VBA 編程繪制GB/T 28043－2019示例E.6大腸桿菌菌落數(shù)量結(jié)果的核密度圖(采用σpt計(jì)算帶寬),同時(shí)與采用R 語(yǔ)言程序繪制的核密度圖[1]進(jìn)行比對(duì)。其中,R 語(yǔ)言程序使用density()函數(shù)計(jì)算各點(diǎn)對(duì)應(yīng)的核密度,帶寬bw采用bw＝0.75σpt計(jì)算,通過(guò)標(biāo)準(zhǔn)正態(tài)分布密度函數(shù)建立核函數(shù),繪圖點(diǎn)數(shù)為512 個(gè),說(shuō)明除了繪圖點(diǎn)數(shù),兩種方法的重要計(jì)算參數(shù)和公式基本一致。根據(jù)GB/T 28043－2019中10.3.2節(jié)規(guī)定,在無(wú)極端離群值存在下,繪圖點(diǎn)數(shù)達(dá)到200個(gè)即可滿足要求,因此認(rèn)為兩種方法的繪圖點(diǎn)數(shù)差異不會(huì)對(duì)核密度圖的繪制造成顯著差異。兩種方法繪制的核密度圖見(jiàn)圖1。

圖1 基于Excel VBA 和R 語(yǔ)言繪制的核密度圖Fig.1 Kernel density plots by Excel VBA and the R programming language

由圖1可知,兩種方法繪制的曲線的形態(tài)以及變化趨勢(shì)幾乎無(wú)差異,兩條曲線基本重合,說(shuō)明基于Excel VBA 編程繪制核密度圖的方法有效且可靠。

4 核密度圖在能力驗(yàn)證結(jié)果統(tǒng)計(jì)分析中的應(yīng)用

合理合規(guī)的數(shù)值結(jié)果是確保能力驗(yàn)證工作順利實(shí)施的關(guān)鍵,而核密度圖有助于在直觀檢查環(huán)節(jié)中確認(rèn)數(shù)據(jù)結(jié)果的預(yù)期分布,識(shí)別意外變異或異常值來(lái)源[1,8]。例如,從上述GB/T 28043－2019 示例E.6大腸桿菌菌落數(shù)量結(jié)果的核密度圖結(jié)果來(lái)看,數(shù)據(jù)呈不對(duì)稱分布,不符合預(yù)期的正態(tài)分布,這說(shuō)明在進(jìn)行后續(xù)能力驗(yàn)證統(tǒng)計(jì)分析計(jì)算指定值時(shí),不宜采用計(jì)算數(shù)據(jù)算術(shù)均值和標(biāo)準(zhǔn)差的經(jīng)典統(tǒng)計(jì)方法或以數(shù)據(jù)服從正態(tài)分布為前提的穩(wěn)健統(tǒng)計(jì)方法,而應(yīng)采用自助法等非參數(shù)估計(jì)方法[9-10]。以大米粉中鎘含量能力驗(yàn)證數(shù)據(jù)結(jié)果(見(jiàn)表1)為基礎(chǔ),采用Excel VBA 編程繪制核密度圖[見(jiàn)圖2,帶寬采用σpt(0.030 mg·kg－1)計(jì)算],所得核密度曲線呈現(xiàn)明顯的雙峰分布(由于不同方法、污染樣本或描述不清楚的操作指令導(dǎo)致的數(shù)據(jù)結(jié)果的混合分布[1])。此時(shí),統(tǒng)計(jì)人員應(yīng)先明確造成數(shù)據(jù)結(jié)果混合分布的原因,再根據(jù)差異來(lái)源進(jìn)行數(shù)據(jù)分組和統(tǒng)計(jì)分析。如果無(wú)法明確原因和完成分組,則建議使用參照值確定能力驗(yàn)證指定值[11]。

5 結(jié)語(yǔ)

本工作使用Excel中VBA 程序處理GB/T 28043－2019中示例E.6的結(jié)果數(shù)據(jù),得到的核密度圖與采用R 語(yǔ)言繪制的完全一致,說(shuō)明了Excel VBA 繪制方法的可靠性和準(zhǔn)確性。提出的核密度圖繪制方法只需要Excel軟件即可完成,對(duì)專業(yè)分析工具依賴性低,降低了獲取數(shù)據(jù)結(jié)果核密度圖的門(mén)檻,有利于推廣核密度圖在數(shù)據(jù)統(tǒng)計(jì)分析中的應(yīng)用。