多重検定を重層的に繰り返したとき(４)
この日の記事の問題点は
『統計量を考えたこと』
独立事象について、生起確率を掛け合わせることと、P値(生起確率の累積)を掛け合わせることとは別物であるから。
以下の記述は、この点について問題があることに留意しつつ、積分その他については、メモとして使うので、残す。
この留意は以降、『複数の多テストスタディの累積』のシリーズの(６)、2007/10/26まで適用されるべきである。
同一の因子について、複数のテストを行い、それを統合するにあたって、Mantel-Haenszelがある(関連記事はこちら)。いくつかの区別できるサブグループに分けられるときなどに用いるし、メタアナリシス(関連記事はこちら)に用いることもある。
基本的なところに戻って、「どうして、複数のテスト(セット)を統合するときに何かの工夫(Mantel-Haenszelとか)をしないといけないのか」というと、個別のデータを単純に足し合わせて、一つのテストを行ったり、個別のデータから出てきた統計量を加減乗除したりすると、その結果として得られるみかけのP値の分布が、０−１に均一に分布しないためである。
単純な例を挙げる。今、２ｘ２分割表が得られるような帰無仮説が成り立つ条件でテストを繰り返しているとする。Nテストを２セット行うものとする。おのおののセットについて、N個の２ｘ２分割表検定が行われたとき、ｐ値は０−１に渡って、均一に分布する。
今、Nテストのそれぞれについて、２セットに得られた２つの２ｘ２分割表を足し合わせて、合算２ｘ２分割表を作ったとしよう。この合算２ｘ２分割表について分割表検定を行って得られるN個のｐ値は、どんな分布をとるだろうか？
- ２つのセットのケース・コントロールの人数が同一で、しかも、テスト対象の因子の母集団での率も同一であるとき、この合算のｐは０−１に渡って均一に分布する。
- ２つのセットのケース・コントロールの総数が異なり、かつ、ケース対コントロールの比率が異なっても、テスト対象の因子の母集団での率が同一であるとき、この合算のｐは０−１に渡って均一に分布する。
- 逆にテスト対象の因子の母集団での率が２セット間で異なるとき、ケース・コントロールの比率が２セットで同一であれば、この合算のｐは０−１に渡って均一に分布する。
- ２セットの母集団における因子の率が異なり、かつ、２セットのケース・コントロールの比率が異なるとき、この合算のｐは０−１に渡って均一に分布せず、小さい側に偏る。
このことからわかるように、複数のテスト(セット)を統合するにあたっては、「２セットの母集団における因子の率が異なり、かつ、２セットのケース・コントロールの比率が異なるようなときに発生する、小ｐ値の増加を補正する」必要が出る。
1の場合については、Mantel-Haenszelテストと２セット単純合算テストとの結果は(ほぼ)同一。
2の場合については、Mantel-Haenszelテストと２セット単純合算テストとの結果は、ある程度の相関はあるものの、ばらつきが大きい。両者のｐ値の大小に一定の傾向もない
3の場合については、Mantel-Haenszelテストと２セット単純合算テストとの結果は、ある程度の相関があり、それは2の場合よりも強い。またMantel-Haenszelのｐ値が単純合算のそれより小さいほうへとずれる
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はサブクラスの長さｘ２ｘ２テーブル
*/
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はサブクラスの長さｘ２ｘ２テーブル
*/
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();
}
}
