Язык R

libraries

ggcats

install.packages("remotes")
remotes::install_github("R-CoderDotCom/ggcats@main")
library(ggcats)

grid <- expand.grid(1:5, 3:1)

ggcats

df <- data.frame(x = grid[, 1],
                 y = grid[, 2],
                 image = c("nyancat", "bongo", "colonel", "grumpy", "hipster", "lil_bub",
                           "maru", "mouth", "pop", "pop_close",  "pusheen", "pusheen_pc",
                           "toast", "venus", "shironeko"))

ggplot(df) +
   geom_cat(aes(x, y, cat = image), size = 5) +
   geom_text(aes(x, y - 0.5, label = image), size = 2.5) +
   xlim(c(0.25, 5.5)) + ylim(c(0.25, 3.5)) + theme_bw()

ggcats

iris$cat <- factor(iris$Species,
                   labels = c("pusheen", "toast", "venus"))
ggplot(iris, aes(Petal.Length, Petal.Width)) +
 geom_cat(aes(cat = cat), size = 4) + theme_bw()

ggdogs

install.packages("remotes")
remotes::install_github("R-CoderDotCom/ggdogs@main")
library(ggdogs)

ggdogs

df <- data.frame(x = grid[, 1],
                 y = grid[, 2],
                 image = c("doge", "doge_strong", "chihuahua", "eyes", "gabe", "glasses",
                           "tail", "surprised", "thisisfine", "hearing", "pug", "ears",
                           "husky", "husky_2", "chilaquil"))
                           
ggplot(df) +
 geom_dog(aes(x, y, dog = image), size = 5) +
 geom_text(aes(x, y - 0.5, label = image), size = 2.5) +
 xlim(c(0.25, 5.5)) + ylim(c(0.25, 3.5)) + theme_bw()

ggdogs

set.seed(1)
df <- data.frame(x = 1:10, y = rnorm(10))

ggplot(df, aes(x = x, y = y)) +
  geom_point(size = 3, color = 4) +
  geom_dog(aes(x = 7, y = -0.5), dog = "thisisfine", size = 5) +
  geom_label(aes(x = 7.75, y = -0.1, label = "This is fine")) + theme_bw()

CatterPlot

install.packages("remotes")
remotes::install_github("Gibbsdavidl/CatterPlots")
library(CatterPlots)

CatterPlot

x <- seq(0, 15, 0.5)
y <- sin(x)

obj <- catplot(x, y, cat = 2, catcolor = 2, type = "line")

cats(obj, 4.5, 0.8, cat = 3, catcolor = 3)

Custom dots

set.seed(123)

d <- data.frame(x = rnorm(10), y = rnorm(10), z=1:10,
  image = 'https://www.svgrepo.com/show/69771/dna.svg')

ggplot(d, aes(x, y)) + 
  geom_image(aes(image=image, color=z), size = 0.1) +
  scale_color_gradient(low='blue', high='red') +
  theme_bw()

Немного биологии

Есть ли различия?

Вопрос

Сколько же генов ассоциированы с заболеванием?

У нас же есть уже p-value, сравненные с уровнем значимости…

Другая ситуация

p-value и ошибка первого рода

Вернемся к больным

Если данные долго анализировать…

… в конце концов они в чем-то сознаются

Scanning Dead Salmon in fMRI Machine Highlights Risk of Red Herrings

Проблема множественного тестирования

Проводим много тестов

Для каждого теста контролируем уровень ошибки первого рода

Получаем ложноположительные значения

Хотелось бы как-то контролировать уровень ошибки первого рода сразу для всех проведенных тестов

Поправки на множественное тестирование

Их очень много

Выделяют два основных семейства:

• FWER - family-wise error rate

• FDR - false discovery rate

FWER

FWER

• Поправка Бонферрони

• Поправка Холм-Бонферрони: выше мощность, меньше ложноотрицательных результатов

Неважна независимость

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

