逐步回归
向前引入法:从一元回归开始,逐步加快变量,使指标值达到最优为止
向后剔除法:从全变量回归方程开始,逐步删去某个变量,使指标值达到最优为止
逐步筛选法:综合上述两种方法
多元线性回归的核心问题:应该选择哪些变量?
- RSS(残差平方和)与R2(相关系数平方)选择法:遍历所有可能的组合,选出使RSS最小,R2最大的模型
- AIC(Akaike information criterion)准则和BIC(Bayesian information criterion)准则
AIC=n×ln(RSSP/n)+2p
n为变量总个数,p为选出的变量个数,AIC越小越好
那么关于逐步回归可使用step()函数,省略自行编写代码的繁琐。
sl=step(s,direction="forward")#向后剔除法 sl=step(s,direction="backward")#向前引入法
例题:某种水泥在凝固时放出的热量Y(卡/克)与水泥中四种化学成分X1,X2,X3,X4有关,现测得13组数据,如下表。希望从中选出主意的变量,建立Y关于它们的线性回归方程
序号 | X1 | X2 | X3 | X4 | Y | 序号 | X1 | X2 | X3 | X4 | Y |
1 | 7 | 26 | 6 | 60 | 78.5 | 8 | 1 | 31 | 22 | 44 | 72.5 |
2 | 1 | 29 | 15 | 52 | 74.3 | 9 | 2 | 54 | 18 | 22 | 93.1 |
3 | 11 | 56 | 8 | 20 | 104.3 | 10 | 21 | 47 | 4 | 26 | 115.9 |
4 | 11 | 31 | 8 | 47 | 87.6 | 11 | 1 | 40 | 23 | 34 | 83.8 |
5 | 7 | 52 | 6 | 33 | 95.9 | 12 | 11 | 66 | 9 | 12 | 113.3 |
6 | 11 | 55 | 9 | 22 | 109.2 | 13 | 10 | 68 | 8 | 12 | 109.4 |
7 | 3 | 71 | 17 | 6 | 102.7 |
首先作多元线性回归方程
> cement<-data.frame( + X1=c( 7,1, 11, 11,7, 11,3,1,2, 21,1, 11, 10), + X2=c(26, 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68), + X3=c( 6, 15,8,8,6,9, 17, 22, 18,4, 23,9,8), + X4=c(60, 52, 20, 47, 33, 22,6, 44, 22, 26, 34, 12, 12), + Y =c(78.5, 74.3, 104.3,87.6,95.9, 109.2, 102.7, 72.5,93.1,115.9,83.8, 113.3, 109.4) + ) > lm.sol<-lm(Y ~ X1+X2+X3+X4, data=cement) > summary(lm.sol) Call: lm(formula = Y ~ X1 + X2 + X3 + X4, data = cement) Residuals: Min 1Q Median 3Q Max -3.1750 -1.6709 0.2508 1.3783 3.9254 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 62.4054 70.0710 0.891 0.3991 X1 1.5511 0.7448 2.083 0.0708 . X2 0.5102 0.7238 0.705 0.5009 X3 0.1019 0.7547 0.135 0.8959 X4 -0.1441 0.7091 -0.203 0.8441 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.446 on 8 degrees of freedom Multiple R-squared: 0.9824, Adjusted R-squared: 0.9736 F-statistic: 111.5 on 4 and 8 DF, p-value: 4.756e-07
从上述计算汇总可以看到,如果选择全部变量作回归方程,效果不太好,因为回归方程系数没有一项通过
下面用step()函数作逐步回归
> lm.step<-step(lm.sol) Start: AIC=26.94 Y ~ X1 + X2 + X3 + X4 Df Sum of Sq RSS AIC - X3 1 0.1091 47.973 24.974 - X4 1 0.2470 48.111 25.011 - X2 1 2.9725 50.836 25.728 <none> 47.864 26.944 - X1 1 25.9509 73.815 30.576 Step: AIC=24.97 Y ~ X1 + X2 + X4 Df Sum of Sq RSS AIC <none> 47.97 24.974 - X4 1 9.93 57.90 25.420 - X2 1 26.79 74.76 28.742 - X1 1 820.91 868.88 60.629
从程序运行结果可看到,用全部变量作回归方程时,AIC值为26.94,接下来显示的数据表告诉你,如果去掉变量X3,得到回归方程的AIC值为24.974.如果去掉变量X4,得到回归方程的AIC值为25.011.后面的类推。由于去掉变量X3可以是AIC达到最小,因此,R自动去掉变量X3,进行下一轮计算。
在下一轮计算中,无论去掉哪一个变量,AIC值均会升高,因此R终止计算,得到“最优”回归方程。
下面分析计算结果
> summary(lm.step) Call: lm(formula = Y ~ X1 + X2 + X4, data = cement) Residuals: Min 1Q Median 3Q Max -3.0919 -1.8016 0.2562 1.2818 3.8982 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 71.6483 14.1424 5.066 0.000675 *** X1 1.4519 0.1170 12.410 5.78e-07 *** X2 0.4161 0.1856 2.242 0.051687 . X4 -0.2365 0.1733 -1.365 0.205395 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.309 on 9 degrees of freedom Multiple R-squared: 0.9823, Adjusted R-squared: 0.9764 F-statistic: 166.8 on 3 and 9 DF, p-value: 3.323e-08
由结果看到:回归系数检验的显著性水平有很大提高,但变量X2,X4系数检验的显著性水平仍不理想。
在R中,还有两个函数可以用作逐步回归。add1()和drop1()。使用格式
add1(object, scope, ...) drop1(object, scope, ...) add1(object, scope, scale=0, test=c("none", "Chisq"),k=2, trace=FALSE, ...) drop1(object, scope, scale=0, test=c("none", "Chisq"),k=2, trace=FALSE, ...) add1(object, scope, scale=0, test=c("none", "Chisq", "F"),x=NULL, k=2, ...) drop1(object, scope, scale=0, all.cols=TRUE,test=c("none", "Chisq", "F"), k=2, ...)
其中object是由拟合模型构成的对象。scope是模型考虑增加或去掉项构成的公式。scale是用于计算Cp的残差的均方估计值,缺省值为0或NULL。其他见help
下面用drop1()计算
> drop1(lm.step) Single term deletions Model: Y ~ X1 + X2 + X4 Df Sum of Sq RSS AIC <none> 47.97 24.974 X1 1 820.91 868.88 60.629 X2 1 26.79 74.76 28.742 X4 1 9.93 57.90 25.420
从运算结果来看,如果去掉变量X4,AIC值从24.97增加到25.42,是增加的最少的。另外,除AIC准则外,残差的平方和也是逐步回归的重要指标之一,从直观来看,拟合越好的方程,残差的平方和应越小。去掉变量X4,残差的平方和上升9.93,也是最少的。因此,从这两项指标来看,应该再去掉变量X4
> lm.opt<-lm(Y ~ X1+X2, data=cement); summary(lm.opt) Call: lm(formula = Y ~ X1 + X2, data = cement) Residuals: Min 1Q Median 3Q Max -2.893 -1.574 -1.302 1.363 4.048 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 52.57735 2.28617 23.00 5.46e-10 *** X1 1.46831 0.12130 12.11 2.69e-07 *** X2 0.66225 0.04585 14.44 5.03e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.406 on 10 degrees of freedom Multiple R-squared: 0.9787, Adjusted R-squared: 0.9744 F-statistic: 229.5 on 2 and 10 DF, p-value: 4.407e-09
这个结果应该是满意的,因为所有的检验均是显著。最后得到“最优”回归方程为
Yˆ=52.58+1.468X1+0.6622X2
(这编辑器真难用)