偏最小二乘回归方法(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,…,y q}和一组自变量X={x1,…,x p}，当数据总体能够满足高斯—马尔科夫假设条件时，根据最小二乘法，有
⋂
Y=X（X T X）-1X T Y
⋂
Y将是Y的一个很好的估计量。

从这个公式容易看出，由于（X T X）必须是可逆矩阵，所以当X中的变量存在严重的多重相关性时，或者在X中的样本点数与变量个数相比显然过少时，
这个最小二乘估计都会失效并将引发一系列应用方面的困难。

考虑到这个问题，偏最小二乘回归分析提出了采用成分提取的方法。

在主成分分析中，对于单张数据表X，为了找到能最好地概括原数据的综合变量，在X中提取了第一主成分F1,使得F1中所包含的原数据变异信息可达到最大，即
Var(F1)→max
在典型相关分析中，为了从整体上研究两个数据表之间的相关关系，分别在X和Y中提取了典型成分F1和G1，它们满足
r(F1,G1) →max
F1T F1=1
G1T G1=1
在能够达到相关度最大的综合变量F1和G1之间，如果存在明显的相关关系，则可以认为，在两个数据表之间亦存在相关关系。

提取成分的做法在数据分析的方法中十分常见，除主成分、典型成分以外，常见到的还有Fisher判别法中的判别成分。

实际上，如果F是X数据表的某种成分，则意味着F是X中变量的某一线性组合F=Xa,而F作为一个综合变量，它在X中所综合提取的信息，将满足我们特殊的分析需要。

2.2 偏最小二乘回归分析的建模方法
设有q个因变量{y1,…,y q}和p个自变量{x1,…,x p},为了研究因变量与自变量的统计关系，观测n个样本点，由此构成了自变量与因变量的数据表X=【x1,…,x p】n*p和Y=【y1,…,y q】n*q。

偏最小二乘法回归分别在X与Y中提取出t1和u1(也就是说，t1是x1,…,x p的线性组合，u1是y1,…,y q的线性组合)。

在提取这两个成分时，为了回归分析的需要，有下列两个要求：
（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,…,t m，偏最小二乘法回归将通过实施Y K 对t1,…,t m的回归，然后再表达成Y K关于原变量x1,…,x p的回归方程，k=1,…,q。

3 计算方法推导
3.1 普遍采用的计算推导过程
为了数学推导方便起见，首先将数据做标准化处理。

X 经标准化处理后的数据矩阵记为E 0=(E 01,…,E 0P )n*p ,Y 经过标准化处理后的数据矩阵记为F 0=(F 01,…,F 0q )n*q 。

第一步，记t 1是E 0的第一个成分，t 1=E 0w 1, w 1是E 0的第一个轴，它是一个单位向量，即||w 1||=1；记u 1是F 0的第一个成分，u 1=F 0c 1, c 1是F 0的第一个轴，它是一个单位向量，即||c 1||=1。

如果要t 1,u 1能分别很好德代表X 与Y 中的数据变异信息，根据主成分分析原理，应该有
Var(t 1)→max Var(u 1)→max
另一方面，由于回归建模的需要，又要求t 1对u 1有最大的解释能力，由典型相关分析的思路，t 1与u 1的相关度应达到最大值，即
r(t 1,u 1)→max
因此综合起来，在偏最小二乘回归中，我们要求t 1与u 1协方差达到最大，即 Cov(t 1,u 1)=
即求解下列优化问题
max<E 0w 1,F 0C 1>
w 1T w 1=1 （3-1） c 1T c 1=1
因此，将在||w 1||=1和||c 1||=1的约束条件下，去求（w 1T E 0T F 0c 1）的最大值。

此种情况下我们就可以用拉格朗日算法求其最优解，记
s= w 1T E 0T F 0c 1-λ1(w 1T w 1-1)- λ2（c 1T c 1-1）
对s 分别求关于w 1、c 1、λ1、λ
2的偏导，并令之为零，有
=∂∂1
s
w E 0T F 0c 1-2λ1 w 1=0 （3-2）
=∂∂1
c s
F 0T E 0w 1-2λ2 c 1=0 （3-3） =∂∂1
s
λ -( w 1T w 1-1)=0 （3-4） =∂∂2
s
λ -( c 1T c 1-1)=0 （3-5）
由（3-2）~（3-5）可以推出
2λ1=2λ2= w 1T E 0T F 0c 1=<E 0w 1,F 0C 1>
记ϴ1=2λ1=2λ2= w 1T E 0T F 0c 1,所以ϴ1是优化问题的目标函数值。

把式（3-2）和式（3-3）写成
E 0T
F 0c 1= ϴ1 w 1 （3-6） F 0T E 0w 1= ϴ1 c 1 （3-7） 将式（3-7）代入式（3-6），有
E 0T
F 0F 0T E 0w 1= ϴ12 w 1 (3-8)
由式（3-8）可知，w 1是矩阵E 0T F 0F 0T E 0特征向量，对应的特征值为ϴ12，ϴ1是目标函数值，要求取得其最大值，所以w 1是对应于矩阵E 0T F 0F 0T E 0最大特征值ϴ12的单位特征向量。

