R语言数据分析与挖掘(第四章):回归分析(3)——变量的选择
R语言数据分析与挖掘(第四章):回归分析(3)——变量的选择
在数据挖掘的实战过程中，经常会遇到变量非常多的情况，即数据的维数很高，也称为“维数灾难”问题。在我们生物医学统计领域，一个数据集中可能存在成百上千个变量，对于回归处模而言，并不是越多变量越好，利用少而精的变量建模显得极为重要，如何选择变量子集就是解决问题的关键。
本文主要介绍几个变量选择的方法:逐步回归、岭回归和lasso回归法。
逐步回归方法
选择变量的最基本方法就是逐步选择，即反复地添加或删除模型中的变量，以达到优化模型的目的，该方法需要确定一个阈值，也就是一个算法停止的标准。逐步回归法就是基于上述思想生成的回归算法，该算法利用AlC准则(赤池信息准则)来度量添加变量的效果，AIC准则的表达式为:
AIC= -2*log(L)+k *edf
其中L代表似然值，edf 代表等效自由度。
R语言中用于实现逐步回归的函数是step()，函数的基本书写格式为:
step(object, scope, scale= 0,direction=c("both", "backward", "forward"),trace= 1, keep = NULL, steps = 1000, k=2, ...)参数介绍:
Object:指定模型的对象，如模型lm;
Scope:指定变量选择的上下界，下界为需要出现在最终模型中的变量组，上界为所有考虑添加到模型中的变量组，若只设置一个公式，则R语言默认其为上界，若需同时设定上下界，则需设置两个公式；
Scale:回归模型和方差分析模型中定义的AIC所需要的值；
Direction:指定变量被添加、移除到模型中或者两者均进行，"forward"即向前法，表示变量被添加，“backward”即向后法，表示变量被移除，“both"表示两者均进行，默认值为“both"；
Trace: 指定是否输出变量选择过程，“0”表示不输出，“I”表示输出，默认值为I;
Keep: 选择从对象中保留的参数的函数，默认值为NULL:
Steps: 指定算法终止的最大迭代次数，默认值为1000;
K:惩罚计算中自由度的倍数，默认值为2。
swiss数据集共有47行观测值，每行有7个变量。
> data(swiss)
> summary(swiss)
Fertility Agriculture Examination Education
Min. :35.00 Min. : 1.20 Min. : 3.00 Min. : 1.00
1st Qu.:64.70 1st Qu.:35.90 1st Qu.:12.00 1st Qu.: 6.00
Median :70.40 Median :54.10 Median :16.00 Median : 8.00
Mean :70.14 Mean :50.66 Mean :16.49 Mean :10.98
3rd Qu.:78.45 3rd Qu.:67.65 3rd Qu.:22.00 3rd Qu.:12.00
Max. :92.50 Max. :89.70 Max. :37.00 Max. :53.00
Catholic Infant.Mortality
Min. : 2.150 Min. :10.80
1st Qu.: 5.195 1st Qu.:18.15
Median : 15.140 Median :20.00
Mean : 41.144 Mean :19.94
3rd Qu.: 93.125 3rd Qu.:21.70
Max. :100.000 Max. :26.60
> corr<-cor(swiss)
> library(corrplot)#如果没有安装，先安装
> corrplot(corr)
上面代码包括Swiss数据集的描述性统计和相关系数计算，并绘制了相关矩阵图。根据数据特征，我们将Fertility作为响应变量，其余变量作为解释变量进行回归分析，然而相关矩阵图显示，解释变量Examination和Education之间的相关性较强，即解释变量之间存在多重共线性，基于这一特点，后续的回归分析将存在一定的问题。
首先对原始数据进行回归分析，将数据中的全部变量用于回归分析，得到的模型称为全模型。
> lm5<-lm(Fertility~.,data=swiss)
> summary(lm5)
Call:
lm(formula = Fertility ~ ., data = swiss)
Residuals:
Min 1Q Median 3Q Max
-15.2743 -5.2617 0.5032 4.1198 15.3213
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 66.91518 10.70604 6.250 1.91e-07 ***
Agriculture -0.17211 0.07030 -2.448 0.01873 *
Examination -0.25801 0.25388 -1.016 0.31546
Education -0.87094 0.18303 -4.758 2.43e-05 ***
Catholic 0.10412 0.03526 2.953 0.00519 **
Infant.Mortality 1.07705 0.38172 2.822 0.00734 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7.165 on 41 degrees of freedom
Multiple R-squared: 0.7067, Adjusted R-squared: 0.671
F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10模型lm5中的变量Examination的参数估计没有通过显著性检验，而与该变量相关性较强的 Education的p值很小，为了解决这一问题，下面利用step（）进行逐步回归。
> tstep<-step(lm5)
Start: AIC=190.69
Fertility ~ Agriculture + Examination + Education + Catholic +
Infant.Mortality
Df Sum of Sq RSS AIC
- Examination 1 53.03 2158.1 189.86
<none> 2105.0 190.69
- Agriculture 1 307.72 2412.8 195.10
- Infant.Mortality 1 408.75 2513.8 197.03
- Catholic 1 447.71 2552.8 197.75
- Education 1 1162.56 3267.6 209.36
Step: AIC=189.86
Fertility ~ Agriculture + Education + Catholic + Infant.Mortality
Df Sum of Sq RSS AIC
<none> 2158.1 189.86
- Agriculture 1 264.18 2422.2 193.29
- Infant.Mortality 1 409.81 2567.9 196.03
- Catholic 1 956.57 3114.6 205.10
- Education 1 2249.97 4408.0 221.43
> summary(tstep)
Call:
lm(formula = Fertility ~ Agriculture + Education + Catholic +
Infant.Mortality, data = swiss)
Residuals:
Min 1Q Median 3Q Max
-14.6765 -6.0522 0.7514 3.1664 16.1422
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 62.10131 9.60489 6.466 8.49e-08 ***
Agriculture -0.15462 0.06819 -2.267 0.02857 *
Education -0.98026 0.14814 -6.617 5.14e-08 ***
Catholic 0.12467 0.02889 4.315 9.50e-05 ***
Infant.Mortality 1.07844 0.38187 2.824 0.00722 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7.168 on 42 degrees of freedom
Multiple R-squared: 0.6993, Adjusted R-squared: 0.6707
F-statistic: 24.42 on 4 and 42 DF, p-value: 1.717e-10上述第一条代码即为逐步回归的操作命令，输出结果展示了变量选择的过程，选择标准是基于AlC值最小:需要注意输出结果的最后一部分，该部分表示逐步回归算法最终选择的变量，可以看出逐步回归在全模型的基础上剔除了变量Examination;利用函数summary()展示逐步回归的具体结果，发现参数估计全部通过了显著性检验，且Adjusted R. squared值为0.6707，说明该模型是有效的。需要注意的是，逐步回归法是基于AIC值最小的选择变量，可能出现挑选变量的参数估计不通过显著性检验的情况，此需要利用函数dropI()进行后线处理，具体操作如下:
> drop1(lm5)
Single term deletions
Model:
Fertility ~ Agriculture + Examination + Education + Catholic +
Infant.Mortality
Df Sum of Sq RSS AIC
<none> 2105.0 190.69
Agriculture 1 307.72 2412.8 195.10
Examination 1 53.03 2158.1 189.86
Education 1 1162.56 3267.6 209.36
Catholic 1 447.71 2552.8 197.75
Infant.Mortality 1 408.75 2513.8 197.03上述输出结果展示了剔除每一个解释变量所对应的AIC值，逐步回归算法选择剔除使AlC值最小的变量，即Examination,若剩余变量的参数估计还是不能通过显著性检验，就根据上述输出结果，选择人为剔除变量Examination外使AIC值最小的变量Agriculture，即利用根据逐步回归选择的变量，然后人为剔除其中的Agriculture,再次进行回归分析。一般来说，此步骤能满足所有参数估计均通过显著性检验:若还是存在参数估计没有通过检验的情况，则需要考虑原始数据可能不适用于线性模型。
岭回归的方法
逐步回归法根据函数lm()来简单拟合模型，缺点在于限定了模型中的变量个数，岭回归就能较好地解决这一问题，下面将详细介绍岭回归法的操作步骤。
岭回归法的思想是:对系数的个数设置约束，并使用不同的算法来拟合模型，以缓解数据内部的多重共线性所带来的方差变大等问题。之前已经介绍了基于最小化残差平方和的参数估计法，即最小二乘法，岭回归则是对每个参数添加一个惩罚项，基于最小化残差平方和与系数的惩罚项总和，一般来说，系数的惩罚项总和是系数平方和的倍数，具体如下:
岭回归的目的就是寻找使RSS最小时的参数估计，在R中，包MASS中的函数lm.ridgc(可以满足要求，函数的基本书写格式为:
Im.ridge(formula, data, subset, na.action, lambda=0, model=FALSE,
x=FALSE, y = FALSE, contrasts = NULL,...)参数介绍:
Formula:指定用于拟合的模型公式，类似于Im中的用法;
Data:指定用于做岭回归的数据对象，可以是数据框、列表或者能强制转换为数据框的其他数据对象:
Subset:一个向量，指定数据中需要包含在模型中的观测值:
Na.action:一个函数，指定当数据中存在缺失值时的处理办法，用法与Im中的一致:
Lambda:指定RSS的表达式中系数平方和的倍数项，默认值为0;
Model:逻辑值，指定是否返回“模型框架”，默认值为FALSE:
X:逻辑值，指定是否返回“模型矩阵”，默认值为FALSE:
Y:逻辑值，制度能够是否返回响应变量，默认值为FALSE:
Contrasts:模型中因子对照的列表。
lm6<-lm(Employed ~., data=longley)
summary(lm6)
library(car)
vif(lm6)
install.packages("MASS")
library(MASS)
plot(lm.ridge(Employed ~ ., longley,lambda = seq(0,0.1,0.001)))
select(lm.ridge(GNP.deflator ~ ., longley,lambda = seq(0,0.1,0.0001)))
plot(lm.ridge(Employed ~ .-GNP.deflator, longley,lambda = seq(0,0.1,0.001)))
select(lm.ridge(GNP.deflator ~ .-GNP.deflator, longley,lambda = seq(0,0.1,0.0001)))
plot(lm.ridge(Employed ~ .-GNP.deflator-GNP, longley,lambda = seq(0,0.1,0.001)))
select(lm.ridge(GNP.deflator ~ .-GNP.deflator-GNP, longley,lambda = seq(0,0.1,0.0001)))
lm7<-lm(Employed ~.-GNP.deflator-GNP, data=longley)
summary(lm7)lasso回归方法
另一种选择变量的回归算法是lasso,与岭回归类似的是，lasso 也添加了针对回由参数的惩罚项，不同之处在于lasso选择的惩罚方式是:用绝对值的平方和取代系数平方和，其RSS的表达式为:
lasso的目的就是寻找使RSS最小时的参数估计，在R语言中，包lars中的函数lasr()可以满足要求，其函数的基本书写格式为:
lars(x, y，type = c("lasso", "lar", "forward.stagewise","stepwise"),
trace = FALSE, normalize = TRUE, intercept = TRUE, Gram,
eps = .Machine$double.eps, max.steps, use.Gram = TRUE)参数介绍:
x:一个矩阵，用于指定预测变量:
y:一个向量，用于指定响应变量;
Type:指定拟合模型的类型，"so"表示进行lasso回归，"lar表示进行最小角回归,"foward. sgewse表示进行极小向前逐段回归，"epis"表示进行遂步回归，默认值为"lasso";
Trace:逻辑值，指定是否打印函数运行过程中的详细信息，默认值为FALSE;
Normalize:逻辑值，指定是否将所有变量，默认值为TRUE;
Intercept:逻辑值，指定是否将解决项包含在模型中,默认值为TRUE;
Gram: 计算过程中的x'x矩阵;
Eps: 有效的0值:
Max.steps:算法迭代的最大次数;
use.Gram: 逻辑位，指定是否预先计算Gram矩阵，默认值为TRUE;
install.packages("lars")
library(lars)
lar.1<-lars(dlongley[,1:6],dlongley[,7])
lar.1
summary(lar.1)
plot(lar.1)
cva <- cv.lars(dlongley[,1:6],dlongley[,7], K = 10, plot.it = TRUE)
best <- cva$index[which.min(cva$cv)]
coef <- coef.lars(lar.1, mode = "fraction", s = best) #使得CV最小时的系数
coef[coef != 0]
best <- cva$index[which.min(cva$cv)]
coef <- coef.lars(lar.1, mode = "fraction", s = best) #使得CV最小时的系数
coef[coef != 0] #通过CV选择的变量本文分享自 MedBioInfoCloud 微信公众号，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体分享计划 ，欢迎热爱写作的你一起参与！