正確確率検定では、P=1となる確率が高い



Fisherの正確確率検定においては、観測された分割表の周辺度数を満足するという条件の下に、とりうるすべての分割表について、その生起確率を求める。この確率が、観測された分割表のそれ以下であるような分割表について、生起確率を足し合わせたものが、P値である。

このことは、言い換えると、ある観測データ分割表があったときに、その分割表を得る確率が、その周辺度数の条件を満たす分割表の中で、もっとも高いときに、P=1となることを意味する。

分割表の数は有限であるので、P=1となるような分割表を得る確率はある値をとる。しかも、その確率とは、とりうる分割表の中で最大であるから、かなり大きめの値となる。さらに、そのような分割表が複数存在することもあるので、ある周辺度数が与えられたときに、かなりの確率でP=1となることがわかる。

2x2分割表から、その分割表がP=1になるかどうかの判断は、その分割表の生起確率が「最大」であるかどうかの判断さえできればつく。その判断は、k実は、生起確率を直接求めなくとも、|ad-bc|が最小になっているかどうかの判断をするだけでよい。ただし、ここで、a,b,c,dは2x2分割表のセルの値を意味するものとする。2x2分割表においては、その中の1つの変数のみを変化させることができるから、|ad-bc|aの関数となっている(f(a)=|ad-bc|。分割表が最大生起確率を与えるかどうかは、|ad-bc|aのときと、a’=a+1,a’’=a-1のときの3通りを比べることで判定可能である。簡単な式変形によりf(a’=a+1)=f(a)+N,f(a’’=a-1)=f(a)-Nであることが分かる(N=a+b+c+d)。この関係示したのが、こちらのエクセルである。

2x2分割表の周辺度数を与えて、そのP=1確率を求めるJavaソースが以下。また、ある生起確率pについて、総サンプル数、ケース数とから、P=1となる確率を算出するソースでもある。

ケース100人、コントロール100人でSNPのアレル頻度比較を行ったときに、Fisherの正確確率検定P値が1となる確率をアレル頻度の関数で表示したのが掲載図(オリジナルファイルはこちら)

関連記事はこちらも。


/*
* 作成日: 2006/09/29
*
* TODO この生成されたファイルのテンプレートを変更するには次へジャンプ:
* ウィンドウ - 設定 - Java - コード・スタイル - コード・テンプレート
*/
package StatUtils;

import DiscreteMathTools.MatrixExec;

/**
* @author ryamada
*
* TODO この生成された型コメントのテンプレートを変更するには次へジャンプ:
* ウィンドウ - 設定 - Java - コード・スタイル - コード・テンプレート
*/
public class Fisher {
/*
* このフィッシャーは両側検定を原則とする
* 片側検定は、どちらかの仮説を事前に知っている、というよりも
* 累積確率がより小さい側について片側検定している
*/
public static void main(String[] args)throws Exception {
/*
* 2x2 分割表の周辺度数(総サンプル数、ケース数、因子保有者数)を与える
* Fisher's exact P=1となる確率を求める
*/
int totalN=21;//総サンプル数
int caseN=10;//ケース数
int factorN=20;//ファクター保有
boolean printMarginal = true;//Prob1Marginalの計算情報を標準出力するかどうか
double prob=Prob1Marginal(totalN,caseN,factorN,printMarginal);
System.out.println("prob(P=1)= "+prob);
System.out.println("\n=========\n");


/*
* p,Ntotal,Ncaseを与え、そのようなスタディにおいてP=1を得る確率を算出する
*
* pとNtotalを与え、観測ファクターごとの確率を返す
* それはpとNtotal=Ncase+Ncontの関数である
* さらに、Ncaseを与え、そのような観測ファクター数のときにP=1
* となるような観測がなされる確率を計算し
* 観測ファクターについて合算する
*/
double p = 0.5;
int Ncase = 100;
int Ncont = 100;
int Ntotal = Ncase+Ncont;
boolean printProbP1=true;
double probPeq1 = ProbP1FromPandN(p,Ntotal,Ncase,printProbP1);
System.out.println("Probability of P=1(p="+p+" Ntotal="+Ntotal+" Ncase="+Ncase+" Ncont="+Ncont+")=\t"+probPeq1);

/*
* Ntotal,Ncaseにて定められるスタディデザインについて
* 複数のpにおけるP=1確率を返す
*/
int Ncase2=100;
int Ncont2=100;
int Ntotal2=Ncase2+Ncont2;
boolean b2=false;//個々のpについて計算過程を出力するか

int kizamiP=100;
double[] ps = new double[kizamiP+1];
for(int i=0;i<ps.length;i++){
ps[i]=i/(double)kizamiP;
}
double[] multiPeq1 = ProbP1FromPandNList(ps,Ntotal2,Ncase2,b2);
for(int i=0;i<multiPeq1.length;i++){
System.out.println(""+ps[i]+"\t"+multiPeq1[i]);
}
}
/*
* ProbP1FromPandNを複数のpについて返す
*/
public static double[] ProbP1FromPandNList(double[] ps, int Ntotal,int Ncase,boolean b2){
double[] ret = new double[ps.length];
for(int i=0;i<ret.length;i++){
ret[i]=ProbP1FromPandN(ps[i],Ntotal,Ncase,b2);
}
return ret;
}
/*
* allele freq,Ntotal,Ncaseを与え、P=1となるような分割表データの
* 得られる確率を返す
*/
public static double ProbP1FromPandN(double p,int Ntotal,int Ncase,boolean b){
double[][] probPNN = ProbContingency(p,Ntotal);
double probPeq1=0;
if(b){
System.out.println("No.factor(+)\t"+"Probablity to observe No.factor\t"+"Probability of P=1 for No.factor\t"
+"Product of both\t");
}
//System.out.println("Ntotal\tNcase\t"+ Ntotal + "\t"+Ncase);
for(int i=0;i<probPNN.length;i++){
double tmpprobPeq = Prob1Marginal(Ntotal,Ncase,i,false);
double tmpprobPeqtimesprob = tmpprobPeq*probPNN[i][1];
if(b){
System.out.println(""+ i + "\t"+probPNN[i][1]+"\t"+tmpprobPeq+"\t"+tmpprobPeqtimesprob);

}
probPeq1+=tmpprobPeqtimesprob;
}
//System.out.println("probP=1\t" + probPeq1);
return probPeq1;
}
/*
* allele freq fについて、Ntotalを与え、
* aごとの観測確率を返す
*/
public static double[][] ProbContingency(double p,int ntotal){
double[][] ret= new double[ntotal+1][2];//a,c,prob
int counter=0;
double numer = logfact(ntotal);
for(int i=0;i<ntotal+1;i++){
//for(int j=0;j<ncont;j++){
ret[counter][0]=i;
//ret[counter][1]=j;
//int[] tmp = {i,ncase-i,j,ncont-j};
double denom1 = logfact(i);
double denom2 = logfact(ntotal-i);
double partP = i*Math.log(p)+(ntotal-i)*Math.log(1-p);

ret[counter][1]=numer+partP-denom1-denom2;
ret[counter][1] = Math.exp(ret[counter][1]);
counter++;
//}
}

return ret;
}
/*
* 2x2分割表データから、P=1(Fisher's exact P)になる確率
*/
public static double Prob1(int[] d){
double p=0;
int sumA = d[0]+d[2];
int ncase = d[0]+d[1];
int ncont = d[2]+d[3];
int nsum = ncase+ncont;

//double expP=(double)sumA/(double)nsum;
//double expcaseA = expP*ncase;
int expcaseAbase= d[0];
int expcaseAplus=expcaseAbase+1;
int expcaseAminus=expcaseAbase-1;
int[] tableBase=new int[4];
int[] tablePlus=new int[4];
int[] tableMinus=new int[4];
tableBase[0]=expcaseAbase;
tableBase[1]=ncase-tableBase[0];
tableBase[2]=sumA-tableBase[0];
tableBase[3]=ncont-tableBase[2];

tablePlus[0]=expcaseAplus;
tablePlus[1]=ncase-tablePlus[0];
tablePlus[2]=sumA-tablePlus[0];
tablePlus[3]=ncont-tablePlus[2];

tableMinus[0]=expcaseAminus;
tableMinus[1]=ncase-tableMinus[0];
tableMinus[2]=sumA-tableMinus[0];
tableMinus[3]=ncont-tableMinus[2];
double probBase = P(tableBase);
double probPlus = P(tablePlus);
double probMinus = P(tableMinus);
String st ="2x2(a):\t";
for(int i=0;i<tableBase.length;i++){
st +=tableBase[i]+"\t";
}
st += "\n2x2(a+1):\t";
for(int i=0;i<tableBase.length;i++){
st +=tablePlus[i]+"\t";
}
st += "\n2x2(a-1):\t";
for(int i=0;i<tableMinus.length;i++){
st +=tableMinus[i]+"\t";
}
System.out.println(st);
System.out.println("probBase="+probBase);
System.out.println("probPlus="+probPlus);
System.out.println("probMinus="+probMinus);
double[] prob={probMinus,probBase,probPlus};
boolean[] max = {true,true,true};//Minus,Base,Plus
boolean trueBase=true;;
if(probBase<probPlus || probBase<probMinus){
trueBase=false;
}
if(probPlus<probBase || probPlus<probMinus){
max[2]=false;
}
if(probBase<probPlus || probBase<probMinus){
max[1]=false;
}
if(probMinus<probPlus || probMinus<probBase){
max[0]=false;
}
if(trueBase){
for(int i=0;i<max.length;i++){
//System.out.println("i " + max[i]);
if(max[i]){
p+=prob[i];
}
}
}

return p;
}
/*
* 2x2分割表周辺度数から、P=1(Fisher's exact P)になる確率
*/

public static double Prob1Marginal(int totalN,int caseN,int factorN,boolean b){
double p=0;
int sumA = factorN;
int ncase = caseN;
int ncont = totalN-ncase;
int nsum = ncase+ncont;

double expP=(double)sumA/(double)nsum;
double expcaseA = expP*ncase;
int expcaseAbase= (int)expcaseA;
int expcaseAplus=expcaseAbase+1;
int expcaseAminus=expcaseAbase-1;
int[] tableBase=new int[4];
int[] tablePlus=new int[4];
int[] tableMinus=new int[4];
double[] absADBC=new double[3];
tableBase[0]=expcaseAbase;
tableBase[1]=ncase-tableBase[0];
tableBase[2]=sumA-tableBase[0];
tableBase[3]=ncont-tableBase[2];
absADBC[1]=Math.abs(tableBase[0]*tableBase[3]-tableBase[1]*tableBase[2]);

tablePlus[0]=expcaseAplus;
tablePlus[1]=ncase-tablePlus[0];
tablePlus[2]=sumA-tablePlus[0];
tablePlus[3]=ncont-tablePlus[2];
absADBC[2]=Math.abs(tablePlus[0]*tablePlus[3]-tablePlus[1]*tablePlus[2]);

tableMinus[0]=expcaseAminus;
tableMinus[1]=ncase-tableMinus[0];
tableMinus[2]=sumA-tableMinus[0];
tableMinus[3]=ncont-tableMinus[2];
absADBC[0]=Math.abs(tableMinus[0]*tableMinus[3]-tableMinus[1]*tableMinus[2]);

double probBase = P(tableBase);
double probPlus = P(tablePlus);
double probMinus = P(tableMinus);
String st ="