基于能量距離法的多維分類變量的分布差異檢驗

彭東海，張留偉

(1.中山職業(yè)技術(shù)學(xué)院數(shù)學(xué)教研室，廣東 中山 528404；2.中山大學(xué)測繪學(xué)院，廣東 珠海 519082)

0 引言

分類數(shù)據(jù)又稱定性數(shù)據(jù)，它是反映對象的屬性特征，是非量化的數(shù)據(jù)。例如，社會調(diào)查數(shù)據(jù)中的性別、職業(yè)、學(xué)歷、婚姻狀況等；醫(yī)學(xué)研究數(shù)據(jù)中病人的血型、吸煙與否、采取的治療方案的類別、檢測結(jié)果的陽性與陰性等?？梢?，分類數(shù)據(jù)是我們實際生活中廣泛存在的一類數(shù)據(jù)類型，因而比較兩組分類數(shù)據(jù)的分布差異也就成了我們經(jīng)常遇到的問題。例如，比較兩個不同行業(yè)從業(yè)人員的學(xué)歷結(jié)構(gòu)差異，或比較不同性別職業(yè)分布差異等。這類問題概括來說，就是要比較兩個獨立的分類變量的分布差異，用數(shù)學(xué)語言表述如下：

假設(shè)X與Y是兩個獨立的離散型隨機變量，概率分布列分別為：

P(X=ai)=pi,i=1,2,…,r,

與

P(Y=ai)=qi,i=1,2,…,r.

其中，pi≥0,qi≥0, ∑ipi=1,∑iqi=1.它們對應(yīng)的分布函數(shù)分別為F與G，考慮檢驗問題：

H0∶F=G

(1)

即

H0∶pi=qi，?i=1,2,…,r.

(2)

假設(shè)X1,X2,…,Xm是來自X的簡單隨機樣本，Y1,Y2,…,Yn是來自Y的簡單隨機樣本。對各樣本進行匯總后，可得頻數(shù)分布如表1、2：

表1 X的樣本頻數(shù)分布表

檢驗問題(1)的經(jīng)典方法為卡方檢驗，卡方檢驗將X的樣本X1,X2,…,Xm與Y的樣本Y1,Y2,…,Yn進行合并，并對每個樣本添加標(biāo)簽，標(biāo)注它是來自X的樣本還是來自Y的樣本。因此，檢驗X與Y是否同分布等價于檢驗標(biāo)簽變量與分類值是否相互獨立，從而，我們可將表1的X的樣本頻數(shù)分布表與表2的Y的樣本頻數(shù)分布表合并成一個新的如表3所示的列聯(lián)表：

表2 Y的樣本頻數(shù)分布表

表3 X與Y的樣本的頻數(shù)列聯(lián)表

基于此表，可構(gòu)造卡方檢驗統(tǒng)計量：

眾所周知，卡方檢驗僅適用于大樣本情形，它要求列聯(lián)表中小于5的理論次數(shù)不能多于20%.因此，當(dāng)X與Y取值的類別比較多，而樣本量又不夠大時，容易出現(xiàn)列聯(lián)表中很多單元格頻數(shù)小于5，甚至為0.此時，如果我們?nèi)圆捎每ǚ綑z驗來檢驗兩個分布的差異，則可能會得到錯誤的結(jié)論。在此情形下，可參考文獻[1]關(guān)于卡方檢驗的一些傳統(tǒng)的修正方法。

能量距離 (Energy distance，ED)的概念由Székely 與Rizzo 在文獻 [2]中提出，它是由兩個獨立隨機向量的特征函數(shù)之差的模構(gòu)造的，在選取一個特殊的權(quán)重函數(shù)進行積分之后，它可以表示成距離函數(shù)的期望形式。ED等于0當(dāng)且僅當(dāng)兩個分布相同，因此，ED是度量兩個獨立隨機向量分布差異的一個很好的測度。而樣本的ED只涉及兩點之間的距離的計算，它可以用來檢驗兩個獨立樣本是否同分布，它檢驗的對象可以是任意有限維數(shù)的隨機向量，只須滿足有限一階矩的條件即可，并且檢驗方法具有一致性。Székely 與Rizzo 在文獻[3～8]中分別展示了將ED用于處理諸如獨立樣本同分檢驗、單樣本分布檢驗、聚類分析等一系列經(jīng)典統(tǒng)計問題的結(jié)果。這些研究結(jié)果表明，相比傳統(tǒng)的方法，ED方法計算簡便、適用于更廣泛分布類型的數(shù)據(jù)，且能處理多變量情況。近年來，ED方法引起了國內(nèi)不少學(xué)者的關(guān)注。文獻[9]討論了基于ED推廣的Ward聚類算法研究，文獻[10]對ED做了幾點補充討論，文獻[11]則將ED用于二維列聯(lián)表的獨立性檢驗問題中，得到了一個比傳統(tǒng)卡方檢驗計算更簡便、在小樣本情形下具有更功效的檢驗統(tǒng)計量。

