多元线性和广义线性回归
在有氧锻炼中，人的耗氧能力是衡量身体状况的重要指标，它可能与一些因素有关：年龄x1，体重x2,1500米所用的时间x3，静止时心速x4，跑步后心速x5,。先对24名志愿者进行了测试，结果如下，根据测得的数据建立耗氧能力y与诸因素之间的回归模型。
4.1 计算相关矩阵
对于多元回归，由于自变量较多，理论回归方差的选择是比较困难的，这里先计算变量间的相关系数矩阵和线性相关性检验矩阵，分析变量间的线性相关性。
data = xlsread('耗氧.xls');
X = data(:,3:7);
y = data(:,2);
[R,P] = corrcoef([y,X])
VarNames1 = {'','y','x1','x2','x3','x4','x5'};
VarNames = {'y','x1','x2','x3','x4','x5'};
newP=[VarNames1;VarNames',num2cell(P)] %转换一下格式，加入标签，便于观察
R =
1.0000 -0.3201 -0.0777 -0.8645 -0.5130 -0.4573
-0.3201 1.0000 -0.1809 0.1845 -0.1092 -0.3757
-0.0777 -0.1809 1.0000 0.1121 0.0520 0.1410
-0.8645 0.1845 0.1121 1.0000 0.6132 0.4383
-0.5130 -0.1092 0.0520 0.6132 1.0000 0.3303
-0.4573 -0.3757 0.1410 0.4383 0.3303 1.0000
P =
1.0000 0.1273 0.7181 0.0000 0.0104 0.0247
0.1273 1.0000 0.3976 0.3882 0.6116 0.0704
0.7181 0.3976 1.0000 0.6022 0.8095 0.5111
0.0000 0.3882 0.6022 1.0000 0.0014 0.0322
0.0104 0.6116 0.8095 0.0014 1.0000 0.1149
0.0247 0.0704 0.5111 0.0322 0.1149 1.0000
newP =
'' 'y' 'x1' 'x2' 'x3' 'x4' 'x5'
'y' [ 1] [0.1273] [0.7181] [5.1258e-08] [0.0104] [0.0247]
'x1' [ 0.1273] [ 1] [0.3976] [ 0.3882] [0.6116] [0.0704]
'x2' [ 0.7181] [0.3976] [ 1] [ 0.6022] [0.8095] [0.5111]
'x3' [5.1258e-08] [0.3882] [0.6022] [ 1] [0.0014] [0.0322]
'x4' [ 0.0104] [0.6116] [0.8095] [ 0.0014] [ 1] [0.1149]
'x5' [ 0.0247] [0.0704] [0.5111] [ 0.0322] [0.1149] [ 1]
从检验的p值矩阵可以看出哪些变量间的线性相关性是显著的，若p值<0.05,则认为变量间的线性相关性是显著的，反之则认为变量间的线性相关性是不显著的。从P矩阵可以看出y与x3，x4，x5的线性相关性是显著的，x3与x4，x5的线性相关性是显著的。
注意：当线性回归模型中有两个或多个自变量高度线性相关时，使用最小二乘法建立回归方程就有可能会失效，甚至会把分析引入歧途，这就是所谓的多重共线性问题。在做多元线性回归分析的时候，应做多重共线性诊断，以期得到较为合理的结果。
4.2 多元线性回归
（1）模型的建立
这里先尝试做5元线性回归，建立y关于x1，x2，....,x5的回归模型：
yi=b0+b1*xi1+b2*xi2+b3*xi3+b4*xi4+b5*xi5+a
a~N(0,aa^2), i=1,2,...,n
（2）调用LinearModel类的fit方法求解模型
调用LinearModel类的fit方法作多元线性回归，返回参数估计结果和显著性检验结果。
mmdl1 = LinearModel.fit(X,y) %5元线性回归拟合
mmdl1 =
Linear regression model:
y ~ 1 + x1 + x2 + x3 + x4 + x5
Estimated Coefficients:
Estimate SE tStat pValue
_________ ________ ________ __________
(Intercept) 121.17 17.406 6.961 1.6743e-06
x1 -0.34712 0.14353 -2.4185 0.026406
x2 -0.016719 0.087353 -0.19139 0.85036
x3 -4.2903 1.0268 -4.1784 0.00056473
x4 -0.039917 0.094237 -0.42357 0.67689
x5 -0.15866 0.078847 -2.0122 0.059407
Number of observations: 24, Error degrees of freedom: 18
Root Mean Squared Error: 2.8
R-squared: 0.816, Adjusted R-Squared 0.765
F-statistic vs. constant model: 16, p-value = 4.46e-06
根据返回值，可写出回归方程为
y=121.17-0.34713*x1-0.016719*x2-4.2903*x3-0.039917*x4-0.15866*x5
回归方差检验的p值=4.46e-6<0.05,可知在显著性水平0.05下，认为回归方差是显著的，但这并不能说明回归方程中的每一项都是显著的。参数估计表中列出了常数项和各项线性进行检验的p值，可以看出，x2，x4和x5所对应的p值大于0.05，说明在显著性水平0.05下，回归方程中的线性项x2，x4，x5都是不显著的，其中x2最不显著，其次是x4，然后x5.
（3）多重共线性诊断
多重共线性诊断的方法很多，这里用的是基于方差膨胀因子(VIF)的多重共线性诊断，
VIFi=1/(1-Ri^2) ; Ri为模型的判定系数
VIF<5时，认为不存在共线性；
5<=VIF<=10,认为存在中等程度共线性；
VIF>10,认为共线性严重，必须设法消除共线性。常用的消除共线性的方法有：去除变量，变量变换，岭回归，主成分回归
可以根据自变量的相关系数矩阵Rx计算各自变量的方差膨胀因子，自变量xi的方差膨胀因子VIFi等于Rx的逆矩阵的对角线上的第i个元素。
Rx = corrcoef(X);
VIF = diag(inv(Rx)) %inv函数求逆，diag函数求对角线上的元素
VIF =
1.5974
1.0657
2.4044
1.7686
1.6985
各自变量的方差膨胀因子均小于5，说明模型中不存在多重共线性。
（4）残差分析与异常值诊断
下面绘制残差直方图和残差正态概率图，并根据学生会残差查找异常值
figure;
subplot(1,2,1);
mmdl1.plotResiduals('histogram'); %绘制残差直方图
title('(a) 残差直方图');
xlabel('残差r');ylabel('f(r)');
subplot(1,2,2);
mmdl1.plotResiduals('probability'); %绘制残差正态该论图
title('(b) 残差正态概率图');
xlabel('残差');ylabel('概率');
Res3 = mmdl1.Residuals; %查询残差值
Res_Stu3 = Res3.Studentized; %学生化残差
id3 = find(abs(Res_Stu3)>2) %查找异常值
id3 =
10
15
从返回结果可以看出，残差基本服从正态分布，有2组数据出现异常，它们的观测序号分别为10和15.
（5）模型改进
下面去掉异常值，并将最不显著的线性项x2，x4去掉，重新建立回归模型
y=b0+b1*xi1+b3*xi3+b5*xi5+a
a~N（0，aa^2）,i=1,2,....,m
说明：在调用LinearModel类对象的fit方法作多元多项式回归时，可通过形如‘polyijk....’的参数指定多项式方差的具体形式，这里的i，j，k,.....,为取值介于0--9的整数，用来指定多项式方程中各变量的最高次数，其中，i用来指定第一个自变量的次数，j用来指定第2个自变量的次数，以此类推。
Model = 'poly10101'; %指定模型的具体形式
mmdl2 = LinearModel.fit(X,y,Model,'Exclude',id3) %去掉异常值和不显著项重新拟合
mmdl2 =
Linear regression model:
y ~ 1 + x1 + x3 + x5
Estimated Coefficients:
Estimate SE tStat pValue
________ _______ _______ __________
(Intercept) 119.5 11.81 10.118 7.4559e-09
x1 -0.36229 0.11272 -3.2141 0.0048108
x3 -4.0411 0.62858 -6.4289 4.7386e-06
x5 -0.17739 0.05977 -2.9678 0.0082426
Number of observations: 22, Error degrees of freedom: 18
Root Mean Squared Error: 2.11
R-squared: 0.862, Adjusted R-Squared 0.84
F-statistic vs. constant model: 37.6, p-value = 5.81e-08
从返回结果可以看出，剔除异常值和线性项x2，x4后的回归方程为
y=119.5-0.36229*x1-4.0411*x3-0.17739*x5
对整个回归方程进行显著性检验p=5.81e-8<0.05,说明该方程是显著的，对常数项和线性项x1，x3，x5的p值均小于0.05，说明常数项和线性项也是显著的。
4.3 多元多项式回归
虽然上面已经剔除了最不显著的线性项x2，x4，并且整个方程是显著的，但是不能认为这样的结果就是最好的回归方程，还应该尝试增加非线性项，作广义线性回归，例如二次多项式。
Model = 'poly22222'; %指定自变量x1,x2,x3,x4,x5的最高次数为2次，
mmdl3 = LinearModel.fit(X,y,Model) %完全二次多项式拟合
mmdl3 =
Linear regression model:
y ~ 1 + x1^2 + x1*x2 + x2^2 + x1*x3 + x2*x3 + x3^2 + x1*x4 + x2*x4 + x3*x4 + x4^2 + x1*x5 + x2*x5 + x3*x5 + x4*x5 + x5^2
Estimated Coefficients:
Estimate SE tStat pValue
__________ _________ ________ _________
(Intercept) 1804.1 176.67 10.211 0.0020018
x1 -26.768 3.3174 -8.069 0.0039765
x2 -16.422 1.4725 -11.153 0.0015449
x3 -7.2417 17.328 -0.41792 0.70412
x4 1.7071 1.5284 1.1169 0.34543
x5 -5.5878 1.2082 -4.6248 0.019034
x1^2 0.034031 0.02233 1.524 0.22489
x1:x2 0.18853 0.014842 12.702 0.0010526
x2^2 -0.0024412 0.0030872 -0.79075 0.48684
x1:x3 0.23808 0.21631 1.1006 0.35145
x2:x3 -0.56157 0.087918 -6.3874 0.0077704
x3^2 0.68822 0.63574 1.0826 0.35825
x1:x4 0.016786 0.015763 1.0649 0.36502
x2:x4 0.0030961 0.0058481 0.52942 0.63319
x3:x4 -0.065623 0.071279 -0.92065 0.42513
x4^2 -0.016381 0.0047701 -3.4342 0.041411
x1:x5 0.03502 0.011535 3.0359 0.056047
x2:x5 0.067888 0.0063552 10.682 0.0017537
x3:x5 0.17506 0.063871 2.7408 0.071288
x4:x5 -0.0016748 0.0056432 -0.29679 0.78599
x5^2 -0.007748 0.0027112 -2.8577 0.064697
Number of observations: 24, Error degrees of freedom: 3
Root Mean Squared Error: 0.557
R-squared: 0.999, Adjusted R-Squared 0.991
F-statistic vs. constant model: 123, p-value = 0.00104
由计算结果可知，整个回归方程显著性检验的p=0.00104<0.05,说明在显著性水平0.05下，y关于x1,x2,x3,x4,x5的完全二多项式回归方差是显著的。
上面调用fit函数分布做了5元线性回归拟合，3元线性回归拟合和完全二次多项式回归拟合得出了3个回归方程，下面绘制出这三个拟合的效果图进行对比
figure;
plot(y,'ko'); %绘制因变量y与观察序号的散点图
hold on
plot(mmdl1.predict(X),':'); %绘制5元线性回归的拟合效果图
plot(mmdl2.predict(X),'r-.'); %绘制3元线性回归的拟合效果图
plot(mmdl3.predict(X),'k'); %绘制完全二次多项式回归的拟合效果图
legend('y的原始散点','5元线性回归拟合',...
'3元线性回归拟合','完全二次回归拟合');
xlabel('y的观测序号');
ylabel('y');
横坐标是因变量的观测序号，纵坐标是因变量的取值，从图中可以看出，完全二次多项式回归拟合的效果最好，5元和3元线性回归拟合的效果差不多。
4.4 逐步回归
逐步回归法是根据自变量对因变量y的影响大小，将他们逐个进入回归方差，影响最显著的变量先引入回归方程，在引入一个变量的同时，对已引入的自变量逐个检验，将不显著的变量再从回归方差中剔除，最不显著的变量线被剔除，直到再也不能向回归方差中引入新的变量，也不能从回归方差中剔除任何一个变量为止。这样就保证了最终得到的回归方程是最优的。
LinearModel类对象的stepwise方法用来作逐步回归，这里在二次多项式的基础上，利用逐步回归方法，建立耗氧能力y与诸多因素之间的二次多项式回归模型。
mmdl4 = LinearModel.stepwise(X,y, 'poly22222') %逐步回归
figure;
plot(y,'ko');
hold on
plot(mmdl4.predict(X),':','linewidth',2)
legend('y的原始散点','逐步回归拟合');
xlabel('y的观测序号');
ylabel('y');
mmdl4 =
Linear regression model:
y ~ 1 + x1^2 + x1*x2 + x2*x3 + x1*x4 + x4^2 + x1*x5 + x2*x5 + x3*x5 + x5^2
Estimated Coefficients:
Estimate SE tStat pValue
__________ _________ _______ __________
(Intercept) 1916.6 106.48 17.999 2.2957e-08
x1 -29.485 1.6156 -18.251 2.0321e-08
x2 -15.841 0.92505 -17.124 3.553e-08
x3 3.3267 4.4986 0.7395 0.47845
x4 0.757 0.43986 1.721 0.11936
x5 -6.547 0.69061 -9.4801 5.5705e-06
x1^2 0.060353 0.0051667 11.681 9.6821e-07
x1:x2 0.17622 0.010126 17.403 3.0846e-08
x2:x3 -0.46789 0.050314 -9.2994 6.5277e-06
x1:x4 0.034115 0.0041517 8.2173 1.7857e-05
x4^2 -0.019258 0.0032306 -5.9612 0.00021239
x1:x5 0.045394 0.0050247 9.0342 8.2768e-06
x2:x5 0.063051 0.0043992 14.332 1.6742e-07
x3:x5 0.165 0.025546 6.4588 0.00011693
x5^2 -0.0052175 0.0016766 -3.1119 0.01248
Number of observations: 24, Error degrees of freedom: 9
Root Mean Squared Error: 0.521
R-squared: 0.997, Adjusted R-Squared 0.992
F-statistic vs. constant model: 201, p-value = 1.82e-09
在二次多项式回归模型的基础上，进行逐步回归，得到的回归方程
y=1916.6-29.485*x1-15.841*x2.......................................
对回归方差进行显著性检验的p=1.82e-09<0.05,说明整个回归方程是显著的。