Язык R

и его применение в биоинформатике

Лекция 9

Анастасия Жарикова

27 октября 2023

Что вас ждет

  • Вы уже знаете статистику

  • При защите проекта мы можем спрашивать теорию, нулевые гипотезы тестов, условия применения тестов, …

  • В качестве дополнительных вопросов могут быть любые вопросы по материалу курса статистики

Будни биоинформатика

  • 70% времени - собрать данные в приемлемый для анализа вид

  • 5% времени - применить готовые функции

  • 20% времени - придумать и реализовать наглядную визуализацию

  • 5% времени - подумать о биологии

Зачем?

У вас есть данные

Их можно визуализировать, что-то посчитать

Зачем статистика? И что это такое?

Зачем нужна статистика?

Инструмент познания окружающего мира

Математическая наука о сборе, анализе, интерпретации и представлении данных

Так зачем?

Как исследовать?

  • Эксперимент: предполагает воздействие

  • Наблюдение: сбор данных без вмешательства и влияния

Пример вопросов

Помогает ли диета?

Новое лекарство лучше старого?

Влияют ли прогулки на продолжительность жизни?

Какой сервис доставки лучше?

Кошек или собак предпочитают в качестве домашних животных в Москве?

Какой средний вес у мопсов?

Как ответить?

Кто больше времени уделяет спорту:

женщины или мужчины?

Как ответить?

Способ 1: выяснить у всех женщин и мужчин, насколько активный образ жизни они ведут

Как ответить?

Способ 2: набрать репрезентативную выборку и выяснить их спортивную активность

  • Сколько людей набрать?

  • Кого брать?

  • Как вообще выбирать людей?

  • Как оценивать спортивную активность?

Генеральная совокупность и выборка

Хотим получить знание о генеральной совокупности, т.е. о всех возможных участниках, которых касается наш опрос или эксперимент

Исследовать генеральную совокупность долго, дорого и часто просто невозможно

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

Генеральная совокупность

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

От этого зависит стратегия составления репрезентативной выборки

Способы создания репрезентативной выборки

  • Простая случайная выборка

  • Кластерная случайная выборка

  • Стратифицированная случайная выборка

  • Систематическая выборка

  • …

Простая случайная выборка

Каждый элемент генеральной совокупности имеет равный шанс попасть в выборку

Кластерная случайная выборка

  • Разбиваем генеральную совокупность на кластеры

  • Случайным образом отбираем несколько кластеров

  • Включаем ВСЕХ участников этого кластера

  • Больше подходит для однородной генеральной совокупности

Стратифицированная случайная выборка

  • Разбиваем генеральную совокупность на кластеры

  • Случайным образом из каждого кластера отбираем случайных участников

  • Больше подходит для неоднородной генеральной совокупности

Систематическая случайная выборка

  • Упорядочиваем всех членов генеральной совокупности в некотором порядке (по алфавиту, по возрасту, …)

  • Выбираем случайным образом начальную точку

  • Отбираем каждого n-ного участника

Невероятностные способы набора выборок

  • “Удобная” выборка

Стоим у дверей факультета и опрашиваем прохожих

  • Метод добровольного опроса

Просим подписчиков заполнить анкету

  • Снежный ком

Отбираем несколько человек для начала исследования и просим их поспособствовать дальнейшему набору

  • Целевой отбор

Решаем, кто нам нужен, и набираем таких людей

Систематическая ошибка отбора

  • Исследовать только интервал, укладывающийся в теорию

  • Ошибка выжившего

  • Некорректная “удобная” выборка

  • Ошибка меткого стрелка

Систематическая ошибка неответа

Участники, попавшие в выборку, значительно отличаются от непопавших.

В итоге выборка не репрезентативна.

Пример: В кофейне решили выяснить, кто больше предпочитает новый сорт кофе: “совы” или “жаворонки”. Опрос устроили в 8 утра.

И еще много других ошибок…

Переменные

Влияет ли занятия спортом на настроение?

  • Независимая переменная: наличие занятий спортом (количество раз в неделю)

  • Зависимая переменная: оценка настроения (по какой-то шкале)

  • Конфаундеры: влияют на значения зависимой переменной, но не контролируются в эксперименте (время года, возраст и пол участников, масса тела, сфера деятельности, ….)

Конфаундеры

При планировании исследования нужно тщательно продумать, что может влиять на результат

По возможности устранить или снизить эти влияния

Ослепление

Проводят исследование эффективности лекарства. Набирают 2 группы: участникам одной группы дают лекарство, другой - пустышку

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

  • Слепое исследование: испытуемый не знает, в какой он группе

  • Двойное слепое исследование: ни испытуемый, ни экспериментатор (врач) не знают, кто в какой группе

  • Тройное слепое исследование: даже тот, кто обрабатывает данные, не знает, какая группа получала плацебо, а какая лекарство

Типы данных

Визуализация категориальных данных

