山東大學(xué)公共衛(wèi)生學(xué)院(250012) 高青松 薛付忠

自從Karl Pearson〔1〕提出主成分分析(PCA)以來(lái)，這種多元統(tǒng)計(jì)分析方法已經(jīng)被廣泛應(yīng)用于各個(gè)領(lǐng)域，如圖像處理、故障診斷、綜合分析等。主成分分析的基本思想是通過(guò)降維(線(xiàn)性變換)的思想，在損失很少信息的前提下把原始數(shù)據(jù)中的多個(gè)變量化為少數(shù)幾個(gè)綜合指標(biāo)(即主成分)，從而達(dá)到簡(jiǎn)化系統(tǒng)結(jié)構(gòu)，抓住問(wèn)題本質(zhì)的目的。

傳統(tǒng)的線(xiàn)性主成分分析可以有效地處理變量間的線(xiàn)性關(guān)系，但是現(xiàn)實(shí)中有很多數(shù)據(jù)是非線(xiàn)性的，而線(xiàn)性主成分分析并不能檢測(cè)出數(shù)據(jù)之間的非線(xiàn)性結(jié)構(gòu)〔2〕。另外，線(xiàn)性主成分分析也不能準(zhǔn)確地代表不滿(mǎn)足高斯分布的數(shù)據(jù)〔3〕。因此，很多文獻(xiàn)提出用非線(xiàn)性的主成分分析方法來(lái)提取數(shù)據(jù)之間的非線(xiàn)性特征，主要包括核主成分分析(KPCA)〔4〕和基于神經(jīng)網(wǎng)絡(luò)(neural network)〔2，5〕的主成分分析等。尤其是核主成分分析，不僅可以提取數(shù)據(jù)的非線(xiàn)性特征，而且其非線(xiàn)性主成分計(jì)算簡(jiǎn)單(求解特征值問(wèn)題)，意義明確(特征空間中的線(xiàn)性主成分)，因此在很多領(lǐng)域得到廣泛應(yīng)用。

在進(jìn)行核主成分分析時(shí)，如果選擇線(xiàn)性核函數(shù)，則核主成分分析等價(jià)于線(xiàn)性主成分分析，即在選擇正確的核函數(shù)的前提下，核主成分分析可用于分析線(xiàn)性數(shù)據(jù)。然而，到目前為止尚未提出一種可以有效地選擇核函數(shù)及其參數(shù)的方法，如果盲目地對(duì)線(xiàn)性數(shù)據(jù)進(jìn)行核主成分分析，不僅核函數(shù)的選擇及其參數(shù)的設(shè)置費(fèi)時(shí)費(fèi)力，而且會(huì)大大增加計(jì)算復(fù)雜度;如果選擇非線(xiàn)性的核函數(shù)，無(wú)法代表原始數(shù)據(jù)的非線(xiàn)性結(jié)構(gòu)，從而降低其效能。所以，在進(jìn)行主成分分析(線(xiàn)性或非線(xiàn)性)之前，必須對(duì)數(shù)據(jù)是線(xiàn)性結(jié)構(gòu)還是非線(xiàn)性結(jié)構(gòu)進(jìn)行檢驗(yàn)。

Uwe Kruger等在過(guò)程控制(process control)領(lǐng)域提出了一個(gè)非線(xiàn)性檢驗(yàn)的度量〔6〕，但此方法需要對(duì)數(shù)據(jù)有一定的先驗(yàn)知識(shí)，或者可以根據(jù)其前幾個(gè)主成分的散點(diǎn)圖來(lái)進(jìn)行分區(qū)，然而我們?cè)谶M(jìn)行數(shù)據(jù)處理時(shí)，絕大部分情況下對(duì)數(shù)據(jù)的了解還處于探索階段，不大可能有明確的先驗(yàn)知識(shí)，而且大部分?jǐn)?shù)據(jù)很難根據(jù)其主成分的散點(diǎn)圖來(lái)進(jìn)行分區(qū)，除非所選擇的樣本有很強(qiáng)的層次結(jié)構(gòu)?；诖?，本文在上述非線(xiàn)性檢驗(yàn)度量的基礎(chǔ)上，提出用主成分聚類(lèi)分析法來(lái)進(jìn)行分區(qū)，減少了主觀(guān)因素對(duì)分區(qū)的影響，并可有效地處理主成分的散點(diǎn)圖無(wú)法分區(qū)的情形，從而得到一個(gè)實(shí)用性更強(qiáng)的非線(xiàn)性檢驗(yàn)度量。

