Мой маленький плот вовсе не так уж плох

Дисклеймер: этот конспект был написан от первого лица, если через пару минут прочтения у вас сложится впечатление что вы сидите в одной комнате с Иваном Поздняковым, это нормально!

Введение

Меня зовут Иван Поздняков, я преподаватель статистики с 10-летним стажем в Бластим, Вышке, РЭШ, Физтехе и Сириусе. Сейчас живу в Киле, Германия.

По образованию я психофизиолог, но со временем понял, что для качественной науки нужно уметь работать с данными, а в кнопочных программах типа Excel или SPSS это делать сложно. Перепробовал Matlab и Python, но остановился на R — считаю его лучшим инструментом для анализа данных. А ggplot2, которым мы сегодня будем активно пользоваться — это вообще лучший инструмент для визуализации.

Этот вебинар — демоверсия примерно середины нашего курса: предполагается, что вы уже немного знакомы с R и RStudio. Буду всё пояснять, но если что-то покажется сложным — для глубокого понимая нужен полноценный курс.

Подготовка данных

В первую очередь мы будем плотно работать с пакетом tidyverse — это наш основной пакет, точнее этот пакет с пакетами, где есть ggplot2 для визуализации данных, где есть dplyr, tidyr и все прочее.

В качестве датасета мы будем использовать пингвинов. Раньше были распространены ирисы, но потом ирисы отменили по политическим причинам (можно прочитать тут), теперь пингвины. Недавно пингвины появились внутри самого R: если у вас достаточно свежая версия R, где-то годовалой давности, то можно просто запустить penguins и у вас будет готовый датасет. У меня версия R более старая, поэтому я воспользуюсь пакетом palmerpenguins, чтобы вытащить этот датасет.

# Подключаем tidyverse
library(tidyverse)

# Подключаем пакет с пингвинами
library(palmerpenguins)

# Загружаем данные
data("penguins")

# Смотрим на данные
View(penguins)

Есть их три вида пингвинов на трех островах, там мальчики и девочки. Тут есть пропущенные значения. Давайте мы их сразу уберем с помощью функции drop_na.

# Удаляем пропущенные значения
my_penguins <- drop_na(penguins)

# Проверяем
View(my_penguins)

Теперь пропущенных значений у нас нет.

Теория квартилей и мусорные журналы

Мы начнем с box-плота. Разберем на примере научных журналов. Многие слышали, что зарплата в Вышке у сотрудников не очень большая, зато дают надбавки за научные статьи. Но казалось бы, эту систему легко обойти, если мы будем просто публиковаться в каких-то мусорных журналах, которые — платишь деньги, а они публикуют твою статью. Но не все так просто, потому что Вышка хочет, чтобы сотрудники печатались именно в самых топовых журналах. И как они определяют эту топовость? Есть большие базы: Web of Science или Scopus, которые оценивают импакт-фактор журнала — это средняя цитируемость на статью. И мы можем отранжировать все научные журналы по импакт-фактору.

Самый высокий импакт-фактор будет у журналов: Science, Nature, а самый низкий — у какого-нибудь вестника Балашихинского университета. Журналы ранжируются по импакт-фактору и делятся на 4 группы. Топ-25% попадают в первый квартиль — это самые топовые журналы, за них дают большую надбавку. Следующие 25% — это второй квартиль. За эти статьи тоже дают денежки, но поменьше. Следующие 25%, то есть от топ-50% до топ-75% — это третий квартиль, за него уже ничего не дают, но хотя бы учитывают, чтобы вас не уволили. А четвертый квартиль — это 25% наименее цитируемых журналов, мусорные журналы, где можно опубликоваться просто заплатив деньги.

Вот квартиль — это как раз разбиение по 25%. В статистике порядок немного другой: первый квартиль — это граница, отделяющая нижние 25% значений. Второй квартиль — это граница 50%, она же медиана. Третий квартиль — граница 75%. Иногда минимальное значение называют нулевым квартилем, максимальное — четвертым квартилем, но это уже эзотерика.

Давайте проверим на наших данных:

# Смотрим описательные статистики для массы тела
summary(my_penguins$body_mass_g)

# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 2700    3550    4050    4207    4775    6300

NB! В новых версиях R колонка body_mass_g называется body_mass

Вот тут можно увидеть минимальное значение, первый квартиль, медиана, среднее, третий квартиль, максимальное значение. Может возникнуть вопрос: первый квартиль есть, третий квартиль есть, куда делся второй квартиль? А это медиана и есть второй квартиль.

