Точный тест Фишера. Python против R
Я понимаю, что для ответа на мой вопрос самым правильным было бы посмотреть исходники 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,;
}
Результаты работы этого теста совпадают с результатами работы стандартной функции из библиотеки 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,;
}