原理與方法

1.原理 首先將原始數(shù)據(jù)中的所有樣本分成若干區(qū)域，然后任意選擇一個(gè)區(qū)域，計(jì)算該區(qū)域內(nèi)變量之間的相關(guān)系數(shù)矩陣，對(duì)于相關(guān)系數(shù)矩陣中的每一個(gè)元素，根據(jù)變量的均值和方差計(jì)算出其置信限;并根據(jù)此置信限估計(jì)出主成分分析的殘差(也就是丟棄的特征值之和)的最大值和最小值，或簡(jiǎn)稱(chēng)為精確邊界;最后計(jì)算每一個(gè)區(qū)域按照選定區(qū)域的均值和方差進(jìn)行標(biāo)準(zhǔn)化之后，進(jìn)行主成分分析時(shí)的殘差。如果所有的殘差都落在估計(jì)出的精確邊界以?xún)?nèi)，則說(shuō)明該數(shù)據(jù)變量之間是線(xiàn)性關(guān)系;當(dāng)至少有一個(gè)殘差落在精確邊界以外時(shí)，就判斷變量之間是非線(xiàn)性關(guān)系。

由于最初是任意選擇一個(gè)區(qū)域來(lái)計(jì)算精確邊界，為了避免可能的偏倚，本文采用交叉驗(yàn)證的原理，對(duì)每個(gè)區(qū)域都估計(jì)其精確邊界并計(jì)算殘差，然后按照上述方法來(lái)判斷變量之間的關(guān)系。

2.方法 具體的檢驗(yàn)過(guò)程分為以下四步:

第一步:定義不連續(xù)區(qū)域 首先對(duì)原始數(shù)據(jù)進(jìn)行主成分分析，選出前幾個(gè)累積貢獻(xiàn)率達(dá)到80%的主成分，然后對(duì)選出的主成分進(jìn)行K均值聚類(lèi)(K-means clustering)。

第二步:計(jì)算相關(guān)系數(shù)矩陣中每個(gè)元素的置信限給定原始數(shù)據(jù)集有N個(gè)變量，則其相關(guān)系數(shù)矩陣為

假設(shè)將原始數(shù)據(jù)集分成m個(gè)區(qū)域，則第h(h=1，…，m)個(gè)區(qū)域的相關(guān)系數(shù)矩陣同樣可根據(jù)該區(qū)域中的數(shù)據(jù)得到，記為R(h)。根據(jù)每個(gè)變量的均值和方差可計(jì)算出R(h)中每個(gè)元素的置信限，用矩陣表示如下

第三步:確定精確邊界 根據(jù)式(1)中相關(guān)系數(shù)矩陣的置信限可直接估計(jì)精確邊界。由主成分分析中殘差矩陣的弗羅貝尼烏斯范數(shù)(frobenius norm)可以證明，殘差σ等價(jià)于主成分模型中丟棄的特征值之和:

其中σj代表第j個(gè)變量的殘差，λk代表相關(guān)系數(shù)矩陣R(h)的第k個(gè)特征值，n(h)是該區(qū)域中的樣本數(shù)，n是保留的主成分的個(gè)數(shù)。

由于特征值 λn+1，…，λN依賴(lài)于相關(guān)系數(shù)矩陣R(h)的元素，所以精確邊界的估計(jì)可轉(zhuǎn)化為求解下面的最優(yōu)化問(wèn)題:

