Изучить возможности языка R для анализа данных и отбора важных признаков. Научиться использовать пакеты CARET, FSelector, arules и Boruta, строить графики, сравнивать полученные результаты и создавать отчёты в R Markdown.
Для выполнения лабораторной работы я использовал язык программирования R и среду разработки RStudio. Сначала я установил необходимые пакеты, которые понадобятся для выполнения заданий. Выполнил в консоли команду установки. (Листинг 1):
Листинг 1. Установка необходимых пакетов.
install.packages(c(
"caret", "FSelector", "arules", "Boruta",
"mlbench", "rmarkdown", "knitr"
))
В процессе установки выяснилось что 3 пакета не доступны к скачиванию, потому что заблокировано для России… Чтобы их установить я открыл репозитории RStudio через VPN в браузере и скачал пакеты которые RStudio не смог скачать на ноутбук. Потом в Tools выбрал Install packages… И выполнил установку из архивов *.tar.gz. После завершения установки я подключил пакеты с помощью функции library() (Листинг 2):
Листинг 2. Подключение пакетов.
library(caret)
library(FSelector)
library(arules)
library(Boruta)
library(mlbench)
Функция install.packages() устанавливает пакеты, а library() позволяет использовать их функции в программе. После этого я перехожу к выполнению практических заданий. (Рисунок 1).
## Пакет Версия
## 1 caret 7.0.1
## 2 FSelector 0.34
## 3 arules 1.7.15
## 4 Boruta 10.0.0
## 5 mlbench 2.1.11
Рисунок 1 - Установка и подключение пакетов в RStudio.
В первом задании мне нужно было посмотреть доступные методы пакета CARET и построить графики для случайных данных. Сначала я вывел список названий моделей с помощью функции getModelInfo() (Листинг 3):
Листинг 3. Получение списка методов CARET.
names(getModelInfo())
В консоли появился большой список названий методов, которые поддерживаются пакетом CARET. Я посмотрел этот список и перешёл к построению графиков (Рисунок 2).
## [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"
Рисунок 2 - Список доступных методов пакета CARET.
Далее я создал данные по примеру из задания. Использовал 50 наблюдений, пять числовых признаков и две группы: A и B. Команда set.seed(123) нужна, чтобы при повторном запуске случайные числа не менялись (Листинг 4).
Листинг 4. Создание случайных данных для графиков.
set.seed(123)
x <- matrix(rnorm(50 * 5), ncol = 5)
y <- factor(rep(c("A", "B"), 25))
dim(x)
table(y)
head(x)
По результату видно, что матрица содержит 50 строк и 5 столбцов. В каждой группе получилось по 25 наблюдений. Команда head(x) показывает первые строки данных.
## [1] 50 5
## y
## A B
## 25 25
## [,1] [,2] [,3] [,4] [,5]
## [1,] -0.56047565 0.25331851 -0.71040656 0.7877388 2.1988103
## [2,] -0.23017749 -0.02854676 0.25688371 0.7690422 1.3124130
## [3,] 1.55870831 -0.04287046 -0.24669188 0.3322026 -0.2651451
## [4,] 0.07050839 1.36860228 -0.34754260 -1.0083766 0.5431941
## [5,] 0.12928774 -0.22577099 -0.95161857 -0.1194526 -0.4143399
## [6,] 1.71506499 1.51647060 -0.04502772 -0.2803953 -0.4762469
Рисунок 3 - Проверка размеров матрицы и количества наблюдений в группах.
Сначала я построил парные диаграммы. Они позволяют посмотреть, как связаны друг с другом разные признаки. Для этого я использовал функцию featurePlot() с параметром pairs (Листинг 5).
Листинг 5. Построение парных диаграмм и сохранение в JPG.
jpeg("01_pairs.jpg", width = 1600, height = 1200)
print(featurePlot(x, y, plot = "pairs"))
dev.off()
После выполнения кода график сохранился в файл 01_pairs.jpg. По нему можно сравнить расположение точек двух групп для разных пар признаков (Рисунок 4).
Рисунок 4 - Парные диаграммы признаков для групп A и B.
Затем я построил диаграммы «ящик с усами». Они нужны, чтобы сравнить разброс значений, медиану и возможные выбросы в двух группах. (Листинг 6).
Листинг 6. Построение диаграмм «ящик с усами» и сохранение в JPG
jpeg("02_boxplot.jpg", width = 1600, height = 1200)
print(featurePlot(x, y, plot = "box"))
dev.off()
Я получил диаграммы для пяти признаков и сравнил по ним распределение значений в группах A и B (Рисунок 5).
Рисунок 5 - Диаграммы «ящик с усами» для групп A и B.
После этого я построил графики плотности распределения. Они показывают, в каких диапазонах чаще встречаются значения признаков. (Листинг 7).
Листинг 7. Построение графиков плотности и сохранение в JPG
jpeg("03_density.jpg", width = 1600, height = 1200)
print(featurePlot(x, y, plot = "density"))
dev.off()
На графиках я увидел кривые распределения для групп A и B. По ним можно сравнить, насколько сильно распределения двух групп отличаются (Рисунок 6).
Рисунок 6 - Графики плотности распределения пяти признаков.
Последним вариантом графика я выбрал strip. Такой график показывает расположение отдельных значений признаков для каждой группы. Для того чтобы точки не накладывались друг на друга, я использовал параметр jitter = TRUE (Листинг 8).
Листинг 8. Построение точечных диаграмм и сохранение в JPG
jpeg("04_strip.jpg", width = 1600, height = 1200)
print(featurePlot(x, y, plot = "strip", jitter = TRUE))
dev.off()
После выполнения программы я получил точечные диаграммы для пяти признаков. По ним можно сравнить разброс значений в группах A и B (Рисунок 7).
Рисунок 7 - Точечные диаграммы распределения признаков для групп A и B.
Вывод по первому заданию: Я познакомился с пакетом CARET и научился строить несколько видов графиков с помощью featurePlot(). Все четыре графика сохранил в файлы JPG. Так как числа были случайными, чёткого разделения групп A и B на графиках не наблюдается.
Во втором задании я использовал набор данных iris. В нём находятся измерения цветков: длина и ширина чашелистика, длина и ширина лепестка, а также вид цветка.
Сначала я загрузил набор данных и посмотрел первые строки. (Листинг 9).
Листинг 9. Загрузка и просмотр набора данных iris.
data(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
Далее мне нужно было определить, какие признаки больше помогают различать виды цветков. Для этого я использовал функцию information.gain() из пакета FSelector (Листинг 10).
Листинг 10. Определение важности признаков iris.
weights <- information.gain(Species ~ ., data = iris)
print(weights)
После выполнения программы появилась таблица со значениями важности для четырёх признаков. Чем больше значение, тем полезнее данный признак для определения вида цветка. Наибольшую важность имеют Petal.Width и Petal.Length, то есть ширина и длина лепестка (Рисунок 8).
## attr_importance
## Sepal.Length 0.4521286
## Sepal.Width 0.2672750
## Petal.Length 0.9402853
## Petal.Width 0.9554360
Рисунок 8 - Результат определения важности признаков набора iris.
Вывод по второму заданию: Я научился оценивать важность признаков с помощью FSelector. По полученной таблице можно определить, какие характеристики лучше использовать для классификации цветков. Наиболее полезными оказались характеристики лепестков.
В третьем задании мне нужно было преобразовать числовые значения в категории. Для этого я взял столбец Sepal.Length из набора iris, в котором записана длина чашелистика.
Я применил четыре метода дискретизации: interval, frequency, cluster и fixed. Во всех случаях я разделял данные на три категории (Листинг 11).
Листинг 11. Дискретизация данных четырьмя методами.
v <- iris$Sepal.Length
interval <- discretize(v, method = "interval", breaks = 3)
frequency <- discretize(v, method = "frequency", breaks = 3)
set.seed(123)
cluster <- discretize(v, method = "cluster", breaks = 3)
fixed <- discretize(v, method = "fixed", breaks = c(4, 5.5, 6.5, 8))
Метод interval делит диапазон значений на интервалы одинаковой ширины. Метод frequency старается распределить наблюдения по категориям примерно поровну. Метод cluster объединяет похожие значения, а в методе fixed я сам указал границы категорий.
Чтобы сравнить полученные результаты, я вывел количество элементов в каждой категории (Листинг 12).
Листинг 12. Сравнение результатов дискретизации.
table(interval)
table(frequency)
table(cluster)
table(fixed)
После выполнения программы появилась таблица со значениями важности для четырёх признаков. Чем больше значение, тем полезнее данный признак для определения вида цветка. Наибольшую важность имеют Petal.Width и Petal.Length, то есть ширина и длина лепестка (Рисунок 9).
## interval:
## interval
## [4.3,5.5) [5.5,6.7) [6.7,7.9]
## 52 70 28
##
## frequency:
## frequency
## [4.3,5.4) [5.4,6.3) [6.3,7.9]
## 46 53 51
##
## cluster:
## cluster
## [4.3,5.37) [5.37,6.36) [6.36,7.9]
## 46 62 42
##
## fixed:
## fixed
## [4,5.5) [5.5,6.5) [6.5,8]
## 52 63 35
Рисунок 9 - Сравнение результатов четырёх методов дискретизации.
Вывод по третьему заданию: Я научился превращать числовые значения в категории с помощью discretize(). Разные способы дают разные результаты, потому что используют разные правила выбора границ. Если нужно самостоятельно задать границы, подходит метод fixed.
В четвёртом задании я использовал набор данных Ozone из пакета mlbench. В нём находятся данные о загрязнении воздуха озоном и погодных условиях. Мне нужно было определить, какие показатели полезны для предсказания уровня озона.
Сначала я загрузил данные и удалил строки, в которых есть пропущенные значения. Это сделал с помощью функции na.omit(). (Листинг 13).
Листинг 13. Загрузка и подготовка данных Ozone.
data("Ozone", package = "mlbench")
ozone <- na.omit(Ozone)
head(ozone)
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13
## 5 1 5 1 5 5760 3 51 54 45.32 1450 25 57.02 60
## 6 1 6 2 6 5720 4 69 35 49.64 1568 15 53.78 60
## 7 1 7 3 4 5790 6 19 45 46.40 2631 -33 54.14 100
## 8 1 8 4 4 5790 3 25 55 52.70 554 -28 64.76 250
## 9 1 9 5 6 5700 3 73 41 48.02 2083 23 52.52 120
## 12 1 12 1 6 5720 3 44 51 54.32 111 9 63.14 150
dim(ozone)
## [1] 203 13
После этого я запустил алгоритм Boruta. В качестве переменной, которую нужно предсказать, использовал V4 - показатель уровня озона. Остальные столбцы использовались как признаки (Листинг 14).
Листинг 14. Определение важных признаков с помощью Boruta.
set.seed(123)
result <- Boruta(V4 ~ ., data = ozone, doTrace = 0)
print(result)
getSelectedAttributes(result)
После выполнения программы я получил результат отбора признаков. Алгоритм проверяет признаки и определяет, какие из них можно считать важными для предсказания. Названия подтверждённых признаков выводятся командой getSelectedAttributes() (Рисунок 10).
## Boruta performed 18 iterations in 0.2186453 secs.
## 9 attributes confirmed important: V1, V10, V11, V12, V13 and 4 more;
## 3 attributes confirmed unimportant: V2, V3, V6;
## [1] "V1" "V5" "V7" "V8" "V9" "V10" "V11" "V12" "V13"
Рисунок 10 - Результаты отбора признаков с помощью Boruta.
Далее я построил график типа boxplot, чтобы наглядно показать важность признаков (Листинг 15).
Листинг 15. Построение графика Boruta и сохранение в JPG.
jpeg("05_boruta.jpg", width = 1600, height = 1100)
plot(result)
dev.off()
На полученном графике показано сравнение признаков по их важности. По нему можно увидеть, какие признаки алгоритм считает полезными, а какие менее значимыми (Рисунок 11).
Рисунок 11 - Диаграмма важности признаков, полученная с помощью Boruta.
Вывод по четвёртому заданию: Я познакомился с алгоритмом Boruta и научился отбирать важные признаки для набора данных Ozone. Результаты отбора я вывел в консоль и показал на графике boxplot. Такой способ помогает понять, какие данные стоит учитывать при построении модели.
После выполнения всех заданий я сохранил содержимое консоли RStudio с командами и результатами в файл LR2_results.txt. Затем через меню Session → «Save Workspace As» сохранил рабочую среду в файл LR2.RData. Далее я создал новый документ R Markdown. После этого открыл его в RStudio заполнил и отредактировал в соответствии с требованиями и форматом из текущего отчета, добавив блоки для вывода выполнения команд и графиков. После чего сохранил и нажал кнопку Knit. Программа выполнила код и создала HTML-документ с текстом и графиками. После этого нажал «Publish» он до установил еще компонентов, выбрал RPubs он долго «крутил» но так и не загрузил мой результат. Видимо по тем же причинам что и не устанавливались компоненты некоторые…
Во время лабораторной работы я изучил основные возможности пакетов CARET, FSelector, arules и Boruta. Я научился создавать наборы данных, строить графики, сохранять их в JPG, определять важность признаков и переводить числовые значения в категории. Также я познакомился с R Markdown и способом публикации отчётов на RPubs. Все задания помогли мне лучше понять, как можно использовать язык R для анализа данных.