Fisher's Exact p value

   Fisher's exact p-value는 다음과 같은 상황에 적합한 값이다. 전체 인구가 10만 5000명인 도시에서 지역 A 에서 무작위로 100명을 골랐는데 그 중 37명이 자전거를 타고 출퇴근을 한다고 했다. 도시 전체 인구 중 자전거로 출퇴근하는 인구가 5000명으로 집계되었을 때, 지역 A는 다른 곳보다 자전거 통근자가 많다고 할 수 있을까? 즉,

전체 인구 10만 5000명
전체 자전거 통근자 수 5,000명
지역 A 에서 무작위로 선택한 사람 수 100명
무작위로 선택된 사람 중 자전거 통근자 37명

인 조건에서 자전거 통근자 37명이 과연 통계적으로 유의미한 값일까? 간단히 생각하면 전체 인구 중 자전거 통근자의 비율이 5,000 / 105,000 = 0.0476 이므로 이 확률이므로, A 지역이 특별하지 않다면 무작위로 선택한 100 명 중 이 비율로 자전거 통근자가 나와야 하므로, 100 x  0.0476 ~ 5 명 정도 나와야 할 것이다. 그런데 37명이 나왔으므로 A 지역은 자전거 통근자의 비율이 매우 높다고 할 수 있다. 그러나 37명에 대한 통계적 유의성을 위와 같이 approximation 하는 것보다는 다음과 같이 정확히 계산할 수 있다.

   개체의 전체 수 A 개 중 C 개를 선택했을 때,  관심있는 특성을 갖는 M 개의 집합 K에서 x 개가 선택되고 나머지는 K 이외에서 선택될 확률은 초기하분포를 따르며 식으로는 다음과 같다.

이 상황에서, K에서 S 개가 선택되었다면, 우연히 K에서 S 보다 많이 선택될 확률을 다음과 같이 계산하여,

이 값이 낮으면 낮을수록 S 개가 선택된 상황은 유의미하다고 할 수 있다. 즉, 우리는 확률적으로 일어나기 힘든 일에 대하여 '유의미하다'라고 말을 하고, 우연히 S개보다 많이 선택될 확률이 매우 낮다면 S 개가 선택된 것은 우연히 일어난 것이라 말하기 힘들고, 따라서 유의미하다고 할 수 있는 것이다.


두 개의 집합이 존재하고, 공통된 개수가 존재할 때, 그 공통된 개수에 대한 통계적 유의미함을 의미하는 것이 Fisher's exact p-value 인 것이다. 이와 같은 값은 위에서 자전거 통근자로 예로 든 것과 같이 binomial distribution 으로 추정할 수도 있고, 값이 커지면 물론 정규분포로 근사화 할 수 있긴 하다.

---

만약 Fisher's p-value를 계산하고자 하는데 A 가 매우 크다면 이항계수값이 커지기 때문에 C++의 built-in type으로는 계산이 되지 않는다. 이럴 땐 GNU Multiprecision (GMP) 라이브러리를 사용하여 계산할 수 있다. GMP는 컴퓨터의 저장공간이 허용하는만큼의 정확도를 갖으며 수를 표현할 수 있도록 해주는 라이브러리로, 2의 1억승, 10000!과 같은 값을 1의 자릿수까지 정확히 계산해 내는 라이브러리이다. 리눅스에서는 컴파일 할 때 -lgmp -lgmpxx option 을 붙여 주면 된다. (역시나 메뉴얼에 자세히 나와 있다). 실제 코드를 보면 다음과 같다.

#include <gmp.h>

#include <gmpxx.h>

 

// Returns

// [target]C[sel] * [total - target]C[sample - sel] / [total]C[sample]

// That is, probability of choosing [sample] numbers out of

// [total] where [sel] is chosen from [target].

mpf_class get_hypergeo_prob(int total, int sample, int target, int sel)

{

        mpz_class ztotal(total);

        mpz_class ztarget(target);

        mpz_class diff = ztotal - ztarget;

        mpz_t denominator; mpz_init(denominator);

 

        mpz_bin_ui(denominator, ztotal.get_mpz_t(), sample);

 

        mpz_t dividend1; mpz_init(dividend1);

        mpz_bin_ui(dividend1, ztarget.get_mpz_t(), sel);

 

        mpz_t dividend2; mpz_init(dividend2);

        mpz_bin_ui(dividend2, diff.get_mpz_t(), sample - sel);

 

        mpz_t dividend; mpz_init(dividend);

        mpz_mul(dividend, dividend1, dividend2);

 

        mpf_t fdividend; mpf_init(fdividend);

        mpf_set_z(fdividend, dividend);

        mpf_t fdenominator; mpf_init(fdenominator);

        mpf_set_z(fdenominator, denominator);

 

        mpf_t fresult; mpf_init(fresult);

        mpf_div(fresult, fdividend, fdenominator);

 

        mpf_class result(fresult);

 

        mpz_clear(denominator);

        mpz_clear(dividend1);

        mpz_clear(dividend2);

        mpz_clear(dividend);

 

        mpf_clear(fdividend);

        mpf_clear(fdenominator);

        mpf_clear(fresult);

 

 

        return result;

}

double get_upper_prob(int total, int sample, int target, int sel)

{

        mpf_class prob_product = 0;

        int i = 0;

        for(i = 0; i<sel; i++){

                prob_product += get_hypergeo_prob(total, sample, target, i);

        }

 

        return 1 - prob_product.get_d();

}

// 사용예

double p = get_upper_prob(all_gene_on_chip_num, sample_size, module_size, com_size);