其中，ΔRmax和ΔRmin中的元素分別是R中非對(duì)角元素的波動(dòng)，用于決定最大值λkmax和最小值λkmin，然后根據(jù)公式(2)可得第h個(gè)區(qū)域的精確邊界為

本文對(duì)上述最優(yōu)化問(wèn)題采用粒子群優(yōu)化算法來(lái)求解(particle swarm optimization，PSO)〔7〕。

第四步:定義非線(xiàn)性度量對(duì)每一個(gè)h(h=1，2，…，m)，按照前面三步得到其精確邊界σhmax，σhmin和殘差以后，對(duì)其進(jìn)行比較，如果所有的殘差都落在精確邊界之內(nèi)，則該數(shù)據(jù)為線(xiàn)性數(shù)據(jù);當(dāng)至少有一個(gè)落在精確邊界之外，則判斷為非線(xiàn)性數(shù)據(jù)。

實(shí)例分析

為驗(yàn)證該非線(xiàn)性檢驗(yàn)度量的有效性，我們首先用該度量分析了文獻(xiàn)〔6〕中的兩份模擬數(shù)據(jù)，包括一個(gè)線(xiàn)性的例子和一個(gè)非線(xiàn)性的例子。接下來(lái)，又用該度量對(duì)北美風(fēng)濕性關(guān)節(jié)炎全基因組關(guān)聯(lián)分析的(SNP)數(shù)據(jù)〔8〕的基因PTPN22基因區(qū)域內(nèi)9個(gè)單核苷酸多態(tài)性(SNPs)之間的線(xiàn)性或非線(xiàn)性特征進(jìn)行判斷，為我們進(jìn)一步選擇核主成分分析構(gòu)建基因組區(qū)域化關(guān)聯(lián)分析檢驗(yàn)統(tǒng)計(jì)量提供依據(jù)。

1．線(xiàn)性例子 首先產(chǎn)生1000個(gè)服從均勻分布的樣本，范圍在－4到4之間，作為隨機(jī)變量的值，然后產(chǎn)生兩個(gè)相互獨(dú)立且均服從均值為0，方差為0.005的正態(tài)分布的隨機(jī)變量和，則線(xiàn)性情況下可根據(jù)如下公式產(chǎn)生兩個(gè)變量z1和z2:

z1=e1+t，z2=t+e2

圖1 線(xiàn)性例子的散點(diǎn)圖及其分區(qū)結(jié)果

圖1給出了線(xiàn)性情況下變量z1和z2的散點(diǎn)圖，并給出了主成分聚類(lèi)分析法的分區(qū)結(jié)果。表1給出了每個(gè)區(qū)域中相關(guān)系數(shù)的值及其置信限，并計(jì)算了每個(gè)區(qū)域的精確邊界。表2給出了每個(gè)區(qū)域中主成分分析的殘差，也就是相關(guān)系數(shù)矩陣的第二個(gè)特征值。其中對(duì)角元素精確邊界是根據(jù)該區(qū)域的數(shù)據(jù)估計(jì)得到的。通過(guò)對(duì)表1和表2的結(jié)果進(jìn)行比較可得，沒(méi)有落在精確邊界以外的數(shù)值，這跟預(yù)期結(jié)果一樣，即我們假設(shè)的變量z1和z2之間的關(guān)系是線(xiàn)性的。

表1 線(xiàn)性例子中每個(gè)區(qū)域內(nèi)的相關(guān)系數(shù)及其置信限和精確邊界

表2 每個(gè)區(qū)域中主成分分析的殘差

2．非線(xiàn)性例子 同樣，首先產(chǎn)生1000個(gè)服從均勻分布的樣本，范圍在－4到4之間，作為隨機(jī)變量t的值，然后產(chǎn)生兩個(gè)相互獨(dú)立且均服從均值為0，方差為0.005的正態(tài)分布的隨機(jī)變量e1和e2，非線(xiàn)性情況下變量z1和z2可根據(jù)下面的公式產(chǎn)生:

