SAS系统和数据分析逐步回归分析

1、第三十三课逐步回归分析一、逐步回归分析在一个多元线性回归模型中，并不是所有的自变量都与因变量有显著关系，有时有些自 变量的作用可以忽略。这就产生了怎样从大量可能有关的自变量中挑选出对因变量有显著影 响的部分自变量的问题。在可能自变量的整个集合有 40到60个，甚至更多的自变量的情况下，使用“最优”子 集算法可能并不行得通。那么，逐步产生回归模型要含有的X变量子集的自动搜索方法，可能是有效的。逐步回归方法可能是应用最广泛的自动搜索方法。这是在求适度“好”的自变 量子集时，同所有可能回归的方法比较，为节省计算工作量而产生的。从本质上说，这种方 法在每一步增加或剔除一个 X变量时，产生一系列回归模型

2、。 增加或剔除一个 X变量的准则， 可以等价地用误差平方和缩减量、偏相关系数或F统计量来表示。无疑选择自变量要靠有关专业知识，但是作为起参谋作用的数学工具，往往是不容轻视 的。通常在多元线性模型中，我们首先从专业角度选择有关的为数众多的因子，然后用数学 方法从中选择适当的子集。本节介绍的逐步回归法就是人们在实际问题中常用的，并且行之 有效的方法。逐步回归的基本思想是，将变量一个一个引入，引入变量的条件是偏回归平方和经检验 是显著的，同时每引入一个新变量后，对已选入的变量要进行逐个检验，将不显著变量剔除， 这样保证最后所得的变量子集中的所有变量都是显著的。这样经若干步以后便得“最优”变 量子集。

3、逐步回归是这样一种方法，使用它时每一步只有一个单独的回归因子引进或从当前的回 归模型中剔除。Efroymoson （1966）编的程序中，有两个F水平，记作Fin和Fout,在每一步时， 只有一个回归因子，比如说人，如果剔除它可能引起 RSS的减少不超过残差均方MSE （即ESS/（N-k-1）的Fut倍，则将它剔除；这就是在当前的回归模型中，用来检验*=0的F比（RSS（x1,X2/ Xi，x-RSSfXjX?, Xi 4）/MSE是小于或等于 Fu。若剔除的变量需要选择，则就选择使RSS减少最少的那一个（或等价的选择 F比最小的）。 用这种方式如果没有变量被剔除，则开始引进一个回归因子，比

4、如Xj,如果引进它后使 RSS的增加，至少是残差均方的Fin倍，则将它引进。即若在当前模型加Xj项后，为了检验 =0的F比，F Fin时，则引进 ,其次，若引进的变量需要选择，则选择F比最大的。程序按照上面的步骤开始拟合，当没有回归因子能够弓I进模型时，该过程停止。二、变量选择的方法若在回归方程中增加自变量 Xi，称为“引入”变量 Xi，将已在回归方程中的自变量 Xj从 回归方程中删除，则称为“剔除”变量 Xj。无论引入变量或剔除变量，都要利用 F检验，将 显著的变量引入回归方程，而将不显著的从回归方程中剔除。记引入变量F检验的临界值为Fin （进），剔除变量F检验的临界值为Fut （出），一

