Описательные статистики

Рассмотрим небольшую выборку из шести наблюдений и сохраним её в вектор v. Посмотрим на описательные статистики по этой выборке:

v <- c(2, 3, 5, 5, 7, 8) 
summary(v)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     2.0     3.5     5.0     5.0     6.5     8.0

Минимальное значение 2, максимальное 8. Среднее значение равно 5, медиана тоже равна 5. Нижний квартиль равен 3.5 (1st Qu.), верхний квартиль равен 6.5 (3rd Qu.). Проинтерпретируем:

Примечание: результаты, полученные в R, а именно значения квартилей, могут отличаться от того, что мы считали вручную, так как в R по умолчанию используется другой алгоритм расчёта.

Теперь посмотрим, как можно вызвать отдельные описательные статистики (по вектору или по столбцу в таблице).

min(v)  # минимум
## [1] 2
max(v) # максимум
## [1] 8
vrange <- max(v) - min(v)  # размах
vrange
## [1] 6
mean(v) # среднее
## [1] 5
median(v) # медиана
## [1] 5

Квантили заслуживают отдельного внимания. Для начала найдём квантиль уровня 0.2: укажем название вектора и через запятую - уровень квантиля.

quantile(v, 0.2)
## 20% 
##   3

Получается, 20% наблюдений в выборке не превышают значение 3. Внутри функции quantile() можно указать сразу несколько уровней, только для этого нужно оформить их в виде вектора:

quantile(v, c(0.1, 0.2, 0.3)) 
## 10% 20% 30% 
## 2.5 3.0 4.0

Получили квантили уровней 0.1, 0.2, 0.3. Если уровней много, и все они отличаются друг от друга на одно и то же число, то их можно задать в виде последовательности. Для этого понадобится функция seq() (от английского sequence), внутри которой нужно указать начало отрезка, конец отрезка и шаг, с которым мы по этому отрезку движемся:

quantile(v, seq(from=0, to=1, by=0.1)) # от 0 до 1 с шагом 0.1
##   0%  10%  20%  30%  40%  50%  60%  70%  80%  90% 100% 
##  2.0  2.5  3.0  4.0  5.0  5.0  5.0  6.0  7.0  7.5  8.0

То, что мы получили выше - это процентили (квантили, только уровень указан в процентах), а если говорить ещё точнее, децили, то есть процентили, уровень которых делится на 10: 10%, 20%, 30%, и так далее. Теперь вернёмся к квартилям.

quantile(v, c(0.25, 0.75)) # нижний и верхний квартили
## 25% 75% 
## 3.5 6.5
quantile(v, 0.75) - quantile(v, 0.25)  # межквартильный размах
## 75% 
##   3
IQR(v) # межквартильный размах - готовая функция R
## [1] 3

Посчитаем выборочную дисперсию и стандартное отклонение:

var(v)  # variance - дисперсия
## [1] 5.2
sd(v) # standard deviation - стандартное отклонение
## [1] 2.280351

А сейчас перейдём к более красочным вещам - к графикам.

Гистограммы (histograms)

Чтобы строить более осмысленные графики, вернёмся к реальным данным и загрузим таблицу из файла survey01.csv с результатами опроса первокурсников ОП «Психология». Для разнообразия загрузим файл по ссылке - просто скопируем её и вставим в кавычках в функцию read.csv():

dat <- read.csv("http://math-info.hse.ru/f/2018-19/psych-ms/survey01.csv",
                dec = ",")

Теперь построим простенькую гистограмму для переменной height (рост студентов):

hist(dat$height)

По умолчанию гистограмма прорисовывается без заливки, получается белой. Исправим это: добавим опцию col и укажем название цвета в кавычках. Список всех цветов в R можно посмотреть здесь.

hist(dat$height, col = "hotpink")