圖2 非線(xiàn)性例子的散點(diǎn)圖及其分區(qū)結(jié)果

圖2給出了非線(xiàn)性情況下變量z1和z2的散點(diǎn)圖，同樣給出了主成分聚類(lèi)分析法的分區(qū)結(jié)果。表3給出了每個(gè)區(qū)域中相關(guān)系數(shù)的值及其置信限，并計(jì)算了每個(gè)區(qū)域的精確邊界。表4給出了每個(gè)區(qū)域中主成分分析的殘差。與線(xiàn)性情況不同的是，每個(gè)區(qū)域都有落在精確邊界以外的數(shù)值，在表4中用粗體表示，這是因?yàn)槲覀兗俣ǖ淖兞縵1和z2之間的關(guān)系是非線(xiàn)性的。

3.實(shí)際數(shù)據(jù) 在全基因組關(guān)聯(lián)分析中，為了減少多重檢驗(yàn)造成的假陽(yáng)性、數(shù)據(jù)共線(xiàn)性特征等所致的困擾，需要構(gòu)建基因組區(qū)域化關(guān)聯(lián)分析檢驗(yàn)統(tǒng)計(jì)量，本課題組已經(jīng)構(gòu)建了一種基于線(xiàn)性主成分bootstrap可信區(qū)間法的檢驗(yàn)統(tǒng)計(jì)量〔8〕。然而，基因組SNP之間的線(xiàn)性關(guān)系的假設(shè)往往不成立，需要進(jìn)一步采用非線(xiàn)性主成分分析(核主成分分析)構(gòu)建統(tǒng)計(jì)效能更高的基因組區(qū)域化關(guān)聯(lián)分析統(tǒng)計(jì)量。為了驗(yàn)證本文所提出的非線(xiàn)性檢驗(yàn)度量的有效性，我們對(duì)北美風(fēng)濕性關(guān)節(jié)炎全基因組關(guān)聯(lián)分析的(SNP)數(shù)據(jù)〔9〕的基因PTPN22基因區(qū)域內(nèi)9個(gè)單核苷酸多態(tài)性(SNPs)之間的線(xiàn)性或非線(xiàn)性特征進(jìn)行判斷，從而為進(jìn)一步選擇核主成分分析構(gòu)建基因組區(qū)域化關(guān)聯(lián)分析檢驗(yàn)統(tǒng)計(jì)量提供依據(jù)。該數(shù)據(jù)包括868個(gè)病例和1 194個(gè)對(duì)照。PTPN22基因區(qū)域中的9個(gè)SNPs(rs114157960-rs114215857)的連鎖不平衡區(qū)域結(jié)構(gòu)(LD block)見(jiàn)圖3。

表3 非線(xiàn)性例子中的相關(guān)系數(shù)及其置信限和精確邊界

表4 每個(gè)區(qū)域中主成分分析的殘差

圖3 PTPN22基因中9個(gè)SNP的LD圖

對(duì)該數(shù)據(jù)進(jìn)行主成分分析，其前三個(gè)主成分的累積貢獻(xiàn)率為80.5%，故將其他6個(gè)主成分舍去。圖4(A)給出了前三個(gè)主成分的散點(diǎn)圖，從該圖中，我們無(wú)法得到任何分區(qū)結(jié)構(gòu)，而根據(jù)主成分聚類(lèi)分析法可以很清晰地將原始數(shù)據(jù)分成四個(gè)區(qū)域(圖4(B))。

表5給出了這9個(gè)SNP的精確邊界及特征值。精確邊界與預(yù)期的一樣，這些值落在區(qū)域之內(nèi)。然而，非對(duì)角元素至少一個(gè)值落在精確邊界之外(用粗體表示)，這就意味著這些SNPs之間是非線(xiàn)性關(guān)系，從而表明在構(gòu)建基因組全區(qū)域化關(guān)聯(lián)分析檢驗(yàn)統(tǒng)計(jì)量時(shí)采用非線(xiàn)性主成分分析(核主成分分析)的必要性。

