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

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

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

同一の因子について、複数のテストを行い、それを統合するにあたって、Mantel-Haenszelがある(関連記事はこちら)。いくつかの区別できるサブグループに分けられるときなどに用いるし、メタアナリシス(関連記事はこちら)に用いることもある。
基本的なところに戻って、「どうして、複数のテスト(セット)を統合するときに何かの工夫(Mantel-Haenszelとか)をしないといけないのか」というと、個別のデータを単純に足し合わせて、一つのテストを行ったり、個別のデータから出てきた統計量を加減乗除したりすると、その結果として得られるみかけのP値の分布が、0−1に均一に分布しないためである。

単純な例を挙げる。今、2x2分割表が得られるような帰無仮説が成り立つ条件でテストを繰り返しているとする。Nテストを2セット行うものとする。おのおののセットについて、N個の2x2分割表検定が行われたとき、p値は0−1に渡って、均一に分布する。
今、Nテストのそれぞれについて、2セットに得られた2つの2x2分割表を足し合わせて、合算2x2分割表を作ったとしよう。この合算2x2分割表について分割表検定を行って得られるN個のp値は、どんな分布をとるだろうか?

  1. 2つのセットのケース・コントロールの人数が同一で、しかも、テスト対象の因子の母集団での率も同一であるとき、この合算のpは0−1に渡って均一に分布する。
  2. 2つのセットのケース・コントロールの総数が異なり、かつ、ケース対コントロールの比率が異なっても、テスト対象の因子の母集団での率が同一であるとき、この合算のpは0−1に渡って均一に分布する。
  3. 逆にテスト対象の因子の母集団での率が2セット間で異なるとき、ケース・コントロールの比率が2セットで同一であれば、この合算のpは0−1に渡って均一に分布する。
  4. 2セットの母集団における因子の率が異なり、かつ、2セットのケース・コントロールの比率が異なるとき、この合算のpは0−1に渡って均一に分布せず、小さい側に偏る。

このことからわかるように、複数のテスト(セット)を統合するにあたっては、「2セットの母集団における因子の率が異なり、かつ、2セットのケース・コントロールの比率が異なるようなときに発生する、小p値の増加を補正する」必要が出る。
1の場合については、Mantel-Haenszelテストと2セット単純合算テストとの結果は(ほぼ)同一。
2の場合については、Mantel-Haenszelテストと2セット単純合算テストとの結果は、ある程度の相関はあるものの、ばらつきが大きい。両者のp値の大小に一定の傾向もない
3の場合については、Mantel-Haenszelテストと2セット単純合算テストとの結果は、ある程度の相関があり、それは2の場合よりも強い。またMantel-Haenszelのp値が単純合算のそれより小さいほうへとずれる
4の場合については、Mantel-Haenszelと単純合算との間に相関はほぼない。

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

public class Dersimonian {
	int[][][][] count;
	public double CIval;
	double[] a;
	double[] or;
	double[] m1;
	double[] m2;
	double[] m3;
	double[] m4;
	double[] n;
	double[] chi;
	double[] p;
	double multp=0;
	double cp1=0;
	double cp2=0;
	double mh1 = 0;//MantelHaenszel stat
	//double mh1_2=0;
	double[][] joint = new double[2][2];
	double chij;
	double pjoint;
	double orjoint;
	double[] orjointCI;
	double pmh;
	
	/*
	 * Dersimonian
	 */
	public int direction;
	double[] logOR;
	double[] varOR;
	double[] logs;
	double[] vars;
	double[] vwts;
	double sumVwts;
	double[] vwtsSq;
	double sumVwtsSq;
	double[] vwtsXlogs;
	double sumVwtsXlogs;
	double logpooled;
	double[] logMinusLogpooled;
	double[] heterog;
	double sumHeterog;
	double tau2;
	double[] wts;
	double sumWts;
	double[] wtsXlogs;
	double sumXtsXlogs;
	double logDSL;
	double varDSL;
	double summaryTable;
	double DLchi;
	double DLp1;
	double DLp2;
	public double DLor;
	double[] DLorCI;
	double seLogDSL;
	
	
	
