偏最小二乘回归是一种新型的多元统计数据分析方法，它与1983年由伍德和阿巴诺等人首次提出。近十年来，它在理论、方法和应用方面都得到了迅速的发展。密西根大学的弗耐尔教授称偏最小二乘回归为第二代回归分析方法。
偏最小二乘回归方法在统计应用中的重要性主要的有以下几个方面：
（1）偏最小二乘回归是一种多因变量对多自变量的回归建模方法。
（2）偏最小二乘回归可以较好地解决许多以往用普通多元回归无法解决的问题。在普通多元线形回归的应用中，我们常受到许多限制。最典型的问题就是自变量之间的多重相关性。如果采用普通的最小二乘方法，这种变量多重相关性就会严重危害参数估计，扩大模型误差，并破坏模型的稳定性。变量多重相关问题十分复杂，长期以来在理论和方法上都未给出满意的答案，这一直困扰着从事实际系统分析的工作人员。在偏最小二乘回归中开辟了一种有效的技术途径，它利用对系统中的数据信息进行分解和筛选的方式，提取对因变量的解释性最强的综合变量，辨识系统中的信息与噪声，从而更好地克服变量多重相关性在系统建模中的不良作用。（3）偏最小二乘回归之所以被称为第二代回归方法，还由于它可以实现多种数据分析方法的综合应用。
由于偏最小二乘回归在建模的同时实现了数据结构的简化，因此，可以在二维平面图上对多维数据的特性进行观察，这使得偏最小二乘回归分析的图形功能十分强大。在一次偏最小二乘回归分析计算后，不但可以得到多因变量对多自变量的回归模型，而且可以在平面图上直接观察两组变量之间的相关关系，以及观察样本点间的相似性结构。这种高维数据多个层面的可视见性，可以使数据系统的分析内容更加丰富，同时又可以对所建立的回归模型给予许多更详细深入的实际解释。
一、 偏最小二乘回归的建模策略\原理\方法 1.1建模原理
设有 q 个因变量{q y y ,...,1}和p 自变量{p x x ,...,1}。为了研究因变量和自变量的统计关系,我们观测了n 个样本点,由此构成了自变量与因变量的数据表X={p x x , (1)
和.Y={q y y ,...,1}。偏最小二乘回归分别在X 与Y 中提取出成分1t 和1u (也就是说, 1t 是
p x x ,...,1 的线形组合, 1u 是q y y ,...,1 的线形组合).在提取这两个成分时,为了回归分析的
需要,有下列两个要求:(1) 1t 和1u 应尽可能大地携带他们各自数据表中的变异信息; (2) 1t 与1u 的相关程度能够达到最大。
这两个要求表明，1t 和1u 应尽可能好的代表数据表X 和Y,同时自变量的成分 1t 对因变量的成分1u 又有最强的解释能力。在第一个成分1t 和 1u 被提取后，偏最小二乘回归分别实施X 对 1t 的回归以及 Y 对1u 的回归。如果回归方程已经达到满意的精度，则算法终止；否则,将利用 X 被1t 解释后的残余信息以及Y 被1t 解释后的残余信息进行第二轮的成分提取。如此往复，直到能达到一个较满意的精度为止。若最终对 X 共提取了 m 个成分1t ，…，m t ，偏最小二乘回归将通过实施 k y 对1t ，…，m t ， 的回归,然后再表达成k y 关于原变量x
1
，…，x m ， 的回
归方程,k=1,2,…,q 。1.2计算方法推导
为了数学推导方便起见,首先将数据做标准化处理。X 经标准化处理后的数据矩阵记为
0E =(E 01，…，E p 0)p n ?，j Y 经标准化处理后的数据矩阵记为0F =(01F ，…，q F 0)p n ?。
第一步 记1t 是0E 的第一个成分，1w 是0E 的第一个轴，它是一个单位向量，既||1w ||=1。
记1u 是0F 的第一个成分，1u =0
F c 1
。c 1
是0
F
的第一个轴，并且||
c 1
||=1。
如果要1t ，1u 能分别很好的代表X 与Y 中的数据变异信息，根据主成分分析原理，应该有
Var(1u )→max Var(1t )→max
另一方面，由于回归建模的需要，又要求1t 对1u 有很大的解释能力，有典型相关分析的思路，1t 与1u 的相关度应达到最大值，既r （1t ，1u ）→max
因此，综合起来，在偏最小二乘回归中，我们要求1t 与1u 的协方差达到最大，既
Cov(1t ，1u )=)()(11u t Var Var r(1t ，1u ) →max
正规的数学表述应该是求解下列优化问题，既
因此，将在||1w ||2
=1和||
c
1
||2
=1的约束条件下，去求(w '
1
E
'0
F
c 1
)的最大值。
如果采用拉格朗日算法，记
s=
w '
1
E
'0
F c 1
－λ1
(w
'
1
1w －1)－λ2 (c '
1c 1－1)
对s 分别求关于1w ，
c 1
，λ1
和λ
2
的偏导并令之为零，有
1
w s
??=E '00F c 1
－λ
1
21w =0 (1 -2)
1
c s
??=F '00E 1w －λ22c 1
=0 (1-3)
1
λ??s =－(w '
11w －1)=0 (1-4) 2
λ??s =－(c '
1c 1－1)=0 (1-5) 由式(1-2)~(1-5),可以推出
>=<==1010100'1'21,22c F w E c F E w λλ
记100'
1'21122c F E w ===λλθ,所以,1θ正是优化问题的目标函数值.
把式(1-2)和式(1-3)写成
11100'w c F E θ= (1-6) 11100'c w E F θ= (1-7)
将式(1-7)代入式(1-6),有
12
1100'00'w w E F F E θ= (1-8)
同理,可得
12
1100'00'c c F E E F θ= (1-9)
可见,1w 是矩阵00'
00'E F F E 的特征向量,对应的特征值为2
1θ.1θ是目标函数值,它要求取最
大值,所以, 1w 是对应于00'
00'E F F E 矩阵最大特征值的单位特征向量.而另一方面, 1c 是对应于矩阵00'
00'F E E F 最大特征值2
1θ的单位特征向量.求得轴1w 和1c 后,即可得到成分
101w E t = 101c F u =
然后,分别求0E 和0F 对1t ,1u 的三个回归方程
11'10E p t E += (1-10) 11'10*+=F q u F (1-11) 11'10F r t F += (1-12)
式中,回归系数向量是
2
11
0'1||||t t E p = (1-13) 2
11
0'1||||u u F q =
(1-14) 2
11
0'1||
||t t F r = (1-15)
而1E ,1*
F ,1F 分别是三个回归方程的残差矩阵.
第二步 用残差矩阵1E 和1F 取代0E 和0F ,然后,求第二个轴2w 和2c 以及第二个成分
2t ,2u ,有
2t =1E 2w 2u =1F 2c
211'2'222,c F E w u t >==<θ
2w 是对应于矩阵11'11'E F F E 最大特征值2
2θ的特征值, 2c 是对应于矩阵1
1'11'F E E F 最大特征值的特征向量.计算回归系数
2
22
1'2||||t t E p = 2
22
1'2||||t t F r =
因此,有回归方程
22'21E p t E += 22'21F r t F +=
如此计算下去,如果X 的秩是A ,则会有
A A p t p t E '1'10++= （1-16）
A A A F r t r t F +++='1'10 （1-17）
由于,A t t ,,1 均可以表示成p E E 001,, 的线性组合,因此,式(1-17)还可以还原成
k k F y 0*
=关于k j E x 0*=的回归方程形式，即
Ak p kp k k F x x y +++=**
11*αα k=1,2,…,q
Ak F 是残差距阵A F 的第k 列。
1.3交叉有效性
下面要讨论的问题是在现有的数据表下,如何确定更好的回归方程。在许多情形下,偏最
小二乘回归方程并不需要选用全部的成分A t t ,,1 进行回归建模,而是可以象在主成分分析一样,采用截尾的方式选择前m 个成分))(,(X A A m 秩=<,仅用这m 个后续的成分
m t t ,,1 就可以得到一个预测性较好的模型。事实上,如果后续的成分已经不能为解释0F 提
供更有意义的信息时,采用过多的成分只会破坏对统计趋势的认识,引导错误的预测结论。在多元回归分析一章中,我们曾在调整复测定系数的内容中讨论过这一观点。下面的问题是怎样来确定所应提取的成分个数。
在多元回归分析中,曾介绍过用抽样测试法来确定回归模型是否适于预测应用。我们把手
中的数据分成两部分:第一部分用于建立回归方程,求出回归系数估计量B b ,拟合值B y ?以及残差均方和2
?B σ
;再用第二部分数据作为实验点,代入刚才所求得的回归方程,由此求出2??T T y σ和。一般地,若有≈2?T σ2?B σ,则回归方程会有更好的预测效果。若 >>2?T σ2
?B σ,则回归
方程不宜用于预测。在偏最小二乘回归建模中,究竟应该选取多少个成分为宜,这可通过考察增加一个新的成分后,能否对模型的预测功能有明显的改进来考虑。采用类似于抽样测试法的工作方式,把所有n 个样本点分成两部分:第一部分除去某个样本点i 的所有样本点集合(共含n-1个样本点),用这部分样本点并使用h 个成分拟合一个回归方程;第二部分是把刚才被排除的样本点i
代入前面拟合的回归方程,得到j y 在样本点i 上的拟合值)(?i hj y -。对于每一个i =1,2,…,n,重复
上述测试,则可以定义j y 的预测误差平方和为hj PRESS ,有∑=--=n
i i hj ij hj y
y PRESS 1
2)()?( (1-18) 定义Y 的预测误差平方和为h PRESS ,有
∑==p
j hj h PRESS PRESS 1
(1-19)
显然,如果回归方程的稳健性不好,误差就很大,它对样本点的变动就会十分敏感,这种扰动误差的作用,就会加大h PRESS 的值。
另外,再采用所有的样本点,拟合含h 个成分的回归方程。这是,记第i 个样本点的预测值
为hji y
?,则可以记j y 的误差平方和为hj SS ,有∑=-=n
i hji ij hj y
y SS 1
2)?( (1-20) 定义Y 的误差平方和为h SS ,有
∑==p
j hj h SS SS 1
(1-21)
一般说来,总是有h PRESS 大于h SS ,而h SS 则总是小于1-h SS 。下面比较1-h SS 和
h PRESS 。1-h SS 是用全部样本点拟合的具有h-1个成分的方程的拟合误差; h PRESS 增加
了一个成分h t ,但却含有样本点的扰动误差。如果h 个成分的回归方程的含扰动误差能在一定程度上小于(h-1)个成分回归方程的拟合误差,则认为增加一个成分h t ,会使预测结果明显提高。因此我们希望)/(1-h h SS PRESS 的比值能越小越好。在SIMCA-P 软件中,指定2195.0)/(≤-h h SS PRESS
即
195.0-≤h h SS PRESS 时,增加成分h t 就是有益的;或者反过来说,当
195.0->h h SS PRESS 时,就认为增加新的成分h t ,对减少方程的预测误差无明显的改善
作用.另有一种等价的定义称为交叉有效性。对每一个变量k y ,定义
k
h hk
hk SS PRESS Q )1(2
1--
= (1-22)
对于全部因变量Y,成分h t 交叉有效性定义为
)
1()1(1
211--=-
=-
=∑∑h h
k
h q
k hk
h
SS PRESS SS
PRESS
Q (1-23)
用交叉有效性测量成分h t 对预测模型精度的边际贡献有如下两个尺度。
(1)
当0975.0)95.01(2
2=-≥h Q 时, h t 成分的边际贡献是显著的。显而易见, 0975.02≥h Q 与2195.0)/(<-h h SS PRESS 是完全等价的决策原则。
(2) 对于k=1,2,…,q,至少有一个k,使得
0975.02
≥h Q
这时增加成分h t ,至少使一个因变量k y 的预测模型得到显著的改善,因此,也可以考虑增加成分h t 是明显有益的。
明确了偏最小二乘回归方法的基本原理、方法及算法步骤后，我们将做实证分析。
附 录
function w=maxdet(A) %求矩阵的最大特征值 [v,d]=eig(A); [n,p]=size(d); d1=d*ones(p,1); d2=max(d1); i=find(d1==d2); w=v(:,i);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [c,m,v]=norm1(C) %对数据进行标准化处理 [n,s]=size(C); for i=1:n for j=1:s
c(i,j)=(C(i,j)-mean(C(:,j)))/sqrt(cov(C(:,j))); end end m=mean(C); for j=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最大特征值所对应的特征向量。
% 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)))';
for i=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
for s=1:px
q(:,s)=p(1,px*(s-1)+1:px*s)';
end
[n,d]=size(q);
for h=1:px
iw=eye(d);
for j=1:h-1
iw=iw*(eye(d)-w(:,j)*q(:,j)');
end
wh(:,h)=iw*w(:,h);
end
for j=1:py
zr(j,:)=(regress1(y(:,j),t))'; %求回归系数
end
for j=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); %normxy 标准化后的数据矩阵
%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)
for m=1:px
RdYtt(1,m)=sum(RdYt(1,1:m)');
end
for j=1:py
for m=1:py
Rdytt(j,m)=sum(Rdyt(j,1:m)');
end
end
for j=1:px
for m=1:px
Rd(j,m)=RdYt(1,1:m)*((w(j,1:m).^2)');
end
end
for j=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);
for m=1:px
RdXtt(1,m)=sum(RdXt(1,1:m)');
end
for j=1:px
for m=1:px
Rdxtt(j,m)=sum(Rdxt(j,1:m)');
end
end
% for j=1:px
% for m=1:px
% Rd(j,m)=RdXt(1,1:m)*((w(j,1:m).^2)');
% end
% end
% for j=1:px
% VIP(j,:)=sqrt((px*ones(1,px)./RdYtt).*Rd(j,:));
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%
function [t,u]=TU(px,py,C)
%t提取的自变量的主成分
%u 提取的因变量的主成分
c=norm1(C);
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);
B=F0'*E0*E0'*F0;
cc(:,1)=maxdet(B);
u(:,1)=F0*cc(:,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%
function drew(px,py,c)
X=c(:,1:px);
Y=c(:,px+1:px+py);
[line,l]=size(Y);
[t,q,w,wh,f0,FF]=fun717(px,py,c);
YY=X*FF(:,2:px+1,3)'+ones(line,1)*FF(:,1,3)';
subplot(1,1,1,1)
bar(f0(:,:,3))
title(' 直方图')
legend('SG','TZBFB','FHL','JK','HPZD','JPZD','TZ','ZG','GPK')
grid on
plot(YY(:,4),Y(:,4),'+');
lsline
for i=1:py
v=mod(i,4);
d=(i-v)/4;
subplot(2,2,v,d+1)
plot(YY(:,i),Y(:,i),'*');
lsline
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%
function [ Qhj,Qh,prey]=crossval7(px,py,c)
%px是自变量的个数；
%py是因量
PRESShj=zeros(px,py);
X=c(:,1:px);
Y=c(:,px+1:px+py);
x=norm1(X);
y=norm1(Y);
[line,row]=size(x);
for h=1:px
for j=1:line
newx=X;
newy=Y;
newx(j,:)=[];
newy(j,:)=[];
[t,p0,w,wh,f0,FF]=fun717(px,py,[newx,newy]);
prey(j,:,h)=X(j,:)*FF(:,2:px+1,h)'+FF(:,1,h)';
end
PRESShj(h,:)=sum((Y-prey(:,:,h)).^2);
end
PRESSh=(sum(PRESShj'))';
for h=1:px
[t1,p0,w,wh,f0,FF]=fun717(px,py,c);
prey2(:,:,h)=X(:,:)*FF(:,2:px+1,h)'+ones(line,1)*FF(:,1,h)';
SShj(h,:)=sum((Y-prey2(:,:,h)).^2);
end
SSh=(sum(SShj'))';
Qhj=ones(px-1,py)-PRESShj(2:px,:)./SShj(1:px-1,:); % 错位
Qh=ones(px-1,1)-PRESSh(2:px,1)./SSh(1:px-1,1);