討 論

圖4 PTPN22基因中9個(gè)SNP前三個(gè)主成分的散點(diǎn)圖(A)及其分區(qū)結(jié)果(B)

表5 PTPN22基因中9個(gè)SNP的非線(xiàn)性檢驗(yàn)結(jié)果

非線(xiàn)性的核主成分分析由于其強(qiáng)大的非線(xiàn)性特征提取能力以及其原理簡(jiǎn)單而被廣泛應(yīng)用于圖像識(shí)別〔10〕、特征提取〔11〕和綜合評(píng)價(jià)〔12〕等多個(gè)領(lǐng)域。然而，核主成分分析的核函數(shù)選擇及其參數(shù)調(diào)試問(wèn)題仍未能得到很好的解決。所以，在進(jìn)行核主成分分析之前，應(yīng)該判斷是否有必要對(duì)數(shù)據(jù)進(jìn)行這種復(fù)雜的分析，也就是說(shuō)，應(yīng)判斷數(shù)據(jù)是否為非線(xiàn)性數(shù)據(jù)。目前，對(duì)數(shù)據(jù)進(jìn)行非線(xiàn)性主成分分析，絕大部分研究都只是簡(jiǎn)單地將數(shù)據(jù)的非線(xiàn)性作為一個(gè)“前提條件(a prior)”〔6〕。為了研究數(shù)據(jù)的結(jié)構(gòu)，Uwe Kruger等人提出了一個(gè)非線(xiàn)性檢驗(yàn)的度量，但該度量并不能很好地應(yīng)用于其他領(lǐng)域，比如說(shuō)全基因組關(guān)聯(lián)研究(GWAS)得到的單核苷酸多態(tài)性(SNP)數(shù)據(jù)，目前對(duì)SNP的研究還處于探索階段，所以不可能根據(jù)其先驗(yàn)知識(shí)來(lái)進(jìn)行分區(qū)，而且大部分的SNP數(shù)據(jù)并不能根據(jù)其前幾個(gè)主成分的散點(diǎn)圖來(lái)進(jìn)行分區(qū)，除非原始數(shù)據(jù)有很好的群體結(jié)構(gòu)。

本文提出的基于主成分聚類(lèi)分析法的非線(xiàn)性檢驗(yàn)度量，有以下優(yōu)點(diǎn):首先，該度量避免了先驗(yàn)知識(shí)匱乏或前幾個(gè)主成分的散點(diǎn)圖不能給出任何分區(qū)結(jié)構(gòu)的不足。主成分聚類(lèi)分析法是對(duì)主成分分析與聚類(lèi)分析方法的綜合利用，利用主成分分析的結(jié)果作為聚類(lèi)分析的樣本矩陣，不僅減少了數(shù)據(jù)的冗余信息，使得聚類(lèi)計(jì)算較為簡(jiǎn)單，而且原理清晰，所得結(jié)論客觀(guān)實(shí)際，可靠性強(qiáng)。其次，在根據(jù)相關(guān)系數(shù)矩陣每一個(gè)元素的置信限估計(jì)其精確邊界時(shí)，本文采用粒子群優(yōu)化算法來(lái)求解最優(yōu)化問(wèn)題〔7〕。粒子群優(yōu)化算法屬于進(jìn)化算法的一種，和遺傳算法相似，它也是從隨機(jī)解出發(fā)，通過(guò)迭代尋找最優(yōu)解，并通過(guò)適應(yīng)度來(lái)評(píng)價(jià)解的優(yōu)良。但粒子群優(yōu)化算法比遺傳算法更容易實(shí)現(xiàn)，并且沒(méi)有許多參數(shù)需要調(diào)整，目前這種算法已廣泛應(yīng)用于函數(shù)優(yōu)化，神經(jīng)網(wǎng)絡(luò)訓(xùn)練，模糊系統(tǒng)控制以及其他遺傳算法的應(yīng)用領(lǐng)域。