5、般取FinFout，它的确定原则一般是对 k上海财经大学经济信息管理系IS/SHUFE* 个自变量的m个(m Fin，则首先将X1引入回归方程，不失一般性，设 Xi就是。接着考虑X1分别与X2,X3,.,Xk与因变量Y组成二元回归方程，对于这k 1个回归方程中 X2,Xk的回归系数进行F检验，计算F值，并选其最大的 F值Fj2,若Fj2 Fin,则接着就将 Xj引入回归方程，不失一般性，设 召就是X2。对已经引入回归方程的变量 X1和X2,如同前面的方法做下去，直至所有未被引入方程的变量的F值均小于Fin时为止。这时的回归方程就是最终选定的回归方程。显然，这种增加法有一定的缺点，主要是，它不能

6、反映后来变化的情况。因为对于某个 自变量，它可能开始是显著的，即将其引入到回归方程，但是，随着以后其他自变量的引入， 它也可能又变为不显著了，但是，并没有将其及时从回归方程中剔除掉。也就是增加变量法， 只考虑引入而不考虑剔除。2. 变量减少法与变量增加法相反，变量减少法是首先建立全部自变量X1,X2,.,Xk对因变量丫的回归方程，然后对k个回归系数进行 F检验，记求得的 F值为 FF；,F：，选其最小的记为Fi1 =min可，,Fk1,若有Fi1 Fout，则可以考虑将自变量人从回归方程中剔除掉，不妨设Xj就取为x1。再对X2,X3,.,Xk对因变量丫建立的回归方程重复上述过程，取最小的F值为

7、Fj2，若有Fj2 Fout，则将Xj也从回归方程中剔除掉。不妨设 Xj就是X2。重复前面的做法，直至在回归方 程中的自变量F检验值均大于Fout，即没有变量可剔除为止。这时的回归方程就是最终的回 归方程。这种减少法也有一个明显的缺点，就是一开始把全部变量都引入回归方程，这样计算量 比较大。若对一些不重要的变量，一开始就不引入，这样就可以减少一些计算。3. 变量增减法前面的两种方法各有其特点，若自变量X1,X2,.,Xk完全是独立的，则可结合这两种方法，但是，在实际的数据中，自变量X1,X2,.,Xk之间往往并不是独立的，而是有一定的相关性存在的，这就会使得随着回归方程中变量的增加和减少，某些

8、自变量对回归方程的贡献也会发 上海财经大学经济信息管理系IS/SHUFE生变化。因此一种很自然的想法是将前两种方法综合起来，也就是对每一个自变量，随着其 对回归方程贡献的变化，它随时可能被引入回归方程或被剔除出去，最终的回归模型是在回 归方程中的自变量均为显著，不在回归方程中的自变量均不显著。三、弓I入变量和剔除变量的依据如果在某一步时，已有 丨个变量被引入到回归方程中，不妨设为X1,X2/ ,X|，即已得回归方程：Y= 1X12X2 ，|X|(33.1)并且有平方和分解式：TSS 二 RSS ESS(33.2)显然，回归平方和 RSS及残差平方和ESS均与引入的变量相关。为了使其意义更清楚起

9、见，将其分别设为 RSS( Xi,X2,，XJ及ESS( Xi,X2,，X|)。下面我们来考虑，又有一个变量Xi (lik)被引入回归方程中，这时对于新的回归方程所对应的平方和分解式为：TSS= RSS( Xi,X2, ,X| , Xi) + ESS( XX?, ,X|, Xi)(33.3)当变量Xj引入后，回归平方和从 RSS( X1,X2/ ,X|)增加到RSS( X1,X2/ ,X| ,Xi)有：，而相应的残差平方和却从 ESS( X1,X2/ ,X|)降到ESS( X1,X2/ ,X|, X)并RS$X1,X2, ,X| , Xi)- RS$X1,X2, ,X|)(33.4)= ESS

10、X1,X2, ,X|)-ESSX1,X2, ,X| , Xi)记 Wi 二RS8Xi,X2，X|,XJ -RSXi,X2，X)它反映了由于引入 Xi 后,Xi对回归平方和的贡献，也等价于引入Xj后残差平方和所减少的量，称其为 Xi对因变量Y的方差贡献，故考虑检验统计量：Lw(Xi,X2，X| )(33.5)Fi 一 Fin，则可以考虑将Fi ESSXi,X2, ,X|,Xi / N -1 -1其中N为样本量，I是已引入回归方程的变量个数，这时若有 自变量Xi引入回归方程，否则不能引入。实际上大于Fin的变量开始时可能同时有几个，那么是否将它们都全部引入呢？实际编程上海财经大学经济信息管理系IS