求得轴w 1和c 1后，即可得到成分
t 1=E 0w 1 u 1=F 0c 1
然后，分别求E 0和F 0对t 1和u 1的回归方程
11101
*1101110,,F r t F F
Q u F E P t E T T +=+=+=
其中，2
1101/t t E P T
=，2
1
101/u u F Q T
=，向量2
1101/t t F r T
=；E 1，F 1*，F 1为回
归方程的残差矩阵。

第2成分t 2的提取，以E 1取代E 0 , F 1取代F 0 , 用上面的方法求第2个轴W 2和第2个成分t 2 ,有
1121
11
12,W E t F E F E W T
T
==
同样，E 1 , F 1分别对t 2做回归, 得到
22212221,F r t F E P t E T T +=+=
同理可推得第h 成分t h , h 的个数可以用交叉有效性原则进行, h 小于X 的秩。

如此计算下去，如果X 的秩为A ，则会有
E 0=t 1P 1T +…+t A P A T
F 0= t 1r 1T +…+t A r A T +F A
由于t 1,…,t A 均可以表示成E 01,…,E 0P 的线性组合，因此，上式可以还原成Y K = F 0K 关于X J =E 0J
的回归方程形式
Y K =b k1X 1+…+ b kP X P +F AK k=1,..,q
3.2一种简洁的计算推导过程
3.1中介绍的推导思路是最为常见的，在3.2中将介绍一种更为简洁的计算方法，即直接在E0,…,E m-1矩阵中提取成分t1,…,t m(m<p)。

要求t h能尽可能多地携带X中的信息，同时，t h对因变量系统F0有最大的解释能力。

这时无需在F0中提取成分u h，并且在迭代算法中也无需使用其残差矩阵，而始终直接用F0进行计算。

这可以使计算过程大为简化，并且对算法结论的解释也更为方便。

下面讨论成分t1,…,t m(m<=A,A=R(X))的一种新原则。

在3.1中推导偏最小二乘法回归算法时，第一步的思路是在因变量F0抽取一个成分u1=F0c1，同时在自变量E0中抽取一个成分
t1=E0w1,成分的抽取原则是max<E0w1,F0C1>。

在这个原则下得知w1，c1，u1，t1的计算方法如下：
（1）w1是矩阵E0T F0F0T E0最大特征值的特征向量，成分t1=E0w1；
（2）c1是矩阵F0T E0E0T F0最大特征值的特征向量，成分u1=F0c1；
在求得成分u1，t1以后，分别实施E0在t1上的回归，并生成残差矩阵E1，以及F0在t1上的回归，得到残差矩阵F1。

再以E1，F1取代E0，F0进行第二轮成分的提取计算，注意到成分u1,…,u m是不参加回归计算的，因此是否可以考虑不提取因变量的成分呢？
为此，用下述原则提取比变量中的成分t2是与3.1中介绍的方法，结果是完全等价的，即
由于F0K是标准化变量，所以
Cov(F0K, E0w1)=r(F0K, E0w1)
因此，该优化原则是求成分t1=E0w1，使得t1能携带尽可能多的E0变异，同时，t1对因变量F0K(k=1,…,q)的解释能力会综合达到最大值。

由于在目标函数上配上常量（n-1）2不影响其求解，即
（n-1）2∑
=
q
1
k
Cov2(F0K, E0w1)=
∑
=
q
1
k
< F0K, E0w1>2
=∑
=
q
1
k
w1T E0T F0K F0K T E0w1= w1T E0T(
∑
=
q
1
k
F0K F0K T)E0w1= w1T E0T F0F0T E0w1
为了求w1采用拉格朗日算法求解，记
s=
∑
=q
1
k < F 0K , E 0w 1>2-λ1(w 1T w 1-1)= w 1T E 0T F 0F 0T E 0w 1-λ1(w 1T w 1-1)
对s 求关于w 1和λ1的偏导，并令之为零，得
=∂∂1
s
w 2 E 0T F 0F 0T E 0w 1-2λ1 w 1=0 （3-9）
=∂∂1
s
λ -( w 1T w 1-1)=0 （3-10） 由式（3-9）可知
E 0T
F 0F 0T E 0w 1=λ1 w 1
可见，最优解w 1应是矩阵E 0T F 0F 0T E 0的一个特征向量，将它代入目标函数，并且由式（3-10）可得
∑
=q
1
k < F 0K , E 0w 1>2= w 1T E 0T F 0F 0T E 0w 1= w 1T （λ1 w 1）=λ
1
因此λ1矩阵E 0T F 0F 0T E 0的最大特征根，w 1则是其相应的特征向量。

由此可见，在新的原则下，w 1仍然是对应于E 0T F 0F 0T E 0最大特征值的特征向量，而这个新的原则完全没有提取到F 0成分u 1提取。

也就是说，t 1=E 0w 1提取可以不依赖对u 1的提取，而这种新的原则又从新的角度说明了t 1的意义。

