Почему я не могу получить значение p меньше, чем 2.2e-16?

Question

Я обнаружил эту проблему с t-тестами и хи-квадратом в R, но я предполагаю, что эта проблема в основном относится к другим тестам. Если я сделаю:

a <- 1:10
b <- 100:110
t.test(a,b) 

Я получаю: t = -64.6472, df = 18.998, p-value < 2.2e-16. Из комментариев я знаю, что 2.2e-16 - это значение .Machine$double.eps - наименьшего числа с плавающей запятой, такого, что 1 + x != 1, но, конечно, R может представлять числа намного меньшие, чем это. Я также знаю из R FAQ, что R должен округлять числа с плавающей точкой до точности 53 двоичных разрядов: R FAQ .

Несколько вопросов: (1) Правильно ли я прочитал, что как 53 двоичные цифры точность или значения в R < .Machine$double.eps не рассчитаны точно? (2) Почему при выполнении таких вычислений R не предоставляет средства для отображения меньшего значения для p-значения даже с некоторой потерей точности? (3) Есть ли способ отобразить меньшее значение p, даже если я потеряю некоторую точность? Для одного теста 2 десятичных значащих числа были бы хороши, для значений, которые я собираюсь исправить в Бонферрони, мне нужно больше. Когда я говорю «потерять некоторую точность», я думаю, что <53 двоичных разряда, но (4) я полностью ошибаюсь, и любое значение p <code>< .Machine$double.eps является крайне неточным? (5) Является ли R просто честным, а другие пакеты статистики - нет?

В моем поле очень маленькие значения p являются нормой, некоторые примеры: http://www.ncbi.nlm.nih.gov/pubmed/20154341, http://www.plosgenetics.org/article/info%3Adoi%2F10.1371%2Fjournal.pgen.1002215, и именно поэтому я хочу представить такие маленькие значения p.

Спасибо за вашу помощь, извините за такой извилистый вопрос.

Answer 1

Я озадачен несколькими вещами при обмене ответами и комментариями здесь.

Прежде всего, когда я пробую оригинальный пример OP, я не получаю значение p , столь же маленькое, как обсуждаемые здесь (несколько разных версий 2.13.x и R-devel) :

a <- 1:10
b <- 10:20
t.test(a,b)
## data:  a and b 
## t = -6.862, df = 18.998, p-value = 1.513e-06

Во-вторых, когда я значительно увеличиваю разницу между группами, я на самом деле получаю результаты, предложенные @eWizardII:

a <- 1:10
b <- 110:120
(t1 <- t.test(a,b))
# data:  a and b 
# t = -79.0935, df = 18.998, p-value < 2.2e-16
#
> t1$p.value
[1] 2.138461e-25

Поведение печатного вывода в t.test определяется его вызовом stats:::print.htest (который также вызывается другими функциями статистического тестирования, такими как chisq.test, как отмечено в OP), который в свою очередь вызывает format.pval, в котором значения p меньше значения eps (по умолчанию .Machine$double.eps) как < eps. Я удивлен, что не согласен с такими обычно проницательными комментаторами ...

Наконец, хотя кажется глупым беспокоиться о точном значении очень маленького значения p , ОП правомерно, что эти значения часто используются в качестве показателей силы доказательств в литературе по биоинформатике - например, можно протестировать 100 000 генов-кандидатов и посмотреть на распределение полученных значений p (найдите «график вулкана», чтобы найти пример такой процедуры).

Answer 2

Два вопроса:

1) Какая возможная разница в статистическом значении будет между значениями p 1e-16 и 1e-32? Если вы действительно можете это оправдать, тогда используйте зарегистрированные значения.

2) Почему вы используете Википедию, когда вы заинтересованы в числовой точности R?

В R-FAQ сказано: «Другие [означающие нецелые] числа должны быть округлены до (обычно) точности 53 двоичных цифр». 16 цифр - это предел. Вот как можно получить пределы точности, когда на консоли:

> .Machine$double.eps
[1] 2.220446e-16

Это число фактически равно нулю при интерпретации в диапазоне [0,1]

Answer 3

Страница Википедии, на которую вы ссылались, была для типа Decimal64, который R не использует - он использует удвоения стандартного выпуска.

Во-первых, некоторые определения со страницы справки .Machine.

double.eps: наименьшее положительное число с плавающей запятой 'x', такое что '1 + x! = 1».... Обычно '2.220446e-16'.

double.xmin: наименьшее ненулевое нормализованное число с плавающей запятой ... Обычно '2.225074e-308'.

Таким образом, вы можете представлять числа меньше, чем 2.2e-16, но их точность уменьшается, и это вызывает проблемы с вычислениями.Попробуйте несколько примеров с числами, близкими к наименьшему представимому значению.

2e-350 - 1e-350
sqrt(1e-350)

Вы упомянули в комментарии, что хотите внести поправки бонферрони.Вместо того, чтобы использовать для этого свой собственный код, я предлагаю вам использовать p.adjust(your_p_value, method = "bonferroni").pairwise.t.test использует это.

Answer 4

Попробуйте что-то вроде этого t.test(a,b)$p.value посмотрите, дает ли это вам необходимую точность.Я считаю, что это больше связано с печатью результата, чем с фактическим сохраненным компьютерным значением, которое должно иметь необходимую точность.

Answer 5

Некоторые пакеты R решают эту проблему.Лучший способ - через пакет pspearman.

source("http://www.bioconductor.org/biocLite.R")
biocLite("pspearman")
library("pspearman")
a=c(1:110,110)
b=1:111
out <- spearman.test(a, b, alternative = "greater", approximation="t-distribution")
out$p.value

[1] 3.819961e-294

Answer 6

Недавно была такая же проблема.Товарищ по статистике рекомендует:

A <- cor.test(…)
p <- 2* pt(A$statistic,  df = A$parameter, lower.tail=FALSE)