变量筛选：逐步回归&最优子集

影响因素比较多的时候，变量间常常存在着各种各样的关系，单因素的分析结果往往不可靠，我们常常采用的解决方案是多因素的回归分析。而在进行多因素的回归分析时，如何处理无统计学意义的变量是避不开的问题：保留还是舍弃？这就涉及到变量筛选的问题。

这需要回到研究的开始：研究的目的是什么？是为了校正主要研究因素的混杂因素，还是为了寻找结局的影响因素？是为了建一个预测模型还是为了因果推断？如果你只是想校正其他因素的影响，显然这个时候不需要进行变量的筛选。我们在进行多因素的方差分析时，即使某个因素没有统计学意义，我们也不会因此把这个没有统计学意义的变量去掉。但如果是为了寻找研究结局的影响因素，变量筛选就很重要的一个过程。如果你是为了建立一个预测模型，最终入选模型的变量也未必都具有统计学意义，能达到当前条件下的最佳的预测效果就可以了。如果是用来推断因果则可能寻求纳入有统计学意义的变量。另外，我觉得样本量也是一个重要的考量。

考虑到样本量，在进行多因素筛选前，常常先进行单因素分析，将单因素分析中具有统计学意义或接近有统计学意义的变量纳入模型，这就是我们常说的先单后多。这也是变量筛选的一种方法，可以算是初步筛选。我们今天说的变量筛选是进行多变量分析时如何舍弃或者保留变量的问题。

变量筛选不是简单地将多因素回归中无统计学意义的变量删除，而是遵循需要一些特殊的方法。变量筛选专业判断第一，统计方法有逐步回归（先前选择法、向后剔除法、双向选择法）、最优子集筛选、lasso回归等。

Lasso回归题目比较大，后面单独笔记。全子集回归我们在<<R笔记：全子集回归 | 最优子集筛选>>已做过演示。今天调用R程序包StepReg来演示逐步回归，该程序包可以实现线性回归、logistic回归、COX风险比例回归的逐步回归，选择变量的方法有向前选择、向后剔除、双向选择以及最优子集选择。

Stepwise Linear Model Regression：

stepwise(formula,data,include = NULL,

selection = c("forward", "backward", "bidirection", "score"),

select = c("AIC", "AICc", "BIC", "CP", "HQ", "HQc", "Rsq", "adjRsq", "SL", "SBC"),

sle = 0.15,sls = 0.15,

multivarStat = c("Pillai", "Wilks", "Hotelling-Lawley", "Roy"),

weights = NULL,best = NULL)

Stepwise Logistic Regression

stepwiseLogit(formula,data,include = NULL,

selection = c("forward", "backward", "bidirection", "score"),

select = c("SL", "AIC", "AICc", "SBC", "HQ", "HQc", "IC(3/2)", "IC(1)"),

sle = 0.15,sls = 0.15,

sigMethod = c("Rao", "LRT"),

weights = NULL,best = NULL)

Stepwise Cox Proportional Hazards Regression

stepwiseCox(formula,data,include= NULL,

selection = c("forward", "backward", "bidirection", "score"),

select = c("SL", "AIC", "AICc", "SBC", "HQ", "HQc", "IC(3/2)", "IC(1)"),

sle = 0.15,sls = 0.15,

method = c("efron", "breslow", "exact"),

weights = NULL,best = NULL)

