Язык R

Что вас ждет

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

Зачем?

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Способ 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

Процентили

a <- rpois(30, 2)
quantile(a)

  0%  25%  50%  75% 100% 
0.00 2.00 2.00 3.75 6.00 

Мода

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

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 

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

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

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

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

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

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

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

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

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

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

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

• Пол ребенка

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

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

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

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

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

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

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

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

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

Нормальное распределение

Принимает значения на \((-\infty; \infty)\)

Параметры: среднее - \(\mu\); дисперсия - \(\sigma^{2}\)

Нормальное распределение

Свойства нормального распределения

Независимые случайные величины Х и Y распределены нормально с параметрами:

\(\mu_{X}; \sigma^{2}_{X}\) \(\mu_{Y}; \sigma^{2}_{Y}\)

Тогда случаная величина Z = X + Y распределена нормально с параметрами:

\(\mu_{X} + \mu_{Y}; \sigma^{2}_{X} + \sigma^{2}_{Y}\)

Свойства нормального распределения

Случайна величина X распределена нормально с параметрами:

\(\mu_{X}; \sigma^{2}_{X}\)

Тогда случайная величина \(Z = \frac{X - \mu}{\sigma}\) распределена нормально с параметрами:

\(\mu = 0; \sigma^{2} = 1\)

Центральная предельная теорема

Есть много (n) независимых случайных величин

Рапределены одинаково, есть матожидание \(\mu\) и дисперсия \(\sigma^{2}\)

\[S_{n} = \sum_{i}X_{i}\]

При очень больших n:

\[\frac{S_{n} - n\mu}{\sigma \sqrt{n}}\]

ЦПТ

Оценки параметров генеральной совокупности

Выборочное среднее: \(\overline{x} = \frac{\sum x_{i}}{n}\)

Выборочная дисперсия: \(s^{2} = \frac{1}{1-n}\sum(x_{i} - \overline{x})^{2}\)

Выборочное стандартное отклонение: \(s = \sqrt{s^{2}}\)

Что дальше?

Есть предположение как устроена генеральная совокупность

Из выборок оцениваем параметры генеральной совокупности

На основании этих оценок хотим делать выводы о генеральной совокупности

Проверять какие-то гипотезы

Нулевая гипотеза

Выдвигается, чтобы быть отвергнутой

Про неинтересный мир

Ничего не отличается друг от друга, эффекта нет, влияния нет, значение параметров равны нулю и т.д.

Нулевая гипотеза: примеры

Есть 2 группы цыплят, одни сидят на диете, другие - нет.

Гипотезы \(H_{0}\):

• среднее веса цыплят из двух групп не отличаются

• в среднем вес цыплят из первой группы отличается от веса цыплят из второй группы не более, чем на 20 г

Альтернативная гипотеза

Как правило является отрицанием \(H_{0}\)

Эффект есть, результаты значимы и т.д.

Нулевая гипотеза и альтернатива: примеры

Есть 2 группы цыплят, одни сядт на диете, другие - нет.

\(H_{0}\): среднее веса цыплят из двух групп не отличаются \(H_{1}\): среднее веса цыплят из двух групп отличаются

\(H_{0}\): в среднем вес цыплят из первой группы отличается от веса цыплят из второй группы не более, чем на 20 г \(H_{1}\): в среднем вес цыплят из первой группы отличается от веса цыплят из второй группы более, чем на 20 г

p-value

Вероятность наблюдать такой же или более критичный результат при условии верности \(H_{0}\)

Расчитываем на основе выборки

Пример

Хотим узнать, влияет ли новая диета на вес людей

Взвешиваем участников в начале исследования и через 60 дней

\(H_{0}\): в среднем вес изменяется не более, чем на 10 кг \(H_{1}\): в среднем вес изменяется более, чем на 10 кг

Получили, что вес изменился на 12 кг

Насколько этот результат значим? Диета влияет?

