1 Цель работы

Изучить графический разведочный анализ данных, оценку информативности признаков, дискретизацию непрерывных переменных и автоматический выбор признаков. Используются пакеты caret, FSelector, arules и Boruta.

В работе четыре вычислительных задания и инструкция публикации отчёта.

2 Подготовка среды

Перед первым нажатием Knit выполнить в консоли RStudio следующий код. Этот блок показан в отчёте, но установка пакетов при сборке отчёта не запускается. Рекомендуется актуальная версия R; для новых версий arules может потребоваться обновление R. FSelector использует Java через свои зависимости: на Windows нужно установить JDK той же разрядности, что и R, затем перезапустить RStudio.

packages <- c("caret", "FSelector", "arules", "Boruta", "mlbench",
              "ranger", "randomForest", "knitr", "rmarkdown")
new_packages <- packages[!vapply(packages, requireNamespace,
                               logical(1), quietly = TRUE)]
if (length(new_packages)) {
  install.packages(new_packages, repos = "https://cloud.r-project.org")
}

Фиксированные значения set.seed() позволяют повторить генерацию данных и случайные процедуры в одинаковой программной среде. Версии пакетов приведены в конце отчёта; при изменении версий отдельные результаты могут отличаться. Каждый рисунок сохраняется в папку graphs в формате JPG и включается в HTML.

3 CARET: модели и графический разведочный анализ

3.1 Список доступных моделей и средства выбора признаков

library(caret)
model_names <- names(getModelInfo())
print(model_names)
##   [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"             "FIR.DM"              "foba"               
##  [53] "FRBCS.CHI"           "FRBCS.W"             "FS.HGD"              "gam"                
##  [57] "gamboost"            "gamLoess"            "gamSpline"           "gaussprLinear"      
##  [61] "gaussprPoly"         "gaussprRadial"       "gbm_h2o"             "gbm"                
##  [65] "gcvEarth"            "GFS.FR.MOGUL"        "GFS.LT.RS"           "GFS.THRIFT"         
##  [69] "glm.nb"              "glm"                 "glmboost"            "glmnet_h2o"         
##  [73] "glmnet"              "glmStepAIC"          "gpls"                "hda"                
##  [77] "hdda"                "hdrda"               "HYFIS"               "icr"                
##  [81] "J48"                 "JRip"                "kernelpls"           "kknn"               
##  [85] "knn"                 "krlsPoly"            "krlsRadial"          "lars"               
##  [89] "lars2"               "lasso"               "lda"                 "lda2"               
##  [93] "leapBackward"        "leapForward"         "leapSeq"             "Linda"              
##  [97] "lm"                  "lmStepAIC"           "LMT"                 "loclda"             
## [101] "logicBag"            "LogitBoost"          "logreg"              "lssvmLinear"        
## [105] "lssvmPoly"           "lssvmRadial"         "lvq"                 "M5"                 
## [109] "M5Rules"             "manb"                "mda"                 "Mlda"               
## [113] "mlp"                 "mlpKerasDecay"       "mlpKerasDecayCost"   "mlpKerasDropout"    
## [117] "mlpKerasDropoutCost" "mlpML"               "mlpSGD"              "mlpWeightDecay"     
## [121] "mlpWeightDecayML"    "monmlp"              "msaenet"             "multinom"           
## [125] "mxnet"               "mxnetAdam"           "naive_bayes"         "nb"                 
## [129] "nbDiscrete"          "nbSearch"            "neuralnet"           "nnet"               
## [133] "nnls"                "nodeHarvest"         "null"                "OneR"               
## [137] "ordinalNet"          "ordinalRF"           "ORFlog"              "ORFpls"             
## [141] "ORFridge"            "ORFsvm"              "ownn"                "pam"                
## [145] "parRF"               "PART"                "partDSA"             "pcaNNet"            
## [149] "pcr"                 "pda"                 "pda2"                "penalized"          
## [153] "PenalizedLDA"        "plr"                 "pls"                 "plsRglm"            
## [157] "polr"                "ppr"                 "pre"                 "PRIM"               
## [161] "protoclass"          "qda"                 "QdaCov"              "qrf"                
## [165] "qrnn"                "randomGLM"           "ranger"              "rbf"                
## [169] "rbfDDA"              "Rborist"             "rda"                 "regLogistic"        
## [173] "relaxo"              "rf"                  "rFerns"              "RFlda"              
## [177] "rfRules"             "ridge"               "rlda"                "rlm"                
## [181] "rmda"                "rocc"                "rotationForest"      "rotationForestCp"   
## [185] "rpart"               "rpart1SE"            "rpart2"              "rpartCost"          
## [189] "rpartScore"          "rqlasso"             "rqnc"                "RRF"                
## [193] "RRFglobal"           "rrlda"               "RSimca"              "rvmLinear"          
## [197] "rvmPoly"             "rvmRadial"           "SBC"                 "sda"                
## [201] "sdwd"                "simpls"              "SLAVE"               "slda"               
## [205] "smda"                "snn"                 "sparseLDA"           "spikeslab"          
## [209] "spls"                "stepLDA"             "stepQDA"             "superpc"            
## [213] "svmBoundrangeString" "svmExpoString"       "svmLinear"           "svmLinear2"         
## [217] "svmLinear3"          "svmLinearWeights"    "svmLinearWeights2"   "svmPoly"            
## [221] "svmRadial"           "svmRadialCost"       "svmRadialSigma"      "svmRadialWeights"   
## [225] "svmSpectrumString"   "tan"                 "tanSearch"           "treebag"            
## [229] "vbmpRadial"          "vglmAdjCat"          "vglmContRatio"       "vglmCumulative"     
## [233] "widekernelpls"       "WM"                  "wsrf"                "xgbDART"            
## [237] "xgbLinear"           "xgbTree"             "xyf"
writeLines(model_names, "results/caret_models.txt")
length(model_names)
## [1] 239

