Когда мы выполняем много статистических тестов, некоторые из них дадут нам статистически значимый результат совершенно случайно. Даже если в действительности никаких различий нет.
К примеру при \(\alpha = 0.05\), из 1000 тестов мы получим примерно 50 результатов, говорящих о достоверных отличиях, даже выборки в действительности не различаются.
Плохо? Плохо!
Familywise error rate – допустимая частота ошибок I рода для всего семейства тестов.
Обычно задается достаточно жестко: например при тестировании лекарств \(FWER \le 0.05\) на всю батарею тестов. То есть вероятность совершить хотя бы одну ошибку I рода составляет 0.05.
Для этого перед сравнением с критическим значением, в результаты каждого отдельного теста вносится поправка.
Исторически первый метод корректировки. Выполняется в один шаг.
Полученные в каждом тесте р-значения (\(p_i\)) нужно умножить на общее число тестов (m) и только после этого сравнить с критическим значением:
\(p_i * m \le \alpha\)
Например:
P-values (sorted): 0.005, 0.011, 0.02, 0.04, 0.13
m = 5; Significance level: 0.05
Corrected p-value: 0.005 * 5 = 0.025 < 0.05 Reject
Corrected p-value: 0.011 * 5 = 0.055 Don’t reject
Corrected p-value: 0.02 * 5 = 0.1 Don’t reject
Corrected p-value: 0.04 * 5 = 0.2 Don’t reject
Corrected p-value: 0.13 * 5 = 0.65 Don’t reject
Плюс - теперь не так-то просто совершить хотя бы одну ошибку I рода на всю группу тестов.
Минус - выросли шансы совершить ошибку II рода, т.е. "потерять" реально существующие различия.
Более щадящий способ это метод Хольма-Бонферрони.
Процедура применения поправки пошаговая. Сначала мы сортируем p-values, полученные в тестах, в порядке возрастания и присваиваем каждой соответствующий ранг j от 1 до m.
Затем применяем поправку \(p_i * (m – j + 1)\), после чего сравниваем с уровнем значимости.
Например:
P-values (sorted): 0.005, 0.011, 0.02, 0.04, 0.13
m = 5; Significance level: 0.05
j = 1 Corrected p-value: 0.005 * (5 – 1 +1) = 0.025 < 0.05: Reject
j = 2 Corrected p-value: 0.011 * (5 – 2 + 1) = 0.044: Reject
j = 3 Corrected p-value: 0.02 * (5 – 3 + 1) = 0.06: Don’t reject
Дальше можно не проверять: все следующие \(p_j\) будут больше \(\alpha\).
j = 4 Corrected p-value: 0.04 * 2 = 0.08: Don’t reject
j = 5 Corrected p-value: 0.13 * 1 = 0.13: Don’t reject
Для случаев, когда нам важнее сохранить истинно-положительные результаты, чем не допустить ложно-положительных, можно регулировать False Discovery Rate.
С помощью FDR мы задаем не количество ошибок первого рода в принципе, а количество ложно-положительных результатов в отношении к истинно-положительным.
Эта цифра обозначается как \(\gamma\) и обычно 10% или меньше считается вполне приемлемой величиной.
Процедура Беньямини-Хохберга пошаговая. Начинаем опять с сортировки и придания рангов р-значениям.
Затем находим такое p-value с наибольшим рангом j чтобы \(p_j \le \gamma * (j/m)\)
Все тесты с рангами меньше j считаем значимыми.
Например:
P-values (sorted): 0.005, 0.011, 0.02, 0.04, 0.13
m = 5; Significance level: 0.05
j=1 0.005 < 0.1*(1/5) = 0.02 Reject
j=2 0.011 < 0.1*(2/5) = 0.04 Reject
j=3 0.02 < 0.1*(3/5) = 0.06 Reject
j=4 0.04 < 0.1*(4/5) = 0.08 Reject
j=5 0.13 > 0.1*(5/5) = 0.1 Don't reject
Как видим, это либеральный метод. Да, здесь могут быть ложноположительные результаты, но мы к этому готовы. Главное – истинно-положительные результаты не пропали.
Прочитаем таблицу симулированных данных и посмотрим как они устроены:
expr <- read.table("data/fake_expression_samples.csv", header=T, sep = ",")
head(expr)
## X gene_ID sample1 sample2 ## 1 1 1 7.83 11.5 ## 2 2 1 11.87 10.5 ## 3 3 1 3.90 11.3 ## 4 4 1 8.56 11.3 ## 5 5 1 4.95 14.9 ## 6 6 1 7.46 10.2
str(expr)
## 'data.frame': 150 obs. of 4 variables: ## $ X : int 1 2 3 4 5 6 7 8 9 10 ... ## $ gene_ID: int 1 1 1 1 1 1 1 1 1 1 ... ## $ sample1: num 7.83 11.87 3.9 8.56 4.95 ... ## $ sample2: num 11.5 10.5 11.3 11.3 14.9 ...
Нам надо сравнить выборки 1 и 2 по каждому из 10 генов. Начнем с первого
t.test(expr$sample1[expr$gene_ID=="1"],
expr$sample2[expr$gene_ID=="1"])
## ## Welch Two Sample t-test ## ## data: expr$sample1[expr$gene_ID == "1"] and expr$sample2[expr$gene_ID == "1"] ## t = -6, df = 20, p-value = 0.000002 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## -5.36 -2.67 ## sample estimates: ## mean of x mean of y ## 7.54 11.56
Нам понадобятся только p-value, поэтому хорошо если мы сможем ее извлечь и записать в новую переменную
t_result <- t.test(expr$sample1[expr$gene_ID=="1"],
expr$sample2[expr$gene_ID=="1"])
t_result$p.value
## [1] 0.00000198
Мы можем добиться того же самого и одной командой:
t1 <- t.test(expr$sample1[expr$gene_ID=="1"],
expr$sample2[expr$gene_ID=="1"])$p.value
t1
## [1] 0.00000198
Теперь надо лишь выполнить эту же операцию для оставшихся девяти сравнений
t2 <- t.test(expr$sample1[expr$gene_ID=="2"],
expr$sample2[expr$gene_ID=="2"])$p.value
t3 <- t.test(expr$sample1[expr$gene_ID=="3"],
expr$sample2[expr$gene_ID=="3"])$p.value
t4 <- t.test(expr$sample1[expr$gene_ID=="4"],
expr$sample2[expr$gene_ID=="4"])$p.value
t5 <- t.test(expr$sample1[expr$gene_ID=="5"],
expr$sample2[expr$gene_ID=="5"])$p.value
t6 <- t.test(expr$sample1[expr$gene_ID=="6"],
expr$sample2[expr$gene_ID=="6"])$p.value
t7 <- t.test(expr$sample1[expr$gene_ID=="7"],
expr$sample2[expr$gene_ID=="7"])$p.value
t8 <- t.test(expr$sample1[expr$gene_ID=="8"],
expr$sample2[expr$gene_ID=="8"])$p.value
t9 <- t.test(expr$sample1[expr$gene_ID=="9"],
expr$sample2[expr$gene_ID=="9"])$p.value
t10 <- t.test(expr$sample1[expr$gene_ID=="10"],
expr$sample2[expr$gene_ID=="10"])$p.value
Мы посчитали десять p-values для всех генов. Соберем их в один вектор и присвоим имена.
pvals <- c(t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
names(pvals) <- c("t1","t2","t3","t4","t5","t6","t7","t8","t9","t10")
pvals
## t1 t2 t3 t4 t5 ## 0.000001985 0.008085760 0.004916486 0.035231583 0.000000128 ## t6 t7 t8 t9 t10 ## 0.000038369 0.478963037 0.240548497 0.028976352 0.564333804
Сколько генов, достоверно меняющих экспрессию, мы нашли в "сырых" t-тестах?
sum(pvals <= 0.05)
## [1] 7
А теперь внесем поправки.
Поправки к p-values в R можно сделать при помощи функции p.adjust() Аргумент method задает тип поправки.
p_bonf <- p.adjust(pvals, method = "bonferroni")
У скольких генов экспрессия достоверно различается после поправки Бонферрони?
sum(p_bonf <= 0.05, na.rm = TRUE)
## [1] 4
После поправки Хольма?
p_holm <- p.adjust(pvals, method = "holm") sum(p_holm <= 0.05, na.rm = TRUE)
## [1] 5
А сколько — после применения процедуры Беньямини-Хохберга?
p_bh <- p.adjust(pvals, method = "BH") sum(p_bh <= 0.05, na.rm = TRUE)
## [1] 6
Так же можем посмотреть результаты каких именно тестов показали значимые различия:
names(p_bh[p_bh <= 0.05])
## [1] "t1" "t2" "t3" "t5" "t6" "t9"
Синдром Кушинга это сборное название для ряда гормональных заболеваний разной этиологии, которые ведут к перевыработке гормона кортизола.
Датасет Cushings содержит данные по объему выведения с мочой двух метаболитов стероидных гормонов.
Cushings <- read_excel("data/Cushings.xlsx", sheet = 1)
head(Cushings)
str(Cushings)
## Classes 'tbl_df', 'tbl' and 'data.frame': 27 obs. of 3 variables: ## $ Tetrahydrocortisone: num 3.1 3 1.9 3.8 4.1 1.9 8.3 3.8 3.9 7.8 ... ## $ Pregnanetriol : num 11.7 1.3 0.1 0.04 1.1 0.4 1 0.2 0.6 1.2 ... ## $ Type : chr "a" "a" "a" "a" ...
Как видим, в таблице данных три переменные:
Tetrahydrocortisone – объем (мг в сутки) выведения тетрагидрокортизона
Pregnanetriol – объем (мг в сутки) выведения прегнантриола
Type – разновидности синдрома: a = adenoma, b = bilateral hyperplasia, c = carcinoma и u = unknown.
Названия переменных слишком замысловатые, переназовем их покороче.
names(Cushings)[names(Cushings)=="Tetrahydrocortisone"] <- "Tetr" names(Cushings)[names(Cushings)=="Pregnanetriol"] <- "Pre"
Ваша задача - выбрать один из метаболитов, и сравнить уровни его экскреции между всеми типами синдромов. Какие типы заболеваний вы стали бы исследовать дальше и почему?
a.b <- t.test(Cushings$Tetr[Cushings$Type=="a"],
Cushings$Tetr[Cushings$Type=="b"])$p.value
a.c <- t.test(Cushings$Tetr[Cushings$Type=="a"],
Cushings$Tetr[Cushings$Type=="c"])$p.value
a.u <- t.test(Cushings$Tetr[Cushings$Type=="a"],
Cushings$Tetr[Cushings$Type=="u"])$p.value
b.c <- t.test(Cushings$Tetr[Cushings$Type=="b"],
Cushings$Tetr[Cushings$Type=="c"])$p.value
b.u <- t.test(Cushings$Tetr[Cushings$Type=="b"],
Cushings$Tetr[Cushings$Type=="u"])$p.value
c.u <- t.test(Cushings$Tetr[Cushings$Type=="c"],
Cushings$Tetr[Cushings$Type=="u"])$p.value
Tetr_pvals <- c(a.b, a.c, a.u, b.c, b.u, c.u)
names(Tetr_pvals) <- c("a.b", "a.c", "a.u", "b.c", "b.u", "c.u")
names(Tetr_pvals[Tetr_pvals <= 0.05])
## [1] "a.b" "a.u"
p_bonf <- p.adjust(Tetr_pvals, method = "bonferroni") names(p_bonf[p_bonf <= 0.05])
## [1] "a.b"
p_holm <- p.adjust(Tetr_pvals, method = "holm") names(p_holm[p_holm <= 0.05])
## [1] "a.b"
p_bh <- p.adjust(Tetr_pvals, method = "BH") names(p_bh[p_bh <= 0.05])
## [1] "a.b"
a_b <- t.test(Cushings$Pre[Cushings$Type=="a"],
Cushings$Pre[Cushings$Type=="b"])$p.value
a_c <- t.test(Cushings$Pre[Cushings$Type=="a"],
Cushings$Pre[Cushings$Type=="c"])$p.value
a_u <- t.test(Cushings$Pre[Cushings$Type=="a"],
Cushings$Pre[Cushings$Type=="u"])$p.value
b_c <- t.test(Cushings$Pre[Cushings$Type=="b"],
Cushings$Pre[Cushings$Type=="c"])$p.value
b_u <- t.test(Cushings$Pre[Cushings$Type=="b"],
Cushings$Pre[Cushings$Type=="u"])$p.value
c_u <- t.test(Cushings$Pre[Cushings$Type=="c"],
Cushings$Pre[Cushings$Type=="u"])$p.value
Pre_pvals <- c(a_b, a_c, a_u, b_c, b_u, c_u)
names(Pre_pvals) <- c("a_b", "a_c", "a_u", "b_c", "b_u", "c_u")
names(Pre_pvals[Pre_pvals <= 0.05])
## [1] "b_c" "c_u"
p_bonf <- p.adjust(Pre_pvals, method = "bonferroni") names(p_bonf[p_bonf <= 0.05])
## character(0)
p_holm <- p.adjust(Pre_pvals, method = "holm") names(p_holm[p_holm <= 0.05])
## character(0)
p_bh <- p.adjust(Pre_pvals, method = "BH") names(p_bh[p_bh <= 0.05])
## [1] "b_c" "c_u"
Нельзя безнаказанно проводить множество тестов! Число ложно-положительных результатов растет очень быстро. Выбирая поправку, решите что вам важнее, избавиться от ложно-положительных, или сберечь истинно-положительные результаты.