MATLAB实现回归分析：从书籍实用数学方法中学习

例1 某公司收集了30个销售周期本公司产品价格与其他品牌价格差，投入广告费用和销售量数据见表6.5.1。试建立产品销售量与价格、广告投入之间的回归模型，预测在不同价格和广告费用下的牙膏销售量。

表6.5.1　产品销售情况

解 设y为公司产品销售量，x1为其他厂家与本公司价格差，x2为公司广告费用。画出分别以x1和x2为自变量的散点图，如图6.5.1和6.5.2所示。

图6.5.1　价格差与销售量散点图

图6.5.2　广告费用与销售量散点图

由散点图中x1,x2与y的关系，可设y=β0+β1 x1+β2 x2+β3 x22+ε，其中β0,β1,β2,β3分别为回归系数，ε为随机误差（均值为零的正态分布随机变量）。下面通过MATLAB统计工具箱对模型求解。直接利用MATLAB统计工具箱中的命令regress求解，使用格式：

输入：y为输入y的数据；x为对应于回归系数ai的数据矩阵；slpha为置信区间α。

输出：b为ai的估计值；bint为b的置信区间；r为残差向量；rint为r的置信区间；stats为回归模型的检验统计量，共有4个值，第1个是回归方程的决定系数R2（R是相关系数），第2个是F统计量值，第3个是与F统计量值对应的概率值P，第4个是剩余方差s2。

结果见表6.5.2。

表6.5.2　计算结果

从表6.5.2可以看出，y的90.54%可由模型确定，p值远小于α=0.05，F值远超过F检验的临界值，模型从整体上看成立，得

例如控制价格差x1=0.2元，投入广告费x2=6.5百万元，代入上式得ˆ8.293 3 y=，销售量预测区间为[7.823 0，8.763 6]（置信度95%）。

模型改进：为随机误差。

在MATLAB命令窗口中输入以下命令：

运行后得到表6.5.3。

表6.5.3　改进模型的计算结果

控制价格差为x1=0.2元，投入广告费x2=6.5百万元，代入上式得=8.327 2，销售量预测区间为[7.895 3，8.759 2]（置信度95%）。yˆ值略有增加，预测区间长度更短。

例2 某公司人事部门为研究员工薪金与他们的资历、管理责任、教育程度等之间的关系，调查了46名员工的档案资料（见表6.5.4）。其中，资历指从事专业工作的年数，管理一列中1表示管理人员，0表示非管理人员，教育一列中1表示中学程度，2表示大学程度，3表示更高程度（研究生）。试由此建立模型研究薪金与资历、管理责任、教育程度之间的关系。

表6.5.4　员工档案资料

续表

解 假设资历每加一年，薪金的增长是常数；管理、教育、资历之间无交互作用。

y表示薪金，1x是资历（年），x2=1表示管理人员，x2=0表示非管理人员，

（1）建立线性回归模型。

薪金y与资历1x、管理责任x2、教育程度3x，x4之间的多元线性回归方程为

其中：a0 ,a1 ,a2 ,a3 ,a4是待估计的回归系数；ε 是随机误差。

（2）模型求解。

直接利用MATLAB统计工具箱中的命令regress求解，在MATLAB命令窗口中输入以下命令：

输入：y为输入y的数据；x为对应于回归系数ai的数据矩阵；slpha为置信水平（0.05）。

输出：b为β的估计值；bint为b的置信区间；r为残差向量；rint为r的置信区间；stats为检验统计量R2，F，p，s2。

运行MATLAB对模型进行求解，可以得到回归系数及其置信区间（置信水平α=0.05），检验统计量R2，F，P，s2的结果，见表6.5.5。第1个是回归方程的决定系数（R是相关系数），第2个是统计量值，第3个是与F统计量值对应的概率值，第4个是剩余方差。

表6.5.5　例2的计算结果

（3）结果分析。

从表6.5.5知2 0.957 R=，即因变量（薪金）的95.7%可由模型确定，F值远远超过F检验的临界值，P远小于α，因而模型从整体看是可用的。

