Дисциплина: Анализ данных Московский политехнический университет, факультет информационных технологий, кафедра «Информатика и информационные технологии».

Цель работы: освоить методы отбора и оценки важности признаков в R с помощью пакетов caret, FSelector, arules и Boruta, выполнить графический разведочный анализ и дискретизацию непрерывных переменных, оформить результаты в R Markdown и опубликовать отчёт на RPubs.

library(caret)      # обучение моделей, featurePlot
library(FSelector)  # оценка важности признаков (фильтры)
library(arules)     # discretize: дискретизация переменных
library(Boruta)     # отбор признаков методом Boruta
library(mlbench)    # набор данных Ozone

set.seed(42)        # фиксируем генератор случайных чисел для воспроизводимости

1 Пакет caret: список методов и графический разведочный анализ

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

Функция getModelInfo() возвращает описание всех моделей, которые можно обучать функцией train(). Названия методов получаем через names().

models <- names(getModelInfo())
length(models)   # число доступных методов
## [1] 239
head(models, 40) # первые 40 названий
##  [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"

В установленной версии caret доступно 239 методов. Часть из них содержит встроенный отбор признаков (embedded-методы), например:

# Методы со встроенным отбором признаков (регуляризация, деревья, ступенчатая регрессия)
fs_like <- c("glmnet", "lasso", "relaxo", "rpart", "rf", "Boruta", "glmStepAIC", "lmStepAIC",
             "leapForward", "leapBackward", "leapSeq", "penalized", "rfRules", "xgbTree")
fs_like[fs_like %in% models]
##  [1] "glmnet"       "lasso"        "relaxo"       "rpart"        "rf"          
##  [6] "glmStepAIC"   "lmStepAIC"    "leapForward"  "leapBackward" "leapSeq"     
## [11] "penalized"    "rfRules"      "xgbTree"

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

1.2 Данные из справочного файла

Используем пример из справки к featurePlot: матрица из 50 наблюдений и 5 случайных нормально распределённых признаков и бинарный класс A/B.

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

dim(x)
## [1] 50  5
table(y)
## y
##  A  B 
## 25 25
summary(x)
##        X1                 X2                X3                X4          
##  Min.   :-2.65646   Min.   :-2.9931   Min.   :-2.0247   Min.   :-1.74506  
##  1st Qu.:-0.63223   1st Qu.:-0.5262   1st Qu.:-0.5884   1st Qu.:-0.57145  
##  Median :-0.10039   Median : 0.2672   Median :-0.2523   Median : 0.03927  
##  Mean   :-0.03567   Mean   : 0.1007   Mean   :-0.1513   Mean   :-0.02372  
##  3rd Qu.: 0.65072   3rd Qu.: 0.6728   3rd Qu.: 0.1341   3rd Qu.: 0.57841  
##  Max.   : 2.28665   Max.   : 1.5757   Max.   : 2.7019   Max.   : 1.81523  
##        X5          
##  Min.   :-2.00093  
##  1st Qu.:-0.69399  
##  Median :-0.15305  
##  Mean   : 0.00794  
##  3rd Qu.: 0.65662  
##  Max.   : 2.05954

1.3 Графики featurePlot()

Функция featurePlot() строит решётчатые (lattice) графики признаков с разбивкой по классу. Каждый график сохраняем в файл *.jpg функцией jpeg() и выводим в отчёт.

# Вспомогательная функция: построить график, сохранить в JPG и показать в отчёте
save_feature_plot <- function(plot_type, file_name, ...) {
  p <- featurePlot(x = x, y = y, plot = plot_type, ...)
  jpeg(file_name, width = 1000, height = 800, quality = 95)
  print(p)
  dev.off()
  p
}

types <- c("pairs", "ellipse", "density", "box", "strip")
files <- paste0("featurePlot_", types, ".jpg")

1.3.1 Матрица диаграмм рассеяния (pairs)

save_feature_plot("pairs", files[1], auto.key = list(columns = 2))

1.3.2 Диаграммы рассеяния с эллипсами (ellipse)

save_feature_plot("ellipse", files[2], auto.key = list(columns = 2))

1.3.3 Плотности распределения (density)

