# Один раз в консоли (не при каждом Knit):
# install.packages(c("caret", "ellipse", "FSelector", "arules", "Boruta", "mlbench"))
library(caret)
library(FSelector)
library(arules)
library(Boruta)
library(mlbench)

1. Пакет caret и разведочный анализ с featurePlot()

Список доступных методов

models <- names(getModelInfo())
length(models)   # сколько методов доступно
## [1] 239
head(models, 50) # первые 50 названий
##  [1] "ada"            "AdaBag"         "AdaBoost.M1"    "adaboost"      
##  [5] "amdai"          "ANFIS"          "avNNet"         "awnb"          
##  [9] "awtan"          "bag"            "bagEarth"       "bagEarthGCV"   
## [13] "bagFDA"         "bagFDAGCV"      "bam"            "bartMachine"   
## [17] "bayesglm"       "binda"          "blackboost"     "blasso"        
## [21] "blassoAveraged" "bridge"         "brnn"           "BstLm"         
## [25] "bstSm"          "bstTree"        "C5.0"           "C5.0Cost"      
## [29] "C5.0Rules"      "C5.0Tree"       "cforest"        "chaid"         
## [33] "CSimca"         "ctree"          "ctree2"         "cubist"        
## [37] "dda"            "deepboost"      "DENFIS"         "dnn"           
## [41] "dwdLinear"      "dwdPoly"        "dwdRadial"      "earth"         
## [45] "elm"            "enet"           "evtree"         "extraTrees"    
## [49] "fda"            "FH.GBML"

Функция getModelInfo() возвращает описания всех моделей, которые поддерживает caret. Среди них есть модели со встроенным отбором признаков (например, glmnet, lasso, rf, rpart). Кроме того, в caret есть отдельные обёртки для отбора признаков: rfe() (рекурсивное исключение признаков), sbf() (фильтрация по одномерным критериям) и gafs() / safs() (генетический алгоритм и имитация отжига).

Данные из справки caret

set.seed(123)
x <- matrix(rnorm(50 * 5), ncol = 5)
y <- factor(rep(c("A", "B"), 25))
colnames(x) <- paste0("X", 1:5)

Графики featurePlot() с сохранением в jpg

plots <- list(
  pairs   = featurePlot(x, y, plot = "pairs",   auto.key = list(columns = 2)),
  ellipse = featurePlot(x, y, plot = "ellipse", auto.key = list(columns = 2)),
  box     = featurePlot(x, y, plot = "box",
                        scales = list(y = list(relation = "free")),
                        layout = c(5, 1)),
  density = featurePlot(x, y, plot = "density",
                        scales = list(x = list(relation = "free"),
                                      y = list(relation = "free")),
                        auto.key = list(columns = 2), layout = c(5, 1)),
  strip   = featurePlot(x, y, plot = "strip", jitter = TRUE,
                        layout = c(5, 1))
)

for (nm in names(plots)) {
  jpeg(paste0("featurePlot_", nm, ".jpg"), width = 1000, height = 700)
  print(plots[[nm]])
  dev.off()
  print(plots[[nm]])   # показать график и в отчёте
}

Выводы. Признаки X1–X5 сгенерированы функцией rnorm() независимо от метки класса, поэтому на всех графиках классы A и B не разделяются: точки на матрице рассеяния перемешаны, эллипсы почти совпадают, медианы и размах на боксплотах близки, кривые плотности накладываются друг на друга. Ни один признак не несёт информации о классе, и для такой задачи классификации отбор признаков ничего не даст. Графики featurePlot() позволяют быстро увидеть это ещё до построения модели.

2. Важность признаков с помощью FSelector (набор iris)

data(iris)

ig  <- information.gain(Species ~ ., iris)
gr  <- gain.ratio(Species ~ ., iris)
su  <- symmetrical.uncertainty(Species ~ ., iris)
chi <- chi.squared(Species ~ ., iris)
set.seed(123)
rfi <- random.forest.importance(Species ~ ., iris, importance.type = 1)

