偏小二乘回归提供一种多对多线性回归建模的方法,特别当两组变量的个数很多,且都存在多重相关性,而观测数据的数量(样本量)又较少时,用偏小二乘回归建立的模型具有传统的经典回归分析等方法所没有的优点。
偏小二乘回归分析在建模过程中集中了主成分分析,典型相关分析和线性回归分析方法的特点,因此在分析结果中,除了可以提供一个更为合理的回归模型外,还可以同时完成一些类似于主成分分析和典型相关分析的研究内容,提供更丰富、深入的一些信息。
例 本节采用兰纳胡德(Linnerud)给出的关于体能训练的数据进行偏小二乘回归建模。在这个数据系统中被测的样本点,是某健身俱乐部的20位中年男子。被测变量分为两组。第一组是身体特征指标X,包括:体重、腰围、脉搏。第二组变量是训练结果指标Y,包括:单杠、弯曲、跳高。原始数据见表1。
表2给出了这6个变量的简单相关系数矩阵。从相关系数矩阵可以看出,体重与腰围是正相关的;体重、腰围与脉搏负相关;而在单杠、弯曲与跳高之间是正相关的。从两组变量间的关系看,单杠、弯曲和跳高的训练成绩与体重、腰围负相关,与脉搏正相关。
表 1 体能训练数据
No |
体重( ) |
腰围( ) |
脉搏( ) |
单杠( ) |
弯曲( ) |
跳高( ) |
1 |
191 |
36 |
50 |
5 |
162 |
60 |
2 |
189 |
37 |
52 |
2 |
110 |
60 |
3 |
193 |
38 |
58 |
12 |
101 |
101 |
4 |
162 |
35 |
62 |
12 |
105 |
37 |
5 |
189 |
35 |
46 |
13 |
155 |
58 |
6 |
182 |
36 |
56 |
4 |
101 |
42 |
7 |
211 |
38 |
56 |
8 |
101 |
38 |
8 |
167 |
34 |
60 |
6 |
125 |
40 |
9 |
176 |
31 |
74 |
15 |
200 |
40 |
10 |
154 |
33 |
56 |
17 |
251 |
250 |
11 |
169 |
34 |
50 |
17 |
120 |
38 |
12 |
166 |
33 |
52 |
13 |
210 |
115 |
13 |
154 |
34 |
64 |
14 |
215 |
105 |
14 |
247 |
46 |
50 |
1 |
50 |
50 |
15 |
193 |
36 |
46 |
6 |
70 |
31 |
16 |
202 |
37 |
62 |
12 |
210 |
120 |
17 |
176 |
37 |
54 |
4 |
60 |
25 |
18 |
157 |
32 |
52 |
11 |
230 |
80 |
19 |
156 |
33 |
54 |
15 |
225 |
73 |
20 |
138 |
33 |
68 |
2 |
110 |
43 |
均值 标准差 |
178.6 24.6905 |
35.4 3.202 |
56.1 7.2104 |
9.45 5.2863 |
145.55 62.5666 |
70.3 51.2775 |
表 2 相关系数矩阵
|
|
|
|
|
|
|
|
1 |
0.8702 |
-0.3658 |
-0.3897 |
-0.4931 |
-0.2263 |
|
0.8702 |
1 |
-0.3529 |
-0.5522 |
-0.6456 |
-0.1915 |
|
-0.3658 |
-0.3529 |
1 |
0.1506 |
0.225 |
0.0349 |
|
-0.3897 |
-0.5522 |
0.1506 |
1 |
0.6957 |
0.4958 |
|
-0.4931 |
-0.6456 |
0.225 |
0.6957 |
1 |
0.6692 |
|
-0.2263 |
-0.1915 |
0.0349 |
0.4958 |
0.6692 |
1 |
MATLAB源代码:
clc,clear
ab0=load('pz.txt'); %原始数据存放在纯文本文件pz.txt中
mu=mean(ab0);sig=std(ab0); %求均值和标准差
rr=corrcoef(ab0); %求相关系数矩阵
ab=zscore(ab0); %数据标准化
a=ab(:,[1:3]);b=ab(:,[4:end]); %提出标准化后的自变量和因变量数据
[XL,YL,XS,YS,BETA,PCTVAR,MSE,stats] =plsregress(a,b);
xw=a\XS; %求自变量提出成分系数,每列对应一个成分,这里xw等于stats.W
yw=b\YS; %求因变量提出成分的系数
a_0=PCTVAR(1,:);b_0=PCTVAR(2,:);
a_1=cumsum(a_0);b_1=cumsum(b_0);
i=1;%赋初始值
%判断提出成分对的个数
while ((a_1(i)<0.9)&(a_0(i)>0.05)&(b_1(i)<0.9)&(b_0(i)>0.05))
i=i+1;
end
ncomp=i;
fprintf('%d对成分分别为:\n',ncomp);
for i=1:ncomp
fprintf('第%d对成分:\n',i);
fprintf('u%d=',i);
for k=1:3%此处为变量x的个数
fprintf('+(%f*x_%d)',xw(k,i),k);
end
fprintf('\n');
fprintf('v%d=',i);
for k=1:3%此处为变量y的个数
fprintf('+(%f*y_%d)',yw(k,i),k);
end
fprintf('\n');
end
[XL2,YL2,XS2,YS2,BETA2,PCTVAR2,MSE2,stats2] =plsregress(a,b,ncomp);
n=size(a,2); m=size(b,2);%n是自变量的个数,m是因变量的个数
beta3(1,:)=mu(n+1:end)-mu(1:n)./sig(1:n)*BETA2([2:end],:).*sig(n+1:end); %原始数据回归方程的常数项
beta3([2:n+1],:)=(1./sig(1:n))'*sig(n+1:end).*BETA2([2:end],:); %计算原始变量x1,...,xn的系数,每一列是一个回归方程
fprintf('最后得出如下回归方程:\n')
for i=1:3%此处为变量y的个数
fprintf('y%d=%f',i,beta3(1,i));
for j=1:3%此处为变量x的个数
fprintf('+(%f*x%d)',beta3(j+1,i),j);
end
fprintf('\n');
end
bar(BETA2','k') %画直方图
fprintf('因变量y的预测值:\n');
yhat=repmat(beta3(1,:),[size(a,1),1])+ab0(:,[1:n])*beta3([2:end],:); %求y1,..,ym的预测值
disp(yhat);
ymax=max([yhat;ab0(:,[n+1:end])]); %求预测值和观测值的最大值
%下面画y1,y2,y3的预测图,并画直线y=x
figure, subplot(2,2,1)
plot(yhat(:,1),ab0(:,n+1),'*',[0:ymax(1)],[0:ymax(1)],'Color','k')
legend('单杠成绩预测图',2)
subplot(2,2,2)
plot(yhat(:,2),ab0(:,n+2),'O',[0:ymax(2)],[0:ymax(2)],'Color','k')
legend('弯曲成绩预测图',2)
subplot(2,2,3)
plot(yhat(:,3),ab0(:,end),'H',[0:ymax(3)],[0:ymax(3)],'Color','k')
legend('跳高成绩预测图',2)
该MATLAB源代码结果分析:根据PCTVAR可知,前两个成分解释自变量的比率为92.13%,只需要取两对成分即可。故ncomp的值应为2,beta3为最后求解出的回归方程。
图1 回归系数的直方图