Строим первый бокс-плот

Давайте теперь нарисуем самый простой бокс-плот для массы тела всех пингвинов.

# Самый простой бокс-плот
ggplot(data = my_penguins) +
  aes(x = "Все пингвины", y = body_mass_g) +
  geom_boxplot()

Здесь x = "Все пингвины" — это просто подпись, потому что у нас одна группа. y = body_mass_g — это масса тела.

Запускаем — и у нас появляется ящик с усами. Видим ящик от первого квартиля до третьего, линию на уровне медианы и усы.

Сравнение нескольких групп

Обычно бокс-плоты используют, чтобы сравнивать распределения в разных группах. Давайте посмотрим, как масса тела различается у трех видов пингвинов.

# Бокс-плот по видам
ggplot(data = my_penguins) +
  aes(x = species, y = body_mass_g) +
  geom_boxplot()

Сразу видно, что пингвины Gentoo заметно тяжелее, а у Adelie и Chinstrap масса похожа, но у Chinstrap медиана чуть выше.

50 оттенков пингвинов

Для наглядности добавим цвета: fill = species — закрашивает ящики разными цветами в зависимости от вида.

# С цветом
ggplot(data = my_penguins) +
  aes(x = species, y = body_mass_g, fill = species) +
  geom_boxplot()

А также меняем тему: theme_bw() — это черно-белая тема, четкие линии.

# Добавляем тему
ggplot(data = my_penguins) +
  aes(x = species, y = body_mass_g, fill = species) +
  geom_boxplot() +
  theme_bw()

Как всё заиграло, а?

Что такое усы?

Вернемся к усам. С ними всё сложно. Многие думают, что усы тянутся до минимума и максимума. Но это не так. В классическом бокс-плоте по Тьюки (а ggplot2 делает именно такие) длина усов ограничена 1.5 интерквартильными размахами.

Интерквартильный размах (IQR) — это расстояние от Q1 до Q3. Верхний ус рисуется до максимального значения, которое не превышает Q3 + 1.5 * IQR. Нижний ус рисуется до минимального значения, которое не меньше Q1 - 1.5 * IQR. Все точки за пределами этих границ считаются выбросами и рисуются отдельными точками.

Давайте проверим на пингвинах вида Adelie: вы видите, что усы не доходят до абсолютного максимума, если есть выбросы.

# Отфильтруем только Adelie
adelie_mass <- my_penguins %>%
  filter(species == "Adelie") %>%
  pull(body_mass_g)

# Считаем квартили
q1 <- quantile(adelie_mass, 0.25)
q3 <- quantile(adelie_mass, 0.75)
iqr <- q3 - q1

# Границы усов
upper_whisker <- q3 + 1.5 * iqr
lower_whisker <- q1 - 1.5 * iqr

# Максимальное и минимальное значение в пределах границ
max_in_whisker <- max(adelie_mass[adelie_mass <= upper_whisker])
min_in_whisker <- min(adelie_mass[adelie_mass >= lower_whisker])

print(paste("Q1:", q1))
# "Q1: 3362.5"
print(paste("Q3:", q3))
# "Q3: 4000"
print(paste("IQR:", iqr))
# "IQR: 637.5"
print(paste("Верхняя граница уса:", upper_whisker))
# "Верхняя граница уса: 4956.25"
print(paste("Фактический верхний ус:", max_in_whisker))
# "Фактический верхний ус: 4775"

Кроме того, есть специальная функция, которая позволяет то же самое всё посчитать. Она называется quantile, потому что квартиль — это только один из видов квантилей. Есть еще процентили, децили, секстили, квинтили и так далее. Все они различаются только тем, на сколько интервалов мы делим: квартили — соответственно на 4, процентили — на 100, децили — на 10 и так далее. По умолчанию quantile считает именно квартили, и это как раз то, что нам и нужно.

Вот, можно даже разобрать на каком-то несложном примере, как это работает. На примере вот таких вот чисел: 1, 2, 4, 5, 8, 12, 14, 30.

Соответственно, 25% снизу — это будет первый квартиль, следующие 25% — это второй квартиль, следующие 25% — третий квартиль, последние 25% — это четвертый квартиль. Ну и соответственно границы между ними — это будут квартили как границы. То есть 3 — это будет Q1, медиана — это будет 6.5 получается, и граница между третьим и четвертым квартилями, то есть третий квартиль — это 13. Получается 3, 6.5, 13.