11、/SHUFE序时并不是一起全部引入，而是选其最大的一个引入回归方程。关于剔除变量，如果已有丨个变量被引入回归方程， 不失一般性，设其为X1,X2 - ,X,，所对应的平方和分解公式为：TSS = RSS(Xi,X2, ,Xi, ,XJ ESS(Xi,X2, ,Xi, XJ(33.6)其中i =1,2，丨为了研究每个变量在回归方程中的作用，我们来考虑分别删掉Xi(i=1,2,.,l后相应的平方和分解公式为：TSS二RSXi,X2, Xi,Xi,X|) ESS(Xi,X2, ,Xi，Xi 1 ,X|)(33.7)这时，回归平方和从 RSSXi,X2, ,Xi, ,X|)降为 RSSXi,X2, X

12、i，Xi 1 ,X|),同时残差也发生相应的变化。残差平方和从ESS(X1,X2，Xi，Xl)增加到ESS1,X2/ ,XiJ,Xi,X|)，xi对回归平方和的贡献，也等价于删除xi后残差平方和所增加的量，同理可表示为：W 二 RSSXi,X2, ,Xi, X|)-RSSXi,X2, ,Xi，Xi 1 X|)(3= ESS(Xi，X2， ,XXiX|)-es&Xi,X2, ,Xi， X|)同理，我们来构造检验统计量：(33.9)W Xi,X2, ,Xi, ,X|ESSXi,X2, ,Xi, X| / N -| -1显然，这时Fi越小，则说明Xj在回归方程中起的作用(对回归方程的贡献)越小，也就

13、是若有Fi F84.01 0.0001F ProbFINTERCEP82.421772683.855303783443.36654076457.050.0001RUNTIME-3.310555360.36119485632.9000998584.010.00011, 1Bounds on condition number:Step 2 Variable AGE Entered R-square = 0.76424693 C(p) = 12.22493455DFSum of SquaresMean SquareF ProbFRegression 2650.66573237325.33286618

14、45.380.0001Error28200.715812477.16842187Total30851.38154484ParameterStandardType IIVariableEstimateErrorSum of SquaresF ProbFINTERCEP 88.462287495.372638851943.41070877271.110.0001AGE -0.150365670.0955146817.765632522.480.1267RUNTIME -3.203950560.35877488571.6775057979.750.0001Bounds on condition nu

15、mber: 1.036941,4.147763Step 3 Variable RUNPULSE Entered R-square = 0.81109446 C(p) = 6.82780371DFSum of SquaresMean SquareF ProbFRegression 3690.55085627230.1836187638.640.0001Error27160.830688575.95669217Total30851.38154484ParameterStandardType IIVariableEstimateErrorSum of SquaresF ProbFINTERCEP11

16、1.71806443 10.23508836709.69013814 119.14 0.0001AGE -0.256398260.0962289242.288674387.100.0129RUNPULSE-0.130908700.0505901139.885123906.70 0.0154RUNTIME -2.825378670.35828041370.4352860762.190.0001Bounds on condition number: 1.354763,11.59745Step 4 Variable MAXPULSE Entered R-square = 0.83681815 C(p

17、) = 4.76608569DF Sum of Squares Mean Square F ProbFRegression 4712.45152692178.1128817333.330.0001在输出结果报告中，提供了进入回归变量逐次改变后回归方差分析和拟合的信息。在报告的最后部分，列出了用逐步回归法挑选自变量过程，四个自变量按run time , age, ru npulse，maxpulse先后次序进入回归模型。所有进入回归的变量在0.15的水平下是显著的，未进入回归的候选变量在 0.15的水平下是不显著的。同时还概要地提供了每个回归模型变化时的 R2增加值、R2值、CP值、相应的F统计

18、量、p值。在逐步回归的每步细节中，还列出了条 件指数的最小值最大值，以及每一个回归变量的类型2平方和。age变量进入模型后，R2值的增加值（Partial R2,称为偏 R2或部分 R2）计算为（650.6657- 632.9001）/ 851.3815= 0.020867。 如果按CP值选择最优子集，随着进入回归模型中的自变量个数P从2增加到5个（包括截上海财经大学经济信息管理系IS/SHUFE距)，相应 CP 值从大到小为 13.51976469、12.22493455、6.82780371 和 4.76608569，按照 Mallows提出的回归模型最优自变量个数的选择准则，CP=4.7

19、6608569是最接近自变量个数P=5 的模型。CP 的计算公式见式(33.11),当 P=5 时，CP=138.93001792/5.39197 (31 - 2 X 5)= 4.76608569。因此，用逐步回归方法及CP值确认的拟合回归模型为：oxyger= 98.14788797 0.佃773470age + 0.27051297maxpulse 0.34810795runpulse 2.76757879runtime条件指数(condition number )为最大特征值和每个特征值之比的平方根。我们看到，当 模型进入第四个自变量 maxpulse时，最大的条件指数从较小的11.59