formula指定进行筛选的模型，data指定包含模型变量的数据集,include可以设定强制纳入模型的变量。当selection选定forward、backward、bidirection表示将对formula指定的模型进行变量选择的方法：向前选择法、向后剔除法、向前向后选择法，选择score时表示将采用全子集回归。当预测变量较多时，全子集回归可生成大量组合模型，可通过best参数设置来控制显示模型的数量。select参数设置的是每一步变量进入模型或从模型剔除的标准，包装各种信息准则和显著水平（SL，Signifificant Levels），线性回归中还会有决定系数和校正的决定系数。sle和sls指定的是进入模型的显著水平（signifificance level for entry）和保留在模型中的显著水平（signifificance level for staying in the model）。线性逐步回归中可以进行多元（多个因变量）回归的逐步回归，可以通过multivarStat指定多元回归的统计量，其对应的formul表达式为d(y1,y2) ~ x1 + x2。Logistic回归中的sigMethod指定变量进入模型的显著性检验方法，包括Rao氏效率得分检验（Rao’s effificient score test）和似然比检验（likelihood ratio test）。Cox回归中的method可以指定处理打结数据的方法。

示例：同<<因变量二分类资料的Probit回归>>中的数据。

library(haven)

fitdata<-read_sav("D:/Temp/logistic_step.sav")

fitdata<-as.data.frame(fitdata)

fitdata<-na.omit(fitdata)

fitdata$race<-factor(fitdata$race,levels=c(1,2,3),labels=c("White","Black","Others"))

fitdata$smoke<-factor(fitdata$smoke,levels=c(0,1),labels=c("No","Yes"))

fitdata$ptl <-ifelse(fitdata$ptl>0,1, 0)

##考虑到样本量的问题，早产次数ptl宜合并为二分类变量。当原ptl值>0时,新的ptl赋值为0，其他情况下(原ptl<=0)时则赋值为0

fitdata$ptl<-factor(fitdata$ptl,levels=c(0,1),labels=c("No","Yes"))

fitdata$ht<-factor(fitdata$ht,levels=c(0,1),labels=c("No","Yes"))

fitdata$ui<-factor(fitdata$ui,levels=c(0,1),labels=c("No","Yes"))

【6】变量筛选

library(StepReg)

stepwise(bwt~age+lwt+race+smoke+ptl+ht+ui+ftv,data=fitdata,selection="bidirection",select="SL",sle=0.05,sls=0.1) 

##双向逐步回归，变量进入标准P=0.05，剔除标准P=0.1

结果解读参见前面的推文，此不赘述！

本例如果采用全子集回归，命令如下：

stepwise(bwt~age+lwt+race+smoke+ptl+ht+ui+ftv,data=fitdata,selection="score",select="AIC",best=2)

全子集回归可以生成256个模型，我们设置了参数best=2，即变量数相同的模型仅显示2个AIC较小的模型。结果显示，当纳入7个变量（age、lwt、race、smoke、ptl、ht和ui）时AIC最小，接下来我们用这7个变量进行线性回归即可。

lmfit<- lm(bwt~age+lwt+race+smoke+ptl+ht+ui,data=fitdata)

summary(lmfit)

如果我们以新生儿是否是低体重出生儿（low：0=正常，1=低体重出生儿）为因变量，可以建立二分类的logistic回归，变量筛选的逐步回归如下：

stepwiseLogit(low~age+lwt+race+smoke+ptl+ht+ui+ftv,data=fitdata,selection="bidirection",select="SL",sle=0.05,sls=0.1,sigMethod="LRT") 

向前向后选择法筛选了3个变量：ptl、ht和lwt。以此三个变量建立logistic回归获取更多信息。

lgrfit<-glm(low~lwt+ptl+ht, family=binomial(link = logit), data=fitdata)

summary(lgrfit)

confint(lgrfit) ##给出模型参数的置信区间

exp(cbind(OR=coef(lgrfit),confint(lgrfit))) ##给出模型参数的OR和置信区间

                   OR     2.5 %     97.5 %

(Intercept) 2.2366523 0.4530008 12.4304398

lwt         0.9846164 0.9711283  0.9968232

ptlYes      4.0705166 1.7788055  9.6045086

htYes       6.1852152 1.6008819 27.3558354

具体解读可参照《R笔记：二分类资料的logistic回归》，此不赘述！

— END —

原创不易，欢迎“在看”

关注“一统浆糊”

获取更多信息