Материал из Linuxformat.

R Свободный инструментарий для статистической обработки данных\

Содержание 1 Интеллектуальный анализ, или Data Mining 1.1 Графический анализ 1.2 Сокращение размерности 1.3 Классификация без обучения 1.4 Кластерный анализ 1.5 Классификация с обучением

Интеллектуальный анализ, или Data Mining

R
Часть 1: Введение в R Часть 2: Данные и графики Часть 3: Мир с точки зрения R Часть 4: Анализ данных с R Часть 5: Работа с двумя переменными Часть 6: Взаимосвязь случайных величин Часть 7: Графические интерфейсы Часть 8: Интеллектуальный анализ, или Data Mining

R на примере
R: Размножаем реальность R: GUI на примере R: Ускоряем на примере R: Поиск закономерностей

ЧАСТЬ 4 Бытует мнение, что сбор невинных данных о вас и их обработка – это дело Врагов Свободы. Так ли эт, каждый решает сам – важно, что мы можем ударить по ним их же собственным оружием. Алексей Шипунов и Евгений Балдин откроют вам арсенал.

Выражение «data mining» всё чаще и чаще встречается на обложках книг по анализу данных, не говоря уж о росте упоминаний об этом методе во всеядном Интернете. Некоторые даже полагают, что эпоха статистики одной-двух переменных закончилась, и наступило новое время – время интеллектуального анализа больших и сверхбольших массивов данных. На самом деле, под «дата-майнингом» подразумеваются любые визуальные или аналитические методы, позволяющие «нащупать» структуру в большом объеме информации – в ином случае ее можно обнаружить и без специальных методов, например, просто поглядев на распределение единственного изучаемого параметра.

Для анализа используются, как правило, многомерные данные, т.е. такие, которые можно представить в виде таблицы из нескольких колонок-переменных. Поэтому более традиционным названием для этих методов является «многомерный анализ» или «многомерная статистика»; но «data mining» звучит, конечно, серьёзнее. Кроме многомерности и большого размера (сотни, а то и тысячи строк и столбцов), используемые данные отличаются ещё и тем, что переменные в них могут быть совершенно разных типов (качественные, балльные, счётные, непрерывные), причём даже «типичные» для статистики непрерывные числовые переменные вполне могут не подчи- няться заранее известным законам распределения, т.е. могут не быть параметрическими.

Сами по себе многомерные методы делятся на визуализационные и методы классификации с обучением. В первом случае результат можно анализировать в основном зрительно, а во втором возможна статистическая проверка результатов. Разумеется, граница между этими группами не резкая, но для удобства мы будем рассматривать их именно в таком порядке.

Графический анализ

Самое простое, что можно сделать с многомерными данными – это построить график. Разумеется, предварительно надо свести всё их разнообразие к двум или, в крайнем случае, трём измерениям. Эта операция называется «сокращением размерности». Если же переменных уже три и все они непрерывные, то с представлением данных замечательно справится RGL, написанный с использованием OpenGL и поддерживающий трёхмерную визуализацию.

Давайте визуализируем при помощи пакета RGL четыре из пяти колонок массива iris (см. врезку):

> library(rgl)
> plot3d(iris$Sepal.Length, iris$Sepal.Width,
+       iris$Petal.Length, col=iris$Species, size=3)

Размер появившегося окна и проекцию отображения данных можно и нужно менять c помощью мыши. Сразу видно, что один из видов (Iris setosa) хорошо отличается от двух других по признаку длины лепестков (Petal.Length). Кстати, визуализированных признаков действительно четыре: один из них (вид ириса) закодирован цветом.

Для трёхмерной визуализации данных можно обойтись и без OpenGL, взяв, например, пакет scatterplot3d. Но, по-хорошему, трёхмерными графиками лучше не злоупотреблять. Очень часто двумерные проекции трёхмерных объектов не раскрывают, а «затемняют» суть явления. Правда, в случае RGL это компенсируется возможностью свободно менять «точку обзора».

Есть и более специализированные системы для визуализации многомерных данных без снижения размерности. Среди них можно отметить пакет rggobi, основанный на системе Ggobi. Сегодня мы его рассматривать не будем, поскольку (строго говоря) использование «стороннего ПО» уже выводит нас за рамки R.

