Мастер-класс: валентинки в R

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

 r
library(ggplot2)
dat <- data.frame(t = seq(0, 2*pi, by = 0.01))
x <-  function(t) 16 * sin(t)^3
y <- function(t) 13*cos(t) - 5*cos(2*t) - 2*cos(3*t) - cos(4*t)
dat$y <- y(dat$t)
dat$x <- x(dat$t)

heart <- ggplot(dat, aes(x,y)) +
  geom_polygon(fill = "#f768a1", col = "firebrick", alpha = 0.9) +
  theme_void()+
  theme(panel.background = element_rect(fill = '#e7e1ef'))

heart + annotate("text", x = -12, y = 16, color = '#ae017e', label = "От кого: ")+
  annotate('text', x = -12, y = 13, color = '#ae017e', label = 'Кому: ') +
  annotate('text', x = -1, y = 0, color = '#49006a', label =  'Люблю тебя больше чем рисовать графики ❤️')+
  annotate('text', x=-5:5,y=15,colour="red", label = '❤️')

С помощью annotate можно добавлять любые дополнительные элементы, чтобы узнать координаты куда, советую временно закомментировать строчку theme_void.

#R #stat_fun

UpSet плоты в R для визуализации пересечений списков

Сегодня разберемся, что это такое, как их строить и в каких случаях они нужны.
UpSet плоты - это остроумная альтернатива диаграммам Венна для большого числа пересечений.
Диаграммы Венна отражают все возможные пересечения, объединения, разности нескольких списков. Если этих множеств не больше пяти, то такие диаграммы позволяют визуально оценить взаимоотношения списков. Однако, диаграммы Венна плохо масштабируются, и в ситуации, когда надо показать больше 5 списков или множеств, неприменимы.

Для визуализации большого количества пересечений Lex с соавторами предложили концепцию UpSet плота [Lex et al., 2014]. Сейчас есть R-пакеты UpSetR и ComplexUpset в репозитории CRAN, есть реализации на Python, например py-upset, UpSetPlot.

Начнем разбор с R-пакета UpSetR. Для отрисовки графиков есть вариант вбить руками количество элементов в списках и прописать все пересечения и количество элементов в них, но очевидно, что это не оптимальный вариант (называется fromExpression()). Намного функциональнее подать на вход список интересующих векторов и использовать функцию fromList(). Эта функция создает матрицу пересечений списков, которая и используется для отрисовки upset-графика. Функция upset() позволяет настроить порядок отображения списков, выбрать только определенные пересечения и выбрать, отображать ли пустые множества. Библиотека опирается только на базовый R, дополнительные библиотеки не нужны. Это как плюс, так и минус, к тому же репозиторий пакета давно не обновлялся, и поэтому сейчас можно использовать более новый пакет ComplexUpset.

Пакет ComplexUpset представляет надстройку над пакетом UpSetR с использованием фишек пакета ComplexHeatmap и ggplot2. Вот здесь разобрано, как работает пакет, возможности кастомизации графиков в нем более гибкие, чем в старом пакете, к тому же интеграция с ggplot2, на мой взгляд, скорее преимущество, чем недостаток. Например, можно сопровождать график боксплотом, подкрашивать уровни пересечения (например, пересечения двух списков обозначить зеленым, трех списков - синим и так далее), выделять цветом наиболее интересные пересечения списков. Еще этот пакет можно использовать в Python, почитать подробнее можно здесь.

Резюмируя - UpSet плоты - это отличный вариант визуализации пересечений множества групп, который часто встречается в биоинформатических статьях.
В следующем посте прикрепляю примеры графиков upset и сравнение с диаграммами Венна.

#R #data_vis

Обновление purrr до версии 1.0

purrr - пакет в R из core tidyverse для замены циклов и функций семейства *apply(). Немного подробнее я писала здесь и скидывала туториал.

Вчера вышло видео Хэдли Викхама про обновление до 1.0 с объяснениями, что изменилось. Однако мне было удобнее прочитать документацию, и оказалось, что обновление вышло в конце прошлого года. Я обновила purrr до версии 1.0.1 и протестировала основные моменты, что изменилось. Давайте разберем.

- Функции для конвертации списков в датафрейм map_dfr()/map_dfc() и их аналоги map2_(), pmap() стали deprecated, то есть разработчики перестали поддерживать эти функции. Объясняется это тем, что они не являются настоящими map-функциями, поскольку в них не сохраняется один к одному соответствие инпута и аутпута. Вместо этого рекомендуется использовать сначала просто map(), а затем эксплицитно прописывать функцию list_rbind()/list_cbind() для соединения датафрейма по строкам/столбцам как в map_dfr()/map_dfc(), соответственно. Вот пример из документации:

