R语言实现偏最小二乘回归

1、第33卷第2期2013年6月数学理论与应月MATHEMATICAL THEOKY AM) APPLICATIONSVol. 33 N<>. 2Jun. 2013偏最小二乘建模在R软件中的实现及实证分析齐琛方秋莲(中南大学数学与统计学院，长沙,410075 )摘 要 通过介绍僞最小二乘(FIS)的建模和显著性检验原理，解决了小样本多变量且变量间存在多重共线 性的回归问题，建立了多变量对多变量的回归模型，并使用K软件(版本为K i386 2. 15. 1 )实现了 PLS建模； 最后基于葡萄和葡萄酒理化指标数据进行了实证分析.关键词 偏最小二乘 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引言人们总能通过普通彊小二乘冋归进行数据的分析和预测；然而半选取的解释变戢过多而 样本很少时无法使用普通最小二乘法进行建模.Herman Wohl在20世纪70年代的经济学研 究中提出偏最小二乘法（PLS）,它能够在小样本的情况下实行多变貳对多变就的冋归建模. 1998年王惠文"对偏彊小二乘冋归方法及其应用进行了详尽的解说在2005年与吴载斌、孟 洁一起，对彊小二乘法再进行了进一步的扩展，并捉出冋归系数

5、的显著性检验方法；高惠 璇°用具体例子对赧小二乘冋归、主成分冋归和偏赧小二乘冋归进行比较分析.并使用SAS 软件实现了 PU建模.R软件是一款具右强大统计分析功能的开源软件，利用K软件进行偏湖南省H然科学基金资助项目 繚号U2JJ5002 ） 收稿日期:2013年4月28日 偏眾小二乘建模在K软件屮的实现及实证分析111啟小二乘冋归建模可以得到理想的模世，并能够对冋归系数进行显著性检验，帮助人们发现 变量的主要影响因素，进行下一步数据探索分析.本文介绍偏最小二乘的冋归原理以及如何 使用K软件中的pls包4实现P於建模，并尝试用PLS对匍萄和匍萄酒理化指标数据进行实 证分析.2偏最小二

6、乘回归及显著性检验原理2.1偏最小二乘回归原理设有g个因变就怙，力，人和p个自变量 尿,%2,灼.观测了 "个样本点，由此构 成了门变量与因变量的数据表X,®和人“偏最小二乘冋归分别在x与y中提取出成分“与 心（即：卩、呦分别是心宀，*p、X2,认的线性组合）在提取这两个成分时，需满足以b 两个条件：（1）人与划应尽可能大地携带它们备自数据表屮的变异信息；Q）/i与山的相关程度达到最大.这两个要求表明，人与妁应尽可能好地代表数据x与y,同时i'i变址的成分“对因变量 的成分心乂有较强的解释能力.在第一个成分林与山被提取后，偏最小二乘回归分别实施X对如的回归及Y对&q

7、uot;的冋 归.如果冋归方程C经达到满意的梢确度则算法终止;否则将利用x与y分别被人解释后的 残余信息进行第二轮的成分提取如此往复玄到能达到一个较满总的精度为止若最终对x 提取F m个成分»一偏最小二乘回归将实行以对勺，儿的回归，然后再表达成九关于 原变量心，久2,秸的回归方程（k = 1 ,2,g ）.2.2偏最小二乘回归步骤Step 1 :先将数据进行标准化处理得到标准化后的X与丫矩阵，记第i对成分为Z,与色， 并且人=Xw, u, = 丫“.于是对第一对成分的提取，即求解以下优化问题：只需求出矩阵M = X YY X的特征值与特征向址，其虽大特征值儿对应的特征向址即为 所求的

8、口，B标函数值等于庆Step 2 :分别做人2,和衍，2,为对"的冋归其中，冋归系数向量从二几，仙），7i '=如，恂）；坨与儿是冋归方程的残差阵' 和仞可山简单最小二乘法的原则求得.Step 3 :用伦与儿代替X与丫进行前两个步骤求得第二对成分，依次循环.设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根据交叉验证结果选择模型的成分个数若选取成分的个数

