. А зависит ли вес от роста в наших искусственных данных? Если да, то какая эта зависимость?
* * *
Ответ к задаче про тест знаков
. Достаточно написать:
> a <- c(1, 2, 3, 4, 5, 6, 7, 8, 9)
> b <- c(5, 5, 5, 5, 5, 5, 5, 5, 5)
> dif <- a-b
> pol.dif <- dif[dif > 0]
> binom.test(length(pol.dif), length(dif))
Exact binomial test data:
length(pol.dif) and length(dif)
Если выборок больше двух
133
number of successes = 4, number of trials = 9, p-value = 1
alternative hypothesis:
true probability of success is not equal to 0.5 95 percent confidence interval:
0.1369957 0.7879915
sample estimates:
probability of success
0.4444444
То есть отличий от половины нет, и тест знаков считает такие вы- борки одинаковыми! Чего и можно было, в общем, ожидать — ведь мы уже писали, что проверяются не выборки как таковые, а лишь их цен- тральные значения.
Ответ к задаче про озон
. Чтобы узнать, нормально ли распреде- ление, применим созданную в предыдущей главе функцию normality3():
> ozone.month <- airquality[,c("Ozone","Month")]
> ozone.month.list <- unstack(ozone.month)
> normality3(ozone.month.list)
$‘5‘
[1] "NOT NORMAL"
$‘6‘
[1] "NORMAL"
$‘7‘
[1] "NORMAL"
$‘8‘
[1] "NORMAL"
$‘9‘
[1] "NOT NORMAL"
Здесь мы применили функцию unstack(), которая «разобрала» дан- ные по месяцам. Разумеется, это можно было сделать и при помощи уже известных вам функций:
> normality3(ozone.month[ozone.month$Month==5,][1])
$Ozone
[1] "NOT NORMAL"
> normality3(ozone.month[ozone.month$Month==6,][1])
134
Анализ связей: двумерные данные
$Ozone
[1] "NORMAL"
и т. д.
Мы использовали здесь две пары квадратных скобок, для того что- бы данные не превратились в вектор, а остались таблицей данных, ведь наша функция normality3() «заточена» именно под таблицы данных.
Если же эти две пары квадратных скобок вам не нравятся, можно по- ступить еще проще (хотя и муторнее) и вернуться к оригинальному shapiro.test()
:
> shapiro.test(ozone.month[ozone.month$Month==5,"Ozone"])
Shapiro-Wilk normality test data:
ozone.month[ozone.month$Month == 5, "Ozone"]
W = 0.714, p-value = 8.294e-06
(Загляните в предыдущую главу, чтобы посмотреть, какая в этом тесте альтернативная гипотеза.)
Как видите, в R одну и ту же задачу всегда можно решить множест- вом способов. И это прекрасно!
Ответ к задаче про оценки двух классов
. Попробуем проана- лизировать данные при помощи статистических тестов:
> otsenki <- read.table("data/otsenki.txt")
> klassy <- split(otsenki$V1, otsenki$V2)
Можно было сделать и по-другому, но функция split() — самое быстрое решение для разбивки данных на три группы, и к
тому же она выдает список, элементы которого можно будет проверить на нор- мальность при помощи созданной нами в предыдущей главе функции normality3()
:
> normality3(klassy)
$A1
[1] "NOT NORMAL"
$A2
[1] "NOT NORMAL"
$B1
[1] "NOT NORMAL"
Увы, параметрические методы неприменимы. Надо работать непа- раметрическими методами:
Если выборок больше двух
135
> lapply(klassy, function(.x) median(.x, na.rm=TRUE))
$A1
[1] 4
$A2
[1] 4
$B1
[1] 5
Мы применили анонимную функцию, для того чтобы избавиться от пропущенных значений (болевшие ученики?). Похоже, что в классе A
результаты первой и второй четвертей похожи, а вот класс B учился в первой четверти лучше A. Проверим это:
> wilcox.test(klassy$A1, klassy$A2, paired=TRUE)
Wilcoxon signed rank test with continuity correction data:
klassy$A1 and klassy$A2
V = 15.5, p-value = 0.8605
alternative hypothesis: true location shift is not equal to 0
Warning messages:
> wilcox.test(klassy$B1, klassy$A1, alt="greater")
Wilcoxon rank sum test with continuity correction data:
klassy$B1 and klassy$A1
W = 306, p-value = 0.02382
alternative hypothesis: true location shift is greater than 0
Warning message:
In wilcox.test.default(klassy$B1, klassy$A1, alt = "greater") :
cannot compute exact p-value with ties
Для четвертей мы применили парный тест: ведь оценки получали одни и те же ученики. А для сравнения разных классов использован од- носторонний тест, поскольку такие тесты обычно чувствительнее дву- сторонних и поскольку здесь как раз и можем проверить, учится ли класс B лучше A, а не просто отличаются ли их результаты.
Итак, результаты класса A за первую и вторую четверти достоверно не отличаются, при этом класс B учился в первую четверть статисти- чески значимо лучше класса A.
136
Анализ связей: двумерные данные
Ответ к задаче про кассиров
. Попробуем сначала выяснить,
можно ли использовать параметрические методы:
> kass <- read.table("data/kass.txt", h=TRUE)
> head(kass)
KASS.1 KASS.2 1
3 12 2
12 12 3
13 9
4 5
6 5
4 2
6 11 9
> normality3(kass)
$KASS.1
[1] "NORMAL"
$KASS.2
[1] "NORMAL"
Отлично! Как насчет средних?
> (kass.m <- sapply(kass, mean))
KASS.1
KASS.2 8.380952 7.809524
Похоже, что у первого кассира очереди побольше. Проверим это:
> with(kass, t.test(KASS.1, KASS.2))
Welch Two Sample t-test data:
KASS.1 and KASS.2
t = 0.4358, df = 39.923, p-value = 0.6653
alternative hypothesis:
true difference in means is not equal to 0 95 percent confidence interval:
-2.078962 3.221819
sample estimates:
mean of x mean of y
8.380952 7.809524
Увы, достоверной разницы нет. Попробуем отобразить это графиче- ски (рис. 27) и проверить, выполняется ли «правило трех сигм»:
Если выборок больше двух
137
> (kass.sd <- sapply(kass, sd))
KASS.1
KASS.2 4.341384 4.154745
> library(gplots)
> barplot2(kass.m, plot.ci=TRUE,
+ ci.l=kass.m, ci.u=(kass.m + (3 * kass.sd)),
+ names.arg=c("Первый кассир","Второй кассир"),
+ ylab="Очередь")
Первый кассир
Второй кассир
О
че ре дь
0 5
10 15 20
Рис. 27. Столбчатые диаграммы с разбросами в три стандартных от- клонения отображают изменчивость очередей к двум кассирам в су- пермаркете
Столбчатые диаграммы с разбросами — это распространенный, хотя и не самый эффективный, способ изображения данных, ящики-с-усами гораздо лучше. Для того чтобы нарисовать разброс на столбике, мы загрузили пакет gplots и вызвали оттуда функцию barplot2(). Кроме того, мы использовали прием с внешними скобками, чтобы результат функции sapply() шел и на экран, и в объект kass.sd.
138
Анализ связей: двумерные данные
«Правило трех сигм» нас не подвело — разбросы шириной в три стандартных отклонения сильно перекрываются. Значит, выборки до- стоверно не отличаются.
Ответ к задаче про
проращивание зараженных грибом се- мян. Загрузим данные, посмотрим на их структуру и применим тест хи-квадрат:
> pr <- read.table("data/prorostki.txt", h=TRUE)
> head(pr)
CID GERM.14 1
63 1
2 63 1
3 63 1
4 63 1
5 63 1
6 63 1
> chisq.test(table(pr[pr$CID %in% c(0,105),]))
Pearson’s Chi-squared test with Yates’ continuity correction data:
table(pr[pr$CID %in% c(0, 105), ])
X-squared = 8.0251, df = 1, p-value = 0.004613
> chisq.test(table(pr[pr$CID %in% c(0,80),]))
Pearson’s Chi-squared test with Yates’ continuity correction data:
table(pr[pr$CID %in% c(0, 80), ])
X-squared = 22.7273, df = 1, p-value = 1.867e-06
> chisq.test(table(pr[pr$CID %in% c(0,63),]))
Pearson’s Chi-squared test with Yates’ continuity correction data:
table(pr[pr$CID %in% c(0, 63), ])
X-squared = 0.2778, df = 1, p-value = 0.5982
Warning message:
In chisq.test(table(pr[pr$CID %in% c(0, 63), ])) :
Chi-squared approximation may be incorrect
Как видим, эффекты двух грибов (CID105 и CID80) отличаются от контроля, а третьего (CID63) — нет.
Если выборок больше двух
139
В ответе использован оператор %in%, который позволяет выбрать из того, что слева, то, что есть справа. Это сэкономило нам немного вре- мени. Можно было написать и table(pr[pr$CID==0 & pr$CID==63, ]),
где & означает логическое «и». Есть, разумеется, и другие способы сде- лать то же самое.
Нужно заметить, что поскольку мы три раза повторили сравнение с одной и той же группой (контролем), нужно пороговое p-значение уменьшить в три раза (это называется «поправка Бонферрони для мно- жественных сравнений»). Кстати, для подобных поправок в R есть спе- циальная функция, p.adjust(). В нашем случае первые два значения все равно останутся значимыми.
Ответ к задаче про вес
. Да и нет:
> anova(lm(data$VES data$CVET))
Analysis of Variance Table
Response: data$VES
Df Sum Sq Mean Sq F value
Pr(>F)
data$CVET
2 1043.5 521.73 31.598 4.837e-11 ***
Residuals 87 1436.5 16.51
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> pairwise.t.test(data$VES, data$CVET)
Pairwise comparisons using t tests with pooled SD
data:
data$VES and data$CVET
bl br br 1.5e-10 - sh 0.15 7.5e-08
P value adjustment method: holm
Отличия между группами есть, но отличаются по весу как раз брю- неты, а не блондины.
Ответ к задаче о показе предметов
. Поступим так же, как и в примере про программистов. Сначала загрузим данные и присоединим их, чтобы было легче работать с переменными:
> pokaz <- read.table("data/pokaz.txt")
> attach(pokaz)
Затем проверим модель:
140
Анализ связей: двумерные данные
> pokaz.logit <- glm(formula=V3 V2, family=binomial)
> summary(pokaz.logit)
Call:
glm(formula = V3 V2, family = binomial)
Deviance Residuals:
Min
1Q
Median
3Q
Max
-2.4029
-0.8701 0.4299 0.7825 1.5197
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)
-1.6776 0.7923
-2.117 0.03423 *
V2 0.9015 0.2922 3.085 0.00203 **
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 62.687
on 49
degrees of freedom
Residual deviance: 49.738
on 48
degrees of freedom
AIC: 53.738
Number of Fisher Scoring iterations: 4
(Вызывая переменные, мы учли тот факт, что R присваивает колон- кам данных без заголовков имена V1, V2, V3 и т. д.)
Модель значима! Так что наша идея, высказанная в первой главе (о том, что человек обучается в ходе попыток), имеет под собой опреде- ленные основания. А вот так можно построить график для этого экс- перимента (рис. 28):
> popytki <- seq(1, 5, 0.081) # ровно 50 значений от 1 до 5
> pokaz.p <- predict(pokaz.logit, list(V2=popytki),
+ type="response")
> plot(V3 jitter(V2, amount=.1))
> lines(popytki, pokaz.p)
> detach(pokaz)
Мы использовали predict(), чтобы рассчитать вероятности успеха в по-разному «расположенных» попытках, а на графике, чтобы точки не накладывались друг на друга, добавили немного случайных откло- нений при помощи функции jitter(). Не забудьте отсоединить данные командой detach().
Если выборок больше двух
141 1
2 3
4 5
0.
0 0.
2 0.
4 0.
6 0.
8 1.
0
Рис. 28. Графическое изображение логистической модели (данные о по- казе предметов)
Ответ к задаче про зависимость веса от роста
. Есть значимая корреляция, но зависимость слабая (рис. 29):
> cor.test(data$VES, data$ROST)
Pearson’s product-moment correlation data:
data$VES and data$ROST
t = 4.8868, df = 88, p-value = 4.565e-06
alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval:
0.2818862 0.6106708
sample estimates:
cor
0.4620071
> ves.rost <- lm(data$VES data$ROST)
142
Анализ связей: двумерные данные
> summary(ves.rost)
Call:
lm(formula = data$VES data$ROST)
Residuals:
Min
1Q
Median
3Q
Max
-10.2435
-3.4745
-0.5642 2.9358 12.3203
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.69580 13.37252 0.800 0.426
data$ROST
0.39742 0.08132 4.887 4.57e-06 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.708 on 88 degrees of freedom
Multiple R-squared: 0.2135,Adjusted R-squared: 0.2045
F-statistic: 23.88 on 1 and 88 DF,
p-value: 4.565e-06
> plot(data$VES data$ROST)
> abline(ves.rost)
Вывод о том, что зависимость слабая, сделан на основании низко- го R
2
и незначимого первого коэффициента. Кроме того, очень много остатков (больше 70% длины интервала между минимумом и максиму- мом) находится снаружи от первого и третьего квартилей. Да и график показывает серьезный разброс.
Если выборок больше двух
143 155 160 165 170 175 180 70 75 80 85 90
Рис. 29. Слабая зависимость
Глава 6
Анализ структуры: data mining
Фразу «data mining» можно все чаще увидеть в Интернете и на обложках книг по анализу данных. Говорят даже, что эпоха статистики одной-двух переменных закончилась, и наступило новое время — время анализа больших и сверхбольших массивов данных.
На самом деле под методами «data mining»
подразумеваются лю- бые методы, как визуальные, так и аналитические, позволяющие «на- щупать» структуру в данных, особенно в данных большого размера.
Данные для такого анализа используются, как правило, многомерные,
то есть такие, которые можно представить в виде таблицы из несколь- ких колонок-переменных. Поэтому более традиционное название для этих методов — «многомерный анализ», или «многомерная статистика».
Но «data mining» звучит, конечно, серьезнее. Кроме многомерности и большого размера (сотни, а то и тысячи строк и столбцов), используе- мые данные отличаются еще и тем, что переменные в них могут быть совершенно разных типов (интервальные, шкальные, номинальные).
Грубо говоря, многомерные методы делятся на методы визуализа- ции и методы классификации с обучением. В первом случае результат можно анализировать в основном зрительно, а во втором — возможна статистическая проверка результатов. Разумеется, граница между эти- ми группами нерезкая, но для удобства мы станем рассматривать их именно в этом порядке.
6.1. Рисуем многомерные данные
Самое простое, что можно сделать с многомерными данными,— это построить график. Разумеется, для того чтобы построить график несколь- ких переменных, надо свести все разнообразие к двум или, в крайнем случае, к трем измерениям. Это называют «сокращением размерности».
Но иногда можно обойтись и без сокращения размерности. Например,
если переменных три.
Рисуем многомерные данные
145 6.1.1. Диаграммы рассеяния
Если все три переменные — непрерывные, то поможет пакет RGL,
который позволяет создавать настоящие трехмерные графики.
Для примера всюду, где возможно, мы будем использовать встро- енные в R данные iris. Эти данные, заимствованные из работы зна- менитого математика (и биолога) Р. Фишера, описывают разнообразие нескольких признаков трех видов ирисов. Соответственно, в них 5 пе- ременных (колонок), причем последняя — это название вида.
Вот как можно изобразить 4 из 5 колонок при помощи RGL (рис. 30):
> library(rgl)
> plot3d(iris$Sepal.Length, iris$Sepal.Width, iris$Petal.Length,
+ col=as.numeric(iris$Species), size=3)
1 2
iris$Sepal.Width
3 4
5 5
6 7
iris$Petal.Length iris$Sepal.Length
6 4.0 3.5 3.0 2.5 2.0
Рис. 30. Пример трехмерного графика RGL
Размер появившегося окна и проекцию можно (и нужно) менять при помощи мышки.
Сразу видно, что один из видов (Iris setosa) хорошо отличается от двух других по признаку длины лепестков (Petal.Length). Кста-
146
Анализ структуры: data mining ти, в том, что мы изобразили именно 4 признака, не было оговорки,
ведь вид ириса — тоже признак, закодированный в данном случае цве- том. Можно обойтись и без RGL, тогда вам может потребоваться пакет scatterplot3d
, содержащий одноименную функцию. Но, вообще гово- ря, трехмерными графиками лучше не злоупотреблять — очень они не раскрывают, а затемняют суть явления. Правда, в случае RGL это ком- пенсируется возможностью свободно менять «точку обзора».
Другой способ визуализации многомерных данных — это построение составных графиков. Здесь у R колоссальные возможности, предостав- ляемые пакетом lattice, который предназначен для так называемой панельной (Trellis) графики. Вот как можно изобразить четыре при- знака ирисов (рис. 31):
> library(lattice)
> xyplot(Sepal.Length Petal.Length + Petal.Width | Species,
+ data=iris, auto.key=TRUE)
Petal.Length + Petal.Width
Se pa l.L
en gt h
5 6
7 8
0 2
4 6
setosa versicolor
5 6
7 8
virginica
Petal.Length
Petal.Width
Рис. 31. Пример Trellis-графика
Рисуем многомерные данные
147
Получилась картинка зависимости длины чашелистиков как от дли- ны, так и от ширины лепестков для каждого из трех видов.
Можно обойтись и без lattice. Несколько подобных графиков до- ступны и в базовой поставке R. Например, coplot() очень похож на xyplot()
(рис. 32):
> coplot(dolja javka | kand, data=vybory2)
0.
0 0.
2 0.
4 0.
6 0.
8 1.
0 0.5 0.6 0.7 0.8 0.9 1.0 0.5 0.6 0.7 0.8 0.9 1.0 0.
0 0.
2 0.
4 0.
6 0.
8 1.
0
javka do lja
KAND.1
KAND.2
KAND.3
Given : kand
Рис. 32. Данные о выборах из раздела о регрессионном анализе, пред- ставленные с помощью функции coplot()
Мы визуализировали данные о выборах из раздела про регресси- онный анализ. Эти данные тоже можно назвать многомерными, ведь каждый кандидат, по сути,— это отдельная переменная.
Матричный график рисуется при помощи функции pairs() (рис. 33):
> pairs(iris[1:4], pch=21, bg=c("red", "green3", "blue")
+ [unclass(iris$Species)])
148
Анализ структуры: data mining
Sepal.Length
2.
0 3.
0 4.
0 4.5 5.5 6.5 7.5 0.
5 1.
5 2.
5 2.0 3.0 4.0
Sepal.Width
Petal.Length
1 2 3 4 5 6 7 0.5 1.5 2.5 4.
5 5.
5 6.
5 7.
5 1
2 3
4 5
6 7
Petal.Width
Рис. 33. Матричный график
Мы получили зависимость
значения каждого признака от каждого признака, причем заодно и «покрасили» точки в цвета видов (для это- го пришлось «деклассировать» фактор iris$Species, превратив его в текстовый вектор). Таким образом, нам удалось отобразить сразу пять переменных.
6.1.2. Пиктограммы
Для визуализации многомерных данных можно использовать раз- нообразные графики-пиктограммы, где каждый элемент представля- ет один объект наблюдений, а его параметры характеризуют значения признаков объекта.
Один из классических графиков-пиктограмм — звезды (рис. 34):
> stars(mtcars[1:9,1:7], cex=1.2)
Здесь каждая пиктограмма — это один тип автомобиля, а длины лу- чей соответствуют значениям разных характеристик этих машин. Легко
Рисуем многомерные данные
149
Mazda RX4
Mazda RX4 Wag
Datsun 710
Hornet 4 Drive
Hornet Sportabout
Valiant
Duster 360
Merc 240D
Merc 230
Рис. 34. Звезды изображают разные марки автомобилей видеть, например, обе «Мазды» схожи с друг с другом больше, чем с остальными марками.
А вот более экзотический график-пиктограмма, так называемые «ли- ца Чернова». Идея графика состоит в том, что лица люди различают очень хорошо (рис. 35):
> library(TeachingDemos)
> faces(mtcars[1:9,1:7])
Здесь каждое лицо соответствует одному исследуемому объекту (ти- пу автомобиля), а черты лица характеризуют значения признаков объ- екта. Сходство «Мазд» хорошо заметно и на этом графике.
«Идеологически» к пиктограммам близок график особого типа —
график параллельных координат
. Его можно построить при помощи функции parcoord() из пакета MASS. Чтобы понять, что он делает, луч- ше всего сразу посмотреть на результат на рис. 36 (для графика ис- пользованы данные об измерениях растений):
150
Анализ структуры: data mining
Mazda RX4
Mazda RX4 Wag
Datsun 710
Hornet 4 Drive
Hornet Sportabout
Valiant
Duster 360
Merc 240D
Merc 230
Рис. 35. Лица Чернова изображают разные марки автомобилей
> measurements <- read.table("data/eq-s.txt", h=TRUE, sep=";")
> library(MASS)
> parcoord(measurements[,-1],
+ col=rep(rainbow(54), table(measurements[,1])))
С одной стороны, это очень простой график, нечто вроде многомер- ного боксплота (или гистограммы). Каждый столбец данных (у нас их всего девять, но первый — это номер местообитания) представлен в виде одной оси, на которой равномерно отмечаются все возможные значения данного признака для каждой строчки данных (в нашем случае — для каждого отдельного растения). Потом для каждого растения эти точ- ки соединяются линиями. В дополнение мы раскрасили эти линии в цвета местообитаний (всего у нас было 54 местообитания). Однако на такой график можно смотреть очень долго и каждый раз находить там интересные детали, рассказывающие о структуре данных. Например,
видно, что некоторые признаки измерены дискретно, даже невзирая на то, что они непрерывные (например, DL.TR.V); видно, что некоторые популяции попадают примерно в одно и то же место на каждом при-
Тени многомерных облаков: анализ главных компонент
151
DL.R
DIA.ST
N.REB
N.ZUB
DL.OSN.Z DL.TR.V
DL.BAZ
DL.PER
Рис. 36.
График параллельных координат знаке; видно, что у некоторых признаков (типа DL.BAZ) изменчивость смещена в сторону меньших значений, и многое другое. Именно поэто- му график параллельных координат (как и «heatmap», описанный в разделе о корреляциях) почти обязательно находится в арсенале любой программы, занимающейся «data mining».
6.2. Тени многомерных облаков: анализ главных компонент
Перейдем теперь к методам сокращения размерности. Самый рас- пространенный из них — это анализ главных компонент. Суть его в том,
что наши объекты можно представить как точки в n-мерном простран- стве, где n — это число анализируемых признаков. Через полученное облако точек проводится прямая, так, чтобы учесть наибольшую до- лю изменчивости признаков, то есть пронизывая это облако вдоль в наиболее вытянутой его части (это можно представить себе как грушу,
152
Анализ структуры: data mining надетую на стальной стержень), так получается первая главная ком- понента. Затем через это облако проводится вторая, перпендикулярная первой прямая, так, чтобы учесть наибольшую оставшуюся долю измен- чивости признаков,— как вы уже догадались, это будет вторая главная компонента. Эти две компоненты образуют плоскость, на которую и проецируются все точки.
В результате все признаки-колонки преобразуются в компоненты,
причем наибольшую информацию о разнообразии объектов несет пер- вая компонента, вторая несет меньше информации, третья — еще мень- ше и т. д. Таким образом, хотя компонент получается столько же, сколь- ко
изначальных признаков, в первых двух-трех из них сосредоточена почти вся нужная нам информация. Поэтому их можно использовать для визуализации данных на плоскости, обычно первой и второй (реже первой и третьей) компоненты. Компоненты часто называют «фактора- ми», и это порождает некоторую путаницу с похожим на анализ глав- ных компонент факторным анализом, преследующим, однако, совсем другие цели (мы его рассматривать не будем).
Вот как делается анализ главных компонент на наших данных про ирисы:
> iris.pca <- princomp(scale(iris[,1:4]))
(Мы употребили функцию scale() для того, чтобы привести все четыре переменные к одному масштабу.)
Теперь посмотрим на сами компоненты (рис. 37):
> plot(iris.pca, main="")
Это служебный график, так называемый «screeplot» (в буквальном переводе «график осыпи»), показывающий относительные вклады каж- дой компоненты в общий разброс данных. Хорошо видно, что компо- нент четыре, как и признаков, но, в отличие от первоначальных призна- ков, наибольший вклад вносят первые две компоненты. Вместо графика можно получить то же самое в текстовом виде, написав:
> 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
Отметьте, что первые две компоненты вместе объясняют почти 96%
разброса (как говорят, их вклады, «proportions of variance» = 73% и
23% соответственно).
Тени многомерных облаков: анализ главных компонент
153
Comp.1
Comp.2
Comp.3
Comp.4 0.
0 0.
5 1.
0 1.
5 2.
0 2.
5
Рис. 37. График, отражающий вклады каждой компоненты в разброс данных
Перейдем к собственно визуализации (рис. 38):
> 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() поз- воляет расположить исходные случаи (строки) в пространстве вновь найденных компонент.
Иногда полезной бывает и функция biplot() (рис. 39):
> biplot(iris.pca)
154
Анализ структуры: data mining
−3
−2
−1 0
1 2
3
−
2
−
1 0
1 2
PC1
P
C
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 s
s 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
a a a
a a
a a
a a
a a
Рис. 38. Разнообразие ирисов на графике первых двух главных компо- нент
Парный график дает возможность понять, насколько силен вклад каждого из четырех исходных признаков в первые две компоненты. В
данном случае хорошо видно, что признаки длины и ширины лепест- ков (в отличие от признаков длины и ширины чашелистиков) вносят гораздо больший вклад в первую компоненту, которая, собственно, и
«различает» виды. К сожалению, графиком, выдаваемым при помощи biplot()
, довольно трудно «управлять». Поэтому часто предпочтитель- нее функция loadings():
> loadings(iris.pca)
Loadings:
Comp.1 Comp.2 Comp.3 Comp.4
Sepal.Length
0.521 -0.377 0.720 0.261
Sepal.Width
-0.269 -0.923 -0.244 -0.124
Petal.Length
0.580
-0.142 -0.801
Petal.Width
0.565
-0.634 0.524
Тени многомерных облаков: анализ главных компонент
155
−0.2
−0.1 0.0 0.1 0.2
−
0.
2
−
0.
1 0.
0 0.
1 0.
2 1
2 3
4 5
6 7
8 9
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
−10
−5 0
5 10
−
10
−
5 0
5 10
Sepal.Length
Sepal.Width
Petal.Length
Petal.Width
Рис. 39. Парный график показывает вклад каждого признака в первые две компоненты
Comp.1 Comp.2 Comp.3 Comp.4
SS loadings
1.00 1.00 1.00 1.00
Proportion Var
0.25 0.25 0.25 0.25
Cumulative Var
0.25 0.50 0.75 1.00
Собственно, показывает она то же самое, что и biplot(), но по- дробнее: в первой части вывода каждая ячейка таблицы соответствует вкладу признака в определенный компонент. Чем
ближе это значение по модулю к единице, тем больше вклад.
Пакеты ade4 и vegan реализуют множество вариаций анализа глав- ных компонент, но самое главное — содержат гораздо больше возможно- стей для визуализации. Вот как можно проанализировать те же данные в пакете ade4 (рис. 40):
> library(ade4)
> iris.dudi <- dudi.pca(iris[,1:4], scannf=FALSE)
156
Анализ структуры: data mining
> s.class(iris.dudi$li, iris[,5])
d = 1 setosa versicolor virginica
Рис. 40. Разнообразие ирисов на графике первых двух главных компо- нент (пакет ade4)
Скачать 3.04 Mb.