paths |> map_dfr(read_csv, .id = "path")
# now
paths |> 
  map(read_csv) |> 
  list_rbind(names_to = "path")

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

- Добавили progress bar для отслеживания времени выполнения функций. Реализуется это подачей аргумента .progress = TRUE в любую map функцию. Однако, если задача выполняется быстро, то прогресс бар не будет появляться (я проверила). Также можно прописать, что именно выполняется, подав строку в .progress, например
map(1:100, \(x) Sys.sleep(0.1), .progress = 'Saving plot')
Это мне очень понравилось, думаю, что буду использовать для долгих маппингов. Всегда приятно знать, что функция не упала, а выполняется и примерно понимать сколько осталось.

- Добавили функцию map_vec() для работы с любым типом векторов. map_vec() является расширением map_lgl(), map_int(), map_dbl() и map_chr() для произвольных типов векторов - например для факторов, даты и времени.

- Уходят от синтаксиса формулы для анонимных функций (~ .x + 1) к другому синтаксису: \(x) x + 1, и от магриттаровского пайпа %>% уходят к нативному |> (хотя это уже давно). Последнее аргументируется тем, что нативный пайп будет работать без подключения дополнительных пакетов, это важно например при разработке своих пакетов, чтобы не получать лишние зависимости. Синтаксис новых анонимных функций я вообще в первый раз вижу, но могу согласиться, что он выглядит более понятно, чем синтаксис формулы. Хотя я только привыкла к формуле, как снова надо переучиваться. Впрочем, старое доброе function(x) x + 1 в качестве анонимных функций тоже никто не отменял, и на мой взгляд это самый понятный синтаксис.

Пример работы нового синтаксиса анонимных функций:

# Previously we wrote
1:10 %>%
  map(~ rnorm(10, .x)) %>%
  map_dbl(mean)
#>  [1]  0.5586355  1.8213041  2.8764412  4.1521664  5.1160393  6.1271905
#>  [7]  6.9109806  8.2808301  9.2373940 10.6269104

# Now we recommend
1:10 |>
  map(\(mu) rnorm(10, mu)) |>
  map_dbl(mean) 
#>  [1]  0.4638639  2.0966712  3.4441928  3.7806185  5.3373228  6.1854820
#>  [7]  6.5873300  8.3116138  9.4824697 10.4590034

Я разобрала наиболее интересные для меня изменения, это не все, поэтому рекомендую далее прочитать документацию самостоятельно.

#R #tidyverse #purrr

Делаем t-тесты или u-тесты в R сразу для многих колонок с помощью `tidyverse`-подхода

Бывают ситуации, когда измерили много количественных переменных для двух групп, например контрольной и с заболеванием, я встречала такие данные у медиков или как результаты масс-спектрометрии. Планируется сравнить все эти количественные переменные тестом Стьюдента или Манна-Уитни, но вручную прописывать 3 или больше раз t-тест кажется не очень хорошей идеей. Что можно сделать? Будем использовать следующий подход: сначала развернем таблицу в длинный формат, соберем в списки значения по каждой группе и количественной переменной, затем таблицу снова превратим в широкий формат, но уже в виде
переменная1 [список значений контрольной группы] [список значений экспериментальной группы]

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

 r
library(dplyr)
library(tidyr)
set.seed(2)
df <- data.frame(lapply(rep(100,15), rnorm), 
                 group = rep(c('control', 'treatment'), each = 50)) # генерируем данные
colnames(df)[1:15] <- paste0('marker', 1:15) # меняем имена колонок на более понятные
df %>% 
  select(where(is.numeric), group) %>% # это на случай, если в исходном датафрейме не только числовые переменные
  pivot_longer(cols = -group, names_to = 'variable') %>% # преобразуем датафрейм в long-формат
  group_by(group, variable) %>% # группируем по типу обработки и типу переменных
  summarise(value = list(value)) %>% # собираем в списки
  pivot_wider(id_cols = c(variable), names_from = group) %>% # разворачиваем обратно
  group_by(variable) %>% # группируем для проведения стат теста
  # запускаем тест Манна-Уитни, сохраняем u-значение и p-value
  mutate(p_value = wilcox.test(unlist(control), unlist(treatment))$p.value,
         u_value = wilcox.test(unlist(control), unlist(treatment))$statistic)