从这个新的原则出发，对c 1,u 1的计算就可以省略。

不过，在偏最小二乘法回归的一些解释技术中，由于u 1可以较好地概括F 0中的信息，因此，它常常也是很有用。

4 应用举例
下面将通过两个具体的案例分析, 以进一步理解偏最小二乘回归的工作过程和它的特点。

4.1 应用举例一
应用举例一将采用Linnerud 给出的关于体能训练的数据进行典型相关分析。

在这个数据系统中被观测样本点，是某健身俱乐部的20位中年男子。

被观测变量分为两组，第一组是身体特征指标X ，包括：体重、腰围、脉搏；第二组变量是训练结果指标Y ，包括：单杠、弯曲、跳高。

原始数据表见表4-1。

表4-1 原始数据表
在简化算法中，对于h=1,2,3时，有
λh =
∑
=q
1
k < F 0K , E h-1w 1>2=（n-1）2
∑
=q
1
k Cov 2(F 0K , t h )
计算可得： λ1/(n-1)2=1.272426
λ2/(n-1)2=0.038763 λ3/(n-1)2=0.026655
而成分t h 的方差，u k 的方差以及t h 与 u k 相关系数的平方r 2(t h , u k )在表4-2中列出。

表4-2 Var(t h ), Var(u k )和r 2(t h , u k )
记第h 个轴是w k ，第h 个成分t k 为
t k =E h-1w h （h=1,2,3）
其中t k 亦可以表示成原自变量E 0的线性组合，即
t k =E 0w h *
则w h *=
∏
-=1
1
j h （1-w j p j T ）w h 。

表4-3给出w h *与w h 的取值。

表4-3 w h *与w h 的取值
在利用E h-1对t h 进行回归时，有回归系数向量p k ，h=1,2,3,见表4-4。

表4-4 回归系数p k
成分t k =E h-1w h 的取值见表4-5。

表4-5 t k 取值表
通过交叉验证的方法可得，之取一个成分t1时，拟合方程的预测性为最佳，不过为了后面作图和解释的方便起见，我们取两个成分t1，t2拟合预测模型。

y k=r1k t1+ r2k t2 k=1,2,3
由于成分t h可以写成自变量x j的函数，即有
t h=w h1*x1+ w h2*x2+ w h3*x3
由此可得两个成分t1，t2所建立的偏最小二乘回归模型为
y k=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 回归系数r k
所以，有
F01=-0.077E01-0.499 E02-0.132 E03
F02=-0.138E01-0.524 E02-0.085 E03
F01=-0.060E01-0.156 E02-0.007 E03
将标准化变量F ok(k=1,2,3)和E oj(j=1,2,3)分别还原成原始变量，y k(k=1,2,3)以及x j(j=1,2,3),
则回归方程为：
Y1*=47.02-0.0166x1-0.824x2-0.097x3
Y2*=612.57-0.351x1-10.52x2-0.741x3
Y3*=183.98-0.125x1-2.497x2-0.052x3
为了快速直观地观察出各个自变量在解释Y k时的作用，可以绘制回归系数图，见图4-1
图4-1 回归系数的直方图
从回归系数图中可以立刻观察到，腰围变量在解释三个回归方程时起到了极为重要的作用，然而，与单杠及弯曲相比，跳高成绩的回归方程显然不够理想，三个自变量对它的解释能力均很低。

因此有必要考虑对自变量做适当的调整。

为了考察这三个回归方程的模型精度，我们以（y ik*,y ik）为坐标值，对所有的样本点绘制预测图。

y ik*是第k个变量，第i个样本点（y ik）的预测值。

在预测直方图上，如果所有样本点都能在图的对角线附近均匀分布，则方程的拟合值与原值差异很小，这个方程的拟合效果就是满意的。

体能训练的预测图如4-2所示。

4.2 应用举例二
这是Cornell在1990年采用的一个化工方面的例子。

此后,偏最小二乘的提出者S.Wold 等人多次引用, 成为单因变量偏最小二乘回归的一个经典案例。

该例中,有个自变量x1~x7, 因变量记为y, 如表4-7所示：
表4-7 自变量和应变量对照表
表4-8给出了12种混合物关于这8个变量的观测数据。

要求建立y对x1~x7,的回归方程, 以确定7种构成元素x1~x7对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个变量的相关系数矩阵
由于q42<0.0975，所以选择h=3, 即采用t1,t2,t3三个成分做偏最小二乘回归模型, 预测效果最好。

从所得到的最终模型看,x6的回归系数值最大, 它与y正相关。

这一点符合我们的基本认识。

x5的回归系数仍然出现反常符号, 但它的取值很低, 几乎可以忽略。

从相关系数表中可以看出,x5与y的相关度不高,并且它与其他自变量之间也没有密切联系。

也就是说, x5是一个相对独立的变量, 它不能直接解释y, 甚至也很难通过其他自变量的传递作用去解释y。

因此, 它在最终模型中的回归系数非常低。

与普通最小二乘回归方程相比, 这个方程的实际含义更加清晰, 也更易于应用。