Квартет Энскомба (Anscombe's) или важность визуализации данных перед началом анализа

Подготовила небольшую статью о таком любопытном наборе данных, который показывает что нужно опираться не только на средние-стандартные отклонения, но и точно знать, что происходит в данных, опираясь на визуализацию.
https://rpubs.com/lena_astr/958484

Бонусом список ссылок на книги и статьи по качественному представлению данных.
Попробовала сделать свой первый документ на quarto, пока полет нормальный, но почему-то на github pages слетело форматирование, поэтому пока залила на rpubs.
Конечно, не всегда можно так легко отразить на одном или нескольких графиках, что происходит в данных. Так что о способах визуализировать многомерные данные, например данные экспрессии генов, поговорим в следующий раз.

#R #base_stat

Супер крутая инструкция по освоению библиотеки purrr для замены циклов

Немного контекста: в R принято обходиться не то чтоб без вложенных циклов, а без циклов вообще. Как это реализуется? Большинство функций в R векторизованы, что означает, что функция будет работать со всеми элементами вектора без необходимости проходить цикл и обрабатывать каждый элемент по одному.
Простейший пример векторизации - возвести числа от 1 до 10 в квадрат. Классическим подходом было бы написать цикл от 1 до 10 и возвести в квадрат каждый элемент. Однако в R это можно сделать гораздо проще!

(1:10)**2
 [1]   1   4   9  16  25  36  49  64  81 100

Мы просто взяли вектор чисел от 1 до 10 и возвели в квадрат одной строчкой без всяких циклов. Это самый простой прием, для более сложных пригодятся функции семейства apply или их более продвинутый аналог функции семейства map_*, map2_* из библиотеки purrr. Именно последним посвящен туториал, скачать который можно по ссылке.

Пример из жизни применения map функции. У меня много однотипных датасетов в определенной директории с одним расширением (например csv). Задача считать все датасеты в один лист для последующих манипуляций с ними. Чтобы не писать цикл, используем purrr::map

library(tidyverse)
csv_file_names <- dir(pattern = 'csv')
list_csv <- map(csv_file_names, ~read_csv(.x))

ИЛИ

list_csv <- map(csv_file_names, function(x) read_csv(x))

, если без синтаксического сахара, а использовать более привычный apply-способ

Вообще на тему функционального программирования в R, векторизации и сравнения *apply и map_* функций можно довольно много что расписать, если есть интерес к этой теме, ставьте реакцию с китом (потому что с котом еще нет реакций, а символ purrr - котик).

#R #tidyverse #purrr

Пределы погрешностей: что это, зачем использовать и как интерпретировать?

Продолжаем тему базовой статистики (кстати базовой не значит простой). Рада представить разбор видов error bar, их различий и правильной интерпретации:

https://ubogoeva.github.io/R4Analytics/posts/types_of_error_bars.html

Тема оказалась сложнее, чем я рассчитывала, например, я сама не знала что пределы погрешностей существуют двух принципиально разных типов. Но думаю получилось интересно, по крайней мере мне точно было интересно писать. В этот раз сверстала в Rmd файл, с quarto пока не разобралась, почему на github pages размещается некрасиво.
Бонусом альтернативные способы визуализации данных. Прикрепляю еще небольшой тизер к посту, очень нравится получившаяся картинка.

Пишите комментарии, какие виды error bar используете для представления своих данных, сталкивались ли с их неправильной интерпретацией и понравились ли предложенные способы визуализации?

#base_stat #R #ggplot2 #data_vis

Как задать хороший вопрос в R чате?

Для русскоязычного R-комьюнити есть два больших чата: R (язык программирования) и Горячая линия R. Вообще правила прописаны в закрепленных постах чатов, но я бы хотела поделиться своим опытом задавания вопросов и создания воспроизводимых примеров (reprex, репрекс).
Reprex (reproducible example) - это минимальный пример кода, воспроизводящий ошибку или описывающий, что требуется сделать.
Такой пример кода удобно копируется в чат с соответствующим форматированием, его легко читать и можно сразу же скопировать себе, чтобы попробовать помочь.

У меня не поместилось все в лимит символов в телеграме, поэтому как всегда ссылка на гитхаб:
https://ubogoeva.github.io/R_question_how_to_ask.html

Надеюсь, будет полезно и повысит культуру общения в чатах)
#R #help

Рисуем дорожные знаки в R

Наткнулась на пост в линкедине, где автор (профессор Nuno Sepúlveda) интересовался, из какого распределения пришел дорожный знак лежачий полицейский. Мне показалось забавным, и я воспроизвела код для этого знака, а еще попробовала нарисовать знак "Неровная дорога" в R, опираясь на похожую идею.

https://ubogoeva.github.io/road_signs_r.html

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

#R #stat_fun

Подборка и небольшой обзор нейросетей для написания кода на R