#> `summarise()` has grouped output by 'group'. You can override using the
#> `.groups` argument.
#> # A tibble: 15 × 5
#> # Groups:   variable [15]
#>    variable control    treatment  p_value u_value
#>    <chr>    <list>     <list>       <dbl>   <dbl>
#>  1 marker1  <dbl [50]> <dbl [50]>  0.293     1403
#>  2 marker10 <dbl [50]> <dbl [50]>  0.0403    1548
#>  3 marker11 <dbl [50]> <dbl [50]>  0.269     1411
#>  4 marker12 <dbl [50]> <dbl [50]>  0.997     1249
#>  5 marker13 <dbl [50]> <dbl [50]>  0.323     1106
#>  6 marker14 <dbl [50]> <dbl [50]>  0.560     1335
#>  7 marker15 <dbl [50]> <dbl [50]>  0.667     1313
#>  8 marker2  <dbl [50]> <dbl [50]>  0.117     1478
#>  9 marker3  <dbl [50]> <dbl [50]>  0.931     1263
#> 10 marker4  <dbl [50]> <dbl [50]>  0.866     1225
#> 11 marker5  <dbl [50]> <dbl [50]>  0.791     1211
#> 12 marker6  <dbl [50]> <dbl [50]>  0.986     1247
#> 13 marker7  <dbl [50]> <dbl [50]>  0.920     1235
#> 14 marker8  <dbl [50]> <dbl [50]>  0.0169    1597
#> 15 marker9  <dbl [50]> <dbl [50]>  0.707     1195

Если понадобится сделать не тест Манна-Уитни, как в примере, а t-test, то надо просто поменять в последней команде wilcox.test() на t.test().

Обращайте внимание на проблему множественных сравнений! Даже в моем примере обнаружились статистические значимые различия там, где их не должно быть, в двух строках (2 и 14). Различий быть не должно, поскольку и "контрольная" и "экспериментальная" группы генерируются из одного распределения (стандартного нормального со средним 0 и стандартным отклонением 1). Следовательно, мы совершили ошибку первого рода, когда нашли различия там, где их нет.
Я сталкивалась с исследованиями у медиков, где они собирали как можно больше анализов, измеряли как можно больше антропометрических данных (рост, вес, толщина плеча, отношение толщины плеча к толщину предплечья), чтобы все-таки найти значимые отличия хоть где-нибудь. И ведь находили! Но на самом деле вопросы даже не к самим студентам или ординаторам, а к их научным руководителям, которые не ставили четкую гипотезу и дизайн эксперимента, и к комиссии, где без заветного p-value < 0.05 можно и не получить диплом.

Мораль такая, что измерить больше параметров не значит лучше, нужно формулировать гипотезу ДО начала эксперимента, чтобы не получилось слишком много сравнений, поскольку после применения поправки на множественное тестирование может пропасть вся значимость.
Поэтому уметь сделать много тестов Стьюдента за раз - это хорошо, а еще лучше понимать, нужно ли это делать, и может быть сделать есть смысл сделать только несколько самых важных тестов.

Так что рекомендую всем подписчикам, формулируйте гипотезы заранее, остерегайтесь проблем с множественным тестированием и помните, что p-value > 0.05 - тоже результат!

#R #stats #tidyverse

Ошибки первого и второго рода — мнемоническое правило

На сайте machinelearning.cards нашла любопытную картинку, как запомнить, что является ошибкой первого и второго рода. Несколько лет назад я путала, где какая ошибка, и приходилось каждый раз гуглить и проверять.

Ошибка I рода — false positive, то есть неверное отклонение нулевой гипотезы (например, нашли значимые различия там, где их на самом деле нет).
Ошибка II рода — false negative, не отклонили нулевую гипотезу, хотя должны (нулевая гипотеза неверна).

Эта картинка нравится графическим отражением ошибок I и II рода в буквах P и N, и на мой взгляд, сделано довольно изящно.
В комментарии скину другие способы запомнить, где какая ошибка.

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

Пишите, не путаете ли вы ошибки I и II рода, как удалось запомнить и какой из предложенных способов понравился больше.

#stats

Как неправильно хранить данные в таблицах?