	/**
	 * @param args
	 */
	public Dersimonian(int[][][] table){
		/*
		 * tableはサブクラスの長さx2x2テーブル
		 */
		CIval=0.95;
		direction=1;
		count = new int[1][][][];
		count[0]=StatUtilsX.MiscUtil.DeepCopyInt3(table);
		//int df1=1;//自由度
		int df1 = 1;
		org.apache.commons.math.distribution.ChiSquaredDistributionImpl chidf1 = 
			new org.apache.commons.math.distribution.ChiSquaredDistributionImpl(df1);
		//int df2=2;//自由度
		int df2=count[0].length-1;
		org.apache.commons.math.distribution.ChiSquaredDistributionImpl chidf2 = 
			new org.apache.commons.math.distribution.ChiSquaredDistributionImpl(df2);
		
		int i=0;
		a=new double[count[0].length];
		or=new double[count[0].length];
		m1=new double[count[0].length];
		m2=new double[count[0].length];
		m3=new double[count[0].length];
		m4=new double[count[0].length];
		n=new double[count[0].length];
		chi=new double[count[0].length];
		p=new double[count[0].length];
		orjointCI=new double[2];
		boolean appropriate=true;
		/*
		 * DL
		 */
		logOR = new double[count[0].length];
		varOR = new double[count[0].length];
		logs = new double[count[0].length];//ORのDLはlogsとしてlogORをvarsとしてvarORを用いる
		vars = new double[count[0].length];
		vwts = new double[count[0].length];
		sumVwts=0;
		vwtsSq = new double[count[0].length];
		sumVwtsSq=0;
		vwtsXlogs = new double[count[0].length];
		sumVwtsXlogs=0;
		logpooled = 0;
		logMinusLogpooled = new double[count[0].length];
		heterog = new double[count[0].length];
		sumHeterog=0;
		tau2=0;
		wts = new double[count[0].length];
		sumWts=0;
		wtsXlogs = new double[count[0].length];
		sumXtsXlogs=0;
		logDSL=0;
		varDSL=0;
		summaryTable=0;
		DLp1=0;
		DLp2=0;
		DLor=0;
		DLorCI = new double[2];
		seLogDSL=0;
		try{
			for(int j=0;j<count[i].length;j++){
				//Pop 1
				
				a[j] = (count[i][j][0][0]*count[i][j][1][1]-count[i][j][0][1]*count[i][j][1][0]);
				or[j] = (double)(count[i][j][0][0]*count[i][j][1][1])/(double)(count[i][j][0][1]*count[i][j][1][0]);
				logOR[j]=Math.log(or[j]);
				varOR[j]=1/(double)count[i][j][0][0]+
				1/(double)count[i][j][0][1]+
				1/(double)count[i][j][1][0]+
				1/(double)count[i][j][1][1];
				//System.out.println("logOR "+logOR[j]);
				//System.out.println("varOR "+varOR[j]);
				logs[j]=logOR[j];
				vars[j]=varOR[j];
				vwts[j]=1/varOR[j];
				//System.out.println("vwts "+vwts[j]);
				sumVwts+=vwts[j];
				//System.out.println("sumVwts++ "+ sumVwts);
				vwtsSq[j]=vwts[j]*vwts[j];
				sumVwtsSq+=vwtsSq[j];
				vwtsXlogs[j]=vwts[j]*logs[j];
				sumVwtsXlogs+=vwtsXlogs[j];
				
				m1[j] = count[i][j][0][0]+count[i][j][0][1];
				m2[j] = count[i][j][1][0]+count[i][j][1][1];
				m3[j] = count[i][j][0][0]+count[i][j][1][0];
				m4[j] = count[i][j][0][1]+count[i][j][1][1];
				n[j] = m3[j]+m4[j];
				double[][] d = {{count[i][j][0][0],count[i][j][0][1]},{count[i][j][1][0],count[i][j][1][1]}};
				chi[j] = StatUtilsX.Chi.simple2x2Chi(d);
				p[j]=1-chidf1.cumulativeProbability(chi[j]);
				joint[0][0]+=count[i][j][0][0];
				joint[0][1]+=count[i][j][0][1];
				joint[1][0]+=count[i][j][1][0];
				joint[1][1]+=count[i][j][1][1];
				if(count[i][j][0][0]==0 ||
						count[i][j][0][1]==0 ||
						count[i][j][1][0]==0 ||
						count[i][j][1][1]==0){
					appropriate=false;
				}
			}
			if(appropriate){
//				System.out.println("sumVwts "+sumVwts);
				//System.out.println("sumVwtsSq "+ sumVwtsSq);
				logpooled=sumVwtsXlogs/sumVwts;
				//System.out.println("logpooled\t"+logpooled);
				for(int j=0;j<count[i].length;j++){
					logMinusLogpooled[j]=logs[j]-logpooled;
					heterog[j]=logMinusLogpooled[j]*logMinusLogpooled[j]/vars[j];
					//System.out.println("heterog "+heterog[j]);
					sumHeterog+=heterog[j];
				}
				//System.out.println("sumHeterog "+sumHeterog);
				
				//System.out.println("tmp "+tmp);
				double tmp2 = (sumHeterog-(double)(count[0].length-1))/(sumVwts-sumVwtsSq/sumVwts);
				//System.out.println("tmp2 "+tmp2);
				tau2=Math.max(0, tmp2);
				for(int j=0;j<count[i].length;j++){
					wts[j]=1/(vars[j]+tau2);
					sumWts+=wts[j];
					wtsXlogs[j]=wts[j]*logs[j];
					sumXtsXlogs+=wtsXlogs[j];
				}
				logDSL=sumXtsXlogs/sumWts;
				varDSL=1/sumWts;
				seLogDSL=Math.sqrt(varDSL);
				summaryTable=logDSL/Math.sqrt(varDSL);
				DLchi=summaryTable*summaryTable;
				DLp1=1-chidf1.cumulativeProbability(DLchi);
				DLp2=1-chidf2.cumulativeProbability(sumHeterog);
				DLor=Math.exp(logDSL);
				DLorCI[0]=Math.exp(logDSL-1.96*seLogDSL);
				DLorCI[1]=Math.exp(logDSL+1.96*seLogDSL);
				
				/*
				
				orjoint=joint[0][0]*joint[1][1]/(joint[0][1]*joint[1][0]);
				if(a[0]*a[1]<0){
					multp=Math.min(p[0]/p[1], p[1]/p[0]);
				}else{
					multp=p[0]*p[1];
				}
				cp1=multp*(1-Math.log(multp));
				cp2=multp*(1-Math.log(multp)/2);
				double mhbunsi=0;
				double mhbunbo=0;
				double jointorbunsi=0;
				double jointorbunbo=0;
				for(int j=0;j<count[0].length;j++){
					mhbunsi+=(a[j]/n[j]);
					mhbunbo+=(m1[j]*m2[j]*m3[j]*m4[j])/(n[j]*n[j]*(n[j]-1));
		
					jointorbunsi+=(count[i][j][0][0]*count[i][j][1][1])/n[j];
					jointorbunbo+=(count[i][j][1][0]*count[i][j][0][1])/n[j];
					
				}
				mhbunsi=Math.pow(mhbunsi,2);
				mh1=mhbunsi/mhbunbo;
				orjoint = jointorbunsi/jointorbunbo;
				
				//mh1 = Math.pow((a[0]/n[0]+a[1]/n[1]),2)/
				//(m1[0]*m2[0]*m3[0]*m4[0]/(n[0]*n[0]*(n[0]-1))+m1[1]*m2[1]*m3[1]*m4[1]/(n[1]*n[1]*(n[1]-1)));
				
				//System.out.println("mh "+ mh1+" "+mh1_2);
				//double[][] joint ={{count[i][0][0][0]+count[i][1][0][0],count[i][0][0][1]+count[i][1][0][1]},
						//{count[i][0][1][0]+count[i][1][1][0],count[i][0][1][1]+count[i][1][1][1]}};
				chij=StatUtilsX.Chi.simple2x2Chi(joint);
				pjoint = 1-chidf1.cumulativeProbability(chij);
				orjointCI[0]=Math.pow(orjoint,(1-1.96/Math.pow(mh1,0.5)));
				orjointCI[1]=Math.pow(orjoint,(1+1.96/Math.pow(mh1,0.5)));
				
				pmh = 1-chidf1.cumulativeProbability(mh1);
				//BioBankPerm2.Run.out3File(bw[0], "i=\t"+i+"\t"+chi[0]+"\t"+chi[1]+"\t"+chij+"\t"+p[0]+"\t"+p[1]+"\t"+pjoint+"\t"+mh1+"\t"+pmh+"\t"+multp+"\t"+cp1+"\t"+cp2+"\t");
				
				//System.out.print("i=\t"+i+"\t"+chi[0]+"\t"+chi[1]+"\t"+p[0]+"\t"+p[1]+"\t"+mh1+"\t"+pmh+"\t"+multp+"\t"+cp1+"\t"+cp2+"\t");
				/*
				for(int j=0;j<count[i].length;j++){
					for(int k=0;k<count[i][j].length;k++){
						for(int l=0;l<count[i][j][k].length;l++){
							//BioBankPerm2.Run.out3File(bw[0], count[i][j][k][l]+ "\t");
							System.out.print(count[i][j][k][l]+ "\t");
						}
					}
				}
				//BioBankPerm2.Run.out3File(bw[0], "\n");
				System.out.println();
				*/
			}
			
		}catch (org.apache.commons.math.MathException ex) {
			
		}
		
	}
	/**
	 * @param args
	 */
	public Dersimonian(int[][][] table,double CIval_){
		/*
		 * tableはサブクラスの長さx2x2テーブル
		 */
		CIval=CIval_;
		direction=1;
		count = new int[1][][][];
		count[0]=StatUtilsX.MiscUtil.DeepCopyInt3(table);
		//int df1=1;//自由度
		int df1 = 1;
		org.apache.commons.math.distribution.ChiSquaredDistributionImpl chidf1 = 
			new org.apache.commons.math.distribution.ChiSquaredDistributionImpl(df1);
		//int df2=2;//自由度
		int df2=count[0].length-1;
		org.apache.commons.math.distribution.ChiSquaredDistributionImpl chidf2 = 
			new org.apache.commons.math.distribution.ChiSquaredDistributionImpl(df2);
		org.apache.commons.math.distribution.NormalDistributionImpl nd = new org.apache.commons.math.distribution.NormalDistributionImpl();
		
		int i=0;
		a=new double[count[0].length];
		or=new double[count[0].length];
		m1=new double[count[0].length];
		m2=new double[count[0].length];
		m3=new double[count[0].length];
		m4=new double[count[0].length];
		n=new double[count[0].length];
		chi=new double[count[0].length];
		p=new double[count[0].length];
		orjointCI=new double[2];
		boolean appropriate=true;
		/*
		 * DL
		 */
		logOR = new double[count[0].length];
		varOR = new double[count[0].length];
		logs = new double[count[0].length];//ORのDLはlogsとしてlogORをvarsとしてvarORを用いる
		vars = new double[count[0].length];
		vwts = new double[count[0].length];
		sumVwts=0;
		vwtsSq = new double[count[0].length];
		sumVwtsSq=0;
		vwtsXlogs = new double[count[0].length];
		sumVwtsXlogs=0;
		logpooled = 0;
		logMinusLogpooled = new double[count[0].length];
		heterog = new double[count[0].length];
		sumHeterog=0;
		tau2=0;
		wts = new double[count[0].length];
		sumWts=0;
		wtsXlogs = new double[count[0].length];
		sumXtsXlogs=0;
		logDSL=0;
		varDSL=0;
		summaryTable=0;
		DLp1=0;
		DLp2=0;
		DLor=0;
		DLorCI = new double[2];
		seLogDSL=0;
		try{
			for(int j=0;j<count[i].length;j++){
				//Pop 1
				
				a[j] = (count[i][j][0][0]*count[i][j][1][1]-count[i][j][0][1]*count[i][j][1][0]);
				or[j] = (double)(count[i][j][0][0]*count[i][j][1][1])/(double)(count[i][j][0][1]*count[i][j][1][0]);
				logOR[j]=Math.log(or[j]);
				varOR[j]=1/(double)count[i][j][0][0]+
				1/(double)count[i][j][0][1]+
				1/(double)count[i][j][1][0]+
				1/(double)count[i][j][1][1];
				//System.out.println("logOR "+logOR[j]);
				//System.out.println("varOR "+varOR[j]);
				logs[j]=logOR[j];
				vars[j]=varOR[j];
				vwts[j]=1/varOR[j];
				//System.out.println("vwts "+vwts[j]);
				sumVwts+=vwts[j];
				//System.out.println("sumVwts++ "+ sumVwts);
				vwtsSq[j]=vwts[j]*vwts[j];
				sumVwtsSq+=vwtsSq[j];
				vwtsXlogs[j]=vwts[j]*logs[j];
				sumVwtsXlogs+=vwtsXlogs[j];
				
				m1[j] = count[i][j][0][0]+count[i][j][0][1];
				m2[j] = count[i][j][1][0]+count[i][j][1][1];
				m3[j] = count[i][j][0][0]+count[i][j][1][0];
				m4[j] = count[i][j][0][1]+count[i][j][1][1];
				n[j] = m3[j]+m4[j];
				double[][] d = {{count[i][j][0][0],count[i][j][0][1]},{count[i][j][1][0],count[i][j][1][1]}};
				chi[j] = StatUtilsX.Chi.simple2x2Chi(d);
				p[j]=1-chidf1.cumulativeProbability(chi[j]);
				joint[0][0]+=count[i][j][0][0];
				joint[0][1]+=count[i][j][0][1];
				joint[1][0]+=count[i][j][1][0];
				joint[1][1]+=count[i][j][1][1];
				if(count[i][j][0][0]==0 ||
						count[i][j][0][1]==0 ||
						count[i][j][1][0]==0 ||
						count[i][j][1][1]==0){
					appropriate=false;
				}
			}
			if(appropriate){
//				System.out.println("sumVwts "+sumVwts);
				//System.out.println("sumVwtsSq "+ sumVwtsSq);
				logpooled=sumVwtsXlogs/sumVwts;
				//System.out.println("logpooled\t"+logpooled);
				for(int j=0;j<count[i].length;j++){
					logMinusLogpooled[j]=logs[j]-logpooled;
					heterog[j]=logMinusLogpooled[j]*logMinusLogpooled[j]/vars[j];
					//System.out.println("heterog "+heterog[j]);
					sumHeterog+=heterog[j];
				}
				//System.out.println("sumHeterog "+sumHeterog);
				
				//System.out.println("tmp "+tmp);
				double tmp2 = (sumHeterog-(double)(count[0].length-1))/(sumVwts-sumVwtsSq/sumVwts);
				//System.out.println("tmp2 "+tmp2);
				tau2=Math.max(0, tmp2);
				for(int j=0;j<count[i].length;j++){
					wts[j]=1/(vars[j]+tau2);
					sumWts+=wts[j];
					wtsXlogs[j]=wts[j]*logs[j];
					sumXtsXlogs+=wtsXlogs[j];
				}
				logDSL=sumXtsXlogs/sumWts;
				varDSL=1/sumWts;
				seLogDSL=Math.sqrt(varDSL);
				summaryTable=logDSL/Math.sqrt(varDSL);
				DLchi=summaryTable*summaryTable;
				DLp1=1-chidf1.cumulativeProbability(DLchi);
				DLp2=1-chidf2.cumulativeProbability(sumHeterog);
				DLor=Math.exp(logDSL);
				double cival2=nd.inverseCumulativeProbability((1-CIval)/2);
				DLorCI[0]=Math.exp(logDSL+cival2*seLogDSL);
				DLorCI[1]=Math.exp(logDSL-cival2*seLogDSL);
				
				/*
				
				orjoint=joint[0][0]*joint[1][1]/(joint[0][1]*joint[1][0]);
				if(a[0]*a[1]<0){
					multp=Math.min(p[0]/p[1], p[1]/p[0]);
				}else{
					multp=p[0]*p[1];
				}
				cp1=multp*(1-Math.log(multp));
				cp2=multp*(1-Math.log(multp)/2);
				double mhbunsi=0;
				double mhbunbo=0;
				double jointorbunsi=0;
				double jointorbunbo=0;
				for(int j=0;j<count[0].length;j++){
					mhbunsi+=(a[j]/n[j]);
					mhbunbo+=(m1[j]*m2[j]*m3[j]*m4[j])/(n[j]*n[j]*(n[j]-1));

					jointorbunsi+=(count[i][j][0][0]*count[i][j][1][1])/n[j];
					jointorbunbo+=(count[i][j][1][0]*count[i][j][0][1])/n[j];
					
				}
				mhbunsi=Math.pow(mhbunsi,2);
				mh1=mhbunsi/mhbunbo;
				orjoint = jointorbunsi/jointorbunbo;
				
				//mh1 = Math.pow((a[0]/n[0]+a[1]/n[1]),2)/
				//(m1[0]*m2[0]*m3[0]*m4[0]/(n[0]*n[0]*(n[0]-1))+m1[1]*m2[1]*m3[1]*m4[1]/(n[1]*n[1]*(n[1]-1)));
				
				//System.out.println("mh "+ mh1+" "+mh1_2);
				//double[][] joint ={{count[i][0][0][0]+count[i][1][0][0],count[i][0][0][1]+count[i][1][0][1]},
						//{count[i][0][1][0]+count[i][1][1][0],count[i][0][1][1]+count[i][1][1][1]}};
				chij=StatUtilsX.Chi.simple2x2Chi(joint);
				pjoint = 1-chidf1.cumulativeProbability(chij);
				orjointCI[0]=Math.pow(orjoint,(1-1.96/Math.pow(mh1,0.5)));
				orjointCI[1]=Math.pow(orjoint,(1+1.96/Math.pow(mh1,0.5)));
				
				pmh = 1-chidf1.cumulativeProbability(mh1);
				//BioBankPerm2.Run.out3File(bw[0], "i=\t"+i+"\t"+chi[0]+"\t"+chi[1]+"\t"+chij+"\t"+p[0]+"\t"+p[1]+"\t"+pjoint+"\t"+mh1+"\t"+pmh+"\t"+multp+"\t"+cp1+"\t"+cp2+"\t");
				
				//System.out.print("i=\t"+i+"\t"+chi[0]+"\t"+chi[1]+"\t"+p[0]+"\t"+p[1]+"\t"+mh1+"\t"+pmh+"\t"+multp+"\t"+cp1+"\t"+cp2+"\t");
				/*
				for(int j=0;j<count[i].length;j++){
					for(int k=0;k<count[i][j].length;k++){
						for(int l=0;l<count[i][j][k].length;l++){
							//BioBankPerm2.Run.out3File(bw[0], count[i][j][k][l]+ "\t");
							System.out.print(count[i][j][k][l]+ "\t");
						}
					}
				}
				//BioBankPerm2.Run.out3File(bw[0], "\n");
				System.out.println();
				*/
			}
			
		}catch (org.apache.commons.math.MathException ex) {
			
		}
		
	}
	public String String(String sep1,String sep2){
		String ret ="";
		ret +=direction+sep1;
		ret += CIval+sep1;
		for(int i=0;i<count.length;i++){
			for(int j=0;j<count[i].length;j++){
				for(int k=0;k<count[i][j].length;k++){
					for(int l=0;l<count[i][j][k].length;l++){
						ret += count[i][j][k][l]+sep1;
					}
				}
			}
		}
		for(int i=0;i<joint.length;i++){
			for(int j=0;j<joint[i].length;j++){
				ret += joint[i][j]+sep1;
			}
		}
		for(int i=0;i<or.length;i++){
			ret += or[i]+sep1;
		}
		ret += orjoint+sep1;
		ret += orjointCI[0]+sep1+orjointCI[1]+sep1;
		for(int i=0;i<chi.length;i++){
			ret += chi[i]+sep1;
		}
		ret += chij+sep1;
		for(int i=0;i<p.length;i++){
			ret += p[i]+sep1;
		}
		ret += pjoint+sep1;
		ret += tau2+sep1;
		ret += DLchi+sep1;
		
		ret += DLp1+sep1;
		ret += sumHeterog+sep1;
		ret += DLp2+sep1;
		ret += DLor+sep1;
		ret += DLorCI[0]+sep1+DLorCI[1]+sep1;
		/*
		ret+=mh1+sep1;
		ret+=pmh+sep1;
		ret+=multp+sep1;
		ret+=cp1+sep1;
		ret+=cp2+sep1;
		*/
		return ret;
	}
	public String StringMin(String sep1,String sep2){
		String ret ="";
		ret +=direction+sep1;
		ret += CIval+sep1;
		/*
		for(int i=0;i<count.length;i++){
			for(int j=0;j<count[i].length;j++){
				for(int k=0;k<count[i][j].length;k++){
					for(int l=0;l<count[i][j][k].length;l++){
						ret += count[i][j][k][l]+sep1;
					}
				}
			}
		}
		for(int i=0;i<joint.length;i++){
			for(int j=0;j<joint[i].length;j++){
				ret += joint[i][j]+sep1;
			}
		}
		for(int i=0;i<or.length;i++){
			ret += or[i]+sep1;
		}
		ret += orjoint+sep1;
		ret += orjointCI[0]+sep1+orjointCI[1]+sep1;
		for(int i=0;i<chi.length;i++){
			ret += chi[i]+sep1;
		}
		ret += chij+sep1;
		for(int i=0;i<p.length;i++){
			ret += p[i]+sep1;
		}
		ret += pjoint+sep1;
		*/
		ret += tau2+sep1;
		ret += DLchi+sep1;
		
		ret += DLp1+sep1;
		ret += sumHeterog+sep1;
		ret += DLp2+sep1;
		ret += DLor+sep1;
		ret += DLorCI[0]+sep1+DLorCI[1]+sep1;
		/*
		ret+=mh1+sep1;
		ret+=pmh+sep1;
		ret+=multp+sep1;
		ret+=cp1+sep1;
		ret+=cp2+sep1;
		*/
		return ret;
	}
	public String Label(String sep1,String sep2){
		String ret="";
		ret+="direction"+sep1;
		ret+="CIval"+sep1;
		//System.out.println("## count len "+count.length);
		for(int i=0;i<count.length;i++){
			//System.out.println("## count[i] len "+count[i].length);
			for(int j=0;j<count[i].length;j++){
				//System.out.println("## count[i][j] len "+count[i][j].length);
				
				for(int k=0;k<count[i][j].length;k++){
					//System.out.println("## count[i][j][k] len "+count[i][j][k].length);
					
					for(int l=0;l<count[i][j][k].length;l++){
						
						//System.out.println("&&&&&");
						ret += "count"+i+","+j+","+k+","+l+sep1;
					}
				}
			}
		}
		for(int i=0;i<joint.length;i++){
			for(int j=0;j<joint[i].length;j++){
				ret += "joint"+i+","+j+sep1;
			}
		}
		for(int i=0;i<or.length;i++){
			ret += "or"+ i + sep1;
		}
		ret += "orjoint"+sep1;
		ret += "orjointCI1"+sep1+"orjointCI2"+sep1;
		for(int i=0;i<chi.length;i++){
			ret += "chi"+i+sep1;
		}
		ret += "chijoint"+sep1;
		for(int i=0;i<p.length;i++){
			ret += "p"+i+sep1;
		}
		ret += "pjoint"+sep1;
	
		
		ret += "tau"+sep1;
		ret += "DLchi"+sep1;
		ret += "DLp"+sep1;
		ret += "heterorg"+sep1;
		ret += "hetero p"+sep1;
		ret += "DLor"+sep1;
		ret +="DLorCI1"+sep1;
		ret +="DLorCI2"+sep2;
		
		/*
		ret+="mh1"+sep1;
		ret+="mhp"+sep1;
		ret+="multp"+sep1;
		ret+="cp1"+sep1;
		ret+="cp2"+sep1;		
		*/
		return ret;
	}
	public String LabelMin(String sep1,String sep2){
		String ret="";
		ret+="direction"+sep1;
		ret+="CIval"+sep2;
		//System.out.println("## count len "+count.length);
		/*
		for(int i=0;i<count.length;i++){
			//System.out.println("## count[i] len "+count[i].length);
			for(int j=0;j<count[i].length;j++){
				//System.out.println("## count[i][j] len "+count[i][j].length);
				
				for(int k=0;k<count[i][j].length;k++){
					//System.out.println("## count[i][j][k] len "+count[i][j][k].length);
					
					for(int l=0;l<count[i][j][k].length;l++){
						
						//System.out.println("&&&&&");
						ret += "count"+i+","+j+","+k+","+l+sep1;
					}
				}
			}
		}
		*/
		/*
		for(int i=0;i<joint.length;i++){
			for(int j=0;j<joint[i].length;j++){
				ret += "joint"+i+","+j+sep1;
			}
		}
		for(int i=0;i<or.length;i++){
			ret += "or"+ i + sep1;
		}
		ret += "orjoint"+sep1;
		ret += "orjointCI1"+sep1+"orjointCI2"+sep1;
		for(int i=0;i<chi.length;i++){
			ret += "chi"+i+sep1;
		}
		ret += "chijoint"+sep1;
		for(int i=0;i<p.length;i++){
			ret += "p"+i+sep1;
		}
		ret += "pjoint"+sep1;

		*/
		ret += "tau"+sep1;
		ret += "DLchi"+sep1;
		ret += "DLp"+sep1;
		ret += "heterorg"+sep1;
		ret += "hetero p"+sep1;
		ret += "DLor"+sep1;
		ret +="DLorCI1"+sep1;
		ret +="DLorCI2"+sep2;
		
		/*
		ret+="mh1"+sep1;
		ret+="mhp"+sep1;
		ret+="multp"+sep1;
		ret+="cp1"+sep1;
		ret+="cp2"+sep1;		
		*/
		return ret;
	}
	public static void main(String[] args) throws IOException{
		// TODO 自動生成されたメソッド・スタブ
		int[][] Nsample = {{200,100},{400,600}};
		int Nmarker=10000;
		BufferedWriter[] bw = new BufferedWriter[1];
		String file ="MAFsameUnbalance.txt";
		bw[0] = new BufferedWriter(new FileWriter(file));
		int x = 10000000;
		int seed = (int)(Math.random()*x);
		Utils.MersenneTwisterFast mz = new Utils.MersenneTwisterFast(seed);
		int df1=1;//自由度
		org.apache.commons.math.distribution.ChiSquaredDistributionImpl chidf1 = 
			new org.apache.commons.math.distribution.ChiSquaredDistributionImpl(df1);
		
		double[][] af = new double[Nmarker][Nsample.length];
	
		for(int i=0;i<af.length;i++){
			af[i][0]=mz.nextDouble()*0.8+0.1;
			af[i][1]=mz.nextDouble()*0.8+0.1;
			af[i][1]=af[i][0];
		}
		int[][][][] count = new int[Nmarker][Nsample.length][Nsample[0].length][2];
		for(int i=0;i<count.length;i++){
			for(int j=0;j<count[i].length;j++){
				for(int k=0;k<count[j].length;k++){
					for(int l=0;l<Nsample[j][k];l++){
						double tmp = mz.nextDouble();
						if(tmp<af[i][j]){
							count[i][j][k][0]++;
						}else{
							count[i][j][k][1]++;
						}
					}
				}
			}
			/*
			int k = mz.nextInt(5)+1;
			int k2 = mz.nextInt(5)+1;
			count[i][1][0][0]=k*count[i][0][0][0];
			count[i][1][0][1]=k*count[i][0][0][1];
			count[i][1][1][0]=k2*count[i][0][1][0];
			count[i][1][1][1]=k2*count[i][0][1][1];
			*/
		}
		
		BioBankPerm2.Run.out3File(bw[0], "i=\tID\tchi1\tchi2\tchisum\tp1\tp2\tpsum\tMHchi\tMHp\tmultp\tcp1\tcp2\n");
		System.out.println("i=\tID\tchi1\tchi2\tp1\tp2\tMHchi\tMHp\tmultp\tcp1\tcp2\t");
		try{
			for(int i=0;i<count.length;i++){
				double[] a = {0,0};
				double[] m1 = {0,0};
				double[] m2 = {0,0};
				double[] m3 = {0,0};
				double[] m4 = {0,0};
				double[] n = {0,0};
				double[] chi = {0,0};
				double[] p = {0,0};
				for(int j=0;j<count[i].length;j++){
					//Pop 1
					a[j] = (count[i][j][0][0]*count[i][j][1][1]-count[i][j][0][1]*count[i][j][1][0]);
					m1[j] = count[i][j][0][0]+count[i][j][0][1];
					m2[j] = count[i][j][1][0]+count[i][j][1][1];
					m3[j] = count[i][j][0][0]+count[i][j][1][0];
					m4[j] = count[i][j][0][1]+count[i][j][1][1];
					n[j] = m3[j]+m4[j];
					double[][] d = {{count[i][j][0][0],count[i][j][0][1]},{count[i][j][1][0],count[i][j][1][1]}};
					chi[j] = StatUtilsX.Chi.simple2x2Chi(d);
					p[j]=1-chidf1.cumulativeProbability(chi[j]);
					
					
				}
				double multp = 0;
				if(a[0]*a[1]<0){
					multp=Math.min(p[0]/p[1], p[1]/p[0]);
				}else{
					multp=p[0]*p[1];
				}
				double cp1=multp*(1-Math.log(multp));
				double cp2=multp*(1-Math.log(multp)/2);
				double mh1 = Math.pow((a[0]/n[0]+a[1]/n[1]),2)/
				(m1[0]*m2[0]*m3[0]*m4[0]/(n[0]*n[0]*(n[0]-1))+m1[1]*m2[1]*m3[1]*m4[1]/(n[1]*n[1]*(n[1]-1)));
				double[][] joint ={{count[i][0][0][0]+count[i][1][0][0],count[i][0][0][1]+count[i][1][0][1]},
						{count[i][0][1][0]+count[i][1][1][0],count[i][0][1][1]+count[i][1][1][1]}};
				double chij=StatUtilsX.Chi.simple2x2Chi(joint);
				double pjoint = 1-chidf1.cumulativeProbability(chij);
				double pmh = 1-chidf1.cumulativeProbability(mh1);
				BioBankPerm2.Run.out3File(bw[0], "i=\t"+i+"\t"+chi[0]+"\t"+chi[1]+"\t"+chij+"\t"+p[0]+"\t"+p[1]+"\t"+pjoint+"\t"+mh1+"\t"+pmh+"\t"+multp+"\t"+cp1+"\t"+cp2+"\t");
				
				System.out.print("i=\t"+i+"\t"+chi[0]+"\t"+chi[1]+"\t"+p[0]+"\t"+p[1]+"\t"+mh1+"\t"+pmh+"\t"+multp+"\t"+cp1+"\t"+cp2+"\t");
				for(int j=0;j<count[i].length;j++){
					for(int k=0;k<count[i][j].length;k++){
						for(int l=0;l<count[i][j][k].length;l++){
							BioBankPerm2.Run.out3File(bw[0], count[i][j][k][l]+ "\t");
							System.out.print(count[i][j][k][l]+ "\t");
						}
					}
				}
				BioBankPerm2.Run.out3File(bw[0], "\n");
				System.out.println();
			}
		}catch (org.apache.commons.math.MathException ex) {
			
		}
		bw[0].close();
	}

}