library(tidyverse)
iris %>%
  sample_n(60) %>%
  ggplot(aes(x = Species, 
             fill = Species)) +
  geom_bar() + theme_bw() +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold"))

Для порядковых категориальных данных важен порядок столбцов

Визуализация непрерывных данных

iris %>%
  ggplot(aes(x = Sepal.Length)) + 
  geom_histogram(fill = 'orange', color = 'brown') + theme_bw() +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold"))

Гистограмма

Многое зависит от ширины бина

iris %>%
  ggplot(aes(x = Sepal.Length)) + 
  geom_histogram(fill = 'orange', color = 'brown', binwidth = 1) + theme_bw() +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold"))

Статистика

  • Описательная статистика

  • Индуктивная статистика

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

Любой набор данных нужно сначала изучить

  • список переменных

  • количество наблюдений

  • наличие пропущенных значений

  • распределения численных переменных

  • композиция категориальных переменных

  • наличие выбросов

  • визуализация

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

Функция от выборки

Характеризуем выборку одним числом

  • медиана

  • среднее

  • мода

  • дисперсия

  • разброс

  • …

Выборочное среднее

\[\overline{X} = {x_{1} + x_{2} + ... + x_{n} \over n}\]

set.seed(123)
x <- rpois(6, 3)
x <- sort(x)
x
[1] 0 2 2 4 5 6
mean(x)
[1] 3.166667

Медиана

В упорядоченном наборе данных такое число, что половина из элементов набора не меньше этого числа, а половина не больше

x = c(1, 4, 8, 9, 11)
x
[1]  1  4  8  9 11
median(x)
[1] 8
x = c(1, 4, 8, 9, 11, 13)
x
[1]  1  4  8  9 11 13
median(x)
[1] 8.5

Процентили

N-ный перцентиль - это такое число, что N% элементов выборки не больше N

a <- rpois(30, 2)
quantile(a)
  0%  25%  50%  75% 100% 
0.00 2.00 2.00 3.75 6.00

  • 25-ый перцентиль - 1ый квартиль

  • 75-ый перцентиль - 3ий квартиль

  • 50-ый перцентиль - медиана

Для описания большинства данных можно использовать 95-ый перцентиль

Мода

Наиболее частое значение

set.seed(123)
a <- rpois(20, 2)
sort(a)
 [1] 0 0 0 1 1 1 2 2 2 2 2 2 3 3 4 4 4 4 5 5

Мода

library(tidyverse)
ggplot(mapping = aes(x = a)) + geom_dotplot() + theme_bw()

Мод может быть несколько

Разброс

Минимум и максимум

sort(a)
 [1] 0 0 0 1 1 1 2 2 2 2 2 2 3 3 4 4 4 4 5 5
range(a)
[1] 0 5

Боксплот - ящик с усами

Представление числовых данных через квартили

Боксплот

quantile(iris$Sepal.Length)
  0%  25%  50%  75% 100% 
 4.3  5.1  5.8  6.4  7.9 
ggplot(iris) +
  geom_boxplot(aes(Sepal.Length)) + theme_bw()

Для одной выборки не очень подходит

Боксплот

ggplot(iris) +
  geom_boxplot(aes(Sepal.Length, fill = Species)) + theme_bw()

Интерквартильный размах

ggplot(iris) +
  geom_boxplot(aes(Sepal.Length, fill = Species)) + theme_bw()

summary(iris[iris$Species == 'setosa',]$Sepal.Length)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.300   4.800   5.000   5.006   5.200   5.800 
IQR(iris[iris$Species == 'setosa',]$Sepal.Length)
[1] 0.4

Выбросы

Есть много методов поиска выбросов

Выброс должен быть объясним с точки зрения структуры эксперимента или исследования

Выброс - это ошибка или интересное наблюдение?

При поиске и удалении выбросов важно вовремя остановиться

Данные про зубы

library(ggpubr)
data("ToothGrowth")
head(ToothGrowth)
   len supp dose
1  4.2   VC  0.5
2 11.5   VC  0.5
3  7.3   VC  0.5
4  5.8   VC  0.5
5  6.4   VC  0.5
6 10.0   VC  0.5

Визуализируем среднее

ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd",
            error.plot = "errorbar",            
            add = "mean"                        
            )

Визуализируем среднее

ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd", color = "black",
            add = "jitter", add.params = list(color = "darkgray")
            )

Визуализируем среднее

ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd", 
            color = "dose", palette = "jco")

Визуализируем среднее

ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd", 
            color = "supp", palette = "jco",
            position = position_dodge(0.3))

Визуализируем среднее

ggbarplot(ToothGrowth, x = "dose", y = "len", 
       add = c("mean_se", "jitter"))

Визуализируем среднее

ggbarplot(ToothGrowth, x = "dose", y = "len", 
          add = c("mean_se", "jitter"),
          color = "supp", palette = "jco",
          position = position_dodge(0.8))

Barplot

