2007-10-19

多重検定を重層的に繰り返したとき

未証明→
この日の記事の問題点は
『統計量 $P_j=\prod_{i=1}^M p_{ij}$ を考えたこと』
独立事象について、生起確率を掛け合わせることと、P値(生起確率の累積)を掛け合わせることとは別物であるから。

以下の記述は、この点について問題があることに留意しつつ、積分その他については、メモとして使うので、残す。
この留意は以降、『複数の多テストスタディの累積』のシリーズの(６)、2007/10/26まで適用されるべきである。

想定はこう。
同一キットを用いて、複数のサンプルセットでスクリーニングをする。
それぞれのキットにはN個の仮説検定がある。
サンプルセット数をMとする。
今、N検定には、諸問題があって、均一分布からのそれとみなせないかもしれないが、それを補正するなどして、均一分布からのP値のセットであるとみなせるだけの処理が適切に行われたものとする。
同一のキットを用いているので、Mサンプルセットのそれぞれから、N検定のそれぞれについて、P値が出てくる。
今、すべてのサンプルセットについて、すべての検定で帰無仮説が仮定できるとする。また、第iサンプルセットにおける、第j検定のP値を $p_{ij}$ とする。
Mサンプルセットについて独立であるなら、 ${p_{1i},p_{2i},...,p_{Mi}$ は、M次元立方体中に均一に分布するはずである。この分布が均一でないとすると、それはMサンプルセットが相互に独立とみなせないことになる(ただし、個々のサンプルセット内でのN個のp値は0-1に均一に分布しているものとみなしていることに注意する)。
その上で、第j検定について、Mサンプルセットについて通算した、統計量 $P_j=\prod_{i=1}^M p_{ij}$ を考える。
この統計量が $q$ 以下になるM次元立方体中の部分体積は１次元のときに
$V_1(q)=q$
２次元のときに
$V_2(q)=q\times (1- log(q))$
３次元のときに
$V_3(q)=q\times (1- log(q) + \frac{1}{2}(log(q))^2)$
４次元のときに
$V_4(q)=q\times (1- log(q) + \frac{1}{2}(log(q))^2-\frac{1}{6}(log(q))^3)$
証明していないが、k次元のときに
$V_k(q)=q\times \sum_{j}^{k} (-1)^{k-1}\frac{1}{(k-1)!}(log(q))^{k-1}$
$V_k(q)=q\times \sum_{j}^{k} \frac{1}{(k-1)!}(-log(q))^{k-1}$

この式は未証明だが、４次元までの移り行きを見た限り、多分大丈夫・・・・。
証明は漸化式を解けばよいはず。考えるにあたっては、 ${1,1,1,...,1}$ なる頂点の側の部分体積を、全体の体積から差し引くこととする、という点と、１次元増やすときには、その増えた次元については $q-1$ の範囲で積分すること、また、その増えた次元のスライスの体積は $q$ と増やした次元の変数とで与えられる、１次元低い次元の体積によって表されることに注意すると、 $V_{k+1}(q)=1-\int_q^1 1-V_{k}(\frac{q}{x}) dx$ が与えられる。これを解く。

これが正しいことは、Rで

rand1<-runif(10000)
rand2<-runif(10000)
rand3<-runif(10000)
rand4<-runif(10000)
multr<-rand1*rand2*rand3*rand4
V4<-multr*(1-log(multr)+1/2*(log(multr))^2-1/6*(log(multr))^3)
plot(sort(V4),seq(from=1/(length(V4)+1),to=length(V4)/(length(V4)+1),by=1/(length(V4)+1)))

により、掲載図のような直線になることから確認できる。
もしこうならなければ、それはセット間で「低いP値を出しがちなテストとそうでないテストの分布に偏りがある」ことを示している(はず)。
javaでは

package multSet;

import java.io.IOException;
import java.io.BufferedWriter;
import java.io.FileWriter;


public class MultSet {

	public int N;//No. tests per set
	public int M;//No. sets
	public double[][] p;//p ncol=M,nrow=N
	public double[] multP;
	public double[] expMultP;
	/**
	 * @param args
	 */
	public static void main(String[] args) throws IOException{
		// TODO 自動生成されたメソッド・スタブ
		int N=500000;
		int M=10;
		double[][] p= new double[N][M];
		BufferedWriter[] bw = new BufferedWriter[1];
		String out = "thomas10000-4_3.txt";
		bw[0] = new BufferedWriter(new FileWriter(out));
		for(int i=0;i<p.length;i++){
			for(int j=0;j<p[i].length;j++){
				p[i][j]=Math.random();
			}
		}
		MultSet ms = new MultSet(p);
		ms.PrintThomas(bw[0]);
		

	}
	
	public MultSet(double[][] p){
		N=p.length;
		M=p[0].length;
		multP=new double[N];
		expMultP=new double[N];
		for(int i=0;i<N;i++){
			for(int j=0;j<M;j++){
				multP[i]+=Math.log(p[i][j]);
			}
			expMultP[i]=ExpMultPLog(multP[i],M);
		}
		
	}
	
	public static double ExpMultP(double p,int M){
		double ret=0;
		double sl[] = StatUtilsX.Fisher.serialLogfact(M);
		for(int i=0;i<M;i++){
			ret += Math.pow((-1),(i-1))/Math.exp(sl[i-1])*Math.pow(Math.log(p),(i-1));
		}
		ret *=p;
		return ret;
	}
	public static double ExpMultPLog(double logp,int M){
		double ret=0;
		double sl[] = StatUtilsX.Fisher.serialLogfact(M);
		for(int i=1;i<=M;i++){
			double tmp =Math.log(-logp);
			//System.out.println(tmp);
			//ret += Math.pow((-1),(i-1))/Math.exp(sl[i-1])*Math.pow(logp,(i-1));
			ret += 1/Math.exp(sl[i-1])*Math.exp(tmp*(i-1));
			
			//ret += Math.pow((-1),(i-1))/Math.exp(sl[i-1])*Math.pow(-1, i-1)*Math.exp(tmp*(i-1));
		}
		ret *=Math.exp(logp);
		return ret;
	}
	public void PrintThomas(BufferedWriter bw)throws IOException{
		double[] exp= StatUtilsX.MiscUtil.DeepCopyDouble1(expMultP);
		StatUtilsX.MiscUtil.quickSortDouble(exp, 0, exp.length-1);
		for(int i=0;i<exp.length;i++){
			double tmp = (double)(i+1)/(double)(exp.length+1);
			BioBankPerm2.Run.out3File(bw, exp[i]+"\t"+tmp+"\n");
		}
		bw.close();
	}

}