R語言實現(xiàn)偏最小二乘回歸

1、第33卷第2期2013年6月數(shù)學理論與應月MATHEMATICAL THEOKY AM) APPLICATIONSVol. 33 N<>. 2Jun. 2013偏最小二乘建模在R軟件中的實現(xiàn)及實證分析齊琛方秋蓮(中南大學數(shù)學與統(tǒng)計學院，長沙,410075 )摘 要 通過介紹僞最小二乘(FIS)的建模和顯著性檢驗原理，解決了小樣本多變量且變量間存在多重共線 性的回歸問題，建立了多變量對多變量的回歸模型，并使用K軟件(版本為K i386 2. 15. 1 )實現(xiàn)了 PLS建模； 最后基于葡萄和葡萄酒理化指標數(shù)據(jù)進行了實證分析.關(guān)鍵詞 偏最小二乘 R語言 jackknife方差 顯著性檢

2、臉Partial Least Squares Modelling with R Softwareand Empirical AnalysisQi Chen Fang Qiulian(Sc hool of Mathrniatics and Statistics, Crntral South Univnity, Changsha 410075, China )Abstract Tliis paper intmduers the Partial Ix-ast Squairs (PIS ) method and its significance test principle for modelling

3、 regression problems in which sample size is small and there is multicollinraiity among observable variablesv and furtheinioit*, illustrates how to set up PLS models with the K software An rxamplr to moilrl the irlation of the pliysi- cochrmical indexes lx-twren gia|)cs and wine is given to deinonst

4、rate the iiKxlvlIing process.Key words FIS K Language Jackknife Vaiian(*e Significance Test1引言人們總能通過普通彊小二乘冋歸進行數(shù)據(jù)的分析和預測；然而半選取的解釋變戢過多而 樣本很少時無法使用普通最小二乘法進行建模.Herman Wohl在20世紀70年代的經(jīng)濟學研 究中提出偏最小二乘法（PLS）,它能夠在小樣本的情況下實行多變貳對多變就的冋歸建模. 1998年王惠文"對偏彊小二乘冋歸方法及其應用進行了詳盡的解說在2005年與吳載斌、孟 潔一起，對彊小二乘法再進行了進一步的擴展，并捉出冋歸系數(shù)

5、的顯著性檢驗方法；高惠 璇°用具體例子對赧小二乘冋歸、主成分冋歸和偏赧小二乘冋歸進行比較分析.并使用SAS 軟件實現(xiàn)了 PU建模.R軟件是一款具右強大統(tǒng)計分析功能的開源軟件，利用K軟件進行偏湖南省H然科學基金資助項目 繚號U2JJ5002 ） 收稿日期:2013年4月28日 偏眾小二乘建模在K軟件屮的實現(xiàn)及實證分析111啟小二乘冋歸建?？梢缘玫嚼硐氲哪Ｊ?，并能夠?qū)諝w系數(shù)進行顯著性檢驗，幫助人們發(fā)現(xiàn) 變量的主要影響因素，進行下一步數(shù)據(jù)探索分析.本文介紹偏最小二乘的冋歸原理以及如何 使用K軟件中的pls包4實現(xiàn)P於建模，并嘗試用PLS對匍萄和匍萄酒理化指標數(shù)據(jù)進行實 證分析.2偏最小二

6、乘回歸及顯著性檢驗原理2.1偏最小二乘回歸原理設(shè)有g(shù)個因變就怙，力，人和p個自變量 尿,%2,灼.觀測了 "個樣本點，由此構(gòu) 成了門變量與因變量的數(shù)據(jù)表X,®和人“偏最小二乘冋歸分別在x與y中提取出成分“與 心（即：卩、呦分別是心宀，*p、X2,認的線性組合）在提取這兩個成分時，需滿足以b 兩個條件：（1）人與劃應盡可能大地攜帶它們備自數(shù)據(jù)表屮的變異信息；Q）/i與山的相關(guān)程度達到最大.這兩個要求表明，人與妁應盡可能好地代表數(shù)據(jù)x與y,同時i'i變址的成分“對因變量 的成分心乂有較強的解釋能力.在第一個成分林與山被提取后，偏最小二乘回歸分別實施X對如的回歸及Y對&q

7、uot;的冋 歸.如果冋歸方程C經(jīng)達到滿意的梢確度則算法終止;否則將利用x與y分別被人解釋后的 殘余信息進行第二輪的成分提取如此往復玄到能達到一個較滿總的精度為止若最終對x 提取F m個成分»一偏最小二乘回歸將實行以對勺，兒的回歸，然后再表達成九關(guān)于 原變量心，久2,秸的回歸方程（k = 1 ,2,g ）.2.2偏最小二乘回歸步驟Step 1 :先將數(shù)據(jù)進行標準化處理得到標準化后的X與丫矩陣，記第i對成分為Z,與色， 并且人=Xw, u, = 丫“.于是對第一對成分的提取，即求解以下優(yōu)化問題：只需求出矩陣M = X YY X的特征值與特征向址，其雖大特征值兒對應的特征向址即為 所求的

