1、1 / 12第三十三课 逐步回归分析一、 逐步回归分析在一个多元线性回归模型中，并不是所有的自变量都与因变量有显著关系，有时有些自变量的作用可以忽略。这就产生了怎样从大量可能有关的自变量中挑选出对因变量有显著影响的部分自变量的问题。在可能自变量的整个集合有 40 到 60 个，甚至更多的自变量的那些情况下，使用“最优”子集算法可能并不行得通。那么，逐步产生回归模型要含有的 X 变量子集的自动搜索方法，可能是有效的。逐步回归方法可能是应用最广泛的自动搜索方法。这是在求适度“好”的自变量子集时，同所有可能回归的方法比较，为节省计算工作量而产生的。本质上说，这种方法在每一步增加或剔除一个 X 变量时
2、，产生一系列回归模型。增加或剔除一个 X 变量的准则，可以等价地用误差平方和缩减量、偏相关系数或 F 统计量来表示。无疑选择自变量要靠有关专业知识，但是作为起参谋作用的数学工具，往往是不容轻视的。通常在多元线性模型中，我们首先从有关专业角度选择有关的为数众多的因子，然后用数学方法从中选择适当的子集。本节介绍的逐步回归法就是人们在实际问题中常用的，并且行之有效的方法。逐步回归的基本思想是，将变量一个一个引入，引入变量的条件是偏回归平方和经检验是显著的，同时每引入一个新变量后，对已选入的变量要进行逐个检验，将不显著变量剔除，这样保证最后所得的变量子集中的所有变量都是显著的。这样经若干步以后便得“最
3、优”变量子集。逐步回归是这样一种方法，使用它时每一步只有一个单独的回归因子引进或从当前的回归模型中剔除。Efroymoson (1966)编的程序中，有两个 F 水平，记作 Fin 和 Fout，在每一步时，只有一个回归因子，比如说 Xi，如果剔除它可能引起 RSS 的减少不超过残差均方MSE（即 ESS/(N-k-1)）的 Fout 倍，则将它剔除；这就是在当前的回归模型中，用来检验 i=0 的 F 比= 是小于或等于 Fout。MSExRSxxRSiii /),(),( 12112 若剔除的变量需要选择，则就选择使 RSS 减少最少的那一个（或等价的选择 F 比最小的） 。用这种方式如果没
4、有变量被剔除，则开始引进一个回归因子，比如 Xj，如果引进它后使 RSS 的增加，至少是残差均方的 Fin 倍，则将它引进。即若在当前模型加 Xj 项后，为了检验 j =0 的 F 比，F Fin 时，则引进 Xj，其次，若引进的变量需要选择，则选择 F 比最大的。程序按照上面的步骤开始拟合，当没有回归因子能够引进模型时，该过程停止。二、 变量选择的方法若在回归方程中增加自变量 Xi，称为“引入”变量 Xi，将已在回归方程中的自变量 Xj从回归方程中删除，则称为“剔除”变量 Xj。无论引入变量或剔除变量，都要利用 F 检验，将显著的变量引入回归方程，而将不显著的从回归方程中剔除。记引入变量 F
5、 检验的临界值为 Fin（进） ，剔除变量 F 检验的临界值为 Fout（出） ，一般取 Fin Fout，它的确定原则一般是对 k 个自变量的 m 个(m k) ，则对显著性水平 df1=1，df 2= 的 F 分布表的值，1mN记为 F*，则取 Fin=Fout= F*。一般来说也可以直接取 Fin=Fout=2.0 或 2.5。当然，为了回归方程中还能够多进入一些自变量，甚至也可以取为 1.0 或 1.5。2 / 121. 变量增加法首先对全部 k 个自变量，分别对因变量 Y 建立一元回归方程，并分别计算这 k 个一元回归方程的 k 个回归系数 F 检验值，记为 ，选其最大的记为 1iF
6、 = max121,kF121,kF,若有 1i F in，则首先将 X1 引入回归方程，不失一般性，设 Xi 就是 X1。接着考虑 X1 分别与 X2,X3,.,Xk 与因变量 Y 二元回归方程，对于这 k 1 个回归方程中X2,.,Xk 的回归系数进行 F 检验，计算得的 F 值，并选其最大的 F 值 2j,若 jF in,则接着就将 Xj 引入回归方程，不失一般性，设 Xj 就是 X2。对已经引入回归方程的变量 X1 和 X2，如同前面的方法做下去，直至所有末被引入方程的变量的 F 值均小于 Fin 时为止。这时的回归方程就是最终选定的回归方程。显然，这种增加法有一定的缺点，主要是，它不
7、能反映后来变化的情况。因为对于某个自变量，它可能开始是显著的，即将其引入到回归方程，但是，随着以后其他自变量的引入，它也可能又变为不显著的了，但是，也并没有将其及时从回归方程中剔除掉。也就是增加变量法，只考虑引入而不考虑剔除。2. 变量减少法与变量增加法相反，变量减少法是首先建立全部自变量 X1,X2,.,Xk 对因变变量 Y 的回归方程，然后对 k 个回归系数进行 F 检验，记求得的 F 值为 1F ，选其最小的记为 1iF=min 121,k ,若有 1iF out，则可以考虑将自变量 Xi 从回归方程中剔除掉，不妨设 Xi 就取为 X1。再对 X2,X3,.,Xk 对因变量 Y 建立的回
8、归方程中重复上述过程，取最小的 F 值为 2j，若有 jF out，则将 Xj 也从回归方程中剔除掉。不妨设 Xj 就是 X2。重复前面的做法，直至在回归方程中的自变量 F 检验值均大于 Fout，即没有变量可剔除为止。这时的回归方程就是最终的回归方程。这种减少法也有一个明显的缺点，就是一开始把全部变量都引入回归方程，这样计算量比较大。若对一些不重要的变量，一开始就不引入，这样就可以减少一些计算。3. 变量增减法前面的二种方法各有其特点，若自变量 X1,X2,.,Xk 完全是独立的，则可结合这二种方法，但是，在实际的数据中，自变量 X1,X2,.,Xk 之间往往并不是独立的，而是有一定的相关性
9、存在的，这就会使得随着回归方程中变量的增加和减少，某些自变量对回归方程的贡献也会发生变化。因此一种很自然的想法是将前二种方法综合起来，也就是对每一个自变量，随着其对回归方程贡献的变化，它随时可能被引入回归方程或被剔除出去，最终的回归模型是在回归方程中的自变量均为显著的，不在回归方程中的自变量均不显著。三、 引入变量和剔除变量的依据如果在某一步时，已有 个变量被引入到回归方程中，不妨设为 ，即已l lX,213 / 12得回归方程 lXXY210(33.1)并且有平方和分解式 ESRTS(33.2)显然，回归平方和 及残差平方和 均与引入的变量相关。为了使其意义更清楚起见，将其分别设为 RSS（
10、 ）及 ESS（ ） 。下面我们来考虑，l,21 l,21又有一个变量 （l;weight 变量 ;by 变量 ;run ;stepwise 至少需要一个 model 语句。by 语句和 weight 语句可以放在任何地方。1) model 语句的。stepwise 中可以有任意多个 model 语句。model 语句中的选项如下： noint不产生一般在模型中自动生成的截距参数。 none请求全回归模型。 forward 或 f请求向前选择法。7 / 12 backward 或 b请求向后淘汰法。 stepwise请求逐步技术，这个任选项是预置的。 maxr请求最大 R2 增量法。 minr
11、请求最小 R2 增量法。 rsquare请求 R2 最大准则法。 adjrsq请求修正 R2 最大准则法。 cp请求 Mallows 的 Cp 统计量法。 slentry= 值指出向前选择和逐步技术中选择变量进入模型的显著水平。如果省略，那么 stepwise 过程便对向前选择技术置 slentry= 0.5，对逐步技术置 slentry0.15。 slstay= 值指出向后淘汰与逐步技术中变量留在模型里的显著水平。如果省略，则逐步技术用 0.15，向后淘汰技术用 0.10。 include=n强迫头 n 个因变量总是在模型中。选择技术由 model 语句中其他变量来完成。 start= s以
12、含有 model 语句中头 s 个自变量的模型为开始，进行比较、选择过程。理所当然地，没有一个被估计的模型含有不足 s 个的变量。此仅应用于 maxr 或 minr模型。 stop= s当它找到“最佳”s 变量模型之后，stepwise 便停止。其中 s 是 stop 的值，此仅应用于 maxr 或 minr 模型。2) 其他语句 weight 语句用于指出含有观察值的权数的变量。分析中仅用具有 weight 变量正值的观察。 by 语句指定的变量值来分组处理某数据集。六、 实例分析例 33.1 例 32.2 续 对 fitness 数据进行逐步回归分析。调用 reg 过程，model 语句中
13、的参数选项使用 selection=stepwise，请求按逐步回归方法挑选自变量子集。程序如下：proc reg data= fitness ;model oxygen = age weight rstpulse maxpulse runpulse runtime /selection=stepwise ;run ;运行后，得到见表 33.1 所示的结果。表 33.1 逐步回归分析结果8 / 12在输出结果报告中，提供了进入回归变量逐次改变后回归方差分析和拟合的信息。在报告的最后部分，列出了用逐步回归法挑选自变量过程，四个自变量按runtime，age， runpulse，maxpulse
14、先后次序进人回归模型。所有进入回归的变量在 0.15的水平下是显著的，未进人回归的侯选变量在 0.15 的水平下是不显著的。同时还概要地提供了每个回归模型变化时的 R2 值增加值、R 2 值、CP 值、相应的 F 统计量、p 值。在逐步回归的每步细节中，还列出了条件指数的最小值最大值，以及每一个回归变量的类型 2 平方和。age 变量进入模型后，R 2 值的增加值（Partial R2，称为偏 R2 或部分 R2）计算为（650.6657632.9001）/ 851.3815= 0.020867。 如果按 CP 值选择最优子集，随着进入回Stepwise Procedure for Depen
15、dent Variable OXYGENStep 1 Variable RUNTIME Entered R-square = 0.74338010 C(p) = 13.51976469DF Sum of Squares Mean Square F ProbFRegression 1 632.90009985 632.90009985 84.01 0.0001Error 29 218.48144499 7.53384293Total 30 851.38154484Parameter Standard Type IIVariable Estimate Error Sum of Squares F
16、ProbFINTERCEP 82.42177268 3.85530378 3443.36654076 457.05 0.0001RUNTIME -3.31055536 0.36119485 632.90009985 84.01 0.0001Bounds on condition number: 1, 1-Step 2 Variable AGE Entered R-square = 0.76424693 C(p) = 12.22493455DF Sum of Squares Mean Square F ProbFRegression 2 650.66573237 325.33286618 45.
17、38 0.0001Error 28 200.71581247 7.16842187Total 30 851.38154484Parameter Standard Type IIVariable Estimate Error Sum of Squares F ProbFINTERCEP 88.46228749 5.37263885 1943.41070877 271.11 0.0001AGE -0.15036567 0.09551468 17.76563252 2.48 0.1267RUNTIME -3.20395056 0.35877488 571.67750579 79.75 0.0001B
18、ounds on condition number: 1.036941, 4.147763-Step 3 Variable RUNPULSE Entered R-square = 0.81109446 C(p) = 6.82780371DF Sum of Squares Mean Square F ProbFRegression 3 690.55085627 230.18361876 38.64 0.0001Error 27 160.83068857 5.95669217Total 30 851.38154484Parameter Standard Type IIVariable Estima
19、te Error Sum of Squares F ProbFINTERCEP 111.71806443 10.23508836 709.69013814 119.14 0.0001AGE -0.25639826 0.09622892 42.28867438 7.10 0.0129RUNPULSE -0.13090870 0.05059011 39.88512390 6.70 0.0154RUNTIME -2.82537867 0.35828041 370.43528607 62.19 0.0001Bounds on condition number: 1.354763, 11.59745-S
20、tep 4 Variable MAXPULSE Entered R-square = 0.83681815 C(p) = 4.76608569DF Sum of Squares Mean Square F ProbFRegression 4 712.45152692 178.11288173 33.33 0.0001Error 26 138.93001792 5.34346223Total 30 851.38154484Parameter Standard Type IIVariable Estimate Error Sum of Squares F ProbFINTERCEP 98.1478
21、8797 11.78569002 370.57373243 69.35 0.0001AGE -0.19773470 0.09563662 22.84231496 4.27 0.0488MAXPULSE 0.27051297 0.13361978 21.90067065 4.10 0.0533RUNPULSE -0.34810795 0.11749917 46.90088674 8.78 0.0064RUNTIME -2.76757879 0.34053642 352.93569605 66.05 0.0001Bounds on condition number: 8.4182, 76.8513
22、5-All variables left in the model are significant at the 0.1500 level.No other variable met the 0.1500 significance level for entry into the model.Summary of Stepwise Procedure for Dependent Variable OXYGEN Variable Number Partial ModelStep Entered Removed In R*2 R*2 C(p) F ProbF1 RUNTIME 1 0.7434 0
23、.7434 13.5198 84.0076 0.00012 AGE 2 0.0209 0.7642 12.2249 2.4783 0.12673 RUNPULSE 3 0.0468 0.8111 6.8278 6.6959 0.01544 MAXPULSE 4 0.0257 0.8368 4.7661 4.0986 0.05339 / 12归模型中的自变量个数 P 从 2 到 5 个（包括截距） ，相应 CP 值从大到小为13.51976469、12.22493455、6.82780371 和 4.76608569，按照 Mallows 提出的回归模型最优自变量个数的选择准则，CP =4.76
24、608569 是最接近自变量个数 P=5 的模型。CP 的计算公式见(33.11)式，当 P=5 时，CP=138.93001792/5.39197 (3125)= 4.76608569。因此，用逐步回归方法及 CP 值确认的拟合回归模型为：oxygen= 98.147887970.19773470age + 0.27051297maxpulse0.34810795 runpulse2.76757879runtime条件指数（condition number）为最大特征值和每个特征值之比的平方根。我们看到当模型进入第四个自变量 maxpulse 时，最大的条件指数从较小 11.59745 变成
25、了较大 76.85135，说明存在一定程度的共线性，根据前面例 33.2 的分析，我们诊断这个共线性方程可能为runpulsemaxpulse= 0。在向前、向后或逐步回归的变量选择过程中，都有一个判断是否可进入或剔除的显著水平，在程序中是分别由 model 语句的选项 slentry=和 slstay设定的，缺省的情况见表 33.2所示。表 33.2 缺省的入选和剔除显著水平forward backward stepwizeslentry 0.50 0.15slstay 0.10 0.15下面我们提供全部可能回归的程序，并且以 R2 值的大到小排序输出。proc reg data= fitn
26、ess ;model oxygen = age weight rstpulse maxpulse runpulse runtime /selection= rsquare b ;run ;在上述程序中，model 语句的选项 selection= rsquare，表示请求 R2 值最大法，选项 b 是表示要输出每种回归的回归系数。程序运行后，得到见表 33.3 所示的结果。10 / 12表 33.3 用 R2 排序全部可能的变量数的逐步回归分析结果程序的输出包括所有只含一个变量的 6 种回归，含 2 个变量的 15 种回归，。总共有 63 种不同形式的回归模型。例如，含 2 个自变量按 R2
27、第二个大值选择回归模型为，R 2 =0.76142381，拟合的回归模型为oxygen= 93.08880.0735runpulse3.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 所示的结果。N = 31 Regress
28、ion Models for Dependent Variable: OXYGENParameterNumber in R-square EstimatesModel Intercept AGE WEIGHT RSTPULSE MAXPULSE RUNPULSE RUNTIME1 0.74338010 82.4218 . . . . . -3.31061 0.15838344 82.4582 . . . . -0.2068 .1 0.11999670 59.3325 . . -0.2225 . . .1 0.09277653 62.2206 -0.3114 . . . . .1 0.05604
29、592 71.2907 . . . -0.1376 . .1 0.02648849 55.4379 . -0.1041 . . . .-2 0.76424693 88.4623 -0.1504 . . . . -3.20402 0.76142381 93.0888 . . . . -0.0735 -3.14022 0.74522106 86.4665 . . . -0.0256 . -3.2723 5 0.78744812 109.8 -0.2528 -0.0481 -0.0255 -0.0836 . -2.90395 0.52421007 116.1 -0.4747 -0.1531 -0.1
30、567 0.3778 -0.5395 .-6 0.84800313 102.2 -0.2199 -0.0724 -0.00084 0.3047 -0.3732 -2.6805-11 / 12表 33.4 只保留 R2 最大两种情况的逐步回归分析结果通过上面的逐步回归分析，我们已经得到回归模型的自变量个数确定时的最优子集或次优子集，但问题是我们到底应该选择几个自变量的回归模型呢？如上表 33.4 中的 3 个自变量、4 个自变量、5 个自变量、6 个自变量的回归模型中哪一个模型呢？一种最简便确定回归模型的自变量个数的方法是 Mallows 的 Cp 方法。确定好模型的自变量个数后，根据上表 33
31、.4 就很容易确定在这个固定自变量数下，最优的自变量组合和相应的参数值估计。以下的程序是对所有可能的回归按 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/s
32、election=cp adjrsq best=5 ;plot cp. * np. /chocking=red cmallows=bluevaxis=0 to 15 by 2haxis=0 to 8 by 1;run ;Model 语句中的 selection=cp 选项请求计算 Mallows 的 Cp 统计量。选项 adjrsq 表示要显示每种回归模型的统计量 Adj-R2。选项 best=5 表示保留 Cp 值最小的前 5 种。plot 语句中的cp. * np.表达式（注意统计量关键字母后的小圆点）表示 Y 轴为 Cp 值 X 轴为 P 值（P 值包N = 31 Regression
33、Models for Dependent Variable: OXYGENParameterNumber in R-square EstimatesModel Intercept AGE WEIGHT RSTPULSE MAXPULSE RUNPULSE RUNTIME1 0.74338010 82.4218 . . . . . -3.31061 0.15838344 82.4582 . . . . -0.2068 .-2 0.76424693 88.4623 -0.1504 . . . . -3.20402 0.76142381 93.0888 . . . . -0.0735 -3.1402
34、-3 0.81109446 111.7 -0.2564 . . . -0.1309 -2.82543 0.80998844 80.9008 . . . 0.3542 -0.3751 -2.9702-4 0.83681815 98.1479 -0.1977 . . 0.2705 -0.3481 -2.76764 0.81649255 115.7 -0.2764 -0.0493 . . -0.1293 -2.7724-5 0.84800181 102.2 -0.2196 -0.0723 . 0.3049 -0.3734 -2.68255 0.83690359 97.9092 -0.1956 . 0
35、.00678 0.2722 -0.3502 -2.7830-6 0.84800313 102.2 -0.2199 -0.0724 -0.00084 0.3047 -0.3732 -2.6805-12 / 12括截距项） 。plot 语句的选项 chocking=red，表示画 Cp=2PP full 红色参考虚线，其中 P 是子模型中含截距的参数个数，P full 是全模型中不含截距的参数个数。Hoching（1976）建议选择满足 Cp2PP full 且 CpP 的模型。plot 语句的选项 cmallows=blue，表示画 Cp=P 蓝色参考实线，其中 P 是子模型中含截距的参数个数。
36、Mallows（1973）建议考虑所有满足 Cp 较小且接近 P 的模型。这一程序的输出结果见表 33.5 和见图 33-1 所示。表 33.5 按 Cp 由小到大进行排序并保留其前 5 种逐步回归分析结果从输出结果可看出，以 Mallows 的建议为标准，age，maxpulse ，runpulse 和 runtime 四个变量进人回归模型时 Cp 最小 (4.76609)，且与 P=4+1=5 最接近，因为54.76609256=4 。而 Cp=5.00021 的模型满足要求，因为 5.000210.83681815)。不同的标准提供不同的选择结果，这是常有的情况。N = 31 Regre
37、ssion Models for Dependent Variable: OXYGENC(p) R-square Adjusted Variables in ModelIn R-square4.76609 0.83681815 4 0.81171325 AGE MAXPULSE RUNPULSE RUNTIME5.00021 0.84800181 5 0.81760218 AGE WEIGHT MAXPULSE RUNPULSE RUNTIME6.75259 0.83690359 5 0.80428431 AGE RSTPULSE MAXPULSE RUNPULSE RUNTIME6.82780 0.81109446 3 0.79010496 AGE RUNPULSE RUNTIME7.00000 0.84800313 6 0.81000391 AGE WEIGHT RSTPULSE MAXPULSE RUNPULSE RUNTIME图 33-1 带有 Mallows 和 Hocking 参考线的 Cp 散点图