20、745变成了较大的76.85135，说明存在一定程度的共线性，根据前面例33.2的分析，我们诊断这个共线性方程可能为 run pulse maxpulse=0。在向前、向后或逐步回归的变量选择过程中，都有一个判断是否可进入或剔除的显著水 平，在程序中是分别由 model语句的选项slentry =和slstay =设定的，缺省的情况见表33.2。表33.2缺省的入选和剔除显著水平forwardbackwardstepwizesle ntry0.500.15slstay0.100.15F面我们提供全部可能回归的程序，并且以R2值由大到小的排序输出。proc reg data= fitness ;

21、model oxygen = age weight rstpulse maxpulse runpulse runtime/selection= rsquare b ;run ;在上述程序中，model语句的选项selection= rsquare,表示请求R2值最大法，选项b是表 示要输出每种回归的回归系数。程序运行后，得到如表33.3所示的结果。表33.3用R2排序全部可能的变量数的逐步回归分析结果N = 31Regression Models for Dependent Variable: OXYGENParameterNumber in R-square EstimatesModeInt

22、ercept AGE WEIGHT RSTPULSE MAXPULSE RUNPULSE RUNTIME10.7433801082.4218-3.310610.1583834482.4582-0.2068 .10.1199967059.3325-0.2225.10.0927765362.2206-0.311410.0560459271.2907-0.1376.10.0264884955.4379-0.1041.20.7642469388.4623-0.1504.-3.204020.7614238193.0888-0.0735 -3.1402程序的输出包括所有只含一个变量的6种回归，含2个变量的

23、15种回归总共有 63种不同形式的回归模型。例如，含2个自变量按R2第二个大值选择回归模型为，R2=0.76142381，拟合的回归模型为：oxyger= 93.0888 0.0735runpulse 3.1402runtime若对每种变量个数，只要保留R2最大的两种情况，可在model语句中加入选项 best=2,即提交proc reg data= fitness ;model oxygen = age weight rstpulse maxpulse runpulse runtime/selection= rsquare b best=2 ;run ;这33.4所示的结果。表33.4只保留

24、R2最大两种情况的逐步回归分析结果N = 31 Regression Models for Dependent Variable: OXYGENParameterNumber in R-square EstimatesModelIntercept AGE WEIGHT RSTPULSE MAXPULSE RUNPULSE RUNTIME1 0.7433801082.4218-3.31061 0.1583834482.4582.-0.2068.2 0.7642469388.4623 -0.1504.-3.20402 0.7614238193.0888.-0.0735 -3.14023 0.811

25、09446111.7 -0.2564.-0.1309 -2.82543 0.8099884480.9008.0.3542 -0.3751 -2.97024 0.8368181598.1479 -0.1977.0.2705 -0.3481 -2.7676通过上面的逐步回归分析，我们已经得到回归模型的自变量个数确定时的最优子集或次优子集，但问题是我们到底应该选择几个自变量的回归模型呢？如表33.4中的3个自变量、4个自变量、5个自变量、6个自变量的回归模型中该选哪一个模型呢？ 一种最简便确定回归 模型的自变量个数的方法是Mallows的Cp方法。确定好模型的自变量个数后，根据表 33.4就很容易确

26、定在这个固定自变量数下，最优的自变量组合和相应的参数值估计。以下的程序 是对所有可能的回归按Cp由小到大进行排序并保留其前5种，并绘制Cp图。goptions reset=global gunit=pct cback=white borderhtitle=6 htext=3 ftext=swissb colors=（back）;title Cp plot with Reference Lines;proc reg data= fitness ;model oxygen = age weight rstpulse maxpulse runpulse runtime/selection=cp ad

27、jrsq best=5 ;plot cp. Model语句中的selection=cp选项请求计算 Mallows的Cp统计量。选项 adjrsq表示要显 示每种回归模型的统计量 Adj- R2。选项best=5表示保留Cp值最小的前5种。plot语句中的cp. np.表达式（注意统计量关键字母后的小圆点）表示 Y轴为Cp值，X轴为P值（P值包括截 np. /chocking=red cmallows=bluevaxis=0 to 15 by 2haxis=0 to 8 by 1;run ;距项)。plot语句的选项chocking=red，表示画Cp=2P Pfuii红色参考虚线，其中P是子

