Лабораторная работа 2
Работа с пакетами CARET и BORUTA на языке R
Редкоус Алексей Алексеевич, группа 233-332
17 октября 2026 г.
Московский политехнический университет
Дисциплина: Анализ данных
Преподаватель: Попов Дмитрий Иванович
1 Цель работы
Я изучил средства графического анализа и отбора признаков в R. Построил графики с помощью caret, оценил важность признаков iris через FSelector, сравнил четыре способа дискретизации и выполнил отбор признаков Ozone алгоритмом Boruta.
2 Графический анализ с caret
2.1 Доступные модели и исходные данные
Я получил список моделей caret командой
names(getModelInfo()). Для отбора признаков в пакете
предусмотрены, в частности, рекурсивное исключение признаков
rfe() и фильтрация sbf(). В этой части я
исследовал признаки графически с помощью featurePlot().
model_names <- names(caret::getModelInfo())
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"
length(model_names)
#> [1] 239
writeLines(model_names, "Выполнение/results/caret_models.txt")Я создал матрицу из 50 наблюдений и пяти случайных признаков. Классы A и B содержат по 25 наблюдений. Перед генерацией зафиксировал начальное состояние генератора, чтобы повторный запуск давал тот же результат.
set.seed(233332)
x <- matrix(rnorm(50*5), ncol=5)
y <- factor(rep(c("A", "B"), 25))
colnames(x) <- paste0("X", seq_len(ncol(x)))
dim(x)
#> [1] 50 5
table(y)
#> y
#> A B
#> 25 25
kable(head(data.frame(x, Class=y)), digits=3,
caption="Первые шесть наблюдений")| X1 | X2 | X3 | X4 | X5 | Class |
|---|---|---|---|---|---|
| 0.377 | -0.646 | -0.331 | 0.166 | 0.346 | A |
| 0.391 | 0.353 | -0.405 | 0.969 | 1.745 | B |
| -2.337 | -1.390 | -0.174 | -0.997 | 0.787 | A |
| -0.554 | -0.918 | -0.930 | 1.264 | -0.091 | B |
| 1.023 | -0.578 | 0.624 | -0.029 | -1.308 | A |
| -0.632 | 1.443 | -1.852 | -1.031 | -0.253 | B |
2.2 Матрица диаграмм рассеяния
print(caret::featurePlot(x=x, y=y, plot="pairs",
auto.key=list(columns=2), pch=16,
col=c("#2463A0", "#D45A32")))
Рисунок 1. Попарные зависимости признаков
2.3 Распределения признаков
print(caret::featurePlot(x=x, y=y, plot="density", layout=c(3,2),
scales=list(x=list(relation="free"), y=list(relation="free")),
auto.key=list(columns=2), adjust=1.2))
Рисунок 2. Плотности распределений в классах A и B
print(caret::featurePlot(x=x, y=y, plot="box", layout=c(3,2),
scales=list(y=list(relation="free"))))
Рисунок 3. Распределения признаков по классам
На диаграммах я сравнил взаимное расположение классов, медианы, разброс и плотности распределений. Области A и B заметно перекрываются. Отдельные различия медиан и формы плотностей возникли в конкретной случайной выборке: признаки сгенерированы независимо от меток классов. Поэтому по этим графикам я не сделал вывода об устойчивой способности признаков разделять A и B.
3 Важность признаков iris с FSelector
Я использовал четыре числовых признака iris для классификации по
столбцу Species. В наборе 150 цветков, по 50 каждого вида.
Для сравнения вычислил информационный выигрыш, отношение информационного
выигрыша, симметричную неопределённость и ReliefF. В каждом методе
большее значение соответствует большей важности; численные шкалы разных
методов напрямую не сравнивал.
data(iris)
stopifnot(!anyNA(iris), nrow(iris)==150)
set.seed(233333)
weights <- list(
InformationGain=FSelector::information.gain(Species~., iris),
GainRatio=FSelector::gain.ratio(Species~., iris),
SymmetricalUncertainty=FSelector::symmetrical.uncertainty(Species~., iris),
ReliefF=FSelector::relief(Species~., iris, neighbours.count=10, sample.size=150)
)
importance <- data.frame(Feature=names(iris)[1:4])
for (method in names(weights)) {
importance[[method]] <- weights[[method]][importance$Feature, 1]
}
rank_matrix <- sapply(importance[-1], function(v) rank(-v, ties.method="average"))
importance$MeanRank <- rowMeans(rank_matrix)
importance <- importance[order(importance$MeanRank), ]
kable(importance, digits=4, row.names=FALSE,
caption="Важность признаков iris и средний ранг")| Feature | InformationGain | GainRatio | SymmetricalUncertainty | ReliefF | MeanRank |
|---|---|---|---|---|---|
| Petal.Width | 0.9554 | 0.8714 | 0.8705 | 0.3548 | 1.00 |
| Petal.Length | 0.9403 | 0.8585 | 0.8572 | 0.3462 | 2.00 |
| Sepal.Length | 0.4521 | 0.4196 | 0.4156 | 0.1339 | 3.25 |
| Sepal.Width | 0.2673 | 0.2473 | 0.2453 | 0.1414 | 3.75 |
selected_iris <- FSelector::cutoff.k(weights$InformationGain, 2)
selected_iris
#> [1] "Petal.Width" "Petal.Length"
write.csv(importance, "Выполнение/results/iris_importance.csv", row.names=FALSE)op <- par(mfrow=c(2,2), mar=c(6,4,3,1))
for (method in names(weights)) {
v <- weights[[method]][,1]
names(v) <- rownames(weights[[method]])
barplot(sort(v, decreasing=TRUE), las=2, col="#437CB4",
main=method, ylab="Importance", cex.names=0.85)
}
Рисунок 4. Сравнение оценок важности iris
По информационному выигрышу я выбрал два признака: Petal.Width и Petal.Length. Сравнение рангов показало, что наибольшую информацию о виде цветка несут размеры лепестков. Длина и ширина чашелистика оказались менее важными. Эти оценки показывают связь признаков с классом на исследуемом наборе; точность отдельной модели классификации в этой части я не оценивал.
4 Дискретизация iris с arules
Я преобразовал длину чашелистика Sepal.Length из
числовой переменной в категориальную. Для сопоставимости использовал три
категории в каждом методе. Для fixed задал границы 5.5 и
6.5 см, крайние интервалы оставил открытыми до бесконечности. В
остальных методах границы определяются по данным.
values <- iris$Sepal.Length
set.seed(233334)
disc <- list(
interval=arules::discretize(values, method="interval", breaks=3),
frequency=arules::discretize(values, method="frequency", breaks=3),
cluster=arules::discretize(values, method="cluster", breaks=3, nstart=50),
fixed=arules::discretize(values, method="fixed", breaks=c(-Inf,5.5,6.5,Inf))
)
breaks <- lapply(disc, attr, "discretized:breaks")
counts <- do.call(rbind, lapply(names(disc), function(method) {
t <- table(disc[[method]])
data.frame(Method=method, Interval=names(t), Count=as.integer(t))
}))
kable(counts, row.names=FALSE, 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 |
kable(data.frame(Value=values[1:10], lapply(disc, as.character) |>
lapply(function(v) head(v,10))),
caption="Первые десять результатов преобразования")| Value | interval | frequency | cluster | fixed |
|---|---|---|---|---|
| 5.1 | [4.3,5.5) | [4.3,5.4) | [4.3,5.45) | [-Inf,5.5) |
| 4.9 | [4.3,5.5) | [4.3,5.4) | [4.3,5.45) | [-Inf,5.5) |
| 4.7 | [4.3,5.5) | [4.3,5.4) | [4.3,5.45) | [-Inf,5.5) |
| 4.6 | [4.3,5.5) | [4.3,5.4) | [4.3,5.45) | [-Inf,5.5) |
| 5.0 | [4.3,5.5) | [4.3,5.4) | [4.3,5.45) | [-Inf,5.5) |
| 5.4 | [4.3,5.5) | [5.4,6.3) | [4.3,5.45) | [-Inf,5.5) |
| 4.6 | [4.3,5.5) | [4.3,5.4) | [4.3,5.45) | [-Inf,5.5) |
| 5.0 | [4.3,5.5) | [4.3,5.4) | [4.3,5.45) | [-Inf,5.5) |
| 4.4 | [4.3,5.5) | [4.3,5.4) | [4.3,5.45) | [-Inf,5.5) |
| 4.9 | [4.3,5.5) | [4.3,5.4) | [4.3,5.45) | [-Inf,5.5) |
stopifnot(all(vapply(disc, function(v) is.factor(v) && !anyNA(v) && length(v)==150, logical(1))))
write.csv(counts, "Выполнение/results/discretization_counts.csv", row.names=FALSE)
write.csv(data.frame(Sepal.Length=values, disc, Species=iris$Species),
"Выполнение/results/discretized_iris.csv", row.names=FALSE)op <- par(mfrow=c(2,2), mar=c(4,4,3,1))
for (method in names(disc)) {
hist(values, breaks=seq(4,8,by=0.2), col="#B6D0E6", border="white",
main=method, xlab="Sepal.Length (cm)", ylab="Count")
abline(v=breaks[[method]][is.finite(breaks[[method]])], col="#BC3E32", lwd=2, lty=2)
}
Рисунок 5. Границы интервалов при четырёх способах дискретизации
count_matrix <- sapply(disc, function(v) as.integer(table(v)))
barplot(count_matrix, beside=TRUE,
col=c("#3977A8", "#82B0D4", "#D5E6F2"),
legend.text=paste("Категория", 1:3),
args.legend=list(x="topright", bty="n"),
xlab="Метод", ylab="Число наблюдений", ylim=c(0,max(count_matrix)+20))
Рисунок 6. Численность полученных категорий
Метод interval разделил диапазон на равные по ширине
интервалы. Метод frequency сформировал группы близкой
численности; точному равенству мешают повторяющиеся значения длины.
Метод cluster нашёл группы с близкими значениями, используя
k-means, поэтому ширина интервалов и их численность различаются. Метод
fixed применил заданные мной границы. Я убедился, что все
150 наблюдений получили категорию без пропусков. Выбор метода зависит от
задачи: для близкой численности групп подходит frequency,
для заранее заданных порогов — fixed.
5 Отбор признаков Ozone с Boruta
5.1 Подготовка данных
Я загрузил набор Ozone из пакета mlbench. Наблюдения содержат
календарные и метеорологические характеристики дней в Лос-Анджелесе в
1976 году. В качестве целевой переменной взял V4 — суточный
максимум средней за один час концентрации озона. Целевая переменная
числовая, поэтому отбор выполняется для задачи регрессии.
data("Ozone", package="mlbench")
missing <- data.frame(Variable=names(Ozone), Missing=colSums(is.na(Ozone)))
kable(missing, row.names=FALSE, caption="Пропуски в исходных данных")| Variable | 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 |
ozone <- na.omit(Ozone)
preparation <- data.frame(Original=nrow(Ozone), Removed=nrow(Ozone)-nrow(ozone),
Used=nrow(ozone), Predictors=ncol(ozone)-1)
kable(preparation, row.names=FALSE, caption="Размер данных до и после удаления неполных строк")| Original | Removed | Used | Predictors |
|---|---|---|---|
| 366 | 163 | 203 | 12 |
stopifnot(!anyNA(ozone), is.numeric(ozone$V4))
write.csv(missing, "Выполнение/results/ozone_missing.csv", row.names=FALSE)Я исключил строки с пропусками, чтобы все признаки оценивались на одной полной выборке. После обработки осталось 203 наблюдений из 366. Такой отбор уменьшает объём данных; результаты относятся к оставшимся полным наблюдениям.
5.2 Запуск и результаты
Алгоритм Boruta сравнивает важность исходных признаков с важностью их случайно перемешанных копий. Для устойчивого повторения я зафиксировал seed, использовал 500 деревьев и один вычислительный поток. Максимальное число итераций задал равным 300, историю важностей сохранил для boxplot.
set.seed(233335)
boruta_result <- Boruta::Boruta(V4~., data=ozone,
maxRuns=300, doTrace=0, holdHistory=TRUE,
getImp=Boruta::getImpRfZ,
num.trees=500, num.threads=1)
print(boruta_result)
#> Boruta performed 54 iterations in 4.0181 secs.
#> 9 attributes confirmed important: V1, V10, V11, V12, V13 and 4 more;
#> 3 attributes confirmed unimportant: V2, V3, V6;
boruta_stats <- Boruta::attStats(boruta_result)
boruta_stats$Variable <- rownames(boruta_stats)
boruta_stats <- boruta_stats[order(-boruta_stats$meanImp),
c("Variable","meanImp","medianImp","minImp","maxImp",
"normHits","decision")]
rownames(boruta_stats) <- NULL
kable(boruta_stats, digits=3, row.names=FALSE,
caption="Результаты отбора признаков Ozone")| Variable | meanImp | medianImp | minImp | maxImp | normHits | decision |
|---|---|---|---|---|---|---|
| V9 | 19.481 | 19.497 | 18.026 | 22.050 | 1.000 | Confirmed |
| V8 | 17.287 | 17.368 | 15.700 | 18.625 | 1.000 | Confirmed |
| V12 | 14.813 | 14.853 | 13.288 | 16.412 | 1.000 | Confirmed |
| V11 | 12.174 | 12.101 | 10.652 | 13.647 | 1.000 | Confirmed |
| V7 | 12.006 | 12.093 | 10.295 | 13.753 | 1.000 | Confirmed |
| V10 | 10.114 | 9.972 | 8.518 | 11.896 | 1.000 | Confirmed |
| V1 | 9.565 | 9.674 | 7.182 | 11.145 | 1.000 | Confirmed |
| V13 | 9.461 | 9.450 | 7.813 | 11.291 | 1.000 | Confirmed |
| V5 | 8.951 | 8.928 | 7.244 | 10.231 | 1.000 | Confirmed |
| V6 | 1.291 | 1.346 | -1.390 | 3.460 | 0.278 | Rejected |
| V2 | 1.241 | 1.244 | -0.690 | 2.671 | 0.056 | Rejected |
| V3 | -1.091 | -0.987 | -2.036 | 0.093 | 0.000 | Rejected |
confirmed <- Boruta::getSelectedAttributes(boruta_result, withTentative=FALSE)
tentative <- names(boruta_result$finalDecision)[boruta_result$finalDecision=="Tentative"]
rejected <- names(boruta_result$finalDecision)[boruta_result$finalDecision=="Rejected"]
write.csv(boruta_stats, "Выполнение/results/boruta_importance.csv", row.names=FALSE)
saveRDS(boruta_result, "Выполнение/results/boruta_result.rds")op <- par(mar=c(8,4,2,1))
plot(boruta_result, las=2, cex.axis=0.85, xlab="", ylab="Importance (Z-score)")
legend("topleft", legend=c("Confirmed", "Tentative", "Rejected", "Shadow"),
fill=c("green", "yellow", "red", "blue"), bty="n", cex=0.8)
Рисунок 7. Важности исходных и теневых признаков Boruta
На boxplot зелёным показаны подтверждённые признаки, красным — отвергнутые, жёлтым — неопределённые, синим — теневые признаки. Я получил 9 подтверждённых, 3 отвергнутых и 0 неопределённых признаков.
Подтверждённые признаки: V1, V5, V7, V8, V9, V10, V11, V12,
V13.
Отвергнутые признаки: V2, V3, V6.
Неопределённые признаки: нет.
| Variable | Meaning | Decision |
|---|---|---|
| V9 | Температура в El Monte | Confirmed |
| V8 | Температура в Sandburg | Confirmed |
| V12 | Температура основания инверсии | Confirmed |
| V11 | Градиент давления | Confirmed |
| V7 | Влажность | Confirmed |
| V10 | Высота основания инверсии | Confirmed |
| V1 | Месяц | Confirmed |
| V13 | Видимость | Confirmed |
| V5 | Высота поверхности 500 мбар | Confirmed |
| V6 | Скорость ветра | Rejected |
| V2 | День месяца | Rejected |
| V3 | День недели | Rejected |
Наибольшую среднюю важность получил признак V9 (Температура в El Monte). Я отобрал признаки, для которых Boruta подтвердил преимущество перед случайными копиями. Важность отражает полезность признака для предсказания в этой модели и не доказывает причинное влияние на концентрацию озона.
6 Вывод
Я выполнил графический анализ случайных данных, сравнил четыре оценки важности iris и четыре способа дискретизации длины чашелистика. Для iris наиболее информативными оказались размеры лепестков. На данных Ozone я обработал пропуски, запустил Boruta и получил набор подтверждённых признаков для прогнозирования концентрации озона. Код, таблицы, графики и выводы объединил в один R Markdown. Все результаты получены при выполнении приведённых блоков кода.