Ещё один способ визуализации многомерных данных – это построение графиков-таблиц. Здесь R обладает колоссальными возможностями, которые предоставляются пакетом lattice, предназначенным для так называемой Trellis-графики. Последняя долгое время была общепризнанной изюминкой S-PLUS, но теперь стала доступна и в R. Вот как можно визуализировать четыре признака ирисов:

> library(lattice)
> xyplot(Sepal.Length ~ Petal.Length + Petal.Width | Species,
+       data = iris, auto.key=TRUE)

Сокращение размерности

Вернёмся к методам сокращения размерности. Самый распространённый из них – это «анализ главных компонент». Суть его заключается в том, что все признаки-колонки преобразуются в компоненты, причём наибольшую информацию о разнообразии объектов несёт первая компонента, вторая – чуть меньше, и так по убыванию.

Таким образом, хотя компонент получается столько же, сколько изначальных признаков, почти все нужные нам сведения сосредоточены в первых двух-трёх из них, поэтому их можно использовать для визуализации данных на плоскости. Обычно используется первая и вторая (реже первая и третья) компоненты. Компоненты часто называют «факторами», и это порождает некоторую путаницу с похожим на анализ главных компонент «факторным анализом», преследующим, однако, совсем иные цели.

Вот как делается анализ главных компонент на наших данных про ирисы:

> iris.pca <- princomp(scale(iris[,1:4]))

Мы использовали функцию scale() для того, чтобы привести все четыре переменные к одному масштабу (по умолчанию она вычитает из данных среднее и делит их на его квадрат), поскольку переменные, варьирующие в разных масштабах, способны исказить результат анализа. В R реализован и другой метод анализа главных компонент, основанный на иных преобразованиях матрицы: он вызывается функцией prcomp().

Выведем служебный график, демонстрирующий относительные вклады каждого компонента в общий разброс данных:

 > plot(iris.pca)

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

> summary(iris.pca)
Importance of components:
                       Comp.1    Comp.2 Comp.3        Comp.4
Standard deviation     1.7026571 0.9528572 0.38180950 0.143445939
Proportion of Variance 0.7296245 0.2285076 0.03668922 0.005178709
Cumulative Proportion 0.7296245 0.9581321 0.99482129 1.000000000

Теперь перейдём собственно к визуализации:

> iris.p <- predict(iris.pca)
> plot(iris.p[,1:2], type="n",xlab="PC1", ylab="PC2")
> text(iris.p[,1:2],
+    labels=abbreviate(iris[,5],1,method="both.sides"))

Как видно из рисунка, Iris setosa (отмечен буквой s) сильно отличается от двух остальных видов, Iris versicolor (v) и Iris virginica (a). Функция predict() позволяет расположить исходные случаи (строки) в пространстве вновь найденных компонент. Функция abbreviate() «умным» образом сокращает названия до одной буквы.

Пакеты ade4 и vegan' реализуют множество вариаций анализа главных компонент, но самое главное то, что они содержат гораздо больше возможностей для визуализации. Например, так можно проанализировать те же данные по ирисам с ade4:

> install.packages("ade4", dependencies=TRUE)
> library(ade4)
> iris.dudi <- dudi.pca(iris[,1:4], scannf=FALSE)
> s.class(iris.dudi$li, iris[,5])

Нетрудно заметить, что различия между ирисами на этом графике видны яснее. Более того, с помощью ade4 можно проверить качество разрешения между классами (в данном случае – видами ирисов):

> iris.between <- between(iris.dudi, iris[,5], scannf=FALSE)
> randtest(iris.between)
Monte-Carlo test
Call: randtest.between(xtest = iris.between)
Observation: 0.7224358
Based on 999 replicates
Simulated p-value: 0.001
Alternative hypothesis: greater
    Std.Obs      Expectation Variance
6.868114e+01 1.374496e-02 1.064728e-04

Из распечатки видно, что классы различаются хорошо (0,7 близко к 1) и стабильно. Этот метод уже действительно ближе к «настоящей статистике», нежели к типичной визуализации данных.

Классификация без обучения

