偏最小二乘回归是一种新型的多元统计数据分析方法,它与

偏最小二乘回归是一种新型的多元统计数据分析方法，它与1983年由伍德和阿巴诺等人

首次提出。近十年来，它在理论、方法和应用方面都得到了迅速的发展。密西根大学的弗耐尔教授称偏最小二乘回归为第二代回归分析方法。

偏最小二乘回归方法在统计应用中的重要性主要的有以下几个方面：

（1）偏最小二乘回归是一种多因变量对多自变量的回归建模方法。

（2）偏最小二乘回归可以较好地解决许多以往用普通多元回归无法解决的问题。在普通多元线形回归的应用中，我们常受到许多限制。最典型的问题就是自变量之间的多重相关性。如果采用普通的最小二乘方法，这种变量多重相关性就会严重危害参数估计，扩大模型误差，并破坏模型的稳定性。变量多重相关问题十分复杂，长期以来在理论和方法上都未给出满意的答案，这一直困扰着从事实际系统分析的工作人员。在偏最小二乘回归中开辟了一种有效的技术途径，它利用对系统中的数据信息进行分解和筛选的方式，提取对因变量的解释性最强的综合变量，辨识系统中的信息与噪声，从而更好地克服变量多重相关性在系统建模中的不良作用。

（3）偏最小二乘回归之所以被称为第二代回归方法，还由于它可以实现多种数据分析方法的综合应用。

偏最小二乘回归=多元线性回归分析+典型相关分析+主成分分析 由于偏最小二乘回归在建模的同时实现了数据结构的简化，因此，可以在二维平面图上对多维数据的特性进行观察，这使得偏最小二乘回归分析的图形功能十分强大。在一次偏最小二乘回归分析计算后，不但可以得到多因变量对多自变量的回归模型，而且可以在平面图上直接观察两组变量之间的相关关系，以及观察样本点间的相似性结构。这种高维数据多个层面的可视见性，可以使数据系统的分析内容更加丰富，同时又可以对所建立的回归模型给予许多更详细深入的实际解释。

一、 偏最小二乘回归的建模策略\\原理\\方法

1.1建模原理

设有 q个因变量{y1,...,yq}和p自变量{x1,...,xp}。为了研究因变量和自变量的统计关系,

我们观测了n个样本点,由此构成了自变量与因变量的数据表X={x1,...,xp}和.Y={y1,...,yq}。偏最小二乘回归分别在X与Y中提取出成分t1 和u1 (也就是说, t1是x1,...,xp 的线形组合,

u1是y1,...,yq 的线形组合).在提取这两个成分时,为了回归分析的需要,有下列两个要求:

(1) t1和u1应尽可能大地携带他们各自数据表中的变异信息; (2) t1与u1 的相关程度能够达到最大。

这两个要求表明，t1和u1 应尽可能好的代表数据表X和Y,同时自变量的成分 t1对因变量的成分u1 又有最强的解释能力。

在第一个成分t1和 u1被提取后，偏最小二乘回归分别实施X 对 t1的回归以及 Y对u1 的回归。如果回归方程已经达到满意的精度，则算法终止；否则,将利用 X被t1解释后的残余信息以及Y 被t1 解释后的残余信息进行第二轮的成分提取。如此往复，直到能达到一个较满意的精度为止。若最终对 X共提取了 m个成分t1yk偏最小二乘回归将通过实施

，…，tm，

对t1，…，tm，

的回归,然后再表达成yk关于原变量

x1，…，

xm，

的回归方

程,k=1,2,…,q 。

1.2计算方法推导

为了数学推导方便起见,首先将数据做标准化处理。X 经标准化处理后的数据矩阵记为

E0=(E01，…，E0p)np，Yj经标准化处理后的数据矩阵记为F0=(F01，…，F0q)np。

第一步 记t1是E0的第一个成分，w1是E0的第一个轴，它是一个单位向量，既||w1||=1。 记u1是F0的第一个成分，u1=F0如果要t1有

