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

Я нашел эту проблему с t-критериями и chi-squared в 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 десятичных значащих цифры будут хорошими, для значений, которые я собираюсь внести в Bonferroni, мне нужно больше. Когда я говорю "потерять некоторую точность", я думаю, 53 двоичных разряда, но (4) я полностью ошибаюсь, и любое значение p < .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-значения.

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

Ответ 1

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

Ответ 2

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

Прежде всего, когда я пытаюсь использовать пример оригинала 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, OP верно, что эти значения часто используются в качестве показателей прочности доказательств в литературе по биоинформатике - например, можно было бы проверить 100 000 генов-кандидатов и посмотреть на распределение результирующих значений p (поиск "графика вулкана" для одного примера такой процедуры).

Ответ 3

Два вопроса:

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

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

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

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

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

Ответ 4

Страница Wikipedia, с которой вы связались, была для типа 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)

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

Ответ 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

Ответ 6

Недавно была та же проблема. Специалист-статистик рекомендует:

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