§2.9 自变量选择

1、.2.9 自变量选择信息时代的一个重要特征是数据便宜信息值钱，我们经常要从海量数据中挖掘有用信息。比如影响产品质量的因素，从生产过程、员工培训过程到原材料供应过程，可能多达几百个，甚至上千个。对这些质量指标和影响因素制造商在日常生产管理过程中都有记录。现在的问题是如何从这众多的影响因素中找出影响产品质量的重要因素。有时只需判断一个自变量对因变量是否有重要影响，而不需要了解它们之间的精确定量关系。比如判断原材料供应对产品质量是否有重要影响比了解它们之间的精确定量关系更重要。线性回归模型的自变量选择就是用于有众多自变量时识别重要自变量的方法。用于线性回归模型自变量选择的方法可分为两类：全局择优法和

2、逐步回归法。一、全局择优法全局择优法就是用衡量回归模型与数据拟合程度的准则，从全部可能的回归模型中选择对数据拟合最优的回归模型。对于一个包含P个自变量的回归问题，全部可能的回归模型有个，全局择优法要求出每个回归模型的准则 值，然后找出最优的回归模型。回归模型对数据的拟合程度可用残差平方和来表示。残差平方和越小，模型拟合的越好。但残差平方和的大小与因变量的计量单位有关，因此我们定义了决定系数。决定系数越大，模型拟合的越好。决定系数不仅与因变量的计量单位无关，而且能说明在因变量的变异中，归功于自变量变化的部分所占比例。但不论是用残差平方和还是用决定系数来度量线性拟合模型拟合程度，都会得出模型中包含

3、越多自变量拟合就越好的结论。但在样本容量给定的情况下，自变量越多，模型就越复杂，模型参数估计就越不精确，导致模型应用的效果就越差。因此我们需要能综合用残差平方和表示的模型拟合精度和用模型中包含的自变量个数表示的模型复杂程度的准则，以便选择出最优的回归模型。回归分析中用于选择自变量的准则很多。由于残差平方和RSSp和决定系数R2只考虑模型拟合精度，因而只能作为自变量个数相同时自变量选择的准则。残差均方s2和修正决定系数是一个综合模型拟合精度和模型复杂程度的准则。综合性准则除了残差均方和修正决定系数外，还有如下一些准则：Mallows Cp准则其中，s2为包含全部自变量的拟合模型的残差均方，RSS

4、p为当前拟合模型的残差平方和，p为当前拟合模型的自变量个数。信息准则信息准则根据公式计算，其中logLik= -nlog(RSS/n)+log(2)+1/2为当前拟合模型的对数似然函数，npar为当前拟合模型的参数个数，当k=2时称为AIC准则，当k=log(n)时称为BIC准则。在小样本情况下，AIC准则的表现不太好，为此人们提出的修正AIC准则AICc，其计算公式为 R中计算当前拟合模型信息准则的函数有（其中fit为当前拟合模型对象）AIC(fit, k = 2) k=2(缺省)时计算 k= log(n)时计算extractAIC(fit, scale, k = 2) 指定scale=s2

5、，计算当前拟合模型的Cp准则 不指定scale，k=2(缺省)时计算 不指定scale，k=log(n)时计算 R的附加程序包qpcR中的函数AICc(fit)可计算当前拟合模型的修正信息准则预测平方和准则其中，表示删除第i个案例后，用剩余的(n-1)个案例估计的拟合模型对第i个案例的预测误差。R的附加程序包qpcR中的函数PRESS(fit)可计算预测平方和。此函数的返回值是一个列表，其中包含三个元素，(1) 名字为stat的预测平方和；(2) 名字为residuals的预测残差向量；(3) 名字为P.square的P2，其计算公式为：R的的附加程序包leaps中的函数leaps( )和re

6、gsubsets( )均可用来完成全局最优的选择。leaps( )依据Cp准则、修正R2准则和R2准则来选择全局最优回归模型；regsubsets( )函数则只能选出不同自变量个数的局部最优的模型，我们再从这些局部的最优模型中选出全局最优的模型。例：高速公路事故数据library(alr3)attach(highway)y=log(Rate)x1=log(Len)x2=log(ADT)x3=log(Trks)x4=log(Sigs*Len+1)/Len)x5=Slimx6=Shldx7=Lanex8=Acpt x9=Itgx10=Lwidx11=as.numeric(Hwy=1)x12=as.