save_feature_plot("density", files[3],
                  scales = list(x = list(relation = "free"), y = list(relation = "free")),
                  adjust = 1.5, pch = "|", layout = c(5, 1), auto.key = list(columns = 2))

1.3.4 Ящики с усами (box)

save_feature_plot("box", files[4],
                  scales = list(y = list(relation = "free")), layout = c(5, 1))

1.3.5 Точечные диаграммы (strip)

save_feature_plot("strip", files[5], jitter = TRUE, layout = c(5, 1))

# Проверяем, что файлы сохранены
data.frame(file = files, exists = file.exists(files), size_kb = round(file.size(files) / 1024, 1))
##                      file exists size_kb
## 1   featurePlot_pairs.jpg   TRUE   153.8
## 2 featurePlot_ellipse.jpg   TRUE   228.6
## 3 featurePlot_density.jpg   TRUE   100.8
## 4     featurePlot_box.jpg   TRUE    86.9
## 5   featurePlot_strip.jpg   TRUE    68.3
# Количественная проверка: различаются ли средние признаков в классах A и B
pv <- sapply(as.data.frame(x), function(col) t.test(col ~ y)$p.value)
round(pv, 3)
##    X1    X2    X3    X4    X5 
## 0.563 0.511 0.293 0.968 0.277

Выводы. Признаки сгенерированы случайно и не связаны с классом, что подтверждают все графики: облака точек классов A и B на матрице рассеяния и эллипсы рассеяния полностью перекрываются, кривые плотности и ящики с усами для двух классов почти совпадают. Корреляции между признаками нет — эллипсы близки к окружностям. P-значения t-теста (0.56; 0.51; 0.29; 0.97; 0.28) в 5 случаях из 5 больше 0,05, то есть различия средних между классами статистически незначимы. Данные случайные, поэтому при другом запуске значения изменятся; отдельное p-значение меньше 0,05 при пяти проверках может появиться просто случайно (ожидаемая доля ложных срабатываний — 5 %). Следовательно, ни один из признаков не обладает разделяющей способностью, и строить по ним классификатор бессмысленно. Пример показывает, как featurePlot() позволяет до обучения модели быстро оценить, есть ли у признаков связь с целевой переменной.

2 Пакет FSelector: важность признаков для классификации ирисов

Набор iris содержит 150 цветков трёх видов и 4 числовых признака. Оценим важность признаков для предсказания вида (Species) разными методами-фильтрами.

data(iris)
str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
# Методы-фильтры FSelector (чем больше вес, тем важнее признак)
w_chi  <- chi.squared(Species ~ ., iris)               # статистика хи-квадрат (коэффициент Крамера)
w_ig   <- information.gain(Species ~ ., iris)          # прирост информации
w_gr   <- gain.ratio(Species ~ ., iris)                # нормированный прирост информации
w_su   <- symmetrical.uncertainty(Species ~ ., iris)   # симметричная неопределённость
w_oner <- oneR(Species ~ ., iris)                      # точность правила одного признака
w_rel  <- relief(Species ~ ., iris, neighbours.count = 5, sample.size = 20)  # алгоритм RReliefF
w_rf   <- random.forest.importance(Species ~ ., iris, importance.type = 1)    # случайный лес

weights <- data.frame(
  chi.squared  = w_chi$attr_importance,
  info.gain    = w_ig$attr_importance,
  gain.ratio   = w_gr$attr_importance,
  sym.uncert   = w_su$attr_importance,
  oneR         = w_oner$attr_importance,
  relief       = w_rel$attr_importance,
  random.forest = w_rf$attr_importance,
  row.names = rownames(w_chi)
)
round(weights, 3)
##              chi.squared info.gain gain.ratio sym.uncert  oneR relief
## Sepal.Length       0.629     0.452      0.420      0.416 0.173  0.173
## Sepal.Width        0.492     0.267      0.247      0.245 0.040  0.154
## Petal.Length       0.935     0.940      0.858      0.857 0.400  0.359
## Petal.Width        0.943     0.955      0.871      0.871 0.407  0.384
##              random.forest
## Sepal.Length        15.590
## Sepal.Width          6.485
## Petal.Length        47.429
## Petal.Width         47.220
# Нормируем веса каждого метода к [0, 1], чтобы сравнить их на одном графике
norm_w <- apply(weights, 2, function(v) (v - min(v)) / (max(v) - min(v)))
barplot(t(norm_w), beside = TRUE, legend.text = colnames(norm_w),
        args.legend = list(x = "topleft", cex = 0.75, bty = "n"),
        col = hcl.colors(ncol(norm_w), "Set 2"),
        main = "Нормированная важность признаков iris по методам FSelector",
        ylab = "Важность (0–1)")