Команда вывела 239 идентификаторов моделей, доступных через train(). Это каталог моделей, а не самостоятельный перечень методов выбора признаков. Наличие модели в каталоге не означает, что установлены все её дополнительные зависимости.

Для выбора признаков в CARET предусмотрены:

Средство Назначение
rfe() Рекурсивное исключение признаков с оценкой качества подмножеств
sbf() Отбор признаков фильтром внутри процедуры ресэмплинга
gafs() Поиск подмножества признаков генетическим алгоритмом
safs() Поиск подмножества методом имитации отжига
varImp() Оценка важности; сама по себе не задаёт правило отбора

Дополнительно выделим модели, отмеченные в каталоге тегом встроенного отбора:

model_info <- getModelInfo()
embedded_names <- names(model_info)[vapply(model_info, function(m) {
  any(grepl("Feature Selection", m$tags, ignore.case = TRUE))
}, logical(1))]
print(embedded_names)
##  [1] "ada"              "AdaBag"           "AdaBoost.M1"      "adaboost"         "bagEarth"        
##  [6] "bagEarthGCV"      "bagFDA"           "bagFDAGCV"        "bartMachine"      "blasso"          
## [11] "BstLm"            "bstSm"            "C5.0"             "C5.0Cost"         "C5.0Rules"       
## [16] "C5.0Tree"         "cforest"          "chaid"            "ctree"            "ctree2"          
## [21] "cubist"           "deepboost"        "earth"            "enet"             "evtree"          
## [26] "extraTrees"       "fda"              "foba"             "gamboost"         "gbm_h2o"         
## [31] "gbm"              "gcvEarth"         "glmnet_h2o"       "glmnet"           "glmStepAIC"      
## [36] "J48"              "JRip"             "lars"             "lars2"            "lasso"           
## [41] "leapBackward"     "leapForward"      "leapSeq"          "lmStepAIC"        "LMT"             
## [46] "LogitBoost"       "M5"               "M5Rules"          "msaenet"          "nodeHarvest"     
## [51] "OneR"             "ordinalNet"       "ordinalRF"        "ORFlog"           "ORFpls"          
## [56] "ORFridge"         "ORFsvm"           "pam"              "parRF"            "PART"            
## [61] "penalized"        "PenalizedLDA"     "qrf"              "ranger"           "Rborist"         
## [66] "relaxo"           "rf"               "rFerns"           "rfRules"          "rotationForest"  
## [71] "rotationForestCp" "rpart"            "rpart1SE"         "rpart2"           "rpartCost"       
## [76] "rpartScore"       "rqlasso"          "rqnc"             "RRF"              "RRFglobal"       
## [81] "sdwd"             "smda"             "sparseLDA"        "spikeslab"        "stepLDA"         
## [86] "stepQDA"          "wsrf"             "xgbDART"          "xgbLinear"        "xgbTree"
writeLines(embedded_names, "results/caret_feature_selection_models.txt")

3.2 Исходные данные