p-value

ref

Что дальше?

Получили “вес” нашего наблюдения

Решаем, что делать с нулевой гипотезой

Ошибка первого и второго рода

Ошибка первого и второго рода

Уровень значимости

Сравниваем p-value с каким-то числом

Заклинание: p-value меньше 0.05

Отвергаем гипотезу \(H_{0}\), поддерживая ошибку первого рода на уровне а

Какую долю ложноположительных результатов мы готовы принять

а - уровень значимости

Итого

СНАЧАЛА формулирует нулевую и альтернативную гипотезу и выбираем уровень значимости

ЗАТЕМ проводим эксперимент

Сравниваем p-value с уровнем значимости

p-value меньше уровня значимости есть основания отвергнуть \(H_{0}\)

Одновыборочный тест Стьюдента

НО!!!! Одна выборка это скучно

Обычно у нас 2 выборки

Выборки

Двувыборочный тест Стьюдента

Формулируем гипотезы

\(H_{0}: \mu_{x} - \mu_{y} = 0\)

\(H_{1}: \mu_{x} - \mu_{y} \not= 0\)

\(H_{1}: \mu_{x} - \mu_{y} > 0\)

\(H_{1}: \mu_{x} - \mu_{y} < 0\)

Формулируем более интересные гипотезы

\(H_{0}: \mu_{x} - \mu_{y} = а\)

\(H_{1}: \mu_{x} - \mu_{y} \not= а\)

\(H_{1}: \mu_{x} - \mu_{y} > а\)

\(H_{1}: \mu_{x} - \mu_{y} < а\)

Тест Велча

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

“Закручиваем гайки” - уменьшеаем число степеней свободы

Далее - аналогичные рассуждения о нулевых и альтернативных гипотезах

Выборки

Парный тест Стьюдента

Дано: две выборки из разных генеральных совокупностей

Наблюдения из этих выборок связаны в пары

Пример: измеряем вес крыс до и после диеты

Не важно знать о дисперсиях генеральных совокупностей

Интересно исследовать среднее разниц!

Парный тест Стьюдента: гипотезы

Нулевая гипотеза: среднее разниц двух выборок не отличается от а

\(H_{0}: \mu_{x-y} = а\)

Альтернативные гипотезы:

\(H_{1}: \mu_{x-y} \not= а\)

\(H_{1}: \mu_{x-y} > а\)

\(H_{1}: \mu_{x-y} < а\)

Задание №1

Смоделируйте скошенное распределение любым способом (положите в вектор)

Визуализируйте распределение, выбрав наиболее оптимальный тип графика

Решение 1.1

set.seed(123)
a <- c(rnorm(10000), rnorm(500, 4, 3), rnorm(500, 6, 5))

Решение 1.2

ggplot(mapping = aes(x = a)) + 
  geom_histogram(fill = '#598392', color = '#01161e') + 
  theme_bw()

Усеченное среднее

mean(a)

[1] 0.4604595

mean(a, trim = .2)

[1] 0.1117887

Асимметрия и эксцесс

skew(a)

[1] 3.24231

kurtosi(a)

[1] 15.81482

Асимметрия и эксцесс: примеры

Сводное описание

summary(a)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-7.9186 -0.6191  0.0928  0.4605  0.8825 20.0108 

describe(a) # library(psych)

   vars     n mean   sd median trimmed  mad   min   max range skew kurtosis
X1    1 11000 0.46 2.18   0.09    0.13 1.11 -7.92 20.01 27.93 3.24    15.81
     se
X1 0.02

describe + group_by/summarise

iris %>%
  group_by(Species) %>%
  summarise(describe(Petal.Width))

# A tibble: 3 × 14
  Species  vars     n  mean    sd median trimmed   mad   min   max range    skew
  <fct>   <dbl> <dbl> <dbl> <dbl>  <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
