真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析.docx
《真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析.docx》由会员分享，可在线阅读，更多相关《真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析.docx（25页珍藏版）》请在冰豆网上搜索。
真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析
偏最小二乘回归是一种新型的多元统计数据分析方法，它与1983年由伍德
和阿巴诺等人首次提出。
近十年来，它在理论、方法和应用方面都得到了迅速的
发展。
密西根大学的弗耐尔教授称偏最小二乘回归为第二代回归分析方法。
偏最小二乘回归方法在统计应用中的重要性主要的有以下几个方面：
（1）偏最小二乘回归是一种多因变量对多自变量的回归建模方法。
（2）偏最小二乘回归可以较好地解决许多以往用普通多元回归无法解决的
问题。
在普通多元线形回归的应用中，我们常受到许多限制。
最典型的问题就是
自变量之间的多重相关性。
如果采用普通的最小二乘方法，这种变量多重相关性
就会严重危害参数估计，扩大模型误差，并破坏模型的稳定性。
变量多重相关问
题十分复杂，长期以来在理论和方法上都未给出满意的答案，这一直困扰着从事
实际系统分析的工作人员。
在偏最小二乘回归中开辟了一种有效的技术途径，它
利用对系统中的数据信息进行分解和筛选的方式，提取对因变量的解释性最强的
综合变量，辨识系统中的信息与噪声，从而更好地克服变量多重相关性在系统建
模中的不良作用。
（3）偏最小二乘回归之所以被称为第二代回归方法，还由于它可以实现多
种数据分析方法的综合应用。
偏最小二乘回归=多元线性回归分析+典型相关分析+主成分分析
由于偏最小二乘回归在建模的同时实现了数据结构的简化，因此，可以在二
维平面图上对多维数据的特性进行观察，这使得偏最小二乘回归分析的图形功能
十分强大。
在一次偏最小二乘回归分析计算后，不但可以得到多因变量对多自变
量的回归模型，而且可以在平面图上直接观察两组变量之间的相关关系，以及观
察样本点间的相似性结构。
这种高维数据多个层面的可视见性，可以使数据系统
的分析内容更加丰富，同时又可以对所建立的回归模型给予许多更详细深入的实
际解释。
一、偏最小二乘回归的建模策略原理方法
1
1.1建模原理
设有q个因变量{
y1,...,}和p自变量{x1,...,xp}。
为了研究因变量和自变量
y
q
的统计关系,我们观测了n个样本点,由此构成了自变量与因变量的数据表
X={x1,...,xp}和.Y={y1,...,yq}。
偏最小二乘回归分别在X与Y中提取出成分t1和
u（也就是说,t1是x1,...,xp的线形组合,u1是y1,...,yq的线形组合）.在提取这
1
两个成分时,为了回归分析的需要,有下列两个要求:
（1）t1和u1应尽可能大地携带他们各自数据表中的变异信息;
（2）t1与u1的相关程度能够达到最大。
这两个要求表明，
t和u1应尽可能好的代表数据表X和Y,同时自变量的成分
1
t对因变量的成分u1又有最强的解释能力。
1
在第一个成分
t和u1被提取后，偏最小二乘回归分别实施X对t1的回归
1
以及Y对
u的回归。
如果回归方程已经达到满意的精度，则算法终止；否则,
1
将利用X被t1解释后的残余信息以及Y被t1解释后的残余信息进行第二轮的
成分提取。
如此往复，直到能达到一个较满意的精度为止。
若最终对X共提取
了m个成分
t
1
，⋯，tm，偏最小二乘回归将通过实施
y对t1
k
，⋯，tm，的
回归,然后再表达成
y关于原变量x
k
m，的回归方程,k=1,2,⋯,q。
1，⋯，x
1.2计算方法推导
为了数学推导方便起见,首先将数据做标准化处理。
X经标准化处理后的数
据矩阵记为
E=（E
0
01
，⋯，E0）
p
n，Yj经标准化处理后的数据矩阵记为
p
F=（F01，⋯，F0q）
0
n。
p
第一步记t是
1
E的第一个成分，w1是
0
E的第一个轴，它是一个单位向量，
0
既||w1||=1。
记1。
c
u是F0的第一个成分，u1=F0c
1是
1
F的第一个轴，并且||c
1||=1。
0
2
如果要
t
1
，u1能分别很好的代表X与Y中的数据变异信息，根据主成分分
析原理，应该有
Var（u1）max
Var（t1）max
另一方面，由于回归建模的需要，又要求t1对
u有很大的解释能力，有典型相关
1
分析的思路，t1与
u的相关度应达到最大值，既
1
r（t1，
u）max
1
因此，综合起来，在偏最小二乘回归中，我们要求
t与u1的协方差达到最大，既
1
Cov（t1
，
u）=Var（）（）r（t1
t1Varu
1
1
，
u）max
1
正规的数学表述应该是求解下列优化问题，既
maxEwFc
w1,c
010
1
1
s.t
w
c
'
1
'
1
w
c
1
1
1
1
'
2=1和||c2=1的约束条件下，去求（w
因此，将在||w||1||
1E
1
'
0F
0c1）的最大
值。
如果采用拉格朗日算法，记
'
s=w
1E
'
'
1－1（w1w1－1）－
'
2（c1c
1－1）
对s分别求关于
w
1
，c
1，1
和
2的偏导并令之为零，有
s
w
1
=E
'
2w1=0（1-2）
1－1
s
c
1
=F
'
0E0w1
2c1=0（1-3）
－2
s
1
'
=－（w
1w1－1）=0（1-4）
3
s
2
'
=－（c1－1）=0（1-5）
1c
由式（1-2）~（1-5）,可以推出
''
212wEFcEw,Fc
10
2010101
记
''
122wEFc,所以,1正是优化问题的目标函数值.
1201
10
把式（1-2）和式（1-3）写成
E（1-6）
'Fcw
0
0111
F（1-7）
'Ewc
0
0111
将式（1-7）代入式（1-6）,有
2
E（1-8）
'FF'Eww
00
00111
同理,可得
2
'
F（1-9）
'EEFcc
00
00111
'FFE
'
可见,w1是矩阵E00的特征向量,对应的特征值为
0
0
2
1.1是目标函数值,它要
求取最大值,所以,
w是对应于E00矩阵最大特征值的单位特征向量.而另
0
'FF'E
10
'EEF
'
一方面,c1是对应于矩阵F00最大特征值
0
0
2
1的单位特征向量.
求得轴
w和c1后,即可得到成分
1
t1Ew
01
u1Fc
01
然后,分别求
E和F0对t1,u1的三个回归方程
0
'
E0tpE（1-10）
1
11
'
F0uqF（1-11）
11
1
'
F0trF（1-12）
1
11
式中,回归系数向量是
4
'
Et
0
1
p（1-13）
1||t||
21
'
Fu
0
1
q（1-14）
1||u||
2
1
'Ft
0
1
r（1-15）
1||t||
2
1
而E1,F1,F1分别是三个回归方程的残差矩阵.
第二步用残差矩阵
E和F1取代
1
E和F0,然后,求第二个轴w2和c2以及第
0
二个成分
t,u2,有
2
t=E1w2
2
u=F1c2
2
''
2t,uwEFc
2212
21
'
w是对应于矩阵E11最大特征值
'FFE
1
21
2
2的特征值,c2是对应于矩阵
F最大特征值的特征向量.计算回归系数
'EEF
'11
11
'
Et
1
2
p
2||t||
22
'Ft
1
2
r
2||t||
2
2
因此,有回归方程
'
E1tpE
2
22
'
F1trF
2
22
如此计算下去,如果X的秩是A,则会有
''
E0t（1-16）
p
tAp
1A
1
''
F0trtArF（1-17）
1A
1A
由于,t1,,tA均可以表示成E01,,E0p的线性组合,因此,式（1-17）还可以还原
5
成
*关于
ykF0
k
xj*E0的回归方程形式，即
k
**
yk*xxFk=1,2,⋯,q
p
k11kpAk
F是残差距阵FA的第k列。
Ak
1.3交叉有效性
下面要讨论的问题是在现有的数据表下,如何确定更好的回归方程。
在许多
情形下,偏最小二乘回归方程并不需要选用全部的成分t1,,tA进行回归建模,而
是可以象在主成分分析一样,采用截尾的方式选择前m个成分
（mA,A秩（X））,仅用这m个后续的成分t1,,tm就可以得到一个预测性较好
的模型。
事实上,如果后续的成分已经不能为解释F提供更有意义的信息时,采用
0
过多的成分只会破坏对统计趋势的认识,引导错误的预测结论。
在多元回归分析
一章中,我们曾在调整复测定系数的内容中讨论过这一观点。
下面的问题是怎样来确定所应提取的成分个数。
在多元回归分析中,曾介绍过用抽样测试法来确定回归模型是否适于预测应用。
我们把手中的数据分成两部分:
第一部分用于建立回归方程,求出回归系数估计量
b,拟合值y?
B以及残差均方和
B
2
?
B;再用第二部分数据作为实验点,代入刚才所求
得的回归方程,由此求出
2
y?
T和?
T。
一般地,若有
2
?
T
2
?
B,则回归方程会有更好的预
测效果。
若
2
?
T
2
?
B,则回归方程不宜用于预测。
在偏最小二乘回归建模中,究竟应该选取多少个成分为宜,这可通过考察增加
一个新的成分后,能否对模型的预测功能有明显的改进来考虑。
采用类似于抽样
测试法的工作方式,把所有n个样本点分成两部分:
第一部分除去某个样本点i的
所有样本点集合（共含n-1个样本点）,用这部分样本点并使用h个成分拟合一个回
归方程;第二部分是把刚才被排除的样本点i代入前面拟合的回归方程,得到yj在
样本点i上的拟合值y?
。
对于每一个i=1,2,⋯,n,重复上述测试,则可以定义
hj（i）
y
j
的预测误差平方和为
PRESS,有
hj
n
2
PRESShj（yy?
（）（1-18）
ijhji）
i1
6
定义Y的预测误差平方和为
PRESS,有
h
p
PRESShPRESS
hj
j1
（1-19）
显然,如果回归方程的稳健性不好,误差就很大,它对样本点的变动就会十分敏感,
这种扰动误差的作用,就会加大
PRESS的值。
h
另外,再采用所有的样本点,拟合含h个成分的回归方程。
这是,记第i个样本
点的预测值为
y?
则可以记
hji
y的误差平方和为SShj,有
j
n
2
SS
hj（yy?
）（1-20）
ijhji
i1
定义Y的误差平方和为
SS,有
h
p
SShSS
hj
j1
（1-21）
一般说来,总是有
PRESS大于SSh,而SSh则总是小于SSh1。
下面比较SSh1和
h
PRESS。
SSh1是用全部样本点拟合的具有h-1个成分的方程的拟合误差;
h
PRESS增加了一个成分th,但却含有样本点的扰动误差。
如果h个成分的回归方
h
程的含扰动误差能在一定程度上小于（h-1）个成分回归方程的拟合误差,则认为增
加一个成分
t,会使预测结果明显提高。
因此我们希望（PRESSh/SSh1）的比值能
h
越小越好。
在SIMCA-P软件中,指定
（PRESSh/SSh
2
1）0.95
即PRESSh0.95SSh1时,增加成分th就是有益的;或者反过来说,当
PRESSh0.95SSh时,就认为增加新的成分th,对减少方程的预测误差无明显
1
的改善作用.
另有一种等价的定义称为交叉有效性。
对每一个变量
y,定义
k
PRESS
21hk（1-22）Q
hkSS
（h1）k
7
对于全部因变量Y,成分
t交叉有效性定义为
h
q
PRESS
hk
PRESS
21k11h
Q（1-23）
hSS
（h1）
SS
（h1）k
用交叉有效性测量成分
t对预测模型精度的边际贡献有如下两个尺度。
h
22
（1）当Q（10.95）0.0975时,th成分的边际贡献是显著的。
显而易
h
2
见,Q0.0975与
h
2
（PRESSh/SSh1）0.95是完全等价的决策原则。
（2）对于k=1,2,⋯,q,至少有一个k,使得
2
Q
h
1.4
这时增加成分
t,至少使一个因变量yk的预测模型得到显著的改善,因此,也
h
可以考虑增加成分
t是明显有益的。
h
明确了偏最小二乘回归方法的基本原理、方法及算法步骤后，我们将做
实证分析。
附录
functionw=maxdet（A）
%求矩阵的最大特征值
[v,d]=eig（A）;
[n,p]=size（d）;
d1=d*ones（p,1）;
d2=max（d1）;
8
i=find（d1==d2）;
w=v（:
i）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[c,m,v]=norm1（C）
%对数据进行标准化处理
[n,s]=size（C）;
fori=1:
n
forj=1:
s
c（i,j）=（C（i,j）-mean（C（:
j）））/sqrt（cov（C（:
j）））;
end
end
m=mean（C）;
forj=1:
s
v（1,j）=sqrt（cov（C（:
j）））;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[t,q,w,wh,f0,FF]=fun717（px,py,C）
%px自变量的输入个数
%py输入因变量的个数。
%C输入的自变量和因变量组成的矩阵
%t提取的主成分
%q为回归系数。
%w最大特征值所对应的特征向量。
9
%wh处理后的特征向量
%f0回归的标准化的方程系数
%FF原始变量的回归方程的系数
c=norm1（C）;%norm1为标准化函数
y=c（:
px+1:
px+py）;%截取标准化的因变量
E0=c（:
1:
px）;
F0=c（:
px+1:
px+py）;
A=E0'*F0*F0'*E0;
w（:
1）=maxdet（A）;%求最大特征向量
t（:
1）=E0*w（:
1）;%提取主成分
E（:
1:
px）=E0-t（:
1）*（E0'*t（:
1）/（t（:
1）'*t（:
1）））';
%获得回归系数
p（:
1:
px）=（E0'*t（:
1）/（t（:
1）'*t（:
1）））';
fori=0:
px-2
B（:
px*i+1:
px*i+px）=E（:
px*i+1:
px*i+px）'*F0*F0'*E（:
px*i+1:
px*i+px）
;
w（:
i+2）=maxdet（B（:
px*i+1:
px*i+px））;
%maxdet为求最大特征值的函数
t（:
i+2）=E（:
px*i+1:
px*i+px）*w（:
i+2）;
p（:
px*i+px+1:
px*i+2*px）=（E（:
px*i+1:
px*i+px）'*t（:
i+2）/（t（:
i+2）'*t（
:
i+2）））';
E（:
px*i+px+1:
px*i+2*px）=E（:
px*i+1:
px*i+px）-t（:
i+2）*（E（:
px*i+1:
px*
i+px）'*t（:
i+2）/（t（:
i+2）'*t（:
i+2）））';
end
fors=1:
px
10
q（:
s）=p（1,px*（s-1）+1:
px*s）';
end
[n,d]=size（q）;
forh=1:
px
iw=eye（d）;
forj=1:
h-1
iw=iw*（eye（d）-w（:
j）*q（:
j）'）;
end
wh（:
h）=iw*w（:
h）;
end
forj=1:
py
zr（j,:
）=（regress1（y（:
j）,t））';%求回归系数
end
forj=1:
px
fori=1:
py%
生成标准化变量的方程的系数矩阵
w1=wh（:
1:
j）;
zr1=（zr（i,1:
j））';
f0（i,:
j）=（w1*zr1）';
end
[normxy,meanxy,covxy]=norm1（C）;%no
rmxy标准化后的数据矩阵
11
%meanxy每一列的均值
%covxy每一列的方差
ccxx=ones（py,1）*meanxy（1,1:
px）;
ccy=（covxy（1,px+1:
px+py））'*ones（1,px）;
ccx=ones（py,1）*（covxy（1,1:
px））;
ff=ccy.*f0（:
:
j）./ccx;
fff=-（sum（（ccy.*ccxx.*f0（:
:
j）./ccx）'）-meanxy（1,px+1:
px+py））';
FF（:
:
j）=[fff,ff];%生成
原始变量方程的常数项和系数矩阵
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[r,Rdyt,RdYt,RdYtt,Rdytt,VIP]=fun8y（px,py,c）
X=c（:
1:
px）;
Y=c（:
px+1:
px+py）;
x=norm1（X）;
y=norm1（Y）;
[t,q,w]=fun717（px,py,[X,Y]）;
r1=corrcoef（[y,t]）;
r=r1（py+1:
px+py,1:
py）';
Rdyt=r.^2;
RdYt=mean（Rdyt）
form=1:
px
RdYtt（1,m）=sum（RdYt（1,1:
m）'）;
end
forj=1:
py
form=1:
py
Rdytt（j,m）=sum（Rdyt（j,1:
m）'）;
12
end
end
forj=1:
px
form=1:
px
Rd（j,m）=RdYt（1,1:
m）*（（w（j,1:
m）.^2）'）;
end
end
forj=1:
px
VIP（j,:
）=sqrt（（px*ones（1,px）./RdYtt）.*Rd（j,:
））;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[r,Rdxt,RdXt,RdXtt,Rdxtt]=fun8x（px,py,c）
X=c（:
1:
px）;
Y=c（:
px+1:
px+py）;
x=norm1（X）;
y=norm1（Y）;
[t,q,w]=fun717（px,py,[X,Y]）;
r1=corrcoef（[x,t]）;
r=r1（px+1:
px+px,1:
px）';
Rdxt=r.^2;
RdXt=mean（Rdxt）;
form=1:
px
RdXtt（1,m）=sum（RdXt（1,1:
m）'）;
end
forj=1:
px
form=1:
px
Rdxtt（j,m）=sum（Rdxt（j,1:
m）'）;
end
end
13
%forj=1:
px
%form=1:
px
%Rd（j,m）=RdXt（1,
升级会员 