2006-08-13

False Discovery Rate

多数の仮説検定を行っている。ある棄却水準を与え、多数の仮説のどれが棄却されどれは棄却されないかを判定する方法のひとつ。

基本となる考え方はこう。

具体的に。

100個の仮説を検定しているとする。100個の仮説のそれぞれに、個別仮説検定P値を得る。独立な100個の仮説があり、それらが独立であるとすると、100個のP値の分布は、０から１までの一様分布をとると期待される。したがって、100個のP値が、100個の仮説検定という条件においても仮説を棄却するかどうかの判断は、仮説が１個の場合よりも厳しくすることが適当であると考える。

Bonferroni補正とFDRを比べるとき、その検出力は、複数の帰無仮説がすべて正しいときには変わらないが、複数の帰無仮説の中に棄却されるべき仮説が多く含まれれば含まれるほど、FDRを用いることによる検出力の改善効果が期待できると言う。この関係は、駆け足で読むシリーズMultivariate permutation testsにおける、combination fucntionのTippettのfunctionとLiptakのfunctionとの間柄に似たところとも言えよう(記事はこちら)(厳密な比較・定義はまだ。印象のみに基づくコメント)

どのくらい厳しくするかの方法がいくつかある。いくつかの方法はこちらの記事に記載があるとおりである。その中で一番簡単なのは、線形なBH法である。

今、全体の棄却水準を $￥alpha=0.05$ としたとき、個別P値として大きいほうからx番目のP値の棄却水準として、 $￥alpha*(100-x+1)/100)$ とする。たとえば、100個の仮説検定において、1番大きなｐ値をとった仮説については、0.05。2番目に大きなｐ値をとった仮説については、0.05-0.05/100=0.0495,...,１番小さなｐ値をとった仮説については0.05-0.05*99/100 = 0.0005。この厳しくした棄却水準との比較を、個別ｐ値の大きい方から順に行い、初めて棄却することとなった仮説を含め、その仮説以下の個別P値を持つ仮説をそろって棄却する。

簡単にその具合を図示した。紺のラインは、個別仮説ｐ値を降順にプロットしたもの。縦軸は、個別仮説ｐ値とBH法補正の棄却水準。右下のプロットは縦軸が対数目盛り。横軸は、100個の仮説。ピンクのラインは、BH法で厳しくした $￥alpha=0.05$ の値。紺ラインがピンクラインより下になった初めての仮説から右側の仮説が棄却されたとする。

Adaptive BH法は、BH法のうち、棄却水準の厳しくするときの補正の仕方を若干緩めたものである。BH法では、補正にあたり、総仮説数と個別検定P値の序列によって線形厳しくした(補正棄却水準は総仮説数で割られている)が、Adaptive BHでは、これを総仮説数ではなく、ある条件を満たした、仮説数に減らしてやることで、補正棄却水準を若干緩めている。

今、全仮説が正しく、それぞれの仮説に得られるp値はランダムに出ているとすると、その値の分布は０から

１にわたって、一様な分布をとる。棄却されるべき仮説が含まれている場合には、その分布は、０がちになるはずである。棄却されるべき仮説群は、小さいｐ値を持つとともに、その値は、０−１一様分布よりも０側に偏った分布となっている。他方、真の仮説群のｐ値は、０−１一様分布になるはずで、今、仮説をその個別検定p値の順に並べたときには、どこかしらにその移行点を示す部分があるはずである。それを示すのに有用と見なせる変数を用いて、『ある条件を満たした、仮説数』を算出する。



public class FdrBH {
	public static void main(String[] args) {

		

		double[] test = {0.5,0.2,0.4,0.1,0.01,0.02,0.6,0.8,0.6,0.0025};

		int[] fdrret = FdrBHtest(test,0.050001);

		System.out.println("rejected list");

		for(int i=0;i<test.length;i++){

			if(fdrret[i]==1){

				System.out.println(i + " p "  + + test[i]);

			}

		}

		

		

	}

	

	public static int[] FdrBHtest(double[] p, double a){

		//for(int i=0;i<p.length;i++){

			//System.out.println("p " + p[i]);

		//}

		int[] rank = MiscUtil.bublSortRetInt(p);

		//for(int i=0;i<p.length;i++){

			//System.out.println("r " + rank[i]);

		//}

		int counter =0;

		int ans =-1;

		

		for(int i=0;i<p.length;i++){

			double thres = (a*(double)(rank[i]+1))/(double)(p.length);

			if(p[i]<thres){

				if(ans<rank[i]){

					ans = rank[i];

				}

			}

		}

		int[] ret = new int[p.length];

		for(int i=0;i<ret.length;i++){

			if(rank[i]<=ans){

				ret[i]=1;

			}else{

				ret[i]=0;

			}

			//System.out.println("p rank judge " + p[i] + " " + rank[i] + " " + ret[i]);

		}

		

		return ret;

		

	}

}