兩份模擬數(shù)據(jù)分別驗(yàn)證了線(xiàn)性模型和非線(xiàn)性模型中該方法的有效性。在線(xiàn)性例子中，所有的殘差都落在估計(jì)出的精確邊界之內(nèi)，說(shuō)明該模型可以很好地檢測(cè)出數(shù)據(jù)之間的線(xiàn)性結(jié)構(gòu);而在非線(xiàn)性的例子中，至少有一個(gè)殘差落在精確邊界之外，從而說(shuō)明該非線(xiàn)性檢驗(yàn)度量對(duì)非線(xiàn)性模型的有效性。對(duì)實(shí)際基因組關(guān)聯(lián)分析數(shù)據(jù)PTPN22基因中的9個(gè)SNP的分析結(jié)果表明，該數(shù)據(jù)是非線(xiàn)性數(shù)據(jù)，所以應(yīng)該采用非線(xiàn)性主成分分析構(gòu)建基因組區(qū)域化關(guān)聯(lián)檢驗(yàn)統(tǒng)計(jì)量;如果采用線(xiàn)性的主成分分析構(gòu)建，則不能代表原來(lái)的數(shù)據(jù)結(jié)構(gòu)，從而不能充分地挖掘SNP數(shù)據(jù)中所蘊(yùn)含的大量信息。

通過(guò)本文提出的非線(xiàn)性檢驗(yàn)度量，我們可以更多地了解數(shù)據(jù)的結(jié)構(gòu)，從而選擇最優(yōu)的模型來(lái)挖掘數(shù)據(jù)中所蘊(yùn)含的信息。

1．Pearson K．Principal components analysis．The London，Edinburgh，and Dublin Philosophical Magazine and Journal of Science，1901，6(2):559．

2．Botelho SSC，Bem RA，Almeida IL，et al．C-NLPCA:extracting non-linear principal components of image dataset．Simpoio Brasileiro de Sensoriamento Remoto，2005:3495-3502．

3．Heo G，Gader P，F(xiàn)rigui H．2009 Special Issue:RKF-PCA:robust kernel fuzzy PCA．Neural Networks，2009，22(5-6):642-650．

4．Schlkopf B，Smola A，Muller K．Kernel principal component analysis．Artificial Neural Networks-ICANN＇97，1997:583-588．

5．Kramer MA．Nonlinear principal component analysis using autoassociative neural networks．AIChE Journal，1991，37(2):233-243．

6．Kruger U，Antory D，Hahn J，et al．Introduction of a nonlinearity measure for principal component models．Computers ＆ Chemical Engineering，2005，29(11-12):2355-2362．

7．Kennedy J，Eberhart R．Particle swarm optimization．IEEE International Conference on Neural Networks，1995:1942-1948．

8．Peng Q，Zhao J，Xue F．PCA-based bootstrap confidence interval tests for gene-disease association involving multiple SNPs．BMC Genet，2010，11:6．

9．Plenge RM，Seielstad M，Padyukov L，et al．TRAF1-C5 as a risk locus for rheumatoid arthritis-A genomewide study．New England Journal of Medicine，2007，357(12):1199-1209．

10．李翊，張靜，吳凌華，等．一種基于改進(jìn)主成分分析的SAR圖像識(shí)別方法研究．海軍航空工程學(xué)院學(xué)報(bào)，2009，24(3):307-312．

11．李岳，溫熙森，呂克洪．基于核主成分分析的鐵譜磨粒特征提取方法研究．國(guó)防科技大學(xué)學(xué)報(bào)，2007，29(2):113-116．

12．祝榮欣，權(quán)龍哲，喬金友，等．核主成分分析在農(nóng)業(yè)機(jī)械性能綜合評(píng)價(jià)中的應(yīng)用．中國(guó)造船，2009(9):187-189．