set.seed(233332)
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                 X5         
##  Min.   :-2.33681   Min.   :-2.2586   Min.   :-1.8520   Min.   :-1.83286   Min.   :-1.8567  
##  1st Qu.:-0.63214   1st Qu.:-0.8816   1st Qu.:-0.9007   1st Qu.:-0.68846   1st Qu.:-0.3613  
##  Median : 0.09398   Median : 0.1130   Median :-0.2494   Median :-0.06545   Median : 0.1450  
##  Mean   : 0.05918   Mean   :-0.1163   Mean   :-0.2176   Mean   :-0.02486   Mean   : 0.1461  
##  3rd Qu.: 0.81534   3rd Qu.: 0.5570   3rd Qu.: 0.3765   3rd Qu.: 0.75455   3rd Qu.: 0.7775  
##  Max.   : 2.42273   Max.   : 1.6827   Max.   : 1.4769   Max.   : 2.39491   Max.   : 2.3676

Матрица содержит 50 наблюдений и пять признаков. Каждый признак сгенерирован из стандартного нормального распределения. Классы A и B содержат по 25 наблюдений и назначены независимо от признаков.

3.3 Диаграммы размаха

p_box <- featurePlot(x = x, y = y, plot = "box",
                     auto.key = list(columns = 2),
                     scales = list(y = list(relation = "free")))
save_lattice(p_box, "caret_box.jpg")
knitr::include_graphics("graphs/caret_box.jpg")

Диаграммы размаха позволяют сравнить медианы, межквартильные интервалы и крайние наблюдения. Перекрытие распределений показывает, что однозначно разделить классы по отдельному признаку трудно. Точки за пределами усов — потенциальные выбросы по правилу диаграммы, а не доказанные ошибки данных.

3.4 Плотности распределений

p_density <- featurePlot(x = x, y = y, plot = "density",
                         auto.key = list(columns = 2),
                         scales = list(x = list(relation = "free"),
                                       y = list(relation = "free")))
save_lattice(p_density, "caret_density.jpg")
knitr::include_graphics("graphs/caret_density.jpg")

Плотности классов могут немного различаться из-за небольшого объёма выборки. Такие различия не подтверждают устойчивую связь с классом: по способу генерации зависимость отсутствует.

3.5 Матрица диаграмм рассеяния

p_pairs <- featurePlot(x = x, y = y, plot = "pairs",
                       auto.key = list(columns = 2))
save_lattice(p_pairs, "caret_pairs.jpg", width = 1800, height = 1800)
knitr::include_graphics("graphs/caret_pairs.jpg")

Матрица рассеяния показывает попарные связи признаков и размещение классов. По конструкции данные не содержат заданных кластеров или линейной связи между признаками. Случайные визуальные закономерности возможны.

3.6 Числовая проверка и вывод

class_means <- aggregate(as.data.frame(x), by = list(Class = y), FUN = mean)
knitr::kable(class_means, digits = 3, caption = "Средние значения по классам")
Средние значения по классам
Class X1 X2 X3 X4 X5
A 0.324 -0.216 -0.169 0.031 -0.059
B -0.206 -0.016 -0.266 -0.080 0.351
write.csv(class_means, "results/caret_class_means.csv", row.names = FALSE)

Вывод. featurePlot() позволяет рассмотреть данные несколькими способами. На этой искусственной выборке наблюдаемые отличия классов обусловлены случайной генерацией. Сделать вывод о полезности признака для новых данных по одним графикам нельзя. Для оценки качества модели нужна независимая проверка. Получены файлы caret_box.jpg, caret_density.jpg, caret_pairs.jpg.

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

4.1 Данные и методы

data(iris, package = "datasets")
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 ...
table(iris$Species)
## 
##     setosa versicolor  virginica 
##         50         50         50

Набор содержит 150 растений трёх видов, по 50 растений каждого вида. Цель классификации — Species. Предикторы: длина и ширина чашелистика (Sepal.Length, Sepal.Width) и лепестка (Petal.Length, Petal.Width).

Используем три энтропийных фильтра FSelector: информационный выигрыш, отношение выигрыша и симметричную неопределённость. Все оценивают отдельные признаки относительно класса; их абсолютные шкалы отличаются.

