Достижение нижнего значения P-значений в R t.test
У меня есть три вектора из ~8000 генов, каждый из которых имеет соответствующую частоту мутаций, например:
alpha = c(0.84, 0.87, 0.91...)
beta = c(0.97, 0.94, 0.99...)
kappa = c(0.72, 0.68, 0.75...)
Я использую функцию t.test в R для вычисления P-значения между этими векторами. Из-за того, как много генов и насколько постоянны различия между каждым вектором, мои значения P всегда очень низкие. Например, когда я запускаю:
ab = t.test(x = alpha, y = beta)
ak = t.test(x = alpha, y = kappa)
bk = t.test(x = beta, y = kappa)
В сводках всех тестов выводится одно и то же значение 2.2e-16, что мне не нравится. Когда я извлекаю значение P напрямую:
t.test(alpha, kappa)$p.value
R просто выводит 0. Я предполагаю, что это связано с каким-то внутренним ограничением размера float/double, но это выглядит не очень хорошо. Мы публикуем эти данные, и независимо от того, возьмем ли мы суммарные P-значения, которые все одинаковы, или равные 0, ни то, ни другое не выглядит хорошо. Есть ли способ обойти это в R, чтобы я мог вычислять сколь угодно низкие значения P? Если бы другой инструмент работал лучше, я бы тоже был в нем заинтересован. Спасибо!
Вы правы, что в сводке выводится 2.2e-16. Однако, когда я извлекаю значение p напрямую с помощью t.test(alpha, kapp)$p.value, R выводит 0. Я добавлю это к вопросу для ясности.
Попробуйте dput(t.test(alpha, kappa)$p.value) или print(t.test(alpha, kappa)$p.value, digits=16). Было бы полезно включить конкретный воспроизводимый пример. Избегайте ... в примерах кода, поскольку это означает, что мы не сможем скопировать/вставить значения в R для тестирования.
...или попробуйте использовать var.equal=TRUE, например. t.test(alpha, kappa, var.equal=TRUE)$p.value на все тесты.
Также возможно, что серия t-тестов не является подходящим статистическим инструментом для этих данных. Размеры вашей выборки очень велики (около 8000 на ген), и вы проводите множественные сравнения без даже грубой поправки. Мы также не знаем, являются ли эти значения взаимозависимыми, ненормальными и т. д. Я бы предложил спросить о целесообразности вашего подхода на CrossValidated.
Ответы 1
Если ваша t-статистика действительно настолько огромна, что ваше значение p меньше нуля (т. е. значение p меньше, чем примерно 1e-308 (!!), и, следовательно, его нельзя отличить от нуля в стандартных 64-битных точность с плавающей запятой), вы можете адаптировать механизм из этого ответа, чтобы получить значение p...
Я не могу (или не хочу) не отметить, что такие экстремальные значения p, безусловно, не имеют смысла в обычном статистическом смысле (Эндрю Гельман язвительно в своем блоге о p-значениях 10^(-246), которые на тысячи порядков больше тех, о которых вы сообщаете в комментариях...). Хотя есть некоторая ценность в соответствии с соглашениями вашей области, возможно, было бы лучше просто сообщить t-статистику, которая даст вам разумную меру разницы между группами (или просто чистую разницу в средних значениях), и оставить вывести (бессмысленные) p-значения?
x <- rep(1:5, 100)
y <- rep(10001:10005, 100)
tt <- t.test(x,y)
Результаты:
Welch Two Sample t-test
data: x and y
t = -111692, df = 998, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-10000.176 -9999.824
sample estimates:
mean of x mean of y
3 10003
tt$p.value ## 0
pvalue.extreme <- function(t, df) {
log.pvalue <- log(2) + pt(abs(t), df = df,
lower.tail = FALSE, log.p = TRUE)
log10.pvalue <- log.pvalue/log(10) ## from natural log to log10
mantissa <- 10^(log10.pvalue %% 1)
exponent <- log10.pvalue %/% 1
## or return(c(mantissa,exponent))
return(sprintf("p value is %1.2f times 10^(%d)",mantissa,exponent))
}
pvalue.extreme(tt$statistic, tt$parameter)
[1] "p value is 1.11 times 10^(-3543)"
Спасибо, это хорошо сработало для моих данных. Как оказалось, выходное значение P этого метода равно 2,20e-8005.
Не знаю, откуда вы взяли форму номер 2.2e-16, но кажется, вы смотрите на резюме, где, вероятно, написано
p-value < 2.2e-16. Если вы явно выберете значение p, это будет действительное число, например.t.test(alpha, kappa)$p.value=>[1] 7.595779e-295, кстати, это бессмысленное различие.