数学实验 12：回归分析

实验 12：回归分析

习题7：

在有氧锻炼中人的耗氧能力y(ml/(min ·kg))是衡量身体状况的重要指标，它可能与以下因素有关：年龄x1,体重x2(kg),1500m 跑的时间x3(min),静止时心跳速度x4（次/min ），跑步后心速x5（次/min ）.对24名40至57岁的志愿者进行了测试，结果如下表（节选），试建立耗氧能力y 与诸因素的之间的回归模型。

（1）若x1~x5中只许选择1个变量，最好的模型是什么 （2）若x1~x5中只许选择2个变量，最好的模型是什么 （3）若不限制变量的个数，最好的模型是什么 （4）对最终模型观察残查，有无异常点，若有，剔除后如何

模型建立

本题不同小问需要建立不同模型，由于专业知识所限，并且提供的数据较少，难以做出精

确符合现实情况的模型，因此这里采用最简单的线性回归法进行拟和,模型基本形式如下：

0111,m m jk j k j k m

y x x x x ββββε

≤≤=++++

+∑

L

事实上，

∑

中的项（高次项和交互项）对于本题目来讲意义不大，因为所给定的5个自变量

和因变量之间关系比较模糊，几个变量彼此之间的联系也很难说清，因此用自变量的一次线性拟和就足以适应本题的要求。但作为练习，还是将每种回归方法都使用到了，可以用于参考。 具体采用的各个模型将在下面单独说明，这里不再重复。

程序设计

由于本题需要建立多组模型，并且要在不断的调试中发现最合理的，很多命令都要在这个过程中不断使用，这里仅仅给出使用的最基本的命令。 数据 clear A=[…];

%数据矩阵，略

n=24; y=A(2,:);

%提取各个数据

x1=A(3,:);x2=A(4,:);x3=A(5,:);x4=A(6,:);x5=A(7,:); 绘制散点图（大致判断影响情况）

for i=1:5

subplot(2,3,i),plot(A(i+2,:),y,'+'),grid pause

end

序号 1 2 3 4 … 21 22

23 24 Y …

X1 44 40 44 42 … 57 54

52 50 X2 … X3 …

X4 62 62 45 40 … 58 62 48 48 X5

178 185 156 166 … 174 156 164 146

pause

单参数回归（第一问）