Что мы можем сказать о распределении роста по этой гистограмме? Во-первых, что распределение роста похоже на нормальное (что никак не противоречит уже существующим биологическим исследованиям), только в отличие от нормального, у этого распределения есть более длинный «хвост» справа. Другими словами, большинство студентов имеет рост, несильно отличающийся от среднего, однако есть студенты, у которых рост довольно высокий, около 190 см. Во-вторых, студентов, у которых рост достаточно маленький, немного.

Теперь давайте приведём график к приличному виду: добавим заголовок и подписи к осям.

# main - заголовок
# xlab - подпись для оси x
# ylab - подпись для оси y
hist(dat$height, col = "hotpink", 
  main = "Моя первая розовая гистограмма", 
  xlab = "Рост (в см)",
  ylab = "Число студентов")

Обратите внимание: все названия, так как они являются текстовыми, вводятся в кавычках.

Чему равен шаг гистограммы на графике выше? Пяти, так как на отрезке в 10 см помещаются ровно два столбца одинаковой ширины. К сожалению, в базовом варианте гистограммы, создаваемой с помощью hist(), изменить величину шага вручную не так просто. Но зато есть опция breaks, которая приблизительно соответствует числу «перегородок», которые необходимо поставить, чтобы получить желаемое число столбцов. То есть, если мы хотим получить пять столбцов, нам потребуется четыре «перегородки»:

hist(dat$height, col = "hotpink", 
  main = "Моя первая розовая гистограмма", 
  xlab = "Рост (в см)", 
  breaks = 4)

Более продвинутые описательные статистики

Теперь посмотрим на более подробную выдачу R с описательными статистиками. Чтобы это сделать, нам понадобится библиотека psych, которая содержит набор функций, часто используемых в психометрических исследованиях. Установим её:

install.packages("psych")

Внимание: если R не хочет устанавливать библиотеку и пишет что-то в стиле is not writable, можно запустить RStudio от имени администратора, обычно помогает. Если нет, стоит прописать путь к папке, куда должны устанавливаться библиотеки.

Обратимся к библиотеке через library(), чтобы R понимал, откуда ему брать последующие функции. Если это не сделать, в ответ на попытки использовать функции из этой библиотеки, R будет писать could not find function....

library(psych)

Теперь запросим описательные статистики с помощью функции describe():

describe(dat$height)
##    vars   n   mean   sd median trimmed  mad min max range skew kurtosis
## X1    1 108 168.76 8.59 168.75  168.55 8.52 145 190    45 0.16    -0.09
##      se
## X1 0.83

Что есть что?

Подробнее про некоторые статистики.

Усечённое среднее, среднее по цензурированной выборке

Коэффициент асимметрии

Коэффициент эксцесса

Библиотека psych удобна тем, что она содержит функцию describeBy(), которая позволяет выводить описательные статистики по группам. Нет необходимости отфильтровывать нужные строки и сохранять их в отдельные датасеты, можно просто указать группирующую переменную. Выведем описательные статистики для переменной height отдельно для групп по R и SPSS:

describeBy(dat$height, group = dat$soft)
## 
##  Descriptive statistics by group 
## group: R
##    vars  n   mean   sd median trimmed mad min max range skew kurtosis   se
## X1    1 93 168.75 8.75    169  168.55 8.9 145 190    45 0.13    -0.11 0.91
## -------------------------------------------------------- 
## group: SPSS
##    vars  n   mean   sd median trimmed mad min max range skew kurtosis se
## X1    1 15 168.87 7.74    168  168.54 8.9 156 186    30 0.36    -0.56  2

Очень удобно!

Ящики с усами (box plots)

Продолжим обсуждать визуализацию распределений данных и рассмотрим ещё один тип графика - ящик с усами (box plot или box-and-whiskers plot). Давайте сначала его построим, а потом посмотрим, что он иллюстрирует.

boxplot(dat$height)