28、模型中含截距的参数个数，Pfuii是全模型中不含截距的参数个数。Hoching( 1976)建议选择满足Cpw 2P PfU|且Cpw P的模型。plot语句的选项 cmallows=blue，表示画Cp=p蓝色参考实线， 其中P是子模型中含截距的参数个数。Mallows (1973)建议考虑所有满足 Cp较小且接近P的模型。这一程序的输出结果如表33.5和图33. 1所示。N = 31C(p) R-square AdjustedIn R-squareRegression Models for Dependent Variable: OXYGENVariables in Model4.7660

29、90.8368181540.811713255.000210.8480018150.817602186.752590.8369035950.804284316.827800.8110944630.79010496AGE MAXPULSE RUNPULSE RUNTIMEAGE WEIGHT MAXPULSE RUNPULSE RUNTIMEAGE RSTPULSE MAXPULSE RUNPULSE RUNTIMEAGE RUNPULSE RUNTIME表33.5按Cp由小到大进行排序并保留其前5种逐步回归分析结果Cp plot wHh Reference LinesPlot i 1 CP*P

30、CP = 2P (P fef M modal) + 1CP = P图33.1 带有Mallows和Hocking参考线的 Cp散点图第三十四课 从输出结果可看出，以Mallows的建议为标准， age, maxpulse, run pulse 和run time四个变量进入回归模型时 Cp最 小(4.76609),且与P=4+1=5最接近，因为5-4.766092X 5 6=4。而 Cp=5.00021 的 模型满足要求，因为 5.000212 X 6 6=6 且 5.000210.83681815)。不同的 标准提供不同的选择结果，这是常有的情况。非线性回归分析现实世界中严格的线性模型并不多

31、见，它们或多或少都带有某种程度的近似；在不少情 况下，非线性模型可能更加符合实际。由于人们在传统上常把“非线性”视为畏途，非线性 回归的应用在国内还不够普及。事实上，在计算机与统计软件十分发达的令天，非线性回归 的基本统计分析已经与线性回归一样切实可行。在常见的软件包中(诸如SAS、SPSS等等)，人们已经可以像线性回归一样，方便的对非线性回归进行统计分析。因此，在国内回归分析 方法的应用中，已经到了“更上一层楼” ，线性回归与非线性回归同时并重的时候。对变量间非线性相关问题的曲线拟合，处理的方法主要有：首先决定非线性模型的函数类型，对于其中可线性化问题则通过变量变换将其线性化，从而归结为前面

32、的多元线性回归问题来解决。若实际问题的曲线类型不易确定时，由于任意曲线皆可由多项式来逼近，故常可 用多项式回归来拟合曲线。若变量间非线性关系式已知(多数未知)，且难以用变量变换法将其线性化，则进行数值迭代的非线性回归分析。六、可变换成线性的非线性回归在实际问题中一些非线性回归模型可通过变量变换的方法化为线性回归问题。例如，对 非线性回归模型2yt -匚 o 亠二ai cosixt - bi sin ixt i亠：t(34.1)即可作变换：Xit = cosxt, X2t 二 sin xt, X3t = cos2xt ,X4t = sin2xt将其化为多元线性回归模型。一般地，若非线性模型的表达