X=[ones(n,1),x4']; %这里检验的是自变量x4，实际操作时要分别检验x1~x5

[b,bint,r,rint,s]=regress(y',X); %回归分析程序（ =）

b,bint,s, %输出回归系数估计值、置信区间、以及统计量

rcoplot(r,rint) %残差图

Polytool(x3',y',2) %检验一元多项式回归的结果，输出交互式画面

双参数回归（第二问）：用逐步回归法找出最合理的两个变量

X5=[x1',x2',x3',x4',x5'];

stepwise(X5,y')；%利用输出的交互式画面，可以选出最佳的两个变量

XX=[x3',x1']; %当得到了最佳的两个变量后（这里假设是x3\x1）

rstool(XX,y','linear') %检验二元情况下的交互项和高次项

全部参数回归（第三问）：

X5=[x1',x2',x3',x4',x5']; %仍然用逐步回归法找出最合理的组合方式

stepwise(X5,y')

第五问要求对残差进行分析，并且剔除异常点，可以在该问得到最终模型后，采用regress得到的残差值和置信区间并根据其绘制残差图，然后再进行剔除操作重新检验。

运行结果及分析

散点图

散点图

从左上到右下的顺序为x1~x5.可以由点的分布大致看出，除了x3自变量呈现比较明显的负相关趋势以外，对于其他的各个自变量都难以直接观测出其对于因变量的影响。根据这种结

果，可以假设自变量x3(1500m 跑后心速)最直接的与锻炼耗氧能力相关，下面通过对各个自变量的单参数回归进行检验。 单参数回归

由单参数回归的结果可以证明X3(1500m 跑后心速)可以最好的反映出y(锻炼耗氧能力)的情

况。由β1置信区间可以看出，x1、x2包含0在内，即y 可能与该参数无关，所以不选择，并且两者的p 值已经明显的大于 =，则不考虑x1、x2。比较x3~x5后发现，x3的2

R -决定系数明显的大于x4、x5的，决定系数反映的是在因变量的总变化中自变量引起的那部分的比例，2

R 大说明x3自变量对因变量起的决定作用最大。并且x3的p 和s^2值也都比较小，所以最终确定x3可以最好的反映出y 的情况。

5

5.5

6

6.5

7

7.5

8

8.5

可以同之前的仅含一次项的结果进行比较，发现各个参量的置信区间都很宽，且β2的置信区间过0。可以认为二次项的引入是不重要的。

因此采用如下单参数模型描述y 是最准确地：

033y x ββ=+ 其中：0383.4438,-5.6682ββ==

双参数回归：

用stepwise 作逐步回归，部分过程和最终结果如下图：

Coeff. t-stat p-val

12

Model History

R M S E

-6

-4

-2

X X

X X X Coefficients with Error Bars

只取x3自变量的回归结果

Coeff. t-stat p-val

123

Model History

R M S E

-6

-4

-2

X X X X X Coefficients with Error Bars

取x3和x5自变量的回归结果

Coeff. t-stat p-val

123

Model History

R M S E

-6

-4

-2

X X X X X Coefficients with Error Bars

取x3和x1自变量的回归结果

根据题目要求，最终得到取双参量时的最佳结果（RMSE 参量最小）是取x3 (1500m 跑后心速)和x1（年龄）自变量。但事实上，实际的逐步回归过程在此时并没有结束，最终的最优结果是只取x3参量。这说明取x3、x1参量同只取x3相比优势并不明显。

通过rstool 命令检验二元情况下的交互项和高次项情况,下图是linear 情况下固定单参数进

行预测的结果：

5.56

6.57

7.584042444648505254

可以看到高次项和相关项的系数都非常小，说明其对于y 的影响不大。根据rmse 的结果进行比较，仍然选择linear 回归方式，即只用二元自变量的一次项。

01133y x x βββ=++ 其中：01390.8529,=-0.1870,-5.4671βββ==

全参数回归

根据以上的分析可以验证模型建立时的猜想，本题中5个字变量和y 的关系都不是很直接

的，除x3外其他变量的影响很小，所以在最终完整模型中，不再考虑高次项和交互项的影响，一方面简化模型，一方面大大节省的筛选的时间。

所以采用stepwise 命令，仅对五元变量x1~x5的一次项进行回归分析，结果如下：

Coeff. t-stat p-val

1234

Model History

R M S E

-6

-5

-4

-3

-2

-1

X X X X X Coefficients with Error Bars

R^2= F= RMSE= P=*e-7

最终取以下三个参数得到最佳回归结果：x3 (1500m 跑后心速)、x1（年龄）以及x5（跑步后心速）。但仍需要进行一般回归分析（regress ）确定常数项并观察残差，结果如下:

5101520

-8

-6-4-202

468

10R e s i d u a l s

Case Number

残差图（全部点）

可以看到10和15号数据异常，剔除，再次观察残差，结果如下

246810121416182022

-6

-4

-2

2

4

6

8

Residual Case Order Plot

R e s i d u a l s

Case Number

残差图（剔除10，15号点）

可以看到4号数据变为异常，再次剔除，结果仍然有异常点。过程从略，最终经过4次剔除，去掉5个点（4,10,15,17,23）以后，得到没有异常点的模型:

2

4

6

8

10

12

14

16

18

-4

-3-2-101

23

4R e s i d u a l s

Case Number

去掉所有异常点之后的结果

最终得到的结果整体上优越于剔除异常点之前的结果（不再粘贴结果）。但是事实上，由于数据点经过剔除不断的结果，模型最终的形式和实际统计到的24组数据的整体情况偏离越来越大，也就是说：剔除异常点虽然能够一应程度上降低其对于整体情况的干扰作用，而剔除的过程也放大了其他原本正常数据点的异常性，所以异常点可能会不断产生，但是剔除的数量增加即采样数据的减少也会削弱模型反省整体性能的能力。是一对矛盾，在数据点较少的时候尤其明显。比较科学的做法是：只进行1次或少次剔除，保证整体性，又去掉了最主要的异常点。 这里的最终结果采用剔除最初两个异常点（10，15号）后的结果，在此也附上完整数据（剔除之前）的结果，作为第3问的答案：

完整数据（第三题结果）：

0113355y x x x ββββ=+++

其中：0135118.0135,=-0.3254,-4.5694,0.1561ββββ===-

一次剔除（最终结果）：

0113355y x x x ββββ=+++

其中：0135119.4955,=-0.03623,-4.0411,0.1774ββββ===-

1500m 跑后心速、年龄以及跑步后心速三个参数最能够反映锻炼耗氧量这个重要的身体状态指标。三种心跳速度越快，说明耗氧量越大；速度越慢，即时间越长，说明耗氧量越小。

习题11：

一个医药公司的新药研究部门为了掌握一种新型止痛剂的疗效，设计了一个药物试验，给24名患有同种痛病的病人使用这种新止痛剂的以下4个剂量中的某一：2、5、7、10,（g ）,并记录每个病人病痛明显减轻的时间（min ）.为了了解新药的疗效与病人性别和血压有什么关系，试验过程中研究人员把柄热男性别集血压的低中高3档平均分配来进行测试。通过比较每个病人血压的历史数据,从低到高分成3组，分别记作,,。实验结束后，公司的记录结果见表格（略）. 请为公司建立一个模型，根据病人用药的计量、性别和血压组别，预测出服

1. 模型建立

本题共提供了三种不同的自变量，设为性别-x1，血压-x2和用药计量-x3。建立回归模型时，

应该充分考虑各个自变量对于因变量的多种影响方式，这里可能涉及到高次项和交互项。因此应该在实验过程中将所有可能的因素列于其中，分别检验，最终得到最佳的回归模型。基本公式如下:

0111,m m jk j k j k m

y x x x x ββββε

≤≤=++++

+∑

L

2. 程序设计 数据

clc,clear

x1=[0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1]; x2=[ ];

x3=[2 2 2 2 2 2 5 5 5 5 5 5 7 7 7 7 7 7 10 10 10 10 10 10];

y=[35 43 55 47 43 57 26 27 28 29 22 29 19 11 14 23 20 22 13 8 3 27 26 5];

用逐步回归法检验所有的一次、二次和交互项 XX=[x1',x2',x3',x2'.^2,x3'.^2,(x1.*x2)',(x2.*x3)',(x3.*x1)']

%注意，由于性别x1是0-1变量，所以齐高次项没有任何意义，才此舍去 stepwise(XX,y')；

根据逐步回归得到的最优模型，进行regress 回归检验，并绘制残差图 XF=[ones(24,1),x2',x3',x2'.^2,x3'.^2,(x2.*x3)',(x3.*x1)']; [b1,bint1,r1,rint1,s1]=regress(y',XF); b1,bint1,s1, rcoplot(r1,rint1);

3. 运行结果及分析

Coeff. t-stat p-val

2 4 6 8 10121416

Model History

R M S E

-50

50

X X X X X X X X Coefficients with Error Bars

上图是逐步回归得到的最优结果，其参数输出结果如下

beta betaci Coeff. t-stat p-val

x1 0 0 0 x2 0 0 0 x3 x2^2 x3^2 x1*x2 0 0 0 x2*x3 0 x3*x1

intercept rmse rsq adjrsq fstat pval

得到模型：

220132233423531y x x x x x x x ββββββ=+++++

其中：

012345=52.8084,=-7.0608,=42.5282,=0.5111,=-7.3746,=0.9551ββββββ

由多项指标可以确定以上模型的合理性：rmse 最小，R^2统计量接近于1，beta(回归系数)置信区间全部部包含0点，pval<<α=,这些都可以说明上模型在数学上是合理的。

从实际意义上来讲，以上模型也是比较有根据的。

首先，在考虑单一变量影响时，x3（用药剂量）的影响是决定性的，同性别和血压相比其显著性是可以将另外二者忽略的；

但是血压因素的2次项（x2^2）却没有被忽略，且其系数很大（40左右），说明在其他因素不变的前提下，x2-y 的函数关系为对称轴为y 轴的抛物线，也就是说所疼痛明显的减轻时间，其增大（时间延长）程度同血压的平方成正比，即血压的增加会明显的降低药效，虽然没有确切的专业知识的证明，但这也是符合常识的；

x3的二次项同样存在，且抛物线的对称轴在正半轴，如下图：

012345678910

这符合一般情况下药物的作用机制，适当用量效果最好，过多过少都会降低药效；但是注意，由于交互项的存在（x1*x3/x2*x3），性别和血压都会影响到上面抛物线的对称轴和上下位置。 两个交互项都包含了药效这个因素，可以说明，性别和血压都会通过用药计量间接的影响药效。

以上模型的实际意义和数学结果都比较合理，下面给出regress 命令和残差图的结果：

5101520

-10

-5

5

10

15

Residual Case Order Plot

R e s i d u a l s

Case Number

最后的两个数据点不合理，剔除，之后结果如下：

2

4

6

8

10

12

14

16

18

20

22

-10

-8-6-4-202

468

10Residual Case Order Plot

R e s i d u a l s

Case Number

根据第7题的分析，不再继续剔除，将这时的结果作为最终结果，参数如下：

220132233423531y x x x x x x x ββββββ=+++++