第六章 线性回归分析 在许多实际问题中，经常会遇到需要同时考虑几个变量的情况，例如，在电路中会遇到电压、电流及电阻的关系，在炼钢过程中会遇到钢水中的碳含量与钢材的物理性能（如强度、延伸率等）之间的关系，医学上经常测量人的身高、体重，研究人的血压与年龄的关系，在制定销售策略时会考虑商品的单价与销售量之间的关系等等．

在微积分中，我们曾提到变量之间的三种关系：无关、相关及函数关系．如电路中会遇到电压、电流及电阻的关系是一种由欧姆定律所描述的确定性的函数关系．而人的身高与体重、血压与年龄的关系则属于相关的关系，也就是说，这些变量相互间有制约关系但又不是确定性的关系．本章则主要利用回归分析的手段，研究这后一种变量之间相关关系的统计规律性．

在回归分析中，把变量分为两类，一类是因变量，它们通常是实际问题中所关心的一些指标，通常用Y表示，而影响因变量取值的另一类变量则称为自变量，又称为解释变量，通常用X1，X2，……，Xp表示．

6-1 一元线性回归 6-1-1 直线回归的概念 首先看一个例子． 例1 这是研究标准混合肥量与马铃薯产量的效应的实例．将硫酸铵、磷酸钙与硫酸钾按1：3：1的比例组成混合肥料，分别按每英亩0，4，8，与12英担（1英担为112磅）给标号为1，2，3，4的四块地施肥，由于自然变异，在给定的肥料量下，马铃薯的产量将因地而异．结果是，1号的地的产量是8.34吨/英亩，2号地的产量是8.89吨/英亩，3号地的产量是9.16吨/英亩，4号地的产量是9.50吨/英亩．

