经典线性回归模型自变量选择.docx
《经典线性回归模型自变量选择.docx》由会员分享，可在线阅读，更多相关《经典线性回归模型自变量选择.docx（14页珍藏版）》请在冰点文库上搜索。
经典线性回归模型自变量选择
§2.9自变量选择
信息时代的一个重要特征是数据便宜信息值钱，我们经常要从海量数据中挖掘有用信息。
比如影响产品质量的因素，从生产过程、员工培训过程到原材料供应过程，可能多达几百个，甚至上千个。
对这些质量指标和影响因素制造商在日常生产管理过程中都有记录。
现在的问题是如何从这众多的影响因素中找出影响产品质量的重要因素。
有时只需判断一个自变量对因变量是否有重要影响，而不需要了解它们之间的精确定量关系。
比如判断原材料供应对产品质量是否有重要影响比了解它们之间的精确定量关系更重要。
线性回归模型的自变量选择就是用于有众多自变量时识别重要自变量的方法。
用于线性回归模型自变量选择的方法可分为两类：
全局择优法和逐步回归法。
一、全局择优法
全局择优法就是用衡量回归模型与数据拟合程度的准则，从全部可能的回归模型中选择对数据拟合最优的回归模型。
对于一个包含P个自变量的回归问题，全部可能的回归模型有
个，全局择优法要求出每个回归模型的准则值，然后找出最优的回归模型。
回归模型对数据的拟合程度可用残差平方和来表示。
残差平方和越小，模型拟合的越好。
但残差平方和的大小与因变量的计量单位有关，因此我们定义了决定系数。
决定系数越大，模型拟合的越好。
决定系数不仅与因变量的计量单位无关，而且能说明在因变量的变异中，归功于自变量变化的部分所占比例。
但不论是用残差平方和还是用决定系数来度量线性拟合模型拟合程度，都会得出模型中包含越多自变量拟合就越好的结论。
但在样本容量给定的情况下，自变量越多，模型就越复杂，模型参数估计就越不精确，导致模型应用的效果就越差。
因此我们需要能综合用残差平方和表示的模型拟合精度和用模型中包含的自变量个数表示的模型复杂程度的准则，以便选择出最优的回归模型。
回归分析中用于选择自变量的准则很多。
由于残差平方和RSSp和决定系数R2只考虑模型拟合精度，因而只能作为自变量个数相同时自变量选择的准则。
残差均方s2和修正决定系数
是一个综合模型拟合精度和模型复杂程度的准则。
综合性准则除了残差均方和修正决定系数外，还有如下一些准则：
·MallowsCp准则
其中，s2为包含全部自变量的拟合模型的残差均方，RSSp为当前拟合模型的残差平方和，p为当前拟合模型的自变量个数。
·信息准则
信息准则根据公式
计算，其中logLik=-n{log（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，计算当前拟合模型的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（）和regsubsets（）均可用来完成全局最优的选择。
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=Slim
x6=Shld
x7=Lane
x8=Acpt
x9=Itg
x10=Lwid
x11=as.numeric（Hwy==1）
x12=as.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：
其它
Rate1973年每百万公里行车的事故率
考虑log（Rate）对log（Len）,log（ADT）,log（Trks）,log（Sigs1）,Slim,Shld,Lane,Acpt,Itg,Lwid,Hwy的回归,其中Sigs1=（Sigs*Len+1）/Len
hwm=lm（y~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13）
summary（hwm）
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）3.9541541.7655412.2400.0342*
x1-0.2144700.099986-2.1450.0419*
x2-0.1546250.111893-1.3820.1792
x3-0.1975600.239812-0.8240.4178
x40.1923220.0753672.5520.0172*
x5-0.0272600.016799-1.6230.1172
x60.0029740.0341590.0870.9313
x7-0.0111330.057021-0.1950.8468
x80.0060490.0081010.7470.4622
x90.0357220.2428180.1470.8842
x100.0421220.1368210.3080.7607
x110.2375450.3998220.5940.5578
x12-0.2857810.273072-1.0470.3053
x13-0.1437290.233458-0.6160.5437
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
0.2607on25degreesoffreedom
MultipleR-squared:
0.7913,AdjustedR-squared:
0.6828
F-statistic:
7.293on13and25DF,p-value:
1.247e-05
kappa（hwm）
[1]1305.216
library（car）
vif（hwm）
x1x2x3x4x5x6x7
2.0330758.3451491.9745914.9948645.3974856.0143523.365510
x8x9x10x11x12x13
3.1854075.5724292.17496110.25124810.6888766.949116
library（leaps）
leaps（x=cbind（x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13）,y=y,nbest=1）
$which
123456789ABCD
1FALSEFALSEFALSEFALSETRUEFALSEFALSEFALSEFALSEFALSEFALSEFALSEFALSE
2TRUEFALSEFALSEFALSETRUEFALSEFALSEFALSEFALSEFALSEFALSEFALSEFALSE
3FALSEFALSEFALSETRUETRUEFALSEFALSEFALSEFALSEFALSEFALSETRUEFALSE
4TRUEFALSEFALSETRUETRUEFALSEFALSEFALSEFALSEFALSEFALSETRUEFALSE
5TRUETRUEFALSETRUETRUEFALSEFALSEFALSEFALSEFALSEFALSETRUEFALSE
6TRUETRUEFALSETRUETRUEFALSEFALSEFALSEFALSEFALSEFALSETRUETRUE
7TRUETRUETRUETRUETRUEFALSEFALSEFALSEFALSEFALSEFALSETRUETRUE
8TRUETRUETRUETRUETRUEFALSEFALSETRUEFALSEFALSEFALSETRUETRUE
9TRUETRUETRUETRUETRUEFALSEFALSETRUEFALSEFALSETRUETRUETRUE
10TRUETRUETRUETRUETRUEFALSEFALSETRUEFALSETRUETRUETRUETRUE
11TRUETRUETRUETRUETRUEFALSETRUETRUEFALSETRUETRUETRUETRUE
12TRUETRUETRUETRUETRUEFALSETRUETRUETRUETRUETRUETRUETRUE
13TRUETRUETRUETRUETRUETRUETRUETRUETRUETRUETRUETRUETRUE
$label
[1]"（Intercept）""1""2""3""4""5"
[7]"6""7""8""9""A""B"
[13]"C""D"
$size
[1]234567891011121314
$Cp
[1]27.72323310.2020837.0867931.6498531.4379302.0455123.0904174.5426566.179744
[10]8.06356810.02710212.00758214.000000
leaps（）函数的一般用法为
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=x3,x4=x4,x5=x5,x6=x6,x7=x7,x8=x8,x9=x9,x10=x10,x11=x11,x12=x12,x13=x13,y=y）
a=regsubsets（y~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13,data=hihway）
summary（a）
Subsetselectionobject
Call:
regsubsets.formula（y~x1+x2+x3+x4+x5+x6+x7+x8+
x9+x10+x11+x12+x13,data=hihway）
13Variables（andintercept）
ForcedinForcedout
x1FALSEFALSE
x2FALSEFALSE
x3FALSEFALSE
x4FALSEFALSE
x5FALSEFALSE
x6FALSEFALSE
x7FALSEFALSE
x8FALSEFALSE
x9FALSEFALSE
x10FALSEFALSE
x11FALSEFALSE
x12FALSEFALSE
x13FALSEFALSE
1subsetsofeachsizeupto8
SelectionAlgorithm:
exhaustive
x1x2x3x4x5x6x7x8x9x10x11x12x13
1
（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","forward","seqrep"））
其中，formula用来指定包含全部自变量的模型；
data=用来指定存放数据的数据框；
nbest=用来指定要报告的不同个数自变量的局部最优模型数。
nvmax=用来指定最大模型的自变量个数；
force.in=用来指定强制进入模型的自变量；
force.out=用来指定强制剔除的自变量；
intercept=表示模型是否包含常数项的逻辑值，缺省表示包含常数项；
method=用来指定选优的方法，包括全局、向后、向前和逐步，缺省为全局。
二、逐步回归法
全局择优法需要大量的运算。
当有5个自变量时，所有可能的回归数为25-1=15个；当有10个自变量时，所有可能的回归数为210-1=1023个；当有50个自变量时，所有可能的回归数为250-1，大约是1015个。
因此在自变量个数较多时，全局择优法是无法实现的，此时需要别一类自变量选择方法，逐步回归法。
逐步回归法分为向前选择、向后剔除和逐步筛选三种。
·向前选择
（a）从不含自变量的回归模型开始；
（b）依据某个标准从候选的自变量中选择一个最优的自变量添加到模型中；
（c）直到候选自变量中没有符合标准的自变量可添加为止。
·向后剔除
（a）从包含全部自变量的回归模型开始；
（b）依据某个标准从模型中剔除一个最差的自变量；
（c）直到模型中没有符合标准的变量可剔除为止。
·逐步筛选
（a）从任意一个回归模型开始；
（b）依据某个标准从候选的自变量中选择一个最优的自变量添加到模型中，或者依
据某个标准从模型中剔除一个最差的自变量；
（c）直到既没有符合标准的候选自变量可添加，模型中也没有符合标准的自变量可
剔除为止。
以上逐步回归法中选择自变量的标准既可以用衡量回归模型与数据拟合程度的准则，也可以用检验系数显著性的t统计量、F统计量或者P值来构造。
R函数step（）可用于逐步回归方法，这个函数的一般用法为
step（object,scope,scale=0,direction=c（"both","backward","forward"）,k=2）
其中，object指定逐步回归的初始模型；
Scope指定逐步回归搜索的模型范围。
如果是包含lower和upper两个公式的列
表，则lower指定强制包含在模型中的自变量（这些自变量必须包含在初始
模型中），upper指定最大的模型。
如果是单个公式，则表示最大的模型。
如果是缺省，则初始模型为最大的模型；
如果scale=s2，则用Cp准则，缺省表示使用信息准则；
direction=指定逐步回归方法，缺省为逐步筛选，"forward"为向前选择、
"backward"为向后剔除；
如果k=log（n），则用BIC准则，缺省表示使用AIC准则。
m0=lm（y~1,data=hihway）
m1<-step（m0,~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13,data=hihway）
summary（m1）
Call:
lm（formula=y~x5+x1+x12+x4+x2,data=hihway）
Residuals:
Min1QMedian3QMax
-0.61322-0.164110.001490.125510.41288
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）3.9010190.4347898.9722.28e-10***
x5-0.0253440.009053-2.8000.00848**
x1-0.2645640.081428-3.2490.00266**
x12-0.3221780.092721-3.4750.00145**
x40.1864710.0527463.5350.00123**
x2-0.0756440.047215-1.6020.11866
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
0.242on33degreesoffreedom
MultipleR-squared:
0.7626,AdjustedR-squared:
0.7267
F-statistic:
21.21on5and33DF,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:
Min1QMedian3QMax
-0.61322-0.164110.001490.125510.41288
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）3.9010190.4347898.9722.28e-10***
x1-0.2645640.081428-3.2490.00266**
x2-0.0756440.047215-1.6020.11866
x40.1864710.0527463.5350.00123**
x5-0.0253440.009053-2.8000.00848**
x12-0.3221780.092721-3.4750.00145**
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
0.242on33degreesoffreedom
MultipleR-squared:
0.7626,AdjustedR-squared:
0.7267
F-statistic:
21.21on5and33DF,p-value:
1.901e-09
m3=step（hwm,scale=summary（hwm）$sigma^2）
summary（m3）
Call:
lm（formula=y~x1+x2+x4+x5+x12）
Residuals:
Min1QMedian3QMax
-0.61322-0.164110.001490.125510.41288
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）3.9010190.4347898.9722.28e-10***
x1-0.2645640.081428-3.2490.00266**
x2-0.0756440.047215-1.6020.11866
x40.1864710.0527463.5350.00123**
x5-0.0253440.009053-2.8000.00848**
x12-0.3221780.092721-3.4750.00145**
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
0.2
升级会员 