穩(wěn)健主成分回歸及其醫(yī)學(xué)應(yīng)用*

劉偉新 郭東星

主成分回歸(principal component regression PCR)可用于處理數(shù)據(jù)中的自變量多重共線性問題，而這個(gè)問題在醫(yī)學(xué)資料的數(shù)據(jù)分析中經(jīng)常會出現(xiàn)。但是，當(dāng)數(shù)據(jù)中存在異常點(diǎn)時(shí)會影響主成分回歸模型的建立，在進(jìn)行主成分回歸分析時(shí)除了可以進(jìn)行異常點(diǎn)的診斷〔1－2〕外，本文將再介紹一種穩(wěn)健主成分回歸方法來解決這個(gè)問題。

原理和方法

1．穩(wěn)健主成分分析(ROBPCA)〔3〕

ROBPCA(robust principal component analysis)是一種穩(wěn)健主成分分析的新方法。其結(jié)合了投影尋蹤(projection pursuit)的思想〔4，5〕和穩(wěn)健的協(xié)方差矩陣〔6，7〕估計(jì)的思想。

步驟一:降低數(shù)據(jù)空間維數(shù)

設(shè)原始數(shù)據(jù)為Xn，p。n表示觀測樣本數(shù)，p表示自變量的原始數(shù)目。通過中心化數(shù)據(jù)矩陣的奇異值分解，產(chǎn)生

令Zn，r0=UD，變成新的數(shù)據(jù)矩陣，這個(gè)奇異值分解就是數(shù)據(jù)在V中的r0列所占據(jù)的子空間中的仿射轉(zhuǎn)置，并且不會失去任何信息。為了方便起見，新的數(shù)據(jù)仍然用xi來表示。

步驟二:找到h(h＜n)個(gè)‘最小異常值’的數(shù)據(jù)點(diǎn)

事先不知道異常點(diǎn)的數(shù)目，就讓h=max{［αn］，［(n+kmax+1)/2］}，kmax代表將要被計(jì)算的主成分的最大數(shù)目，默認(rèn)為10。參數(shù)α在0.5和1之間選擇。默認(rèn)時(shí)α=0.75。

(1)對每一個(gè)數(shù)據(jù)點(diǎn)xi，計(jì)算它的異常值(outlyingness)，從中找到h個(gè)具有最小異常值的數(shù)據(jù)點(diǎn)，

(2)用^μ1和S0來代表H0中的h個(gè)觀測點(diǎn)的均數(shù)和協(xié)方差矩陣。協(xié)方差矩陣的特征值按降序排列，而且特征向量被相應(yīng)地標(biāo)記。則

L=diag(~l1，…，~lr)，是特征值的對角矩陣，r≤ri

協(xié)方差矩陣S0決定在未來的分析中，將要保留的主成分的個(gè)數(shù)k0(k0≤r)。在這過程可以用多種方法達(dá)到這個(gè)目的，例如，可以觀察特征值的單調(diào)遞減的一個(gè)斜線圖，或者能利用累積貢獻(xiàn)率的選擇標(biāo)準(zhǔn)。

(3)將數(shù)據(jù)點(diǎn)投影到S0的前k0個(gè)特征向量所在的子空間上。即令:

這里Pr1，k0包括了(3)中的P0的前k0列。步驟三:利用MCD估計(jì)值〔7〕，穩(wěn)健地估計(jì)X*

n，k0中的數(shù)據(jù)點(diǎn)的方差-協(xié)方差矩陣。

需要找到h個(gè)數(shù)據(jù)點(diǎn)，使它們的協(xié)方差矩陣有著最小的行列式。利用由步驟二得到的異常值測量(2)

(4)將P2的列轉(zhuǎn)換回原始空間，產(chǎn)生最后的穩(wěn)健特征向量矩陣Pp，k。最后穩(wěn)健中心^μ通過將^μ5轉(zhuǎn)換回原始空間來獲得，而最后的p維秩為k的穩(wěn)健分散矩陣 S 由 S=Pp×kLk×kP'p×k給出。公式(8)中的得分在Rp中可以寫為Tn×k=(Xn×p－1n^μ')Pp×k。