Помним, что потеряем много истинноположительных

FDR

• Поправка Беньямини — Хохберга

Готовы к “мусорным” данным в результате

Порог на FDR зависит от амбиций исследователя

Задание 1

Есть ли зависимость между длиной лепестка и длиной чашелистика у ирисов?

Рисуем

Рисуем

ggplot(data = iris) + 
  geom_point(aes(x = Sepal.Length, y = Petal.Length)) + theme_bw()

cor(iris$Sepal.Length, iris$Petal.Length)

[1] 0.8717538

Пирсона или Спирмена?

cor(iris$Sepal.Length, iris$Petal.Length)

[1] 0.8717538

cor(iris$Sepal.Length, iris$Petal.Length, method = 'spearman')

[1] 0.8818981

А значимость?

cor.test(iris$Sepal.Length, iris$Petal.Length)

    Pearson's product-moment correlation

data:  iris$Sepal.Length and iris$Petal.Length
t = 21.646, df = 148, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.8270363 0.9055080
sample estimates:
      cor 
0.8717538 

cor.test(iris$Sepal.Length, iris$Petal.Length, method = 'spearman')

    Spearman's rank correlation rho

data:  iris$Sepal.Length and iris$Petal.Length
S = 66429, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.8818981 

Вопрос

Что делать, если коэффициент корреляции = 2.3 (p.value ~ -0.04)?

Матрица корреляций

cor(mtcars) -> cor.m
cor.m[1:5, 1:5]

            mpg        cyl       disp         hp       drat
mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684  0.6811719
cyl  -0.8521620  1.0000000  0.9020329  0.8324475 -0.6999381
disp -0.8475514  0.9020329  1.0000000  0.7909486 -0.7102139
hp   -0.7761684  0.8324475  0.7909486  1.0000000 -0.4487591
drat  0.6811719 -0.6999381 -0.7102139 -0.4487591  1.0000000

Рисуем

corrplot(cor(mtcars), 
         col=colorRampPalette(c("#1b4965","white","#c1121f"))(200))

Значимость?

psych::corr.test(mtcars) 

Call:psych::corr.test(x = mtcars)
Correlation matrix 
       mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
mpg   1.00 -0.85 -0.85 -0.78  0.68 -0.87  0.42  0.66  0.60  0.48 -0.55
cyl  -0.85  1.00  0.90  0.83 -0.70  0.78 -0.59 -0.81 -0.52 -0.49  0.53
disp -0.85  0.90  1.00  0.79 -0.71  0.89 -0.43 -0.71 -0.59 -0.56  0.39
hp   -0.78  0.83  0.79  1.00 -0.45  0.66 -0.71 -0.72 -0.24 -0.13  0.75
drat  0.68 -0.70 -0.71 -0.45  1.00 -0.71  0.09  0.44  0.71  0.70 -0.09
wt   -0.87  0.78  0.89  0.66 -0.71  1.00 -0.17 -0.55 -0.69 -0.58  0.43
qsec  0.42 -0.59 -0.43 -0.71  0.09 -0.17  1.00  0.74 -0.23 -0.21 -0.66
vs    0.66 -0.81 -0.71 -0.72  0.44 -0.55  0.74  1.00  0.17  0.21 -0.57
am    0.60 -0.52 -0.59 -0.24  0.71 -0.69 -0.23  0.17  1.00  0.79  0.06
gear  0.48 -0.49 -0.56 -0.13  0.70 -0.58 -0.21  0.21  0.79  1.00  0.27
carb -0.55  0.53  0.39  0.75 -0.09  0.43 -0.66 -0.57  0.06  0.27  1.00
Sample Size 
[1] 32
Probability values (Entries above the diagonal are adjusted for multiple tests.) 
      mpg cyl disp   hp drat   wt qsec   vs   am gear carb