Var(u1)max Var(t1)max

另一方面，由于回归建模的需要，又要求t1对u1有很大的解释能力，有典型相关分析的思路，

c。c 是F的第一个轴，并且||c||=1。

1101，

u1能分别很好的代表X与Y中的数据变异信息，根据主成分分析原理，应该

t1与u1的相关度应达到最大值，既

r（t1，u1）max

因此，综合起来，在偏最小二乘回归中，我们要求t1与u1的协方差达到最大，既

Cov(t1，

u1)=Var(t1)Var(u1)r(t1u1) max

，

正规的数学表述应该是求解下列优化问题，既

maxE0w1,F0c1 w1,c1'w1w11 s.t 'cc111因此，将在||w1||2=1和||c1||2=1的约束条件下，去求(w1如果采用拉格朗日算法，记

s=w1对s分别求关于w1，c1，1''EF0'0c)的最大值。

1E'0F0c－ (w11'1w1－1)－2 (c1c1－1)

'和2的偏导并令之为零，有

s=E'0F0 w1

c －211w1=0 (1 -2)

s=F'0E0w12 －2c1c=0 (1-3)

1

s'=－(w1w1－1)=0 (1-4) 1s'=－(c1c1－1)=0 (1-5) 2

由式(1-2)~(1-5),可以推出

2122w'1E'0F0c1E0w1,F0c1

记12122w'1E'0F0c1,所以,1正是优化问题的目标函数值. 把式(1-2)和式(1-3)写成

E'0F0c11w1 (1-6)

F'0E0w11c1 (1-7)

将式(1-7)代入式(1-6),有

E'0F0F'0E0w11w1 (1-8)

同理,可得

2F'0E0E'0F0c11c1 (1-9)

可见,w1是矩阵E'0F0F'0E0的特征向量,对应的特征值为1.1是目标函数值,它要求取最大值,所以, w1是对应于E'0F0F'0E0矩阵最大特征值的单位特征向量.而另一方面, c1是对应于矩阵

2F'0E0E'0F0最大特征值1的单位特征向量.

22求得轴w1和c1后,即可得到成分

t1E0w1 u1F0c1 然后,分别求E0和F0对t1,u1的三个回归方程

E0t1p'1E1 (1-10) F0u1q'1F1 (1-11) F0t1r'1F1 (1-12) 式中,回归系数向量是

E'0t1 p1 (1-13)

||t1||2F'0u1 q1 (1-14) 2||u1||F'0t1 (1-15) r1||t1||2而E1,F1,F1分别是三个回归方程的残差矩阵.

第二步 用残差矩阵E1和F1取代E0和F0,然后,求第二个轴w2和c2以及第二个成分

t2,u2,有

t2=E1w2 u2=F1c2

2t2,u2w'2E'1F1c2

2w2是对应于矩阵E'1F1F'1E1最大特征值2的特征值, c2是对应于矩阵F'1E1E'1F1最大特

征值的特征向量.计算回归系数

E'1t2 p2 2||t2||F'1t2 r2 2||t2||因此,有回归方程

E1t2p'2E2 F1t2r'2F2 如此计算下去,如果X的秩是A,则会有

E0t1p'1tAp'A （1-16）

F0t1r'1tAr'AFA （1-17）

由于,t1,,tA均可以表示成E01,,E0p的线性组合,因此,式(1-17)还可以还原成ykF0k关于xj*E0k的回归方程形式，即

*yk*k1x1kpx*pFAk k=1,2,…,q

FAk是残差距阵FA的第k列。 1.3交叉有效性

下面要讨论的问题是在现有的数据表下,如何确定更好的回归方程。在许多情形下,偏最小二乘回归方程并不需要选用全部的成分t1,,tA进行回归建模,而是可以象在主成分分析一样,采用截尾的方式选择前m 个成分(mA,A秩(X)),仅用这m 个后续的成分t1,,tm就可以得到一个预测性较好的模型。事实上,如果后续的成分已经不能为解释F0提供更有意义的信息时,采用过多的成分只会破坏对统计趋势的认识,引导错误的预测结论。在多元回归分析一章中,