Ещё одним способом снижения размерности является ординация (упорядочение или классификация без обучения), проводимая на основании заранее вычисленных значений сходства между всеми парами объектов (строк). В результате этой процедуры получается квадратная матрица расстояний, диагональ которой обычно состоит из нулей (расстояние между объектом и им же самим равно нулю). За десятилетия развития этой области статистики придуманы сотни коэффициентов сходства, из которых наиболее употребимыми являются евклидово и квартальное (манхэттеновское). Эти коэффициенты применимы в основном к непрерывным переменным. Балльные и бинарные переменные в общем случае требуют других коэффициентов, но в пакете cluster реализована функция daisy(), которая способна распознавать тип переменной и применять соответствующие коэффициенты, а в пакете vegan доступно множество дополнительных коэффициентов сходства.

Вот как можно построить матрицу сходства (лучше её всё-таки называть матрицей различий, поскольку в её ячейках стоят именно расстояния) для наших ирисов:

> library(cluster)
> iris.dist <- daisy(iris[,1:4], metric="manhattan")

С полученной таким образом матрицей можно делать довольно многое. Одним из самых простых её применений является «многомерное шкалирование» или, иначе, «анализ главных координат» (это название применяют в основном для метрических вариантов данного метода).

Попробуем разобраться в сути «многомерного шкалирования». Допустим, в результате долгой и изнурительной последовательности действий над картой были измерены расстояния между десятком городов, результаты сохранены, а сама карта была неожиданно потеряна. Теперь перед исследователями стоит задача: восстановить карту взаимного расположения городов, зная только расстояния между ними. Именно её и решает многомерное шкалирование, причём это отнюдь не метафора. Можно набрать example(cmdscale) и на примере 21 европейского города посмотреть, как подобное вычисляется на самом деле. Для наших же ирисов многомерное шкалирование можно применить так:

> iris.c <- cmdscale(iris.dist)
> plot(iris.c[,1:2], type="n", xlab="Dim. 1", ylab="Dim. 2")
> text(iris.c[,1:2],
+      labels=abbreviate(iris[,5],1, method="both.sides"))

Как видно из графика, результат очень напоминает полученный при анализе главных компонент, что неудивительно, так как внутренняя структура данных (которую нам и надо найти в процессе «data mining») не изменилась. Кроме cmdscale(), советуем обратить внимание на непараметрический вариант этого метода, осуществляемый при помощи функции isoMDS().

Кластерный анализ

Ещё одним вариантом работы с матрицей различий является «кластерный анализ». Существует множество его разновидностей, причём наиболее употребимыми являются иерархические методы, которые вместо уже привычных нам двумерных графиков производят «полуторамерные» деревья классификации или дендрограммы.

> iriss <- iris[seq(1,nrow(iris),5),]
> iriss.dist <- daisy(iriss[,1:4])
> iriss.h <- hclust(iriss.dist, method="ward")
> plot(iriss.h,
+     labels=abbreviate(iriss[,5],1, method="both.sides"))

Для построения «дерева» была взята только каждая пятая строка данных, иначе ветки сидели бы слишком плотно. Это, кстати, недостаток иерархической кластеризации как метода визуализации данных. Метод Уорда [Ward] даёт очень хорошо очерченные кластеры – при условии, естественно, что их удаётся найти; поэтому неудивительно, что в нашем случае все три вида разделились. При этом отлично видно, что виды на «v» (versicilor и virginica) разделяются на более низком уровне (Height ~ 12), то есть сходство между ними сильнее, чем каждого из них с третьим видом.

Кластерный анализ этого типа весьма привлекателен тем, что даёт готовую классификацию. Однако не стоит забывать, что это «всего лишь визуализация». Насколько «хороши» получившиеся кластеры, проверить порой непросто, хотя и здесь существует множество методов. Один из них, так называемый «silhouette plot», реализован в пакете cluster. Чтобы увидеть его в действии, достаточно набрать example(agnes). Ещё один, очень модный сейчас метод, основанный на boostrap-репликации, реализован в пакете pvclust:

> install.packages("pvclust", dependencies=TRUE)
> library(pvclust)
> iriss.pv <- pvclust(t(iriss[,1:4]),
+ method.dist="manhattan", method.hclust="ward", nboot=100)
> plot(iriss.pv, print.num=FALSE)