9、过多，会很容易;I；现过度拟合的问题，因此我们需耍一个右效的原则來 确定成分的个数.采用类似抽样测试的丁作方式把所有样本点分成两部分:第一部分用来重 新拟合一个偏最小二乘模型，第二部分的样本点作为测试数据;带入拟合模型中求得预测值误 差平方和PRESS二X匕-必尺再以这种方式重复g次,直到所有的样本都被预测了一次, 赧后把每个样本的预测误差平方和加总称为PKESS.PRESS = £ PRESS（2.4 ）常见的交叉验证方法有“留一验证'，“K折交叉验证'，"Holdout验证”等方法，选取一种 方法分别求出第1 r个成分对应的PRESS值取PRESS赧小的

10、或者PKESS几乎不再变化的 成分个数作为最终模型选取的成分个数m.2.4回归系数的显著性检验在R软件的pls包（package pls冲提供了函数jack, test,用来检验回归系数的显著性.由 于偏赧小二乘法不同于一般最小二乘法它的冋归系数方差无法得到准确的无偏估计.Miller 提出 Quenouille - Tukey jackknife方法来估计冋归系数的方差：先抽出g个样本子集，然后用 只去除一个子集的样本做一次偏彊小二乘的冋归系数估计记去除第i个样本子集对应的冋 归系数为则jackknife方差估计为var © ） =£ J -p）20.5 ）S <

11、* I其中0是心的均值，皿常见的巫抽样法是jackknife法，即每次选一个样本点，于是共右“ 个样本子集（即g = n）.在估计出方差后，类似于普通彊小二乘法求出仅对应的t统计就再 进行均值是否为零的假设检脸.山于偏赧小二乘冋归妾数确定的复杂公式，我们至今无法确定 准确的/分布在R软件中默认为服从白由度为加的/分布（JU为建模使用成分的个数）.3在R软件中的实现首先;需要在加载R的程庁包pls； pls包是由Bjoni - Helge Mevik. Kon Welnens和Kristian Hovde IJland创建，专门用来做偏最小二乘冋归的程序包.代码如下：> libran (&

12、quot; pls" Jib. loc = " C： /Program Files/R/K -2. 15. 1 /library")再导人H变戢和因变就的样木数据，并且使用scale ()函数将数据进行标准化消除戢纲的 影响记标准化之后的门变量为X,因变量为Y,进行PLS冋归的代码如下：> plsl < - plsr (Y X,validation = " LOO" Jackknife 二TRUE )#进行偏最小二乘冋归,模型存为对象plsl> summary (plsl, what = HallH ) #显示回归结果(包括P

13、KESS与变异解释度)苴中validation = * L00-表示使用昭一交叉验证计算PRESS, jackknife = TRUE表示使用 jackknife法估计冋归系数方差(为后面的显著性检验做准备).在没给定成分个数的情况下，会 默认使用所冇的主成分进行冋归，因此我们需要在选择的成分个数尽可能小的前提下,选择使 PRESS最小或几乎不变的成分个数.假设选定了成分个数为新进行冋归，并对冋归系数 假设检验代码如下：> pls2 < - plsr (Y X, neon屮二 m, validation = " LOO" , jackknife = TRUE )

14、# “ncomp二m”表示模型成分个数为m> jack, test (pls2 )另外还可以使用coef ()函数得到冋归系数scores ()得到得分矩阵.loadings ()得到载荷矩 predict ()得到对应样本的预测值，以及plot ()函数将结果以图的形式展现.4基于葡萄和葡萄酒理化指标的PLS实证分析匍萄酒是山制萄精细酿造而成，因此二者的理化指标之间必然存在一定的联系本文采用 中国2012年数学建模大赛A题中提供的数据对红匍萄酒的理化指标与酿酒红匍萄的理化指 标进行最小二乘法建模分析(以下的葡萄酒与酿酒匍萄均指红匍萄酒与酿酒红匍萄).4.1建模过程Step 1 :导入数

