在实际问题中,影响因变量Y的因素很多,人们可以从中挑选若干个变量建立回归方程,这便涉及变量选择的问题。
一般来讲,如果在一个回归方程中忽略了对Y有显著影响的自变量,那么所建立的方程必与实际有较大的偏离,但变量选得过多,使用就不方便,特别是方程中含有对Y影响不大的变量时。因此适当地选择变量以建立一个“最优”的回归方程是十分重要的。
什么是“最优”回归方程呢?对于这个问题有许多不同的准则,在不同的准则下,“最优”回归方程也可能不同,这讲的“最优”是指从可供选择的所有变量中选出对Y有显著影响的变量建立方程,且在方程中不含有对Y无显著影响的变量。
有许多方法准则可以获得“最优”回归方程,如“一切子集回归法”、"前进法"、“后退法”、“逐步回归法”等,其中“逐步回归法”由于计算程序简便,因而使用较为普遍。
R语言提供了较为方便的“逐步回归”计算函数step(),它是以AIC信息统计量为准则,通过选择最小的AIC信息统计量,来达到删除或增加变量的目的。step()函数的调用格式为:
其中object是回归模型,scope是确定逐步搜索的区域,scale是用于AIC统计量,direction确定逐步搜索的方向,缺省值为“both”,是“一切子集回归法”,“backward”是“后退法”,“forward”是“前进法”。下面通过一个简单的例子,介绍如何使用R语言来完成逐步回归的过程,从而达到选择“最优”方程的目的。
例子:某种水泥在凝固时放出的热量Y与水泥中四种化学成分X1、X2、X3、X4有关,现测得13组数据,求从中选出主要的变量,建立Y关于它们的线性回归方程。
首先作多元线性回归方程
> 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);
> cement<-data.frame(X1,X2,X3,X4,Y)
> 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,进行下一轮计算。
> 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(),它们的调用格式为:
其中object是由拟合模型构成的对象,scope是模型考虑增加或去掉项构成的公式,scale是用于计算残差的均方估计值,缺省值为0或NULL。下面用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.974增加到25.420,是增加最少的。另外,除了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.468*X1 0.6623*X2
,