- Начнем с всем известной chatGPT3. Вообще история общения с этой нейросетью получилась интересная, потому что сначала я допытывала определение p-value, на что сеть отвечала набором неправильных интерпретаций почти как из справочника "чем p-value не является". После нескольких попыток и формулировок "write a correct p-value definition" сеть обучилась и стала выдавать-таки правильную формулировку, хотя по-моему студенты обучались быстрее. Теперь сеть уже не ошибается в этом вопросе, потому что она локально обучилась на моих вопросах и реакции на ответы. При этом для всех остальных юзеров это не так.

С R кодом ненамного лучше, часть кода еще похожа на правду, но где-то используются устаревшие подходы, а где-то и вообще несуществующие названия функций. Нужно стараться как можно более полно и конкретно описать, что нужно сделать, тогда повышается шанс корректного ответа. Самое плохое, что сеть не использует правильные паттерны, для нее норма писать один и тот же код 10 раз (например mean1 = mean(variable1), mean2 = mean(variable2) и так далее). Поэтому неудивительно, что нейросеть забанили на stackoverflow. Общаться с нейросетью весело, но, на мой взгляд, не стоит ожидать, что она станет заменой программистов в ближайшее время.

- Rtutor - shiny app, основанная на нейросети OpenAI's Davinci, схожей с chatGPT. Заточена специально на написание R кода. Можно загружать свои данные и описывать текстом, что нужно сделать, далее нейросеть подсказывает решения. Я сильно много не экспериментировала с этой сетью, потому что каждый запрос платный для автора сервиса. Если планируется активно использовать сеть, то предлагается создать свой аккаунт. Сам разработчик говорит, что результат работы нейросети нужно использовать лишь как стартовую точку, полученный код нужно проверять и редактировать.

- Интеграция chatGPT в RStudio. Дополнение в арстудию, которое позволяет сразу из IDE отправлять запросы на сервер chatGPT3. Я сама не тестила, но звучит интересно.

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

#R #neural_network #stat

Подборка источников о ggplot2 и обзор некоторых интересных расширений

Начнем с того, что такое ggplot2.
ggplot2 - пакет в R для визуализации данных, входит в ядро tidyverse. На самом деле это не просто "какой-то очередной пакет для графиков", а целая философия визуализации! gg в названии означает grammar of graphics - язык для построения графиков, описанный в книге Леланда Вилкинсона (Leland Wilkinson).

Делюсь материалами по освоению ggplot2:

- Англоязычные
• Замечательная инструкция как работает пакет, изучение с самых основ. Далее разбирается визуализация конкретных данных, с самого начала, то есть не сразу готовое красивое решение, а постепенно, с улучшением и дополнением на каждом этапе (а иногда и нужно вернуться на предыдущий шаг);
• Интересный разбор использования stat_summary для создания графиков с error bar (пределами погрешностей);
• Очень объемный туториал, разбирается много важных вещей (сама целиком не прочитала).

- Русскоязычные
• глава книги Ивана Позднякова про устройство пакета
• воркшоп по ggplot2 Татьяны Балтыжаковой (я все не читала, но судя по заголовкам очень хороший материал)
• вебинар Филиппа Управителева, есть также запись на ютуб

Материалов много, поэтому я решила пока не писать свой интродакшн в ggplot, а вместо этого рассмотреть некоторые интересные расширения. Для пакета существует очень много расширений (на момент написания поста 120), с полным списком можно ознакомиться здесь.

• Расширение patchwork для объединения нескольких графиков в один. Пример работы можно посмотреть у меня в посте про пределы погрешностей. Очень крутая вещь, позволяет в том числе сразу добавлять аннотации на плоты (например А, B, C), что делает их отлично отформатированными для научной публикации.
• По поводу публикаций - здесь есть несколько расширений, чтобы например добавлять результаты статистических тестов на графики: ggstatplot, ggpubr, ggsignif. Минус этих расширений в том, что они поддерживают заданный набор базовых плотов: барплоты, боксплоты и другое, но не все что угодно, как обычный ggplot2 без расширений. Кроме того, с коллегой обнаружили, что ggpubr не поддерживает односторонние тесты, конечно можно залезть в код и отредактировать, но немного неприятно.
• Еще по поводу статей: пакет hrbrthemes делает очень красивые графики для публикаций с нужными шрифтами, фоном, выравниванием заголовков.
• Для генетиков и биоинформатиков: gggenes позволяет визуализировать гены и их относительное расположение, направление цепи, думаю экзон-интронную структуру также можно воспроизвести.
• ggpattern, чтобы раскрашивать графики котятами (и не только).
• ggheatmap для отрисовки хитмапов

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

Также рекомендуйте в комментариях, какие есть еще интересные расширения.

#R #data_vis #ggplot2

Мастер-класс: валентинки в 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