8、口，B標函數(shù)值等于慶Step 2 :分別做人2,和衍，2,為對"的冋歸其中，冋歸系數(shù)向量從二幾，仙），7i '=如，恂）；坨與兒是冋歸方程的殘差陣' 和仞可山簡單最小二乘法的原則求得.Step 3 :用倫與兒代替X與丫進行前兩個步驟求得第二對成分，依次循環(huán).設(shè)X的秩為 r （r Wp）,則存在個主成分，使得(2.2 )Q.3)X 二 /|/>|' + E,丫 =佃+ *M/ + 人再把人=X叭帶入Q. 2 ）即可得到兒6 兒分別對 “2，已 的冋歸方程y, = Bj0 + 時p G = l,2,g）2.3根據(jù)交叉驗證結(jié)果選擇模型的成分個數(shù)若選取成分的個數(shù)

9、過多，會很容易;I；現(xiàn)過度擬合的問題，因此我們需耍一個右效的原則來 確定成分的個數(shù).采用類似抽樣測試的丁作方式把所有樣本點分成兩部分:第一部分用來重 新擬合一個偏最小二乘模型，第二部分的樣本點作為測試數(shù)據(jù);帶入擬合模型中求得預測值誤 差平方和PRESS二X匕-必尺再以這種方式重復g次,直到所有的樣本都被預測了一次, 赧后把每個樣本的預測誤差平方和加總稱為PKESS.PRESS = £ PRESS（2.4 ）常見的交叉驗證方法有“留一驗證'，“K折交叉驗證'，"Holdout驗證”等方法，選取一種 方法分別求出第1 r個成分對應的PRESS值取PRESS赧小的

10、或者PKESS幾乎不再變化的 成分個數(shù)作為最終模型選取的成分個數(shù)m.2.4回歸系數(shù)的顯著性檢驗在R軟件的pls包（package pls沖提供了函數(shù)jack, test,用來檢驗回歸系數(shù)的顯著性.由 于偏赧小二乘法不同于一般最小二乘法它的冋歸系數(shù)方差無法得到準確的無偏估計.Miller 提出 Quenouille - Tukey jackknife方法來估計冋歸系數(shù)的方差：先抽出g個樣本子集，然后用 只去除一個子集的樣本做一次偏彊小二乘的冋歸系數(shù)估計記去除第i個樣本子集對應的冋 歸系數(shù)為則jackknife方差估計為var © ） =£ J -p）20.5 ）S <

11、* I其中0是心的均值，皿常見的巫抽樣法是jackknife法，即每次選一個樣本點，于是共右“ 個樣本子集（即g = n）.在估計出方差后，類似于普通彊小二乘法求出僅對應的t統(tǒng)計就再 進行均值是否為零的假設(shè)檢臉.山于偏赧小二乘冋歸妾數(shù)確定的復雜公式，我們至今無法確定 準確的/分布在R軟件中默認為服從白由度為加的/分布（JU為建模使用成分的個數(shù)）.3在R軟件中的實現(xiàn)首先;需要在加載R的程庁包pls； pls包是由Bjoni - Helge Mevik. Kon Welnens和Kristian Hovde IJland創(chuàng)建，專門用來做偏最小二乘冋歸的程序包.代碼如下：> libran (&

12、quot; pls" Jib. loc = " C： /Program Files/R/K -2. 15. 1 /library")再導人H變戢和因變就的樣木數(shù)據(jù)，并且使用scale ()函數(shù)將數(shù)據(jù)進行標準化消除戢綱的 影響記標準化之后的門變量為X,因變量為Y,進行PLS冋歸的代碼如下：> plsl < - plsr (Y X,validation = " LOO" Jackknife 二TRUE )#進行偏最小二乘冋歸,模型存為對象plsl> summary (plsl, what = HallH ) #顯示回歸結(jié)果(包括P

13、KESS與變異解釋度)苴中validation = * L00-表示使用昭一交叉驗證計算PRESS, jackknife = TRUE表示使用 jackknife法估計冋歸系數(shù)方差(為后面的顯著性檢驗做準備).在沒給定成分個數(shù)的情況下，會 默認使用所冇的主成分進行冋歸，因此我們需要在選擇的成分個數(shù)盡可能小的前提下,選擇使 PRESS最小或幾乎不變的成分個數(shù).假設(shè)選定了成分個數(shù)為新進行冋歸，并對冋歸系數(shù) 假設(shè)檢驗代碼如下：> pls2 < - plsr (Y X, neon屮二 m, validation = " LOO" , jackknife = TRUE )

