А. Б. Шипунов, Е. М. Балдин, П. А. Волкова, А. И. Коробейников, С. А. Назарова
 Скачать 3.04 Mb.
|
Кроме того, можно проверить качество разрешения между классами (в данном случае видами ирисов): > iris.between <- bca(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 Классификация без обучения, или Кластерный анализ 157 Классы (виды ирисов) различаются хорошо. Об этом говорит осно- ванное на 999 повторениях (со слегка различающимися параметрами) анализа значение Observation. Если бы это значение было меньше 0.5 (то есть 50%), то нам пришлось бы говорить о нечетких различиях. Как видно, использованный метод уже ближе к «настоящей статисти- ке», нежели к типичной визуализации данных. 6.3. Классификация без обучения, или Кластерный анализ Другим способом снижения размерности является классификация без обучения (упорядочение, или ординация), проводимая на основании заранее вычисленных значений сходства между всеми парами объектов (строк). В результате этой процедуры получается квадратная матрица расстояний, диагональ которой обычно составлена нулями (ведь рас- стояние между объектом и им же самим равно нулю). За десятилетия развития этой области статистики придуманы сотни коэффициентов сходства, из которых наиболее употребительными являются эвклидово и квартальное (манхеттеновское), применимые в основном к непрерыв- ным переменным (объяснение см. на рис. 41). Коэффициент корреляции тоже может быть мерой сходства. Балльные и бинарные переменные в общем случае требуют других коэффициентов, но в пакете cluster реа- лизована функция daisy(), способная распознавать тип переменной и применять соответствующие коэффициенты, а в пакете vegan реализо- вано множество дополнительных коэффициентов сходства. Рис. 41. Эвклидово (1) и манхеттеновское (2) расстояния Вот как можно построить матрицу сходства (лучше ее все-таки на- зывать матрицей различий, поскольку в ее ячейках стоят именно рас- стояния) для наших ирисов: > library(cluster) > iris.dist <- daisy(iris[,1:4], metric="manhattan") 158 Анализ структуры: data mining (Указанное нами в качестве аргумента манхеттеновское расстояние будет вычислено только для интервальных данных, для данных других типов функция daisy() самостоятельно рассчитает более универсаль- ное расстояние Говера, «Gower’s distance»). С полученной матрицей можно делать самые разные вещи. Одно из самых простых применений — сделать многомерное шкалирование (или, как его еще иногда называют, «анализ главных координат»). Суть метода можно объяснить так. Допустим, мы измерили по карте расстоя- ние между десятком городов, а карту потеряли. Задача — восстановить карту взаимного расположения городов, зная только расстояния между ними. Такую же задачу решает многомерное шкалирование. Причем это не метафора — вы можете запустить example(cmdscale) и посмотреть (на примере 21 европейского города), как это происходит на самом деле. А вот для наших ирисов многомерное шкалирование можно применить так (рис. 42): > 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(). Результат многомерного шкалирования тоже можно украсить. Вот как мы попытались выделить получившиеся группы (рис. 43): > library(KernSmooth) > est <- bkde2D(iris.c[,1:2], bandwidth=c(.7,1.5)) > 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")) > contour(est$x1, est$x2, est$fhat, add=TRUE, + drawlabels=FALSE, lty=3) Функция bkde2D() из пакета KernSmooth предсказывает для каж- дой точки графика плотность распределения данных (в нашем случае — букв), а функция contour() рисует нечто вроде карты, где пиками явля- ются места с максимальной плотностью букв. Чтобы график выглядел красиво, параметры bkde2d() были подобраны вручную. Другой вариант работы с матрицей различий — это кластерный ана- лиз. Существует множество его разновидностей. Наиболее употреби- Классификация без обучения, или Кластерный анализ 159 −4 −2 0 2 4 6 − 3 − 2 − 1 0 1 2 Dim. 1 D im . 2 s s s s s s s s s s s s s s s s s s s s s s s s s s s ss s s s s s s s s s s s s s s s s s s s s s v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a aa a a a a a a a a a Рис. 42. Результат многомерного шкалирования данных об ирисах тельными являются иерархические методы, которые вместо уже при- вычных нам двумерных графиков производят «полуторамерные» де- ревья классификации, или дендрограммы. Вот как это делается (рис. 44): > iriss <- iris[seq(1,nrow(iris),5),] > iriss.dist <- daisy(iriss[,1:4]) > iriss.h <- hclust(iriss.dist, method="ward.D") > plot(iriss.h, labels=abbreviate(iriss[,5],1, + method="both.sides"), main="") Поскольку на выходе получается «дерево», мы взяли для его по- строения каждую пятую строку данных, иначе ветки сидели бы слиш- ком плотно (это, кстати, недостаток иерархической кластеризации как метода визуализации данных). Метод Уорда (ward.D) дает очень хоро- шо очерченные кластеры (конечно, если их удается найти), и поэтому неудивительно, что в нашем случае все три вида разделились, причем отлично видно, что виды на «v» (Iris versicilor и I. virginica) разделяют- 160 Анализ структуры: data mining −4 −2 0 2 4 6 − 3 − 2 − 1 0 1 2 Dim. 1 D im . 2 s s s s s s s s s s s s s s s s s s s s s s s s s s s ss s s s s s s s s s s s s s s s s s s s s s v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a aa a a a a a a a a a Рис. 43. Визуализация результата многомерного шкалирования при по- мощи функции плотности ся на более низком уровне (Height≈ 12), то есть сходство между ними сильнее, чем каждого из них с третьим видом. Задача . Попробуйте построить иерархическую классификацию сор- тов пива. Данные содержатся в файле pivo-s.txt, а расшифровка при- знаков — в файле pivo-c.txt. Сведения собирались в 2001 году в Москве (это надо учесть, если ваш опыт говорит о чем-то ином). Кластерный анализ этого типа весьма привлекателен тем, что дает готовую классификацию. Однако не ст 0 оит забывать, что это — все- го лишь визуализация, и кластерный анализ может «навязать» дан- ным структуру, которой у них на самом деле нет. Насколько «хороши» получившиеся кластеры, проверить порой непросто, хотя и здесь со- здано множество методов. Один метод реализует пакет cluster — это так называемый «silhouette plot» (за примером можно обратиться к example(agnes) ). Другой, очень «модный» метод, основанный на так называемой bootstrap-репликации, реализует пакет pvclust (рис. 45): > library(pvclust) Классификация без обучения, или Кластерный анализ 161 s s s s s s s s s s v v v v v v v v v v a a a a a a a a a a 0 10 20 30 40 Рис. 44. Дендрограмма, отражающая иерархическую классификацию ирисов > irisst <- t(iriss[,1:4]) > colnames(irisst) <- paste(abbreviate(iriss[,5], 3), + colnames(irisst)) > iriss.pv <- pvclust(irisst, method.dist="manhattan", + method.hclust="ward.D", nboot=100) Bootstrap (r = 0.5)... Done. Bootstrap (r = 0.5)... Done. Bootstrap (r = 0.5)... Done. Bootstrap (r = 0.75)... Done. Bootstrap (r = 0.75)... Done. > plot(iriss.pv, col.pv=c(au=0, bp="darkgreen", edge=0), main="") Над каждым узлом печатаются p-значения (au), связанные с устой- чивостью кластеров в процессе репликации исходных данных. Значе- ния, близкие к 100, считаются «хорошими». В нашем случае мы видим как раз «хорошую» устойчивость получившихся основных трех класте- 162 Анализ структуры: data mining st s 46 st s 26 st s 31 st s 36 st s 1 st s 41 st s 16 st s 6 st s 11 st s 21 vr g 12 6 vr g 13 1 vr g 10 6 vr g 13 6 vr g 14 6 vr g 11 1 vr g 11 6 vr g 10 1 vr g 12 1 vr g 14 1 vr s 71 vr s 86 vr s 51 vr s 66 vr s 76 vr s 96 vr s 56 vr s 91 vr s 61 vr s 81 0 10 20 30 40 50 60 70 84 61 61 31 43 58 54 42 24 47 44 42 61 54 53 59 66 20 46 43 48 45 20 38 90 46 46 77 bp Рис. 45. Дендрограмма со значениями поддержки кластеров (пакет pvclust ) ров (разных видов ирисов) и неплохую устойчивость кластера, состоя- щего из двух видов на «v». Кроме иерархических методов кластеризации, существуют и другие. Есть, например, «методы средних», которые стараются получить зара- нее заданное количество кластеров. Приведем пример. В файле eq.txt записаны данные измерений разных частей растений у двух сильно от- личающихся видов хвощей. Попробуем проверить, сможет ли такой ме- тод распознать эти виды, зная лишь признаки. Используем функцию kmeans() : > eq <- read.table("data/eq.txt", h=TRUE) > eq.k <- kmeans(eq[,-1], centers=2) > table(eq.k$cluster, eq$SPECIES) arvense fluviatile 1 37 5 2 1 41 Классификация без обучения, или Кластерный анализ 163 Получилось! Ошибка классификации довольно низкая. Интересно будет потом посмотреть, на каких именно растениях метод ошибся и почему. Отметьте, кстати, что функция kmeans() (как и следующие две) берет на входе не матрицу расстояний, а оригинальные данные. Очень интересны так называемые нечеткие методы кластеризации, основанные на идее того, что каждый объект может принадлежать к нескольким кластерам сразу — но с разной «силой». Вот как реализу- ется такой метод в пакете cluster (рис. 46): > iris.f <- fanny(iris[,1:4], 3) > plot(iris.f, which=1, main="") > 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(), но в этом случае вместо количества кластеров можно указать предполагаемые центры, вокруг которых будут группи- роваться элементы. * * * Как вы уже поняли, методы классификации без обучения могут ра- ботать не только с интервальными и шкальными, но и с номинальными данными. Для этого надо лишь закодировать номинальные признаки в виде нулей и единиц. Дальше умные функции типа daisy() сами выбе- рут коэффициент сходства, а можно и указать его вручную (например, так: dist(..., method="binary"). А есть ли многомерные методы, ко- торые могут работать с таблицами сопряженности? Оказывается, такие не просто существуют, но широко используются в разных областях, на- пример в экологии. 164 Анализ структуры: data mining −3 −2 −1 0 1 2 3 − 3 − 2 − 1 0 1 2 3 Рис. 46. Результат кластеризации функцией fanny() Анализ связей (correspondence analysis) — один из них. Как и любой многомерный метод, он позволяет визуализировать структуру данных, причем работает он через составление таблиц сопряженности. Простой вариант анализа связей реализован в пакете MASS функцией corresp() (рис. 47): > library(MASS) > caith.ru <- caith > row.names(caith.ru) <- abbreviate(c("голубоглазые", + "сероглазые","кареглазые","черноглазые"), 10, method="both") > names(caith.ru) <- abbreviate(c("блондины","рыжие", + "русоволосые","шатены","брюнеты"), 10, method="both") > caith.ru блонд рыжие русов шатен брюне голуб 326 38 241 110 3 серог 688 116 584 188 4 карег 343 84 909 412 26 черно 98 48 403 681 85 Классификация без обучения, или Кластерный анализ 165 > biplot(corresp(caith.ru, nf = 2)) −0.4 −0.2 0.0 0.2 0.4 0.6 − 0. 4 − 0. 2 0. 0 0. 2 0. 4 0. 6 голубоглаз сероглазые кареглазые черноглазы −0.5 0.0 0.5 1.0 − 0. 5 0. 0 0. 5 1. 0 блондины рыжие русоволосы шатены брюнеты Рис. 47. График взаимосвязей, полученный из таблицы сопряженности Данные представляют собой статистическую таблицу, созданную в результате переписи населения, где отмечались, среди прочего, цвета глаз и волос. На получившемся графике удалось визуализовать обе пе- ременные, причем чем чаще две признака (скажем, голубоглазость и светловолосость) встречаются вместе у одного человека, тем ближе на графике эти надписи. Эта возможность (визуализация одновременно нескольких наборов признаков) широко используется в экологии при анализе сообществ жи- вотных и растений, обитающих на одной территории. Если есть, напри- мер, данные по 10 озерам, где указаны их названия и то, какие виды рыб в них обитают, то можно будет построить график, где названия озер бу- дут расположены ближе к тем названиям рыб, которые для этих самых озер более характерны. Более того, если в данные добавить еще и абио- тические признаки озер (скажем, максимальную глубину, прозрачность воды, придонную температуру летом), то можно будет визуализовать 166 Анализ структуры: data mining сразу три набора признаков! Пакет vegan содержит несколько таких функций, например cca() и decorana(). 6.4. Классификация с обучением, или Дискриминантный анализ Перейдем теперь к методам, которые лишь частично могут назы- ваться «визуализацией». В зарубежной литературе именно их приня- то называть «методами классификации». Для того чтобы работать с методами классификации с обучением, надо сначала освоить технику «обучения». Как правило, выбирается часть данных с известной груп- повой принадлежностью. На основании анализа этой части (трениро- вочной выборки) строится гипотеза о том, как должны распределять- ся по группам остальные, неклассифицированные данные. В результате анализа, как правило, можно узнать, насколько хорошо работает та или иная гипотеза. Кроме того, методы классификации с обучением мож- но с успехом применять и для других целей, например для выяснения важности признаков для классификации. Один из самых простых методов в этой группе — линейный дискри- минантный анализ (linear discriminant analysis). Его основная идея — создание функций, которые на основании линейных комбинаций значе- ний признаков (это и есть классификационная гипотеза) «сообщают», куда нужно отнести данную особь. Вот как можно применить такой анализ: > library(MASS) > iris.train <- iris[seq(1,nrow(iris),5),] > iris.unknown <- iris[-seq(1,nrow(iris),5),] > iris.lda <- lda(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 40 0 0 versicolor 0 40 7 virginica 0 0 33 Наша тренировочная выборка привела к построению гипотезы, по которой все виды (за исключением части virginica) попали в «свою» группу. Заметьте, что линейный дискриминантный анализ не требу- ет шкалирования (функция scale()) переменных. Можно более точно подсчитать ошибки классификации. Для этого мы написали функцию misclass() : Классификация с обучением, или Дискриминантный анализ 167 > misclass <- function(pred, obs) { + tbl <- table(pred, obs) + sum <- colSums(tbl) + dia <- diag(tbl) + msc <- (sum - dia)/sum * 100 + m.m <- mean(msc) + cat("Classification table:", "\n") + print(tbl) + cat("Misclassification errors:", "\n") + print(round(msc, 1)) + } > misclass(iris.ldap, iris.unknown[,5]) Classification table: obs pred setosa versicolor virginica setosa 40 0 0 versicolor 0 40 7 virginica 0 0 33 Misclassification errors: setosa versicolor virginica 0.0 0.0 17.5 Как видим, ошибка классификации — 17.5%. Результат дискриминантного анализа можно проверить статистиче- ски. Для этого используется многомерный дисперсионный анализ (MA- NOVA), который позволяет выяснить соответствие между данными и моделью (классификацией, полученной при помощи дискриминацион- ного анализа): > ldam <- manova(as.matrix(iris.unknown[,1:4]) iris.ldap) > summary(ldam, test="Wilks") Df Wilks approx F num Df den Df Pr(>F) iris.ldap 2 0.026878 145.34 8 228 < 2.2e-16 *** Residuals 117 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Здесь имеет значение и значимость статистики Фишера (F) (см. гла- ву про анализ связей), и само значение «Wilks», поскольку это — «likeli- hood ratio», то есть отношение правдоподобия, в данном случае вероят- ность того, что группы не различаются. Чем ближе статистика Wilks к нулю, тем лучше классификация. В данном случае есть и достоверные различия между группами, и качественная классификация. 168 Анализ структуры: data mining Линейный дискриминантный анализ можно использовать и для ви- зуализации данных, например так (рис. 48): > iris.lda2 <- lda(iris[,1:4], iris[,5]) > iris.ldap2 <- predict(iris.lda2, dimen=2)$x > plot(iris.ldap2, type="n", xlab="LD1", ylab="LD2") > text(iris.ldap2, labels=abbreviate(iris[,5], 1, + method="both.sides")) −5 0 5 10 − 2 − 1 0 1 2 LD1 LD 2 s s s s s s s s s s s s s s s s s s s s s s s s s s s ss s s s s s s s s s s s s s s s s s s s s s v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v v vv v v a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a Рис. 48. Графическое представление результатов дискриминантного анализа Здесь мы в качестве тренировочной использовали всю выборку це- ликом. Как видно, виды на «v» разделились даже лучше, чем в преды- дущих примерах, поскольку дискриминантный анализ склонен часто переоценивать различия между группами. Это свойство, а также па- раметричность метода привели к тому, что в настоящее время этот тип анализа используется все реже. На замену дискриминантному анализу придумано множество мето- дов с похожим принципом работы. Один из самых оригинальных — Классификация с обучением, или Дискриминантный анализ 169 это так называемые «деревья классификации», или «деревья решений» («classification trees», или «decision trees»). Они позволяют выяснить, какие именно показатели могут быть использованы для разделения объектов на заранее заданные группы. В результате строится ключ, в котором на каждой ступени объекты делятся на две группы (рис. 49): > library(tree) > iris.tree <- tree(iris[,5] ., iris[,-5]) > plot(iris.tree) > text(iris.tree) | Petal.Length < 2.45 Petal.Width < 1.75 Petal.Length < 4.95 Sepal.Length < 5.15 Petal.Length < 4.95 setosa versicolor versicolor virginica virginica virginica Рис. 49. Дерево классификации для данных об ирисах Здесь мы опять использовали всю выборку в качестве тренировоч- ной (чтобы посмотреть на пример частичной тренировочной выборки, можно набрать ?predict.tree). Получился график, очень похожий на так называемые определительные таблицы, по которым биологи опре- деляют виды организмов. Из графика легко понять, что к setosa от- носятся все ирисы, у которых длина лепестков меньше 2.45 (читать 170 Анализ структуры: data mining график надо влево), а из оставшихся ирисов те, у которых ширина ле- пестков меньше 1.75, а длина — меньше 4.95, относятся к versicolor. Все остальные ирисы относятся к virginica. Деревья классификации реализованы и в другом пакете — rpart. Задача . Попробуйте выяснить, какие именно значения признаков различают два вида растений, которые мы использовали как пример к функции kmeans() (данные eq.txt, расшифровка признаков — в файле eq-c.txt ). Еще один, набирающий все большую популярность метод тоже идео- логически близок к деревьям классификации. Он называется «Random Forest», поскольку основой метода является производство большого ко- личества классификационных «деревьев». > library(randomForest) > set.seed(17) > iris.rf <- randomForest(iris.train[,5] ., + data=iris.train[,1:4]) > iris.rfp <- predict(iris.rf, iris.unknown[,1:4]) > table(iris.rfp, iris.unknown[,5]) iris.rfp setosa versicolor virginica setosa 40 0 0 versicolor 0 39 9 virginica 0 1 31 Заметна значительно более высокая эффективность этого метода по сравнению с линейным дискриминантным анализом. Кроме того, Random Forest позволяет выяснить значимость («importance») каждо- го признака, а также дистанции между всеми объектами тренировочной выборки («proximity»), которые затем можно использовать для класте- ризации или многомерного шкалирования. Наконец, этот метод позво- ляет производить «чистую визуализацию» данных, то есть он может работать как метод классификации без обучения (рис. 50): > set.seed(17) > iris.urf <- randomForest(iris[,-5]) > MDSplot(iris.urf, iris[,5]) Великое множество методов «data mining», разумеется, нельзя охва- тить в одной главе. Однако нельзя не упомянуть еще об одном совре- менном методе классификации с обучением, основанном на идее вычис- ления параметров гиперплоскости, разделяющей различные группы в многомерном пространстве признаков, «Support Vector Machines». Классификация с обучением, или Дискриминантный анализ 171 −0.4 −0.2 0.0 0.2 − 0. 4 − 0. 2 0. 0 0. 2 0. 4 Рис. 50. Визуализация данных при помощи «Random Forest» > 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 его можно использовать (как мы здесь видим) и для большего числа групп. Работает он эффективнее дискриминант- ного анализа, хотя и не так точно, как Random Forest. В заключение приведем краткий обзор применимости основных из обсуждавшихся в этой главе методов многомерного анализа данных (рис. 51). При чтении диаграммы стоит помнить о том, что исходные данные можно трансформировать (например, создать из таблицы дан- |