ig <- FSelector::information.gain(Species ~ ., iris, unit = "log2")
gr <- FSelector::gain.ratio(Species ~ ., iris)
su <- FSelector::symmetrical.uncertainty(Species ~ ., iris)
features <- rownames(ig)
importance <- data.frame(
  Feature = features,
  InformationGain_bits = ig$attr_importance,
  GainRatio = gr[features, "attr_importance"],
  SymmetricalUncertainty = su[features, "attr_importance"]
)
importance <- importance[order(-importance$InformationGain_bits), ]
rownames(importance) <- NULL
knitr::kable(importance, digits = 4, caption = "Оценки важности признаков iris")
Оценки важности признаков iris
Feature InformationGain_bits GainRatio SymmetricalUncertainty
Petal.Width 1.3784 0.8714 0.8705
Petal.Length 1.3565 0.8585 0.8572
Sepal.Length 0.6523 0.4196 0.4156
Sepal.Width 0.3856 0.2473 0.2453
write.csv(importance, "results/iris_importance.csv", row.names = FALSE)
top2 <- FSelector::cutoff.k(ig, 2)
selected_formula <- FSelector::as.simple.formula(top2, "Species")
top2
## [1] "Petal.Width"  "Petal.Length"
selected_formula
## Species ~ Petal.Width + Petal.Length
## <environment: 0x5561124c16c8>

4.2 График оценок

save_base(function() {
  par(mfrow = c(1, 3), mar = c(8, 4, 3, 1))
  for (column in names(importance)[-1]) {
    barplot(importance[[column]], names.arg = importance$Feature,
            las = 2, col = "steelblue", main = column, ylab = "Importance")
  }
}, "iris_importance.jpg", width = 2160)
knitr::include_graphics("graphs/iris_importance.jpg")

4.3 Вывод

Максимальный информационный выигрыш получил Petal.Width: 1.3784 бит. По cutoff.k(ig, 2) выбраны Petal.Width, Petal.Length. Полный порядок по информационному выигрышу: Petal.Width → Petal.Length → Sepal.Length → Sepal.Width.

Большая оценка означает более сильную индивидуальную связь признака с видом растения. Длина и ширина лепестка обычно хорошо различают виды iris, что следует сопоставить с полученной таблицей. Фильтр не оценивает все сочетания признаков; два высоко оценённых признака могут содержать общую информацию. Отбор здесь описательный, по полной выборке. При оценке классификатора отбор нужно выполнять только на обучающей части каждого разбиения, иначе возникает утечка информации. Повышение точности после отбора в данной работе не измерялось.

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

5.1 Четыре метода

В качестве непрерывной переменной выбрана длина чашелистика Sepal.Length. Во всех методах задаются три категории. Для fixed используются заранее выбранные границы 5.5 и 6.5 см; это учебные пороги, а не границы видов.

z <- iris$Sepal.Length
set.seed(233332)
disc <- list(
  interval = arules::discretize(z, method = "interval", breaks = 3),
  frequency = arules::discretize(z, method = "frequency", breaks = 3),
  cluster = arules::discretize(z, method = "cluster", breaks = 3, nstart = 25),
  fixed = arules::discretize(z, method = "fixed",
                            breaks = c(-Inf, 5.5, 6.5, Inf))
)
stopifnot(all(vapply(disc, function(v) length(v) == 150L &&
                      is.factor(v) && !anyNA(v), logical(1))))
disc_table <- do.call(rbind, lapply(names(disc), function(m) {
  freq <- table(disc[[m]])
  data.frame(Method = m, Interval = names(freq), Count = as.integer(freq))
}))
knitr::kable(disc_table, caption = "Интервалы и частоты после дискретизации")
Интервалы и частоты после дискретизации
Method Interval Count
interval [4.3,5.5) 52
interval [5.5,6.7) 70
interval [6.7,7.9] 28
frequency [4.3,5.4) 46
frequency [5.4,6.3) 53
frequency [6.3,7.9] 51
cluster [4.3,5.45) 52
cluster [5.45,6.46) 63
cluster [6.46,7.9] 35
fixed [-Inf,5.5) 52
fixed [5.5,6.5) 63
fixed [6.5, Inf] 35
write.csv(disc_table, "results/discretization_counts.csv", row.names = FALSE)
disc_values <- data.frame(Sepal.Length = z, as.data.frame(disc),
                         Species = iris$Species)