importance <- data.frame(
  information.gain        = ig$attr_importance,
  gain.ratio              = gr$attr_importance,
  symmetrical.uncertainty = su$attr_importance,
  chi.squared             = chi$attr_importance,
  random.forest           = rfi$attr_importance,
  row.names = rownames(ig)
)
round(importance, 3)
##              information.gain gain.ratio symmetrical.uncertainty chi.squared
## Sepal.Length            0.452      0.420                   0.416       0.629
## Sepal.Width             0.267      0.247                   0.245       0.492
## Petal.Length            0.940      0.858                   0.857       0.935
## Petal.Width             0.955      0.871                   0.871       0.943
##              random.forest
## Sepal.Length        15.233
## Sepal.Width          6.772
## Petal.Length        47.478
## Petal.Width         46.266
# Два лучших признака по информационному выигрышу
cutoff.k(ig, 2)
## [1] "Petal.Width"  "Petal.Length"
# Итоговая формула для модели
as.simple.formula(cutoff.k(ig, 2), "Species")
## Species ~ Petal.Width + Petal.Length
## <environment: 0x0000020cc60a4660>
barplot(sort(setNames(ig$attr_importance, rownames(ig))),
        horiz = TRUE, las = 1, col = "steelblue",
        main = "Information gain, iris")

Выводы. Все методы дают один и тот же порядок: самые важные признаки — Petal.Width и Petal.Length, затем Sepal.Length, наименее информативен Sepal.Width. Для классификации ирисов достаточно двух признаков лепестка, и признаки чашелистика можно исключить почти без потери качества.

3. Дискретизация непрерывной переменной (arules::discretize)

Преобразуем Petal.Length в категориальную переменную с тремя уровнями.

pl <- iris$Petal.Length

d_interval  <- discretize(pl, method = "interval",  breaks = 3)
d_frequency <- discretize(pl, method = "frequency", breaks = 3)
set.seed(123)
d_cluster   <- discretize(pl, method = "cluster",   breaks = 3)
d_fixed     <- discretize(pl, method = "fixed",
                          breaks = c(-Inf, 2, 5, Inf),
                          labels = c("short", "medium", "long"))

table(d_interval)
## d_interval
##    [1,2.97) [2.97,4.93)  [4.93,6.9] 
##          50          54          46
table(d_frequency)
## d_frequency
##   [1,2.63) [2.63,4.9)  [4.9,6.9] 
##         50         49         51
table(d_cluster)
## d_cluster
##    [1,2.85) [2.85,4.89)  [4.89,6.9] 
##          50          49          51
table(d_fixed)
## d_fixed
##  short medium   long 
##     50     54     46

Насколько каждая дискретизация совпадает с видом цветка:

table(d_interval,  iris$Species)
##              
## d_interval    setosa versicolor virginica
##   [1,2.97)        50          0         0
##   [2.97,4.93)      0         48         6
##   [4.93,6.9]       0          2        44
table(d_frequency, iris$Species)
##             
## d_frequency  setosa versicolor virginica
##   [1,2.63)       50          0         0
##   [2.63,4.9)      0         46         3
##   [4.9,6.9]       0          4        47
table(d_cluster,   iris$Species)
##              
## d_cluster     setosa versicolor virginica
##   [1,2.85)        50          0         0
##   [2.85,4.89)      0         46         3
##   [4.89,6.9]       0          4        47
table(d_fixed,     iris$Species)
##         
## d_fixed  setosa versicolor virginica
##   short      50          0         0
##   medium      0         48         6
##   long        0          2        44
par(mfrow = c(2, 2))
hist(pl, breaks = 20, main = "interval",  xlab = "Petal.Length")
abline(v = attr(d_interval,  "discretized:breaks"), col = "red",  lwd = 2)
hist(pl, breaks = 20, main = "frequency", xlab = "Petal.Length")
abline(v = attr(d_frequency, "discretized:breaks"), col = "blue", lwd = 2)
hist(pl, breaks = 20, main = "cluster",   xlab = "Petal.Length")
abline(v = attr(d_cluster,   "discretized:breaks"), col = "darkgreen", lwd = 2)
hist(pl, breaks = 20, main = "fixed",     xlab = "Petal.Length")
abline(v = c(2, 5), col = "purple", lwd = 2)

par(mfrow = c(1, 1))

