А. Б. Шипунов, Е. М. Балдин, П. А. Волкова, А. И. Коробейников, С. А. Назарова
 Скачать 3.04 Mb.
|
|
. Создавая матрицу ma, мы задали byrow=TRUE , то есть указали, что элементы вектора будут соединяться в матрицу построчно. Если бы мы указали byrow=FALSE, то получили бы точно такую же матрицу, как mb: > ma <- matrix(m, ncol=2, byrow=FALSE) > ma [,1] [,2] [1,] 1 3 [2,] 2 4 Ответ к задаче про сортировку . Для того чтобы что-то сделать с колонками, надо использовать квадратные скобки с запятой в центре, а команды помещать справа от запятой. Логично использовать ту же команду order(): > d.sorted <- d[order(d$sex, d$height), ] > d.sorted[, order(names(d.sorted))] height sex size weight Женя 162 female S 68 Катя 165 female S 59 Вася 168 male M 82 Жора 172.5 male L 72 Коля 174 male L 69 Петя 188 male XL 93 Саша 192 male XXL 87 Заметьте, что нельзя просто написать order() после запятой, эта команда должна выдать новый порядок колонок. Поэтому мы «скорми- ли» ей имена этих колонок, которые для таблицы данных всегда можно получить, вызвав команду names(). Кстати, вместо order() здесь подо- шла бы и команда sort(), потому что нам надо отсортировать всего лишь один вектор. Данные гадания на ромашке . Вот они (проценты от общего ко- личества исходов): 72 Типы данных > rev(sort(romashka)) плюнет к сердцу прижмет поцелует 24 17 16 не любит любит к черту пошлет 15 15 13 (Мы использовали rev(), потому что sort() сортирует по возрас- танию.) Как видно, гадание было не очень-то удачное... Глава 4 Великое в малом: одномерные данные Теперь, наконец, можно обратиться к статистике. Начнем с самых элементарных приемов анализа — вычисления общих характеристик одной-единственной выборки. 4.1. Как оценивать общую тенденцию У любой выборки есть две самые общие характеристики: центр (центральная тенденция) и разброс (размах). В качестве центра ча- ще всего используются среднее и медиана, а в качестве разброса — стандартное отклонение и квартили. Среднее отличается от медиа- ны прежде всего тем, что оно хорошо работает в основном тогда, когда распределение данных близко к нормальному (мы еще поговорим об этом ниже). Медиана не так зависит от характеристик распределения, как говорят статистики, она более робастна (устойчива). Понять разни- цу легче всего на таком примере. Возьмем опять наших гипотетических сотрудников. Вот их зарплаты (в тыс. руб.): > salary <- c(21, 19, 27, 11, 102, 25, 21) Разница в зарплатах обусловлена, в частности, тем, что Саша — экспедитор, а Катя — глава фирмы. > mean(salary); median(salary) [1] 32.28571 [1] 21 Получается, что из-за высокой Катиной зарплаты среднее гораз- до хуже отражает «типичную», центральную зарплату, чем медиана. Отчего же так получается? Дело в том, что медиана вычисляется со- вершенно иначе, чем среднее. Медиана — это значение, которое отсекает половину упорядоченной выборки. Для того чтобы лучше это показать, вернемся к тем двум векторам, на примере которых в предыдущей главе было показано, как присваиваются ранги: 74 Великое в малом: одномерные данные > a1 <- c(1,2,3,4,4,5,7,7,7,9,15,17) > a2 <- c(1,2,3,4,5,7,7,7,9,15,17) > median(a1) [1] 6 > median(a2) [1] 7 В векторе a1 всего двенадцать значений, то есть четное число. В этом случае медиана — среднее между двумя центральными числами. У век- тора a2 все проще, там одиннадцать значений, поэтому для медианы просто берется середина. Кроме медианы, для оценки свойств выборки очень полезны квар- тили , то есть те значения, которые отсекают соответственно 0%, 25%, 50%, 75% и 100% от всего распределения данных. Если вы читали преды- дущий абзац внимательно, то, наверное, уже поняли, что медиана — это просто третий квартиль (50%). Первый и пятый квартили — это соответственно минимум и максимум, а второй и четвертый кварти- ли используют для робастного вычисления разброса (см. ниже). Мож- но понятие «квартиль» расширить и ввести специальный термин для значения, отсекающего любой процент упорядоченного распределения (не обязательно по четвертям),— это называется «квантиль». Кванти- ли используются, например, при анализе данных на нормальность (см. ниже). Для характеристики разброса часто используют и параметрическую величину — стандартное отклонение. Широко известно «правило трех сигм», которое утверждает, что если средние значения двух выборок различаются больше чем на тройное стандартное отклонение, то эти выборки разные, то есть взяты из разных генеральных совокупностей. Это правило очень удобно, но, к сожалению, подразумевает, что обе выборки должны подчиняться нормальному распределению. Для вы- числения стандартного отклонения в R предусмотрена функция sd(). Кроме среднего и медианы, есть еще одна центральная характерис- тика распределения, так называемая мода, самое часто встречающееся в выборке значение. Мода применяется редко и в основном для но- минальных данных. Вот как посчитать ее в R (мы использовали для подсчета переменную sex из предыдущей главы): > sex <- c("male", "female", "male", "male", "female", "male", + "male") > t.sex <- table(sex) > mode <- t.sex[which.max(t.sex)] > mode male Как оценивать общую тенденцию 75 5 Таким образом, мода нашей выборки — male. Часто стоит задача посчитать среднее (или медиану) для целой таб- лицы данных. Есть несколько облегчающих жизнь приемов. Покажем их на примере встроенных данных trees: > attach(trees) # Первый способ > mean(Girth) [1] 13.24839 > mean(Height) [1] 76 > mean(Volume/Height) [1] 0.3890012 > detach(trees) > with(trees, mean(Volume/Height)) # Второй способ [1] 0.3890012 > lapply(trees, mean) # Третий способ $Girth [1] 13.24839 $Height [1] 76 $Volume [1] 30.17097 Первый способ (при помощи attach()) позволяет присоединить ко- лонки таблицы данных к списку текущих переменных. После этого к переменным можно обращаться по именам, не упоминая имени табли- цы. Важно не забыть сделать в конце detach(), потому что велика опасность запутаться в том, что вы присоединили, а что — нет. Если присоединенные переменные были как-то модифицированы, на самой таблице это не скажется. Второй способ, в сущности, аналогичен первому, только присоедине- ние происходит внутри круглых скобок функции with(). Третий способ использует тот факт, что таблицы данных — это списки из колонок. Для строк такой прием не сработает, надо будет запустить apply(). (Если вам пришел в голову четвертый способ, то напоминаем, что цикличе- ские конструкции типа for в R без необходимости не приветствуются). Стандартное отклонение, дисперсия (его квадрат) и так называемый межквартильный разброс вызываются аналогично среднему: > sd(salary); var(salary); IQR(salary) [1] 31.15934 76 Великое в малом: одномерные данные [1] 970.9048 [1] 6 Последнее выражение, дистанция между вторым и четвертым квар- тилями IQR (или межквартильный разброс), робастен и лучше подхо- дит для примера с зарплатой, чем стандартное отклонение. Применим эти функции к встроенным данным trees: > attach(trees) > mean(Height) [1] 76 > median(Height) [1] 76 > sd(Height) [1] 6.371813 > IQR(Height) [1] 8 > detach(trees) Видно, что для деревьев эти характеристики значительно ближе друг к другу. Разумно предположить, что распределение высоты де- ревьев близко к нормальному. Мы проверим это ниже. В наших данных по зарплате — всего 7 цифр. А как понять, есть ли какие-то «выдающиеся» цифры, типа Катиной зарплаты, в данных большого, «тысячного» размера? Для этого есть графические функ- ции. Самая простая — так называемый «ящик-с-усами», или боксплот. Для начала добавим к нашим данным еще тысячу гипотетических ра- ботников с зарплатой, случайно взятой из межквартильного разброса исходных данных (рис. 9): > new.1000 <- sample((median(salary) - IQR(salary)) : + (median(salary) + IQR(salary)), 1000, replace=TRUE) > salary2 <- c(salary, new.1000) > boxplot(salary2, log="y") Это интересный пример еще и потому, что в нем впервые пред- ставлена техника получения случайных значений. Функция sample() способна выбирать случайным образом данные из выборки. В данном случае мы использовали replace=TRUE, поскольку нам нужно было вы- брать много чисел из гораздо меньшей выборки. Если писать на R ими- тацию карточных игр (а такие программы написаны!), то надо исполь- зовать replace=FALSE, потому что из колоды нельзя достать опять ту же самую карту. Кстати говоря, из того, что значения случайные, сле- дует, что результаты последующих вычислений могут отличаться, если Как оценивать общую тенденцию 77 20 50 10 0 Рис. 9. «Ящик-с-усами», или боксплот их воспроизвести еще раз, поэтому ваш график может выглядеть чуть иначе. Но вернемся к боксплоту. Как видно, Катина зарплата представ- лена высоко расположенной точкой (настолько высоко, что нам даже пришлось вписать параметр log="y", чтобы нижележащие точки ста- ли видны лучше). Сам бокс, то есть главный прямоугольник, ограничен сверху и снизу квартилями, так что высота прямоугольника — это IQR. Так называемые «усы» по умолчанию обозначают точки, удаленные на полтора IQR. Линия посередине прямоугольника — это, как легко до- гадаться, медиана. Точки, лежащие вне «усов», рассматриваются как выбросы и поэтому рисуются отдельно. Боксплоты были специально придуманы известным статистиком Джоном Тьюки, для того чтобы быстро, эффективно и устойчиво отражать основные робастные харак- теристики выборки. R может рисовать несколько боксплотов сразу (то есть эта команда векторизована, см. результат на рис. 10): > boxplot(trees) 78 Великое в малом: одномерные данные Girth Height Volume 20 40 60 80 Рис. 10. Три боксплота, каждый отражает одну колонку из таблицы данных Есть две функции, которые связаны с боксплотами. Функция quantile() по умолчанию выдает все пять квартилей, а функция fivenum() — ос- новные характеристики распределения по Тьюки. Другой способ графического изображения — это гистограмма, то есть линии столбиков, высота которых соответствует встречаемости дан- ных, попавших в определенный диапазон (рис. 11): > hist(salary2, breaks=20, main="") В нашем случае hist() по умолчанию разбивает переменную на 10 интервалов, но их количество можно указать вручную, как в предло- женном примере. Численным аналогом гистограммы является функция cut() . При помощи этой функции можно выяснить, сколько данных ка- кого типа у нас имеется: > table(cut(salary2, 20)) (10.9,15.5] (15.5,20] (20,24.6] (24.6,29.1] (29.1,33.7] Как оценивать общую тенденцию 79 20 40 60 80 100 0 10 0 20 0 30 0 40 0 Рис. 11. Гистограмма зарплат 1007 гипотетических сотрудников 76 391 295 244 0 (33.7,38.3] 0 (38.3,42.8] (42.8,47.4] (47.4,51.9] (51.9,56.5] (56.5,61.1] 0 0 0 0 0 (61.1,65.6] 0 (65.6,70.2] (70.2,74.7] (74.7,79.3] (79.3,83.9] (83.9,88.4] 0 0 0 0 0 (88.4,93] (93,97.5] (97.5,102] 0 0 1 Есть еще две графические функции, «идеологически близкие» к гис- тограмме. Во-первых, это stem() — псевдографическая (текстовая) гис- тограмма: > stem(salary, scale=2) The decimal point is 1 digit(s) to the right of the | 80 Великое в малом: одномерные данные 1 | 19 2 | 1157 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 2 Это очень просто — значения данных изображаются не точками, а цифрами, соответствующими самим этим значениям. Таким образом, видно, что в интервале от 10 до 20 есть две зарплаты (11 и 19), в ин- тервале от 20 до 30 — четыре и т. д. Другая функция тоже близка к гистограмме, но требует гораздо более изощренных вычислений. Это график плотности распределения (рис. 12): > plot(density(salary2, adjust=2), main="") > rug(salary2) (Мы использовали «добавляющую» графическую функцию rug(), чтобы выделить места с наиболее высокой плотностью значений.) По сути, перед нами сглаживание гистограммы — попытка превра- тить ее в непрерывную гладкую функцию. Насколько гладкой она бу- дет, зависит от параметра adjust (по умолчанию он равен единице). Результат сглаживания называют еще графиком распределения. Кроме боксплотов и различных графиков «семейства» гистограмм, в R много и других одномерных графиков. График-«улей», например, отражает не только плотность распределения значений выборки, но и то, как расположены сами эти значения (точки). Для того чтобы по- строить график-улей, потребуется загрузить (а возможно, еще и уста- новить сначала) пакет beeswarm. После этого можно поглядеть на сам «улей» (рис. 13): > library("beeswarm") > beeswarm(trees) > boxplot(trees, add=TRUE) Мы здесь не просто построили график-улей, но еще и добавили ту- да боксплот, чтобы стали видны квартили и медиана. Для этого нам понадобился аргумент add=TRUE. Ну и, наконец, самая главная функция, summary(): Как оценивать общую тенденцию 81 20 40 60 80 100 0. 00 0. 02 0. 04 0. 06 0. 08 Рис. 12. Плотность распределения зарплат 1007 гипотетических сотруд- ников > lapply(list(salary, salary2), summary) [[1]] Min. 1st Qu. Median Mean 3rd Qu. 11.00 20.00 21.00 32.29 26.00 Max. 102.00 [[2]] Min. 1st Qu. Median Mean 3rd Qu. 11.00 18.00 21.00 21.09 24.00 Max. 102.00 Фактически она возвращает те же самые данные, что и fivenum() с добавлением среднего значения (Mean). Заметьте, кстати, что у обеих «зарплат» медианы одинаковы, тогда как средние существенно отлича- ются. Это еще один пример неустойчивости средних значений — ведь с 82 Великое в малом: одномерные данные 20 40 60 80 Girth Height Volume Girth Height Volume 20 40 60 80 Рис. 13. График-«улей» с наложенными боксплотами для трех харак- теристик деревьев добавлением случайно взятых «зарплат» вид распределения не должен был существенно поменяться. Функция summary() — общая, и по законам объект-ориентированно- го подхода она возвращает разные значения для объектов разного типа. Вы только что увидели, она работает для числовых векторов. Для спис- ков она работает немного иначе. Вывод может быть, например, таким (на примере встроенных данных attenu о 23 землетрясениях в Кали- форнии): > summary(attenu) event mag station dist Min. : 1.00 Min. :5.000 117 : 5 Min. : 0.50 1st Qu.: 9.00 1st Qu.:5.300 1028 : 4 1st Qu.: 11.32 Median :18.00 Median :6.100 113 : 4 Median : 23.40 Mean :14.74 Mean :6.084 112 : 3 Mean : 45.60 3rd Qu.:20.00 3rd Qu.:6.600 135 : 3 3rd Qu.: 47.55 Ошибочные данные 83 Max. :23.00 Max. :7.700 (Other):147 Max. :370.00 NA’s : 16 accel Min. :0.00300 1st Qu.:0.04425 Median :0.11300 Mean :0.15422 3rd Qu.:0.21925 Max. :0.81000 Переменная station (номер станции наблюдений) — фактор, и к тому же с пропущенными данными, поэтому отображается иначе. Перед тем как завершить рассказ об основных характеристиках вы- борки, надо упомянуть еще об одной характеристике разброса. Для сравнения изменчивости признаков (особенно таких, которые измере- ны в разных единицах измерения) часто применяют безразмерную ве- личину — коэффициент вариации. Это просто отношение стандартного отклонения к среднему, взятое в процентах. Вот так можно сравнить коэффициент вариации для разных признаков деревьев (встроенные данные trees): > 100*sapply(trees, sd)/colMeans(trees) Girth Height Volume 23.686948 8.383964 54.482331 Здесь мы для быстроты применили sapply() — вариант lapply() с упрощенным выводом, и colMeans(), которая просто вычисляет сред- нее для каждой колонки. Сразу отметим, что подобных функций в R несколько. Например, широко используются функции colSums() и rowSums() , которые выдают итоговые суммы соответственно по колон- кам и по строкам (главная функция электронных таблиц!). Есть, разу- меется, еще и rowMeans(). 4.2. Ошибочные данные Способность функции summary() указывать пропущенные данные, максимумы и минимумы служит очень хорошим подспорьем на самом раннем этапе анализа данных — проверке качества. Предположим, у нас есть данные, набранные с ошибками, и они находятся в директории data внутри текущей директории: > dir("data") [1] "errors.txt" ... 84 Великое в малом: одномерные данные > err <- read.table("data/errors.txt", h=TRUE, sep="\t", + stringsAsFactors=TRUE) > str(err) ’data.frame’: 7 obs. of 3 variables: $ AGE : Factor w/ 6 levels "12","22","23",..: 3 4 3 5 1 6 2 $ NAME : Factor w/ 6 levels "","John","Kate",..: 2 3 1 4 5 6 2 $ HEIGHT: num 172 163 161 16.1 132 155 183 > summary(err) AGE NAME HEIGHT 12:1 :1 Min. : 16.1 22:1 John :2 1st Qu.:143.5 23:2 Kate :1 Median :161.0 24:1 Lucy :1 Mean :140.3 56:1 Penny:1 3rd Qu.:167.5 a :1 Sasha:1 Max. :183.0 Обработка начинается с проверки наличия нужного файла. Кроме команды summary(), здесь использована также очень полезная команда str() . Как видно, переменная AGE (возраст) почему-то стала фактором, и summary() показывает, почему: в одну из ячеек закралась буква a. Кроме того, одно из имен пустое, скорее всего, потому, что в ячейку забыли поставить NA. Наконец, минимальный рост — 16.1 см! Такого не бывает обычно даже у новорожденных, так что можно с уверенностью утверждать, что наборщик просто случайно поставил точку. 4.3. Одномерные статистические тесты Закончив разбираться с описательными статистиками, перейдем к простейшим статистическим тестам (подробнее тесты рассмотрены в следующей главе). Начнем с так называемых «одномерных», которые позволяют проверять утверждения относительно того, как распределе- ны исходные данные. Предположим, мы знаем, что средняя зарплата в нашем первом при- мере — около 32 тыс. руб. Проверим теперь, насколько эта цифра до- стоверна: > t.test(salary, mu=mean(salary)) One Sample t-test data: salary t = 0, df = 6, p-value = 1 alternative hypothesis: true mean is not equal to 32.28571 95 percent confidence interval: |