knitr::kable(head(disc_values, 12), caption = "Первые 12 преобразованных наблюдений")
Первые 12 преобразованных наблюдений
Sepal.Length interval frequency cluster fixed Species
5.1 [4.3,5.5) [4.3,5.4) [4.3,5.45) [-Inf,5.5) setosa
4.9 [4.3,5.5) [4.3,5.4) [4.3,5.45) [-Inf,5.5) setosa
4.7 [4.3,5.5) [4.3,5.4) [4.3,5.45) [-Inf,5.5) setosa
4.6 [4.3,5.5) [4.3,5.4) [4.3,5.45) [-Inf,5.5) setosa
5.0 [4.3,5.5) [4.3,5.4) [4.3,5.45) [-Inf,5.5) setosa
5.4 [4.3,5.5) [5.4,6.3) [4.3,5.45) [-Inf,5.5) setosa
4.6 [4.3,5.5) [4.3,5.4) [4.3,5.45) [-Inf,5.5) setosa
5.0 [4.3,5.5) [4.3,5.4) [4.3,5.45) [-Inf,5.5) setosa
4.4 [4.3,5.5) [4.3,5.4) [4.3,5.45) [-Inf,5.5) setosa
4.9 [4.3,5.5) [4.3,5.4) [4.3,5.45) [-Inf,5.5) setosa
5.4 [4.3,5.5) [5.4,6.3) [4.3,5.45) [-Inf,5.5) setosa
4.8 [4.3,5.5) [4.3,5.4) [4.3,5.45) [-Inf,5.5) setosa
write.csv(disc_values, "results/iris_discretized.csv", row.names = FALSE)

breaks — актуальное имя аргумента: прежний categories устарел. interval делит диапазон на равные отрезки; frequency использует квантили; cluster строит разбиение по k-means; fixed применяет заданные пороги. У frequency одинаковые значения не разделяются между категориями, поэтому численности могут отличаться. Бесконечные крайние границы fixed охватывают все исходные значения.

5.2 Распределения и границы

save_base(function() {
  par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
  for (m in names(disc)) {
    hist(z, breaks = 15, col = "lightblue", border = "white",
         main = m, xlab = "Sepal.Length (cm)", ylab = "Count")
    bounds <- attr(disc[[m]], "discretized:breaks")
    abline(v = bounds[is.finite(bounds)], col = "firebrick", lwd = 2, lty = 2)
  }
}, "iris_discretization.jpg", height = 1440)
knitr::include_graphics("graphs/iris_discretization.jpg")

Границы берутся из уже выполненной дискретизации, поэтому повторный запуск случайной кластеризации для рисунка не требуется.

for (m in names(disc)) {
  cat("\nМетод:", m, "\n")
  print(table(Category = disc[[m]], Species = iris$Species))
}
## 
## Метод: interval 
##            Species
## Category    setosa versicolor virginica
##   [4.3,5.5)     45          6         1
##   [5.5,6.7)      5         38        27
##   [6.7,7.9]      0          6        22
## 
## Метод: frequency 
##            Species
## Category    setosa versicolor virginica
##   [4.3,5.4)     40          5         1
##   [5.4,6.3)     10         31        12
##   [6.3,7.9]      0         14        37
## 
## Метод: cluster 
##              Species
## Category      setosa versicolor virginica
##   [4.3,5.45)      45          6         1
##   [5.45,6.46)      5         35        23
##   [6.46,7.9]       0          9        26
## 
## Метод: fixed 
##             Species
## Category     setosa versicolor virginica
##   [-Inf,5.5)     45          6         1
##   [5.5,6.5)       5         35        23
##   [6.5, Inf]      0          9        26

Таблицы по видам нужны для интерпретации: вид растения не используется при расчёте границ этих четырёх методов.

5.3 Вывод

Метод Что сохраняет Ограничение
interval Одинаковую ширину интервалов Численность категорий может сильно различаться
frequency Приблизительно равные численности Повторы мешают получить строго равные группы
cluster Группировку близких значений Зависит от распределения и случайной инициализации
fixed Заданные интерпретируемые пороги Пороги требуют содержательного обоснования

В данном запуске численности категорий составили:

  • interval: 52, 70, 28 наблюдений.
  • frequency: 46, 53, 51 наблюдений.
  • cluster: 52, 63, 35 наблюдений.
  • fixed: 52, 63, 35 наблюдений.

Дискретизация преобразовала значения в факторы и уменьшила детализацию. Ни один из методов сам по себе не оптимизирует качество классификации. Для новых наблюдений нужно применять границы, рассчитанные по обучающей выборке, а не вычислять новые границы отдельно на тестовой выборке.

6 Boruta: выбор признаков Ozone

6.1 Подготовка данных

