Для начала подключим пакет caret и посмотрим, какие модели в нём доступны.
library(caret)
## Loading required package: ggplot2
## Loading required package: lattice
names(getModelInfo())
## [1] "ada" "AdaBag" "AdaBoost.M1"
## [4] "adaboost" "amdai" "ANFIS"
## [7] "avNNet" "awnb" "awtan"
## [10] "bag" "bagEarth" "bagEarthGCV"
## [13] "bagFDA" "bagFDAGCV" "bam"
## [16] "bartMachine" "bayesglm" "binda"
## [19] "blackboost" "blasso" "blassoAveraged"
## [22] "bridge" "brnn" "BstLm"
## [25] "bstSm" "bstTree" "C5.0"
## [28] "C5.0Cost" "C5.0Rules" "C5.0Tree"
## [31] "cforest" "chaid" "CSimca"
## [34] "ctree" "ctree2" "cubist"
## [37] "dda" "deepboost" "DENFIS"
## [40] "dnn" "dwdLinear" "dwdPoly"
## [43] "dwdRadial" "earth" "elm"
## [46] "enet" "evtree" "extraTrees"
## [49] "fda" "FH.GBML" "FIR.DM"
## [52] "foba" "FRBCS.CHI" "FRBCS.W"
## [55] "FS.HGD" "gam" "gamboost"
## [58] "gamLoess" "gamSpline" "gaussprLinear"
## [61] "gaussprPoly" "gaussprRadial" "gbm_h2o"
## [64] "gbm" "gcvEarth" "GFS.FR.MOGUL"
## [67] "GFS.LT.RS" "GFS.THRIFT" "glm.nb"
## [70] "glm" "glmboost" "glmnet_h2o"
## [73] "glmnet" "glmStepAIC" "gpls"
## [76] "hda" "hdda" "hdrda"
## [79] "HYFIS" "icr" "J48"
## [82] "JRip" "kernelpls" "kknn"
## [85] "knn" "krlsPoly" "krlsRadial"
## [88] "lars" "lars2" "lasso"
## [91] "lda" "lda2" "leapBackward"
## [94] "leapForward" "leapSeq" "Linda"
## [97] "lm" "lmStepAIC" "LMT"
## [100] "loclda" "logicBag" "LogitBoost"
## [103] "logreg" "lssvmLinear" "lssvmPoly"
## [106] "lssvmRadial" "lvq" "M5"
## [109] "M5Rules" "manb" "mda"
## [112] "Mlda" "mlp" "mlpKerasDecay"
## [115] "mlpKerasDecayCost" "mlpKerasDropout" "mlpKerasDropoutCost"
## [118] "mlpML" "mlpSGD" "mlpWeightDecay"
## [121] "mlpWeightDecayML" "monmlp" "msaenet"
## [124] "multinom" "mxnet" "mxnetAdam"
## [127] "naive_bayes" "nb" "nbDiscrete"
## [130] "nbSearch" "neuralnet" "nnet"
## [133] "nnls" "nodeHarvest" "null"
## [136] "OneR" "ordinalNet" "ordinalRF"
## [139] "ORFlog" "ORFpls" "ORFridge"
## [142] "ORFsvm" "ownn" "pam"
## [145] "parRF" "PART" "partDSA"
## [148] "pcaNNet" "pcr" "pda"
## [151] "pda2" "penalized" "PenalizedLDA"
## [154] "plr" "pls" "plsRglm"
## [157] "polr" "ppr" "pre"
## [160] "PRIM" "protoclass" "qda"
## [163] "QdaCov" "qrf" "qrnn"
## [166] "randomGLM" "ranger" "rbf"
## [169] "rbfDDA" "Rborist" "rda"
## [172] "regLogistic" "relaxo" "rf"
## [175] "rFerns" "RFlda" "rfRules"
## [178] "ridge" "rlda" "rlm"
## [181] "rmda" "rocc" "rotationForest"
## [184] "rotationForestCp" "rpart" "rpart1SE"
## [187] "rpart2" "rpartCost" "rpartScore"
## [190] "rqlasso" "rqnc" "RRF"
## [193] "RRFglobal" "rrlda" "RSimca"
## [196] "rvmLinear" "rvmPoly" "rvmRadial"
## [199] "SBC" "sda" "sdwd"
## [202] "simpls" "SLAVE" "slda"
## [205] "smda" "snn" "sparseLDA"
## [208] "spikeslab" "spls" "stepLDA"
## [211] "stepQDA" "superpc" "svmBoundrangeString"
## [214] "svmExpoString" "svmLinear" "svmLinear2"
## [217] "svmLinear3" "svmLinearWeights" "svmLinearWeights2"
## [220] "svmPoly" "svmRadial" "svmRadialCost"
## [223] "svmRadialSigma" "svmRadialWeights" "svmSpectrumString"
## [226] "tan" "tanSearch" "treebag"
## [229] "vbmpRadial" "vglmAdjCat" "vglmContRatio"
## [232] "vglmCumulative" "widekernelpls" "WM"
## [235] "wsrf" "xgbDART" "xgbLinear"
## [238] "xgbTree" "xyf"
Сохраним в переменные x и y случайные данные согласно примеру
в задании. Значения генерируются случайно, поэтому
set.seed(1) фиксирует их, чтобы графики не менялись при
каждой сборке отчёта. Столбцам x даём имена, чтобы они были подписаны на
графиках.
set.seed(1)
x <- matrix(rnorm(50*5), ncol = 5)
colnames(x) <- paste0("V", 1:5)
y <- factor(rep(c("A", "B"), 25))
Построим диаграммы на основе данных x и y. Каждый
график выводится в отчёт и сохраняется в .jpg файл. Пока открыт
jpeg(), график рисуется в файл, а dev.off()
закрывает файл.
проверим, куда будут сохраняться изображения:
getwd()
## [1] "C:/Users/mixai/Downloads"
boxplot:
p <- featurePlot(x, y, plot = "box")
print(p) # в отчёт
jpeg("box.jpg", width = 800, height = 600)
print(p) # в файл
invisible(dev.off())
densityPlot
p <- featurePlot(x, y, plot = "density",
scales = list(x = list(relation = "free"),
y = list(relation = "free")),
auto.key = list(columns = 2))
print(p)
jpeg("density.jpg", width = 800, height = 600)
print(p)
invisible(dev.off())
stripPlot
p <- featurePlot(x, y, plot = "strip", jitter = TRUE)
print(p)
jpeg("strip.jpg", width = 800, height = 600)
print(p)
invisible(dev.off())
pairsPlot
p <- featurePlot(x, y, plot = "pairs", auto.key = list(columns = 2))
print(p)
jpeg("pairs.jpg", width = 800, height = 600)
print(p)
invisible(dev.off())
ellipsePlot
p <- featurePlot(x, y, plot = "ellipse", auto.key = list(columns = 2))
print(p)
jpeg("ellipse.jpg", width = 800, height = 600)
print(p)
invisible(dev.off())
Выводы: Список names(getModelInfo())
содержит несколько сотен методов, доступных в caret; у многих из них
отбор признаков встроен в сам алгоритм (например, у lasso и случайного
леса). На всех пяти графиках распределения признаков в классах A и B
сильно перекрываются: на boxplot медианы и размахи почти совпадают, на
графиках плотности кривые двух классов накладываются друг на друга, на
strip-графике точки классов перемешаны, а на pairs и ellipse облака
точек и эллипсы классов почти совпадают. Ни один из пяти признаков не
разделяет классы. Это ожидаемо, поскольку данные сгенерированы случайно
и никак не связаны с классом. Небольшие различия (например, чуть более
высокая медиана признака у одного из классов) случайны и при такой малой
выборке (50 наблюдений) значимыми считать их нельзя.
Датасет “iris”
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
Подключим пакет FSelector и посчитаем важность признаков разными методами:
library(FSelector)
w_ig <- information.gain(Species ~ ., data = iris) #Прирост информации
w_gr <- gain.ratio(Species ~ ., data = iris)
w_su <- symmetrical.uncertainty(Species ~ ., data = iris)
w_chi <- chi.squared(Species ~ ., data = iris)
w_one <- oneR(Species ~ ., data = iris)
set.seed(1) # случайный лес использует случайные числа
w_rf <- random.forest.importance(Species ~ ., data = iris, importance.type = 1)
print(w_ig)
## attr_importance
## Sepal.Length 0.4521286
## Sepal.Width 0.2672750
## Petal.Length 0.9402853
## Petal.Width 0.9554360
Сравним в таблице важность четырех признаков с помощью всех 6-ти методов
fs_methods <- list(IG = information.gain, GR = gain.ratio,
SU = symmetrical.uncertainty, Chi2 = chi.squared,
OneR = oneR, RF = random.forest.importance)
set.seed(1)
res <- sapply(fs_methods, function(f) {
w <- f(Species ~ ., data = iris)
setNames(w$attr_importance, rownames(w))
})
res <- res[names(iris)[1:4], ] # одинаковый порядок признаков
round(res, 3) # исходные значения
## IG GR SU Chi2 OneR RF
## Sepal.Length 0.452 0.420 0.416 0.629 0.173 15.593
## Sepal.Width 0.267 0.247 0.245 0.492 0.040 7.283
## Petal.Length 0.940 0.858 0.857 0.935 0.400 48.076
## Petal.Width 0.955 0.871 0.871 0.943 0.407 45.889
res_norm <- apply(res, 2, function(v) v / max(v))
round(res_norm, 2) # нормированные (максимум = 1)
## IG GR SU Chi2 OneR RF
## Sepal.Length 0.47 0.48 0.48 0.67 0.43 0.32
## Sepal.Width 0.28 0.28 0.28 0.52 0.10 0.15
## Petal.Length 0.98 0.99 0.98 0.99 0.98 1.00
## Petal.Width 1.00 1.00 1.00 1.00 1.00 0.95
res
## IG GR SU Chi2 OneR RF
## Sepal.Length 0.4521286 0.4196464 0.4155563 0.6288067 0.1733333 15.593303
## Sepal.Width 0.2672750 0.2472972 0.2452743 0.4922162 0.0400000 7.282986
## Petal.Length 0.9402853 0.8584937 0.8571872 0.9346311 0.4000000 48.075563
## Petal.Width 0.9554360 0.8713692 0.8705214 0.9432359 0.4066667 45.888977
barplot(t(res_norm), beside = TRUE, las = 1,
legend.text = TRUE, args.legend = list(x = "topleft", cex = 0.7),
main = "Важность признаков iris по разным методам")
Выводы: Все шесть методов (information gain, gain ratio, symmetrical uncertainty, хи-квадрат, OneR, случайный лес) дали один и тот же порядок важности: признаки лепестка Petal.Width и Petal.Length значительно важнее признаков чашелистика Sepal.Length и Sepal.Width. Например, по information gain лепестковые признаки получили 0,955 и 0,940, а чашелистиковые 0,452 и 0,267. Самым слабым признаком все методы называют Sepal.Width. Внутри пары лепестка различия минимальны (пять методов ставят выше Petal.Width, случайный лес чуть выше оценивает Petal.Length), поэтому один из них безусловно лучшим считать нельзя. Значения у методов в разных шкалах, поэтому сравнивать можно только порядок признаков, а не сами числа. Для различения видов ириса достаточно признаков лепестка. Поскольку они сильно коррелируют, использованные методы оценивают их по отдельности и избыточности не видят, поэтому на практике можно оставить один из них.
Применим функцию разными методами к одному признаку датасета iris,
Sepal.Length, и разобьём его на 3 категории. Функция
table() покажет, сколько наблюдений попало в каждую
категорию.
library(arules)
## Loading required package: Matrix
##
## Attaching package: 'arules'
## The following objects are masked from 'package:base':
##
## abbreviate, write
x <- iris$Sepal.Length
set.seed(1) # cluster использует k-means со случайным стартом
d_interval <- discretize(x, method = "interval", breaks = 3)
d_frequency <- discretize(x, method = "frequency", breaks = 3)
d_cluster <- discretize(x, method = "cluster", breaks = 3)
d_fixed <- discretize(x, method = "fixed",
breaks = c(-Inf, 5.5, 6.5, Inf),
labels = c("short", "medium", "long"))
lapply(list(interval = d_interval, frequency = d_frequency,
cluster = d_cluster, fixed = d_fixed), table)
## $interval
##
## [4.3,5.5) [5.5,6.7) [6.7,7.9]
## 52 70 28
##
## $frequency
##
## [4.3,5.4) [5.4,6.3) [6.3,7.9]
## 46 53 51
##
## $cluster
##
## [4.3,5.33) [5.33,6.27) [6.27,7.9]
## 46 53 51
##
## $fixed
##
## short medium long
## 52 63 35
Выводы: Способ дискретизации влияет на результат,
универсально лучшего метода нет. Метод interval делит
диапазон значений на отрезки равной ширины и не учитывает, как данные
распределены внутри диапазона, поэтому группы получились неравными (52,
70 и 28 наблюдений). Метод frequency специально выравнивает
количество наблюдений в группах (46, 53 и 51), но границы ставит там,
где это нужно для равного счёта, а не там, где данные естественно
делятся. Метод cluster ищет границы по кластерам значений;
для Sepal.Length явных «провалов» в данных нет, поэтому его
результат близок к frequency (46, 53 и 51). Метод
fixed целиком зависит от границ, заданных вручную (здесь
5.5 и 6.5), и уместен там, где есть предметные пороги, например
возрастные группы. Это связано с частью 2: значения information gain
зависят от того, как нарезаны числовые признаки.
Подключим пакет Boruta
library(Boruta)
library(mlbench)
Загружаем данные
data("Ozone")
nrow(Ozone)
## [1] 366
Убираем строки, в которых есть пропуски
oz <-na.omit(Ozone)
Количество строк в новом множестве
nrow(oz)
## [1] 203
Запускаем алгоритм Boruta для определения важных свойств.
set.seed(123)
b <- Boruta(V4 ~ ., data = oz)
print(b)
## Boruta performed 18 iterations in 0.9217548 secs.
## 9 attributes confirmed important: V1, V10, V11, V12, V13 and 4 more;
## 3 attributes confirmed unimportant: V2, V3, V6;
строим график
plot(b, las = 2, xlab = "", cex.axis = 0.8) # boxplot важности
attStats(b)
## meanImp medianImp minImp maxImp normHits decision
## V1 0.40660229 0.39874710 0.37208192 0.50361585 1.00000000 Confirmed
## V2 0.04391890 0.03169822 -0.02604472 0.13137821 0.11111111 Rejected
## V3 -0.03235907 -0.04528594 -0.08039941 0.03178173 0.00000000 Rejected
## V5 0.37961385 0.37663333 0.33809868 0.42147363 1.00000000 Confirmed
## V6 0.04740489 0.05979747 -0.01572443 0.12718160 0.05555556 Rejected
## V7 0.51883205 0.50588169 0.46725744 0.59114026 1.00000000 Confirmed
## V8 0.77086858 0.76889989 0.73251688 0.81327051 1.00000000 Confirmed
## V9 0.86915268 0.86432645 0.84273725 0.92413075 1.00000000 Confirmed
## V10 0.42389058 0.41739463 0.38164901 0.49245967 1.00000000 Confirmed
## V11 0.55143690 0.54721914 0.49482770 0.61161304 1.00000000 Confirmed
## V12 0.66051279 0.65832222 0.60234621 0.73958611 1.00000000 Confirmed
## V13 0.42761454 0.42057856 0.38618408 0.46920282 1.00000000 Confirmed
Выводы: После удаления строк с пропусками для
анализа осталось 203 наблюдений из 366. Алгоритм Boruta завершил работу
за 18 итераций, не дойдя до лимита, так как по всем признакам удалось
принять решение: 9 признаков подтверждены как важные для предсказания
максимума концентрации озона (V1, V5, V7, V8, V9, V10, V11, V12, V13), 3
признака отвергнуты (V2, V3, V6), сомнительных признаков нет. Отвергнуты
число месяца (V2), день недели (V3) и скорость ветра (V6). Первые два
признака не имеют причинной связи с озоном, что согласуется со здравым
смыслом. Среди подтверждённых есть месяц (V1, сезонность) и
температурные показатели (V8, V9, V12): озон образуется на солнце и в
тепле. Ограничение результата в том, что из-за пропусков пришлось убрать
почти половину строк, а метод случайный, поэтому при другом
set.seed статусы пограничных признаков могли бы
измениться.