本文受文獻[1]啟發(fā)，基于獨立變量的ED概念，通過引入虛擬向量的方式，給出了兩個分類變量的ED的總體形式與樣本估計，基于該樣本估計，我們得到了一種新的檢驗方法，并采用置換程序計算該方法的檢驗p值。數(shù)值模擬顯示，相比傳統(tǒng)的卡方檢驗法，我們的方法不受列聯(lián)表中單元格的頻數(shù)影響，因此，大部分時候具有更高的功效。最后，我們將方法應(yīng)用于一份涉及兩個多維分類向量的分布差異檢驗的數(shù)據(jù)中。

1 能量距離(ED)方法簡介

我們先回顧一下文獻[1]提出的兩個獨立的隨機變量的能量距離(ED)的定義及主要結(jié)果。

定義1[1]設(shè)X、Y為取值于d的兩個獨立的隨機向量，并且E|X|+E|Y|<∞，則X與Y的能量距離(ED)定義為：

ε(X,Y)∶=2E|X-Y|-E|X-X′|-E|Y-Y′|

其中|.|表示歐氏距離，X′是與X獨立同分布的隨機變量，Y′是與Y獨立同分布的隨機變量。

假設(shè)X1,X2,…,Xm是來自X的簡單隨機樣本，Y1,Y2,…,Yn是來自Y的簡單隨機樣本，則ε(X,Y)的樣本估計為

(7)

文獻[1]證明了εm,n(X,Y)如下的大樣本性質(zhì)：

定理2[1]設(shè)X與Y為取值于d的兩個獨立的隨機變量，X,Y分別有分布函數(shù)F,G，且E|X|+E|Y|<∞.則當(dāng)X與Y的分布相同時，有

(8)

步驟1：由原始樣本{X1,X2,…,Xm}和{Y1,Y2,…,Yn}計算檢驗統(tǒng)計量Tm,n；

步驟2：對聯(lián)合樣本{Z1,Z2,…,Zm,Zm+1,…,Zm+n}={X1,X2,…,Xm,Y1,Y2,…,Yn}隨機地重新排列得

{Zπ(1),Zπ(2),…,Zπ(m),Zπ(m+1),…,Zπ(m+n)}

(9)

2 分類變量的能量距離

定義1中給出的能量距離的概念針對的兩個隨機向量X、Y要求為數(shù)值型變量，若X與Y為分類變量(向量)，我們無法直接用(8)式給出的Tm,n統(tǒng)計量來檢驗X與Y的分布差異。為此，我們首先對X與Y構(gòu)造相應(yīng)的虛擬變量，將X與Y的分布差異檢驗問題等價地轉(zhuǎn)換成對兩個虛擬向量的分布差異的檢驗。具體地，記ei,i=1,2,…,r為r空間中第i個分量元素為1，其他位置元素均為0的單位向量，引入虛擬向量eX與eY如下：

eX=ei,若X=ai,i=1,2,…,r,

及

eY=ei,若Y=ai,i=1,2,…,r.

顯然，P(eX=ei)=P(X=ai)=pi及P(eY=ei)=P(Y=ai)=qi,i=1,2,…,r.因此，檢驗X與Y是否同分布等價于檢驗eX與eY是否同分布。即我們將分類變量X與Y的分布差異檢驗問題等價地轉(zhuǎn)換成了Rr空間中兩個隨機向量的eX與eY的分布差異的檢驗問題。此時，eX與eY的能量距離有表達式：

ε(eX,eY)=2E|eX-eY|-E|eX-eX′|-E|eY-eY′|

其中，eX′是與eX獨立同分布的隨機向量，eY′是與eY獨立同分布的隨機向量。ε(eX,eY)的樣本估計相應(yīng)為：

3 數(shù)值模擬

在這一節(jié)中，我們利用R軟件通過幾個簡單的模擬實驗來比較本文提出的方法與傳統(tǒng)的卡方檢驗的表現(xiàn)?？紤]以下6個例子：

例3：X與Y獨立，均以q，q2，…，qr-1, 1-(q+q2+…+qr-1) 的概率分別取a1，a2，…，ar，其中

q=1/2；

1-(q+q2+…+qr-1)的概率分別取a1,a2,…,ar，其中q=1/2；

顯然，上面例子中，例1至例3為X與Y同分布的情形，而例4至例6為X與Y不同分類的情形。類別個數(shù)考慮r=3, 5, 10, 20, 40五種情況，設(shè)定顯著性水平α=0.05，Tm,n通過499次置換計算p值，樣本容量考慮：1)m=30,n=50; 2)m=50,n=80; 3)m=80,n=100; 4)m=100,n=100; 5)m=100,n=150五種情況。所有的模擬都通過1 000次重復(fù)試驗來近似相應(yīng)的第一類錯誤與功效。表4給出了各種情形下卡方檢驗與Tm,n的模擬結(jié)果。