data("Ozone", package = "mlbench")
dim(Ozone)
## [1] 366  13
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 ...
missing_counts <- colSums(is.na(Ozone))
knitr::kable(data.frame(Feature = names(missing_counts),
                        Missing = as.integer(missing_counts)),
             caption = "Пропущенные значения Ozone")
Пропущенные значения Ozone
Feature Missing
V1 0
V2 0
V3 0
V4 5
V5 12
V6 0
V7 15
V8 2
V9 139
V10 15
V11 1
V12 14
V13 0
ozo <- na.omit(Ozone)
removed_n <- nrow(Ozone) - nrow(ozo)
c(original = nrow(Ozone), complete = nrow(ozo), removed = removed_n)
## original complete  removed 
##      366      203      163
stopifnot(is.numeric(ozo$V4), !anyNA(ozo))

Ozone содержит 366 ежедневных наблюдений за загрязнением воздуха в Лос-Анджелесе в 1976 году и 13 переменных. Целевая переменная V4 — максимальное за день среднее часовое значение озона. В отличие от предыдущей задачи это регрессия; преобразовывать V4 в фактор без дополнительного условия не нужно.

После удаления строк с пропусками осталось 203 наблюдений, удалено 163. Метод прост, но уменьшает выборку и может внести смещение, если пропуски неслучайны.

Переменная Смысл
V1–V3 Месяц, день месяца, день недели
V4 Целевое значение озона
V5 Высота поверхности давления 500 мбар
V6–V7 Скорость ветра и влажность
V8–V9 Температура на двух станциях
V10 Высота основания инверсии
V11 Градиент давления
V12 Температура основания инверсии
V13 Видимость

Сохраняется исходное числовое кодирование календарных признаков, как в справочном примере Boruta. Специальные календарные преобразования не выполняются.

6.2 Запуск алгоритма

set.seed(233332)
boruta_fit <- Boruta::Boruta(V4 ~ ., data = ozo,
                             maxRuns = 300, doTrace = 0,
                             holdHistory = TRUE,
                             getImp = Boruta::getImpRfZ, num.threads = 1)
print(boruta_fit)
## Boruta performed 30 iterations in 3.657959 secs.
##  9 attributes confirmed important: V1, V10, V11, V12, V13 and 4 more;
##  3 attributes confirmed unimportant: V2, V3, V6;
initial_stats <- Boruta::attStats(boruta_fit)
knitr::kable(initial_stats, digits = 3,
             caption = "Статистика важности и исходные решения Boruta")
Статистика важности и исходные решения Boruta
meanImp medianImp minImp maxImp normHits decision
V1 9.591 9.497 8.266 11.353 1.000 Confirmed
V2 1.018 0.827 -0.743 3.845 0.200 Rejected
V3 -1.746 -2.081 -2.734 -0.483 0.000 Rejected
V5 8.923 8.914 7.222 10.419 1.000 Confirmed
V6 1.426 1.373 0.069 2.937 0.067 Rejected
V7 11.646 11.628 9.903 12.827 1.000 Confirmed
V8 17.079 16.961 16.150 18.389 1.000 Confirmed
V9 19.287 19.420 17.348 20.662 1.000 Confirmed
V10 10.036 9.928 8.952 11.330 1.000 Confirmed
V11 12.026 11.956 10.333 13.585 1.000 Confirmed
V12 14.681 14.697 13.750 15.565 1.000 Confirmed
V13 9.671 9.760 8.467 10.643 1.000 Confirmed

Boruta сравнивает важность признаков со случайными теневыми признаками, которые получаются перестановкой значений. Источник важности явно задан как getImpRfZ: он использует случайный лес из ranger. Решения: Confirmed — подтверждён, Rejected — отклонён, Tentative — пока не определён.

6.3 Диаграмма размаха важности

save_base(function() {
  par(mar = c(8, 4, 3, 1))
  plot(boruta_fit, las = 2, cex.axis = 0.8,
       main = "Boruta: Ozone (response V4)", xlab = "")
}, "ozone_boruta_boxplot.jpg", width = 2160, height = 1440)
knitr::include_graphics("graphs/ozone_boruta_boxplot.jpg")

plot() для объекта Boruta строит именно boxplot распределений важности по итерациям. Зелёный цвет обозначает подтверждённые признаки, красный — отклонённые, жёлтый — неопределённые, синий — теневые. Этот рисунок отражает исходное статистическое решение до дополнительной обработки.

6.4 Неопределённые признаки и итоговое подмножество