模型中1x的系数为546，说明资历增加1年薪金增长546；2x的系数为6883，说明管理人员薪金多6 883；3x的系数为-2 994，说明中学程度薪金比更高的少2 994；x4的系数为148，说明大学程度薪金比更高的多148，但是注意到a4置信区间包含零点，说明这个系数的解释不可靠。

需要注意，以上解释是就平均值来说，而且一个因素改变引起的因变量的变化量，是在其他因素不变的条件下成立的。

a4的置信区间包含零点，说明基本模型存在缺点。为了寻找改进的方向，常用残差ε（薪金的实际值y与用模型估计的薪金之差）分析方法，将影响因素分成资历与管理—教育组合两类，管理—教育组合见表6.5.6。

表6.5.6　管理—教育组合

为了对残差进行分析，图6.5.3给出了ε与资历1x的关系，图6.5.4给出了ε与管理、教育组合间的关系。

(www.xing528.com)

图6.5.3　ε与1x的关系

图6.5.4　ε与管理—教育组合间的关系

从图6.5.3看，残差大概分成三个水平，这是由于6种管理—教育组合混合在一起，在模型中未被正确反映的结果；从图6.5.4看，对于前4个管理—教育组合，残差或者全为正，或者全为负，也表明管理—教育组合在模型中处理不当。

在模型中国管理责任和教育程度是分别起作用的，事实上，二者可能起着交互作用，如大学程度的管理人员的薪金会比二者分别的薪金制和高一点。因此，考虑在基本模型中增加2x与3x，4x的交互项，建立新的回归模型。

模型改进：

增加管理2x与教育3x，4x的交互项，令

利用MATLAB的统计工具得到的结果见表6.5.7。

由表6.5.7可知，改进模型的2R和F值均优于原模型，并且所有回归系数的置信区间都不含零点，表明改进模型是完全可用的。与原模型类似，做改进模型的两个残差分析图（图6.5.5、图6.5.6），可以看出，已经消除了图6.5.3、图6.5.4中的不正常现象，这也说明了改进模型的适用性。

表6.5.7　改进模型的计算结果

图6.5.5　改进模型（2）ε与1x的关系

图6.5.6　改进模型ε与x2-x3 ,x4组合的关系

从图6.5.5、图6.5.6发现一个异常点：具有10年资历、大学程度的管理人员（从表6.5.4可以查出是33号），他的实际薪金明显低于模型的估计值，也明显低于与他有类似经历的其他人的薪金。这可能是由我们未知的原因造成的，为了使个别的数据不致影响整个模型，剔除这个异常数据，对改进模型重新估计回归系数，得到的结果见表6.5.8，残差分析图见图6.5.7、图6.5.8，可以看出，剔除异常数据后的结果又有所改善。

表6.5.8　改进模型剔除异常数据后的计算结果

续表

图6.5.7　改进模型剔除异常数据后ε与1x的关系

图6.5.8　改进模型ε与x2-x3 ,x4组合的关系

对于改进回归模型，用剔除模型异常数据（33号）后估计出的系数，得到的结果是满意的。作为这个模型的应用之一，不妨用它来“制订”6种管理—教育组合人员的“基础”薪金（即资历为零的薪金，当然，这也是平均意义上的薪金）。利用改进模型和表6.5.8容易得到表6.5.9。

表6.5.9　6种管理—教育组合人员的“基础”薪金

从表6.5.9可以看出，大学程度的管理人员的薪金比研究生程度的管理人员的薪金高，大学程度的非管理人员的薪金比研究生程度的管理人员的薪金高，而大学程度的非管理人员的薪金比研究生的非管理人员的薪金略低。当然，这是根据这家公司实际数据建立的模型得到的结果，并不具普遍性。

例3 为了考评教师的教学质量，教学研究部门对学生进行问卷调查，得到12位教师的15门课程（其中3位教师有2门课）各项评分的平均值，按下面7项内容打分，分值从低到高为1~5分，见表6.5.10。

其中：X1表示内容组织的合理性；X2表示问题展开的逻辑性；X3表示回答学生的有效性；X4表示课下交流的有助性；X5表示教材的帮助性；X6表示考试的公正性；Y表示总体评价。

