偏最小二乘回归方法(PLS)

偏最小二乘回归方法1 偏最小二乘回归方法(PLS)背景介绍在经济管理、教育学、农业、社会科学、工程技术、医学和生物学中，多元线性回归分析是一种普遍应用的统计分析与预测技术。多元线性回归中，一般采用最小二乘方法(Ordinary Least Squares :OLS)估计回归系数，以使残差平方和达到最小，但当自变量之间存在多重相关性时，最小二乘估计方法往往失效。而这种变量之间多重相关性问题在多元线性回归分析中危害非常严重，但又普遍存在。为消除这种影响，常采用主成分分析(principal Components Analysis :PCA)的方法，但采用主成分分析提取的主成分，虽然能较好地概括自变量系统中的信息，却带进了许多无用的噪声，从而对因变量缺乏解释能力。最小偏二乘回归方法(Partial Least Squares Regression：PLS)就是应这种实际需要而产生和发展的一种有广泛适用性的多元统计分析方法。它于1983年由S.Wold和C.Albano等人首次提出并成功地应用在化学领域。近十年来，偏最小二乘回归方法在理论、方法和应用方面都得到了迅速的发展，己经广泛地应用在许多领域，如生物信息学、机器学习和文本分类等领域。偏最小二乘回归方法主要的研究焦点是多因变量对多自变量的回归建模，它与普通多元回归方法在思路上的主要区别是它在回归建模过程中采用了信息综合与筛选技术。它不再是直接考虑因变量集合与自变量集合的回归建模，而是在变量系统中提取若干对系统具有最佳解释能力的新综合变量(又称成分)，然后对它们进行回归建模。偏最小二乘回归可以将建模类型的预测分析方法与非模型式的数据内涵分析方法有机地结合起来，可以同时实现回归建模、数据结构简化(主成分分析)以及两组变量间的相关性分析(典型性关分析)，即集多元线性回归分析、典型相关分析和主成分分析的基本功能为一体。下面将简单地叙述偏最小二乘回归的基本原理。2 偏最小二乘法的工作目标2.1 偏最小二乘法的工作目标在一般的多元线性回归模型中，如果有一组因变量Y=y1,yq和一组自变量X=x1,xp，当数据总体能够满足高斯马尔科夫假设条件时，根据最小二乘法，有=X（XTX）-1XTY将是Y的一个很好的估计量。从这个公式容易看出，由于（XTX）必须是可逆矩阵，所以当X中的变量存在严重的多重相关性时，或者在X中的样本点数与变量个数相比显然过少时，这个最小二乘估计都会失效并将引发一系列应用方面的困难。考虑到这个问题，偏最小二乘回归分析提出了采用成分提取的方法。在主成分分析中，对于单张数据表X，为了找到能最好地概括原数据的综合变量，在X中提取了第一主成分F1,使得F1中所包含的原数据变异信息可达到最大，即Var(F1)max在典型相关分析中，为了从整体上研究两个数据表之间的相关关系，分别在X和Y中提取了典型成分F1和G1，它们满足r(F1,G1) maxF1T F1=1G1T G1=1在能够达到相关度最大的综合变量F1和G1之间，如果存在明显的相关关系，则可以认为，在两个数据表之间亦存在相关关系。提取成分的做法在数据分析的方法中十分常见，除主成分、典型成分以外，常见到的还有Fisher判别法中的判别成分。实际上，如果F是X数据表的某种成分，则意味着F是X中变量的某一线性组合F=Xa,而F作为一个综合变量，它在X中所综合提取的信息，将满足我们特殊的分析需要。2.2 偏最小二乘回归分析的建模方法设有q个因变量y1,yq和p个自变量x1,xp,为了研究因变量与自变量的统计关系，观测n个样本点，由此构成了自变量与因变量的数据表X=【x1,xp】n*p和Y=【y1,yq】n*q。偏最小二乘法回归分别在X与Y中提取出t1和u1(也就是说，t1是x1,xp的线性组合，u1是y1,yq的线性组合)。在提取这两个成分时，为了回归分析的需要，有下列两个要求：（1） t1和u1应尽可能大地携带它们各自数据表中的变异信息（2） t1和u1的相关程度能达到最大这两个要求表明，t1和u1应尽可能好地代表数据表X和Y，同时自变量的成分t1对因变量的成分u1又有最强的解释能力。在第一个成分t1和u1被提取后，偏最小二乘法回归分别实施X对t1的回归以及Y对t1的回归。如果方程达到了满意的精度，则算法终止；否则，将利用X被t1解释后的残余信息以及Y被t1解释后的残余信息进行第二轮的成分提取。如此递推，直到能达到一个较为满意的精度为止。若最终对X共提取了m个成分t1,tm，偏最小二乘法回归将通过实施YK对t1,tm的回归，然后再表达成YK关于原变量x1,xp的回归方程，k=1,q。3 计算方法推导3.1 普遍采用的计算推导过程为了数学推导方便起见，首先将数据做标准化处理。X经标准化处理后的数据矩阵记为E0=(E01,E0P)n*p,Y经过标准化处理后的数据矩阵记为F0=(F01,F0q)n*q。第一步，记t1是E0的第一个成分，t1=E0w1, w1是E0的第一个轴，它是一个单位向量，即|w1|=1；记u1是F0的第一个成分，u1=F0c1, c1是F0的第一个轴，它是一个单位向量，即|c1|=1。如果要t1,u1能分别很好德代表X与Y中的数据变异信息，根据主成分分析原理，应该有Var(t1)maxVar(u1)max另一方面，由于回归建模的需要，又要求t1对u1有最大的解释能力，由典型相关分析的思路，t1与u1的相关度应达到最大值，即r(t1,u1)max因此综合起来，在偏最小二乘回归中，我们要求t1与u1协方差达到最大，即Cov(t1,u1)=Var(t1)Var(u1)r(t1,u1)max即求解下列优化问题maxw1T w1=1 （3-1）c1T c1=1因此，将在|w1|=1和|c1|=1的约束条件下，去求（w1TE0TF0c1）的最大值。此种情况下我们就可以用拉格朗日算法求其最优解，记s= w1TE0TF0c1-1(w1T w1-1)- 2（c1T c1-1）对s分别求关于w1、c1、1、2的偏导，并令之为零，有 E0TF0c1-21 w1=0 （3-2） F0TE0w1-22 c1=0 （3-3） -( w1T w1-1)=0 （3-4） -( c1T c1-1)=0 （3-5）由（3-2）（3-5）可以推出21=22= w1TE0TF0c1=记1=21=22= w1TE0TF0c1,所以1是优化问题的目标函数值。把式（3-2）和式（3-3）写成 E0TF0c1= 1 w1 （3-6） F0TE0w1= 1 c1 （3-7） 将式（3-7）代入式（3-6），有 E0TF0F0TE0w1= 12 w1 (3-8)由式（3-8）可知，w1是矩阵E0TF0F0TE0特征向量，对应的特征值为12，1是目标函数值，要求取得其最大值，所以w1是对应于矩阵E0TF0F0TE0最大特征值12的单位特征向量。求得轴w1和c1后，即可得到成分t1=E0w1u1=F0c1然后，分别求E0和F0对t1和u1的回归方程其中，向量；E1，F1*，F1为回归方程的残差矩阵。 第2成分t2的提取，以E1取代E0 , F1取代F0 , 用上面的方法求第2个轴W2和第2个成分t2 ,有 同样，E1 , F1分别对t2做回归, 得到同理可推得第h 成分th , h 的个数可以用交叉有效性原则进行, h 小于X 的秩。如此计算下去，如果X的秩为A，则会有E0=t1P1T+tAPATF0= t1r1T+tArAT+FA由于t1,tA均可以表示成E01,E0P的线性组合，因此，上式可以还原成YK= F0K关于XJ=E0J的回归方程形式YK=bk1X1+ bkPXP+FAK k=1,.,q3.2一种简洁的计算推导过程3.1中介绍的推导思路是最为常见的，在3.2中将介绍一种更为简洁的计算方法，即直接在E0,Em-1矩阵中提取成分t1,tm(mp)。要求th能尽可能多地携带X中的信息，同时，th对因变量系统F0有最大的解释能力。这时无需在F0中提取成分uh，并且在迭代算法中也无需使用其残差矩阵，而始终直接用F0进行计算。这可以使计算过程大为简化，并且对算法结论的解释也更为方便。下面讨论成分t1,tm(m=A,A=R(X)的一种新原则。在3.1中推导偏最小二乘法回归算法时，第一步的思路是在因变量F0抽取一个成分u1=F0c1，同时在自变量E0中抽取一个成分t1=E0w1,成分的抽取原则是max。在这个原则下得知w1，c1，u1，t1的计算方法如下：（1）w1是矩阵E0TF0F0TE0最大特征值的特征向量，成分t1=E0w1；（2）c1是矩阵F0TE0E0TF0最大特征值的特征向量，成分u1=F0c1；在求得成分u1，t1以后，分别实施E0在t1上的回归，并生成残差矩阵E1，以及F0在t1上的回归，得到残差矩阵F1。再以E1，F1取代E0，F0进行第二轮成分的提取计算，注意到成分u1,um是不参加回归计算的，因此是否可以考虑不提取因变量的成分呢？为此，用下述原则提取比变量中的成分t2是与3.1中介绍的方法，结果是完全等价的，即由于F0K是标准化变量，所以Cov(F0K, E0w1)=Var(EOw1) r(F0K, E0w1)因此，该优化原则是求成分t1=E0w1，使得t1能携带尽可能多的E0变异，同时，t1对因变量F0K(k=1,q)的解释能力会综合达到最大值。由于在目标函数上配上常量（n-1）2不影响其求解，即（n-1）2Cov2(F0K, E0w1)= 2= w1T E0TF0KF0KTE0w1= w1T E0T(F0KF0KT)E0w1= w1T E0TF0F0TE0w1为了求w1采用拉格朗日算法求解，记s= 2-1(w1T w1-1)= w1T E0TF0F0TE0w1-1(w1T w1-1)对s求关于w1和1的偏导，并令之为零，得 2 E0TF0F0TE0w1-21 w1=0 （3-9） -( w1T w1-1)=0 （3-10）由式（3-9）可知E0TF0F0TE0w1=1 w1可见，最优解w1应是矩阵E0TF0F0TE0的一个特征向量，将它代入目标函数，并且由式（3-10）可得2= w1T E0TF0F0TE0w1= w1T （1 w1）=1因此1矩阵E0TF0F0TE0的最大特征根，w1则是其相应的特征向量。由此可见，在新的原则下，w1仍然是对应于E0TF0F0TE0最大特征值的特征向量，而这个新的原则完全没有提取到F0成分u1提取。也就是说，t1=E0w1提取可以不依赖对u1的提取，而这种新的原则又从新的角度说明了t1的意义。从这个新的原则出发，对c1,u1的计算就可以省略。不过，在偏最小二乘法回归的一些解释技术中，由于u1可以较好地概括F0中的信息，因此，它常常也是很有用。4 应用举例下面将通过两个具体的案例分析, 以进一步理解偏最小二乘回归的工作过程和它的特点。4.1 应用举例一应用举例一将采用Linnerud给出的关于体能训练的数据进行典型相关分析。在这个数据系统中被观测样本点，是某健身俱乐部的20位中年男子。被观测变量分为两组，第一组是身体特征指标X，包括：体重、腰围、脉搏；第二组变量是训练结果指标Y，包括：单杠、弯曲、跳高。原始数据表见表4-1。表4-1 原始数据表 在简化算法中，对于h=1,2,3时，有h=2=（n-1）2Cov2(F0K, th)计算可得： 1/(n-1)2=1.2724262/(n-1)2=0.0387633/(n-1)2=0.026655而成分th的方差，uk的方差以及th与 uk相关系数的平方r2(th, uk)在表4-2中列出。表4-2 Var(th), Var(uk)和r2(th, uk)hVar(th)Var(uk)r2(th, uk)12.02522.05030.306620.43811.90710.046430.23551.15050.0983记第h个轴是wk，第h个成分tk为tk=Eh-1wh （h=1,2,3）其中tk亦可以表示成原自变量E0的线性组合，即tk=E0wh*则wh*=（1-wjpjT）wh。表4-3给出wh*与wh的取值。表4-3 wh*与wh的取值自变量w1w2w3w1*w2*w3*X1-0.589890.46879-0.65747-0.589890.36793-0.93459X2-0.77134-0.568010.28706-0.77134-0.699890.80231X30.23888-0.67647-0.696660.23888-0.63562-0.22282在利用Eh-1对th进行回归时，有回归系数向量pk，h=1,2,3,见表4-4。表4-4 回归系数pkp1P2P3-0.6659-0.0197-0.6574-0.6760-0.35460.28700.3589-1.1942-0.6966成分tk=Eh-1wh的取值见表4-5。表4-5 tk取值表NOt1t2t31-0.6430.591-0.1312-0.7700.1670.1343-0.907-0.5210.04840.688-0.6800.3465-0.4871.133-0.1826-0.229-0.0720.0257-1.404-0.077-0.57280.744-0.211-0.03291.715-0.655-1.557101.1630.1670.333110.3650.7010.201120.7430.6980.002131.187-0.7570.36614-4.390-0.7600.25515-0.8230.974-0.08316-0.749-0.521-0.66717-0.393-0.2030.564181.1990.7830.092191.0480.3730.319201.942-1.1290.568通过交叉验证的方法可得，之取一个成分t1时，拟合方程的预测性为最佳，不过为了后面作图和解释的方便起见，我们取两个成分t1，t2拟合预测模型。yk=r1kt1+ r2kt2 k=1,2,3由于成分th可以写成自变量xj的函数，即有th=wh1*x1+ wh2*x2+ wh3*x3由此可得两个成分t1，t2所建立的偏最小二乘回归模型为yk=r1k（w11*x1+ w12*x2+ w13*x3）+ r2k（w21*x1+ w22*x2+ w23*x3） =（r1k w11*+ r2k w21*）x1+（r1k w12*+ r2k w22*）x2+（r1k w13*+ r2k w23*）x3回归系数的计算结果见表4-6。表4-6 回归系数rkKr1r2r310.34160.33630.477220.41600.29070.455430.14290.0651-0.2125所以，有F01=-0.077E01-0.499 E02-0.132 E03F02=-0.138E01-0.524 E02-0.085 E03F01=-0.060E01-0.156 E02-0.007 E03将标准化变量Fok(k=1,2,3)和Eoj(j=1,2,3)分别还原成原始变量，yk(k=1,2,3)以及xj(j=1,2,3),则回归方程为：Y1*=47.02-0.0166x1-0.824x2-0.097x3Y2*=612.57-0.351x1-10.52x2-0.741x3Y3*=183.98-0.125x1-2.497x2-0.052x3为了快速直观地观察出各个自变量在解释Yk时的作用，可以绘制回归系数图，见图4-1图4-1 回归系数的直方图从回归系数图中可以立刻观察到，腰围变量在解释三个回归方程时起到了极为重要的作用，然而，与单杠及弯曲相比，跳高成绩的回归方程显然不够理想，三个自变量对它的解释能力均很低。因此有必要考虑对自变量做适当的调整。为了考察这三个回归方程的模型精度，我们以（yik*,yik）为坐标值，对所有的样本点绘制预测图。yik*是第k个变量，第i个样本点（yik）的预测值。在预测直方图上，如果所有样本点都能在图的对角线附近均匀分布，则方程的拟合值与原值差异很小，这个方程的拟合效果就是满意的。体能训练的预测图如4-2所示。4.2 应用举例二这是Cornell在1990年采用的一个化工方面的例子。此后,偏最小二乘的提出者S.Wold等人多次引用, 成为单因变量偏最小二乘回归的一个经典案例。该例中,有个自变量x1x7, 因变量记为y, 如表4-7所示：表4-7 自变量和应变量对照表-直接蒸馏成分-重整汽池;-原油热裂化油-原油催化裂化油;-聚合物-烷基化物-天然香精-原辛烷值表4-8给出了12种混合物关于这8个变量的观测数据。要求建立y对x1x7,的回归方程, 以确定7种构成元素x1x7对y的影响。表4-8 12种混合物关于8个变量的观测数据表这8个变量的相关系数矩阵见表4-9。从相关系数矩阵中可以看出，在自变量之间存在严重的多重相关性，例如r(x1,x3)=0.999, r(x4,x7)=0.92, r(x1,x6)=-0.80。实际上，这7个自变量之间有如下关系： x1+x2+x7=1表4-9 8个变量的相关系数矩阵由于q420.0975，所以选择h=3, 即采用t1,t2,t3三个成分做偏最小二乘回归模型, 预测效果最好。从所得到的最终模型看,x6的回归系数值最大, 它与y正相关。这一点符合我们的基本认识。x5的回归系数仍然出现反常符号, 但它的取值很低, 几乎可以忽略。从相关系数表中可以看出,x5与y的相关度不高,并且它与其他自变量之间也没有密切联系。也就是说, x5是一个相对独立的变量, 它不能直接解释y, 甚至也很难通过其他自变量的传递作用去解释y。因此, 它在最终模型中的回归系数非常低。与普通最小二乘回归方程相比, 这个方程的实际含义更加清晰, 也更易于应用。