最小廣義特征值在多元方差分析中的應用探討

江忠偉

（中國人民銀行南通市中心支行，江蘇 南通 226007）

0 引言

多元方差分析是一元方差分析的推廣，在選擇檢驗統計量方面，通常的做法是：考慮到組內差異是由隨機誤差造成的，組間差異可能是由隨機誤差和系統誤差共同引起的，與一元方差分析的基本思想相同。在一元方差分析中，若各個總體之間沒有顯著差異，則組間離差平方和與組內離差平方和近似相等?？梢宰C明組間離差平和與組內離差平方和的比值服從F分布，給定顯著性水平后，就可以算出臨界值即得出拒絕域。與一元方差分析不同的是：多元統計分析需要將一元方差分析中的組間離差平方和、組內離差平方和推廣為組間離差陣以及組內離差陣。然后基于組間離差陣與組內離差陣的比值構建檢驗統計量，可以證明該統計量為wilks統計量，給定顯著性水平后，就可以算出臨界值即得出拒絕域[1]。另外還有一些其他的檢驗統計量，例如Hotelling跡檢驗統計量[2,3]、Pil?lai-Bartlett準則檢驗統計量（Pillai-Bartlett criterion）[4,5]Roy最大特征值檢驗統計量（Roy’s Largest Root）[6]，具體表達形式見表1。

表1 四種檢驗統計量匯總

通過推導證明，四個檢驗統計量經過適當的變形均轉化成服從F分布的檢驗統計量[7]。在進行多元方差分析時選擇哪個檢驗統計量，是一個很有實際意義的問題。Stevens[7]對上述四個檢驗統計量拒絕原假設能力進行了比較，結果表明：在相同條件下，Roy最大特征值檢驗統計量拒絕能力最強。Olson[8]對上述四種檢驗統計量的檢驗穩健性進行了比較，結果表明：通常，Pillai-Bartlett準則檢驗統計量的穩健性好。

綜上所述，四個檢驗統計量經過適當的變形均可以形成一個以F分布為漸近分布的隨機變量，據此可以在給定的顯著性水平下，設置一個小概率事件：當原假設成立時，檢驗統計量的取值落入構建的小概率事件中，則拒絕原假設。例如，利用wilks檢驗統計量進行檢驗的思路為：首先利用似然比原則導出服從wilks分布的檢驗統計量；由于對wilks檢驗統計量不夠熟悉，通常將wilks檢驗統計量轉換成F檢驗統計量；最后結合一個給定的顯著性水平，就確定了拒絕域，即檢驗法則。其三個檢驗統計量也是按照這種思路：先利用樣本資料導出一個統計量，再將該檢驗統計量轉換成F檢驗統計量，最后結合一個給定的顯著性水平確定拒絕域。有一個很自然的想法是：能否先對樣本資料進行變換，然后再根據變換后的樣本資料構建F檢驗統計量進行方差分析？

1 基本思路

多元方差分析的主要任務是檢驗因子的不同處理（類型變量）對不同處理下得到的樣本觀測值（數值變量）有無顯著影響，即分類自變量對數值因變量有無顯著影響。該模型可以表述為：設分類自變量有K個處理，可以將每個處理看成一個總體，則有總體：

從這K個總體抽取如下樣本：

其中是相互獨立的。

檢驗：

H0:至少有一組i≠j，使得μi≠μj，H1:μ1=…=μK，可以對m個總體中的所有樣品做同一變換即選擇一個p維行向量與所有的樣品進行線性組合，顯然：若H0:至少有一組i≠j，使得μi≠μj成立，則選取任意一個p維行向量，必有H0:至少有一組i≠j，使得≠成立；反之也是如此。