Недели три назад в интернете форсился мем про хранение "параметров" людей в Excel-таблицах (примеры в комментариях). Я не поддерживаю такой утилитарный подход к окружающим, но оставляя за скобками морально-этические аспекты, хотелось бы рассказать, почему с технической стороны это было организовано неверно.
Основная претензия в том, что переменные (параметры) записаны в строках, а наблюдения (в данном случае люди) в столбцах. Это влечет за собой массу неудобств.

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

Во-вторых, для анализа подобных таблиц средствами языков программирования критично иметь одинаковый тип данных в столбцах. К примеру, в R это вытекает из структуры датафрейма (поскольку столбцы – это вектора, в которых требуется одинаковый тип данных).
Для переменных (признаков, параметров) обычно гарантируется одинаковый тип данных, в то время как для наблюдений – вовсе нет. Следовательно, в R при загрузке такой таблицы, все столбцы приведутся к наиболее общему типу – скорее всего строковому, что также усложнит или сделает невозможным анализ.

В-третьих, не в контексте конкретно таких таблиц, можно отметить, что в Excel-таблицах формата .xlsx максимальное количество столбцов: 16384 (2^14), в то время как строк: 1048576 (2^20). Следовательно, при заполнении наблюдений по столбцам, место закончится гораздо раньше.
В качестве иллюстрации я хотела привести пример, как в 2020 году британские медики заполняли новые случаи Covid-19 в эксель таблицу по столбцам, в результате чего новые случаи ковид просто не влезли! Однако я почитала еще, и оказалось, что проблема была немного другая: данные хранили в строках (ура!), но вместо .xlsx использовали устаревший формат .xls, который поддерживает только 65536 (2^16) строк. Похоже, что медицинские центры присылали данные в формате .csv, где нет ограничения на количество строк, в то время как министерство здравоохранения (Public Health England) открывало их в .xls формате, и все строки после 65536 для них исчезали. Таким образом потерялось около 16 тысяч положительных результатов ковид-тестов (а сколько всего было потерянных строк, не сообщается).

Вообще хранение большого количества строк в Excel в целом вызывает сомнение, но это отдельная тема, пока что отмечу, как у генетиков некоторые названия генов (MARCH1, SEPT2) превращались в даты (например описано здесь).

UPD: в комментариях дополнили про недостатки экселя:
- Превращение чисел в даты, что сильно раздражает и с этим сталкиваются не только генетики
- Привязка к региональным разделителям (по умолчанию запятая как десятичный разделитель, что приходится менять)
- Отсутствие стандартного обозначения пропущенных значений

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

#stat_humor #data #excel

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

Я давно составляла и размещала на канале рекомендуемые источники по статистике, однако кое-что устарело, плюс нашла еще классные материалы, так что решила обновить.

Книги
1) Статистика и котики из старой подборки пусть остаются, книга критикуется где-то за излишние упрощения, но для старта все еще неплохо.

2) Медико-биологическая статистика Гланца пусть тоже остается, хотя я сама недолюбливаю эту книгу, по крайней мере перевод на русский язык. В нем много небольших неточностей, например достоверность вместо статистической значимости, плюс переведено только старое издание, в котором нет упоминаний поправок на множественное тестирование FDR. Однако во многих чатах его советуют, поэтому оставляю в подборке.

3) Статистика для всех (Statistics in a nutshell) Сары Бослаф неплохая книга, однако в переводе большинство терминов насильно переведены на русский (иногда довольно неудачно), и плохо, что не дано исходного термина на английском, поскольку гуглить придется все равно так. Читать рекомендую самое свежее издание, благо найти в интернете легко.

4) Анализ данных и статистика в R Ивана Позднякова, что хорошо, автор постоянно обновляет часть по R (особенно по tidyverse) и дополняет часть по статистике, например многомерные методы. Материалы вплоть до обобщенной линейной модели доступны, написано хорошо, в общем рекомендую.

5) Рандомизация и бутстреп: статистический анализ в биологии и экологии с использованием R. В.К. Шитиков, Г.С. Розенберг. Остаются с предыдущей подборки.

Курсы

У меня в старой подборке есть однозначная рекомендация курсов Анатолия Карпова по статистике, однако после длительных дискуссий на тему качества этих курсов, я пересмотрела первую часть (Основы статистики).
Соглашусь, что действительно много неточностей как мелких, так и иногда довольно существенных в этом курсе. Следовательно, сейчас я не могу рекомендовать его, поскольку вышло много бесплатных и более качественных материалов (однако для 2013го года это был хороший продукт, так как русскоязычных курсов по статистике практически не было). Возможно, я сделаю более подробный разбор неточностей, которые я собрала, но пока не классифицировала и не оформила в пост.

