Точный тест Фишера, основы реализации на Python

Таблицы сопряженности 2 х 2

Представим себе исследование влияния какого-либо воздействия на явление. Часть экспериментов выполняется в условиях заданного воздействия, часть — вне его, т.е. в обычных условиях. При этом регистрируется факт обнаружения явления как при наличии воздействия, так и без него. Задача состоит в том, чтобы выяснить, можно ли считать, что воздействие содействует явлению, или наоборот, воздействие, позволяет избежать явления. Такие задачи встречаются повсеместно — это и оценка влияния лечения на течение болезней, исследование процедур обработки семян, направленных на повышение их всхожести и многое многое другое.

Если привести задачу с проращиванием семян, то в качестве воздействия можно принять, например, факт обработки семян раствором марганцовки, а в роли явления — результат — проросло семя или нет в течение, допустим, 10 дней. Если эксперимент проводится над N — общим числом семян, то его результаты можно записать в виде таблицы сопряженности 2 x 2:

- П Н -
М a b a+b
В c d c+d
- a+c b+d N

 В таблице слева приняты следующие сокращения: П, Н — факт прорастания (П) или непрорастания (Н) семени; М, В — факт обработки марганцовкой (М) или обработки водой (В). a,b,c,d —  наблюдаемые количества взаимных сочетаний воздействия и явления, N = a+b+c+d.

Очевидно, если доля проросших семян в строке, соответствующей обработке марганцовкой М, будет значительно больше непроросших, и более того, при обработке водой, эти доли, например, будут практически равны, то можно говорить, что обработка марганцовкой, скорее всего влияет на прорастание семян. В этой попытке решить задачу мы использовали слишком много слов, допускающих неточности, это прежде всего "значительно больше", "практически равны".  Попробуем формализовать данные понятия, опираясь на комбинаторные представления при формировании таблиц сопряженности. Поставим такой вопрос: какова вероятность наблюдать конкретную таблицу сопряженности (a, b, c, d), если размещение в ячейках осуществляется случайно; иными словами, если N общее число семян (в некоторой корзине), которые далее размещаются по ячейкам в таблице сопряженности. Какова вероятность при этом наблюдать ситуацию, что ячейки будут заполнены по количеству семян в порядке a, b, c, d?  Эта вероятность определяется гипергеометрическим распределением:
$$ P(a,b; c,d) = \dfrac{C_{a+c}^a\cdot C_{b+d}^b}{C_N^{a+b}}$$

Формула гипергеометрического распределения определяется отношением возможных вариантов выбора a+b элементов из N (в данном случае — подлежащих обработке марганцовкой), и далее — последующим выбором "П" — проросших "элементов" (семян) и "Н" — непроросших соответственно. Поскольку выбор проросших и непроросших осуществляется независимо, общее число возможных вариантов — определяется произведением (числитель формулы гипергеометрического распределения).

Вероятность, задаваемая формулой гипергеометрического распределения — это вероятность наблюдать заданную конфигурацию — таблицу сопряженности — в условиях случайности. Это можно использовать, чтобы построить статистический тест. Пусть из эксперимента мы наблюдаем некоторую конфигурацию; по формуле гипергеометрического распределения мы можем вычислить вероятность наблюдать такую конфигурацию; если эта вероятность оказалась очень мала, то наблюдаемая конфигурация — маловероятна в условиях случайного размещения, что, в свою очередь, свидетельствует в пользу какого-то правила или порядка в организации таблицы. Таким образом, малая вероятность наблюдать заданную конфигурацию в таблице сопряженности свидетельствует в пользу значимого влияния тестируемого воздействия на результаты эксперимента. 

В пакете научных и инженерных расчетов SciPy гипергеометрическое распределение определяется функцией:

pmf = st.hypergeom(M, n, N) # обозначения M, n, N -- взяты из официальной документации функции.

$$pmf(k) = P(k, M, n, N) = \dfrac{C_n^k\cdot C_{M-n}^{N-k}}{C_M^N}$$

В терминах обозначений таблицы сопряженности: M = a+b+c+d; k=a, n = a+c, N = a + b.

Пример вычислений.

Исследуется вопрос об эффективности обработки с целью последующего проращивания жёлудей. В результате эксперимента построена следующая таблица сопряженности:

- не взошел взошел -
обработано 1 10 11
не обработано 4 3 7
- 5 13 18

 Нужно вычислить вероятность реализации  таблицы [1, 10; 4, 3], а также более "худшего" ("худшего" — по отношению к нулевой гипотезе, которая предполагает, что обработка не оказывает влияния) варианта, [0,11; 5,2], т.е. когда после обработки вообще все семена взошли. Если сумма этих вероятностей будет мала, то понятно, что обработка (а не случайность) определяет исход прорастания. Вычислить данные вероятности можно используя либо функцию scipy.stats.hypergeom, либо непосредственно выражения для \(C_n^k\).

В данном конкретном случае, вычисления приводят к следующему результату: $$ \begin{array}{l}\dfrac{C^1_{5}\cdot C^{10}_{13}}{C^{11}_{18}}+\dfrac{C^0_5C^{11}_{13}}{C^{11}_{18}} = \\\dfrac{1430}{31824} + \dfrac{78}{31824}\approx 0.047\end{array}$$

Таким образом, вероятность наблюдать текущую таблицу сопряженности, а также ту, которая в большей мере подтверждает эффективность обработки желудей, равна 0.047. Что весьма мало. В условиях случайности, проращивая много раз по 18 желудей, мы бы приблизительно в 5 случаях из 100 получали данную конфигурацию (считая 11 желудей условно обработанными, хотя реально — не подвергавшиеся никакой обработке). Следовательно, если мы вдруг наблюдаем данную конфигурацию, то скорее всего, обработка оказывает значимое влияние на всхожесть. Ошибка такого суждения, - или вероятность отвергнуть нулевую гипотезу (нулевой здесь является гипотеза о том, что обработка не влияет на всхожесть), когда она верна,  или - ошибка первого рода - равна 0.047 (прибл.). Это и есть так называемое p-value теста, вероятность наблюдать в условиях нулевой гипотезы текущую конфигурацию, или еще более — в данном случае — выраженную в эффективности действия обработки.

from math import factorial

def Cab(k, n):
    return factorial(n)/factorial(n-k)/factorial(k)


def P(a,b,c,d):
    return Cab(a, a + c) * Cab(b, b + d) / Cab(a + b, a + b + c + d)



if __name__ == '__main__':
    print('P(1,10; 4,3) + P(0,11;5,2) = ', P(1,10,4,3) + P(0,11,5,2))
    # P(1,10;4,3) + P(0,11;5,2)= 0.04738562091503268

Единственный вопрос с тестом Фишера — определение наиболее худших ситуаций по сравнению с наблюдаемой в таблице сопряженности. Если речь идет о проверке влияния действия какой-либо обработки на результаты эксперимента, то рассмотренный подход — вполне приемлем, он позволяет оценить вероятность наблюдать текущую конфигурацию, а также те, в которых это влияние еще более выражено. В случаях, когда понятие "наихудших" сценариев затруднено, прозводят суммирование всех вероятностей других возможных конфигураций, которые по значению меньше вероятности наблюдаемой конфигурации.

blog comments powered by Disqus