А. Б. Шипунов, Е. М. Балдин, П. А. Волкова, А. И. Коробейников, С. А. Назарова
 Скачать 3.04 Mb.
|
96 Анализ связей: двумерные данные Вполне возможно, что 1001-ый (или 1 000 001-ый) факт опровергнет эту теорию. Поэтому в любом статистическом тесте выдвигаются две про- тивоположные гипотезы. Одна — это то, что мы хотим доказать (но не можем!),— называется альтернативной гипотезой (ее обозначают H 1 ). Другая — противоречащая альтернативной — называется нулевой гипотезой (обозначается H 0 ). Нулевая гипотеза всегда является пред- положением об отсутствии чего-либо (например, зависимости одной переменной от другой или различия между двумя выборками). Мы не можем доказать альтернативную гипотезу, а можем лишь опроверг- нуть нулевую гипотезу и принять альтернативную (чувствуете раз- ницу?). Если же мы не можем опровергнуть нулевую гипотезу, то мы вынуждены принять ее. 5.1.2. Статистические ошибки Естественно, что когда мы делаем любые предположения (в нашем случае — выдвигаем статистические гипотезы), мы можем ошибаться (в нашем случае — делать статистические ошибки). Всего возможно четыре исхода (табл. 5.1). Выборка Популяция Верна H 0 Верна H 1 Принимаем H 0 Правильно! Статистическая ошибка второго рода Принимаем H 1 Статистическая ошибка первого рода Правильно! Таблица 5.1. Статистические ошибки Если мы приняли для выборки H 0 (нулевую гипотезу) и она вер- на для генеральной совокупности (популяции), то мы правы, и все в порядке. Аналогично и для H 1 (альтернативной гипотезы). Ясно, что мы не можем знать, что в действительности верно для генеральной со- вокупности, и сейчас просто рассматриваем все логически возможные варианты. А вот если мы приняли для выборки альтернативную гипотезу, а она оказалась неверна для генеральной совокупности, то мы соверши- ли так называемую статистическую ошибку первого рода (нашли не- существующую закономерность). Вероятность того, что мы соверши- ли эту ошибку (это и есть p-значение, «p-value»), всегда отображается Есть ли различие, или Тестирование двух выборок 97 при проведении любых статистических тестов. Очевидно, что если ве- роятность этой ошибки достаточно высока, то мы должны отвергнуть альтернативную гипотезу. Возникает естественный вопрос: какую веро- ятность считать достаточно высокой? Так же как и с размером выбор- ки (см. первую главу), однозначного ответа на этот вопрос нет. Более или менее общепринято, что пороговым значением надо считать 0.05 (то есть альтернативная гипотеза отвергается, если вероятность ошиб- ки при ее принятии больше или равна 5%). В медицине ценой ошибки нередко являются человеческие жизни, поэтому там пороговое значе- ние часто принимают равным 0.01 или даже 0.001 (то есть решение о существовании закономерности принимается, если вероятность ошибки ничтожна). Таким образом, решение о результатах статистических тестов при- нимается главным образом на основании вероятности статистической ошибки первого рода. Степень уверенности исследователя в том, что за- ключение, сделанное на основании статистической выборки, будет спра- ведливо и для генеральной совокупности, отражает статистическая до- стоверность. Допустим, если вероятность статистической ошибки пер- вого рода равна 3%, то говорят, что найденная закономерность досто- верна с вероятностью 97%. А если вероятность статистической ошибки первого рода равна, например, 23%, то говорят, что достоверной зако- номерности не найдено. В случае, если мы принимаем нулевую гипотезу для выборки, в то время как для генеральной совокупности справедлива альтернативная гипотеза, то мы совершаем статистическую ошибку второго рода (не за- мечаем существующей закономерности). Этот параметр характеризует так называемую мощность (power) статистического теста. Чем меньше вероятность статистической ошибки второго рода (то есть чем меньше вероятность не заметить несуществующую закономерность), тем более мощным является тест. Для полной ясности перерисуем нашу таблицу для конкретного слу- чая — исследования зависимости между двумя признаками (табл. 5.2). 5.2. Есть ли различие, или Тестирование двух выборок При ответе на этот вопрос нужно всегда помнить, что упомянутые ниже тесты проверяют различия только по центральным значениям (например, средним) и подразумевают, что разброс данных в выборках примерно одинаков. Например, выборки с одинаковыми параметрами средней тенденции и разными показателями разброса данных, c(1, 2, 98 Анализ связей: двумерные данные Выборка Популяция Верна H 0 Верна H 1 Принимаем H 0 Принимаем H 1 Таблица 5.2. Статистические ошибки при исследовании зависимости между двумя признаками. Мелкие точки отражают объекты генераль- ной совокупности (популяции), крупные точки — объекты из нашей выборки 3, 4, 5, 6, 7, 8, 9) и c(5, 5, 5, 5, 5, 5, 5, 5, 5) различаться не будут. Как вы помните, для проведения статистического теста нужно вы- двинуть две статистические гипотезы. Нулевая гипотеза здесь: «разли- чий между (двумя) выборками нет» (то есть обе выборки взяты из од- ной генеральной совокупности). Альтернативная гипотеза: «различия между (двумя) выборками есть». Напоминаем, что ваши данные должны быть организованы в виде двух векторов — отдельных или объединенных в лист или таблицу дан- ных. Например, если нужно узнать, различается ли достоверно рост мужчин и женщин, то в одном векторе должен быть указан рост муж- чин, а в другом — рост женщин (каждая строчка — это один обследо- ванный человек). Если данные параметрические, то нужно провести параметрический «t-тест» (или «тест Стьюдента»). Если при этом переменные, которые мы хотим сравнить, были получены на разных объектах, мы будем ис- пользовать двухвыборочный t-тест для независимых переменных (two sample t-test), который запускается при помощи команды t.test(). На- пример, если две сравниваемые выборки записаны в первом и втором столбцах таблицы данных data, то команда t.test(data[,1], data[,2]) Есть ли различие, или Тестирование двух выборок 99 по умолчанию выдаст нам результаты двухвыборочного t-теста для независимых переменных. Если же пары сравниваемых характеристик были получены на од- ном объекте , то есть переменные зависимы (например, частота пульса до и после физической нагрузки измерялась у одного и того же челове- ка), надо использовать парный t-тест (paired t-test). Для этого в коман- де t.test() надо указать параметр paired=TRUE. В нашем примере, ес- ли данные зависимы, надо использовать команду вида t.test(data[,1], data[,2], paired=TRUE) Второй тест, тест для зависимых переменных, более мощный. Пред- ставьте себе, что мы измеряли пульс до нагрузки у одного человека, а после нагрузки — у другого. Тогда было бы не ясно, как объяснить полученную разницу: может быть, частота пульса увеличилась после нагрузки, а может быть, этим двум людям вообще свойственна разная частота пульса. В случае же «двойного» измерения пульса каждый че- ловек как бы является своим собственным контролем, и разница между сравниваемыми переменными (до и после нагрузки) обусловливается только тем фактором, на основе которого они выделены (наличием на- грузки). Если мы имеем дело с непараметрическими данными, то нужно провести непараметрический двухвыборочный тест Вилкоксона, «Wil- coxon test» (он известен еще и как как тест Манна-Уитни, «Mann- Whitney test»). Для этого надо использовать команду wilcox.test(). В случае с зависимыми выборками, аналогично t-тесту, надо использо- вать параметр paired=TRUE. Приведем несколько примеров. Для начала воспользуемся классиче- ским набором данных, который использовался в оригинальной работе Стьюдента (псевдоним ирландского математика Уильяма Сили Госсе- та). В упомянутой работе производилось сравнение влияния двух раз- личных снотворных на увеличение продолжительности сна (рис. 15). В R эти данные доступны под названием sleep. В столбце extra со- держится среднее приращение продолжительности сна после начала приема лекарства (по отношению к контрольной группе), а в столбце group — код лекарства (первое или второе). > plot(extra group, data = sleep) Здесь использована так называемая «формула модели»: в рассмат- риваемом случае extra group означает, что group используется для разбивки extra. Влияние лекарства на каждого человека индивидуально, но сред- нее увеличение продолжительности сна можно считать вполне логич- ным показателем «силы» лекарства. Основываясь на этом предположе- 100 Анализ связей: двумерные данные 1 2 − 1 0 1 2 3 4 5 Рис. 15. Среднее приращение продолжительности сна после начала при- ема разных лекарств в двух группах по отношению к контрольной нии, попробуем проверить при помощи t-теста, значимо ли различие в средних для этих двух выборок (соответствующих двум разным лекар- ствам): > with(sleep, t.test(extra[group == 1], extra[group == 2], + var.equal = FALSE)) Welch Two Sample t-test data: extra[group == 1] and extra[group == 2] t = -1.8608, df = 17.776, p-value = 0.0794 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -3.3654832 0.2054832 sample estimates: mean of x mean of y Есть ли различие, или Тестирование двух выборок 101 0.75 2.33 Параметр var.equal позволяет выбрать желаемый вариант крите- рия: оригинальный t-критерий Стьюдента в предположении, что раз- бросы данных одинаковы(var.equal = TRUE), или же t-тест в модифи- кации Уэлча (Welch), свободный от этого предположения (var.equal = FALSE ). Хотя формально мы не можем отвергнуть нулевую гипотезу (о ра- венстве средних), p-значение (0.0794) все же достаточно маленькое, что- бы попробовать другие методы для проверки гипотезы, увеличить ко- личество наблюдений, еще раз убедиться в нормальности распределе- ний и т. д. Может прийти в голову провести односторонний тест — ведь он обычно чувствительнее. Так вот, этого делать нельзя! Нельзя пото- му, что большинство статистических тестов рассчитаны на то, что они будут проводиться ad hoc, то есть без знания какой-либо дополнитель- ной информации. Существуют и post hoc тесты (скажем, обсуждаемый ниже «Tukey Honest Significant Differences test»), но их немного. Для сравнения двух выборок существуют, разумеется, и непарамет- рические тесты. Один из них, тест знаков, настолько прост, что его нет в R. Можно, однако, легко сделать его самому. Тест знаков вычисляет разницу между всеми парами элементов двух одинаковых по размеру выборок (то есть это парный тест). Затем можно отбросить все отри- цательные значения и оставить только положительные. Если выборки взяты из одной генеральной совокупности, то положительных разниц будет примерно половина, и тогда уже знакомый нам биномиальный тест не найдет отличий между 50% и пропорцией положительных раз- ниц. Если выборки разные, то пропорция положительных разниц будет значимо больше или меньше половины. Задача . Придумайте, при помощи какого R-выражения сделать та- кой тест, и протестируйте две выборки, которые мы привели в начале раздела. Перейдем тем временем к более сложным непараметрическим тес- там. Рассмотрим пример. Стандартный набор данных airquality со- держит информацию о величине озона в воздухе города Нью-Йорка с мая по сентябрь 1973 года. Концентрация озона представлена округ- ленными средними значениями за день, поэтому анализировать ее «от греха подальше» лучше непараметрическими методами. Впрочем, по- пробуйте выяснить самостоятельно (это задача!), насколько близки к нормальным помесячные распределения значений концентраций озона (см. ответ в конце главы). Проверим, например, гипотезу о том, что распределение уровня озо- на в мае и в августе было одинаковым: > wilcox.test(Ozone Month, data = airquality, 102 Анализ связей: двумерные данные + subset = Month %in% c(5, 8)) Wilcoxon rank sum test with continuity correction data: Ozone by Month W = 127.5, p-value = 0.0001208 alternative hypothesis: true location shift is not equal to 0 Поскольку Month дискретен (это просто номер месяца), то значения Ozone будут сгруппированы помесячно. Кроме того, мы использовали параметр subset вместе с оператором %in%, который выбирает изо всех месяцев май и август (пятый и восьмой месяцы). Критерий отвергает гипотезу о согласии распределений уровня озо- на в мае и в августе с достаточно большой надежностью. В это до- статочно легко поверить, потому что уровень озона в воздухе сильно зависит от солнечной активности, температуры и ветра. Различия между выборками хорошо иллюстрировать при помощи графика «ящик-с-усами», например (рис. 16): > boxplot(Ozone Month, data = airquality, + subset = Month %in% c(5, 8)) (Обратите внимание на то, что здесь для boxplot() мы использова- ли ту же самую формулу модели.) Часто считают, что если «ящики» перекрываются более чем на треть своей высоты, то выборки достоверно не различаются. Можно использовать t-тест и тест Вилкоксона и для одной выбор- ки, если стоит задача сравнить ее с неким «эталоном». Поскольку вы- борка одна, то и соответствующие тесты называют одновыборочными. Нулевую гипотезу в этом случае надо формулировать как равенство выборочной средней или медианы (для t-теста и теста Вилкоксона, со- ответственно) заданному числу µ (см. предыдущую главу). Задача . В файле данных otsenki.txt записаны оценки одних и тех же учеников некоторого школьного класса за первую четверть (значе- ние A1 во второй колонке) и за вторую четверть (A2), а также оценки учеников другого класса за первую четверть (B1). Отличаются ли ре- зультаты класса A за первую и вторую четверть? Какой класс учился в первой четверти лучше — A или B? Ответ см. в конце главы. Задача . В супермаркете, среди прочих, работают два кассира. Для того чтобы проанализировать эффективность их работы, несколько раз в день в середине недели подсчитывали очереди к каждой из них. Дан- ные записаны в файле kass.txt. Какой кассир работает лучше? Есть ли соответствие, или Анализ таблиц 103 5 8 0 50 10 0 15 0 Рис. 16. Распределение озона в воздухе в мае и июне 5.3. Есть ли соответствие, или Анализ таблиц А как сравнить между собой две выборки из номинальных (кате- гориальных) данных? Для этого часто используют таблицы сопряжен- ности (contigency tables). Построить таблицу сопряженности можно с помощью функции table(): > with(airquality, table(cut(Temp, quantile(Temp)), Month)) Month 5 6 7 8 9 (56,72] 24 3 0 1 10 (72,79] 5 15 2 9 10 (79,85] 1 7 19 7 5 (85,97] 0 5 10 14 5 Строки этой таблицы — четыре интервала температур (по Фаренгей- ту), а столбцы — месяц года. На пересечении каждой строки и столбца 104 Анализ связей: двумерные данные находится число, которое показывает, сколько значений в данном меся- це попадает в данный интервал температур. Если факторов больше двух, то R будет строить многомерную табли- цу и выводить ее в виде серии двумерных таблиц, что не всегда удобно. Можно, однако, построить «плоскую» таблицу сопряженности, когда все факторы, кроме одного, объединяются в один «многомерный» фак- тор, чьи градации используются при построении таблицы. Построить такую таблицу можно с помощью функции ftable(): > ftable(Titanic, row.vars = 1:3) Survived No Yes Class Sex Age 1st Male Child 0 5 Adult 118 57 Female Child 0 1 Adult 4 140 2nd Male Child 0 11 Adult 154 14 Female Child 0 13 Adult 13 80 3rd Male Child 35 13 Adult 387 75 Female Child 17 14 Adult 89 76 Crew Male Child 0 0 Adult 670 192 Female Child 0 0 Adult 3 20 Параметр row.vars позволяет указать номера переменных в набо- ре данных, которые следует объединить в один-единственный фактор, градации которого и будут индексировать строки таблицы сопряжен- ности. Параметр col.vars проделывает то же самое, но для столбцов таблицы. Функцию table можно использовать и для других целей. Самое простое — это подсчет частот. Например, можно считать пропуски: > d <- factor(rep(c("A","B","C"), 10), levels=c("A","B","C","D", + "E")) > is.na(d) <- 3:4 > table(factor(d, exclude = NULL)) Есть ли соответствие, или Анализ таблиц 105 A B C 9 10 9 2 Функция mosaicplot позволяет получить графическое изображение таблицы сопряженности (рис. 17): > titanic <- apply(Titanic, c(1, 4), sum) > titanic > titanic Survived Class No Yes 1st 122 203 2nd 167 118 3rd 528 178 Crew 673 212 > mosaicplot(titanic, col = c("red", "green"), main = "", + cex.axis=1) При помощи функции chisq.test() можно статистически прове- рить гипотезу о независимости двух факторов. К данным будет приме- нен тест хи-квадрат (Chi-squared test). Например, проверим гипотезу о независимости цвета глаз и волос (встроенные данные HairEyeColor): > x <- margin.table(HairEyeColor, c(1, 2)) > chisq.test(x) Pearson’s Chi-squared test data: x X-squared = 138.2898, df = 9, p-value < 2.2e-16 (Того же эффекта можно добиться, если функции summary() пере- дать таблицу сопряженности в качестве аргумента.) Набор данных HairEyeColor — это многомерная таблица сопряжен- ности. Для суммирования частот по всем «измерениям», кроме двух, использовалась функция margin.table. Таким образом, в результате была получена двумерная таблица сопряженности. Тест хи-квадрат в качестве нулевой гипотезы принимает независимость факторов, так что в нашем примере (поскольку мы отвергаем нулевую гипотезу) следует принять, что факторы обнаруживают соответствие. Чтобы графически изобразить эти соответствия, можно воспользо- ваться функцией assocplot() (рис. 18): > x <- margin.table(HairEyeColor, c(1,2)) > assocplot(x) 106 Анализ связей: двумерные данные Class Su rv iv ed 1st 2nd 3rd Crew N o Y es Рис. 17. Доля выживших на «Титанике» пассажиров и членов команды в зависимости от класса их билета: мозаичный график На графике видны отклонения ожидаемых частот от наблюдаемых величин. Высота прямоугольника показывает абсолютную величину это- го отклонения, а положение — знак отклонения. Отчетливо видно, что для людей со светлыми волосами характерен голубой цвет глаз и совсем не характерен карий цвет, а для обладателей черных волос ситуация об- ратная. Чтобы закрепить изученное о таблицах сопряженности, рассмотрим еще один интересный пример. Однажды большая компания статистиков- эпидемиологов собралась на банкет. На следующее утро после банкета, в воскресенье, многие встали с симптомами пищевого отравления, а в понедельник в институте много сотрудников отсутствовало по болезни. Поскольку это были статистики, и не просто статистики, а эпидемио- логи, то они решили как следует вспомнить, кто что ел на банкете, и тем самым выяснить причину отравления. Получившиеся данные име- ют следующий формат: > tox <- read.table("data/otravlenie.txt", h=TRUE) Есть ли соответствие, или Анализ таблиц 107 Black Brown Red Blond G re en H az el B lu e B ro w n Рис. 18. Связь между цветом волос и цветом глаз: график сопряженно- сти > head(tox) ILL CHEESE CRABDIP CRISPS BREAD CHICKEN RICE CAESAR TOMATO 1 1 1 1 1 2 1 1 1 1 2 2 1 1 1 2 1 2 2 2 3 1 2 2 1 2 1 2 1 2 4 1 1 2 1 1 1 2 1 2 5 1 1 1 1 2 1 1 1 1 6 1 1 1 1 1 1 2 1 1 ICECREAM CAKE JUICE WINE COFFEE 1 1 1 1 1 1 2 1 1 1 1 2 3 1 1 2 1 2 4 1 1 2 1 2 5 2 1 1 1 1 6 2 1 1 2 2 108 Анализ связей: двумерные данные Первая переменная (ILL) показывает, отравился участник банкета (1) или нет (2), остальные переменные соответствуют разным блюдам. Простой взгляд на данные ничего не даст, ведь на банкете было 45 человек и 13 блюд! Значит, надо попробовать статистические методы. Так как данные номинальные, то можно попробовать проанализировать таблицы сопряженности 1 : > for (m in 2:ncol(tox)) + { + tmp <- chisq.test(tox$ILL, tox[,m]) + print(paste(names(tox)[m], tmp$p.value)) + } [1] "CHEESE 0.840899679390882" [1] "CRABDIP 0.94931385140737" [1] "CRISPS 0.869479670886473" [1] "BREAD 0.349817724258644" [1] "CHICKEN 0.311482217451896" [1] "RICE 0.546434435905853" [1] "CAESAR 0.000203410168460333" [1] "TOMATO 0.00591250292451728" [1] "ICECREAM 0.597712594782716" [1] "CAKE 0.869479670886473" [1] "JUICE 0.933074267280188" [1] "WINE 0.765772843686273" [1] "COFFEE 0.726555246056369" Небольшая хитрость (цикл for) позволила нам не писать тест три- надцать раз подряд (хотя copy-paste никто не отменял, и можно было бы сделать так, как сделано в разделе «Прогноз» главы про временные ряды). Без функции print() цикл бы ничего не напечатал, а paste() помогла оформить результат. Ну а результат таков, что есть два зна- чимых соответствия — с цезарским салатом (CAESAR) и с помидорами (TOMATO). Посмотрим, что нам покажет график ассоциаций (рис. 19): > assocplot(table(ILL=tox$ILL, CAESAR=tox$CAESAR)) Практически такая же картина и по тем, кто ел помидоры. Винов- ник найден! Почти. Ведь вряд ли испорченными были сразу оба блюда. Надо теперь попробовать разобраться, что же было главной причиной. Для этого вам придется заглянуть ниже, туда, где мы знакомимся с логистической регрессией. 1 На предупреждения можно не обращать внимания. Полученные значения лучше дополнительно обработать с помощью p.adjust(). Есть ли соответствие, или Анализ таблиц 109 1 2 2 1 ILL C A E SA R Рис. 19. Соответствие между отравившимися и теми, кто ели цезарский салат * * * Кроме методов работы с таблицами сопряженности (таких, напри- мер, как тест хи-квадрат), номинальные и шкальные данные можно обрабатывать различными, более специализированными методами. На- пример, для сравнения результатов экспертных оценок популярны так называемые тесты согласия (concordance). Среди этих тестов широко распространен тест Коэна (Cohen), который вычисляет так называе- мую каппу Коэна (Cohen’s kappa) — меру согласия, изменяющуюся от 0 до 1, и вдобавок вычисляет p-значение для нулевой гипотезы (о том, что каппа равна 0). Приведем такой пример: в 2003 году две группы независимо (и примерно в одно и то же время года) обследовали один и тот же остров Белого моря. Целью было составить список всех видов растений, встреченных на острове. Таким образом, данные были бинар- ными (0 — вида на острове нет, 1 — вид на острове есть). Результаты обследований записаны в файл pokorm03.dat: 110 Анализ связей: двумерные данные > pok <- read.table("data/pokorm_03.dat", h=TRUE, sep=";") > library(psych) > cohen.kappa(pok) Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries lower estimate upper unweighted kappa 0.55 0.67 0.8 weighted kappa 0.55 0.67 0.8 Number of subjects = 193 (Мы загрузили пакет psych, потому что именно в нем находится нужная нам функция cohen.kappa().) Каппа довольно высока, нижний предел 95% доверительного интер- вала выше 0.5. Это значит, что результаты исследования можно считать согласными друг с другом. Задача . В файле данных prorostki.txt находятся результаты экс- перимента по проращиванию семян васильков, зараженных различны- ми грибами (колонка CID, CID=0 — это контроль, то есть незараженные семена). Всего исследовали по 20 семян на каждый гриб, тестировали три гриба, итого с контролем — 80 семян. Отличается ли прорастание семян, зараженных грибами, от контроля? Ответ см. в конце главы. 5.4. Есть ли взаимосвязь, или Анализ корреляций Мерой линейной взаимосвязи между переменными является коэф- фициент корреляции Пирсона (обозначается латинской буквой r). Зна- чения коэффициента корреляции могут изменяться по модулю от ну- ля до единицы. Нулевой коэффициент корреляции говорит о том, что значения одной переменной не связаны со значениями другой перемен- ной, а коэффициент корреляции, равный единице (или минус единице), свидетельствует о четкой линейной связи между переменными. Поло- жительный коэффициент корреляции говорит о положительной взаи- мосвязи (чем больше, тем больше), отрицательный — об отрицательной (чем больше, тем меньше). Степень взаимосвязи между переменными отражает коэффициент детерминации : это коэффициент корреляции, возведенный в квадрат. Эта величина показывает, какая доля изменений значений одной пере- менной сопряжена с изменением значений другой переменной. Напри- мер, если коэффициент корреляции увеличится вдвое (0.8), то степень взаимосвязи между переменными возрастет в четыре раза (0.8 2 = 0.64 ). Есть ли взаимосвязь, или Анализ корреляций 111 Повторим еще раз, что коэффициент корреляции характеризует ме- ру линейной связи между переменными. Две переменные могут быть очень четко взаимосвязаны, но если эта связь не линейная, а, допу- стим, параболическая, то коэффициент корреляции будет близок к ну- лю. Примером такой связи может служить связь между степенью воз- бужденности человека и качеством решения им математических задач. Очень слабо возбужденный человек (скажем, засыпающий) и очень сильно возбужденный (во время футбольного матча) будет решать за- дачи гораздо хуже, чем умеренно возбужденный человек. Поэтому, пе- ред тем как оценить взаимосвязь численно (вычислить коэффициент корреляции), нужно посмотреть на ее графическое выражение. Луч- ше всего здесь использовать диаграмму рассеяния, или коррелограмму (scatterplot) — именно она вызывается в R командой plot() с двумя аргументами-векторами. Очень важно также, что всюду в этом разделе речь идет о наличии и силе взаимосвязи между переменными, а не о характере этой взаимо- связи. Если мы нашли достоверную корреляцию между переменными А и Б, то это может значить, что: • А зависит от Б; • Б зависит от А; • А и Б зависят друг от друга; • А и Б зависят от какой-то третьей переменной В, а между собой не имеют ничего общего. Например, хорошо известно, что объем продаж мороженого и чис- ло пожаров четко коррелируют. Странно было бы предположить, что поедание мороженого располагает людей к небрежному обращению с огнем или что созерцание пожаров возбуждает тягу к мороженому. Все гораздо проще — оба этих параметра зависят от температуры воздуха! Для вычисления коэффициента корреляции в R используется функ- ция cor(): > cor(5:15, 7:17) [1] 1 > cor(5:15, c(7:16, 23)) [1] 0.9375093 В простейшем случае ей передаются два аргумента (векторы одина- ковой длины). Кроме того, можно вызвать ее с одним аргументом, если это — матрица или таблица данных. В этом последнем случае функ- 112 Анализ связей: двумерные данные ция cor() вычисляет так называемую корреляционную матрицу, со- ставленную из коэффициентов корреляций между столбцами матрицы или набора данных, взятых попарно, например: > cor(trees) Girth Height Volume Girth 1.0000000 0.5192801 0.9671194 Height 0.5192801 1.0000000 0.5982497 Volume 0.9671194 0.5982497 1.0000000 Если все данные присутствуют, то все просто, но что делать, когда есть пропущенные наблюдения? Для этого в команде cor есть пара- метр use. По умолчанию он равен all.obs, что при наличии хотя бы одного пропущенного наблюдения приводит к сообщению об ошибке. Если use приравнять к значению complete.obs, то из данных авто- матически удаляются все наблюдения с пропусками. Может оказаться так, что пропуски раскиданы по исходному набору данных достаточно хаотично и их много, так что после построчного удаления от матрицы фактически ничего не остается. В таком случае поможет попарное уда- ление пропусков, то есть удаляются строчки с пропусками не из всей матрицы сразу, а только лишь из двух столбцов непосредственно перед вычислением коэффициента корреляции. Для этого опцию use следует приравнять к значению pairwise.complete.obs (надо иметь в виду, что в в этом коэффициенты корреляции вычисляются по разному количе- ству наблюдений и сравнивать их друг с другом может быть опасно). Если данные непараметрические, нужно использовать ранговый ко- эффициент корреляции Спирмена (Spearman) ρ. Он менее подвержен влиянию случайных «выбросов» в данных, чем коэффициент Пирсона. Для подсчета ρ достаточно приравнять параметр method к значению spearman : > x <- rexp(50) > cor(x, log(x), method="spearman") [1] 1 Если корреляционная матрица большая, то читать ее довольно труд- но. Поэтому существует несколько способов визуального представления таких матриц. Можно, например, использовать функцию symnum(), ко- торая выведет матрицу в текстовом виде, но с заменой чисел на буквы в зависимости от того, какому диапазону принадлежало значение: > symnum(cor(longley)) GNP. GNP U A P Y E GNP.deflator 1 Есть ли взаимосвязь, или Анализ корреляций 113 GNP B 1 Unemployed , , 1 Armed.Forces . 1 Population B B , . 1 Year B B , . B 1 Employed B B . . B B 1 attr(,"legend") [1] 0 ‘ ’ 0.3 ‘.’ 0.6 ‘,’ 0.8 ‘+’ 0.9 ‘*’ 0.95 ‘B’ 1 Эта функция имеет большое количество разнообразных настроек, но по умолчанию они все выставлены в значения, оптимальные для отображения корреляционных матриц. Второй способ — это графическое представление корреляционных коэффициентов. Идея проста: нужно разбить область от −1 до +1 на отдельные диапазоны, назначить каждому свой цвет, а затем все это отобразить. Для этого можно воспользоваться функциями image() и axis() (рис. 20): > cor.l <- cor(longley) > image(1:ncol(cor.l), 1:nrow(cor.l), cor.l, + col=heat.colors(22), axes=FALSE, xlab="", ylab="") # Подписи к осям: > axis(1, at=1:ncol(cor.l), labels=abbreviate(colnames(cor.l))) > axis(2, at=1:nrow(cor.l), labels=abbreviate(rownames(cor.l)), + las = 2) (Мы сократили здесь длинные названия строк и столбцов при по- мощи команды abbreviate().) Полученный график часто называют «heatmap» («карта темпера- туры»). Еще один интересный способ представления корреляционной мат- рицы предоставляется пакетом ellipse. В этом случае значения коэф- фициентов корреляции рисуются в виде эллипсов. Чем ближе значение коэффициента корреляции к +1 или −1 — тем более вытянутым стано- вится эллипс. Наклон эллипса отражает знак. Для получения изобра- жения необходимо вызвать функцию plotcorr (рис. 21): > library(ellipse) > cor.l <- cor(longley) > colnames(cor.l) <- abbreviate(colnames(cor.l)) > rownames(cor.l) <- abbreviate(rownames(cor.l)) > plotcorr(cor.l, type="lower", mar=c(0,0,0,0)) 114 Анализ связей: двумерные данные GNP. GNP Unmp Ar.F Pplt Year Empl GNP. GNP Unmp Ar.F Pplt Year Empl Рис. 20. Графическое представление корреляционной матрицы А как проверить статистическую значимость коэффициента корре- ляции? Это равносильно проверке статистической гипотезы о равенстве нулю коэффициента корреляции. Если гипотеза отвергается, то связь одного признака с другим считается значимой. Для проверки такой ги- потезы используется функция cor.test(): > with(trees, cor.test(Girth, Height)) Pearson’s product-moment correlation data: Girth and Height t = 3.2722, df = 29, p-value = 0.002758 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.2021327 0.7378538 sample estimates: cor 0.5192801 Какая связь, или Регрессионный анализ 115 GNP Unmp Ar.F Pplt Year Empl G N P G N P U nm p A r. F P pl t Y ea r Рис. 21. Коэффициенты корреляции в виде эллипсов Логика рассуждений здесь абсолютно такая же, как и в рассмотрен- ных выше тестах. В данном случае нам нужно принять альтернативную гипотезу о том, что корреляция действительно существует. Обратите внимание на доверительный интервал — тест показывает, что реальное значение корреляции может лежать в интервале 0.2–0.7. 5.5. Какая связь, или Регрессионный анализ Корреляционный анализ позволяет определить, зависимы ли пере- менные, и вычислить силу этой зависимости. Чтобы определить тип зависимости и вычислить ее параметры, используется регрессионный анализ. Наиболее известна двумерная линейная регрессионная модель: m = b 0 + b 1 × x, 116 Анализ связей: двумерные данные где m — модельное (predicted) значение зависимой переменной y, x — независимая переменная, а b 0 и b 1 — параметры модели (коэффи- циенты). Такая модель позволяет оценить среднее значение переменной y при известном значении x. Разность между истинным значением и модель- ным называется ошибкой («error») или остатком («residual»): E = y − m В идеальном варианте остатки имеют нормальное распределение с нулевым средним и неизвестным, но постоянным разбросом σ 2 , который не зависит от значений x и y. В этом случае говорят о гомогенности остатков. Если же разброс зависит еще от каких-либо параметров, то остатки считаются гетерогенными. Кроме того, если обнаружилось, что средние значения остатков зависят от x, то это значит, что между y и x имеется нелинейная зависимость. Для того чтобы узнать значения параметров b 0 и b 1 , применяют метод наименьших квадратов. Затем вычисляют среднеквадратичное отклонение R 2 : R 2 = 1 − σ 2 m /σ 2 y , где σ 2 y — разброс переменной y. Для проверки гипотезы о том, что модель значимо отличается от нуля, используется так называемая F-статистика (статистика Фишера). Как обычно, если p-значение меньше уровня значимости (обычно 0.05), то модель считается значимой. Вот пример. Во встроенной таблице данных women содержатся 15 наблюдений о росте (дюймы) и весе (фунты) женщин в возрасте от 30 до 39 лет. Попробуем узнать, как вес зависит от роста (рис. 22): # Преобразование в метрическую систему: > women.metr <- women > women.metr$height <- 0.0254*women.metr$height > women.metr$weight <- 0.45359237*women.metr$weight # Вычисление параметров уравнения регрессии: > lm.women<-lm(formula = weight height, data = women.metr) > lm.women$coefficients (Intercept) height -39.69689 61.60999 # Вывод модельных значений: > b0 <- lm.women$coefficient[1] > b1 <- lm.women$coefficient[2] > x1 <- min(women.metr$height) Какая связь, или Регрессионный анализ 117 > x2 <- max(women.metr$height) > x <- seq(from = x1, to = x2, length.out =100) > y <- b0 + b1*x # Вывод графика зависимости: > plot(women.metr$height, women.metr$weight, main="", + xlab="Рост (м)", ylab="Вес (кг)") > grid() > lines(x, y, col="red") 1.50 1.55 1.60 1.65 1.70 1.75 1.80 55 60 65 70 75 Рост (м) В ес ( кг ) Рис. 22. Зависимость между ростом и весом Для просмотра результатов линейной аппроксимации можно исполь- зовать функцию > summary(lm.women) Call: lm(formula = weight height, data = women.metr) 118 Анализ связей: двумерные данные Residuals: Min 1Q Median 3Q Max -0.7862 -0.5141 -0.1739 0.3364 1.4137 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -39.697 2.693 -14.74 1.71e-09 *** height 61.610 1.628 37.85 1.09e-14 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.6917 on 13 degrees of freedom Multiple R-squared: 0.991,Adjusted R-squared: 0.9903 F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14 on 1 and 13 DF, p-value: 1.091e-14 Отсюда следует, что: 1. Получена линейная регрессионная модель вида Веc(мод) = -39.697 + 61.61 * Рост , то есть увеличение роста на 10 см соответствует увеличению веса примерно на 6 кг. 2. Наибольшее положительное отклонение истинного значения от- клика от модельного составляет 1.4137 кг, наибольшее отрица- тельное −0.7862 кг. 3. Почти половина остатков находится в пределах от первой кварти- ли (1Q = −0.5141 кг) до третьей (3Q = 0.3364 кг). 4. Все коэффициенты значимы на уровне p-value < 0.001; это по- казывают *** в строке Coefficients) и сами значения p-value Pr(>|t|) : 1.71e-09 для b 0 («(Intercept)») и 1.09e-14 («height») для b 1 5. Среднеквадратичное отклонение (Adjusted R-squared) для дан- ной модели составляет R 2 = 0.9903 . Его значение близко к 1, что говорит о высокой значимости. 6. Значимость среднеквадратичного отклонения подтверждает и вы- сокое значение F-статистики, равное 1433, и общий уровень зна- чимости (определяемый по этой статистике): p-value: 1.091e-14, что много меньше 0.001. Какая связь, или Регрессионный анализ 119 7. Если на основе этого анализа будет составлен отчет, то надо будет указать еще и значения степеней свободы df: 1 и 13 (по колонкам и по строкам данных соответственно). На самом деле это еще далеко не конец, а только начало долго- го пути анализа и подгонки модели. В этом случае, например, пара- болическая модель (где используется не рост, а квадрат роста), будет еще лучше вписываться в данные (как выражаются, будет «лучше объ- яснять» вес). Можно попробовать и другие модели, ведь о линейных, обобщенных линейных и нелинейных моделях написаны горы книг! Но нам сейчас было важно сделать первый шаг. Время от времени требуется не просто получить регрессию, а срав- нить несколько регрессий, основанных на одних и тех же данных. Вот пример. Известно, что в методике анализа выборов большую роль иг- рает соответствие между процентом избирателей, проголосовавших за данного кандидата, и процентом явки на данном избирательном участ- ке. Эти соответствия для разных кандидатов могут выглядеть по-разно- му, что говорит о различии их электоратов. Данные из файла vybory.txt содержат результаты голосования за трех кандидатов по более чем сотне избирательных участков. Посмотрим, различается ли для разных кандидатов зависимость их результатов от явки избирателей. Сначала загрузим и проверим данные: > vybory <- read.table("data/vybory.txt", h=TRUE) > str(vybory) ’data.frame’: 153 obs. of 7 variables: $ UCHASTOK: int 1 2 3 4 5 6 7 8 9 10 ... $ IZBIR : int 329786 144483 709903 696114 ... $ NEDEJSTV : int 2623 859 5656 4392 3837 4715 ... $ DEJSTV: int 198354 97863 664195 619629 ... $ KAND.1 : int 24565 7884 30491 54999 36880 ... $ KAND.2 : int 11786 6364 11152 11519 10002 ... $ KAND.3 : int 142627 68573 599105 525814 ... > head(vybory) UCHASTOK IZBIR NEDEJSTV DEJSTV KAND.1 KAND.2 KAND.3 1 1 329786 2623 198354 24565 11786 142627 2 2 144483 859 97863 7884 6364 68573 3 3 709903 5656 664195 30491 11152 599105 4 4 696114 4392 619629 54999 11519 525814 5 5 717095 3837 653133 36880 10002 559653 6 6 787593 4715 655486 72166 25204 485669 Теперь присоединим таблицу данных, для того чтобы с ее колонками можно было работать как с независимыми переменными: 120 Анализ связей: двумерные данные > attach(vybory) Вычислим доли проголосовавших за каждого кандидата и явку: > DOLJA <- cbind(KAND.1, KAND.2, KAND.3) / IZBIR > JAVKA <- (DEJSTV + NEDEJSTV) / IZBIR Посмотрим, есть ли корреляция: > cor(JAVKA, DOLJA) KAND.1 KAND.2 KAND.3 [1,] -0.3898262 -0.41416 0.9721124 Корреляция есть, хотя и невысокая для первых двух кандидатов. Похоже, что голосование по третьему кандидату серьезно отличалось. Проверим это: > lm.1 <- lm(KAND.1/IZBIR JAVKA) > lm.2 <- lm(KAND.2/IZBIR JAVKA) > lm.3 <- lm(KAND.3/IZBIR JAVKA) > lapply(list(lm.1, lm.2, lm.3), summary) Call: lm(formula = KAND.1/IZBIR JAVKA) Residuals: Min 1Q Median 3Q Max -0.046133 -0.013211 -0.001493 0.012742 0.057943 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.117115 0.008912 13.141 < 2e-16 *** JAVKA -0.070726 0.013597 -5.202 6.33e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.01929 on 151 degrees of freedom Multiple R-squared: 0.152,Adjusted R-squared: 0.1463 F-statistic: 27.06 on 1 and 151 DF, p-value: 6.331e-07 Call: lm(formula = KAND.2/IZBIR JAVKA) Residuals: Какая связь, или Регрессионный анализ 121 Min 1Q Median 3Q Max -0.036542 -0.014191 -0.000162 0.011782 0.050139 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.093487 0.007807 11.974 < 2e-16 *** JAVKA -0.066597 0.011911 -5.591 1.03e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.0169 on 151 degrees of freedom Multiple R-squared: 0.1715,Adjusted R-squared: 0.166 F-statistic: 31.26 on 1 and 151 DF, p-value: 1.027e-07 Call: lm(formula = KAND.3/IZBIR JAVKA) Residuals: Min 1Q Median 3Q Max -0.116483 -0.023570 0.000518 0.025714 0.102810 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.43006 0.01702 -25.27 <2e-16 *** JAVKA 1.32261 0.02597 50.94 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.03684 on 151 degrees of freedom Multiple R-squared: 0.945,Adjusted R-squared: 0.9446 F-statistic: 2595 on 1 and 151 DF, p-value: < 2.2e-16 Коэффициенты регрессии у третьего кандидата сильно отличают- ся. Да и R 2 гораздо выше. Посмотрим, как это выглядит графически (рис. 23): > plot(KAND.3/IZBIR JAVKA, xlim=c(0,1), ylim=c(0,1), + xlab="Явка", ylab="Доля проголосовавших за кандидата") > points(KAND.1/IZBIR JAVKA, pch=2) > points(KAND.2/IZBIR JAVKA, pch=3) > abline(lm.3) > abline(lm.2, lty=2) > abline(lm.1, lty=3) 122 Анализ связей: двумерные данные > legend("topleft", lty=c(3,2,1), + legend=c("Кандидат 1","Кандидат 2","Кандидат 3")) > detach(vybory) 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Явка Д ол я пр ог ол ос ов ав ш их з а ка нд ид ат а Кандидат 1 Кандидат 2 Кандидат 3 Рис. 23. Зависимость результатов трех кандидатов от явки на выборы Итак, третий кандидат имеет электорат, который значительно отли- чается по своему поведению от электоратов двух первых кандидатов. Разумеется, хочется проверить это статистически, а не только визу- ально. Для того чтобы сравнивать регрессии, существует специальный метод, анализ ковариаций (ANCOVA). Подробное рассмотрение этого метода выходит за рамки настоящей книги, но мы сейчас покажем, как можно применить анализ ковариаций для данных о выборах: > vybory2 <- cbind(JAVKA, stack(data.frame(DOLJA))) > names(vybory2) <- c("javka","dolja","kand") > str(vybory2) ’data.frame’: 459 obs. of 3 variables: $ javka: num 0.609 0.683 0.944 0.896 0.916 ... Какая связь, или Регрессионный анализ 123 $ dolja: num 0.0745 0.0546 0.043 0.079 0.0514 ... $ kand : Factor w/ 3 levels "KAND.1","KAND.2",..: 1 1 1 ... > head(vybory2, 3) javka dolja kand 1 0.6094164 0.07448770 KAND.1 2 0.6832776 0.05456697 KAND.1 3 0.9435810 0.04295094 KAND.1 Мы создали и проверили новую таблицу данных. В ней две интере- сующие нас переменные (доля проголосовавших за каждого кандидата и явка) расположены теперь в один столбец, а третий столбец — это фактор с именами кандидатов. Чтобы получить такую таблицу, мы ис- пользовали функцию stack(). > ancova.v <- lm(dolja javka * kand, data=vybory2) > summary(ancova.v) Call: lm(formula = dolja javka * kand, data = vybory2) Residuals: Min 1Q Median 3Q Max -0.116483 -0.015470 -0.000699 0.014825 0.102810 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.117115 0.011973 9.781 < 2e-16 *** javka -0.070726 0.018266 -3.872 0.000124 *** kandKAND.2 -0.023627 0.016933 -1.395 0.163591 kandKAND.3 -0.547179 0.016933 -32.315 < 2e-16 *** javka:kandKAND.2 0.004129 0.025832 0.160 0.873074 javka:kandKAND.3 1.393336 0.025832 53.938 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.02591 on 453 degrees of freedom Multiple R-squared: 0.9824,Adjusted R-squared: 0.9822 F-statistic: 5057 on 5 and 453 DF, p-value: < 2.2e-16 Анализ ковариаций использует формулу модели «отклик воздей- ствие * фактор », где звездочка (*) обозначает, что надо проверить од- новременно отклик на воздействие, отклик на фактор и взаимодействие между откликом и фактором. (Напомним, что линейная регрессия ис- 124 Анализ связей: двумерные данные пользует более простую формулу, «отклик воздействие»). В нашем случае откликом была доля проголосовавших, воздействием — явка, а фактором — имя кандидата. Больше всего нас интересовало, правда ли, что имеется сильное взаимодействие между явкой и третьим кандида- том. Полученные результаты (см. строку javka:kandKAND.3) не остав- ляют в этом сомнений, взаимодействие значимо. * * * Довольно сложно анализировать данные, если зависимость нели- нейная. Здесь на помощь могут прийти методы визуализации. Более того, часто их бывает достаточно, чтобы сделать выводы, пусть и пред- варительные. Рассмотрим такой пример. По берегу Черного моря, от Новороссийска до Сочи, растут первоцветы, или примулы. Самая заме- чательная их черта — это то, что окраска цветков постепенно меняется при продвижении от Новороссийска на юго-запад, вдоль берега: сначала большинство цветков белые или желтые, а ближе к Дагомысу появля- ется все больше и больше красных, малиновых, фиолетовых и розовых тонов. В файле данных primula.txt записаны результаты десятилет- него исследования береговых популяций первоцветов. Попробуем выяс- нить, как именно меняется окраска цветков по мере продвижения на юго-запад (рис. 24): > prp.coast <- read.table("data/primula.txt", + as.is=TRUE, h=TRUE) > plot(yfrac nwse, data=prp.coast, type="n", + xlab="Дистанция от Новороссийска, км", + ylab="Пропорция светлых цветков, %") > rect(129, -10, 189, 110, col=gray(.8), border=NA); box() > mtext("129", at=128, side=3, line=0, cex=.8) > mtext("189", at=189, side=3, line=0, cex=.8) > points(yfrac nwse, data=prp.coast) > abline(lm(yfrac nwse, data=prp.coast), lty=2) > lines(loess.smooth(prp.coast$nwse, prp.coast$yfrac), lty=1) В этом небольшом анализе есть несколько интересных для нас мо- ментов. Во-первых, окраска цветков закодирована пропорцией, для того чтобы сделать ее интервальной. Во-вторых, линейная регрессия, оче- видно, не отражает реальной ситуации. В-третьих, для изучения нели- нейного взаимоотношения мы применили здесь так называемое сгла- живание, «LOESS» (Locally weighted Scatterplot Smoothing), которое помогло увидеть нам общий вид возможной кривой. Вероятность успеха, или Логистическая регрессия 125 50 100 150 200 250 0 20 40 60 80 10 0 Дистанция от Новороссийска, км П ро по рц ия с ве тл ы х цв ет ко в, % Рис. 24. Изменение пропорций светлых цветков в популяциях первоцве- та 5.6. Вероятность успеха, или Логистическая регрессия Номинальные данные очень трудно обрабатывать статистически. Тест пропорций да хи-квадрат — практически все, что имеется в на- шем арсенале для одно- и двумерных номинальных данных. Однако часто встречаются задачи, в которых требуется выяснить не просто достоверность соответствия, а его характер, то есть требуется приме- нить к номинальным данным нечто вроде регрессионного анализа. Вот, например, данные о тестировании программистов с разным опытом ра- боты. Их просили за небольшое время (скажем, двадцать минут) напи- сать несложную программу, не тестируя ее на компьютере. Потом напи- санные программы проверяли, и если программа работала, писали «S» («success», успех), а если не работала — писали «F» («failure», неудача): > l <- read.table("data/logit.txt", stringsAsFactors=TRUE) > head(l) 126 Анализ связей: двумерные данные V1 V2 1 14 F 2 29 F 3 6 F 4 25 S 5 18 S 6 4 F Как узнать, влияет ли стаж на успех? Таблицы сопряженности тут подходят мало, ведь стаж (V1) — интервальная переменная. Линейная регрессия не подходит, потому что для нее требуется не только интер- вальное воздействие (в нашем случае стаж), но еще и интервальный отклик. У нас же отклик — бинарный. Оказывается, существует вы- ход. Для этого надо в качестве отклика исследовать не успех/неуспех, а вероятность успеха, которая, как всякая вероятность, непрерывно изменяется в интервале от 0 до 1. Мы не будем здесь вдаваться в мате- матические подробности, а просто покажем, как это работает для дан- ных о программистах: > l.logit <- glm(formula=V2 V1, family=binomial, data=l) > summary(l.logit) Call: glm(formula = V2 V1, family = binomial, data = l) Deviance Residuals: Min 1Q Median 3Q Max -1.9987 -0.4584 -0.2245 0.4837 1.5005 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.9638 2.4597 -2.018 0.0436 * V1 0.2350 0.1163 2.021 0.0432 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 18.249 on 13 degrees of freedom Residual deviance: 10.301 on 12 degrees of freedom AIC: 14.301 Number of Fisher Scoring iterations: 5 Вероятность успеха, или Логистическая регрессия 127 Итак, оба параметра логистической регрессии значимы, и p-значения меньше пороговых. Этого достаточно, чтобы сказать: опыт програм- миста значимо влияет на успешное решение задачи. Попробуйте теперь решить такую задачу. В файле pokaz.txt нахо- дятся результаты эксперимента с показом предметов на доли секунды, описанного в первой главе. Первая колонка содержит номер испытуе- мого (их было десять), вторая — номер испытания (пять на человека, испытания проводились последовательно), а третья — результат в би- нарной форме: 0 — не угадал(а), 1 — угадал(а). Выясните, есть ли связь между номером испытания и результатом. В выдаче функции summary.glm() стоит обратить внимание на пред- последнюю строчку, AIC (Akaike’s Information Criterion) — информаци- онный критерий Акаике. Это критерий, который позволяет сравнивать модели между собой (несколько подробнее про AIC рассказано в главе про временные ряды). Чем меньше AIC, тем лучше модель. Покажем это на примере. Вернемся к нашим данным про «виновника» отравле- ния. Помидоры или салат? Вот так можно попробовать решить пробле- му: > tox.logit <- glm(formula=I(2-ILL) CAESAR + TOMATO, + family=binomial, data=tox) > tox.logit2 <- update(tox.logit, . . - TOMATO) > tox.logit3 <- update(tox.logit, . . - CAESAR) Сначала мы построили логистическую модель. Поскольку логисти- ческая модель «хочет» в качестве отклика (слева от тильды) получить 0 или 1, то мы просто вычли значения переменной ILL из 2, при этом заболевшие стали закодированы как 0, а здоровые — как 1. Функция I() позволяет нам избежать интерпретации минуса как части формулы модели, то есть возвращает минусу обычные арифметические свойства. Затем при помощи update() мы видоизменили исходную модель, убрав из нее помидоры, а потом убрали из исходной модели салат (точки в формуле модели означают, что используются все воздействия и все от- клики). Теперь посмотрим на AIC: > tox.logit$aic [1] 47.40782 > tox.logit2$aic [1] 45.94004 > tox.logit3$aic [1] 53.11957 Что же, похоже, что настоящим виновником отравления является салат, ведь модель с салатом обладает наименьшим AIC. Почему же 128 Анализ связей: двумерные данные тест хи-квадрат указал нам еще и на помидоры? Дело, скорее всего, в том, что между салатом и помидорами существовало взаимодействие (interaction). Говоря человеческим языком, те, кто брали салат, обычно брали к нему и помидоры. 5.7. Если выборок больше двух А что если теперь мы захотим узнать, есть ли различия между тремя выборками? Первое, что приходит в голову (предположим, что это параметрические данные),— это провести серию тестов Стьюдента: между первой и второй выборками, между первой и третьей и, нако- нец, между второй и третьей — всего три теста. К сожалению, чис- ло необходимых тестов Стьюдента будет расти чрезвычайно быстро с увеличением числа интересующих нас выборок. Например, для сравне- ния попарно шести выборок нам понадобится провести уже 15 тестов! А представляете, как обидно будет провести все эти 15 тестов толь- ко для того, чтобы узнать, что все выборки не различаются между собой! Но главная проблема заключена не в сбережении труда исследо- вателя (все-таки обычно нам нужно сравнить не больше 3–4 выборок). Дело в том, что при повторном проведении статистических тестов, ос- нованных на вероятностных понятиях, на одной и той же выборке ве- роятность обнаружить достоверную закономерность по ошибке воз- растает . Допустим, мы считаем различия достоверными при p-value < 0.05 , при этом мы будем ошибаться (находить различия там, где их нет) в 4 случаях из 100 (в 1 случае из 25). Понятно, что если мы про- ведем 25 статистических тестов на одной и той же выборке, то, скорее всего, однажды мы найдем различия просто по ошибке. Аналогичные рассуждения могут быть применены и к экстремальным видам спорта. Например, вероятность того, что парашют не раскроется при прыжке, довольно мала (допустим, 1/100), и странно бы было ожидать, что па- рашют не раскроется как раз тогда, когда человек прыгает впервые. При этом любой десантник, имеющий опыт нескольких сотен прыж- ков, может рассказать несколько захватывающих историй о том, как ему пришлось использовать запасной парашют. Поэтому для сравнения трех и более выборок используется специ- альный метод — однофакторный дисперсионный анализ (ANOVA, от английского «ANalysis Of VAriance»). Нулевая гипотеза здесь будет утверждать, что выборки не различаются между собой, а альтерна- тивная гипотеза — что хотя бы одна пара выборок различается между собой. Обратите внимание на формулировку альтернативной гипотезы! Результаты этого теста будут одинаковыми и в случае, если различается только одна пара выборок, и в случае, если различаются все выборки. Если выборок больше двух 129 Чтобы узнать, какие именно выборки отличаются, надо будет прово- дить специальные тесты (о них рассказано ниже). Данные для ANOVA лучше организовать следующим образом: за- дать две переменные, и в одной из них указать все значения всех срав- ниваемых выборок (например, рост брюнетов, блондинов и шатенов), а во второй — коды выборок, к которым принадлежат значения первой переменной (например, будем ставить напротив значения роста брюне- та «br», напротив роста блондина — «bl» и напротив роста шатена — «sh»). Тогда мы сможем задействовать уже знакомую нам логическую регрессию (при этом отклик будет количественный, а вот воздействие — качественное). Сам анализ проводится при помощи двойной функции anova(lm()) (рис. 25): > set.seed(1683) > VES.BR <- sample(70:90, 30, replace=TRUE) > VES.BL <- sample(69:79, 30, replace=TRUE) > VES.SH <- sample(70:80, 30, replace=TRUE) > ROST.BR <- sample(160:180, 30, replace=TRUE) > ROST.BL <- sample(155:160, 30, replace=TRUE) > ROST.SH <- sample(160:170, 30, replace=TRUE) > data <- data.frame(CVET=rep(c("br", "bl", "sh"), each=30), + VES=c(VES.BR, VES.BL, VES.SH), + ROST=c(ROST.BR, ROST.BL, ROST.SH)) > boxplot(data$ROST data$CVET) > anova(lm(data$ROST data$CVET)) Analysis of Variance Table Response: data$ROST Df Sum Sq Mean Sq F value Pr(>F) data$CVET 2 2185.2 1092.58 81.487 < 2.2e-16 *** Residuals 87 1166.5 13.41 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Надо теперь объяснить, что именно мы сделали. В R дисперсион- ный анализ проводится при помощи той же самой функции линейной модели lm(), что и регрессия, но формула модели здесь другая: отклик фактор , где фактором является индикатор группы (в нашем случае c("br", "bl", "sh") ). Такой анализ сначала строит три модели. Каж- дую модель можно условно представить в виде вертикальной линии, по которой распределены точки; вертикальной потому, что все точки со- ответствуют одному значению индикатора группы (скажем, bl). Потом функция anova() сравнивает эти модели и определяет, есть ли между 130 Анализ связей: двумерные данные bl br sh 15 5 16 0 16 5 17 0 17 5 18 0 Рис. 25. Различается ли вес между тремя выборками людей с разным цветом волос? (Данные получены искусственно) ними какие-нибудь значимые отличия. Мы создали наши данные ис- кусственно, как и данные о зарплате в главе про одномерные данные. Ну и, естественно, получили ровно то, что заложили,— надо принять альтернативную гипотезу о том, что хотя бы одна выборка отличается от прочих. Глядя на боксплот, можно предположить, что это — блон- дины. Проверим предположение при помощи множественного t-теста (pairwise t-test): > pairwise.t.test(data$ROST, data$CVET) Pairwise comparisons using t tests with pooled SD data: data$ROST and data$CVET bl br br < 2e-16 - sh 1.1e-11 1.2e-05 P value adjustment method: holm Если выборок больше двух 131 И вправду, из полученной таблички видно, что блондины значимо (оба p-значения много меньше 0.05) отличаются от остальных групп. Есть и еще один похожий тест. Это «Tukey Honest Significant Differences test» (у него есть и специальный график, см. рис. 26): > rost.cvet <- aov(lm(data$ROST data$CVET)) > (rost.cvet.hsd <- TukeyHSD(rost.cvet)) Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = lm(data$ROST data$CVET)) $‘data$CVET‘ diff lwr upr p adj br-bl 11.933333 9.678935 14.187732 0.00e+00 sh-bl 7.533333 5.278935 9.787732 0.00e+00 sh-br -4.400000 -6.654399 -2.145601 3.43e-05 > plot(rost.cvet.hsd) Результаты обоих тестов, разумеется, похожи. Важно помнить, что ANOVA — параметрический метод, то есть на- ши данные должны иметь нормальное (или близкое к нормальному) распределение. Кроме того, стандартные отклонения выборок должны быть, по крайней мере, похожи. Последнее условие можно обойти при помощи функции onеway.test(), которая не предполагает по умолча- нию равенства разбросов. Но если есть сомнения в параметричности ис- ходных данных, то лучше всего использовать непараметрический ана- лог ANOVA, тест Краскала-Уоллиса (Kruskal-Wallis test). Можно по- пробовать его на наших искусственных данных: > kruskal.test(data$ROST data$CVET) Kruskal-Wallis rank sum test data: data$ROST by data$CVET Kruskal-Wallis chi-squared = 63.9153, df = 2, p-value = 1.321e-14 Результат аналогичен «классическому» ANOVA, а p-значение чуть больше, как и положено для непараметрического теста. Задача про вес . Попробуйте самостоятельно выяснить, отлича- ются ли наши группы еще и по весу? Если да, то отличаются ли от остальных по-прежнему блондины? 132 Анализ связей: двумерные данные −5 0 5 10 sh − br sh − bl br − bl Рис. 26. Результаты post hoc теста Тьюки. Обе разницы между средними для роста блондинов и двух остальных групп лежат внутри доверитель- ного интервала (справа от пунктирной линии) Задача про зависимость роста и веса |