15、据，并进行数据的标准化：> 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年数学建模大赛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个主成分对应的各项结果此处Ll达到最大个数25 )选取部分结果如表1所示：表1初步模型拟合结果(部分)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即为不同主成分个数对应的PRESS,adjcv为调整后的PRESS, TRAINING： % variance explained *'一栏为主成分对各变就的累积贡献率.山结果可知，主成分个数为3个时.模型在经过留一交叉验证法后求得的PRESS总和最 小,随着成分个数的增加PRESS值也没有太大改变，并且3个成分对各个因变就的累积贡献 率也基本达到了 85%，因此定下冋归的成分个数m = 3.Step 3:根据成分数皿=3,建立最终模型:> pls2 < - plsr (Y X, ncomp = 3 , validation = &qu

19、ot; LOO" Jackknife = TRUE )>coef (pls2 ) #得到回归系数得到冋归系数后,便能耳出各因变址对所右解释变呈的冋归方程，下面耳出对各解释 变杲的冋归方程(由于变昴太多，因此中间有省略)：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 ()函数可以画出PLS模型在不同主成分数下对应的RMSEP

20、 (山留一交 叉验证法算得的均方预测误差根)，对初始模型的结果进行画图如图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不同成分数对应的均方误差图眦处仅裁取前9个因变量)图1屮纵坐标“1M1ESP”表示均方预测误差根，横坐标为不同模型的成分个数；由图1可 知大部分因变竝在成分数为3时对应的均方误差根赧小证明选择3个成分参与建模是正确 的.使用predplot (pls2 )函数画出最终模型的预测效果图，如图2所示.图2

21、屮纵坐标为各因变量的预测值，横坐标为各因变呈的实际测址值;散点集中分布在主对角线上則说明预测效果很好.图2中15张预测图射应15个因变就)的散点大致都分布在 对角线上，说明最终模型的拟合效果较好然而,Y8对应的散点图几乎是一条垂直的线预测 很糟糕;査阅原始数据得他第26个样本点的Y8原始值是1.6239,远大于其他只有0.02左右 的数据，因此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回归参数的显著性检验使用jack, test ()函数进行检验挑出与通过检验的冋归系数对应的|'|变址，整理结果如表2所示（由于篇幅限制，只列出前2个因变量对应的结果）：表2回归系数显著性统计表葡萄酒理朋酒筍萄系数筍萄ifl理廉酒葡萄系数化指标¥理化指标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+注：显甘性符号衷示*'极英卫誉*'非常显暮佯很掘富校显蔷.通过显著性检验可以知道各因变就（葡萄酒理化指标）受哪些H变竝（酿酒匍萄理化指 标）的影响较大，及其受影响程度.山表2可血，对因变量Y1花色昔而言，对其右显著影响的 H变就有X19蛋白质、X21花色背、X27 DPPH白由墓、X28总酚、X29单宁、X37杨梅黄酮、X40 杲鼠7索，并均足対具有止向的影响.服酒制匍的花色仲址刈匍萄酒的花色口伉有

25、止向影响 是显.然成立的，此处也通过显著性检验，因此可初步判断此模型与实际相符，对于其他影响 显著的变誠便是给了我们个探索点，可以通过别的方法深入探讨同理于苴他14个因变量, 宦著性检验可以比我们初步了解到各因变最的受影响因索.5总结偏赧小二乘冋归能够解决许多以往用普通多元线性【叫归不能解决的问题，在解释变讯:个 数大于样本个数的情况下也能建立出很有效的模型本文主要介绍的是使用R软件的|血包 进行偏处小二乘建模，成功地对超萄酒和酿酒匍萄的理化指标进行偏最小二乘的建模然而 由于円$公式的复杂性，对于冋归系数的方差估计至今没有特别完善的方法，因此我们需要 辩证地看jack, test显著性检验的结果，以它作为一个研究的参考再进一步进行深入分析.偏眾小二乘建模在K软件屮的实现及实证分析#参考文献1 王惠文.偏最小二乘冋归方法及其应用LMJ.北京：国防丁业出版社，1999.2 王惠文.吴载斌.孟洁.備最小二乘回归的线性与非线