第三十章偏最小二乘回归 在实际问题中,经常遇到需要研究两组多重相关变量间的相互依赖关系,并研究用 组变量(常称为自变量或预测变量)去预测另一组变量(常称为因变量或响应变量), 除了最小二乘准则下的经典多元线性回归分析(MLR),提取自变量组主成分的主成 分回归分析(PCR)等方法外,还有近年发展起来的偏最小二乘(PLS)回归方法 偏最小二乘回归提供一种多对多线性回归建模的方法,特别当两组变量的个数很 多,且都存在多重相关性,而观测数据的数量(样本量)又较少时,用偏最小二乘回归 建立的模型具有传统的经典回归分析等方法所没有的优点 偏最小二乘回归分析在建模过程中集中了主成分分析,典型相关分析和线性回归分 析方法的特点,因此在分析结果中,除了可以提供一个更为合理的回归模型外,还可以 同时完成一些类似于主成分分析和典型相关分析的研究内容,提供更丰富、深入的一些 信息。 本章介绍偏最小二乘回归分析的建模方法;通过例子从预测角度对所建立的回归模 型进行比较。 §1偏最小二乘回归 考虑P个变量y,y2…,yn与m个自变量x,x2…,xm的建模问题。偏最小二乘 回归的基本作法是首先在自变量集中提出第一成分1(1是x1,…,xm的线性组合,且 尽可能多地提取原自变量集中的变异信息);同时在因变量集中也提取第一成分l1 并要求h与相关程度达到最大。然后建立因变量y,…,y与4的回归,如果回归方 程已达到满意的精度,则算法中止。否则继续第二对成分的提取,直到能达到满意的精 度为止。若最终对自变量集提取厂个成分1,12…1,偏最小二乘回归将通过建立 n,…y与4,42…的回归式,然后再表示为y…,y与原自变量的回归方程式 即偏最小二乘回归方程式 为了方便起见,不妨假定P个因变量y,…,yn与m个自变量x,…,xm均为标准 化变量。因变量组和自变量组的n次标准化观测数据阵分别记为 y F Eo 偏最小二乘回归分析建模的具体步骤如下:
-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)分别提取两变量组的第一对成分,并使之相关性达最大 假设从两组变量分别提出第一对成分为1和l1,1是自变量集X=(x1,…,xn)的 线性组合:1=W1x+…+1nxn=w1X,1是因变量集Y=(y,…,yn)的线性组 合:4=1y+…+Vpy=v1Y。为了回归分析的需要,要求 ①t1和1各自尽可能多地提取所在变量组的变异信息 ②1和1的相关程度达到最大。 由两组变量集的标准化观测数据阵E0和F0,可以计算第一对成分的得分向量,记 为1和 xInWu E Fov 第一对成分1和1的协方差Cov(1,1)可用第一对成分的得分向量1和的内积 来计算。故而以上两个要求可化为数学上的条件极值问题: <i1>=<E01,0”1>=WEFx→max =1,rn=|n 利用 Lagrange乘数法,问题化为求单位向量w1和v1,使B1=wEFV1→最大。问 题的求解只须通过计算m×m矩阵M=E6F0FE的特征值和特征向量,且M的最大特 征值为,相应的单位特征向量就是所求的解,而v可由w计算得到n=FE (2)建立y,…,yn对1的回归及x1…,xm对t1的回归。 假定回归模型为
-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 的回归。 假定回归模型为
alte Fo=uB:+F 其中a1=(a1…a1n),B1=(1…,月n)分别是多对一的回归模型中的参数向量, E1和F1是残差阵。回归系数向量a1,B1的最小二乘估计为 a=E2/ [B=Foi,/l 称a1,B1为模型效应负荷量 (3)用残差阵E和F代替E和F重复以上步骤 记E0=1a1,F0=iB,则残差阵E1=E6-E0,F=F-F0。如果残差阵F 中元素的绝对值近似为0,则认为用第一个成分建立的回归式精度已满足需要了,可以 停止抽取成分。否则用残差阵E1和F代替E和F重复以上步骤即得 W2=(21…,w2m);v2=(n2,…,2n)分别为第二对成分的权数。而 2=E12,i2=Fv2为第二对成分的得分向量。 分别为X,Y的第二对成分的负荷量。这时有 E0=1a1 E F=iB+2月2+F2 (4)设n×m数据阵E0的秩为r≤min(n-1,m),则存在r个成分1,12… 使得 t a+e F=i1B1+…+iB+F 把l=Wk1x1+…+mxn(k=1,2…r),代入Y=1月1+…+1,B,即得P个 因变量的偏最小二乘回归方程式
-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 个 因变量的偏最小二乘回归方程式
y (5)交叉有效性检验。 一般情况下,偏最小二乘法并不需要选用存在的r个成分12…来建立回归 式,而像主成分分析一样,只选用前l个成分(l≤r),即可得到预测能力较好的回归 模型。对于建模所需提取的主成分个数l,可以通过交叉有效性检验来确定 每次舍去第i个观测(i=1,2,…,n),用余下的n-1个观测值按偏最小二乘回归 方法建模,并考虑抽取h个成分后拟合的回归式,然后把舍去的第i个观测点代入所拟 合的回归方程式,得到y(j=1,2…,P)在第i个观测点上的预测值yy(h)。对 i=1,2,…,n重复以上的验证,即得抽取h个成分时第j个因变量y,(=1,2,…,p)的 预测误差平方和为 PRESS/(h)=∑(y-jn/(h)(j=12…P) =(y1…,yn)的预测误差平方和为 PRESSO(h)=∑ PRESS/(h) 另外,再采用所有的样本点,拟合含h个成分的回归方程。这时,记第i个样本点 的预测值为y(h),则可以定义y的误差平方和为 S8(b)=∑(yn-j(h)2 定义Y的误差平方和为 SSh)=∑SS(h) 当 PRESS(h)达到最小值时,对应的h即为所求的成分个数。通常,总有 PRESS(h)大于SS(h),而SS(h)则小于SS(h-1)。因此,在提取成分时,总希望比 值 PRESS(h)/SS(h-1)越小越好;一般可设定限制值为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，即当
PRESSO(h)/SS(h-1)≤(1-0.05)2=0952 ,增加成分L有利于模型精度的提高。或者反过来说,当 PRESS(h)/SS(h-1)>0.95 时,就认为增加新的成分b,对减少方程的预测误差无明显的改善作用。 为此,定义交叉有效性为Q2=1- PRESS(h)/SSh-1),这样,在建模的每一步 计算结束前,均进行交叉有效性检验,如果在第h步有Q2<1-0.952=00985,则模 型达到精度要求,可停止提取成分:若Q2≥00975,表示第h步提取的tn成分的边际 贡献显著,应继续第h+1步计算。 §2一种更简洁的计算方法 上节介绍的算法原则和推导过程的思路在目前的文献中是最为常见的。然而,还有 一种更为简洁的计算方法,即直接在E0,…E-1矩阵中提取成分l1…,L1(r≤m)。要 求L能尽可能多地携带X中的信息,同时,t对因变量系统F有最大的解释能力。注 意,无需在F中提取成分得分,这可以使计算过程大为简化,并且对算法结论的解 释也更为方便。 偏最小二乘法的简记算法的步骤如下: (1)求矩阵 EFFE0最大特征值所对应的特征向量1,求得成分1=v1X 计算成分得分向量=E,和残差矩阵E1=E0-ia1,其中a=E/ (2)求矩阵EFFE1最大特征值所对应的特征向量v2,求得成分t2=2X 计算成分得分向量=E,和残差矩阵E2=E1-1a,其中a3=EF/ (r)至第r步,求矩阵EFFE最大特征值所对应的特征向量v,,求得成分 ln=WX,计算成分得分向量,=Ew
-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 ˆ = −