Обязательно указывайте, что именно отображает ваш барплот:

  • количество наблюдений в группе

  • количество каких-то конкретных значений в группе

  • среднее значение параметра в группе

И обязательно указывайте, что именно отображают “усы”

Еще раз про NA

Есть ли хотя бы одно пропущенное значение?

anyNA(airquality)
[1] TRUE

Сколько всего NA?

colSums(is.na(airquality))
  Ozone Solar.R    Wind    Temp   Month     Day 
     37       7       0       0       0       0

В каких столбцах и в каких колонках находятся пропущенные значения?

Как вы это узнали?

А если столбцов 1000?

Если вы что-то оцениваете методом глаза, то вы делаете это неправильно

Случайная величина

Величина, принимающая в ходе эксперимента то или иное значение

Численные исходы некого случайного эксперимента

Точно сказать, какое это будет значение, мы не можем

Случайные величины

  • Число, выпавшее на шестигранном кубике

  • Пол ребенка

  • Время прибытия автобуса

  • Время работы лампочки

  • Количество клиентов в день

  • Число бракованных товаров в партии

Бывают дискретными и непрырывными

Дискретные случайные величины

Принимает не более, чем счетное множество значений

Распределение можно задавать при помощи таблицы: каждому значению сопоставлена его вероятность.

В сумме вероятности равны единице

Задача

Есть 2 случайные величины

Насколько согласовано их изменение?

Ковариация

Мера совместной изменчивости двух случайных величин

Большие значения одной переменной в основном соответствуют большим значениям в другой переменной, аналогично для меньших значений

Ковариация

library(tidyverse)
x = rnorm(10)
y = rnorm(10)
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cov(x, y)
[1] -0.5481997

Ковариация

x = 1:10
y = 4:13
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cov(x, y)
[1] 9.166667

Ковариация

x = 1:10
y = 13:4
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cov(x, y)
[1] -9.166667

Ковариация

Не нормированная величина

Зависит от величин переменных

Сложно интерпретировать

Ковариация

x = 1:10
y = x*0.1
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cov(x, y)
[1] 0.9166667

Ковариация

x = 1:10
y = x*10
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cov(x, y)
[1] 91.66667

Корреляция

Можно отнормировать

Получим корреляцию Пирсона

Безразмерная величина

Изменяется от -1 до 1

Можно применять не только к нормально распределенным данным

Корреляция = 0

Из этого не следует независимость переменных

Корреляция != 0

Нельзя рассуждать о причинно-следственной связи

Корреляция 0,66%

ref

Корреляция != 0

Корреляция 0,99%

Визуализируй!

set.seed(123)
x = rnorm(1000)
y = rnorm(1000)
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cor(x, y)
[1] 0.08647944

Визуализируй!

set.seed(123)
x = c(rnorm(100), rpois(5, 5))
y = c(rnorm(100), rpois(5, 6))
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cor(x, y)
[1] 0.5275666

Визуализируй!

set.seed(123)
x = c(rnorm(100), rnorm(100, 6))
y = c(rnorm(100), rnorm(100, 6))
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cor(x, y)
[1] 0.8977794

Корреляция

Если зависимость не линейная, но монотонная?

Корреляция Спирмена

Считаем корреляцию между рангами наблюдений

Корреляционная матрица

mt_filt <- mtcars[, c(1,3,4,5,6,7)]

cor(mt_filt$mpg, mt_filt$disp)
[1] -0.8475514
cor(mt_filt$mpg, mt_filt$disp, method = 'spearman')
[1] -0.9088824

Визуализация корреляционной матрицы

library(ggcorrplot)
correlation_matrix <- round(cor(mt_filt),1)
head(correlation_matrix[, 1:4])
      mpg disp   hp drat
mpg   1.0 -0.8 -0.8  0.7
disp -0.8  1.0  0.8 -0.7
hp   -0.8  0.8  1.0 -0.4
drat  0.7 -0.7 -0.4  1.0
wt   -0.9  0.9  0.7 -0.7
qsec  0.4 -0.4 -0.7  0.1
ggcorrplot(correlation_matrix, method ="square")

Визуализация корреляционной матрицы

ggcorrplot(correlation_matrix, method ="circle")

Визуализация корреляционной матрицы

ggcorrplot(correlation_matrix, lab =TRUE)

Визуализация корреляционной матрицы

ggcorrplot(correlation_matrix, type ="lower", lab =TRUE)

Визуализация корреляционной матрицы

ggcorrplot(correlation_matrix, type ="lower", lab =TRUE, hc.order = T)

Визуализация корреляционной матрицы

library(corrplot)
corr_mat <- cor(mt_filt)
corrplot(corr_mat)

Визуализация корреляционной матрицы

corrplot(corr_mat, col=colorRampPalette(c("blue","white","red"))(200))

Визуализация корреляционной матрицы

corrplot(corr_mat, method = 'number')

Визуализация корреляционной матрицы

corrplot(corr_mat, method = 'ellipse')

Конец!