基于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軟件即可完成,對專業分析工具依賴性低,降低了獲取數據結果核密度圖的門檻,有利于推廣核密度圖在數據統計分析中的應用。