Выводы.

  • interval (равная ширина, ≈1,97 см) дал группы 50 / 54 / 46; первая группа (1–2,97) точно совпадает с 50 цветками setosa.
  • frequency дал почти равные группы 50 / 49 / 51, как и задумано методом, но границы (2,63 и 4,9) выбраны по числу наблюдений, а не по структуре данных.
  • cluster (k-means) дал группы 50 / 66 / 34: setosa отделена точно, но вторая граница стоит на 5,13, и средняя группа получилась заметно больше крайней.
  • fixed с границами 2 и 5 дал 50 / 54 / 46 (short / medium / long), что близко к реальному делению на три вида по 50 цветков.

Все методы уверенно отделяют setosa (короткие лепестки до 2 см), а граница между versicolor и virginica размыта, и методы проводят её по-разному. Выбор метода зависит от задачи: frequency — для равных по размеру групп, cluster — для поиска естественных групп, fixed — когда границы известны из предметной области.

4. Отбор признаков с Boruta (набор Ozone)

data("Ozone", package = "mlbench")
str(Ozone)
## 'data.frame':    366 obs. of  13 variables:
##  $ V1 : Factor w/ 12 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ V2 : Factor w/ 31 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ V3 : Factor w/ 7 levels "1","2","3","4",..: 4 5 6 7 1 2 3 4 5 6 ...
##  $ V4 : num  3 3 3 5 5 6 4 4 6 7 ...
##  $ V5 : num  5480 5660 5710 5700 5760 5720 5790 5790 5700 5700 ...
##  $ V6 : num  8 6 4 3 3 4 6 3 3 3 ...
##  $ V7 : num  20 NA 28 37 51 69 19 25 73 59 ...
##  $ V8 : num  NA 38 40 45 54 35 45 55 41 44 ...
##  $ V9 : num  NA NA NA NA 45.3 ...
##  $ V10: num  5000 NA 2693 590 1450 ...
##  $ V11: num  -15 -14 -25 -24 25 15 -33 -28 23 -2 ...
##  $ V12: num  30.6 NA 47.7 55 57 ...
##  $ V13: num  200 300 250 100 60 60 100 250 120 120 ...
# V4 — суточный максимум озона (целевая переменная). Boruta не работает с NA.
ozone <- na.omit(Ozone)

set.seed(123)
boruta_res <- Boruta(V4 ~ ., data = ozone, doTrace = 0)
print(boruta_res)
## Boruta performed 18 iterations in 0.4642251 secs.
##  9 attributes confirmed important: V1, V10, V11, V12, V13 and 4 more;
##  3 attributes confirmed unimportant: V2, V3, V6;
final <- TentativeRoughFix(boruta_res)
getSelectedAttributes(final)
## [1] "V1"  "V5"  "V7"  "V8"  "V9"  "V10" "V11" "V12" "V13"
attStats(final)
##         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
plot(boruta_res, las = 2, cex.axis = 0.8, xlab = "",
     main = "Важность признаков Boruta, набор Ozone")

Выводы. На графике зелёные боксплоты — подтверждённо важные признаки, красные — отклонённые, синие — «теневые» признаки (перемешанные копии, эталон случайного шума). Признак считается важным, если его важность устойчиво выше лучшего теневого (shadowMax).

Boruta завершилась за 18 итераций и подтвердила важность 9 признаков: V1 (месяц), V5 (высота изобары 500 мб), V7 (влажность), V8 и V9 (температуры), V10 (высота основания инверсии), V11 (градиент давления), V12 (температура основания инверсии), V13 (видимость). Самые важные — температуры V9 и V8, затем V12. Отклонены 3 признака: V2 (день месяца), V3 (день недели) и V6 (скорость ветра): их важность не выше, чем у случайного шума. Уровень озона определяется метеоусловиями и сезоном (месяц важен), а день месяца и день недели на него не влияют.

Общие выводы

В работе рассмотрены три подхода к подготовке признаков: визуальный разведочный анализ (caret::featurePlot), фильтрационные методы оценки важности (FSelector) и обёрточный метод на основе случайного леса (Boruta), а также дискретизация непрерывных переменных (arules::discretize). Фильтры работают быстро и дают понятный рейтинг признаков, а Boruta сравнивает признаки со случайным шумом и отбирает все признаки, действительно связанные с целевой переменной.