另一方面，由于服從p維多元正態分布的向量的分量的線性組合仍然服從正態分布，所以變換之后的樣品數據仍然服從正態分布。據此可以構建F檢驗統計量進行一元方差分析。但F檢驗統計量的取值是無法確定的，雖然樣本觀測值是已知的，但p維行向量是未知的。如何求出？假設檢驗的目的是尋找證據支持本文的觀點。通常的做法是設置兩個對立事件，然后尋找一個特例拒絕與本文觀點對立的觀點，這樣可以從一定置信水平上認為本文觀點是正確的。因為拒絕一個觀點只需要找到一個特例就行了，而接受一個觀點需要考慮所有的情況（通常是做不到的），因此只需尋找特例來拒絕原假設。利用矩陣的譜分解以及向量的線性表出等知識，可以解出上述F檢驗統計量的最小值以及相對應l′的具體形式。如何利用這個極端值？一般的，對于假設檢驗中的原假設H0，可以認為H0是根據實際問題提出來的，往往是從過去經驗中總結出來的，沒有充分理由不能拒絕它。所以在多元方差分析中，當原假設為：H0:至少有一組i≠j，使得μi≠μj，若原假設為真，即各個總體的均值向量有顯著差異，此時各水平的系統誤差不為零，此時F檢驗統計量（為組間離差平方和與組內離差平方和的比值）會很大。但若由樣本計算出的F檢驗統計量的值小到可以將其看成一個小概率事件，則可以認為原假設是不正確的，此時有較大把握拒絕原假設H0，接受備擇假設H1。

2 依據樣本資料直接構造F檢驗統計量

設分類自變量有K個處理，可以將每個處理看成一個子總體，則有總體：

從這K個子總體抽取如下樣本：

其中是 相 互 獨 立的。按照上文的內容，選擇一個p維向量l′與所有樣品相乘，得出線性組合后的樣本：

樣本數據經過線性組合后均變成了一維數據，由上文可知，檢驗H0:至少有一組i≠j，使得μi≠μj與檢驗H0:至少有一組i≠j，使得≠是等價的。這樣就將多元方差分析轉換為一元方差分析?？梢詷嬙霧檢驗統計量進行一元方差分析。這里存在兩個問題：第一個問題是該樣本數據經歷線性組合之后是否仍然服從正態分布；第二個問題是變換后的樣本數據的組間離差平方和與組內離差平方和是否仍然獨立。接下來分別論證這兩個問題。

2.1 樣本數據線性組合后正態性證明

在一元正態分布中，若Z~N(0 ，1) ，則X=μ+σ Z~N(μ，σ2)。類似的在多元正態分布中，可以類似的定義多元正態分布。設相互獨立且有相同的分布N(0 ，1)，μ為p維常數向量，A為p階常數矩陣，則稱：x=μ+的分布為多元正態分布，記作

可以利用上述定義證明樣本數據進行線性組合后仍然服從正態分布。具體過程如下：

協差陣∑可以分解為：∑=

則可以寫成μj+

則

故得證。

2.2 線性組合后的數據組間離差與組內離差平方和獨立性證明

由上知樣本數據進行線性組合后仍然服從正態分布，可以計算出變換后的樣本數據的總離差平方和SST、組間離差平方SSB和組內離差平方和SSE，經過適當變形之后總離差平方和SST、組間離差平方SSB和組內離差平方和SSE均服從卡方分布，若組間離差平方SSB和組內離差平方和SSE相互獨立，則可以構造出F檢驗統計量進行方差分析。下面證明組間離差平方SSB和組內離差平方和SSE相互獨立。

變換后樣本數據的總離差平方和SST、組間離差平方SSB和組內離差平方和SSE為：

p維行向量l′為一個常數向量，要證明組間離差平方SSB和組內離差平方和SSE之間相互獨立，即證明組間離差陣B和組內離差陣E相互獨立。隨機矩陣的獨立性可以利用的科克朗（Cochran）定理來證明：設X~Nn×p(M，In?Σ )，C和D為n階對稱矩陣，X′CX與X′DX獨立，當且僅當CD=0。另外，若A是投影陣則I-A也是投影陣并且有A(I-A)=0成立。利用科克朗（Cochran）定理以及投影陣的性質，可以很方便地證明組間離差平方SSB和組內離差平方和SSE之間相互獨立。具體證明過程如下：