mpg  0.00   0 0.00 0.00 0.00 0.00 0.22 0.00 0.01 0.10 0.02
cyl  0.00   0 0.00 0.00 0.00 0.00 0.01 0.00 0.04 0.08 0.04
disp 0.00   0 0.00 0.00 0.00 0.00 0.20 0.00 0.01 0.02 0.30
hp   0.00   0 0.00 0.00 0.17 0.00 0.00 0.00 1.00 1.00 0.00
drat 0.00   0 0.00 0.01 0.00 0.00 1.00 0.19 0.00 0.00 1.00
wt   0.00   0 0.00 0.00 0.00 0.00 1.00 0.02 0.00 0.01 0.20
qsec 0.02   0 0.01 0.00 0.62 0.34 0.00 0.00 1.00 1.00 0.00
vs   0.00   0 0.00 0.00 0.01 0.00 0.00 0.00 1.00 1.00 0.02
am   0.00   0 0.00 0.18 0.00 0.00 0.21 0.36 0.00 0.00 1.00
gear 0.01   0 0.00 0.49 0.00 0.00 0.24 0.26 0.00 0.00 1.00
carb 0.00   0 0.03 0.00 0.62 0.01 0.00 0.00 0.75 0.13 0.00

 To see confidence intervals of the correlations, print with the short=FALSE option

Значимость

psych::corr.test(mtcars, adjust = 'BH') -> mt.c

Рисуем

corrplot(corr = mt.c$r,
p.mat = mt.c$p,
method = "color",
order = "hclust")

Проблема

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

Что же делать?

Категориальные признаки

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

\(\chi^{2}\) - тест на сопряженность

\(H_{0}\): связи нет

\(H_{1}\): связь есть

\(\chi^{2}\) - тест на сопряженность

\(\chi^{2}\) - тест на сопряженность

\(\chi^{2}\) - тест на сопряженность

При условии верности нулевой гипотезы

Рассчитаем ожидаемое количество людей

Возьмем выборку из 116 человек

\(\chi^{2}\) - тест на сопряженность

Можно посчитать некоторую статистику - \(\chi^{2}\)

\(\chi^{2}\)

Статистика \(\chi^{2}\) распределена согласно распределению \(\chi^{2}\) с числом степеней свободы: 1

\(\chi^{2}\) - тест на сопряженность

Односторонний критерий

Интересны случаи, когда статистика больше критического значения

ref

\(\chi^{2}\) - тест на сопряженность

Рассмотрели случай таблицы 2х2

Можно проводить и на более сложных таблицах

Связаны ли количество заболеваний ОРВИ (их много: грипп, риновирусные, коронавирусные инфекции и т.д.) с временем года (4 категории)

Проблема

Ну нет у нас столько данных, чтобы получить ожидаемые нужного размера!!!

Точный тест Фишера

Наблюдения должны быть независимы

\[p(table) = \frac{(a + b)!(c + d)!(a + c)!(b + d)!}{a!b!c!d!}\]

Точный тест Фишера

Вспомним t-test

Усложним входные данные

Как поступить?

Вообще говоря…

One-Way ANOVA

head(PlantGrowth)

  weight group
1   4.17  ctrl
2   5.58  ctrl
3   5.18  ctrl
4   6.11  ctrl
5   4.50  ctrl
6   4.61  ctrl

One-Way ANOVA

library(ggpubr)
ggboxplot(PlantGrowth, x = "group", y = "weight", 
          color = "group", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("ctrl", "trt1", "trt2"),
          ylab = "Weight", xlab = "Treatment")

One-Way ANOVA

ggline(PlantGrowth, x = "group", y = "weight", 
       add = c("mean_se", "jitter"), 
       order = c("ctrl", "trt1", "trt2"),
       ylab = "Weight", xlab = "Treatment")

One-Way ANOVA

Предположение ANOVA

• Каждая выборка была взята из нормально распределенной популяции.

• Дисперсии совокупностей, из которых взяты выборки, равны

• Шум в данных распределен нормально

Конец!