had_tentative <- any(boruta_fit$finalDecision == "Tentative")
boruta_final <- if (had_tentative) Boruta::TentativeRoughFix(boruta_fit) else boruta_fit
final_stats <- Boruta::attStats(boruta_final)
final_table <- data.frame(Feature = rownames(final_stats), final_stats,
                          InitialDecision = initial_stats[rownames(final_stats), "decision"],
                          row.names = NULL)
final_table <- final_table[order(-final_table$medianImp), ]
knitr::kable(final_table, digits = 3,
             caption = "Итоговые оценки; исходные решения сохранены для сравнения")
Итоговые оценки; исходные решения сохранены для сравнения
Feature meanImp medianImp minImp maxImp normHits decision InitialDecision
8 V9 19.287 19.420 17.348 20.662 1.000 Confirmed Confirmed
7 V8 17.079 16.961 16.150 18.389 1.000 Confirmed Confirmed
11 V12 14.681 14.697 13.750 15.565 1.000 Confirmed Confirmed
10 V11 12.026 11.956 10.333 13.585 1.000 Confirmed Confirmed
6 V7 11.646 11.628 9.903 12.827 1.000 Confirmed Confirmed
9 V10 10.036 9.928 8.952 11.330 1.000 Confirmed Confirmed
12 V13 9.671 9.760 8.467 10.643 1.000 Confirmed Confirmed
1 V1 9.591 9.497 8.266 11.353 1.000 Confirmed Confirmed
4 V5 8.923 8.914 7.222 10.419 1.000 Confirmed Confirmed
5 V6 1.426 1.373 0.069 2.937 0.067 Rejected Rejected
2 V2 1.018 0.827 -0.743 3.845 0.200 Rejected Rejected
3 V3 -1.746 -2.081 -2.734 -0.483 0.000 Rejected Rejected
confirmed_initial <- Boruta::getSelectedAttributes(boruta_fit, withTentative = FALSE)
confirmed <- Boruta::getSelectedAttributes(boruta_final, withTentative = FALSE)
rejected <- names(boruta_final$finalDecision)[boruta_final$finalDecision == "Rejected"]
unresolved_initial <- names(boruta_fit$finalDecision)[boruta_fit$finalDecision == "Tentative"]
write.csv(final_table, "results/ozone_boruta_stats.csv", row.names = FALSE)
writeLines(confirmed, "results/ozone_selected_features.txt")
saveRDS(boruta_fit, "results/ozone_boruta_original.rds")
saveRDS(boruta_final, "results/ozone_boruta_final.rds")

Вывод по выполненному запуску. Исходно подтверждены V1, V5, V7, V8, V9, V10, V11, V12, V13. Неопределёнными остались нет. TentativeRoughFix() не потребовался. Это дополнительное эвристическое решение по медианам важности, а не эквивалент продолжения основного статистического теста.

Итоговое подмножество: V1, V5, V7, V8, V9, V10, V11, V12, V13 (9 признаков); отклонены V2, V3, V6. Лидирует по медиане важности V9. Нужно опираться на таблицу и boxplot, а не заранее назначать важные признаки. Алгоритм ищет все релевантные признаки, поэтому может сохранять несколько связанных между собой предикторов. Важность не доказывает причинную связь.

Выбор признаков выполнен на всём очищенном наборе для изучения алгоритма. Изменение ошибки прогноза не измерялось. При последующей проверке модели отбор следует повторять внутри обучающего разбиения; поскольку наблюдения ежедневные, для прогноза будущих дней разумно учитывать порядок времени.

7 Общий вывод

В работе получены три графика разведочного анализа случайных данных CARET, таблица и график важности признаков iris, четыре разбиения длины чашелистика на категории и boxplot выбора признаков Ozone. Все рисунки сохранены в JPG.

Разведочный анализ показывает распределения и возможные связи, фильтры FSelector оценивают отдельные признаки, а Boruta сравнивает их важность со случайным ориентиром. Дискретизация меняет представление данных и не заменяет выбор признаков. Результаты отбора требуют проверки предсказательной полезности на новых данных.

8 Полученные файлы

Сохранённые графики
File Size_KB
caret_box.jpg 166.5
caret_density.jpg 233.7
caret_pairs.jpg 534.8
iris_discretization.jpg 255.1
iris_importance.jpg 198.1
ozone_boruta_boxplot.jpg 156.8

