第三十章偏最小二乘回归 在实际问题中，经常遇到需要研究两组多重相关变量间的相互依赖关系，并研究用 一组变量（常称为自变量或预测变量）去预测另一组变量（常称为因变量或响应变量） 除了最小二乘准则下的经典多元线性回归分析(MLR),提取自变量组主成分的主成 分回归分析(PCR)等方法外，还有近年发展起来的偏最小二乘(PLS)回归方法。 偏最小二乘回归提供一种多对多线性回归建模的方法，特别当两组变量的个数很 多，日都存在多重相关性，而观测数据的量（样本量）又较少时，用偏最小一乘回归 建立的模型具有传统的经典回归分析等方法所没有的优点 偏最小二乘回归分析在建模过程中集中了主成分分析，典型相关分析和线性回归分 析方法的特点，因此在分析结果中，除了可以提供一个更为合理的回归模型外，还可以 同时完成一些类似于主成分分析和典型相关分析的研究内容，提供更丰富、深入的一些 信息。 本章介绍偏最小二乘回归分析的建模方法：通过例子从预测角度对所建立的回归模 型进行比较。 §1偏最小二乘回归 考虑p个变量乃片，y,与m个自变量x,2,x的建模问题。偏最小二乘 回归的基本作法是首先在自变量集中提出第一成分4(是x,x的线性组合，且 尽可能多地提取原自变量集中的变异信息)：同时在因变量集中也提取第一成分山， 并要求4与山相关程度达到最大。然后建立因变量，y。与1的回归，如果回归方 程已达到满意的精度，则算法中止。否测继续第二对成分的提取，直到能达到满意的精 度为止。若最终对自变量集提取r个成分4,2，.，4，偏最小二乘回归将通过建立 片，.，。与4，山2，.，4，的回归式，然后再表示为，.，y,与原自变量的回归方程式， 即偏最小二乘回归方程式。 为了方便起见，不妨假定p个因变量y,.,y,与m个自变量x,.,x均为标准 化变量。因变量组和自变量组的n次标准化观测数据阵分别记为 .p X1.Xm Fo= ,E。=:: yn.ypx.m 偏最小二乘回归分析建模的具体步骤如下： -531
-531- 第三十章 偏最小二乘回归 在实际问题中，经常遇到需要研究两组多重相关变量间的相互依赖关系，并研究用 一组变量（常称为自变量或预测变量）去预测另一组变量（常称为因变量或响应变量）， 除了最小二乘准则下的经典多元线性回归分析（MLR），提取自变量组主成分的主成 分回归分析（PCR）等方法外，还有近年发展起来的偏最小二乘（PLS）回归方法。 偏最小二乘回归提供一种多对多线性回归建模的方法，特别当两组变量的个数很 多，且都存在多重相关性，而观测数据的数量（样本量）又较少时，用偏最小二乘回归 建立的模型具有传统的经典回归分析等方法所没有的优点。 偏最小二乘回归分析在建模过程中集中了主成分分析，典型相关分析和线性回归分 析方法的特点，因此在分析结果中，除了可以提供一个更为合理的回归模型外，还可以 同时完成一些类似于主成分分析和典型相关分析的研究内容，提供更丰富、深入的一些 信息。 本章介绍偏最小二乘回归分析的建模方法；通过例子从预测角度对所建立的回归模 型进行比较。 §1 偏最小二乘回归 考虑 p 个变量 p y , y , , y 1 2 " 与 m 个自变量 m x , x , , x 1 2 " 的建模问题。偏最小二乘 回归的基本作法是首先在自变量集中提出第一成分 1t （ 1t 是 m x , , x 1 " 的线性组合，且 尽可能多地提取原自变量集中的变异信息）；同时在因变量集中也提取第一成分 1 u ， 并要求 1t 与 1 u 相关程度达到最大。然后建立因变量 p y , , y 1 " 与 1t 的回归，如果回归方 程已达到满意的精度，则算法中止。否则继续第二对成分的提取，直到能达到满意的精 度为止。若最终对自变量集提取 r 个成分 r t ,t , ,t 1 2 " ，偏最小二乘回归将通过建立 p y , , y 1 " 与 r t ,t , ,t 1 2 " 的回归式，然后再表示为 p y , , y 1 " 与原自变量的回归方程式， 即偏最小二乘回归方程式。 为了方便起见，不妨假定 p 个因变量 p y , , y 1 " 与 m 个自变量 m x , , x 1 " 均为标准 化变量。因变量组和自变量组的n 次标准化观测数据阵分别记为 ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = n np p y y y y F " # # " 1 11 1 0 ， ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = n nm m x x x x E " # # " 1 11 1 0 偏最小二乘回归分析建模的具体步骤如下：
(1)分别提取两变量组的第一对成分，并使之相关性达最大。 假设从两组变量分别提出第一对成分为马，和4,4是自变量集X=(x,.,x)了的 线性组合：=+.+xm=X,出是因变量集了=(，.，y,)了的线性组 合：4=M出++y。=Y。为了回归分析的需要，要求： ①4和4各自尽可能多地提取所在变量组的变异信息： ②1，和山，的相关程度达到最大。 由两组变量集的标准化观测数据阵E,和F。,可以计算第一对成分的得分向量，记 为i和立 a.ypp 第一对成分1和山的协方差Cov(化，4)可用第一对成分的得分向量和元，的内积 来计算。故而以上两个要求可化为数学上的条件极值问题 iin >=<Eow,Yov>=wEFox=max w=hiP=1 wy=h=1 利用Lagrange乘数法，问题化为求单位向量和，使日=Fov→最大.问 题的求解只须通过计算m×m矩阵M=EFF。E。的特征值和特征向量，且M的最大特 征值为，相应的单位特征向量就是所求的解州，面%可由斯计算得到=日？E,(。 (2)建立，.，y。对的回归及x,x对的回归。 假定回归模型为 52
-532- （1）分别提取两变量组的第一对成分，并使之相关性达最大。 假设从两组变量分别提出第一对成分为 1t 和 1 u ，1t 是自变量集 T m X (x , , x ) = 1 " 的 线性组合：t w x w x w XT 1 = 11 1 +"+ 1m m = 1 ， 1 u 是因变量集 T p Y ( y , , y ) = 1 " 的线性组 合：u v y v y v Y T 1 = 11 1 +"+ 1p p = 1 。为了回归分析的需要，要求： ① 1t 和 1 u 各自尽可能多地提取所在变量组的变异信息； ② 1t 和 1 u 的相关程度达到最大。 由两组变量集的标准化观测数据阵 E0 和 F0 ，可以计算第一对成分的得分向量，记 为 1t ˆ 和 1 uˆ ： ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = = 1 11 1 11 1 11 1 1 0 1 ˆ n nm m n m t t w w x x x x t E w # # " # # " ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = = 1 11 1 11 1 11 1 1 0 1 ˆ n np p n p u u v v y y y y u F v # # " # # " 第一对成分 1t 和 1 u 的协方差Cov( , ) 1 1 t u 可用第一对成分的得分向量 1t ˆ 和 1 uˆ 的内积 来计算。故而以上两个要求可化为数学上的条件极值问题： ⎪⎩ ⎪ ⎨ ⎧ = = = = < >=< >= ⇒ 1, 1 ˆ , ˆ , max 2 1 1 1 2 1 1 1 1 0 1 0 1 1 0 0 1 w w w v v v t u E w Y v w E F x T T T T 利用Lagrange乘数法，问题化为求单位向量 w1和 1 v ，使θ1 = w1 T E0 T F0v1 ⇒ 最大。问 题的求解只须通过计算 m× m 矩阵 M E0 F0F0 E0 T T = 的特征值和特征向量，且 M 的最大特 征值为 2 θ1 ，相应的单位特征向量就是所求的解 w1，而 1 v 可由 w1计算得到 0 0 1 1 1 1 v F E w T θ = 。 （2）建立 p y , , y 1 " 对 1t 的回归及 m x , , x 1 " 对 1t 的回归。 假定回归模型为
[E。=ia+E Fo=inB'+F 其中a,=(C,.,a1m)了，月=(B,.,Bp)》了分别是多对一的回归模型中的参数向量 E,和F是残差阵。回归系数向量a,B的最小二乘估计为 a=E/作 B=Fi 称a,B为模型效应负荷量 (3)用残差阵E和F代替E。和F。重复以上步骤。 记E。=,户。=B,则残差阵E,=E。-E。,F=F-户。如果残差阵F 中元素的绝对值近似为0，则认为用第一个成分建立的回归式精度已满足需要了，可以 停止抽取成分。否则用残差阵E和F代替E和F重复以上步骤即得： W=(21,.,m)了；=(2，.，2p)了分别为第二对成分的权数。而 2=Ew2,立2=Fy2为第二对成分的得分向量。 么爵么骨物y隆防 Eli, Eo=iaf +iga+Ez F。-R+iB+F (4)设n×m数据阵E的秩为r≤min(n-l,m),则存在r个成分1,2， 使得 E。=ia+.+i,a+E F=i+.+i,B+F 把=w出+.+wmxn(k=l2,.,r),代入Y=1B+.+1,B,即得p个 因变量的偏最小二乘回归方程式 -533
-533- ⎪⎩ ⎪ ⎨ ⎧ = + = + 0 1 1 1 0 1 1 1 ˆ ˆ F u F E t E T T β α 其中 T m ( , , ) α1 = α11 " α1 ， T p ( , , ) β1 = β11 " β1 分别是多对一的回归模型中的参数向量， E1和 F1是残差阵。回归系数向量 1 1 α , β 的最小二乘估计为 ⎪ ⎩ ⎪ ⎨ ⎧ = = 2 1 0 1 1 2 1 0 1 1 ˆ ˆ ˆ ˆ F t t E t t T T β α ， 称 1 1 α , β 为模型效应负荷量。 （3）用残差阵 E1和 F1代替 E0 和 F0 重复以上步骤。 记 T E t 0 1 1 ˆ = ˆα ， T F t 0 1 1 ˆ = ˆ β ，则残差阵 1 0 0 E = E − Eˆ ， 1 0 0 F = F − Fˆ 。如果残差阵 F1 中元素的绝对值近似为 0，则认为用第一个成分建立的回归式精度已满足需要了，可以 停止抽取成分。否则用残差阵 E1和 F1代替 E0 和 F0 重复以上步骤即得： T w w w m ( , , ) 2 = 21 " 2 ； T p v (v , , v ) 2 = 21 " 2 分别为第二对成分的权数。而 2 1 2 t ˆ = E w ， 2 1 2 uˆ = F v 为第二对成分的得分向量。 2 2 1 2 2 ˆ ˆ t E t T α = ， 2 2 1 2 2 ˆ ˆ t F t T β = 分别为 X ,Y 的第二对成分的负荷量。这时有 ⎪⎩ ⎪ ⎨ ⎧ = + + = + + 0 1 1 2 2 2 0 1 1 2 2 2 ˆ ˆ ˆ ˆ F t t F E t t E T T T T β β α α （4）设 n × m 数据阵 E0 的秩为 r ≤ min(n −1,m) ，则存在 r 个成分 r t ,t , ,t 1 2 " ， 使得 ⎪⎩ ⎪ ⎨ ⎧ = + + + = + + + r T r r T r T r r T F t t F E t t E β β α α ˆ ˆ ˆ ˆ 0 1 1 0 1 1 " " 把 k k km m t = w x +"+ w x 1 1 （ k = 1,2,",r ），代入 r r Y = t1β1 +"+ t β ，即得 p 个 因变量的偏最小二乘回归方程式
,=a+.+amxm,（j=l,2,.,m) (5)交叉有效性检验。 一般情况下，偏最小二乘法并不需要选用存在的个成分4，2，4，来建立回归 式，而像主成分分析一样，只选用前1个成分(1≤r),即可得到预测能力较好的回归 模型。对于建模所需提取的主成分个数1，可以通过交叉有效性检验来确定。 每次舍去第i个观测(i=1,2.,n),用余下的n-1个观测值按偏最小二乘回归 方法建模，并考虑抽取h个成分后拟合的回归式，然后把舍去的第i个观测点代入所批 合的回归方程式，得到y,U=1,2.,P)在第i个观测点上的预测值戈，()。对 i=1,2,.,n重复以上的验证，即得抽取h个成分时第j个因变量y,（U=1,2,.,p)的 预测误差平方和为 PRESS(h)=(Y((J=1.2.P) Y=(y,.,y)'的预测误差平方和为 PRESS()=左PRESS,(h: 另外，再采用所有的样本点，拟合含h个成分的回归方程。这时，记第i个样本点 的预测值为少，()，则可以定义y,的误差平方和为 SS,()=∑0y,-,(h》 定义Y的误差平方和为 ss()ss,(h) 当PRESS(h)达到最小值时，对应的h即为所求的成分个数。通常，总有 PRESS(h)大于SS(h),而SS(h)则小于SS(h-I)。因此，在提取成分时，总希望比 值PRESS(h)/SS(h-)越小越好：一般可设定限制值为0.05，即当 -534-
-534- j j jm m y = a x +"+ a x 1 1 ，（ j = 1,2,",m） （5）交叉有效性检验。 一般情况下，偏最小二乘法并不需要选用存在的 r 个成分 r t ,t , ,t 1 2 " 来建立回归 式，而像主成分分析一样，只选用前l 个成分（l ≤ r ），即可得到预测能力较好的回归 模型。对于建模所需提取的主成分个数l ，可以通过交叉有效性检验来确定。 每次舍去第i 个观测（i = 1,2,", n ），用余下的 n −1个观测值按偏最小二乘回归 方法建模，并考虑抽取h 个成分后拟合的回归式，然后把舍去的第i 个观测点代入所拟 合的回归方程式，得到 y ( j 1,2, , p) j = " 在第i 个观测点上的预测值 ˆ ( ) y(i) j h 。对 i = 1,2,", n 重复以上的验证，即得抽取 h 个成分时第 j 个因变量 y ( j 1,2, , p) j = " 的 预测误差平方和为 ∑= = − n i j h yij y i j h 1 2 ( ) PRESS ( ) ( ˆ ( )) （ j = 1,2,", p ） T p Y ( y , , y ) = 1 " 的预测误差平方和为 ∑= = p i h j h 1 PRESS( ) PRESS ( ) 。 另外，再采用所有的样本点，拟合含h 个成分的回归方程。这时，记第i 个样本点 的预测值为 yˆ (h) ij ，则可以定义 j y 的误差平方和为 ∑= = − n i j h yij yij h 1 2 SS ( ) ( ˆ ( )) 定义Y 的误差平方和为 ∑= = p j h j h 1 SS( ) SS ( ) 当 PRESS(h) 达到最小值时，对应的 h 即为所求的成分个数。通常，总有 PRESS(h) 大于SS(h) ，而SS(h) 则小于SS(h −1)。因此，在提取成分时，总希望比 值 PRESS(h) SS(h −1) 越小越好；一般可设定限制值为 0.05，即当
PRESS(h)/SSh-1)≤1-0.05)2=0.952 时，增加成分1，有利于模型精度的提高。或者反过来说，当 PRESS(h)/SS(h-1)>0.95 时，就认为增加新的成分1，对减少方程的预测误差无明显的改善作用。 为此，定义交叉有效性为Q=1-PRESS(h)/SS(h-),这样，在建模的每一步 计算结束前，均进行交叉有效性检验，如果在第h步有Q<1-0.952=0.0985,则模 型达到精度要求，可停止提取成分：若Q2之0.0975，表示第h步提取的，成分的边际 贡献显著，应继续第h+1步计算 §2一种更简洁的计算方法 上节介绍的算法原则和推导过程的思路在目前的文献中是最为常见的。然而，还有 一种更为简洁的计算方法，即直接在E0,.,E,矩阵中提取成分4，.，4(r≤m)。要 求1，能尽可能多地携带X中的信息，同时，对因变量系统F有最大的解释能力。注 意，无需在F中提取成分得分4，这可以使计算过程大为简化，并且对算法结论的解 (1)求矩阵EFFE,最大特征值所对应的特征向量州，求得成分4=wX, 计算成分得分向量i=E”,和残差矩阵E=E。-ia,其中a=E/作。 (2)求矩阵EFFE,最大特征值所对应的特征向量W2,求得成分12=wX 计算成分得分向量，=Ew2,和残差矩阵E2=E-i,其中a=Ei/作。 ()至第r步，求矩阵EF。FE,-最大特征值所对应的特征向量w,求得成分 4,=wX,计算成分得分向量，=E,w,。 -535-
-535- 2 2 PRESS(h) SS(h −1) ≤ (1− 0.05) = 0.95 时，增加成分 ht 有利于模型精度的提高。或者反过来说，当 2 PRESS(h) SS(h −1) > 0.95 时，就认为增加新的成分 ht ，对减少方程的预测误差无明显的改善作用。 为此，定义交叉有效性为 1 PRESS( ) SS( 1) 2 Qh = − h h − ，这样，在建模的每一步 计算结束前，均进行交叉有效性检验，如果在第 h 步有 1 0.95 0.0985 2 2 Qh < − = ，则模 型达到精度要求，可停止提取成分；若 0.0975 2 Qh ≥ ，表示第h 步提取的 ht 成分的边际 贡献显著，应继续第 h +1步计算。 §2 一种更简洁的计算方法 上节介绍的算法原则和推导过程的思路在目前的文献中是最为常见的。然而，还有 一种更为简洁的计算方法，即直接在 0 1 , , E " Er− 矩阵中提取成分 r t , ,t 1 " (r ≤ m )。要 求 ht 能尽可能多地携带 X 中的信息，同时， ht 对因变量系统 F0 有最大的解释能力。注 意，无需在 F0 中提取成分得分uh ，这可以使计算过程大为简化，并且对算法结论的解 释也更为方便。 偏最小二乘法的简记算法的步骤如下： （1）求矩阵 E0 F0F0 E0 T T 最大特征值所对应的特征向量 w1，求得成分t w XT 1 = 1 ， 计算成分得分向量 1 0 1 t ˆ = E w ，和残差矩阵 T E E t 1 0 1 1 = − ˆα ，其中 2 1 0 1 1 E t ˆ t ˆ T α = 。 （2）求矩阵 E1 F0F0 E1 T T 最大特征值所对应的特征向量 w2 ，求得成分t w XT 2 = 2 ， 计算成分得分向量 2 1 2 t ˆ = E w ，和残差矩阵 T E E t 2 1 2 2 = − ˆα ，其中 2 2 1 2 2 E t ˆ t ˆ T α = 。 # （r）至第 r 步，求矩阵 −1 0 0 r−1 T T Er F F E 最大特征值所对应的特征向量 wr ，求得成分 t w XT r = r ，计算成分得分向量 r Er wr t 1 ˆ = −