33、式为：yt=bb bgiXtb2g2 Xbmgm xt(34.2)则可作变量变换：X1tgiXt, X2t g 2Xt, ,Xmt = gmXt(34.3)将其化为线性回归模型的表达式，从而用前面线性模型的方法来解决，其中式(34.3)中的Xt也可为自变量构成的向量。这种变量变换法也适用于因变量和待定参数bi。如：yt =aexpbxit +b2X?t + b3(XitX/ 1 卩(34.4)时上式两边取对数得：In yt = In a bt b?X2t b3 XitX2t -1(34.5)现作变换：yt =lnyt,b =ln a,X3t = XitX2t -1(34.6)则可得线性表达式：

34、y；二 bo biXit b2X2t b3X3t(34.7)利用前面方法确定了 ！?，i=0,1,2,3，并由0?=exp(?)得到？的值。变量变换的线性化方法可推广到下列形式的非线性模型：h(yj = c0(b。) G(bjgi 人Cm(bm)gm 人(34.8)其中x= (X1,X2,xp),而h (yt)、G ( b )、gi (xj则分别化为新的因变量、线性回归 参数和自变量，即可归结为线性回归模型来解。表34.1给出了一些常见的可线性化的非线性模型。表34.1典型的函数及线性化方法函数名称函数表达式线性化方法双曲线函数1b=a 十一 yx11v =u =yx幂函数by =axv =

35、l n y u = l n x指数函数bxy =aev = l n y u = xb/ xy =ae1v = l n y u =x对数函数y = a + bl n xv = yu = l n xS型函数1y =a + be1丄v =一u = ey当曲线的函数类型未确定时，我们常采用上述非线性模型作为其拟合曲线，即将自变量 的各种初等函数的组合作为新自变量，用逐步回归法（或正交筛选法等）对新变量进行筛选, 以确定一个项数不多的线性函数表达式。该方法对表达式形式没限制且精度要求不高的问题 颇为有效。七、多项式回归分析在式（34.2）中，若取gi x = xi，则为多项式回归模型。由数学分析知识可知

36、，一般函数都可用多项式来逼近，故多项式回归分析可用来处理相当广泛的非线性问题。 对观测数据（xt, yt） （t= 1，N）,多项式回归模型为：yt二bRxtbzx2亠亠bmXtm；t，t=1,2，Nyily2Jn1X12X12X2X2JXn2XnmX-imX2blbiPm?N则模型可表示为:当X列满秩时，由前面的讨论知，其最小二乘估计为:齐 XX XY由此即可求得其多项式回归方程。但由于xX J的计算既复杂又不稳定，故我们一般采用正交多项式法来进行多项式回归。八、不可变换成线性的非线性回归分析假设因变量y与自变量(xi, X2,，xp )之间满足非线性模型：y = F XiK, ,Xp； 1

37、；(34.9)其中，卩：丄1,辽，m为未知参数，F为已知表达式，；为误差项。现将观察数据：yt,Xit,X2t, ,Xpt ,t=1,2，N代人式(34.9)得非线性回归模型：% 二 F Xit,X2t, ,Xpt；t,t=l,2, ,N常记为：Y = F C ) E其中，丫二yi, y2/ , yN为y的观察向量， r,，：m为非线性回归系数，E=；1,2，；N为观察误差向量，F为未知参数1的函数向量。非线性回归分析就是利用最小二乘准则来估计回归系数1 ,即求?使得残差平方和：11Q EEY-F： Y-F：22在-_ ?处达到最小。非线性回归分析一般用数值迭代法来进行，其共同特点是：由选定1

38、的初值1 0出发，通过逐步迭代： = -0 t . ：(34.10)即选择适当的步长t ( 0 )及确定搜索方向向量 心=(鸟，心2,，Am),使得：Q :：: Q - 0(34.11)再由取代1 0，重复上述迭代过程，直至Q( -)可认为达到最小值为止，即可将所得的1作为其最小二乘估计？，从而得到非线性回归方程? = F X1,X2, ,Xp； ?1.下降方向和步长的选择44*首先考察Q 一： =EE = 1 丫一 F iii iY 一 Fl的梯度向量（即导数）22Y-F ：G Y F 1 3丿其中,G 1,，m为F的梯度矩阵。（34.11）。现考虑一元函数为使1 迭代收敛到?，其迭代公式应

39、满足下降性质t二Q，它从一：0出发以 厶为方向的射线上取值。由复合求导公式得：d 八 t 宀=乎.Y F 1 G :可以证明，当 d0时，在以 厶为方向向量的射线上可以找到r- : T0 -1，使得Q 1 : Q -0。我们将满足 d0的厶称为下降方向，Bard于1974年给出了厶为下降方向的 充要条件为：:二 PG Y - F 其中，P为对称正定阵，由此我们可得下降算法的迭代公式为：- = -0 tPG Y -F.（34.12）其中，P为任意正定阵，G为F的梯度，t为满足Q 1 :：: Q 0的正实数，即步长。如何计算厶以便修改参数向量1有五种常用的非线性回归迭代方法：高斯一牛顿法（Gauss-Newton）、最速下降法（梯度法，Gradient）、牛顿法（Newton ）、麦夸特法（Marquardt ）、 正割法（DUD ）。以下我们介绍高斯牛顿法。2. 咼斯一牛顿法首先选取1的一