資料陣Y~Nn×p(M，In?Σ )，其中M的各行是各個子總體的均值向量的轉置按照各個子總體的觀測次數重復排列而成。

可以將總離差陣改寫成：

其中：

可以驗證：

故C為投影陣且rank(C)=n-1；

類似的有：

也可以將組內離差陣E寫成：E=Y′C*Y

其中，C*=diag(C2，…，CK)

顯然C*也是投影陣并且rank(C*)=rank(C1)+rank(C2)+…+rank(CK)=n-K；

組間離差陣B可以改寫成：

其中，

顯然有，C**=C**′， (C**)2=C**，故C**是投影陣并且rank(C**)=trC**=trC+trC*=K-1。

C，C*，C**均為投影陣，并且有C=C*+C**，所以C*C**=0，由科克朗（Cochran）定理知組間離差陣B和組內離差陣E是相互獨立的，故組間離差平方SSB和組內離差平方和SSE之間相互獨立。

綜上所述，本文可以構造出F檢驗統計量：

3 構建檢驗法則

可以將原假設和備擇假設設為：

H0:H0:存在μi≠μj，i≠j；H1:μ1=…=μK

由上文知，可以將原假設和備擇假設改寫成：

H0:存在l′μi≠l′μj，i≠j；H1:l′μ1= … =l′μK

并且這兩組原假設和備擇假設的檢驗結果是等價的。檢驗統計量為:

是一個已知分布的統計量，只需要給出顯著性水平α就可以確定拒絕域的臨界值Fα即得出檢驗法則。

該F檢驗統計量與一般的F統計量有所不同，其中的p維行向量l′事先并不知道，所以無法計算出檢驗統計量的具體數值。但考慮到檢驗的初衷：拒絕與本文觀點對立的觀點，從而證明本文的觀點是正確的。故只需要找到一個特例說明與本文觀點對立的觀點是錯誤的。原假設H0:存在l′μi≠l′μj，i≠j成立時，即系統誤差不為零。所以組間離差平方和與組內離差平方和應該相差很大。若將樣本觀測值帶入檢驗統計量F，計算得出的結果很小，小到可以看成是一個小概率事件，則我們有充分的理由拒絕原假設。所以上述的假設檢驗問題就轉化為已知樣本數據的條件下求解F檢驗統計量的最小值，再與臨界值Fα（下分為數）做出比較。F檢驗統計量的最小值的計算過程如下：

組內離差陣組間離差陣顯然E、B為正定矩陣并且是對稱矩陣，檢驗統計量F可以改寫成：

其中是p階對稱矩陣，故其特征值是實數；又因為為正定矩陣，故其特征值全部大于零。

由矩陣的譜分解知：

其中λ1≥λ2≥…≥λp為B相對于E的廣義特征值，β1，β2，…，βp為B相對于E的廣義特征值λ1≥λ2≥…≥λp所對應的標準化特征向量。β2，…，βp為一組線性無關的p維向量，對β2，…，βp做適當變換后，可以將其看成p維向量空間中的一組標準正交基，該正交基仍然記作β2，…，βp。

由向量的線性表出知：

其中a2，…，ap為常數。

將式（2）、式（3）帶入式（1）得：

當l=βp時，等號成立。

綜上所述，檢驗法則為：當時，有充分理由拒絕原假設，接受備擇假設；當時，不拒絕原假設。

4 利用投影思想進行多元方差分析的優點

傳統的構造檢驗統計量的步驟為：先構造出一個統計量，該統計量的分布是不為我們所熟悉的，為此一般的做法是將該統計量做適當變換使得變換后的統計量的分布漸近服從一個我們熟悉的分布即F分布。這一過程通常計算量較大，并且理論性強不易理解。若直接從投影后的樣本資料出發構建F分布。首先，從推導過程中可以發現，所使用的都是基本的統計知識以及一些線性代數知識，推導過程也十分簡單，可以方便大家理解以及運用該分析方法；其次，隨著計算機的普及以及儲存技術的發展，所研究的數據往往是海量、高維的數據，這是挖掘數據中有價值信息的一個障礙，利用投影思想可以將高維度數據變換成低維度，這種思想的應用無疑帶來了巨大的便利。