Красным цветом на графике печатаются p-значения (p-values), связанные с устойчивостью кластеров в процессе репликации исходных данных. Величины, близкие к 100, считаются «хорошими». В данном случае мы видим как раз «хорошую» устойчивость получившихся основных трёх кластеров (разных видов ирисов), и неплохую устойчивость кластера, состоящего из двух видов на «v».

Кроме иерархических методов кластеризации, существуют и другие. Из них наиболее интересны так называемые fuzzy-методы, основанные на той идее, что каждый объект может принадлежать к нескольким кластерам сразу, но с разной «силой». Вот как реализуется такой метод в пакете cluster:

> library(cluster)
> iris.f <- fanny(iris[,1:4], 3)
> plot(iris.f, which=1)
> head(data.frame(sp=iris[,5], iris.f$membership))
    sp              X1                X2           X3
1 setosa             0.9142273         0.03603116   0.04974153
2 setosa            0.8594576         0.05854637   0.08199602
3 setosa            0.8700857         0.05463714   0.07527719
4 setosa            0.8426296         0.06555926   0.09181118
5 setosa            0.9044503         0.04025288   0.05529687
6 setosa            0.7680227         0.09717445   0.13480286

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

Подобный метод реализован также и в пакете e1071. Функция называется cmeans(), но в этом случае вместо количества кластеров можно указать предполагаемые центры, вокруг которых будут группироваться элементы.

Классификация с обучением

Теперь обратимся к методам, которые могут называться «визуализацией» лишь частично. В зарубежной литературе именно их принято называть «методами классификации». Чтобы с ними работать, надо освоить технику «обучения». Как правило, выбирается часть данных с известной групповой принадлежностью. На основании анализа этой части, называемой «тренировочной выборкой», строится гипотеза о том, как должны распределяться по группам остальные, неклассифицированные данные. При этом обычно можно узнать, насколько хорошо работает та или иная гипотеза. Кроме того, методы классификации с обучением можно с успехом применять и для других целей: например, для выяснения важности признаков.

Один из самых простых методов в данной группе – это «линейный дискриминантный анализ». Его основной идеей является создание функций, которые на основании линейных комбинаций значений признаков (это и есть классификационная гипотеза) «сообщают», куда нужно отнести данную особь. Воспользуемся им для выяснения структуры данных ирисов:

> iris.train <- iris[seq(1,nrow(iris),5),]
> iris.unknown <- iris[-seq(1,nrow(iris),5),]
> library(lda.cv)
> iris.lda <- lda(scale(iris.train[,1:4]), iris.train[,5])
> iris.ldap <- predict(iris.lda, iris.unknown[,1:4])$class
> table(iris.ldap, iris.unknown[,5])
iris.ldap             setosa           versicolor         virginica
  setosa                 0                  0                 0
  versicolor            34                  0                 0
  virginica              6                 40                40

На выходе получился довольно забавный результат: наша тренировочная выборка привела к построению гипотезы, по которой все virginica и versicolor (а также часть setosa) попали в одну группу. Это говорит не только о близости видов на «v», но также и о недостаточной «тренировке».

Великое множество методов «data mining», разумеется, нереально охватить в небольшой статье. Однако нельзя и не упомянуть ещё об одном современном методе, основанном на идее вычисления параметров гиперплоскости, разделяющей различные группы в многомерном пространстве признаков – Support Vector Machines.

> library(e1071)
> iris.svm <- svm(Species ~ ., data = iris.train)
> iris.svmp <- predict(iris.svm, iris[,1:4])
> table(iris.svmp, iris[,5])
iris.svmp              setosa           versicolor            virginica
  setosa                 50                  0                    0
  versicolor              0                 50                    14
  virginica               0                  0                    36

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

На этом мы завершаем наше «углубленное введение» в мир R вообще и в технологии интеллектуального анализа – в частности. В заключение хотелось бы отметить, что настоящая статья ни в коем случае не заменяет больших книг и справочников, написанных по теме «data mining». Мы ставили целью лишь дать примерное представление о многообразии методов анализа многомерных данных и наиболее распространенных способах решения основных проблем, возникающих при поиске порядка в больших массивах информации. LXF