7、numeric(Hwy=2)x13=as.numeric(Hwy=3) 考虑汽车意外事故率（事故数/百万行车）与一些可能的相关之间的关系。数据包括1973年在明尼苏达州的39段高速公路。ADT 以千计的平均是流量（估计）Trks 卡车容量在全部容量中的百分比Lane 在两个方向上的交通车道总数Acpt 路段中每英里的进入点Sigs 路段中每英里信号交换数Itg 路段中每英里的快车道类型交换数Slim 时速限制（在1973年）Len 段的长度（英里）Lwid 道路宽度（英尺）Shld 道路的外侧路肩宽度Hwy 公路类型的因子变量，0：州际高速公路、1：首要干道高速公路、2：主干道高速公路、3：其

8、它Rate 1973年每百万公里行车的事故率考虑log(Rate)对log(Len), log(ADT), log(Trks), log(Sigs1), Slim, Shld, Lane, Acpt, Itg, Lwid, Hwy的回归, 其中Sigs1=(Sigs*Len+1)/Lenhwm=lm(yx1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13)summary(hwm) Estimate Std. Error t value Pr(|t|) (Intercept) 3.954154 1.765541 2.240 0.0342 *x1 -0.214470

9、 0.099986 -2.145 0.0419 *x2 -0.154625 0.111893 -1.382 0.1792 x3 -0.197560 0.239812 -0.824 0.4178 x4 0.192322 0.075367 2.552 0.0172 *x5 -0.027260 0.016799 -1.623 0.1172 x6 0.002974 0.034159 0.087 0.9313 x7 -0.011133 0.057021 -0.195 0.8468 x8 0.006049 0.008101 0.747 0.4622 x9 0.035722 0.242818 0.147 0

10、.8842 x10 0.042122 0.136821 0.308 0.7607 x11 0.237545 0.399822 0.594 0.5578 x12 -0.285781 0.273072 -1.047 0.3053 x13 -0.143729 0.233458 -0.616 0.5437 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 Residual standard error: 0.2607 on 25 degrees of freedomMultiple R-squared: 0.7913, Adjusted R-squared

11、: 0.6828 F-statistic: 7.293 on 13 and 25 DF, p-value: 1.247e-05kappa(hwm)1 1305.216library(car)vif(hwm) x1 x2 x3 x4 x5 x6 x7 2.033075 8.345149 1.974591 4.994864 5.397485 6.014352 3.365510 x8 x9 x10 x11 x12 x13 3.185407 5.572429 2.174961 10.251248 10.688876 6.949116library(leaps)leaps(x=cbind(x1,x2,x

12、3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13), y=y, nbest=1)$which 1 2 3 4 5 6 7 8 9 A B C D1 FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE2 TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE3 FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FA

13、LSE TRUE FALSE4 TRUE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE5 TRUE TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE6 TRUE TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE7 TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE T

14、RUE TRUE8 TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE FALSE FALSE FALSE TRUE TRUE9 TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE FALSE FALSE TRUE TRUE TRUE10 TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE FALSE TRUE TRUE TRUE TRUE11 TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE12TRUE TRUE

15、TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE13TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE$label 1 (Intercept) 1 2 3 4 5 7 6 7 8 9 A B 13 C D $size 1 2 3 4 5 6 7 8 9 10 11 12 13 14$Cp 1 27.723233 10.202083 7.086793 1.649853 1.437930 2.045512 3.090417 4.542656 6.1797441

16、0 8.063568 10.027102 12.007582 14.000000leaps( )函数的一般用法为leaps(x=, y=, method=c(Cp, adjr2, r2),int=TRUE, nbest=10, names=NULL)其中，x=用来指定自变量的矩阵；y=用来指定因变量的向量； method=用来指定准则，缺省为Cp准则； int=表示模型是否包含常数项的逻辑值，缺省为TRUE，表示包含常数项。nbest=用来指定要报告的不同个数自变量的局部最优模型数。 names=用来指定自变量名称的字符向量。hihway=data.frame(x1=x1,x2=x2,x3=x

17、3,x4=x4,x5=x5,x6=x6,x7=x7,x8=x8,x9=x9,x10=x10,x11=x11,x12=x12,x13=x13,y=y)a=regsubsets(yx1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13,data=hihway)summary(a)Subset selection objectCall: regsubsets.formula(y x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + x12 + x13, data = hihway)13 Variables (a

18、nd intercept) Forced in Forced outx1 FALSE FALSEx2 FALSE FALSEx3 FALSE FALSEx4 FALSE FALSEx5 FALSE FALSEx6 FALSE FALSEx7 FALSE FALSEx8 FALSE FALSEx9 FALSE FALSEx10 FALSE FALSEx11 FALSE FALSEx12 FALSE FALSEx13 FALSE FALSE1 subsets of each size up to 8Selection Algorithm: exhaustive x1 x2 x3 x4 x5 x6