# Отбор подмножеств признаков
cutoff.k(w_ig, 2)                       # два лучших признака по приросту информации
## [1] "Petal.Width"  "Petal.Length"
cutoff.k.percent(w_chi, 0.5)            # лучшие 50 % по хи-квадрат
## [1] "Petal.Width"  "Petal.Length"
cfs(Species ~ ., iris)                  # CFS: корреляционный отбор подмножества
## [1] "Petal.Length" "Petal.Width"
consistency(Species ~ ., iris)          # отбор по критерию согласованности
## [1] "Petal.Width"

Выводы. Все семь методов единодушно ставят на первые места признаки лепестка — Petal.Width и Petal.Length: их веса близки к максимальным (прирост информации около 0,94–0,96 против 0,27–0,45 у чашелистика). Наименее информативен Sepal.Width. Функция cutoff.k выбирает два лучших признака — Petal.Width и Petal.Length, а корреляционный отбор CFS оставляет подмножество Petal.Length, Petal.Width. Это согласуется с известной структурой набора: вид setosa полностью отделяется по размерам лепестка, а versicolor и virginica различаются в основном по ним же. Для классификации ирисов достаточно признаков лепестка, а признаки чашелистика можно исключить почти без потери качества.

3 Пакет arules: дискретизация непрерывной переменной

Функция discretize() преобразует непрерывную переменную в категориальную (фактор). Применим её к длине лепестка iris$Petal.Length четырьмя методами, разбивая на 3 категории.