首先绘制单位施肥量X与单位产量Y的数据散点图 > X=c(0,4,8,12);Y=c(8. 34,8. 89,9. 16,9 首先绘制单位施肥量X与单位产量Y的数据散点图 > X=c(0,4,8,12);Y=c(8.34,8.89,9.16,9.50) > plot(X,Y)

显然，产量是随着肥料量的增加呈现上升趋势，并近于直线，从而可认为Y与X之间的关系基本上是线性的，但并非4个点恰好全都在一条件直线上，这与两变量之间的严格对应的函数关系不同，称之为直线回归．可作一条直线，使得这4点到该直线的垂直偏离（称之为残差）的积累最小，称该直线为回归直线．

式中的β0与β1是决定直线的两个系数，β0是回归直线在Y轴上的截距，β1称为回归系数，即直线的斜率. 为区别于一般的直线方程，称该直线的方程为直线回归方程．记之为 式中的β0与β1是决定直线的两个系数，β0是回归直线在Y轴上的截距，β1称为回归系数，即直线的斜率.

β1>0，表明Y值随着X的增大而增大，β1<0，表明Y值随着X的增大而减小，β1=0，表明回归直线与X轴平行，即Y与X无直线关系．

5-1-2 直线回归方程的求法及假设检验 根据调查或实验，得到解释变量与因变量的n对对应值 X： x1 x2 x3 …… xn Y： y1 y2 y3 …… yn 设直线回归方程是，将X的n个值代入回归方程，得到n个估计值，构造量Q如下

刚才提及，各散点到回归直线的垂直偏离的积累最小，也就是说β0与β1的应使上述Q值最小．求Q的最小值如下

实际分析可知，β0与β1的适当取值，能够使得Q值最小，因而在这一唯一的驻点处Q取得最小值，从而直线回归方程中的系数β0与β1计算按上式即可，只不过这是一个非常繁锁的过程，在后面我们将对例1用R软件处理．

此时计算所得的系数β0与β1仅仅是根据含有n对数据这样一个样本而得的结果，那么总体的情况如何呢？比如，如果计算所得的系数β1不等于零，是否意谓着总体的系数β1也不为零呢，对于系数β0也有同样的问题．于是，用上述方法计算而得的系数β0与β1的值仅仅是对总体情况的一个点估计，这种估计有没有统计意义则需要进行假设检验．回归分析中的假设检验包括回归方程总的假设检验与各个回归系数的假设检验．

R软件的命令函数lm( )提供系数β0与β1的计算及假设检验功能，下面用例1的数据操作如下： > x=c(0,4,8,12);y=c(8.34,8.89,9.16,9.50) > lm.sol<-lm(y~1+x) > summary(lm.sol) 这里summary( )的功能是将lm( )的计算结果提取到显示终端．

显示如下： Call: lm(formula = y ~ 1 + x) Residuals: 1 2 3 4 -7.000e-02 1.050e-01 1.735e-17 -3.500e-02 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.41000 0.07748 108.550 8.49e-05 *** X 0.09375 0.01035 9.055 0.0120 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.0926 on 2 degrees of freedom Multiple R-Squared: 0.9762, Adjusted R-squared: 0.9643 F-statistic: 82 on 1 and 2 DF, p-value: 0.01198

显示共有四个部分 第一部分（Call）列出了相应的回归模型公式，其模型公式y ~ 1 + x表示回归方程是． 第二部分（Residuals）列出了4对数据各自的残差即(i=1,2,3,4)．

第三部分（Coefficients）中，Estimate表示回归方程中系数β0与β1的估计值，在该例中β0=8. 41000，β1=0 第三部分（Coefficients）中，Estimate表示回归方程中系数β0与β1的估计值，在该例中β0=8.41000，β1=0.09375，Std. Error表示回归系数的标准差，t value为t值，Pr(>|t|)表示回归系数假设检验双尾概率，右面跟的就是星号属于显著性标记，三个星号表明极其显著，二个星号说明高度显著，一个星号说明显著，点号说明不太显著． 第四部分中，Residual standard error表示残差的标准差，Multiple R-Squared表示相关系数的平方，计算结果的最后一行是对回归方程作F检验的概率．

从结果来看，例1的回归方程通过了回归系数的假设检验与回归方程的检验，由此得到回归方程是．

例2 现调查了生产某种产品的个同类企业的月产量与生产成本的统计资料如下： 1 2 3 4 5 6 7 8 月产量 （千吨） 1.2 2.0 3.1 3.8 5.0 6.1 7.2 8.0 生产成本 （万元） 62 86 80 110 115 132 135 160 试作回归分析．

> x=c(1.2, 2.0, 3.1, 3.8, 5.0, 6.1, 7.2, 8.0) > y=c(62, 86, 80, 110, 115, 132, 135, 160) > lm.sol<-lm(y~1+x) > summary(lm.sol)

结果输出 Call: lm(formula = y ~ 1 + x) Residuals: Min 1Q Median 3Q Max -11.301 -5.892 0.604 6.353 9.672 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 51.323 6.755 7.598 0.000271 *** x 12.896 1.326 9.723 6.8e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 8.587 on 6 degrees of freedom Multiple R-Squared: 0.9403, Adjusted R-squared: 0.9304 F-statistic: 94.55 on 1 and 6 DF, p-value: 6.795e-05

从结果来看，回归方程是 并且通过了回归系数的假设检验与回归方程的检验．其中51.323相当于固定成本，12.896相当于边际成本．

6-1-3 线性预测与控制 线性回归的主要应用，是由相应的X信息预报Y，即测得一个新的X的值后，预报该总体中新的Y的值，确切地说，是根据已有的X的值，得到Y期望值的点估计值，也可得到Y期望值的区间估计．

> predict(object,newdata=data.frame, level=0.95 ) R软件中用于线性预测的命令函数是predict( )，其使用格式如下： > predict(object,newdata=data.frame, level=0.95 ) 其中，object是由lm构成的对象，newdata是预测点的数据，它必须由数据框的形式输入，level用于输入置信度，默认是95%，返回Y期望值与区间估计．

例如在例1中，若标准混合肥料的施肥量是10英担，则预测马铃薯单位产量如下： > X=c(0,4,8,12);Y=c(8.34,8.89,9.16,9.50) > lm.sol<-lm(Y~1+X) > new<-data.frame(X=10) > lm.pred<-predict(lm.sol,new,interval="prediction",level=0.95);lm.pred 命令“new<-data.frame(X=10)”的目的是以数据框的形式输入新数据new,interval="prediction"的目的是要求输出期望值的区间估计．

软件输出是： fit lwr upr [1,] 9.3475 8.867725 9.827275 说明，当标准混合肥料的施肥量是10英担，则预测马铃薯单位产量期望值点估计是9.3475吨/英亩，期望值95%的区间估计是[8.867725，9.827275]吨/英亩．

线性回归的另一应用则是质量控制．例如，在例2 中，我们已经得到了生产某种产品的同类企业产量与成本之间的回归方程是 并极其显著地通过了假设检验．

假设通过调查得到另一个同类企业的产量7. 5千吨，容易得到成本期望值95%的区间估计是[123. 7875 172 假设通过调查得到另一个同类企业的产量7.5千吨，容易得到成本期望值95%的区间估计是[123.7875 172.2988] 万元，若实际成本数据小于123.7875，则说明该企业在成本控制上取得了重大突破，如果实际成本大于172.2988，则说明该企业在成本控制中出现重大问题．

一般而言，如果实际数据超出了区间估计的范围，则说明在质量控制上出现异常．

6-1-4 一个著名的实例 例3 十九世纪四五十年代，当时已经可以通过气压计较为精确地测量大气压，并通过大气压估计海拔高度，海拔越高，气压越低．但在当时运输精密的气压计相当困难，于是苏格兰物理学家James D. Forbes试图通过水的沸点来估计大气压进而估计海拔高度．该项研究将给旅行者提供通过测量水的沸点快速估计海拔高度的方法．Forbes在阿尔卑斯山及苏格兰山区收集数据．数据计量单位，沸点单位采用的是华氏，气压单位采用的是英寸汞柱．从他1857年的论文中选取了17对沸点与气压数据，如下表，试用回归分析

观察点序号 沸点（华氏） 气压（英寸汞柱） 1 194.5 20.79 2 194.3 3 197.9 22.40 4 198.4 22.67 5 199.4 23.15 6 199.9 23.35 7 200.9 23.89 8 201.1 23.99 9 201.4 24.02 10 201.3 24.01 11 203.6 25.14 12 204.6 26.57 13 209.5 28.49 14 208.6 27.76 15 210.7 29.04 16 211.9 29.88 17 212.2 30.06

下面是R操作过程： 首先输入数据 > x=c(194.5, 194.3, 197.9, 198.4, 199.4, 199.9, 200.9, 201.1, 201.4, 201.3, 203.6, 204.6, 209.5, 208.6, 210.7, 211.9, 212.2) > y=c(20.79, 20.79, 22.4, 22.67, 23.15, 23.35, 23.89, 23.99, 24.02, 24.01, 25.14, 26.57, 28.49, 27.76, 29.04, 29.88, 30.06)

作散点图，如图下图. > plot(x,y)

对x与y，线性回归分析如下： > lm.sol<-lm(y~1+x);summary(lm.sol) 软件输出如后

Call: lm(formula = y ~ 1 + x) Residuals: Min 1Q Median 3Q Max -0.25717 -0.11246 -0.05102 0.14283 0.64994 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -81.06373 2.05182 -39.51 <2e-16 *** x 0.52289 0.01011 51.74 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.2328 on 15 degrees of freedom Multiple R-Squared: 0.9944, Adjusted R-squared: 0.9941 F-statistic: 2677 on 1 and 15 DF, p-value: < 2.2e-16

结果显示回归方程是 并通过了回归系数的假设检验与回归方程的检验，残差标准差是0.2328．

在散点图中作出回归直线，命令如下： > abline(lm.sol1) 输出如图

观察可知，第12个数据点（204.6，26.57）残差较大，该样本点的数据采集可能存在问题，剔除这一点，重新作回归分析如下： > x=c(194.5, 194.3, 197.9, 198.4, 199.4, 199.9, 200.9, 201.1, 201.4, 201.3, 203.6, 209.5, 208.6, 210.7, 211.9, 212.2) > y=c(20.79, 20.79, 22.4, 22.67, 23.15, 23.35, 23.89, 23.99, 24.02, 24.01, 25.14, 28.49, 27.76, 29.04, 29.88, 30.06) > lm.sol<-lm(y~1+x);summary(lm.sol)

软件输出为： Call: lm(formula = y ~ 1 + x) Residuals: Min 1Q Median 3Q Max -0.21493 -0.09546 -0.01500 0.09049 0.27793 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -80.667294 1.419984 -56.81 <2e-16 *** x 0.520738 0.006997 74.42 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1608 on 14 degrees of freedom Multiple R-Squared: 0.9975, Adjusted R-squared: 0.9973 F-statistic: 5538 on 1 and 14 DF, p-value: < 2.2e-16

结果显示回归方程是 并通过了回归系数的假设检验与回归方程的检验，残差标准差是0.1608，比刚才的0.2328要小．

6-2 多元线性回归分析 6-2-1 多元线性回归模型及软件操作 上一小节介绍的直线回归分析，由于只有一个解释变量，描述了一个因变量Y与一个解释变量X之间的线性依存或伴随关系，因而又称之为一元线性回归分析．而事物之间的联系常常是多方面的，实际问题中影响因变量Y取值的解释变量往往不止一个，通常设为p个，记为X1，X2，…，Xp，也可记之为向量X=（X1 X2 … Xp）．由于一般无法借助于图形来确定模型类型，所以仅讨论一种最为简单但又应用广泛的模型，即多元线性回归模型．

多元线性回归模型的一般表达式是 式中，是因变量的估计值，X1，X2，…，Xp是解释变量，p是自变量的个数，β0是回归方程常数项，β1，β2，…，βp是偏回归系数，将其合计为向量β=(β0 β1 … βp)T．

β= (β0，β1，β2，…，βp)T=A-1·B 其中，xij是第i个变量的第j个观察值，计算矩阵A与矩阵B如下 A=XT·X，B=XT·Y 其中，XT是矩阵X的转置矩阵，也就是行列互换的矩阵. 最后计算出回归方程的常数项与偏回归系数向量如下： β= (β0，β1，β2，…，βp)T=A-1·B 上述结论也是根据最小二乘法推导而得，推导过程从略．

在R软件中，多元线性回归分析的命令函数与一元线性回归分析一样，都是lm( )，现举例如下．

例1 根据经验，在人的身高相等的情况下，血压的收缩压Y与体重X1(kg)、年龄X2(岁)有关．现集了13名成人男子的数据，试作线性回归分析． 序号 1 2 3 4 5 6 7 8 9 10 11 12 13 X1 76.0 91.5 85.5 82.5 79.0 80.5 74.5 85.0 76.5 82.0 95.0 92.5 X2 50 20 30 60 40 55 Y 120 141 124 126 117 125 123 132 155 147

R软件操作如下： > x1=c(76.0, 91.5, 85.5, 82.5, 79.0, 80.5, 74.5, 79.0, 85.0, 76.5, 82.0, 95.0, 92.5) > x2=c(50, 20, 20, 30, 30, 50, 60, 50, 40, 55, 40, 40, 20) > y=c(120, 141, 124, 126, 117, 125, 123, 125, 132, 123, 132, 155, 147) > lm.sol<-lm(y~1+x1+x2);summary(lm.sol)

软件输出如下： Call: lm(formula = y ~ 1 + x1 + x2) Residuals: Min 1Q Median 3Q Max -4.0404 -1.0183 0.4640 0.6908 4.3274 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -62.96336 16.99976 -3.704 0.004083 ** x1 2.13656 0.17534 12.185 2.53e-07 *** x2 0.40022 0.08321 4.810 0.000713 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.854 on 10 degrees of freedom Multiple R-Squared: 0.9461, Adjusted R-squared: 0.9354 F-statistic: 87.84 on 2 and 10 DF, p-value: 4.531e-07

结果表明，多元线性回归方程是 该方程通过了回归系数的假设检验与回归方程的检验．

6-2-2 多元线性预测 当多元线性回归方程经过假设检验是显著的，且每一个回归系数均具有统计意义时，可用此方程作线性预测，包括因变量的期望值预测与置信区间的预测．

例2 在例1中，经过回归分析，我们得到了多元线性回归方程，并且该方程通过了回归系数的假设检验与回归方程的检验，那么与一元线性回归分析过程中的预测一样，也可通过命令predict( )，求得X =（75，40）时收缩压的点估计值及期望值区间估计．R操作如下： > new<-data.frame(x1=75,x2=40) > lm.pred<-predict(lm.sol,new,interval="prediction",level=0.95);lm.pred

软件输出 fit lwr upr [1,] 113.2871 106.0412 120.5331 结果显示，因变量的估计值是113.2871，95%的期望值区间估计是 [106.0412 ，120.5331]