Давайте сопоставим наши расчеты с тем, что выдает функция quantile. Передадим функции quantile
вектор: 1, 2, 4, 5, 8, 12, 14, 30.

# Считаем квартили для нашего примера
quantile(c(1, 2, 4, 5, 8, 12, 14, 30))

# 0%  25%  50%  75% 100% 
# 1.0  3.5  6.5 12.5 30.0

Вот и получим немножко другие результаты. Почему?

9 видов квантилей. Чтоа?

На самом деле есть всего 9 типов квантилей в R, и тот, который мы сейчас сделали — это второй тип, а по дефолту используется седьмой. Разница на самом деле между ними не очень большая, особенно когда у вас достаточно много данных. Давайте уточним тип, чтобы полностью воспроизвести наши расчеты:

# Считаем квартили для нашего примера
quantile(c(1, 2, 4, 5, 8, 12, 14, 30), type = 2)

# 0%  25%  50%  75% 100% 
# 1.0  3  6.5 13 30.0

Вот, мы полностью воспроизвели: первый квартиль 3, второй квартиль (медиана) 6,5, третий квартиль 13 — вот как раз ровно то же самое получили.

Если у вас достаточно много данных, то разница между типами становится всё меньше и меньше. В основном они просто по-разному решают проблему, что делать с пограничными значениями: считать между ними просто среднее арифметическое или что-то более сложное.

Внутренность усиковой диаграммы

Итак, теперь, после того как мы ознакомились с квартилями, мы готовы разбирать, как устроены эти самые boxplot. Мы будем горизонтально рисовать boxplot. Вот эта нижняя и верхняя граница коробки — это первый и третий квартили:

Что мы можем сказать по этой коробке сразу? Что 50% значений нашей выборки находится внутри коробки, а 50% (по 25% слева и по 25% справа) находится вне этой коробки. В общем-то, это довольно полезно, когда вы смотрите на boxplot, понимать эту информацию.

Соответственно, медиана, то есть второй квартиль. Как видите, она далеко не обязательно находится посередине, и как раз указывает на асимметрию нашей выборки. Мы можем увидеть, что она чуть левее центра — действительно, у нас такой длинный жирный хвост справа получается.

Ну с этим всё более-менее понятно. А что касается усов — насчет усов на самом деле всё сложно, потому что, к сожалению, нет единого какого-то правила, единого закона, единого стандарта для boxplot, как именно рисовать усы. Есть более-менее стандартные варианты, которые можно прочитать на Википедии, которые используются в большинстве программ по умолчанию, в том числе в том же самом R. Именно его мы и разберем.

Довольно часто говорят, что усы обозначают минимальное и максимальное значение. Но тогда возникает вопрос: откуда появляются точки? Вот эти вот, которые называют многие "outliers" (выбросы). Я не очень люблю, когда их так называют, потому что это как будто наделяет их какой-то математической сущностью, которой они не являются.

Boxplot в базовом R

Давайте нарисуем нормальный boxplot с нашими данными в базовом R. Базовый R, мне кажется, для визуализации стандартно используется, когда нужно что-то быстренько посмотреть, не очень заморачиваться с украшательствами и настройками, а просто для исследования. Потому что базовые инструменты для визуализации в R — они простые, дубовые, но не очень красивые.

Давайте это сделаем. Тут используется синтаксис формул — вот этот значок тильда (~). Он находится там же, где у вас буква ё или backtick, скорее всего, слева сверху на клавиатуре. И это условно обозначает такое "равно" в уравнении: слева зависимая переменная, справа — какая-то независимая переменная, которой мы манипулируем, или какая-то группирующая переменная.

# Boxplot в базовом R
boxplot(body_mass_g ~ species, data = my_penguins)

...и в ggplot2

Давайте то же самое теперь нарисукм с помощью ggplot2. У ggplot гораздо более сложный синтаксис, чем у базового R. GG (grammar of graphics)— это специальный язык для описания графики, который был реализован в этом пакете. Это не просто пакет, который делает красивые картинки, а за ним стоит целая система, целая теория. Но про это у нас есть отдельно в курсе.

Сейчас просто про основные элементы синтаксиса. Мы инициализируем картинку функцией ggplot(). Можно туда сразу передать данные, которые мы будем использовать. И график у нас состоит из одного или нескольких слоев. Для каждого слоя у нас должны быть:

Данные
Геометрия — во что мы превращаем эти данные
Эстетики — которые соединяют данные и геометрию

А также могут быть еще статистические преобразования или position adjustment