建立Y与X1~X6间简单、有效的模型，给教师提出建议。

表6.5.10　课程评分的平均值

问题分析：

从X1~X6中挑选出对Y影响显著的变量建立回归模型。将所有对Y影响显著的X选入模型，而影响不显著的X不选入模型，使模型中自变量个数尽可能少，为此利用逐步回归的方法。

逐步回归的基本思路是，确定一个包含若干X的初始集合S0，从S0外的X中引入一个对Y影响最大的，S0 →S1，对S1中的X进行检验，移出一个影响最小的，S1 →S2，继续进行，直到不能引入和移出为止，引入和移出都以给定的显著性水平为标准。

MATLAB统计工具箱中的逐步回归命令是stepwise（逐步回归），通常用法是

其中：x是自变量集合的n×k数据矩阵（n是数据容量，k是变量数目）；y是因变量数据向量（n维）；inmodel是初始模型S0中包括的自变量集合的指标（即矩阵X的列序数），缺省时为无自变量；penter是引入变量的显著性水平（缺省时为0.05）；premove为移出变量的显著性水平（缺省时为0.10），premove值不能小于penter值。

回归模型的建立：

将表6.5.10中的数据利用stepwise（x, y）命令得到stepwise regression的初始界面窗口，如图6.5.9所示。

图6.5.9　stepwise regression初始界面

初始界面的左上方用水平线给出了所有6个变量的回归系数的估计及误差界，其中在计算机屏幕上彩色的水平线表示置信度为90%的置信区间，灰色的为其95%的置信区间。水平线若为红色，表示该水平线所对应的变量没有被选入模型中。初始界面的上方中间部分的表格显示的是一旦某个变量被选入模型中，该变量的回归系数的估计值、检验的t统计量值以及p值。一般来说，在每一步中选入的是具有最小p值或最大统计量值的项目，在图6.5.9中即自变量X2，一般只需按Next step按钮进行下一步操作，程序就会自动选择所要引入或移出的自变量，并在界面的右上方给出相应的结果。也可以手动操作，用鼠标点击其中的一行，改变其状态，即目前不在模型中的一个变量（红色的行）被引入（变蓝），目前在模型中的一个变量（蓝色的行）被移出（变红），直到界面提示“Move no terms”为止。通常我们可以直接按All steps按钮来完成整个模型的逐步回归过程。如在图6.5.9中按All steps按钮，可得到逐步回归的最终结果，如图6.5.10所示。

图6.5.10中，界面的上方部分同图6.5.9。界面的中间部分的表格给出了该回归模型的所有计算结果，包括Intercept（截距，即回归常数）、决定系数R2，检验的F值，RMSE（剩余标准差）、调整的决定系数R2以及检验的p值，这里调整的决定系数R2定义为其中k为选入模型的自变量个数。

界面的下方部分Model History给出了逐步回归中每步所对应的模型的剩余标准差的点图（蓝色的点），将鼠标移至某步对应的蓝色点，会显示此步对应的模型中所包含的自变量，点击该蓝色的点，可追逐当前模型所对应的界面。

图6.5.10　stepwise regression最终界面

从图6.5.9可以看出，最终入选的自变量只有X1，X3（显示为蓝色），分析可知，随着逐步回归的进行，每步模型所对应剩余标准差（RMSE）逐渐减少到0.1，虽然在个别步骤中略有下降，但模型的F值却大大提高。这些都表明，仅含X1，X3的模型是合适的，由X1，X3的回归系数和回归常数，利用逐步回归最终得到的模型为

模型解释：

最终X1（内容组织的合理性）与X3（回答学生的有效性）列入教师考评的重点。得到的最终模型表明，X1的分值每增加1分，对教师的总体评价增加约0.5分；X1的分值每增加1分，对教师的总体评价增加约0.77分。

初步看来如果影响因变量的因素较多，并得到了大量的数据，为了建立一个有效的、便于应用的模型，可以利用逐步回归只选择那些影响显著的变量。如果认为原有变量的平方项、交叉项等也对因变量有显著影响，可以将这些项作为新的自变量加入候选行列，再用逐步回归处理。