14、# “ncomp二m”表示模型成分個數(shù)為m> jack, test (pls2 )另外還可以使用coef ()函數(shù)得到冋歸系數(shù)scores ()得到得分矩陣.loadings ()得到載荷矩 predict ()得到對應樣本的預測值，以及plot ()函數(shù)將結(jié)果以圖的形式展現(xiàn).4基于葡萄和葡萄酒理化指標的PLS實證分析匍萄酒是山制萄精細釀造而成，因此二者的理化指標之間必然存在一定的聯(lián)系本文采用 中國2012年數(shù)學建模大賽A題中提供的數(shù)據(jù)對紅匍萄酒的理化指標與釀酒紅匍萄的理化指 標進行最小二乘法建模分析(以下的葡萄酒與釀酒匍萄均指紅匍萄酒與釀酒紅匍萄).4.1建模過程Step 1 :導入數(shù)

15、據(jù)，并進行數(shù)據(jù)的標準化：> Gl < - read, csv (" K: W0RK 論文 R grape, csv'* )> W1 < - read, csv (*'K： WWORKW論文WRWwine. csvM )> X < - scale (G1 )> Y< -scale (W1 )得到的IT變戢X是27 x59的矩陣，兒他9依次代表釀酒匍萄理化指標如下(在此只列 出部分名稱，具體參見P012年數(shù)學建模大賽A題附件2”)：氨基酸總就、天門冬氨酸、蘇氨酸、絲氨酸、谷氨酸果皮顏色H、果皮顏色C.得到的因變量Y是27x1

16、5的矩陣,X、 ri5依次代表葡萄酒理化指標如下:花色背、單寧、總酚、酒總黃酮色澤b*、色澤H、色澤C.Step 2:進行初步偏域小二乘冋歸：> plsl < - plsr (Y X, ncomp = 10. validation = " LOO" , jackknife = TRUE)> suimnaiy (plsl , what = *'all")#注汕中默認赧多只能顯示25個主成分對應的各項結(jié)果此處Ll達到最大個數(shù)25 )選取部分結(jié)果如表1所示：表1初步模型擬合結(jié)果(部分)VALIDATION：Cross - validated u

17、sing 27 leavKMSEPe 一 one 一 out segmentsRoponx-:Y1Qntrrcrpt )1 comps2 comps3 comps4 compw.CV1.0190.75190.7CO20.71060.6757adjCV1.0190.74900.75400.70580.6759.THAINING： % variance explained1 coinp>2 c!<>inp>3 coinps4 compsX16.1729924.84535.87246. 10Y165.3819874.14781.72882.46Y283.1518783.729

18、84.34384.3&苴中CV即為不同主成分個數(shù)對應的PRESS,adjcv為調(diào)整后的PRESS, TRAINING： % variance explained *'一欄為主成分對各變就的累積貢獻率.山結(jié)果可知，主成分個數(shù)為3個時.模型在經(jīng)過留一交叉驗證法后求得的PRESS總和最 小,隨著成分個數(shù)的增加PRESS值也沒有太大改變，并且3個成分對各個因變就的累積貢獻 率也基本達到了 85%，因此定下冋歸的成分個數(shù)m = 3.Step 3:根據(jù)成分數(shù)皿=3,建立最終模型:> pls2 < - plsr (Y X, ncomp = 3 , validation = &qu

19、ot; LOO" Jackknife = TRUE )>coef (pls2 ) #得到回歸系數(shù)得到冋歸系數(shù)后,便能耳出各因變址對所右解釋變呈的冋歸方程，下面耳出對各解釋 變杲的冋歸方程(由于變昴太多，因此中間有省略)：Y二一 0. 0066人 + 0. 046禺 + 0. 0087X, - 0. 04IX+ 0. 0057九-0. 012X6+ -0. 036*55 - 0. 018X56 - 0. 000495, + 0. 012X58 - 0. 017X594.2模型擬合效果分析使用valiilationplot ()函數(shù)可以畫出PLS模型在不同主成分數(shù)下對應的RMSEP

20、 (山留一交 叉驗證法算得的均方預測誤差根)，對初始模型的結(jié)果進行畫圖如圖1所示：02468 10Y5Y302468 10Y40 24 6 8 10Y70 2 4 6 6 10Y8Y60 24 6 8 100 2 4 6 8 10 number of components02 4 6 S 10圖1不同成分數(shù)對應的均方誤差圖眥處僅裁取前9個因變量)圖1屮縱坐標“1M1ESP”表示均方預測誤差根，橫坐標為不同模型的成分個數(shù)；由圖1可 知大部分因變竝在成分數(shù)為3時對應的均方誤差根赧小證明選擇3個成分參與建模是正確 的.使用predplot (pls2 )函數(shù)畫出最終模型的預測效果圖，如圖2所示.圖2