Что показывает данный график? Квартили и минимальное/максимальное значения или границы типичных значений. Нижняя и верхняя граница «ящика» соответствуют нижнему и верхнему квартилю (Q1 и Q3), толстая горизонтальная линия внутри «ящика» соответствует медиане. Горизонтальная линия не всегда будет лежать в середине «ящика», в нашем случае так совпало, так как распределение роста студентов примерно симметричное. Если бы распределение было скошено (слишком много людей высокого роста или наоборот, низкого), то медиана была бы ближе к нижнему или верхнему квартилю, то есть к нижней или верхней границе «ящика». «Усы» графика (похожи на антенны сверху и снизу) строятся по-разному, и их построение зависит от наличия нетипичных/нехарактерных наблюдений или выбросов (outliers). Выбросы на графике изображаются отдельными точками, которые находятся всегда за рамками «усов». В нашем случае есть один выброс, в нижней части графика, то есть нетипично маленькое значение.

Если выбросов ни с какой стороны нет, то границы «усов» совпадают с минимальным и максимальным значением в выборке. Если выбросы есть с обеих сторон (слишком маленькие и слишком большие значения), то границы «усов» совпадают с границами типичных значений, которые вычисляются следующим образом:

\[[\text{Q}_1 - 1.5 \times \text{IQR} ; \text{Q}_3 + 1.5 \times \text{IQR}].\] Все значения, которые выходят за границы этого интервала, считаются нетипичными значениями. Иногда среди нетипичных значений выделяют совсем нехарактерные (в SPSS они обозначаются звёздочкой) - значения, которые выходят за границы следующего интервала:

\[[\text{Q}_1 - 3 \times \text{IQR} ; \text{Q}_3 + 3 \times \text{IQR}].\]

Если выбросы есть только сверху, слишком большие значения, то нижний «ус» совпадает с минимальным значением, а верхний - со значением \(\text{Q}_3 + 1.5 \times \text{IQR}\). Если выбросы есть только снизу, слишком маленькие значения, то нижний «ус» соответствует значению \(\text{Q}_1 - 1.5 \times \text{IQR}\), а верхний - максимальному.

Дополнительно: скрипичные диаграммы (violin plots)

Дополнительно рассмотрим ещё один тип графика, менее типичный, который называется скрипичной диаграммой (на английском есть целый ряд терминов: violin plot, bean plot, vase plot, кто-то видит скрипки, кто-то бобы, а кто-то вазы). Давайте опять сначала на него посмотрим, а потом разберёмся, что к чему. Чтобы построить такой график, нам потребуется библиотека vioplot.

install.packages("vioplot")

Удалим пропущенные значения и построим скрипичную диаграмму для баллов за ЕГЭ по биологии:

library(vioplot)

dat <- na.omit(dat)
vioplot(dat$bio)

Как можно догадаться, внутри этого графика, внутри «скрипки», находится «ящик» с усами, на котором белой точкой обозначена медиана. А что снаружи? Этот график состоит из двух симметричных половинок, каждая из которых представляет собой график плотности распределения, о котором можно думать как о сглаженной гистограмме. Этот график хорошо характеризует форму распределения; по нему видно, например, что у нас есть две группы студентов («два горба»): те, у которых баллы достаточно средние, около 75, и те, у кого баллы высокие, около 90. Чтобы было совсем ясно, из чего этот график складывается, построим график плотности, определённым образом «сглаженную» гистограмму:

d <- density(dat$bio)
plot(d, xlim = c(60, 100))
polygon(d, col="magenta")

Скрипичные диаграммы довольно часто используются для визуализации данных, особенно, если исследователя интересует форма распределения и распределение отличается от нормального, но стоит иметь в виду, что такой график корректно строить в случае, когда выборка достаточно большая. Если выборка маленькая, то сглаживание, которое используется для построения графика, сработает очень грубо, и картина исказится. (Попробуйте построить гистограмму с маленьким шагом для маленькой выборки вручную и «сгладить» - схематично обвести её контур, вы поймёте, о чём идёт речь).