穩(wěn)健主成分分析的部分至此完成。

2．穩(wěn)健回歸法則-LTS(least trimmed squares〔8－10〕)

在穩(wěn)健主成分回歸中，用重新加權(quán)的LTS法，眾所周知，最小二乘法回歸是將殘差平方和最小化。對于LTS方法，殘差平方和被殘差平方的修剪的和所代替。殘差平方由低到高排列，然后從最低的殘差到排

^y－i，k是主成分?jǐn)?shù)目為k時(shí)，觀測點(diǎn)i暫時(shí)作為驗(yàn)證樣本時(shí)從已建模型中獲得的預(yù)測值。這時(shí)具有最小的RMSECV值所對應(yīng)的k被認(rèn)為是最優(yōu)的數(shù)目。

在穩(wěn)健主成分回歸中，因?yàn)椴幌朐赗MSECV值中包括異常點(diǎn)的預(yù)測誤差，所以我們用穩(wěn)健的RMSECV值

w－i由上面穩(wěn)健回歸時(shí)的公式所得到。

由此，具有最小的R-RMSECV值所對應(yīng)的k被認(rèn)為是最優(yōu)的數(shù)目。

(2)穩(wěn)健的R2值

n

實(shí)例與分析

采用17所醫(yī)院的人力利用情況及有關(guān)醫(yī)院任務(wù)的資料〔1〕，其中，X1為平均每天住院人數(shù);X2為每月X線照光人數(shù);X3為每月占病床天數(shù);X4為服務(wù)范圍內(nèi)人口數(shù)(千人);X5為每名病人平均住院天數(shù);Y為每月使用人力(小時(shí))。以下結(jié)果采用SAS 8.0和MATLAB 7.1軟件編程實(shí)現(xiàn)。

1.常規(guī)描述分析和共線性診斷見表1，表2。

表1 17所醫(yī)院的人力利用及醫(yī)院任務(wù)情況的簡相關(guān)矩陣

表1中可見，X1與 X2，X3，X4的相關(guān)系數(shù)、X2與X3，X4、X3與X4的相關(guān)系數(shù)均大于90%，提示自變量間可能有多重共線性存在。

表2 17所醫(yī)院的人力利用及醫(yī)院任務(wù)情況的共線性診斷

在表2中可以看出，條件指數(shù)11.416，33.881，390.423均大于 10，對應(yīng)的方差比 0.728，0.944，0.999，0.998大于0.50，因此確定自變量之間存在多重共線性。進(jìn)行穩(wěn)健主成分回歸分析。

2.進(jìn)行主成分個(gè)數(shù)選擇

圖1，2中表示，選擇第一，第二主成分后穩(wěn)健的R-RMSECV值最小，且穩(wěn)健的R2值最高，所以選擇第一和第二兩個(gè)主成分。

3．穩(wěn)健主成分分析ROBPCA結(jié)果:穩(wěn)健均數(shù)估計(jì)值^μ=

［81.2328 10595.3629 2435.8449 65.4358 5.4438］，穩(wěn)健特征值為［84131696.2849 517807.0308］，穩(wěn)健

數(shù)據(jù)中心化，可以計(jì)算得到穩(wěn)健的主成分。

圖 1 穩(wěn)健的 RMSECV 值(1912，757．4，851，1023，639．2)

圖 2 穩(wěn)健的值(0．9871，0．9960，0．9951，0．9963，0．9964)

4．兩個(gè)主成分與因變量的穩(wěn)健回歸結(jié)果

參數(shù)估計(jì)值^φ1=0.2256，^φ2=0.6436，截距=2683.029，R2=0.99435，所以方程為 ^Y=2683.029+0.2256T1+0.6436T2

反代回原始自變量得:

^Y=34.9719+0.0219X1+0.0946X2+0.6751X3+0.0018X4+0.0005X5。

5．穩(wěn)健主成分回歸還將產(chǎn)生可以診斷異常點(diǎn)的診斷圖:

圖3 穩(wěn)健主成分分析的主成分診斷圖