21、屮縱坐標為各因變量的預測值，橫坐標為各因變呈的實際測址值;散點集中分布在主對角線上則說明預測效果很好.圖2中15張預測圖射應15個因變就)的散點大致都分布在 對角線上，說明最終模型的擬合效果較好然而,Y8對應的散點圖幾乎是一條垂直的線預測 很糟糕;査閱原始數(shù)據(jù)得他第26個樣本點的Y8原始值是1.6239,遠大于其他只有0.02左右 的數(shù)據(jù)，因此26號樣本為模型的強影響點，它使得對Y8的擬合效果很差；可以占慮將苴剔除 重新進行預測對于其他的變戢,可以直觀的看出效果是不錯的.Y1 J eompt. /UMfttonY2.3ctKnpf vaiaioonY3 3 wnw wMJiOwY8. 3 co

22、mos.vaiiaaOonA 0 bQb %0OLs .- o«o Oa oo?：°。oOo° o O。瞽r9 9 <? oV 0特9 Q V o 0, 1 Y4.3 comps, adat“r$. 3 comps vahUBon1MlclplgV7. 3 8mt efcdBi"ri.Sconce EMatknW 3comps.va«aanonmeasuredV11.3axn(n, valiOatnv> 0 O OAo e o° o Q <?Oo； 爲O °公°aq %0 Y10. 3omps Q

23、OstionV12,3 comps vahdatonY113 compt *3onY15.3 comps vUidatonr>9 -2： e 圖2預測效果圖4.3回歸參數(shù)的顯著性檢驗使用jack, test ()函數(shù)進行檢驗挑出與通過檢驗的冋歸系數(shù)對應的|'|變址，整理結(jié)果如表2所示（由于篇幅限制，只列出前2個因變量對應的結(jié)果）：表2回歸系數(shù)顯著性統(tǒng)計表葡萄酒理朋酒筍萄系數(shù)筍萄ifl理廉酒葡萄系數(shù)化指標¥理化指標X符號程度化指標Y理化指標X符石程度X19+X28+* * *X21+X29+* *X27+* *X30+Y1X28+*X34+* *X29+* *X36+X3

24、7+*X38+*X40+*Y2X40+* *XI+*X41+X6+X45+XII+X49+* *Y2X18X51X21+* * *X52+X26+*X53+X27+* * *X55+注：顯甘性符號衷示*'極英衛(wèi)譽*'非常顯暮佯很掘富校顯薔.通過顯著性檢驗可以知道各因變就（葡萄酒理化指標）受哪些H變竝（釀酒匍萄理化指 標）的影響較大，及其受影響程度.山表2可血，對因變量Y1花色昔而言，對其右顯著影響的 H變就有X19蛋白質(zhì)、X21花色背、X27 DPPH白由墓、X28總酚、X29單寧、X37楊梅黃酮、X40 杲鼠7索，并均足対具有止向的影響.服酒制匍的花色仲址刈匍萄酒的花色口伉有

25、止向影響 是顯.然成立的，此處也通過顯著性檢驗，因此可初步判斷此模型與實際相符，對于其他影響 顯著的變誠便是給了我們個探索點，可以通過別的方法深入探討同理于苴他14個因變量, 宦著性檢驗可以比我們初步了解到各因變最的受影響因索.5總結(jié)偏赧小二乘冋歸能夠解決許多以往用普通多元線性【叫歸不能解決的問題，在解釋變訊:個 數(shù)大于樣本個數(shù)的情況下也能建立出很有效的模型本文主要介紹的是使用R軟件的|血包 進行偏處小二乘建模，成功地對超萄酒和釀酒匍萄的理化指標進行偏最小二乘的建模然而 由于円$公式的復雜性，對于冋歸系數(shù)的方差估計至今沒有特別完善的方法，因此我們需要 辯證地看jack, test顯著性檢驗的結(jié)果，以它作為一個研究的參考再進一步進行深入分析.偏眾小二乘建模在K軟件屮的實現(xiàn)及實證分析#參考文獻1 王惠文.偏最小二乘冋歸方法及其應用LMJ.北京：國防丁業(yè)出版社，1999.2 王惠文.吳載斌.孟潔.備最小二乘回歸的線性與非線