Следовательно, рекомендую:
1) Курсы Марины Варфоломеевой: Линейные модели, дисперсионный и регрессионный анализ с использованием R, с еще несколькими курсами можно ознакомиться здесь. Наткнулась на материалы практически случайно, увидела что в других чатах тоже рекомендуют, сама кое-какие идеи оттуда брала для своих лекций (конечно с ссылкой на источник). Очень качественно и подробно разобрано практически все, с чем можно столкнуться в статистике, супер круто, что такие материалы есть в открытом доступе.

2) Курс Data Analysis with R Specialization на курсере, я сама не проходила, но много где советовали, программа очень солидная.

3) Платные курсы от Института биоинформатики и бластима. Ссылки кидать не буду, но они легко гуглятся, если кто-то не найдет, пишите в личку или в комменты. На бластиме я преподавала, как многие тут уже знают, могу ручаться за качество лекций и поддержки ассистентов. Курсы или точнее программы повышения квалификации от института биоинформатики тоже однозначно рекомендую, много общалась с людьми оттуда, делают очень качественный продукт.

Ютуб-каналы

1) StatQuest - отличный канал с короткими, но очень четкими разборами конкретных тем по статистике. У автора забавный стиль изложения, с приколами и всякими фразочками, сам он биоинформатик (а по первому образованию музыкант).
2) TileStats - пожалуй еще лучше чем первый канал (хоть и малоизвестный). Тоже формат коротких видео, где наглядно показывается применение формул, как это выглядит на графике и что означает. Абсолютно незаменим для подготовки собственных лекций, очень удачные идеи по визуализации сложных концептов. Стиль изложения более строгий, без забавных фразочек, но все равно максимально понятно.

Если что-то забыла из хороших материалов, пишите комментарии, буду стараться держать подборку лучших книг и курсов всегда обновленной!

#recommendation

Пятничный мем.
Еще вариант: tidyverse/data.table

Идея мема взята у Хэдли Викхама

#R #tidyverse #stat_fun

ПоздRавление от подписчика
У админа канала вчера был день рождения! А еще завтра (а у кого-то уже сегодня) день знаний, поэтому мы поздравляем админа и вас таким тортиком, сделанным в R, а заодно ещё раз исследуем базовую графику!
Вы можете тоже сделать такой тортик для себя и друзей:

install.packages('DescTools')
library('DescTools')
candle = function(pos)
{
x=pos[1]
y=pos[2]
rect(x,y,x+.2,y+2,col="red")
polygon(c(x+.05,x-.1,x+.1,x+.3,x+.15,x+0.05), c(y+2,y+2.3,y+2.6,y+2.3,y+2,y+2),col="#EECF08")

}

plot(c(0,10), c(0,10),type="n", bty="n",xaxt="n",yaxt="n", main="Cake", xlab="",ylab="")
DrawEllipse(x = 5, y = 2,col="seagreen",radius.x=4.4,radius.y=1.7)
DrawEllipse(x = 5,y = 2,col="violet",radius.x=4,radius.y=1.4)
rect(1,2,9,5,col="violet",border="violet")
lines(c(1,1),c(2,5))
lines(c(9,9),c(2,5))
DrawEllipse(x = 5,y = 5,col="brown",radius.x=4,radius.y=1.4)

candle(c(2.5,4.5))
candle(c(3,5))
candle(c(4,4.5))
candle(c(5,5))
candle(c(6,4.5))
candle(c(7,5.2))

Источник кода: https://qiuyiwu.github.io/2019/01/29/Birthday/ (там можно увидеть, что некоторые функции в аре устаревают и приходится искать аналоги)

В комментариях получившийся торт)

Неявное приведение типов - иллюстрация

Пятничный мем! Сделала еще давно, но вспомнила только сейчас.

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

#stat_fun #R

Опубликовала первое видео на ютубе, посвященное истории языка R, и ответу на животрепещущий вопрос*, почему в сообществе принято использовать оператор <-, а не =. Планируется записать ~15 видео по R и примерно столько же по статистике, так что подписывайтесь, все будет опубликовано здесь.

По результатам опроса 66% подписчиков канала используют оператор <-, что в целом неудивительно, скорее меня озадачило использование обоих операторов одновременно (почему - смотрите в видео).

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

*примерно раз в месяц в биоинформатических чатах появляется обсуждение, чем же <- лучше =

#R #youtube