由圖3 可見，點(diǎn)10，14，15，16，17 為異常點(diǎn)。其中點(diǎn)14，15，17為無影響PCA杠桿點(diǎn)，點(diǎn)10為正交異常點(diǎn)，點(diǎn)16為有影響PCA杠桿點(diǎn)。

圖4 穩(wěn)健主成分回歸的回歸方面的異常點(diǎn)診斷圖

由圖4可見，點(diǎn)15，17為無影響異常點(diǎn)，點(diǎn)9，10為垂直異常點(diǎn)，點(diǎn)14，16為有影響異常點(diǎn)。

圖5 經(jīng)典的主成分分析的主成分診斷圖

由圖5可見，只有點(diǎn)10，17被診斷為異常點(diǎn)。

圖6 經(jīng)典主成分回歸的回歸方面的異常點(diǎn)診斷圖

由圖6可見，只有點(diǎn)9和點(diǎn)16，17被診斷為異常點(diǎn)，其他三個(gè)異常點(diǎn)并沒有被診斷出來。

討 論

本文所介紹的穩(wěn)健主成分回歸方法是由兩部分穩(wěn)健方法組成:穩(wěn)健主成分分析方法ROBPCA和穩(wěn)健回歸方法即重新加權(quán)的LTS法。這兩種方法均為目前最新的非常穩(wěn)健的方法。這兩部分方法的失效點(diǎn)均達(dá)到50%。這種穩(wěn)健主成分回歸方法計(jì)算速度快，并且對于低維和高維的數(shù)據(jù)都能夠處理。而且既可以用在有異常點(diǎn)的數(shù)據(jù)，也可以用于沒有異常點(diǎn)的數(shù)據(jù)。在實(shí)例分析中表明，當(dāng)數(shù)據(jù)中包含異常點(diǎn)時(shí)，與經(jīng)典的主成分回歸相比較，此方法得到了穩(wěn)健的估計(jì)值，并且診斷圖對于確定異常點(diǎn)也非常有用。

通常在診斷出異常點(diǎn)以后，不能簡單地將異常點(diǎn)刪除，因?yàn)檫@樣做可能將異常點(diǎn)攜帶的一些有用的信息丟失，所以應(yīng)該對不同情況的異常點(diǎn)給予不同處理。如果證實(shí)是數(shù)據(jù)錄入錯(cuò)誤，可以刪除。而多數(shù)情況下，剔除只是一種識別數(shù)據(jù)是否異常的方法，不是診斷分析的最終目的。對于處理的方法，除了有本文提到的穩(wěn)健估計(jì)方法外，還需要以后進(jìn)一步的研究和探討。

1．郭東星，劉偉新．主成分回歸中異常點(diǎn)的穩(wěn)健診斷．中國衛(wèi)生統(tǒng)計(jì)，2008，25(1):31-34．

2．劉偉新，郭東星．主成分回歸中異常點(diǎn)的二步診斷法及其醫(yī)學(xué)應(yīng)用．現(xiàn)代預(yù)防醫(yī)學(xué)，2007，34(13):2423-2425．

3．Mia H，Peter JR．ROBPCA:a new approach to robust principal component analysis．Technometrics，2005，47:64-79．

4．Jolliffe IT．Principal component analysis．New York:Springer，1986．

5．Li G，Chen Z．Projection-Pursuit approach to robust dispersion matrices and principal components:primary theory and Monte Carlo．Journal of A-merican statistical association，1982，80:759-766．

6．Rousseeuw PJ，Van DK．A fast algorithm for the minimum covariance determinant estimator．Technometrics，1999，41:212-223．

7．Croux C，Haesbroeck G．Influence function and efficiency of the minimum covariance determinant scatter matrix estimator．Journal of Multivariate A-nalysis，1999，71:161-190．

8．Pell RJ．Multiple outlier detection for multivariate calibration using robust statistical techniques．Chemometrics and Intelligent Laboratory Systems，2000，52:87-104．

9．Rousseeuw PJ，Leroy A．Robust regression and outlier detection．New York:John wiley，1987．

10．Walczak B．Outlier detection in multivariate calibration．Chemometrics and Intelligent Laboratory Systems，1998，28:259-272．