1、偏最小二乘回归方法PLS偏最小二乘回归方法1偏最小二乘回归方法 （PLS） 背景介绍在经济管理、教育学、农业、社会科学、工程技术、医学和生物学中，多元线性回归分 析是一种普遍应用的统计分析与预测技术。多元线性回归中，一般采用最小二乘方法 （Ordinary Least Squares :OLS） 估计回归系数，以使残差平方和达到最小，但当自变量之 间存在多重相关性时， 最小二乘估计方法往往失效。 而这种变量之间多重相关性问题在多元 线性回归分析中危害非常严重，但又普遍存在。为消除这种影响，常采用主成分分析 （principal Components Analysis :PCA） 的方法，但采用
2、主成分分析提取的主成分，虽然能 较好地概括自变量系统中的信息，却带进了许多无用的噪声，从而对因变量缺乏解释能力。最小偏二乘回归方法 （Partial Least Squares Regression ： PLS）就是应这种实际需要 而产生和发展的一种有广泛适用性的多元统计分析方法。 它于 1983年由 S.Wold 和 C.Albano 等人首次提出并成功地应用在化学领域。 近十年来， 偏最小二乘回归方法在理论、 方法和应 用方面都得到了迅速的发展， 己经广泛地应用在许多领域， 如生物信息学、 机器学习和文本 分类等领域。偏最小二乘回归方法主要的研究焦点是多因变量对多自变量的回归建模， 它与普
3、通多元 回归方法在思路上的主要区别是它在回归建模过程中采用了信息综合与筛选技术。 它不再是 直接考虑因变量集合与自变量集合的回归建模， 而是在变量系统中提取若干对系统具有最佳 解释能力的新综合变量 （又称成分 ） ，然后对它们进行回归建模。 偏最小二乘回归可以将建模 类型的预测分析方法与非模型式的数据内涵分析方法有机地结合起来， 可以同时实现回归建 模、数据结构简化 （主成分分析 ）以及两组变量间的相关性分析 （典型性关分析 ） ，即集多元线 性回归分析、 典型相关分析和主成分分析的基本功能为一体。 下面将简单地叙述偏最小二乘 回归的基本原理。2偏最小二乘法的工作目标2.1偏最小二乘法的工作目
4、标在一般的多元线性回归模型中，如果有一组因变量 Y=y 1, ,yq 和一组自变量X=x 1, ,x p ，当数据总体能够满足高斯马尔科夫假设条件时，根据最小二乘法，有Y =X（XTX）-1 XTYY 将是 Y 的一个很好的估计量。从这个公式容易看出，由于（ XTX）必须是可逆矩阵，所以当 X 中的变量存在严重的多重相关性时， 或者在 X 中的样本点数与变量个数相比显然过少时， 这个最小二乘估计都会失效并将引发一系列应用方面的困难。 考虑到这个问题，偏最小二乘回归分析提出了采用成分提取的方法。在主成分分析中， 对于单张数据表 X，为了找到能最好地概括原数据的综合变量， 在 X 中提取了第一主成
5、分 F1, 使得 F1 中所包含的原数据变异信息可达到最大，即Var（F 1） max在典型相关分析中， 为了从整体上研究两个数据表之间的相关关系， 分别在 X 和 Y中提取了 典型成分 F1 和 G1，它们满足r（F 1,G1） maxF1T F 1=1G1 G1=1在能够达到相关度最大的综合变量 F1 和 G1 之间，如果存在明显的相关关系，则可以认为，在两个数据表之间亦存在相关关系。提取成分的做法在数据分析的方法中十分常见， 除主成分、 典型成分以外， 常见到的还 有 Fisher 判别法中的判别成分。实际上，如果 F 是 X 数据表的某种成分，则意味着 F 是 X中变量的某一线性组合
6、F=Xa,而 F 作为一个综合变量，它在 X中所综合提取的信息，将满足 我们特殊的分析需要。2.2偏最小二乘回归分析的建模方法设有 q个因变量 y1, ,y q和 p 个自变量 x 1, ,x p, 为了研究因变量与自变量的统计关系，观测 n个样本点，由此构成了自变量与因变量的数据表 X=【x1, ,x p】n*p和 Y=【 y1, ,y q】n*q。偏最小二乘法回归分别在 X与 Y中提取出 t1和 u1（也就是说， t 1是 x1,x p的线性组合， u1 是 y1, ,y q的线性组合 ） 。在提取这两个成分时，为了回归分析的需要，有下列两个要求： （1） t1和 u1应尽可能大地携带它们
7、各自数据表中的变异信息（ 2） t 1 和 u1 的相关程度能达到最大这两个要求表明， t1和 u1应尽可能好地代表数据表 X和 Y，同时自变量的成分 t 1对因变 量的成分 u1 又有最强的解释能力。在第一个成分 t1和 u1被提取后，偏最小二乘法回归分别实施 X对 t 1的回归以及 Y对 t1的回归。如果方程达到了满意的精度，则算法终止；否则，将利用 X 被 t 1解释后的残余信 息以及 Y 被 t1 解释后的残余信息进行第二轮的成分提取。如此递推，直到能达到一个较为 满意的精度为止。 若最终对 X共提取了 m个成分 t1, ,t m，偏最小二乘法回归将通过实施 YK对 t1,t m的回归
8、，然后再表达成 YK关于原变量 x1,x p的回归方程， k=1, ,q 。3计算方法推导3.1普遍采用的计算推导过程 为了数学推导方便起见，首先将数据做标准化处理。 X 经标准化处理后的数据矩阵记为E0=(E01, ,E 0P) n*p,Y 经过标准化处理后的数据矩阵记为 F0=(F 01, ,F 0q) n*q。第一步，记 t 1是E0的第一个成分， t 1=E0w1, w1是 E0的第一个轴，它是一个单位向量，即|w 1|=1 ；记 u1是 F0的第一个成分， u1=F0c1, c1是 F0的第一个轴，它是一个单位向量， 即|c 1|=1 如果要 t 1,u 1能分别很好德代表 X 与
9、Y 中的数据变异信息，根据主成分分析原理，应该 有Var(t 1) maxVar(u 1) max 另一方面，由于回归建模的需要，又要求 t 1对 u1有最大的解释能力，由典型相关分析的思路， t1与 u1的相关度应达到最大值，即r(t 1,u 1) max 因此综合起来，在偏最小二乘回归中，我们要求 t 1 与 u1 协方差达到最大，即Cov(t 1,u 1)=即求解下列优化问题maxw1T w 1=1 ( 3-1 )Tc1 c 1=1因此，将在 |w 1|=1 和 |c 1|=1 的约束条件下，去求( w1TE0TF0c1)的最大值。此种情况 下我们就可以用拉格朗日算法求其最优解，记s=
10、w 1TE0TF0c1- 1(w1 T w 1-1)- 2(c1T c 1-1 )对 s 分别求关于 w1、 c1、 1、 2的偏导，并令之为零，有E0F0c1-21 w1=0 （3-2 ）F 0TE0w1-2 2 c 1=0 （3-3）-（ w 1T w 1-1）=0 （3-4 ）-（ c 1T c 1-1）=0 （3-5 ）由（ 3-2 ） （ 3-5 ）可以推出TT2 1=2 2= w 1 E0 F0c1=1=2 1=2 2= w1TE0TF0c1,所以 ? 1是优化问题的目标函数值。求得轴 w1和 c1 后，即可得到成分t 1=E0w1u1=F0c1然后，分别求 E0和F0对 t1和
11、u1的回归方程归方程的残差矩阵。F2, h 小于 X 的秩。E1 t 2P2T E2 , F1 t 2r T2同理可推得第 h 成分 th , h 的个数可以用交叉有效性原则进行如此计算下去，如果 X 的秩为 A，则会有E0=t 1P1 + +t APAF0= t 1r 1 +t Ar A +FA由于 t 1, ,t A均可以表示成 E01, ,E0P的线性组合，因此，上式可以还原成 YK= F0K关于XJ=E0J 的回归方程形式YK=bk1X1+ + b kPXP+FAK k=1,.,q3.2一种简洁的计算推导过程3.1中介绍的推导思路是最为常见的，在 3.2 中将介绍一种更为简洁的计算方法
12、，即直 接在 E0, ,E m-1矩阵中提取成分 t 1, ,t m(mp) 。要求 t h能尽可能多地携带 X中的信息，同时， t h 对因变量系统 F0 有最大的解释能力。这时无需在 F0 中提取成分 uh，并且在迭代算法中也无需使用其残差矩阵，而始终直接用 F0 进行计算。这可以使计算过程大为简化，并且对算法结论的解释也更为方便。下面讨论成分 t1,t m(m=A,A=R(X)的一种新原则。在 3.1 中推导偏最小二乘法回归 算法时，第一步的思路是在因变量 F0 抽取一个成分 u1=F0c1，同时在自变量 E0中抽取一个成分 t 1=E0w1, 成分的抽取原则是 max。在这个原则下得知
13、 w1， c1， u1，t 1的计算方法如下：(1)w1是矩阵 E0TF0F0TE0最大特征值的特征向量，成分 t 1=E0w1；( 2) c1是矩阵 F0TE0E0TF0最大特征值的特征向量，成分 u1=F0c1；在求得成分 u1， t 1以后，分别实施 E0在t1上的回归，并生成残差矩阵 E1，以及 F0在 t1 上的回归，得到残差矩阵 F1。再以 E1，F1 取代 E0，F0进行第二轮成分的提取计算，注意到成分 u1,u m是不参加回归计算的，因此是否可以考虑不提取因变量的成分呢？ 为此，用下述原则提取比变量中的成分 t2是与 3.1 中介绍的方法， 结果是完全等价的，即由于 F0K 是
14、标准化变量，所以Cov(F0K, E 0w1)= r(F 0K, E 0w1)因此，该优化原则是求成分 t 1=E0w1，使得 t 1能携带尽可能多的 E0变异，同时， t 1对因变量 F0K(k=1, ,q) 的解释能力会综合达到最大值。由于在目标函数上配上常量( n-1)2 不影响其求解，即为了求 w1 采用拉格朗日算法求解，记可得对 s 求关于 w1和 1的偏导，并令之为零，得由式（ 3-9 ）可知TT2 E 0 F0F0 E0w1-2 1 w1=0 ( 3-9 )w1E0 F0F0 E0w1=1 w12 T T T T = w 1 E 0 F0F0 E0w1= w 1因此 1矩阵 E0
15、TF0F0TE0 的最大特征根， w1则是其相应的特征向量。由此可见，在新的原则下， w1仍然是对应于 E0TF0F0TE0 最大特征值的特征向量，而这个新 的原则完全没有提取到 F0成分 u1提取。也就是说， t 1=E0w1 提取可以不依赖对 u1 的提取，而 这种新的原则又从新的角度说明了 t 1的意义。从这个新的原则出发，对 c1,u 1的计算就可以省略。不过，在偏最小二乘法回归的一些解释技术中， 由于 u1 可以较好地概括 F0中的信息， 因此，它常常也是很有用。4应用举例下面将通过两个具体的案例分析 , 以进一步理解偏最小二乘回归的工作过程和它的特 点。4.1应用举例一应用举例一将
16、采用 Linnerud 给出的关于体能训练的数据进行典型相关分析。在这个数 据系统中被观测样本点，是某健身俱乐部的 20 位中年男子。被观测变量分为两组，第一组 是身体特征指标 X，包括：体重、 腰围、脉搏；第二组变量是训练结果指标 Y，包括：单杠、 弯曲、跳高。原始数据表见表 4-1 。表 4-1 原始数据表在简化算法中，对于 h=1,2,3 时，有计算可得：2 1/(n-1) 2=1.272426 2/(n-1) 2=0.038763 3/(n-1) 2=0.026655而成分 t h的方差， uk的方差以及 th与 u k相关系数的平方 r2(t h, u k)在表 4-2 中列出。2表
17、 4-2 Var(t h), Var(u k) 和 r 2(t h, u k)hVar(t h)Var(u k)r 2(t h, u k)12.02522.05030.306620.43811.90710.046430.23551.15050.0983记第 h 个轴是 wk，第 h 个成分 t k 为t k=Eh-1 wh ( h=1,2,3 ) 其中 t k亦可以表示成原自变量 E0的线性组合，即t k=E0wh*h1* T * 则 wh*= (1-wj pjT)wh。表 4-3 给出 wh*与 wh的取值。j1表 4-3 w h* 与 wh的取值自变量w1w2w3w1*w2*w3*X1-0
18、.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 对 t h进行回归时，有回归系数向量 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成分 t k=Eh-1 wh 的取值见表 4-5 。表 4-5
19、 t k 取值表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
20、.66717-0.393-0.2030.564181.1990.7830.092191.0480.3730.319201.942-1.1290.568通过交叉验证的方法可得，之取一个成分 t1 时，拟合方程的预测性为最佳，不过为了后面作图和解释的方便起见，我们取两个成分 t 1，t 2拟合预测模型。y k=r 1kt 1+ r 2kt 2 k=1,2,3 由于成分 t h 可以写成自变量 xj 的函数，即有t h=wh1*x1+ wh2*x2+ w h3*x3 由此可得两个成分 t 1，t 2所建立的偏最小二乘回归模型为 yk=r 1k( w11 x1+ w 12 x2+ w13 x3)+ r
21、 2k( w21 x1+ w22 x2+ w23 x3) = ( r 1k w11 + r 2k w21 ) x1+( r 1k w12 + r 2k w 22 ) x2+( r 1k w13 + r 2k w23 ) x3 回归系数的计算结果见表 4-6 。表 4-6 回归系数 r kKr1r2r310.34160.33630.477220.41600.29070.455430.14290.0651-0.2125所以，有F01=-0.077E 01-0.499 E 02- 0.132 E 03F02=-0.138E 01-0.524 E 02- 0.085 E 03F01=-0.060E 0
22、1-0.156 E 02- 0.007 E 03 将标准化变量 Fok(k=1,2,3) 和 Eoj (j=1,2,3) 分别还原成原始变量， yk(k=1,2,3) 以及xj (j=1,2,3), 则回归方程为：为了快速直观地观察出各个自变量在解释 Yk时的作用，可以绘制回归系数图， 见图 4-1图 4-1 回归系数的直方图从回归系数图中可以立刻观察到， 腰围变量在解释三个回归方程时起到了极为重要的作 用，然而， 与单杠及弯曲相比， 跳高成绩的回归方程显然不够理想，三个自变量对它的解释能力均很低。因此有必要考虑对自变量做适当的调整。 为了考察这三个回归方程的模型精度，我们以( yik*,y
23、ik )为坐标值，对所有的样本点绘制预测图。 yik*是第 k 个变量，第 i 个样本点( yik )的预测值。在预测直方图上，如果所有样本 点都能在图的对角线附近均匀分布， 则方程的拟合值与原值差异很小， 这个方程的拟合效果 就是满意的。体能训练的预测图如 4-2 所示。4.2 应用举例二这是 Cornell 在 1990年采用的一个化工方面的例子。 此后, 偏最小二乘的提出者 S.Wold等人多次引用 , 成为单因变量偏最小二乘回归的一个经典案例。该例中 , 有个自变量 x1x7,因变量记为 y, 如表 4-7 所示：表 4-7 自变量和应变量对照表x1- 直接蒸馏成分x2 - 重整汽池
24、;x3- 原油热裂化油x4- 原油催化裂化油 ;x5- 聚合物x6 - 烷基化物x7- 天然香精y - 原辛烷值表 4-8 给出了 12 种混合物关于这 8 个变量的观测数据。要求建立 y 对 x1x7, 的回归方 程, 以确定 7 种构成元素 x1x7对 y 的影响。表 4-8 12 种混合物关于 8 个变量的观测数据表这 8 个变量的相关系数矩阵见表 4-9 。从相关系数矩阵中可以看出，在自变量之间存在 严重的多重相关性，例如 r(x 1,x 3)=0.999, r(x 4,x 7)=0.92, r(x 1,x 6)=-0.80 。实际上，这 7 个自变量之间有如下关系： x 1+x2+
25、+x7=1表 4-9 8 个变量的相关系数矩阵由于 q420.0975，所以选择 h=3, 即采用 t1,t 2,t 3三个成分做偏最小二乘回归模型 , 预 测效果最好。从所得到的最终模型看 ,x 6 的回归系数值最大 , 它与 y 正相关。这一点符合我们的基本 认识。 x5的回归系数仍然出现反常符号 , 但它的取值很低 , 几乎可以忽略。 从相关系数表中 可以看出 ,x 5 与 y 的相关度不高 , 并且它与其他自变量之间也没有密切联系。也就是说 , x 5是一个相对独立的变量 , 它不能直接解释 y, 甚至也很难通过其他自变量的传递作用去解 释 y 。因此 , 它在最终模型中的回归系数非常低。与普通最小二乘回归方程相比 , 这个方程的实际含义更加清晰 , 也更易于应用。
