А. Б. Шипунов, Е. М. Балдин, П. А. Волкова, А. И. Коробейников, С. А. Назарова
 Скачать 3.04 Mb.
|
|
3.468127 61.103302 Одномерные статистические тесты 85 sample estimates: mean of x 32.28571 Это вариант теста Стьюдента для одномерных данных. Статисти- ческие тесты (в том числе и этот) пытаются высчитать так называе- мую тестовую статистику, в данном случае статистику Стьюдента (t- статистику). Затем на основании этой статистики рассчитывается «p- величина» (p-value), отражающая вероятность ошибки первого рода. А ошибкой первого рода (ее еще называют «ложной тревогой»), в свою очередь, называется ситуация, когда мы принимаем так называемую альтернативную гипотезу, в то время как на самом деле верна нулевая (гипотеза «по умолчанию»). Наконец, вычисленная p-величина исполь- зуется для сравнения с заранее заданным порогом (уровнем) значимо- сти . Если p-величина ниже порога, нулевая гипотеза отвергается, если выше — принимается. Подробнее о статистических гипотезах можно прочесть в следующей главе. В нашем случае нулевая гипотеза состоит в том, что истинное сред- нее (то есть среднее генеральной совокупности) равно вычисленному нами среднему (то есть 32.28571). Перейдем к анализу вывода функции. Статистика Стьюдента при шести степенях свободы (df=6, поскольку у нас всего 7 значений) дает единичное p-значение, то есть 100%. Какой бы распространенный по- рог мы не приняли (0.1%, 1% или 5%), это значение все равно больше. Следовательно, мы принимаем нулевую гипотезу. Поскольку альтерна- тивная гипотеза в нашем случае — это то, что «настоящее» среднее (среднее исходной выборки) не равно вычисленному среднему, то по- лучается, что «на самом деле» эти цифры статистически не отличают- ся. Кроме всего этого, функция выдает еще и доверительный интервал (confidence interval), в котором, по ее «мнению», может находиться на- стоящее среднее. Здесь он очень широк — от трех с половиной тысяч до 61 тысячи рублей. Непараметрический (то есть не связанный предположениями о рас- пределении) аналог этого теста тоже существует. Это так называемый ранговый тест Уилкоксона: > wilcox.test(salary2, mu=median(salary2), conf.int=TRUE) Wilcoxon signed rank test with continuity correction data: salary2 V = 221949, p-value = 0.8321 alternative hypothesis: true location is not equal to 21 95 percent confidence interval: 20.99999 21.00007 86 Великое в малом: одномерные данные sample estimates: (pseudo)median 21.00004 Эта функция и выводит практически то же самое, что и t.test() выше. Обратите, однако, внимание, что этот тест связан не со средним, а с медианой. Вычисляется (если задать conf.int=TRUE) и доверитель- ный интервал. Здесь он значительно 0 уже, потому что медиана намного устойчивее среднего значения. Понять, соответствует ли распределение данных нормальному (или хотя бы близко ли оно к нормальному), очень и очень важно. Напри- мер, все параметрические статистические методы основаны на предпо- ложении о том, что данные имеют нормальное распределение. Поэто- му в R реализовано несколько разных техник, отвечающих на вопрос о нормальности данных. Во-первых, это статистические тесты. Самый простой из них — тест Шапиро-Уилкса (попробуйте провести этот тест самостоятельно): > shapiro.test(salary) > shapiro.test(salary2) Но что же он показывает? Здесь функция выводит гораздо мень- ше, чем в предыдущих случаях. Более того, даже встроенная справка не содержит объяснений того, какая здесь, например, альтернативная гипотеза. Разумеется, можно обратиться к литературе, благо справка дает ссылки на публикации. А можно просто поставить эксперимент: > set.seed(1638) > shapiro.test(rnorm(100)) Shapiro-Wilk normality test data: rnorm(100) W = 0.9934, p-value = 0.9094 rnorm() генерирует столько случайных чисел, распределенных по нормальному закону, сколько указано в его аргументе. Это аналог функ- ции sample(). Раз мы получили высокое p-значение, то это свидетель- ствует о том, что альтернативная гипотеза в данном случае: «распре- деление не соответствует нормальному». Кроме того, чтобы результаты при вторичном воспроизведении были теми же, использована функция set.seed() , регулирующая встроенный в R генератор случайных чисел так, чтобы числа в следующей команде были созданы по одному и тому же «закону». Одномерные статистические тесты 87 Таким образом, распределение данных в salary и salary2 отлича- ется от нормального. Другой популярный способ проверить, насколько распределение по- хоже на нормальное,— графический. Вот как это делается (рис. 14): > qqnorm(salary2, main="") > qqline(salary2, col=2) −3 −2 −1 0 1 2 3 20 40 60 80 10 0 Рис. 14. Графическая проверка нормальности распределения Для каждого элемента вычисляется, какое место он занимает в сор- тированных данных (так называемый «выборочный квантиль») и какое место он должен был бы занять, если распределение нормальное (тео- ретический квантиль). Прямая проводится через квартили. Если точки лежат на прямой, то распределение нормальное. В нашем случае мно- гие точки лежат достаточно далеко от красной прямой, а значит, не похожи на распределенные нормально. Для проверки нормальности можно использовать и более универ- сальный тест Колмогорова—Смирнова, который сравнивает любые два 88 Великое в малом: одномерные данные распределения, поэтому для сравнения с нормальным распределением ему надо прямо указать «pnorm», то есть так называемую накопленную функцию нормального распределения (она встроена в R): > ks.test(salary2, "pnorm") One-sample Kolmogorov-Smirnov test data: salary2 D = 1, p-value < 2.2e-16 alternative hypothesis: two-sided Он выдает примерно то же, что и текст Шапиро-Уилкса. 4.4. Как создавать свои функции Тест Шапиро-Уилкса всем хорош, но не векторизован, как и многие другие тесты в R, поэтому применить его сразу к нескольким колонкам таблицы данных не получится. Сделано это нарочно, для того чтобы подчеркнуть нежелательность множественных парных сравнений (об этом подробнее написано в главе про двумерные данные). Но в нашем случае парных сравнений нет, а сэкономить время хочется. Можно, ко- нечно, нажимая «стрелку вверх», аккуратно повторить тест для каж- дой колонки, но более правильный подход — создать пользовательскую функцию. Вот пример такой функции: > normality <- function(data.f) + { + result <- data.frame(var=names(data.f), p.value=rep(0, + ncol(data.f)), normality=is.numeric(names(data.f))) + for (i in 1:ncol(data.f)) + { + data.sh <- shapiro.test(data.f[, i])$p.value + result[i, 2] <- round(data.sh, 5) + result[i, 3] <- (data.sh > .05) + } + return(result) + } Чтобы функция заработала, надо скопировать эти строчки в окно консоли или записать их в отдельный файл (желательно с расширени- ем *.r), а потом загрузить командой source(). После этого ее можно вызвать: Как создавать свои функции 89 > normality(trees) var p.value normality 1 Girth 0.08893 TRUE 2 Height 0.40342 TRUE 3 Volume 0.00358 FALSE Функция не только запускает тест Шапиро-Уилкса несколько раз, но еще и разборчиво оформляет результат выполнения. Разберем функ- цию чуть подробнее. В первой строчке указан ее аргумент — data.f. Дальше, в окружении фигурных скобок, находится тело функции. На третьей строчке формируется пустая таблица данных такой размер- ности, какая потребуется нам в конце. Дальше начинается цикл: для каждой колонки выполняется тест, а потом (это важно!) из теста из- влекается p-значение. Процедура эта основана на знании структуры вывода теста — списка, где элемент p-value содержит p-значение. Про- верить это можно, заглянув в справку, а можно и экспериментально (как? — см. ответ в конце главы). Все p-значения извлекаются, округ- ляются, сравниваются с пороговым уровнем значимости (в данном слу- чае 0.05) и записываются в таблицу. Затем таблица выдается «наружу». Предложенная функция совершенно не оптимизирована. Ее легко мож- но сделать чуть короче и к тому же несколько «смышленее», скажем так: > normality2 <- function(data.f, p=.05) + { + nn <- ncol(data.f) + result <- data.frame(var=names(data.f), p.value=numeric(nn), + normality=logical(nn)) + for (i in 1:nn) + { + data.sh <- shapiro.test(data.f[, i])$p.value + result[i, 2:3] <- list(round(data.sh, 5), data.sh > p) + } + return(result) + } > normality2(trees) Результаты, разумеется, не отличаются. Зато видно, как можно до- бавить аргумент, причем сразу со значением по умолчанию. Теперь можно писать: > normality2(trees, 0.1) 90 Великое в малом: одномерные данные var p.value normality 1 Girth 0.08893 FALSE 2 Height 0.40341 TRUE 3 Volume 0.00358 FALSE То есть если вместо 5% взять десятипроцентный порог, то уже и для первой колонки можно отвергнуть нормальное распределение. Уже не раз говорилось, что циклов в R следует избегать. Можно ли сделать это в нашем случае? Оказывается, да: > lapply(trees, shapiro.test) $Girth Shapiro-Wilk normality test data: X[[1L]] W = 0.9412, p-value = 0.08893 $Height Как видите, все еще проще! Если мы хотим улучшить зрительный эффект, можно сделать так: > lapply(trees, function(.x) ifelse(shapiro.test(.x)$p.value > + .05, "NORMAL", "NOT NORMAL")) $Girth [1] "NORMAL" $Height [1] "NORMAL" $Volume [1] "NOT NORMAL" Здесь применена так называемая анонимная функция, функция без названия, обычно употребляемая в качестве последнего аргумента ко- манд типа apply(). Кроме того, используется логическая конструкция ifelse() И наконец, на этой основе можно сделать третью пользовательскую функцию (проверьте самостоятельно, как она работает): Всегда ли точны проценты 91 > normality3 <- function(df, p=.05) + { + lapply(df, function(.x) ifelse(shapiro.test(.x)$p.value > + p, "NORMAL", "NOT NORMAL")) + } > normality3(list(salary, salary2)) > normality3(log(trees+1)) Примеры тоже интересны. Во-первых, нашу третью функцию мож- но применять не только к таблицам данных, но и к «настоящим» спис- кам, с неравной длиной элементов. Во-вторых, простейшее логарифми- ческое преобразование сразу же изменило «нормальность» колонок. 4.5. Всегда ли точны проценты Полезной характеристикой при исследовании данных является про- порция (доля). В статистике под пропорцией понимают отношение объ- ектов с исследуемой особенностью к общему числу наблюдений. По- скольку доля — это часть целого, то отношение части к целому нахо- дится в пределах от 0 до 1. Для удобства восприятия доли ее умножают на 100% и получают процент — число в пределах от 0% до 100%. Сле- дует внимательно относиться к вычислению долей и не забывать про исходные данные. В 1960-е годы в районном отделе по сельскому хозяй- ству было обнаружено, что падеж лошадей в одном сельском поселении составил 50%. Пришлось составить компетентную комиссию для про- верки причин такого ужасного явления. Прибыв на место, комиссия обнаружила, что в данном поселении было всего две лошади, в том числе и сдохшая недавно старая кляча, которая и послужила причиной формирования и выезда на место комиссии! Рассмотрим проблему, которая часто встречается при проведении статистических исследований. Как узнать, отличается ли вычисленный нами процент от «истинного» процента, то есть доли интересующих нас объектов в генеральной совокупности? Вот пример. В больнице есть группа из 476 пациентов, среди кото- рых 356 курят. Мы знаем, что в среднем по больнице доля курящих составляет 0.7 (70%). А вот в нашей группе она чуть побольше — при- мерно 75%. Для того чтобы проверить гипотезу о том, что доля куря- щих в рассматриваемой группе пациентов отличается от средней доли по больнице, мы можем задействовать так называемый биномиальный тест: > binom.test(x=356, n=476, p=0.7, alternative="two.sided") 92 Великое в малом: одномерные данные Exact binomial test data: 356 and 476 number of successes = 356, number of trials = 476, p-value = 0.02429 alternative hypothesis: true probability of success is not equal to 0.7 95 percent confidence interval: 0.7063733 0.7863138 sample estimates: probability of success 0.7478992 Поскольку p-значение значительно меньше 0.05 и альтернативная гипотеза состоит в том, что доля курильщиков (она здесь довольно издевательски называется «probability of success», «вероятность успе- ха») не равна 0.7, то мы можем отвергнуть нулевую гипотезу и при- нять альтернативную, то есть решить, что наши 74% отличаются от средних по больнице 70% не случайно. В качестве опции мы исполь- зовали alternative="two.sided", но можно было поступить иначе — проверить гипотезу о том, что доля курящих в рассматриваемой груп- пе пациентов превышает среднюю долю курящих по больнице. Тогда альтернативную гипотезу надо было бы записать как alt="greater". Кроме биномиального, мы можем применить здесь и так называе- мый тест пропорций. Надо заметить, что он применяется шире, потому что более универсален: > prop.test(x=356, n=476, p=0.7, alternative="two.sided") 1-sample proportions test with continuity correction data: 356 out of 476, null probability 0.7 X-squared = 4.9749, df = 1, p-value = 0.02572 alternative hypothesis: true p is not equal to 0.7 95 percent confidence interval: 0.7059174 0.7858054 sample estimates: p 0.7478992 Как видим, результат практически такой же. Тест пропорций можно проводить и с двумя выборками, для это- го используется все та же функция prop.test() (для двухвыборочного Всегда ли точны проценты 93 текста пропорций), а также mcnemar.test() (для теста Мак-Немара, ко- торый проводится тогда, когда выборки связаны друг с другом). Узнать, как их использовать, можно, прочитав справку (и особенно примеры!) по обеим функциям. * * * Ответ к задаче про функцию normality() . Чтобы узнать, отку- да взять значения p-value, надо сначала вспомнить, что в R почти все, что выводится на экран,— это результат «печати» различных списков при помощи невидимой команды print(). А из списка легко «достать» нужное либо по имени, либо по номеру (если, скажем, имени у эле- ментов нет — так иногда бывает). Сначала выясним, из чего состоит список-вывод функции shapiro.test(): > str(shapiro.test(rnorm(100))) List of 4 $ statistic: Named num 0.992 ..- attr(*, "names")= chr "W" $ p.value : num 0.842 $ method : chr "Shapiro-Wilk normality test" $ data.name: chr "rnorm(100)" - attr(*, "class")= chr "htest" В этом списке из 4 элементов есть элемент под названием p.value, что нам и требовалось. Проверим на всякий случай, то ли это: > set.seed(1683) > shapiro.test(rnorm(100))$p.value [1] 0.8424077 Именно то, что надо. Осталось только вставить это название в нашу функцию. Ответ к задаче про выборы из предисловия. Для того чтобы рассчитать p-value для нашей пропорции (48% про- голосовавших за кандидата B), можно воспользоваться описанным вы- ше prop.test(): > prop.test(0.48*262, 262) 1-sample proportions test with continuity correction data: 0.48 * 262 out of 262, null probability 0.5 X-squared = 0.343, df = 1, p-value = 0.5581 94 Великое в малом: одномерные данные alternative hypothesis: true p is not equal to 0.5 95 percent confidence interval: 0.4183606 0.5422355 sample estimates: p 0.48 Получается, что отвергнуть нулевую гипотезу о равенстве пропор- ций мы не можем, слишком велико p-value (гораздо больше «положен- ных» 0.05 и даже «либеральных» 0.1). Процент же проголосовавших может лежать в интервале от 42% до 54%! Итак, по результатам опроса нельзя сказать, что кандидат А победил. А вот так можно рассчитать количество человек, которых надо опро- сить, чтобы быть уверенным в том, что 48% и 52% и вправду отражают реальную ситуацию (генеральную совокупность): > power.prop.test(p1=0.48, p2=0.52, power=0.8) Two-sample comparison of proportions power calculation n = 2451.596 p1 = 0.48 p2 = 0.52 sig.level = 0.05 power = 0.8 alternative = two.sided NOTE: n is number in *each* group Получается, что опросить надо было примерно 5 тысяч человек! Мы использовали здесь так называемый power-тест, который часто применяется для прогнозирования эксперимента. Мы задали значение power = 0.8 , потому что именно такое значение является общеприня- тым порогом значимости (оно соответствует p-value < 0.1). Чуть бо- лее подробно про мощность (power) и ее связь с другими характеристи- ками теста можно прочитать в следующей главе. Глава 5 Анализ связей: двумерные данные Здесь речь пойдет о том, как работать с двумя выборками. Если у нас два набора цифр, то первое, что приходит в голову,— сравнить их. Для этого нам потребуются статистические тесты. 5.1. Что такое статистический тест Настала пора подробнее познакомиться с ядром статистики — тес- тами и гипотезами. В предыдущей главе мы уже коротко объясняли, как строятся статистические гипотезы, но если анализ одной выборки можно сделать, не вдаваясь подробно в их смысл, то в анализе связей без понимания этих принципов не обойтись. 5.1.1. Статистические гипотезы Из первой главы мы знаем, что статистическая выборка должна быть репрезентативной (то есть адекватно характеризовать генераль- ную совокупность, или, как ее еще называют, популяцию). Но как же мы можем знать, репрезентативна ли выборка, если мы не исследовали всю генеральную совокупность? Этот логический тупик называют пара- доксом выборки. Хотя мы и обеспечиваем репрезентативность выборки соблюдением двух основных принципов ее создания (рандомизации и повторности), неопределенность все же остается. Кроме того, если мы принимаем вероятностную точку зрения на происхождение наших дан- ных (которые получены путем случайного выбора), то все дальнейшие суждения, основанные на этих данных, будут иметь вероятностный ха- рактер. Таким образом, мы никогда не сможем на основании нашей выборки со стопроцентной уверенностью судить о свойствах генераль- ной совокупности! Мы можем лишь выдвигать гипотезы и вычислять их вероятность Великие философы науки (например, Карл Поппер) постулировали, что мы ничего не можем доказать, мы можем лишь что-нибудь опро- вергнуть. Действительно, если мы соберем 1000 фактов, подтверждаю- щих какую-нибудь теорию, это не будет значить, что мы ее доказали. |