Точный тест Фишера. Python против R

[apximhd]

Я понимаю, что для ответа на мой вопрос самым правильным было бы посмотреть исходники R, но прежде чем погрузиться в разбор чужого кода, я хотел бы задать вопрос здесь: а вдруг кто-нибудь из моих друзей уже знает ответ? Проблема возникла из-за того, что мне понадобилось заменить скрипт на R скриптом на Perl. Реализация функции, вычисляющей точный тест Фишера, есть в библиотеке Python и в библиотеке R, но её нет в Perl. Поэтому я написал следующий тестовый код на Perl (он приводится под катом).
Результаты работы этого теста совпадают с результатами работы стандартной функции из библиотеки Python, но расходятся с результатами, выдаваемыми стандартной функцией из R.
Тестовый вывод приведенной под катом программы:

A = |30 60|
    |10 70|
Fisher exact test, Perl
P two-sided =         0.00187
Odds ratio  =         3.500
Confidence Interval:  1.581...7.746
P one-sided less =    0.99972
P one sided greater = 0.00109
P hypergeom =         0.00081

Для сравнения привожу тестовый вывод с теми же данными результатов работы стандартной библиотечной функции Python stats.fisher_exact:

Python's library function stats.fisher_exact
P two-sided =         0.00187
Odds ratio  =         3.500
P one-sided less =    0.99972
P one-sided greater = 0.00109

Однако результаты работы функции fisher.test из библиотеки R несколько отличаются от реализации на Perl, выполненной в соответствии с описанием в учебниках, и от реализации в статистической библиотеке Python. Отличающиеся значения выделены красным цветом:

R: Fisher's Exact Test for Count Data
P two-sided =         0.001871
Odds ratio  =         3.475
Confidence Interval:  1.503...8.655
P one-sided less =    0.99972
P one-sided greater = 0.00109

В документации R написано, что при расчете отношения шансов и доверительного интервала используется дополнительная «подгонка» результатов с использованием метода наилучшего соответствия. Есть ещё вот такая статья по расчету отношения шансов: http://www.people.fas.harvard.edu/~mparzen/published/parzen17.pdf, но в R всё равно используется что-то другое. Может ли кто-нибудь вкратце описать алгоритм работы функции R или всё же придётся лезть в исходники?

Upd.: Вкратце: мне нужно понять, как «руками» получить такое же отношение шансов (odds ratio), которое возвращает функция R, потому что сами вероятности я считаю правильно.

#!/usr/bin/env perl
use Time::HiRes qw(time);
use Math::Trig qw(pi);
use List::Util qw(min max);


my @A = ([30,60],[10,70]);
print "A = |@$_|\n" for $A[0];
print "    |@$_|\n" for $A[1];
printf
"Fisher exact test, Perl
P two-sided =         %.5f
Odds ratio  =         %.3f
Confidence Interval:  %.3f...%.3f
P one-sided less =    %.5f
P one sided greater = %.5f
P hypergeom =         %.5f
\n",
fisher_exact_test(@A);

sub stirling { # Вычисление логарифма факториала с использованием приближения Стирлинга
    my $n = $_[0]; 
    if ($n == 0 or $n == 1) {
        return 0;
    }
    if ($n <= 16) { # Если n <= 16, вычисляем факториал точно
        my $factorial = 1; my $i = 1;
        $factorial *= ++$i while $i < $n;
        return log($factorial);
    }
    my $n2 = $n**2;
    return $n*log($n) - $n + 1/2*log(2*pi*$n) + (1/12-1/360/$n2)/$n;
}


sub logbinomial { # Вычисление логарифма биномиального коэффициента
    my ($n, $k) = @_;
    return stirling($n) - stirling($k) - stirling($n-$k);
}


sub hypergeom { # Вычисление вероятностного параметра гипергеометрического распределения
    my ($a, $b, $c, $d) = @_;
    return exp(logbinomial($a+$b, $a)
             + logbinomial($c+$d, $c)
             - logbinomial($a+$b+$c+$d, $a+$c));
}


sub fisher_exact_test { # Основная функция, вычисляющая параметры теста Фишера
                        # Входные параметры -- матрица 2x2
                        # Выходные параметры:
                        # $p_value -- двухсторонний вероятностный критерий 
                        # $odds_ratio -- отношение шансов
                        # $ci95_low -- нижняя граница 95% доверительного интервала
                        # $ci95_high -- верхняя граница 95% доверительного интервала
                        # $l_value -- односторонний вероятностный критерий «слева» или «меньше»
                        # $r_value -- односторонний вероятностный критерий «справа» или «больше»
                        # $m_value -- «гипергеометрическая» вероятность
    my ($a, $b, $c, $d) = map @$_, @_; # Разбор переданной матрицы на отдельные переменные
    my $i_min = max (0, $a-$d);
    my $i_max = min ($a+$b, $a+$c);
    my $p_value = 0; # Двухсторонний вероятностный критерий
    if ($i_min == $i_max) {
        my p_value = 1;
    }
    my $m_value = hypergeom($a, $c, $b, $d); # «Гипергеометрическая» вероятность
    my $l_value = 0; # Односторонний вероятностный критерий «слева» или «меньше»
    my $r_value = 0; # Односторонний вероятностный критерий «справа» или «больше»
    for (my $i = $i_min; $i <= $i_max; $i = $i + 1) {
        $p = hypergeom($i, $a+$c-$i, $a+$b-$i, $d-$a+$i);
        if ($i <= $a) {
            $l_value += $p;
        }
        if ($i >= $a) {
            $r_value += $p;
        }
        if ($p <= $m_value) {
             $p_value += $p;
        }
    }
    my $eps = 0.00000001; # Добавляем ко всем переменным, чтобы избежать деления на ноль
    my $se = sqrt(1/($a+$eps) + 1/($b+$eps) + 1/($c+$eps) + 1/($d+$eps)); # Станд. откл.
    my $odds_ratio = ($a+$eps)*($d+$eps)/(($b+$eps)*($c+$eps)); # Отношение шансов
    my $ci95_low = exp(log($odds_ratio) - 1.96*$se); # Нижняя граница 95% доверительного интервала
    my $ci95_high = exp(log($odds_ratio) + 1.96*$se); # Верхняя граница 95% доверительного интервала
    return $p_value, $odds_ratio, $ci95_low, $ci95_high, $l_value, $r_value, $m_value,;
}

Comments

[al-pas.livejournal.com]

Спасибо, но моя задача, к сожалению, состоит не в том, чтобы заставить библиотечную функцию R работать так же, как мой код, а в том, чтобы написать код, работающий так же, как библиотечная функция R именно в таком варианте её вызова:

h <- fisher.test(matrix(c(a_11, a_12, a_21, a_22), nrow=2)
pvalues[i] <- round(h$p.value, 5)
oddratio[i] <- round(h$estimate, 5)

[ogn-slon.livejournal.com]

А-а-а, понял. Да, тогда совет по моей ссылке вряд ли релевантен поставленной задаче.