Язык R

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

Лекция 11

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

15 ноября 2024

Проблема

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

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

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

Курит Не курит
Больной 25 11
Здоровый 32 48

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

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

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

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

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

Курит Не курит Вероятность встретить
Больной 25 11 Больного: 36 / 116
Здоровый 32 48 Здорового 80 / 116
Вероятность встретить Курящего: 57 / 116 Некурящего: 59 / 116 Всего: 116

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

Курит Не курит Вероятность встретить
Больной 0.31 * 0.59 0.31 * 0.51 Больного: 0.31
Здоровый 0.68 * 0.49 0.68 * 0.51 Здорового 0.68
Вероятность встретить Курящего: 0.49 Некурящего: 0.51 Всего: 116

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

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

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

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

Курит Не курит
Больной 0.18 * 116 0.16 * 116
Здоровый 0.33 * 116 0.35 * 116

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

Наблюдаемое Курит Не курит
Больной 25 11
Здоровый 32 48
Ожидаемое Курит Не курит
Больной 20.88 18.56
Здоровый 46.4 40.6

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

\(\chi^{2}\)

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

Ожидаемое число в каждой ячейке должно быть не менее 5

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

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

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

ref

Таблица 2х2

library(vcd)

data1 <- as.table(matrix(c(25, 32, 11, 48), 2, 2))

data1
   A  B
A 25 11
B 32 48
dimnames(data1) <- list(Measurement_1 = c("Больной", "Здоровый"),
                       Measurement_2 = c("Курит", "Не курит"))

data1
             Measurement_2
Measurement_1 Курит Не курит
     Больной     25       11
     Здоровый    32       48

chisq.test

chisq.test(data1)

    Pearson's Chi-squared test with Yates' continuity correction

data:  data1
X-squared = 7.4747, df = 1, p-value = 0.006257

Таблица 2х2 (%)

percentages <- round(100*prop.table(data1), 2)

percentages
             Measurement_2
Measurement_1 Курит Не курит
     Больной  21.55     9.48
     Здоровый 27.59    41.38

Таблица 2х2 (%)

etiquettes <- as.table(matrix(paste0(data1, "; ", percentages, "%"), 2, 2))

dimnames(etiquettes) <- dimnames(data1)

etiquettes
             Measurement_2
Measurement_1 Курит      Не курит  
     Больной  25; 21.55% 11; 9.48% 
     Здоровый 32; 27.59% 48; 41.38%

Мозаичный график

vcd::mosaic(data1, 
       pop = F, shade = T, colorize = T, 
       gp = gpar(fill = matrix(c("grey", "#756bb1", "#756bb1", "grey"), 2, 2))) 

labeling_cells(text = etiquettes, margin = 0)(data1)

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

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

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

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

Проблема

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

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

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

fisher.test()

Курит Не курит Всего
Больной A B A + B
Здоровый C D C + D
Всего A + C B + D A + B + C + D

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

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

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

Можно рассмотреть всевозможные таблицы, перекошенные также как наша таблица или еще больше в ту или другую сторону

Для каждой такой таблицы считаем вероятность

Складываем вероятности

Получаем p-value

Решаем, есть ли основание отклонить нулевую гипотезу

Форма связи между переменными

cars <- mtcars[, 1:7]
head(cars)
                   mpg cyl disp  hp drat    wt  qsec
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02
Datsun 710        22.8   4  108  93 3.85 2.320 18.61
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02
Valiant           18.1   6  225 105 2.76 3.460 20.22

pairs(cars)

pairs(cars, panel = panel.smooth)

library(GGally) 
ggpairs( cars )

ggpairs(iris)

Вспомним t-test

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

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

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

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

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

Мы пришли в спортзал, где есть групповые занятия по:

  • йоге

  • аэробике

  • тренировке на тренажерах

Хотим понять, какой вид тренировок больше способствует потере веса.

Уже три группы

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

Давайте посравниваем попарно все группы между собой

  • йога vs аэробика

  • йога vs тренировка

  • …

Чем сравнить? t-test

Много t-test`ов

А еще поправка…

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

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

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

Однофакторый дисперсионный анализ - One-Way ANOVA

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

Насколько эти отличия значимы?

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
PG <- PlantGrowth

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

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

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

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

PG$group <- factor(PG$group,
                         levels = c("ctrl", "trt1", "trt2"))
levels(PG$group)
[1] "ctrl" "trt1" "trt2"

library(tidyverse)
group_by(PG, group) %>%
  summarise(
    count = n(),
    mean = mean(weight, na.rm = TRUE),
    sd = sd(weight, na.rm = TRUE)
  )
# A tibble: 3 × 4
  group count  mean    sd
  <fct> <int> <dbl> <dbl>
1 ctrl     10  5.03 0.583
2 trt1     10  4.66 0.794
3 trt2     10  5.53 0.443
res.aov <- aov(weight ~ group, data = PG)
summary(res.aov)
            Df Sum Sq Mean Sq F value Pr(>F)  
group        2  3.766  1.8832   4.846 0.0159 *
Residuals   27 10.492  0.3886                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(res.aov)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = weight ~ group, data = PG)

$group
            diff        lwr       upr     p adj
trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
trt2-trt1  0.865  0.1737839 1.5562161 0.0120064

pairwise.t.test(PG$weight, PG$group,
                 p.adjust.method = "BH")

    Pairwise comparisons using t tests with pooled SD 

data:  PG$weight and PG$group 

     ctrl  trt1 
trt1 0.194 -    
trt2 0.132 0.013

P value adjustment method: BH

Однородность дисперсий

plot(res.aov, 1)

plot(res.aov, 2)

aov_residuals <- residuals(object = res.aov )
shapiro.test(x = aov_residuals )

    Shapiro-Wilk normality test

data:  aov_residuals
W = 0.96607, p-value = 0.4379

kruskal.test(weight ~ group, data = PG)

    Kruskal-Wallis rank sum test

data:  weight by group
Kruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842

Конец!