1、创建虚拟变量
caret包通过model.matrix()、dummyVars()或其他方式将因子转换为虚拟变量。
dummyVars()可用于从一个或多个因子变量生成一个完整的虚拟变量集。该函数获取一个公式和一个数据集，并输出一个对象，该对象创建虚拟变量可用于预测。
> library(pacman)
> p_load(earth, caret)
> data(etitanic)
> str(etitanic)
## 'data.frame': 1046 obs. of 6 variables:
## $ pclass : Factor w/ 3 levels "1st","2nd","3rd": 1 1 1 1 1 1 1 1 1 1 ...
## $ survived: int 1 1 0 0 0 1 1 0 1 0 ...
## $ sex : Factor w/ 2 levels "female","male": 1 2 1 2 1 2 1 2 1 2 ...
## $ age : num 29 0.917 2 30 25 ...
## $ sibsp : int 0 1 1 1 1 0 1 0 2 0 ...
## $ parch : int 0 2 2 2 2 0 0 0 0 0 ...
该数据集包括两个因子变量，pclass表示乘客等级，sex表示性别。
看看R基本函数model.matrix()生成的结果：
> lm.mat <- model.matrix(survived ~ ., data = etitanic)
> head(lm.mat, 3)
## (Intercept) pclass2nd pclass3rd sexmale age sibsp parch
## 1 1 0 0 0 29.0000 0 0
## 2 1 0 0 1 0.9167 1 2
## 3 1 0 0 0 2.0000 1 2
再看看dummyVars()函数的结果：
> dummies <- dummyVars(survived ~ ., data = etitanic)
> pred.dummies <- predict(dummies, newdata = etitanic)
> head(pred.dummies, 3)
## pclass.1st pclass.2nd pclass.3rd sex.female sex.male age sibsp parch
## 1 1 0 0 1 0 29.0000 0 0
## 2 1 0 0 0 1 0.9167 1 2
## 3 1 0 0 1 0 2.0000 1 2
注意没有截距项，每个级别的因子都有一个虚拟变量，所以这个参量化的函数对于某些模型可能没有用处，比如lm()。
2、零方差和接近零方差变量的预测
在某些情况下，数据生成机制可以创建只有一个唯一值的预测因子(即“零方差预测器”)。 对于许多模型(不包括基于树的模型) ，这可能会导致模型崩溃或拟合得不健壮。
> data(mdrr)
> table(mdrrDescr$nR11)
##
## 0 1 2
## 501 4 23
数据集中可能只有少数值出现的频率很低，非常的不平衡。当使用交叉验证建模时，这些样本可能对模型的影响很小而成为零方差预测因子。
nearZeroVar()函数用来识别常数自变量，或者是方差极小的自变量：
> # saveMetrics默认为F，用来显示细节
> nzv <- nearZeroVar(mdrrDescr,
+ # 第一众数与第二众数的比率的cutoff(临界值)
+ freqCut = 95/5,
+ # 剔重后的唯一值 与 样本总数量的百分比,大于这个值不会被剔除
+ uniqueCut = 10,
+ saveMetrics = T)
> str(nzv)
## 'data.frame': 342 obs. of 4 variables:
## $ freqRatio : num 1.25 1.12 1 1.25 1.25 ...
## $ percentUnique: num 90 42.6 83 84.3 82.8 ...
## $ zeroVar : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ nzv : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
> # 显示前10个接近零方差的变量
> nzv[nzv$nzv,][1:10,]
## freqRatio percentUnique zeroVar nzv
## nTB 23.00000 0.3787879 FALSE TRUE
## nBR 131.00000 0.3787879 FALSE TRUE
## nI 527.00000 0.3787879 FALSE TRUE
## nR03 527.00000 0.3787879 FALSE TRUE
## nR08 527.00000 0.3787879 FALSE TRUE
## nR11 21.78261 0.5681818 FALSE TRUE
## nR12 57.66667 0.3787879 FALSE TRUE
## D.Dr03 527.00000 0.3787879 FALSE TRUE
## D.Dr07 123.50000 5.8712121 FALSE TRUE
## D.Dr08 527.00000 0.3787879 FALSE TRUE
> dim(mdrrDescr)
## [1] 528 342
> index <- nearZeroVar(mdrrDescr)
> filter.Descr <- mdrrDescr[, -index]
> dim(filter.Descr)
## [1] 528 297
可以看到，数据由342个变量删减为297个变量了。
3、识别与其他变量有强相关性的变量
> # 建立相关性矩阵
> descr.cor <- cor(filter.Descr)
> summary(descr.cor[upper.tri(descr.cor)])
查看相关性统计信息：
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.99607 -0.05373 0.25006 0.26078 0.65527 1.00000
> # 找出相关性高于0.75的变量
> high.cor <- findCorrelation(descr.cor, cutoff = 0.75)
> length(high.cor)
## [1] 247
一共247个变量，与其他变量的相关性高于0.75。
> low.cor <- filter.Descr[, -high.cor]
> descr.cor.2 <- cor(low.cor)
> summary(descr.cor.2[upper.tri(descr.cor.2)])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.70728 -0.05378 0.04418 0.06692 0.18858 0.74458
剔除后再看看相关性统计信息，全部在0.75以下了。
4、识别变量共线性
> mat <- matrix(0, nrow = 6, ncol = 6)
> mat[, 1] <- c(1, 1, 1, 1, 1, 1)
> mat[, 2] <- c(1, 1, 1, 0, 0, 0)
> mat[, 3] <- c(0, 0, 0, 1, 1, 1)
> mat[, 4] <- c(1, 0, 0, 1, 0, 0)
> mat[, 5] <- c(0, 1, 0, 0, 1, 0)
> mat[, 6] <- c(0, 0, 1, 0, 0, 1)
请注意，第二列和第三列相加就等于第一列。同样，第四、五和六列加起来就是第一列。Findlinearcombos()函数将返回一个列表，列举这些依赖项。对于每个线性组合，它都会递增地从矩阵中删除列，并测试依赖关系是否已被解析。还将返回一个列位置的向量，以消除线性依赖:
> combo.info <- findLinearCombos(mat)
> combo.info
## $linearCombos
## $linearCombos[[1]]
## [1] 3 1 2
##
## $linearCombos[[2]]
## [1] 6 1 4 5
##
##
## $remove
## [1] 3 6
> mat <- mat[, -combo.info$remove]
> mat
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 0
## [2,] 1 1 0 1
## [3,] 1 1 0 0
## [4,] 1 0 1 0
## [5,] 1 0 0 1
## [6,] 1 0 0 0
5、preProcess函数
缺省参数是标准化数据，其高级功能还包括用K近邻和装袋决策树两种方法来预测缺失值。此外它还可以进行cox幂变换和主成分提取。
method为处理方法，可能的值为"BoxCox", "YeoJohnson", "expoTrans", "center", "scale", "range", "knnImpute", "bagImpute", "medianImpute", "pca", "ica", "spatialSign", "corr", "zv", "nzv", "conditionalX"。
BoxCox 变换用于响应变量的转换，对于估计幂变换也同样有效
YeoJohnson 转换类似于 Box-Cox 模型，但是可以容纳零和 / 或负值的预测器(BoxCox 转换的预测器值必须严格为正)
center 预测值减去预测值数据的平均值
scale 预测值除以预测值数据的标准差
range 转换扩展了范围内的数据。 如果新样本的值大于或小于训练集中的值，则值将超出此范围。在计算中将忽略非数值型预测值。
zv 用单个值（即方差为零）标识数值型预测列，并将其排除在进一步计算之外
nzv 通过应用nearZeroVar函数标识方差接近0的变量，并将其排除在进一步计算之外
除非使用zv或者nzv方法，否则任何数值型变量少于2个唯一值时函数将报错
corr 通过findCorrelation函数过滤掉高度相关的预测因子
conditionalX 对于分类，根据结果检查每个预测因子的分布情况。如果某类中只有一个唯一值，则预测器将被排除在进一步计算之外。如果预测的不是因子变量，则不执行此计算。 这个操作在通过train函数进行重采样时会很耗时。
运算按如下顺序进行: zero-variance filter, near-zero variance filter, correlation filter, Box-Cox/Yeo-Johnson/exponential transformation, centering, scaling, range, imputation, PCA, ICA then spatial sign.
使用PCA和ICA方法时，数据将被自动中心化和标准化。如果同时使用PCA和ICA方法，则会抛出警告。Ica由fastICA包实现，在查找ICA得分之前自动进行PCA分解
注意，当使用PCA或ICA时，非数值型变量在预测时可能位于不同的位置。非数值型数据不会被预处理，但是他们的值会保留在预测值中。
knnImpute 通过训练集的K最近邻（欧式距离）填充缺失值
bagImpute 通过装袋法填充缺失值，该方法简单、准确、可接受缺失值，但是计算量大
medianImpute 通过训练集中预测值的中位数填充缺失值，该方法简单、快速、可接受缺失值，但可能不准确
5.1 中心化和标准化
> set.seed(120)
> ind <- sample(2, nrow(low.cor), replace = T, prob = c(0.7, 0.3))
>
> train <- low.cor[ind == 1, ]
> test <- low.cor[ind == 2, ]
>
> dim(train)
## [1] 367 50
> dim(test)
## [1] 161 50
> # 同时使用多种预处理方法
> preprocess.values <- preProcess(train,
+ method = c("center", "scale", "YeoJohnson", "nzv"))
> preprocess.values
## Created from 367 samples and 50 variables
##
## Pre-processing:
## - centered (50)
## - ignored (0)
## - scaled (50)
## - Yeo-Johnson transformation (31)
##
## Lambda estimates for Yeo-Johnson transformation:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.70101 -0.41917 0.04850 -0.07775 0.35785 1.88606
> train.transformed <- predict(preprocess.values, train)
> test.transformed <- predict(preprocess.values, test)
>
> train[1:3, 1:3]
## AMW Mp Ms
## METHOPROMAZINE 7.15 0.67 1.92
## ACEPROMAZINE 7.26 0.68 2.06
## TRIMEPRAZINE 6.94 0.68 1.85
> train.transformed[1:3, 1:3]
## AMW Mp Ms
## METHOPROMAZINE 0.322580832 0.5994524 -0.9359891
## ACEPROMAZINE 0.480970227 0.9717011 -0.3131278
## TRIMEPRAZINE 0.007559779 0.9717011 -1.2474198
5.2 特征选择
在进行数据挖掘时，我们并不需要将所有的自变量用来建模，而是从中选择若干最重要的变量，这称为特征选择（feature selection）。一种算法是后向选择，即先将所有的变量都包括在模型中，然后计算其效能（如误差、预测精度）和变量重要排序，然后保留最重要的若干变量，再次计算效能，这样反复迭代，找出合适的自变量数目。这种算法的一个缺点在于可能会存在过度拟合，所以需要在此算法外再套上一个样本划分的循环。在caret包中的rfe()函数可以完成这项任务。
> # functions确定用什么样的模型进行自变量排序，本例选择的模型是随机森林即rfFuncs
> # 可以选择的还有lmFuncs（线性回归），nbFuncs（朴素贝叶斯），treebagFuncs（装袋决策树），caretFuncs(使用caret包的训练模型）
> # method确定用什么样的抽样方法，本例使用cv即交叉检验
> # 还有提升boot以及留一交叉检验LOOCV
> ctrl <- rfeControl(functions = rfFuncs,
+ # 确定抽样方法，默认为“boot”
+ method = "repeatedcv",
+ verbose = F,
+ # 交叉验证折数或抽样迭代次数，默认10和25
+ number = 20,
+ returnResamp = "final")
>
> # 使用随机森林算法，设置多线程（linux系统）
> if (require("multicore", quietly = TRUE, warn.conflicts = FALSE)) {
+ ctrl$workers <- multicore:::detectCores()
+ ctrl$computeFunction <- mclapply
+ ctrl$computeArgs <- list(mc.preschedule = FALSE, mc.set.seed = FALSE)
+ }
>
> # 特征选择
> profile <- rfe(
+ # 训练集自变量矩阵或数据库，注意，列名必须唯一
+ train.transformed,
+ # 训练集的结果向量（数值型或因子型）
+ mdrrClass[ind==1],
+ # 指定选择最优模型的指标，默认分类用Acc，回归用RMSE
+ metric = "Accuracy",
+ # metric是否最大化，"Accuracy"为T，RMSE为F
+ maximize = T,
+ # 控制项列表，包括拟合预测的函数
+ rfeControl = ctrl,
+ rerank=F)
> profile
##
## Recursive feature selection
##
## Outer resampling method: Cross-Validated (20 fold, repeated 1 times)
##
## Resampling performance over subset size:
##
## Variables Accuracy Kappa AccuracySD KappaSD Selected
## 4 0.7744 0.5374 0.07025 0.1434
## 8 0.8234 0.6376 0.07440 0.1559
## 16 0.8311 0.6544 0.06564 0.1292 *
## 50 0.8311 0.6539 0.05745 0.1134
##
## The top 5 variables (out of 16):
## X5v, VRA1, H3D, G1, QXXm
最终选择了16个变量。画图看看重要性排序：
> plot(profile)
可以看到16个变量时效能最高。返回最终保留的自变量：
> profile$optVariables
## [1] "X5v" "VRA1" "H3D" "G1" "QXXm" "Xt" "nBnz" "Jhetm" "SPAM"
## [10] "TIE" "nAB" "FDI" "VEA1" "BLI" "SPI" "SPH"
