Расчет неизвестного p-значения

Недавно я отлаживал R-скрипт и обнаружил что-то очень странное, автор определил их собственную функцию p-значения

pval <- function(x, y){
    if (x+y<20) { # x + y is small, requires R.basic
        p1<- nChooseK(x+y,x) * 2^-(x+y+1);
        p2<- nChooseK(x+y,y) * 2^-(x+y+1);
        pvalue = max(p1, p2)
    }
    else { # if x+y is large, use approximation
        log_p1 <- (x+y)*log(x+y) - x*log(x) - y*log(y) - (x+y+1)*log(2);
        pvalue<-exp(log_p1);
    }
    return(pvalue)
}

Где X и Y - положительные значения, значения которых больше 0. Случай <20, кажется, является вычислением для некоторого гипергеометрического распределения (что-то похожее на критерий Фишера?), И кто-нибудь знает, что такое другое вычисление? В качестве идентификатора я пытаюсь оптимизировать этот код, поэтому пытаюсь выяснить правильную функцию R для вызова и заменить ее на.

Редактировать: Формула детализации бумаги для расчета p-значения может быть найдена здесь (необходимо щелкнуть pdf, чтобы увидеть формулы). Методы начинаются на странице 8 pdf, а соответствующая формула может быть найдена на странице 9 в (1). Распределение, которое они предполагают, является пуассоновским.

Вторая вещь выглядит так, как будто это приближение к вычислению, используемому для x+y < 20случая, но основанное на приближении Стирлинга .

Обычно, когда он используется для такого рода аппроксимации, люди используют по крайней мере следующий дополнительный член (коэффициент в аппроксимации для ), Что значительно улучшит относительную аппроксимацию для малых , n! $\sqrt{2\pi n}$$n!$$n$

Например, если и равны 10, первое вычисление дает около 0,088, а приближение, когда коэффициент включен во все члены, составляет около 0,089, что достаточно близко для большинства целей ... но опускание этого члена в приближении дает 0,5 - что на самом деле недостаточно близко! Автор этой функции явно не удосужился проверить точность своего приближения в граничном случае.у $x$$y$$\sqrt{2\pi n}$

Для этого автору, вероятно, следовало бы просто вызвать встроенную lgammaфункцию, в частности, используя это вместо того, что он имеет для log_p1:

log_p1 <- lgamma(x+y+1)-lgamma(x+1)-lgamma(y+1)-(x+y+1)*log(2)

что приводит к ответу, который он пытается приблизить (поскольку на lgamma(x+1)самом деле возвращает , то самое, что он пытается приблизить - плохо - через приближение Стирлинга).$\log(x!)$

Точно так же я не уверен, почему автор не использует встроенную chooseфункцию в первой части, функцию, которая входит в стандартное распределение R. В этом отношении соответствующая функция распределения, вероятно, также встроена.

Вам действительно не нужны два отдельных случая; lgammaодин работает нормально , вплоть до самых маленьких значений. С другой стороны, chooseфункция работает для довольно больших значений (например, choose(1000,500)работает просто отлично). Возможно, более безопасный вариант lgamma, хотя вам нужно было бы иметь достаточно большие и прежде чем это станет проблемой.$x$$y$

Имея больше информации, можно определить источник теста. Я предполагаю, что автор взял это откуда-то, поэтому должна быть возможность отследить его. У вас есть какой-то контекст для этого?

Когда вы говорите «оптимизировать», вы имеете в виду сделать его быстрее, короче, удобнее в обслуживании или что-то еще?

Редактировать после быстрого прочтения на бумаге:

Авторы, похоже, ошибаются по ряду пунктов. Точный тест Фишера не предполагает, что поля являются фиксированными, он просто обуславливает их, что совсем не то же самое, что обсуждалось, например, здесь со ссылками. Действительно, они, кажется, в значительной степени совершенно неосведомлены о дебатах по поводу условий наценки и почему это делается. Ссылки там стоит прочитать.

[Они переходят от «теста Фишера всегда более консервативно, чем у нас» к утверждению, что тест Фишера слишком консервативен … что не обязательно следует, если оно не является условным . Они должны были бы это установить, но, учитывая, что об этом спорят статистики уже около 80 лет, и эти авторы, похоже, не знают, почему это делается, я не думаю, что эти ребята достаточно хорошо разбираются в этом вопросе. .]

Авторы статьи, по крайней мере, похоже, понимают, что вероятности, которые они дают, должны быть суммированы, чтобы дать p-значения; например около середины первого столбца страницы 5 (выделено мое):

Статистическая значимость в соответствии с точным критерием Фишера для такого результата составляет 4,6% (значение P с двумя хвостами, т. Е. Вероятность появления такой таблицы в гипотезе о том, что частоты актинового EST не зависят от библиотек кДНК). Для сравнения, значение P, рассчитанное по кумулятивной форме (уравнение 9, см. Методы) уравнения 2 (т. Е. Относительная частота актиновых EST одинакова в обеих библиотеках, учитывая, что по меньшей мере 11 родственных EST наблюдаются в библиотека печени после двух наблюдалась в библиотеке мозга) составляет 1,6%.

(хотя я не уверен, что согласен с их расчетом стоимости там; мне пришлось бы тщательно проверить, что они на самом деле делают с другим хвостом.)

Я не думаю, что программа делает это.

Остерегайтесь, однако, что их анализ не является стандартным биномиальным тестом; они используют байесовский аргумент для получения p-значения в тесте, который часто используется. Они также кажутся - как ни странно, на мой взгляд - условием для , а не для . Это означает, что у них должно получиться что-то вроде отрицательного бинома, а не бинома, но я считаю, что статья действительно плохо организована и ужасно плохо объяснена (и я привык выяснять, что происходит в статистических работах), поэтому я не собираюсь быть уверенным, если я не пройду внимательно.х + $x$$x+y$

Я даже не уверен, что сумма их вероятностей равна 1 на данный момент.

Здесь можно сказать гораздо больше, но вопрос не в бумаге, а в реализации в программе.

-

В любом случае, результат, по крайней мере, в документе правильно определяет, что p-значения состоят из суммы вероятностей, подобных тем, что приведены в уравнении 2, но программа этого не делает . (См. Уравнения 9a и 9b в разделе «Методы»).

Код просто неверен в этом.

[Вы можете использовать pbinom, как подразумевает комментарий @ whuber, для определения индивидуальных вероятностей (но не хвоста, поскольку это не биномиальный тест, как они его структурируют), но тогда в их уравнении 2 есть дополнительный коэффициент 1/2, так что если вы хотите воспроизвести результаты в документе, вам нужно изменить их.]

Вы можете получить его, немного поигравшись, из pnbinom-

Обычные формы отрицательного биномиального являются либо числа испытаний в успеха или количество отказов в успеха. Два эквивалентны; Википедия дает вторую форму здесь . Функция вероятности: к т $k^\rm{th}$$k^\rm{th}$

Уравнение 2 на p4 (а также уравнение 1 на p3) является отрицательным биномом, но смещено на 1. Пусть , и .$p = N_1/(N_1+N_2)$$k=x$$r = y+1$

Это заставляет меня беспокоиться о том, что, поскольку пределы по не были смещены аналогичным образом, их вероятности могут даже не увеличиваться до 1.$y$

Это было бы плохо.