Точный тест Фишера. 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

[ogn-slon.livejournal.com]

На StackOverflow задавали аналогичный вопрос по немного другому (но сходному) стат. поводу. Вот как на этот вопрос там ответили: Basically whenever you're having trouble matching output between statistical software, start looking at the default arguments to see if you should be enabling or disabling these adjustements. То есть перед тем, как лезть глубоко в код стандартных библиотек, можно попробовать поискать в документации/комментариях, какими аргументами можно выключить корректирующие примочки данного конкретного софта чтобы использовать просто стандартный метод (типа, "correction=False" в примере по ссылке).

P.S. Сам я в теме не разбираюсь, но этот коммент все же рискну запостить. :)

[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]

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

[yazon.livejournal.com]

С прошедшим!