1 setosa      1    50 0.246 0.105    0.2   0.238 0       0.1   0.6   0.5  1.18  
2 versic…     1    50 1.33  0.198    1.3   1.32  0.222   1     1.8   0.8 -0.0293
3 virgin…     1    50 2.03  0.275    2     2.03  0.297   1.4   2.5   1.1 -0.122 
# ℹ 2 more variables: kurtosis <dbl>, se <dbl>

library(skimr)

qqplot

set.seed(1)
my_df <- data.frame(y=rexp(200, rate=3))
ggplot(my_df, mapping = aes(x = y)) + geom_histogram(fill = '#598392', color = '#01161e') + theme_bw()

qqplot

ggplot(my_df, aes(sample=y)) +
  stat_qq() + 
  stat_qq_line() + theme_bw()

qqplot

Задание 2

Отличается ли в среднем длина лепестков у ирисов от 4 см?

Сначала рисуйте!

Нарисуем!

ggplot(iris) +
  geom_histogram(aes(x = Petal.Length, fill = Species)) + theme_bw()

Какие идеи?

Подготовим данные

df <- iris %>%
  filter(Species %in% c('versicolor', 'virginica'))

Решение 2.1

t.test(df$Petal.Length)

    One Sample t-test

data:  df$Petal.Length
t = 59.425, df = 99, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 4.742187 5.069813
sample estimates:
mean of x 
    4.906 

Решение 2.2

t.test(df$Petal.Length, mu = 4)

    One Sample t-test

data:  df$Petal.Length
t = 10.974, df = 99, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 4
95 percent confidence interval:
 4.742187 5.069813
sample estimates:
mean of x 
    4.906 

Без интерпретации результата - ответа нет!

Неверная интерпретация - ответа нет!

Вывод результата теста

t.test(df$Petal.Length, mu = 4) -> x
x$p.value

[1] 8.307922e-19

Задание 3

Больше ли в среднем длина лепестков у ирисов, чем 3 см?

Какая нулевая гипотеза?

Решение 3.1

t.test(df$Petal.Length, mu = 3, alternative = 'greater')

    One Sample t-test

data:  df$Petal.Length
t = 23.087, df = 99, p-value < 2.2e-16
alternative hypothesis: true mean is greater than 3
95 percent confidence interval:
 4.768922      Inf
sample estimates:
mean of x 
    4.906 

Задание 4

Различаются ли в среднем длины лепестков у ирисов видов versicolor и virginica?

Какая нулевая гипотеза?

Решение 4.1

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

Теперь визуализируйте правильно

Решение 4.2

ggerrorplot(df, x = "Species", y = "Petal.Length", 
            desc_stat = "mean_sd", color = "black",
            add = "jitter", add.params = list(color = "darkgray"), size = 0.8
            )

Решение 4.3

a <- filter(df, Species == 'versicolor')
b <- filter(df, Species == 'virginica')

t.test(a$Petal.Length, b$Petal.Length)

    Welch Two Sample t-test

data:  a$Petal.Length and b$Petal.Length
t = -12.604, df = 95.57, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.49549 -1.08851
sample estimates:
mean of x mean of y 
    4.260     5.552 

Решение 4.4

Лучше

t.test(Petal.Length ~ Species, df)

    Welch Two Sample t-test

data:  Petal.Length by Species
t = -12.604, df = 95.57, p-value < 2.2e-16
alternative hypothesis: true difference in means between group versicolor and group virginica is not equal to 0
95 percent confidence interval:
 -1.49549 -1.08851
sample estimates:
mean in group versicolor  mean in group virginica 
                   4.260                    5.552 

Решение 4.5

Теперь визуализируем еще лучше

ggbarplot(df, x = "Species", y = "Petal.Length", 
          fill = ("Species"), 
          add = c('mean', 'jitter'),
          palette = c("#9e2a2b", "#62929e")) +
  stat_compare_means(method = "t.test")

Конец!