從表4中我們可以看到，對于例1至例3，Tm,n在各種情形下都能將第一類錯誤很好地控制在名義水平0.05附近，而卡方檢驗在r=20與r=40的情形下，當(dāng)樣本量較小(如m=30,n=50與m=50,n=80)的時候，第一類錯誤遠小于名義水平0.05.這是因為，當(dāng)分類的類別個數(shù)較多，而樣本量較小時，列聯(lián)表中可能會出現(xiàn)很多單元格的期望頻數(shù)小于5，從而導(dǎo)致卡方檢驗統(tǒng)計量近似卡方分布的效果不好，因此第一類錯誤不準(zhǔn)確。很明顯，在樣本量增大的過程中，這種情況稍微得到了改善，卡方檢驗的第一類錯誤逐漸靠近名義水平。對于例4至例6，在類別個數(shù)r≤10，且樣本量較大(如m=80,n=100及以上)時，大部分時候卡方檢驗的功效稍微高于Tm,n，但在樣本量較小(如m=30,n=50)時，或類別個數(shù)較多(如r=40)時，卡方檢驗的功效均不如Tm,n.總的來說，我們的方法Tm,n在小樣本或類別個數(shù)較多時相比卡方檢驗更有優(yōu)勢。

表4 6個例子在不同樣本容量下卡方檢驗與Tm,n檢驗法的第一類錯誤與功效比較(α=0.05)

4 實例分析

本小節(jié)中，我們利用本文提出的方法來分析一份實際數(shù)據(jù)。數(shù)據(jù)來源于R軟件的MASS包中的birthwt數(shù)據(jù)集，該數(shù)據(jù)集記錄了1986年在馬薩諸塞州斯普林菲爾德的Baystate醫(yī)療中心產(chǎn)婦與新生兒的信息，共189條觀測與10個變量。各變量的意義及取值如表5所示。

我們提取其中的race、smoke、ht、ui、ftv這5個分類變量與離散變量來做分析，看看低體重產(chǎn)兒(low=0)的母親與正常產(chǎn)兒(low=1)的母親的這些指標(biāo)之間是否存在差異。下面的圖1至圖5分別給出了兩組產(chǎn)婦之間的這5個指標(biāo)的分布差異對比圖。從這些圖中，我們可以看到，低體重產(chǎn)兒的產(chǎn)婦與正常產(chǎn)兒的產(chǎn)婦之間的這些指標(biāo)的分布有明顯的差異。比如，相比正常產(chǎn)兒的母親，低產(chǎn)兒母親中黑人比例、吸煙比例、有高血壓史的比例、存在子宮刺激的比例以及在懷孕前三個月未看過醫(yī)生的比例都偏高。下面我們采用本文提出的方法與經(jīng)典的卡方檢驗來檢驗兩組樣本的這5個分類變量的聯(lián)合分布的差異。

表5 birthwt數(shù)據(jù)集中各變量名稱及其意義與取值

顯然，兩組的race、smoke、ht、ui與ftv這5個變量的聯(lián)合分布都是取值在3×2×2×2×7=168個類別上的離散分布，理論上的聯(lián)合列聯(lián)表為2×168的二維表，但由于單元格中大部分元素為0，去掉整行或整列為0的行或列后，只剩下2×46的列聯(lián)表。此時，我們采用卡方檢驗對此列聯(lián)表做獨立性檢驗，可得卡方統(tǒng)計量的值為54.688, 相應(yīng)的p值為0.1527.這說明，卡方檢驗結(jié)果顯示，低體重產(chǎn)兒的產(chǎn)婦與正常產(chǎn)兒的產(chǎn)婦之間的這5個分類指標(biāo)的聯(lián)合分布沒有明顯的差異。接下來，我們采用本文提出的Tm,n重新對兩組的5個分類指標(biāo)的聯(lián)合分布做差異性檢驗，可得統(tǒng)計量的值為 28 447.99，在499次置換樣本的情況下，得到相應(yīng)的p值為0.024.因此，我們的方法說明，低體重產(chǎn)兒的產(chǎn)婦與正常產(chǎn)兒的產(chǎn)婦之間的這5個分類指標(biāo)的聯(lián)合分布是有明顯的差異的，這跟圖1至圖5展示的結(jié)果是一致的。

圖1 低體重產(chǎn)兒與正常產(chǎn)兒的產(chǎn)婦的種族分布對比圖

圖2 低體重產(chǎn)兒與正常產(chǎn)兒的產(chǎn)婦是否吸煙的分布對比圖

圖5 低體重產(chǎn)兒與正常產(chǎn)兒的產(chǎn)婦是在懷孕前三個月看醫(yī)生次數(shù)的分布對比圖

5 總結(jié)

本文基于獨立隨機向量的能量距離(ED)的概念給出了分類變量(向量)的能量距離(ED)的定義，由此給出了一種新的檢驗兩個分類變量(向量)的分布差異的檢驗統(tǒng)計量。該檢驗方法的檢驗統(tǒng)計量只涉及兩組樣本頻數(shù)的差距，因此計算非常簡便，且適用于任意有限維的分類的隨機向量。數(shù)值模擬結(jié)果說明，在小樣本或類別變量較多的情形下，該方法比經(jīng)典的卡方檢驗法能更好地控制第一類錯誤，同時具有更高的功效。