19、x7 x8 x9 x10 x11 x12 x131 ( 1 ) * 2 ( 1 ) * * 3 ( 1 ) * * * 4 ( 1 ) * * * * 5 ( 1 ) * * * * * 6 ( 1 ) * * * * * *7 ( 1 ) * * * * * * *8 ( 1 ) * * * * * * * *regsubsets( )函数的一般用法为regsubsets(formula, data=, nbest=1, nvmax=8, force.in=NULL, force.out=NULL, intercept=TRUE, method=c(exhaustive, backward,

20、 forward, seqrep)其中，formula用来指定包含全部自变量的模型； data=用来指定存放数据的数据框； nbest=用来指定要报告的不同个数自变量的局部最优模型数。 nvmax=用来指定最大模型的自变量个数； force.in=用来指定强制进入模型的自变量； force.out=用来指定强制剔除的自变量； intercept=表示模型是否包含常数项的逻辑值，缺省表示包含常数项； method=用来指定选优的方法，包括全局、向后、向前和逐步，缺省为全局。二、逐步回归法全局择优法需要大量的运算。当有5个自变量时，所有可能的回归数为25-1=15个；当有10个自变量时，所有可能的

21、回归数为210-1=1023个；当有50个自变量时，所有可能的回归数为250-1，大约是1015个。因此在自变量个数较多时，全局择优法是无法实现的，此时需要别一类自变量选择方法，逐步回归法。逐步回归法分为向前选择、向后剔除和逐步筛选三种。向前选择(a) 从不含自变量的回归模型开始；(b) 依据某个标准从候选的自变量中选择一个最优的自变量添加到模型中；(c) 直到候选自变量中没有符合标准的自变量可添加为止。向后剔除(a) 从包含全部自变量的回归模型开始；(b) 依据某个标准从模型中剔除一个最差的自变量；(c) 直到模型中没有符合标准的变量可剔除为止。逐步筛选(a) 从任意一个回归模型开始；(b)

22、 依据某个标准从候选的自变量中选择一个最优的自变量添加到模型中，或者依据某个标准从模型中剔除一个最差的自变量；(c) 直到既没有符合标准的候选自变量可添加，模型中也没有符合标准的自变量可剔除为止。 以上逐步回归法中选择自变量的标准既可以用衡量回归模型与数据拟合程度的准则，也可以用检验系数显著性的t统计量、F统计量或者P值来构造。R函数step( )可用于逐步回归方法，这个函数的一般用法为step(object, scope, scale = 0,direction = c(both, backward, forward), k = 2)其中，object指定逐步回归的初始模型；Scope指定逐

23、步回归搜索的模型范围。如果是包含lower和upper两个公式的列表，则lower指定强制包含在模型中的自变量（这些自变量必须包含在初始模型中），upper指定最大的模型。如果是单个公式，则表示最大的模型。如果是缺省，则初始模型为最大的模型； 如果scale =s2，则用Cp准则，缺省表示使用信息准则；direction =指定逐步回归方法，缺省为逐步筛选，forward为向前选择、backward为向后剔除；如果k= log(n)，则用BIC准则，缺省表示使用AIC准则。m0=lm(y1,data=hihway)m1|t|) (Intercept) 3.901019 0.434789 8.9

24、72 2.28e-10 *x5 -0.025344 0.009053 -2.800 0.00848 * x1 -0.264564 0.081428 -3.249 0.00266 * x12 -0.322178 0.092721 -3.475 0.00145 * x4 0.186471 0.052746 3.535 0.00123 * x2 -0.075644 0.047215 -1.602 0.11866 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 Residual standard error: 0.242 on 33 degrees of

25、 freedomMultiple R-squared: 0.7626, Adjusted R-squared: 0.7267 F-statistic: 21.21 on 5 and 33 DF, p-value: 1.901e-09 m=lm(y.,data=insu1)m2=step(hwm)summary(m2)Call:lm(formula = y x1 + x2 + x4 + x5 + x12)Residuals: Min 1Q Median 3Q Max -0.61322 -0.16411 0.00149 0.12551 0.41288 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 3.901019 0.434789 8.972 2.28e-10 *x1 -0.264564 0.081428 -3.249 0.00266 * x2 -0.075644 0.047215 -1.602 0.11866 x4 0.186471 0.052746 3.535 0.00123 * x5 -0.025344 0.009053 -2.800 0.00848 * x12 -0.322178 0.092721 -3.475 0.00145 * -Signif. codes: 0 * 0.