pl <- iris$Petal.Length
summary(pl)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   1.600   4.350   3.758   5.100   6.900
d_interval  <- discretize(pl, method = "interval",  breaks = 3)                 # равная ширина
d_frequency <- discretize(pl, method = "frequency", breaks = 3)                 # равная частота
d_cluster   <- discretize(pl, method = "cluster",   breaks = 3)                 # k-средних
d_fixed     <- discretize(pl, method = "fixed",     breaks = c(-Inf, 2.5, 4.8, 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     45     55
# Гистограммы с границами интервалов каждого метода
par(mfrow = c(2, 2))
draw_cuts <- function(d, title) {
  hist(pl, breaks = 30, col = "grey85", border = "white", main = title, xlab = "Petal.Length")
  abline(v = attr(d, "discretized:breaks"), col = "red", lwd = 2, lty = 2)
}
draw_cuts(d_interval,  "interval (равная ширина)")
draw_cuts(d_frequency, "frequency (равная частота)")
draw_cuts(d_cluster,   "cluster (k-средних)")
draw_cuts(d_fixed,     "fixed (заданные границы)")

par(mfrow = c(1, 1))
# Насколько полученные категории соответствуют видам ирисов
table(interval  = d_interval,  Species = iris$Species)
##              Species
## 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(frequency = d_frequency, Species = iris$Species)
##             Species
## 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(cluster   = d_cluster,   Species = iris$Species)
##              Species
## 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(fixed     = d_fixed,     Species = iris$Species)
##         Species
## fixed    setosa versicolor virginica
##   short      50          0         0
##   medium      0         44         1
##   long        0          6        49
# Дискретизация всех числовых признаков набора сразу (метод frequency)
iris_disc <- discretizeDF(iris, default = list(method = "frequency", breaks = 3))
head(iris_disc)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1    [4.3,5.4)   [3.2,4.4]     [1,2.63) [0.1,0.867)  setosa
## 2    [4.3,5.4)   [2.9,3.2)     [1,2.63) [0.1,0.867)  setosa
## 3    [4.3,5.4)   [3.2,4.4]     [1,2.63) [0.1,0.867)  setosa
## 4    [4.3,5.4)   [2.9,3.2)     [1,2.63) [0.1,0.867)  setosa
## 5    [4.3,5.4)   [3.2,4.4]     [1,2.63) [0.1,0.867)  setosa
## 6    [5.4,6.3)   [3.2,4.4]     [1,2.63) [0.1,0.867)  setosa
# Сколько цветков попадает в категорию, где преобладает другой вид (чем меньше — тем лучше категории совпадают с видами)
mismatch <- function(d) { t <- table(d, iris$Species); sum(t) - sum(apply(t, 1, max)) }
quality <- data.frame(method = c("interval", "frequency", "cluster", "fixed"),
                      breaks = sapply(list(d_interval, d_frequency, d_cluster, d_fixed),
                                      function(d) paste(round(attr(d, "discretized:breaks"), 2), collapse = " | ")),
                      sizes = sapply(list(d_interval, d_frequency, d_cluster, d_fixed),
                                     function(d) paste(table(d), collapse = "/")),
                      mismatches = sapply(list(d_interval, d_frequency, d_cluster, d_fixed), mismatch))
quality
##      method                 breaks    sizes mismatches
## 1  interval  1 | 2.97 | 4.93 | 6.9 50/54/46          8
## 2 frequency   1 | 2.63 | 4.9 | 6.9 50/49/51          7
## 3   cluster  1 | 2.85 | 4.89 | 6.9 50/49/51          7
## 4     fixed -Inf | 2.5 | 4.8 | Inf 50/45/55          7

Выводы.

  • interval делит диапазон [1; 6,9] на три интервала одинаковой ширины; численность категорий получается неравной (50/54/46), так как распределение длины лепестка двухвершинное.
  • frequency строит границы по квантилям, поэтому категории содержат почти одинаковое число наблюдений (50/49/51).
  • cluster подбирает границы методом k-средних по «провалам» распределения (1 | 2.85 | 4.89 | 6.9), то есть по естественной группировке данных.
  • fixed использует границы, заданные исследователем по знанию предметной области (2,5 и 4,8 см); категория short точно совпадает с видом setosa.

Все методы безошибочно выделяют setosa, у которой лепестки заметно короче. Виды versicolor и virginica по длине лепестка частично перекрываются, поэтому при любом методе 7–8 цветков попадают в «чужую» категорию. Методы различаются положением границ и численностью категорий: frequency удобен, когда нужны сбалансированные категории, cluster — когда границы должны отражать естественные группы, а fixed — когда границы заданы предметной областью и категории должны иметь понятные названия.

4 Пакет Boruta: отбор признаков для набора Ozone

Набор Ozone (пакет mlbench) содержит ежедневные измерения в районе Лос-Анджелеса за 1976 год. Целевая переменная V4 — среднесуточная максимальная концентрация озона; признаки — календарные (V1 месяц, V2 день месяца, V3 день недели) и метеорологические (V5–V13: давление, ветер, влажность, температура, высота инверсии, видимость и др.).

data("Ozone", package = "mlbench")
dim(Ozone)
## [1] 366  13
colSums(is.na(Ozone))          # пропуски по столбцам
##  V1  V2  V3  V4  V5  V6  V7  V8  V9 V10 V11 V12 V13 
##   0   0   0   5  12   0  15   2 139  15   1  14   0
oz <- na.omit(Ozone)           # Boruta не работает с пропусками — удаляем неполные строки
dim(oz)
## [1] 203  13

Алгоритм Boruta добавляет к данным «теневые» признаки (перемешанные копии исходных), многократно обучает случайный лес и сравнивает важность каждого реального признака с наибольшей важностью теневых. Признак подтверждается (Confirmed), если стабильно важнее теневых, и отвергается (Rejected), если стабильно хуже.

set.seed(42)
boruta_oz <- Boruta(V4 ~ ., data = oz, doTrace = 0, maxRuns = 200)
print(boruta_oz)
## Boruta performed 21 iterations in 1.497149 secs.
##  9 attributes confirmed important: V1, V10, V11, V12, V13 and 4 more;
##  3 attributes confirmed unimportant: V2, V3, V6;
# Boxplot важности признаков: зелёные — подтверждённые, жёлтые — неопределённые,
# красные — отвергнутые, синие — теневые признаки (min, mean, max)
plot(boruta_oz, las = 2, cex.axis = 0.8, xlab = "",
     main = "Boruta: важность признаков для предсказания V4 (Ozone)")

# История изменения важности по итерациям
plotImpHistory(boruta_oz, main = "Boruta: история важности по итерациям")

final_oz <- TentativeRoughFix(boruta_oz)   # окончательное решение по неопределённым признакам
stats <- attStats(final_oz)
stats[order(-stats$meanImp), ]
##         meanImp   medianImp      minImp     maxImp  normHits  decision
## V9   0.88789089  0.88686438  0.82668819 0.99602188 1.0000000 Confirmed
## V8   0.76671886  0.76754267  0.70941659 0.83272127 1.0000000 Confirmed
## V12  0.65481180  0.66243636  0.58700258 0.69820162 1.0000000 Confirmed
## V11  0.54787436  0.54192563  0.50000214 0.62146819 1.0000000 Confirmed
## V7   0.52271613  0.51879561  0.47630171 0.56514800 1.0000000 Confirmed
## V10  0.43906302  0.43962045  0.36477347 0.51909248 1.0000000 Confirmed
## V13  0.41297531  0.41782410  0.35739084 0.45129002 1.0000000 Confirmed
## V1   0.40832700  0.41286596  0.32226876 0.47443026 1.0000000 Confirmed
## V5   0.39237797  0.39139823  0.31671328 0.44197223 1.0000000 Confirmed
## V6   0.05371338  0.06176711 -0.04041718 0.12062160 0.1428571  Rejected
## V2   0.03630079  0.04288927 -0.03301760 0.10983565 0.0000000  Rejected
## V3  -0.04668555 -0.05121816 -0.12973322 0.05836646 0.0000000  Rejected
getSelectedAttributes(final_oz, withTentative = FALSE)
## [1] "V1"  "V5"  "V7"  "V8"  "V9"  "V10" "V11" "V12" "V13"
# Проверка: сравним качество случайного леса на всех признаках и на отобранных
library(randomForest)
selected <- getSelectedAttributes(final_oz)
set.seed(42); rf_all <- randomForest(V4 ~ ., data = oz)
set.seed(42); rf_sel <- randomForest(reformulate(selected, "V4"), data = oz)
data.frame(model = c("Все признаки", "Отобранные Boruta"),
           n_features = c(ncol(oz) - 1, length(selected)),
           var_explained = round(c(tail(rf_all$rsq, 1), tail(rf_sel$rsq, 1)) * 100, 1))
##               model n_features var_explained
## 1      Все признаки         12          67.8
## 2 Отобранные Boruta          9          74.7

Выводы. Boruta подтвердил важность признаков V1, V5, V7, V8, V9, V10, V11, V12, V13 и отверг V2, V3, V6. На boxplot подтверждённые признаки (зелёные) расположены значительно выше максимального теневого признака (shadowMax), а отвергнутые (красные) находятся на уровне теневых, то есть их важность не отличается от важности случайного шума.

Наибольшую важность имеют V9, V8, V12 — метеорологические показатели: температуры в двух точках наблюдения (V8, V9), температура и высота инверсионного слоя (V12, V10), влажность (V7), градиент давления (V11). Важен и месяц (V1) — он отражает сезонность образования озона. День месяца (V2) и день недели (V3) отвергнуты, что физически объяснимо: концентрация озона определяется погодой и сезоном, а не календарной датой.

Проверка случайным лесом показала, что модель на 9 отобранных признаках объясняет 74.7 % дисперсии против 67.8 % на всех 12 признаках. Удаление отвергнутых признаков не ухудшает модель, а избавляет её от шума.

5 Общие выводы

  1. Пакет caret предоставляет единый интерфейс к сотням методов машинного обучения и инструменты разведочного анализа; featurePlot() позволяет визуально оценить разделяющую способность признаков ещё до обучения модели. На случайных данных из справки признаки не разделяют классы.
  2. Методы-фильтры FSelector на наборе iris согласованно выделяют признаки лепестка (Petal.Width, Petal.Length) как наиболее важные для классификации видов.
  3. Функция discretize() пакета arules реализует четыре стратегии дискретизации; результат существенно зависит от метода, и для данных с групповой структурой лучше подходят кластерный метод или заданные экспертом границы.
  4. Метод-обёртка Boruta на наборе Ozone отбирает метеорологические и сезонный признаки и отвергает нерелевантные (день месяца, день недели), не ухудшая качество модели.