第四章 回归分析.docx
《第四章 回归分析.docx》由会员分享，可在线阅读，更多相关《第四章 回归分析.docx（35页珍藏版）》请在冰点文库上搜索。
第四章回归分析
§4.2回归变量的选择与逐步回归
二、逐步回归（stepwise）
逐步回归分三种：
向前选择法，从模型中无自变量开始，根据给定的条件，每次将一个最符合条件的变量进入模型，直至所有符合条件的变量都进入模型为止。
向后剔除法，先建立全模型，根据给定的条件，每次剔除一个最不符合条件的自变量，直到回归方程不在含有不符合条件的自变量为止。
Stepwise法，即前面两种方法的结合，从含有某几个变量或没有自变量开始，根据给定的条件，将一个最符合条件的变量进入模型，再剔除新老变量中不符合条件的变量，接着再选入符合条件的变量，再剔除新老变量不符合条件的变量。
如此反复选入变、剔除变量，直到没有一个变量可选入和剔除为止。
命令：
stepwise（X,y）
stepwise（X,y,inmode）
stepwise（X,y,inmodel,penter,premove）
stepwise（X,y）
X为不包括全为1列向量n×m，n为样本容量，m为自变量个数。
y为因变量n×1列向量。
stepwise（X,y,inmode）
Inmode为逐步回归时，最初所包括的自变量。
如果n=4，如果inmode为[1,3]，则表明最初所包括的自变量为X矩阵第1列和第3列所对应的自变量。
Inmode缺失时，表明最初没有包括自变量，只包括n×1全为1的列向量。
stepwise（X,y,inmodel,penter,premove）
逐步回归时，为了了解增加和剔除变量的原则,以增加一个变量为例：
相应的P值：
当相应的P值小于等于penter时，新的变量将被引进时。
同理，删除一个变量x时：
当相应的P值大于等于premove时，相应的变量x将被删除。
如果最小的P值小于等于给定penter，或最大的P值大于等于给定的premove，则每一步都是选择最大的F值（或的P值最小的）变量引进模型。
将最小的F值（或最大的P值）对应的变量删除。
penter一定小于等于premove
缺失的情况下，penter为0.05，premove为0.1。
值得注意的是，以增加一个变量为例，新模型中F值等于新模型中增加变量对应的t值的平方，新模型中F值对应的P值等于新模型中增加变量对应t值的P值。
下面看一个例子：
序号
推销开支（x1）
实际帐目数（x2）
同类产品竞争数（x3）
地区销售潜力（x4）
建材销售量（千方）（y）
1
5.5
31
10
8
79.3
2
2.5
55
8
6
200.1
3
8
67
12
9
163.2
4
3
50
7
16
200.1
5
3
38
8
15
146
6
2.9
71
12
17
177.7
7
8
30
12
8
30.9
8
9
56
5
10
291.9
9
4
42
8
4
160
10
6.5
73
5
16
339.4
11
5.5
60
11
7
159.6
12
5
44
12
12
86.3
13
6
50
6
6
237.5
14
5
39
10
4
107.2
15
3.5
55
10
4
155
16
8
70
6
14
201.4
17
6
40
11
6
100.2
18
4
50
11
8
135.8
19
7.5
62
9
13
223.3
20
7
59
9
11
195
z=[5.500031.000010.00008.000079.3000
2.500055.00008.00006.0000200.1000
8.000067.000012.00009.0000163.2000
3.000050.00007.000016.0000200.1000
3.000038.00008.000015.0000146.0000
2.900071.000012.000017.0000177.7000
8.000030.000012.00008.000030.9000
9.000056.00005.000010.0000291.9000
4.000042.00008.00004.0000160.0000
6.500073.00005.000016.0000339.4000
5.500060.000011.00007.0000159.6000
5.000044.000012.000012.000086.3000
6.000050.00006.00006.0000237.5000
5.000039.000010.00004.0000107.2000
3.500055.000010.00004.0000155.0000
8.000070.00006.000014.0000201.4000
6.000040.000011.00006.0000100.2000
4.000050.000011.00008.0000135.8000
7.500062.00009.000013.0000223.3000
7.000059.00009.000011.0000195.0000]
x=z（:
[1:
4]）;y=z（:
5）;
stepwise（x,y）%回车得：
解释一下上面这个对话框，同四个部分组成：
左上角
右上角
中间
最低端
第一部分，彩色水平柱状图是回归系数90%的置信区间，黑色水平柱状图是回归系数95%的置信区间。
如果柱状图穿过中间虚线（横坐标为0），则在相应的显著性水平下，回归系数为0。
柱状图中间的红点，为对应回归系数的值。
第二部分，红色字体表示在原始模型上加上相应变量时，对应变量的回归系数，对应的t统计量值和对应的p值。
蓝色模型为原始模型的变量的回归系数，对应的t统计量值和对应p值。
在此例中，全为红色，说明原始模型自变量是包括只有全为1列向量。
y=c1+6.53444×x1回归系数t值：
0.7768对应的p值0.4473
y=c2+4.02871×x2回归系数t值：
0.44192对应的p值0.0003
另外两个意义也一样。
注意，四个P值，其中x3对应的P值最小，或x3的t统计量的绝对值最大。
x3用阴影注明。
最右端Nextstep:
Movex3in，即把x3引进原始模型。
AllStep可以马上显示最后的逐步回归结果。
第三部分，显示原始模型各统计量的值，本例中：
y=169.495
不包括x1,x2,x3,x4变量，RMSE是模型的标准误差72.5666
第四部分，x轴中1表示第一次逐步回归，y轴是原始模型的RMSE。
在上面的对话框中,点击NextStep出现如下对话框:
简要说明以上对话框：
第二部分，右上角，蓝色对应下面方程的t值和p值。
y=c3-23.935×x3回归系数t值：
-5.3877对应的p值0.00
红色对应的方程为在上面这个蓝色方程上，加以相应的变量所得方程。
如x1对应的方程：
y=c13+2.24413×x1+c3×x3
x1回归系数t值：
0.4094对应的p值0.6873
同理x2对应的方程：
y=c23+3.09067×x2+c3×x3
x2回归系数t值：
7.0497对应的p值0.000
第三部分，为新的原始模型的相应统计量，即：
y=387.303-23.935×x3的可决系数，F统计量值
第四部分，为第一步、第二步的模型的标准误差图。
点击Export，将输出beta、betaci、coeftab、history、in、out、stats。
这些统计量值存在workspace。
beta为原始模型的自变量回归系数，这里为-23.935。
betaci为原始模型的自变量回归系数的95%置信区间。
这里为[-33.268，-14.602]
coeftab为上图的第二部分的数据，加上对应的标准误差。
history：
包括rmse、nvars、in。
rmse为逐步回归各步对应的模型的标准误差。
这里例子是两步，第一步的标准误差为72.567，第二步的标准误差为46.125。
nvars为各步对应的原始模型的自变量个数。
本例中，第一步的变量个数为0，第二步的变量个数为1。
in说明有几个自变量在相应的原始模型，更详细的说明是哪些变量在原始模型中，1表示在原始模型，0表示不在原始模型。
本例中：
[00000;10010]其中红色表示原始模型的变量数。
第一步原始模型没有自变量，第一行全为0。
第二步原始模型中x3在原始模型中。
因此第二行第四列为1。
in表示哪个变量在原始模型中，本例中，in显示为3，即表示x3出现在原始模型中。
Out表示哪些变量不在原始模型中，本例中，out显示为[124]，即表示x1,x2,x4不在原始模型中。
Stats显示原始模型的各统计量值。
具体为上面说的对话框中第三部分值。
即包括intercept、rmse、rsq、adjrsq、fstat、pval。
上面的对话框中，注意到右上角NextStep下方显示Movex2in，即把x2引进到原始模型中。
因为对话框中，三个红色的p值x2的最小（或t值的绝对值最大）。
在上面的对话框中，再点击NextStep出现如下对话框:
对话框中第二个部分，即右上角NextStep下方显示Moveno
terms逐步回归过程全部完成。
右上角的蓝色所对应的自变量为最后模型所选择的变量。
此例中，最值模型是：
y=186.0484+3.09067x2-19.514x3
也可以对penter和premove进行设定，如：
stepwise（x,y,[1,4],0.00009,0.1）
几点说明：
最原始模型的自变量是x1和x4，在右上角相应的统计量值显示为蓝色。
红色（x2、x3）所对应的最小的P值为0.0002，但是它大于0.00009，故不能引进x3模型。
原始模型中P值最大的是x1对应的P值0.4507，大于0.1，故删除x1，NextStep下方显示Movex1out
如果改变penter和premove的大小，则出现的对话框是不一样的。
如：
stepwise（x,y,[1,4],0.001,0.1）
0.001大于0.0002，则要引要x3
如果双击对话框左上角或右上角（第一部分或第二部分）相应的柱状线或统计量值，则蓝色的线或蓝色的数字变为红色，或者，红色的线或蓝色的数字变为蓝色。
蓝色的线或蓝色的数字变为红色，表明双击的线或数字对应的变量删除；红色的线或蓝色的数字变为蓝色，表明双击的线或数字对应的变量添加。
还有，在对话框中，点击菜单Stepwise，选择scaleinputs，则是先把所有自变量标准化，再进行逐步回归。
我们验证一下：
x=z（:
[1:
4]）;y=z（:
5）;b=zscore（z）;x1=b（:
[1:
4]）;stepwise（x,y,[1,4]）再选择scaleinputs得：
stepwise（x1,y,[1,4]）得：
两对话框右上角相同，即第二部分。
点Stepwise选择AddedVariablePlot得：
上图不知道什么意思？
显然，选择变量的方法不止stepwise这一种方法，比如有m个自变量，则所有可能的模型个数是2m－1个。
选择变量的原则有:
模型的标准误差最小，Cp最小，AIC最小等，具体可参见
《统计手册》茆诗松主编科学出版社2003.1p470-p472
《线性模型引论》王松桂等编著科学出版社2004.5p159-p164
§4.3多因变量的多元线性回归
例4.3.1
xy=[0.90.80.146.630.241.477.31;12.10.157.070.461.257.42;2.96.30.337.61.022.0511.13;54.40.7812.881.612.4916.08;8.213.31.1815.861.633.1622.86;13.116.81.5618.791.933.8729.52;23.817.82.1114.632.314.534.54;34.827.83.0919.793.326.0941.22;35.422.13.5816.54.446.7847.54;4732.27.3126.227.1810.7360;62.633.29.61288.7717.6578;6855.612.8527.569.8926.8496.2;35.324.46.7610.955.5824.252.37;31.317.95.0810.156.0320.0837.77;35.224.85.5414.237.1819.2840.07;45.337.87.1420.388.822.8950.36;49.578.811.226.5610.4528.9465.33;59.7101.615.8933.1812.5139.0583.64;47.874.910.8623.911.4239.0968.16;17.740.25.117.569.0326.8141.64;3673.313.1427.28.0537.1967.3;62138.625.5436.2810.354.09103.57;9724731.3141.5314.1877.39135.8;95.227028.7940.2415.1984.02118.1;118.4233.528.0338.215.7788.39119.62;99.920526.531.5412.2986.32112.39;15128838.6146.8717.36107.94144.41;108262.231.4638.6215.1102.76130.66;162.5358.646.2152.4820.48118.84175.1;238.2454.855.8655.9626.4139.3214.44]
x=xy（:
[1:
5]）,y=xy（:
[6,7]）,c=[ones（30,1）x],b=inv（c'*c）*c'*y
b=
8.99124.3224
-0.16750.2757
0.1724-0.1341
1.70362.2313
-0.76220.9880
1.97561.0502
pf=c*b-repmat（mean（y）,30,1）;R2=[norm（pf（:
1））^2norm（pf（:
2））^2]./（30*std（y,1）.^2）
R2=
0.98040.9867
复相关系数：
sqrt（R2）
ans=
0.99010.9933
rmse=sqrt（（（30*std（y,1）.^2）-[norm（pf（:
1））^2norm（pf（:
2））^2]）./（30-6））
rmse=
6.25366.5627
或者：
stats=regstats（y（:
1）,x,'linear',{'beta','mse','tstat','fstat'}）
stats=
source:
'regstats'
beta:
[6x1double]
mse:
39.1069
tstat:
[1x1struct]
fstat:
[1x1struct]
stats.beta
ans=
8.9912
-0.1675
0.1724
1.7036
-0.7622
1.9756
sqrt（stats.mse）
ans=
6.2536
stats.tstat.pval
ans=
0.0343
0.0799
0.0066
0.0129
0.0179
0.0067
stats.fstat.pval
ans=
0
[B,BINT,R,RINT,STATS]=regress（y（:
1）,[ones（30,1）,x]）;
R2=STATS
（1）
R2=
0.9804
第二个方程的计算相似。
4.31多变量检验统计量的计算，P137
c=[ones（30,1）x],Q=y'*（eye（30）-c*inv（c'*c）*c'）*y,c1=[ones（30,1）x（:
[1,2]）];D=x（:
[3:
5]）'*（eye（30）-c1*inv（c1'*c1）*c1'）*x（:
[3:
5]）,U=det（Q）/det（Q+B2'*D*B2）
U=
0.173********358
看p65
n1=30-5-1;n2=3;F=（n1-1）/n2*（1-sqrt（U））/sqrt（U）
F=
10.71758928235678
P=1-fcdf（F,2*n2,2*（n1-1））
P=
1.985914258595400e-007
§4.4多因变量的逐步回归
xy=[0.90.80.146.630.241.477.31;12.10.157.070.461.257.42;2.96.30.337.61.022.0511.13;54.40.7812.881.612.4916.08;8.213.31.1815.861.633.1622.86;13.116.81.5618.791.933.8729.52;23.817.82.1114.632.314.534.54;34.827.83.0919.793.326.0941.22;35.422.13.5816.54.446.7847.54;4732.27.3126.227.1810.7360;62.633.29.61288.7717.6578;6855.612.8527.569.8926.8496.2;35.324.46.7610.955.5824.252.37;31.317.95.0810.156.0320.0837.77;35.224.85.5414.237.1819.2840.07;45.337.87.1420.388.822.8950.36;49.578.811.226.5610.4528.9465.33;59.7101.615.8933.1812.5139.0583.64;47.874.910.8623.911.4239.0968.16;17.740.25.117.569.0326.8141.64;3673.313.1427.28.0537.1967.3;62138.625.5436.2810.354.09103.57;9724731.3141.5314.1877.39135.8;95.227028.7940.2415.1984.02118.1;118.4233.528.0338.215.7788.39119.62;99.920526.531.5412.2986.32112.39;15128838.6146.8717.36107.94144.41;108262.231.4638.6215.1102.76130.66;162.5358.646.2152.4820.48118.84175.1;238.2454.855.8655.9626.4139.3214.44];
x=xy（:
[1:
5]）;y=xy（:
[6,7]）;
L0=（xy-repmat（mean（xy）,length（xy）,1））'*（xy-repmat（mean（xy）,length（xy）,1））;
LYY=（y-repmat（mean（y）,length（y）,1））'*（y-repmat（mean（y）,length（y）,1））;
V=zeros（1,5）;
fori=1:
5
V（i）=[L0（i,6）L0（i,7）]*inv（LYY）*[L0（i,6）L0（i,7）]'/L0（i,i）
end%得V=0.94840.96860.98570.93550.9338
n=length（xy）;p=2;P1=1-fcdf（（n-p-1）*max（V）/（p*（1-max（V）））,p,n-p-1）%得：
P1=0，引入变量x3,再L1
L1=zeros（7,7）;
fori=1:
7,j=1:
7,
L1（i,j）=L0（i,j）-L0（i,3）*L0（3,j）/L0（3,3）
end
L1（:
3）=-L0（:
3）/L0（3,3）;L1（3,:
）=L0（3,:
）/L0（3,3）;L1（3,3）=1/L0（3,3）;
formatlong
inv（[ones（30,1）x（:
3）]'*[ones（30,1）x（:
3）]）*[ones（30,1）x（:
3）]'*y%它的第三行与L1中Lxy
（1）第三行相等于，见p154页最下面,再考虑是否要引入新的变量：
L1YY=L1（[6,7],:
）;L1YY=L1YY（:
[6,7]）;
V1=zeros（1,4）;
fori=[1245]
V1（i）=[L1（i,6）L1（i,7）]*inv（L1YY）*[L1（i,6）L1（i,7）]'/L1（i,i）
end%得：
V1=0.23360.27800.98570.35850.4231
r=1;P2=1-fcdf（（n-p-r-1）*max（V1）/（p*（1-max（V1）））,p,n-p-r-1）;%P2=7.8446e-004,引入x4
L2=zeros（7,7）;
fori=1:
7,j=1:
7,
L2（i,j）=L1（i,j）-L1（i,4）*L1（4,j）/L1（4,4）
end
L2（:
4）=-L1（:
4）/L1（4,4）;L2（4,:
）=L1（4,:
）/L1（4,4）;L2（4,4）=1/L1（4,4）
inv（[ones（30,1）x（:
[3,4]）]'*[ones（30,1）x（:
[3,4]）]）*[ones（30,1）x（:
[3,4]）]'*y%它的第三、四行与L2中Lxy
（2）第三、四行相等，见p154页最下面,再考虑是否要删除变量：
L2YY=L2（[6,7],:
）;L2YY=L2YY（:
[6,7]）
V2=zeros（1,5）;
fori=1:
5
V2（i）=[L2（i,6）L2（i,7）]*in
升级会员 