Но пока начнем с простого варианта, когда у нас один слой. Идем к boxplot. Что нам нужно? Мы выбрали геометрию, у нас есть данные my_penguins. Нам нужно прописать эстетики, которые указывают, какие именно колонки в наших данных превращаются в какие особенности графика.

Для boxplot это будет группировка по x и по y. y — количественная переменная body_mass_g. Всё, что у нас как-то зависит от данных, все какие-то характеристики, которые зависят от данных, они прописываются внутри функции aes(). aes — сокращение от aesthetics. Эстетика — от греческого "то, как абстрактные идеи превращаются в конкретные чувственные ощущения". Поэтому они назвали это эстетикой — не потому что это делает красиво, а вот такая идея: как абстрактные данные превращаются в конкретные характеристики геометрии.

# Boxplot в ggplot2
ggplot(data = my_penguins) +
  aes(x = species, y = body_mass_g) +
  geom_boxplot()

Давайте сделаем небольшое изменение — чуть-чуть поуже сделаем, потому что вот такие широкие boxplot выглядят не очень симпатично.

# Boxplot с настройкой ширины
ggplot(data = my_penguins) +
  aes(x = species, y = body_mass_g) +
  geom_boxplot(width = 0.5)

Гистограмма

Давайте сразу нарисуем гистограмму по тем же самым данным. Важно обратить внимание: гистограмма — это принципиально одномерная визуализация. Внешне гистограмма похожа на барплоты, на столбики, но по сути сильно от них различается. Потому что здесь у нас есть просто одна переменная, числовая, которую дальше алгоритм как-то разделяет на эти интервалы, и всё. Нет никакой второй переменной!

# Гистограмма в ggplot2
ggplot(data = my_penguins) +
  aes(x = body_mass_g) +
  geom_histogram()

Если мы сделаем bins (количество интервалов), например, равным 4, то обратите внимание, не обязательно ровно 4 интервала получится.

# Гистограмма с 4 интервалами
ggplot(data = my_penguins) +
  aes(x = body_mass_g) +
  geom_histogram(bins = 4)

Там используется алгоритм, который задает интервалы таким образом, чтобы это число примерно соответствовало тому, которое вы задаете, и при этом интервалы, по возможности, были бы ровными, круглыми числами. Например, если вы сделаете 4, то у нас интервалы по 1000 ровно. Не просто так получилось, что границы у нас получились 2000, 3000, 4000, 5000, 6000.

Давайте то же самое теперь сделаем в ggplot, но с группировкой по видам:

# Гистограмма с группировкой по видам
ggplot(data = my_penguins) +
  aes(x = body_mass_g, fill = species) +
  geom_histogram()

Когда используете гистограмму, обязательно, обязательно пробуйте разное количество интервалов, потому что картинка при этом может довольно сильно меняться.

Density plot

Есть еще альтернатива гистограмме, которая при этом весьма и весьма популярна — а именно density plot. Давайте здесь пойдем по-другому: сначала нарисуем это в ggplot, а потом разберем, как это работает. Потому что по сути density — это просто такая сглаженная гистограмма. Ее так часто и объясняют.

# Density plot
ggplot(data = my_penguins) +
  aes(x = body_mass_g, fill = species) +
  geom_density()

Вот, и у нас появились вот такие вот сглаженные гистограммы. Правда, здесь они уже не накладываются, не стакуются друг на друга, а перекрывают друг друга. Это можно увидеть, если мы добавим прозрачность через alpha. Поскольку мы хотим одинаковую прозрачность для всех, то мы это прописываем вне функции aes().

# Density plot с прозрачностью
ggplot(data = my_penguins) +
  aes(x = body_mass_g, fill = species) +
  geom_density(alpha = 0.5)

Видите, они как бы пересекаются. Это, кстати, тоже можно поменять. Если мы сделаем position = "stack":

# Density plot со стакингом
ggplot(data = my_penguins) +
  aes(x = body_mass_g, fill = species) +
  geom_density(alpha = 0.5, position = "stack")

Продолжение следует...

Как именно получается сглаженная гистограмма, какие бывают ядра, а также море других визуализаций можно получить в полной версии конспекта. Иван строит violin-плоты, strip-плоты, beeswarm, raincloud и еще некоторые графики, для которых просто нет какого-то специального названия.

Забрать у Варвары ли Евы

Бластим. Статистика, R и анализ данных

2026-03-27 13:15 Статистика, R и анализ данных Полезные советы