5 模擬

為了證實方法的正確性，分兩步進行模擬。

第一步利用R軟件產生9個子總體，每個子總體有20個樣品，這9個子總體的均值向量和協方差陣相同，所有樣品均為5維向量（見表2），其中均值向量和協方差陣是隨機選取的，在此基礎上進行多元方差分析，驗證檢驗結果是否能夠拒絕原假設。

表2 第一個子總體前十個樣品的5維向量

利用計算出的組間離差陣相對于組內離差陣最小廣義特征值為λp=0.0129，故F檢驗統計量的取值為F=，該分位點對應的p=0.02703，非常接近0，因此有充分理由拒絕原假設。

第二步繼續利用R軟件產生9個子總體，每個子總體有20個樣品，與第一步不同的是，這9個子總體的均值向量不相同，所有樣品均為5維向量（見表3），其中均值向量和協差陣是隨機選取的，在此基礎上進行多元方差分析，驗證檢驗結果是否為不能拒絕原假設。

表3 第一個子總體前十個樣品的5維向量

利用計算出的組間離差陣相對于組內離差陣最小廣義特征值為λp=0.0228，故F檢驗統計量的取值為F=，該分位點對應的p=0.136，不是一個非常小的數值，因此沒有充分理由拒絕原假設。

6 結論

本文首先利用投影思想構建的F檢驗統計量在模擬試驗中，當各個子總體均值向量之間不存在差異時，檢驗結果拒絕原假設，接受備擇假設；當各個子總體均值向量之間確實存在差異，F檢驗統計量的取值不能夠拒絕原假設，故可以達到多元方差分析的目的。在假設檢驗過程中，當沒有充分理由拒絕原假設時，這時很多人便認為原假設是正確的。贊同這個觀點的人并沒有考慮原假設錯誤但檢驗統計量取值沒有落入拒絕域中的概率的大小（即納偽的概率），若原假設錯誤時建議統計量取值沒有落入拒絕域中的概率很大，這時認為原假設是正確的顯然是不可信的。此時可以認為檢驗工作并沒有取得實質進展。如何有效克服這個問題有待更進一步的探討。

參考文獻：

[1] Finch H.Comparison of the Performance of Nonparametric and Para?metric MANOVA Test Statistics When Assumptions Are Violated[J].Methodology,2005,1(1).

[2] Kapstad H,Hanestad B R,Langeland N,et al.Cutpoints for Mild,Moderate and Severe Pain in Patients With Osteoarthritis of the Hip or Knee Ready for Joint Replacement Surgery[J].BMC Musculoskele?tal Disorders,2008,9(1).

[3] Hatlen M A,Arora K,Vacic V,et al.Integrative Genetic Analysis of Mouse and Human AML Identifies Cooperating Disease Alleles[J].The Journal of Experimental Medicine,2016,213(1).

[4] Ullah I,Jones B.Regularised Manova for High-Dimensional Data[J].Australian&New Zealand Journal of Statistics,2015,57(3).

[5] Chiani M.Distribution of the Largest Root of a Matrix for Roy’s Test in Multivariate Analysis of Variance[J].Journal of Multivariate Analy?sis,2016,(143).

[6] Haase R F,Ellis M V.Multivariate Analysis of Variance[J].Journal of Counseling Psychology,1987,34(4).

[7] Stevens J P.Power of the Multivariate Analysis of Variance Tests[J].Psychological Bulletin,1980,88(3).

[8] Olson C L.On Choosing a Test Statistic in Multivariate Analysis of Variance[J].Psychological Bulletin,1976,83(4).