9 Публикация на RPubs

  1. Открыть этот .Rmd в RStudio и установить пакеты из раздела подготовки.
  2. Нажать Knit → Knit to HTML. Проверить, что код завершился без ошибок, в документе видны таблицы и все шесть рисунков.
  3. В окне предварительного просмотра нажать Publish → RPubs.
  4. Войти в свой аккаунт RPubs, указать название отчёта и завершить публикацию.
  5. Скопировать адрес опубликованного документа и отправить преподавателю.

HTML содержит встроенные изображения и подходит для публикации. Ссылка появляется после фактической публикации: локальный HTML её не заменяет.

10 Источники

  1. Документация CARET: визуализация.
  2. Справочник FSelector.
  3. Справка arules: discretize.
  4. Справка mlbench: Ozone.
  5. Справка Boruta.
  6. Kursa M. B., Rudnicki W. R. Feature Selection with the Boruta Package. Journal of Statistical Software, 2010, 36(11).
  7. Инструкция RPubs по публикации.

11 Версии программной среды

versions <- data.frame(Package = required,
                       Version = vapply(required, function(p)
                         as.character(packageVersion(p)), character(1)))
knitr::kable(versions)
Package Version
caret caret 6.0.94
FSelector FSelector 0.34
arules arules 1.7.7
Boruta Boruta 8.0.0
mlbench mlbench 2.1.3.1
ranger ranger 0.16.0
randomForest randomForest 4.7.1.1
knitr knitr 1.45
rmarkdown rmarkdown 2.25
sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 24.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8          LC_NUMERIC=C              LC_TIME=C.UTF-8          
##  [4] LC_COLLATE=C.UTF-8        LC_MONETARY=C.UTF-8       LC_MESSAGES=C.UTF-8      
##  [7] LC_PAPER=C.UTF-8          LC_NAME=C.UTF-8           LC_ADDRESS=C.UTF-8       
## [10] LC_TELEPHONE=C.UTF-8      LC_MEASUREMENT=C.UTF-8    LC_IDENTIFICATION=C.UTF-8
## 
## time zone: Europe/Paris
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] caret_6.0-94   lattice_0.22-5 ggplot2_3.4.4 
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.4         xfun_0.41            bslib_0.6.1          recipes_1.0.9       
##  [5] rJava_1.0-11         vctrs_0.6.5          tools_4.3.3          generics_0.1.3      
##  [9] stats4_4.3.3         parallel_4.3.3       tibble_3.2.1         fansi_1.0.5         
## [13] highr_0.10           FSelector_0.34       pkgconfig_2.0.3      ModelMetrics_1.2.2.2
## [17] Matrix_1.6-5         data.table_1.14.10   lifecycle_1.0.4      compiler_4.3.3      
## [21] stringr_1.5.1        mlbench_2.1-3.1      munsell_0.5.0        codetools_0.2-19    
## [25] Boruta_8.0.0         arules_1.7-7         htmltools_0.5.7      class_7.3-22        
## [29] sass_0.4.8           yaml_2.3.8           prodlim_2023.08.28   pillar_1.9.0        
## [33] jquerylib_0.1.4      MASS_7.3-60.0.1      entropy_1.3.2        cachem_1.0.8        
## [37] gower_1.0.1          iterators_1.0.14     rpart_4.1.23         foreach_1.5.2       
## [41] nlme_3.1-164         parallelly_1.37.1    lava_1.7.3           tidyselect_1.2.0    
## [45] digest_0.6.34        stringi_1.8.3        future_1.33.1        dplyr_1.1.4         
## [49] reshape2_1.4.4       purrr_1.0.2          listenv_0.9.1        RWekajars_3.9.3-2   
## [53] splines_4.3.3        fastmap_1.1.1        grid_4.3.3           colorspace_2.1-0    
## [57] cli_3.6.2            magrittr_2.0.3       randomForest_4.7-1.1 survival_3.5-8      
## [61] utf8_1.2.4           future.apply_1.11.1  withr_2.5.0          scales_1.3.0        
## [65] lubridate_1.9.3      timechange_0.3.0     rmarkdown_2.25       globals_0.16.2      
## [69] nnet_7.3-19          RWeka_0.4-50         timeDate_4032.109    ranger_0.16.0       
## [73] evaluate_0.23        knitr_1.45           hardhat_1.3.1        rlang_1.1.3         
## [77] Rcpp_1.0.12          glue_1.7.0           pROC_1.18.5          ipred_0.9-14        
## [81] jsonlite_1.8.8       R6_2.5.1             plyr_1.8.9