*

我们曾在调整复测定系数的内容中讨论过这一观点。 下面的问题是怎样来确定所应提取的成分个数。

在多元回归分析中,曾介绍过用抽样测试法来确定回归模型是否适于预测应用。我们把手中

ˆB以及残差的数据分成两部分:第一部分用于建立回归方程,求出回归系数估计量bB,拟合值y22ˆBˆT和ˆT均方和;再用第二部分数据作为实验点,代入刚才所求得的回归方程,由此求出y。一2222ˆTˆBˆTˆB般地,若有,则回归方程会有更好的预测效果。若 ,则回归方程不宜用于预测。

在偏最小二乘回归建模中,究竟应该选取多少个成分为宜,这可通过考察增加一个新的成分后,能否对模型的预测功能有明显的改进来考虑。采用类似于抽样测试法的工作方式,把所有n个样本点分成两部分:第一部分除去某个样本点i的所有样本点集合(共含n-1个样本点),用这部分样本点并使用h个成分拟合一个回归方程;第二部分是把刚才被排除的样本点i代入前面

ˆhj(i)。对于每一个i=1,2,…,n,重复上述测试,拟合的回归方程,得到yj在样本点i上的拟合值y则可以定义yj的预测误差平方和为PRESShj,有

ˆhj(i))2 (1-18) PRESShj(yijyi1n定义Y 的预测误差平方和为PRESSh,有

PRESShPRESShj (1-19)

j1p显然,如果回归方程的稳健性不好,误差就很大,它对样本点的变动就会十分敏感,这种扰动误差的作用,就会加大PRESSh的值。

另外,再采用所有的样本点,拟合含h 个成分的回归方程。这是,记第i个样本点的预测值为

ˆhji,则可以记yj的误差平方和为SShj,有 y

ˆhji)2 (1-20) SShj(yijyi1n定义Y的误差平方和为SSh,有

SShSShj (1-21)

j1p一般说来,总是有PRESSh大于SSh,而SSh则总是小于SSh1。下面比较SSh1和PRESSh。

SSh1是用全部样本点拟合的具有h-1个成分的方程的拟合误差; PRESSh增加了一个成分th,但却含有样本点的扰动误差。如果h个成分的回归方程的含扰动误差能在一定程度上小于(h-1)个成分回归方程的拟合误差,则认为增加一个成分th,会使预测结果明显提高。因此我们希望

(PRESSh/SSh1)的比值能越小越好。在SIMCA-P软件中,指定

(PRESSh/SSh1)0.952

即PRESSh0.95SSh1时,增加成分th就是有益的;或者反过来说,当PRESSh0.95SSh1时,就认为增加新的成分th,对减少方程的预测误差无明显的改善作用.

另有一种等价的定义称为交叉有效性。对每一个变量yk,定义

2 Qhk1PRESShk (1-22)

SS(h1)k对于全部因变量Y,成分th交叉有效性定义为

Q12hPRESSk1qhk1SSPRESSh (1-23)

SS(h1)(h1)k用交叉有效性测量成分th对预测模型精度的边际贡献有如下两个尺度。 (1)

2 当Qh(10.952)0.0975时, th成分的边际贡献是显著的。显而易见, 2Qh0.0975与(PRESSh/SSh1)0.952是完全等价的决策原则。

(2) 对于k=1,2,…,q,至少有一个k,使得

2Qh0.0975

这时增加成分th,至少使一个因变量yk的预测模型得到显著的改善,因此,也可以考虑增加成分th是明显有益的。

明确了偏最小二乘回归方法的基本原理、方法及算法步骤后，我们将做实证分析。

附 录

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);