第７章 ノンパラメトリック組み合わせ解析の例 7.2 Multivariate paired Observation
- 7.2 Multivariate paired Observation
- ２群に属するn個体についてq種類の変数をk時刻に観測する、というようなデータ構造
- 各時刻において、２群の値に大小の傾向がないことを検定
- 教科書に記載されている個別時刻のP値が少しずれているように思えるので少々不安だが。。。このソースで、６章までと整合性がとれていることから、擬似乱数列のぶれによる値の乖離なのかもしれない。ちなみに、個別仮説１，２，３に対するPは0.15605 0.01785 0.01155 、Combining fucntiosとしてFisher,LiptakLogit,Tippettのぞれぞれにて
corrP Fisher 0.0114
corrPLiptakLogit 0.0114
corrP Tippett 0.0354
public class TensorDataVec {
public int rank;
public int[] dimension;
public double[][] datavec;
public double[] data;
public int[][] index;public static void main(String[] args) {
int rank_= 4;
int[] dim_ = {3,10,2,1};
//double [] data_ = {8.5,8.6,9.8,9.4,14.9,14.1,6.1,6.2,7.7,7.3,8.0,7.4,9.8,9.6,11.5,11.0,13.0,12.9,
//14.5,14.1,16.7,15.9,19.0,18.6,6.1,6.3,6.4,6.5,7.7,7.9,5.1,4.9,8.7,7.6,9.1,8.7,6.4,6.7,9.2,8.9,13.3,12.6,
//7.9,7.6,7.7,8.0,10.2,10.5,5.1,5.1,6.3,6.5,9.6,9.2,10.7,10.2,12.8,12.1,14.7,13.8};
double[] data_ = {8.5,8.6,6.1,6.2,9.8,9.6,14.5,14.1,6.1,6.3,5.1,4.9,6.4,6.7,7.9,7.6,5.1,5.1,10.7,10.2,
9.8,9.4,7.7,7.3,11.5,11.0,16.7,15.9,6.4,6.5,8.7,7.6,9.2,8.9,7.7,8.0,6.3,6.5,12.8,12.1,
14.9,14.1,8.0,7.4,13.0,12.9,19.0,18.6,7.7,7.9,9.1,8.7,13.3,12.6,10.2,10.5,9.6,9.2,14.7,13.8};
TensorDataVec t = new TensorDataVec(rank_,dim_,data_);
//Create Tensor Object and print out its features
//System.out.println("Tensor object is consisted of int rank, int[] dimension, double[] data" +
//" and int[][] index");
//PrintTensor(t,"\t","\n");
/*
* 個々の仮説は、第１次元のそれぞれについて、第３次元をパーミュテーションしても
*/
double DEV= 0.00000001;
int perm =10000;
//異なるdelta値での統計量を算出
//double[] delta = {-1000,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,1000};
double[] delta = {0};
double[][][] permT = new double[t.dimension[0]][2][perm];
String mval = "";
String mval2 ="";
int c = 1000000;
int seed = (int)(Math.random()*c);
//int seed =100;
for(int m=0;m<delta.length;m++){
System.out.println("\n\n=================\nd value_" + m + "\n=================\n");
double[][] x = new double[t.dimension[0]][t.dimension[1]];
// int seed = 10000;
MersenneTwisterFast mt = new MersenneTwisterFast(seed);
MersenneTwisterFast mt2 = new MersenneTwisterFast(seed);
for(int i=0;i<x.length;i++){
for(int j=0;j<x[i].length;j++){
int[] a = {i,j,0,0};
int[] b = {i,j,1,0};
double aElem = ReturnElem(t,a);
double bElem = ReturnElem(t,b);
//System.out.println("aElem bElem delta " + aElem + " " + bElem + " " + delta[m]);
x[i][j] = aElem-bElem+delta[m];
//System.out.println("x " + x[i][j]);
}
}
double[][] to = new double[x.length][2];
//System.out.println("to.len " + to.length);
for(int i=0;i<to.length;i++){
to[i][0] = signPerm.sum(x[i]);
to[i][1] = signPerm.K(x[i]);
}
//複数の統計量を比較できるように
permT = signPerm.SumAndKAndWilcoxonSignPermVector(x,perm,mt);
//permutation P(個別仮説)の算出
double[] Pvalue = new double[t.dimension[0]];
for(int i=0;i<Pvalue.length;i++){
Pvalue[i] = signPerm.PermTestP(to[i][0],permT[i][0]);
}
//permutation P(個別仮説)をオリジナルデータと全パーミュテーション試行に算出
double[][] PAllPerm = new double[t.dimension[0]][perm];
for(int i=0;i<PAllPerm.length;i++){
PAllPerm[i] = signPerm.PermTestPPerm(permT[i][0]);
}
//PermTestP(to[0],permT[0]);
String StpermP = "PermP\t";
for(int i=0;i<Pvalue.length;i++){
StpermP += Pvalue[i] + "\t";
}
System.out.println(StpermP);
// T''算出
//Direct Combining Function
double tppDirect =0;
for(int i=0;i<to.length;i++){
tppDirect+=to[i][0];
}
double[] tppDirectPerm = new double[perm];
for(int i=0;i<perm;i++){
tppDirectPerm[i]=0;
for(int j=0;j<permT.length;j++){
tppDirectPerm[i]+=permT[j][0][i];
}
}
double corrPDirect = signPerm.PermTestP(tppDirect,tppDirectPerm);
//Fisher's omnibus combining function
double tppFisher = CombFxFisher(Pvalue);
double[] tppFisherPerm = new double[perm];
// LiptakLogit's omnibus combining function
double tppLiptakLogit = CombFxLiptakLogit(Pvalue);
double[] tppLiptakLogitPerm = new double[perm];
// Tippett's omnibus combining function
double tppTippett = CombFxTippett(Pvalue);
double[] tppTippettPerm = new double[perm];
for(int i=0;i<tppFisherPerm.length;i++){
double[] tmp = new double[PAllPerm.length];
for(int j=0;j<tmp.length;j++){
tmp[j]=PAllPerm[j][i];
}
tppFisherPerm[i] = CombFxFisher(tmp);
tppLiptakLogitPerm[i] = CombFxLiptakLogit(tmp);
tppTippettPerm[i] = CombFxTippett(tmp);
}
double corrPFisher = signPerm.PermTestP(tppFisher,tppFisherPerm);
System.out.println("corrP Fisher " + corrPFisher);
double corrPLiptakLogit = signPerm.PermTestP(tppLiptakLogit,tppLiptakLogitPerm);
System.out.println("corrPLiptakLogit " + corrPLiptakLogit);
double corrPTippett = signPerm.PermTestP(tppTippett,tppTippettPerm);
System.out.println("corrP Tippett " + corrPTippett);
System.out.println("corrP Direct " + corrPDirect);/*
for(int i=0;i<permT.length;i++){
// 度数分布・信頼区間など
int[] h=signPerm.hist(permT[i][0],20);
double[][] hdouble = signPerm.histAccum(permT[i][0],1000);
// 統計量の信頼区間を算出表示
double confA = 0.000001;
double[] confInterval = signPerm.confInterval(hdouble,confA);
System.out.println("confident interval\t" + confInterval[0] + "\t" + confInterval[1]);
// 度数分布を標準出力
// System.out.println(histSt);
System.out.println("Histogram\n
第６章 ノンパラメトリック組み合わせ方法
- 補助記事はこちらとこちら
- 6.1 イントロダクション
- 6.2 組み合わせ方、Combining frunctions
- あるデータセットがある。そのデータセットから、q個の変数それぞれについて、帰無仮説を棄却するか否かのp値
が得られる。q個の変数のすべてについてのglobal 帰無仮説に関するパーミュテーション統計量
は
の関数であり(
)、Combininb functions と呼ばれる。これらは次の条件を満たす
について、
:q個の変数のすべてにおいて、変数を単独で評価したときに、あるデータセットがより得られにくいときには、
- q個の変数のいずれか１つでもその帰無仮説のp値が無限に小さくなるとき、
- 0より大なるあらゆる棄却水準について前項で示した上限値未満の
- その他、必須条件ではないが、次のような特徴は有用である。
- 変数ごとのp値を交換しても、global統計量は変わらない
- このような
- ２段階の検定プロセス
- あるデータセットがある。q個の変数について、統計量が得られる(１つずつ、q種類)。
- ついで、パーミュテーションデータセットにつき、同様に１つずつ、q種類の統計量が得られる。
- 観測データセットについてq個の棄却確率が得られる
- また、パーミュテーションデータセットのそれぞれについて、q個の変数のそれぞれについてそのデータセットで帰無仮説が棄却される確率も得られる。これは、全部でqｘパーミュテーション回数である。
- 組み合わせ(複合)統計量は、このq個の棄却確率の関数として定義でき、観測データセットについて、その値が得られるとともに、パーミュテーションのそれぞれについても、その値が得られる。
- 組み合わせ統計量は後述のとおり、複数のそれが知られているが、そのそれぞれについて、観測データに対する値と、パーミュテーションのそれぞれに対する値とが得られている。global棄却確率は、この組み合わせ統計量が、分布上どのくらい偏っているかで決まる。
- 組み合わせ統計量が複数知られ、それぞれに、特徴があるので、使い分けたり、複数のそれを用いて総合判定したりする。
- Combining functionsの例
- Fisher omnibus combining function
- Liptak combining function
- Liptak function with logistic transformation
- Tippett combining function
- Lancaster combining function
- Others
- Direct combining functions
- 複数の個別Tが同質であるときは、思い切って、それらの和をとって、その和について、観測データでの値がパーミュテーション分布での位置をｐとする方法もある。この方法が許される条件であると、圧倒的に計算量が減る
- Direct combining functions
- Fisher omnibus combining function
- あるデータセットがある。そのデータセットから、q個の変数それぞれについて、帰無仮説を棄却するか否かのp値
public static double CombFxFisher(double[] a){
double ret = 0;
for(int i=0;i<a.length;i++){
ret += -2 * Math.log(a[i]);
}
return ret;
}
public static double CombFxLiptakLogit(double[] a){
double ret=0;
for(int i=0;i<a.length;i++){
ret += Math.log((1-a[i])/a[i]); }
return ret;
}
public static double CombFxTippett(double[] a){
double ret=0;
double min =a[0];
for(int i=1;i<a.length;i++){
if(min>a[i]){
min = a[i];
}
}
ret =1-min;
return ret;
}
- 変数の数が複数なだけでなく、その変数について検定したい内容が複数化したり、変数の組み合わせ方を特定のしかたで定めたりする場合の取り扱いについての記載が少しある
- 6.3,6.4 Unbiasedness, consistency, asymptotic properties
- Unbiasedness と asymptotic properties については単変量での記述の拡張
- Consistency
- 個々の変数の
から
- 前述のcombining functionsはいずれも、グラフで描くと、globalな棄却水準はグラフが区切る空間のうち原点を含む側の面積(体積)で表される。このとき、１つの変数についての棄却水準が０に近づくにつれ、グラフはそれ以外の変数については、0-1のどこにあっても、棄却域に入るような形をしている。この事情はglobalな棄却水準を無限小にしても成り立つこれは、変数のうち１つでも棄却されれば、その他の変数での棄却の可能性によらず、global な仮説全体も棄却されるからである
- このようなグラフは凹なグラフとなっている
- combining fucntionのグラフが凸なとき、上の条件を満たすことはできない。たとえば
1は、凸グラフ(凹でない）ので、上の条件を満たさない
- このように、上の条件を満たすことをconsistentと呼び、combining functionとして、multivariate permutation testsで用いることができるための条件となっている
- 個々の変数の
- 6.5 本章の結語
- 解析データにうまく用いることのできるcombining fucntion(s)がある場合は、それを用いることを考慮するとして、そのようなcombining function(s)がないときには、変数の数の次元で
- 解析データにうまく用いることのできるcombining fucntion(s)がある場合は、それを用いることを考慮するとして、そのようなcombining function(s)がないときには、変数の数の次元で
*1：Hotelling's T2 Wikipediaの記事
第５章 多標本パーミュテーションテストの理論
- 5.1 基本事項の復習
- パーミュテーションテストでは、全パーミュテーションを調べる場合と、モンテカルロで分布を作る場合とがある。前者では、観測データの正確なパーセンタイルがでる。後者では推定分布がでる。棄却水準を考えるとき、前者では、棄却水準を与える統計量以上の統計量が観測データから得られている場合には、棄却し、そうでないときは棄却しない。他方、後者では、観測データ統計量が棄却水準統計量と等しい場合には、棄却するかしないかは確率的に決まる。
- 5.2 多標本パーミュテーションで検定するいろいろ(第４章であつかったものを列挙、整理)
- ２群の平均値を比較
』もしくは
- One-way ANOVA
- Goodness Of Fit (分割表)
- グループ数
で序列のあるカテゴリカル変数の場合
、ただし、
は、グループjにおけるカテゴリ1-i累積人数の割合。
、
- グループ数
- グループ数
- Two-way ANOVAのフレーム
- ２群の平均値を比較
- 5.3 推定量の不偏性と検定力
- 5.4 漸近性
- 5.5 パーミュテーションの中心極限定理
第４章（３）Repeated Observations サンプルは対等・multiple量的変数 多標本１変数のパーミュテーションテストの例
- 4.5
- 複数の標本(n)があり、それらの違いはないものとできるとき、それぞれの標本に複数の観測(k)がなされて、n x k 行列状のデータが得られている。これについて、ｋ変数は偏りがないことを示したいとする。具体的には、ある測定量について経時データをとりつつ、測定時が影響していないことを示したいときなど
- 順列は、
である
- Friedman's Rankテストが統計量として知られている。
- 各サンプルのデータをサンプル内で順序づけし、nサンプルにわたって、各計測時において平均値をとり
とする
- 各サンプルのデータをサンプル内で順序づけし、nサンプルにわたって、各計測時において平均値をとり
- 対応するPermutation 統計量
を次の式で与える
- ただし、
は全データの平均、
は、それぞれ、サンプルごとの平均、変数ごとの平均。
- ただし、
- これらの一致はこの程度(掲載図、およびこちらの図)
- ソースはこちら
public class RepeatedData {public static void main(String[] args) {
//2D data
double[][] data ={
{320,278,236,222,232},
{478,513,415,359,292},
{921,701,645,526,458},
{213,230,261,253,199},
{273,338,323,332,222},
{392,302,289,305,172},
{469,443,292,235,233},
{422,389,359,331,185},
{613,649,626,588,636},
{395,318,298,269,328},
{462,400,360,247,284}
};
int numUnit = data.length;
int numVar = data[0].length;//No. categories
// 異なるdelta値での統計量を算出
int numseg = 1;
double[] delta = new double[numseg];
double initial=0;
double dif = -initial*2/(double)numseg;
for(int i=0;i<delta.length;i++){
delta[i] = initial+dif*i;
}
//double[] delta = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
//double[] delta = {0};
//for(int i=0;i<delta.length;i++){
//delta[i] /= 0.1;
//}
/*
int[] categoryNum = new int[data.length];
for(int i=0;i<categoryNum.length;i++){
categoryNum[i] = data[i].length;
}
double[] concat = concatMult(data);
*/
int c = 100000000;
int seed = (int)(Math.random()*c);
//int seed = 10000;
String mval ="";
String mval2 ="";
for(int w=0;w<delta.length;w++){
MersenneTwisterFast mt = new MersenneTwisterFast(seed);
double to = GoFTR(data);
double to2 = Friedman(data);
/*
double[] toTd = GoFTypT(data,group);
double to = toTd[0];
double to2 = toTd[1];
double to3 = GoFT(data,group);
double to4 = GoFT2(data, group);
double t2ADo = GoFT2AD(data,group);
double t2fADo = GoFT2fAD(data,group);
double tchio = Chisq(data,group);
*/
double DEV= 0.00000001;
int perm =100;
System.out.println("to to2 " + to + " " + to2);
double[] t = new double[perm];
double[] t2 = new double[perm];
//double[] t2AD = new double[perm];
//double[] t2fAD = new double[perm];
//double[] tchi = new double[perm];
for(int i=0;i<perm;i++){
for(int j=0;j<data.length;j++){
shuffleDoubleMZ(data[j],mt);
}
//shuffleDoubleMZ(data[0],mt);
t[i] = GoFTR(data);
t2[i] = Friedman(data);
/*
double[] tmpt = GoFTypT(data,group);
t[i] = tmpt[0];
t2[i] = tmpt[1];
double tmpt2 = GoFT(data,group);
double tmpt3 = GoFT2(data,group);
t2AD[i] = GoFT2AD(data,group);
t2fAD[i] = GoFT2fAD(data,group);
tchi[i] = Chisq(data,group);
*/
System.out.println("t\tmwt\t" + t[i] +"\t" + t2[i]);
}
// 度数分布・信頼区間など
int[] h=signPerm.hist(t,20);
double[][] hdouble = signPerm.histAccum(t,1000);
// 統計量の信頼区間を算出表示
double confA = 0.000001;
double[] confInterval = signPerm.confInterval(hdouble,confA);
System.out.println("confident interval\t" + confInterval[0] + "\t" + confInterval[1]);
// 度数分布を標準出力
// System.out.println(histSt);
System.out.println("Histogram\n
第４章（３）Goodness-of-fit (分割表検定) 多標本１変数順序のあるカテゴリカル変数のパーミュテーションテストの例
- 4.4.1
- 分割表は、そのままだと、パーミュテーション処理に向かないので、unit-by-unit(１サンプルずつ）が帰属カテゴリを持つようなデータにしたうえで、サンプルについてモンテカルロパーミュテーション処理をする
- 4.4.2 スコアづけによる検定
- カテゴリには順序があるので、個々のカテゴリに値を与える
- 量的変数データのときに、個々のサンプルのデータに順位をつける解析(Mann-Whitney)は、サンプルごとに異なるカテゴリを与え、そのカテゴリに差１の値を付与したこととみなせる。ただし、分割表検定では、同点サンプルが多数ある点が量的変数データのランクづけと異なる点である
- 統計量は
。これは、
とpermutationally eqivalent である (それを示すプロットはこちら)。それを算出するソースは末尾
- 4.4.3 典型的なGoodness-of-Fit テスト
- 前項で与えたカテゴリ値は、唐突な感じをぬぐえない。観測データに即して、カテゴリ値を調整してみた統計量が教科書に提示されており、さらにそのpermutational equivalenceが示されている。
- 前項のカテゴリ値の与え方で問題となるのは、次のような例を考えるとわかりやすい。
- ２群に３カテゴリがある。{5,10,5},{10,10,10}なる観測データだったとする。３カテゴリにそれぞれ、1,2,3の値を与えたとき、
- ２群に３カテゴリがある。{5,10,5},{10,10,10}なる観測データだったとする。３カテゴリにそれぞれ、1,2,3の値を与えたとき、
- 次の統計量は、観測度数から、カテゴリに与える値を調整している。その数式表現は次の通り
- ここで、
は、グループ２のk=1,2,...,iまでの累積観測度数。
は２グループについての累積観測度数の和。
は総サンプル数。
はグループ１のサンプル数となっている。
- この式だと、前項で、各カテゴリに値を与えた式と形が異なるのでわかりにくい。この式を変形すると、各カテゴリiには、カテゴリ数kについて
を与えたことになっている。
の項は、カテゴリ順行でのサンプルの増え方と逆行でのサンプルの増え方の両方を加味して、サンプル全体がカテゴリについて均一か否かに対応する調整項である。
の項は、グループ１とグループ２とのサンプル数の違いについての調整項である
- 前項でもpermutational equivalentとして簡易な統計量を示したが、こちらにも同様に簡易統計量が存在する
- カテゴリ値に単純な値1,2,3などを与えた場合と、観測データ準拠のカテゴリ値を与えた場合とでは、パーミュテーション統計量の大小に逆転はおきないが、単純なカテゴリ値の場合に同点例が多発してくることを図示したのが、掲載図および、こちらの図
- 4.4.5 カテゴリに序列がない場合と、グループが３以上の場合への拡張
- グループが２個でカテゴリに順序があるときには、カテゴリの値は、１つの尺度で比較可能であったので、その差のみを数量化すればよく、また、グループ同士の比較も、１次元の比較であるので、その差を数量化すればよかった
- カテゴリに順序がなくなると、おのおののカテゴリについて、数値化し、最終的にすべてのカテゴリで合算する必要がある。これは、カテゴリ数の次元を持つ空間について、個々のカテゴリが表す次元における値を定量し、それを合算することになり、これをカテゴリ数空間内の距離と考えてよい
- カテゴリ数が２であり、片側検定するときには、カテゴリに本来、大小関係になくても、ダミー変数を用いて１次元的に取り扱うことができる。カテゴリ数が３以上になるとそれができなくなる
- グループ数が３以上に増えることも、グループが作る次元が多次元となる点は、カテゴリの多次元化と同様である。
- それを反映して、カテゴリ数もしくはグループ数が３以上にも対応した統計量は、観測度数の自乗によって、表現されている。
- また、グループ数が３以上のときのカテゴリが順序がないときと、あるときとでは、カテゴリ別の観測度数をそのまま用いるか、累積観測度数を用いるかの違いがあることもわかる。カテゴリに順序があるとき、その順序情報は観測度数を累積すること・カテゴリに与える値に累積観測度数を用いた値を適用することで、順序情報を取り込んでいる。逆に、カテゴリに順序がないときは、累積しない観測度数を用いていることも式からわかる
- 順序のないカテゴリカル変数の分割表データに対する統計量としては、Pearson's chi square を含め、複数存在するが、順序のある場合には、適用できない。順序のあるカテゴリの分割表検定統計量については、こちらの記事も参考に。
- 順序のないカテゴリカル変数におけるカイ自乗統計量と以下で示す、
とは、相関はあるが、permutationally equivalent ではないことに注意。それを示すプロットはこちら。それをエクセルで確かめるためのファイルはこちら
- グループ数
- グループ数
- グループ数
第４章（２）One-way ANOVA 多標本１変数パーミュテーションテストの例
- C>=２群のすべてが等しいか、そうでないかの検定
- 全nデータが、
個ずつC群に分かれている。
- パーミュテーションは
通り
- 統計量としては、
がわかりやすいが、それとpermutationally にequivalentな統計量は
がある。また、one-way ANOVAで用いるFもT,Sと１対１の単調な関係にある。他方、ランクテストの統計量である、Kruskal-Wallisもあるが、これは、permutationally にequivalentではない。それを示したのが掲載図とこの図(リンク)である。
- これまでと同様に、教科書中のデータサンプルについて、実行し、掲載図を作るにあたって用いたJavaソースはこちら
public class oneWayANOVAunivarMultiSamplePerm {public static void main(String[] args) {
// 異なるdelta値での統計量を算出
int numseg = 1;
double[] delta = new double[numseg];
double initial=0;
double dif = -initial*2/(double)numseg;
for(int i=0;i<delta.length;i++){
delta[i] = initial+dif*i;
}
//double[] delta = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
//double[] delta = {0};
//for(int i=0;i<delta.length;i++){
//delta[i] /= 0.1;
//}
//Data
double[][] data = new double[3][];
double[] tmp1 = {10.2,8.2,8.9,8.0,8.3,8.0};
double[] tmp2 = {12.2,10.6,9.9,13.0,8.1,10.8,11.5};
double[] tmp3 = {9.2,10.5,9.2,8.7,9.0};
//double[] tmp1 = {3.42,3.84,3.96,3.76};
//double[] tmp2 = {3.17,3.63,3.47,3.44,3.39};
//double[] tmp3 = {3.64,3.72,3.91};
/*
double[][] data = new double[4][];
double[] tmp1 = {9,9.5,9.75,10,13,9.5};
double[] tmp2 = {11,10,10,11.75,10.5,15};
double[] tmp3 = {11.5,12,9,11.5,13.25,13};
double[] tmp4 = {13.25,11.5,12,13.5,11.5,12.35};
*/
data[0] = MiscUtil.DeepCopyDouble1(tmp1);
data[1] = MiscUtil.DeepCopyDouble1(tmp2);
data[2] = MiscUtil.DeepCopyDouble1(tmp3);
//data[3] = MiscUtil.DeepCopyDouble1(tmp4);
int[] categoryNum = new int[data.length];
for(int i=0;i<categoryNum.length;i++){
categoryNum[i] = data[i].length;
}
double[] concat = concatMult(data);
//double[] x1 = {1,2,3,4,5};
//double[] x2 = {11,12,13,14,15};
//double[] x1 = {66,57,81,62,61,60,73,59,80,55,67,70};
//double[] x2 = {64,58,45,43,37,56,44,42};
//double[] concat = concat(x1,x2);
int c = 100000000;
int seed = (int)(Math.random()*c);
//int seed = 10000;
String mval ="";
String mval2 ="";
for(int w=0;w<delta.length;w++){
MersenneTwisterFast mt = new MersenneTwisterFast(seed);
//double[] concatMu = plusMu(concat,delta[w],x1.length);
double[] concatMu = MiscUtil.DeepCopyDouble1(concat);
double kw = KruskalWallis(data);
//System.out.println("mw " + mw);
double DEV= 0.00000001;
int perm =10000;
//double to = sumUpto(concatMu,x1.length);
double to = tAnova(data);
double to2 = sAnova(data);
double to3 = fRatio(data);
System.out.println("fRatio " + to3);
System.out.println("to to2 " + to + " " + to2);
//System.out.println("mu " + delta[w]);
//for(int i=0;i<concat.length;i++){
//System.out.println("concat concatMu " + concat[i] + " " + concatMu[i]);
//}
//double[][] tmp0 = segment(concatMu,x1.length);
//double to2 = difMean(tmp0[0],tmp0[1]);
double[] t = new double[perm];
double[] t2 = new double[perm];
double[] fRatio = new double[perm];
double[] kwt = new double[perm];
for(int i=0;i<perm;i++){
//shuffleDoubleMZ(x1,mt);
//shuffleDoubleMZ(x2,mt);
shuffleDoubleMZ(concatMu,mt);
//concatMu = plusMu(concat,delta[w],x1.length);
//concat = concat(x1,x2);
//double[][] tmp = segment(concatMu,x1.length);
double[][] tmp = segmentMult(concatMu,categoryNum);
//t[i] = sumUpto(concatMu,x1.length);
//t2[i] = difMean(tmp[0],tmp[1]);
t[i]=tAnova(tmp);
t2[i]=sAnova(tmp);
fRatio[i]=fRatio(tmp);
kwt[i]=KruskalWallis(tmp);
//mwt[i] = MannWhitney(concat,x1.length);
System.out.println("t\tmwt\t" + t[i] +"\t" + t2[i] + "\t" + kwt[i] + "\t" + fRatio[i]);
}
// 度数分布・信頼区間など
int[] h=signPerm.hist(t,20);
double[][] hdouble = signPerm.histAccum(t,1000);
// 統計量の信頼区間を算出表示
double confA = 0.000001;
double[] confInterval = signPerm.confInterval(hdouble,confA);
System.out.println("confident interval\t" + confInterval[0] + "\t" + confInterval[1]);
// 度数分布を標準出力
// System.out.println(histSt);
System.out.println("Histogram\n
第４章 多標本１変数パーミュテーションテストの例
- 4.1 イントロダクション
- 対象となるのは、(１)２群の平均や位置の比較、（２）one-way ANOVA、（３）goodness-of-fit fir irdered categirucak varuabkes、（４）without-interaction two-way ANOVA
- 4.2 ２群の平均値の差の検定
- 全部で20人のデータ。２群に分けられて、群１は12人、群２は８人。２０人はある点数を持ち、群１と群２との間に点数の分布の違いがあるかどうかを検定したい。
x1 = {66,57,81,62,61,60,73,59,80,55,67,70};
x2 = {64,58,45,43,37,56,44,42};
- ２群間の平均値が異なるかの検定というのはよくやられる。
- Permutation testでは、20人をパーミュテーションして
することもできるし、それよりは場合数を減らしてやって、20人を12人と8人の２群に分ける場合わけとして
というパーミュテーションもできる
- Permutation 統計量は、『平均
がpermutationally equivalentなので、この方が簡単。
- この統計量で計算して、ｐとして、
オーダーの値が返る
public class univariateMultiSamplePerm {public static void main(String[] args) {
// 異なるdelta値での統計量を算出
int numseg = 100;
double[] delta = new double[numseg];
double initial=-20;
double dif = -initial*2/(double)numseg;
for(int i=0;i<delta.length;i++){
delta[i] = initial+dif*i;
}
//double[] delta = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
//double[] delta = {0};
//for(int i=0;i<delta.length;i++){
//delta[i] /= 0.1;
//}
double[] x1 = {1,2,3,4,5};
double[] x2 = {11,12,13,14,15};
//double[] x1 = {66,57,81,62,61,60,73,59,80,55,67,70};
//double[] x2 = {64,58,45,43,37,56,44,42};
double[] concat = concat(x1,x2);
int c = 100000000;
int seed = (int)(Math.random()*c);
//int seed = 10000;
String mval ="";
String mval2 ="";
for(int w=0;w<delta.length;w++){
MersenneTwisterFast mt = new MersenneTwisterFast(seed);
double[] concatMu = plusMu(concat,delta[w],x1.length);
double mw = MannWhitney(concatMu,x1.length);
//System.out.println("mw " + mw);
double DEV= 0.00000001;
int perm =10000;
double to = sumUpto(concatMu,x1.length);
//System.out.println("mu " + delta[w]);
//for(int i=0;i<concat.length;i++){
//System.out.println("concat concatMu " + concat[i] + " " + concatMu[i]);
//}
double[][] tmp0 = segment(concatMu,x1.length);
double to2 = difMean(tmp0[0],tmp0[1]);
double[] t = new double[perm];
double[] t2 = new double[perm];
double[] mwt = new double[perm];
for(int i=0;i<perm;i++){
//shuffleDoubleMZ(x1,mt);
//shuffleDoubleMZ(x2,mt);
shuffleDoubleMZ(concatMu,mt);
//concatMu = plusMu(concat,delta[w],x1.length);
//concat = concat(x1,x2);
double[][] tmp = segment(concatMu,x1.length);
t[i] = sumUpto(concatMu,x1.length);
t2[i] = difMean(tmp[0],tmp[1]);
mwt[i] = MannWhitney(concat,x1.length);
//System.out.println("t\tmwt\t" + t[i] +"\t" + t2[i] + "\t" + mwt[i]);
}
// 度数分布・信頼区間など
int[] h=signPerm.hist(t,20);
double[][] hdouble = signPerm.histAccum(t,1000);
// 統計量の信頼区間を算出表示
double confA = 0.000001;
double[] confInterval = signPerm.confInterval(hdouble,confA);
System.out.println("confident interval\t" + confInterval[0] + "\t" + confInterval[1]);
// 度数分布を標準出力
// System.out.println(histSt);
System.out.println("Histogram\n
疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)の目次
教科書
京都大学大学院医学研究科集中講義 全８限Part I(2007)
京都大学大学院医学研究科集中講義 全８限Part II(2008)
京都大学大学院医学研究科統計遺伝学講義 全８限(2009)
以下の関連サイトへのリンク道案内ページ
遺伝子多型解析Wiki
東大医科研HGCゲノム機能解析分野(旧)
理研リウマチチーム解析資料ページ
疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)
遺伝統計学のための統計学基礎 全９限
Haploviewを用いた連鎖不平衡マッピング解析実習 全９限
Arlequinを用いた集団遺伝学実習 全７限
一から始める遺伝子多型・形質関連解析
多数のGWAS/GWASメタアナリシスとその周辺解析用のアプリケーションへのリンク
複合遺伝性疾患とは
----------News----------
20万アクセス
2008/07/20
40万アクセス
2010/09/12
-------------------------
