トレンドテスト



SNPジェノタイプデータの、Armitage's trend testについての記事を3月に書いた(こちら)。

式はY^2=¥frac{N(N(r_1+2r_2)-R(n_1+2n_2))^2}{RS(N(n_1+4n_2)-(n_1+2n_2)^2)}

ただし、2x3の分割表が

r_0,r_1,r_2,R

s_0,s_1,s_2,S

n_0,n_1,n_2,N

のような場合

少し修正。

3月の記事の計算式は、Trend-Chiテストと別称する方がわかりやすい。

Cochran-Armitage trendテストでは、上式の1箇所だけ(分子の先頭のNをN-1へ)、検体数 N の扱いを N-1 に変えるので、値が少し変わってくる。カイ自乗値のYates 補正に相当。

補正式は

Y^2=¥frac{(N-1)(N(r_1+2r_2)-R(n_1+2n_2))^2}{RS(N(n_1+4n_2)-(n_1+2n_2)^2)}

ただし、あちらこちらの記事でこのあたりは曖昧に扱われているようである。

再度、こちらのサイトの比較例を見ると、SASのPatrick Royston’s ptrend commandは「Trend-Chi」の値を返すらしい(筆者は未検証)。

Cytel社のStatXactでは、Cochran-Armitage trendテストの近似統計量とより近似度を高めた『正確統計量』(こちら)を返すが、この近似統計量では、「補正あり」のCochran-Armitage trendテストがなされている(ようだ)。

Trend-Chiテストと「補正あり」のCochran-Armitage trendテストを載せたエクセルはこちら

さらに補足で、群馬大青木先生のRの関数(こちら) Cochran.ArmitageはTrend-Chiテストの値を返します。


public static double armitage(double[][] data){
double N=0;
double R=0;
double S=0;
double[] n = new double[data[0].length];
for(int i=0;i<data.length;i++){
for(int j=0;j<data[i].length;j++){
N+=data[i][j];
if(i==0){
R+=data[i][j];
}else if(i==1){
S+=data[i][j];
}
n[j]+=data[i][j];
}
}
//System.out.println("data[0][1] data[0][2] "+data[0][1] + " " + data[0][2]);
//System.out.println("N="+N+" R="+R+" S="+S);
//System.out.println("n1n2n3="+n[0]+" "+n[1]+" "+n[2]);
double Y=((N-1)*Math.pow((N*(data[0][1]+2*data[0][2])-R*(n[1]+2*n[2])),2))/
(R*S*(N*(n[1]+4*n[2])-Math.pow((n[1]+2*n[2]),2)));
//System.out.println("Y="+Y);
double s=Math.pow(3,2);
//System.out.println("3^2="+s);
return Y;
}