第２章　単純な検定問題 - ryamadaの遺伝学・遺伝統計学メモ

2.0.0
- 帰無仮説は、統計量としてどういうものを選ぶかを定め、対立仮説は、統計量がどうなっていたら、確率的に起こりにくいと判断するべきかの情報を与える
2.0 サンプル数20の対応のある観測データでの例
- 各種検定(permutation testを含む)で算出してみた(一部、適用することの妥当性を欠くものも)結果は次の通り



Student's t	0.0001134

Peason's product moment correlation	0.00001001

Wilcoxson's Signed Rank test	0.04776

McNemar test (Binomial test)	0.001288

Permutation test	0.0003

Permutation testでは、下記で出るXの和と、Xの標準化した和とを用いることができる。今、パーミュテーションしたデータセットについて、それぞれの統計量を算出し、その値の大小で順序づけすると、『正しい統計量』にが作る順序は一致するはずである。実際に、2.1以下で用いたデータについて、permutation testで用いたXの和、Xの標準化した和、Wilcoxon、Student's tを算出し、その順序関係を調べると、『和』『標準化した和』『Student't t』では一致するが、『Wilcoxon』では、かなり高い一致を示すが、多少の順序のずれが出る。これを図で説明したものはこちら。
以下、詳細

2.1 対応のある観測データ
- データの例
  - 20個の独立サンプルについて、あるインターベンションの前後に観測値YA,YBをとる。XはYA-YBである。この観測値は、アンケートの質問項目の応答数に相当しており、その分布には適当な分布が仮定し得ない



i	YA	YB	X

1	19	14	5

2	22	23	-1

3	18	13	5

4	18	17	1

5	24	20	4

6	30	22	8

7	26	30	-4

8	28	21	7

9	15	11	4

10	30	29	1

11	16	17	-1

12	25	20	5

13	22	18	4

14	19	17	2

15	27	22	5

16	23	21	2

17	24	21	3

18	18	15	3

19	28	24	4

20	27	22	5

Rにてデータ処理するために



> ya<-c(19,22,18,18,24,30,26,28,15,30,16,25,22,19,27,23,24,18,28,27)

> yb<-c(14,23,13,17,20,22,30,21,11,29,17,20,18,17,22,21> x<-ya-yb

> x

 [1]  5 -1  5  1  4  8 -4  7  4  1 -1  5  4  2  5  2  3  3  4  5,21,15,24,22)

> mean(ya)

[1] 22.95

> mean(yb)

[1] 19.85

> mean(x)

[1] 3.1

> var(x)

[1] 8.2

どのような因子が考慮されるか
- インターベンションのもたらす影響
- インターベンションとは無関係であるが、個人に依存するずれ
- 個人特異的なインターベンション反応性
- 個人に依存しかつ、その他の因子の影響を受けるランダムなぶれ

2.2 Student't t-testを適用してみる(相関検定もしてみる)
- RでStudent't t-test
  - p<0.001,p>0.0001である

> t.test(x)

One Sample t-test

data: x
t = 4.8414, df = 19, p-value = 0.0001134
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
1.759811 4.440189
sample estimates:
mean of x
3.1

RでPeasonの積率相関係数算出とその検定
- Student's t-testではX=YA-YBの値の分布について、パラメトリックな分布(平均値やそこからの分散がパラメタ化されている)を仮定した。これがそもそも、怪しいが、相関係数の方では、YAとYBの両方にパラメトリックな分布を仮定している。これはXに仮定するよりさらに怪しいことに留意しながら、先へ進む

> cor.test(ya,yb)

Pearson's product-moment correlation

data: ya and yb
t = 6.0578, df = 18, p-value = 1.001e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.5906741 0.9259794
sample estimates:
cor
0.8190952

2.3 Signed rank test を適用してみる
- X が連続変数ながら、それがもたらす確率分布は不明であるとする。Xの期待値については、有限であるとの制約はない。

> wilcox.test(ya,yb)

Wilcoxon rank sum test with continuity correction

data: ya and yb
W = 273.5, p-value = 0.04776
alternative hypothesis: true mu is not equal to 0

Warning message:
タイがあるため、正確な p 値を計算することができません in: wilcox.test.default(ya, yb)

2.4 McNemar test を適用してみる
Xはもはや、連続である必要もなくなった
- 今回のデータは、YAが大きい場合が17、小さい場合が3、合わせて20ペアなので、２項検定に相当する(McNemar testは２ｘ２のときは２項検定の近似)。

> binom.test(17,20,alternative = "greater")

Exact binomial test

data: 17 and 20
number of successes = 17, number of trials = 20, p-value = 0.001288
alternative hypothesis: true probability of success is greater than 0.5
95 percent confidence interval:
0.6563362 1.0000000
sample estimates:
probability of success
0.85

2.5 Permutation test を適用してみる
- 帰無仮説においては、YA と YB とは等しいから、その差であるXの分布は、０を中心に偏りがない。Xiについて、その符号を付け替えた観測データのパターン数は、 $2^n; n=20$ 。ただし、今回のサンプルデータには存在していないが、Xi の値が０と観測されているときには、符号の付け替えは何も違いをもたらさないので、一般的には、観測サンプル数nに対して、そのXiが０のサンプル数vを引いて、 $2^{n-v}$ 通りの値の割り付けパターンがある。このすべての場合について、データセットが与える値の分布を取る。たとえば、 $￥sum_{i=0}^{n-v}X_i$ は、YA>YB なる仮設において、大きな正の値をとればとるほど、起こりにくい値となるので、この $￥sum_{i=0}^{n-v}X_i$ について、全通り(permutations)のTの値における観測Tの相対的順位を持って、棄却率が得られる。
- $2^{20}=1048576$
2.6 Conditional Monte Carlo (CMC) (Permutation) Method を適用してみる
- CMC
  - すべての順列(permutation)をしらみつぶしに計算する代わりに、観測データの順列をランダムに繰り返し発生させ(Random permutations)、それが作る統計量の分布を全順列が作る分布の代用とする
- 試しに符号割付けCMCを実行すると、教科書で0.0003が得られたとあるようにそれに近い数値が得られる
- Permutation 分布の近似の実際(目安)
  - サンプルサイズが小さいとき(n<25)
    - 全permutationsを算出するのがよい
  - nが200を越えるような大サイズの場合
    - 中心極限定理が使えて、正規化した統計量を正規分布より棄却確率が算出できる
    - $Ko=￥frac{￥sum_{i=1}^{20}X_i}{￥sqrt{￥sum_{i=1}^{20}X_i^2}}=3.324$ となり、これは、平均０、分散１の正規分布検定でそのP値は、両側検定で0.000444である。Permutation test 的には、符号入れ替えを行って、 $K^*=￥frac{￥sum_{i=1}^{20}X_i}{￥sqrt{￥sum_{i=1}^{20}X_i^2}}$ の分布をとってやると、これが平均０、分散１の正規分布になるということ(やらなくても分布がわかっているので、計算せずにPを出してやる。提示サンプルデータでは、サンプル数が20と条件を満たさないが、それでも与えられるP値は、それほど悪くない数値である。
  - この中間の場合
    - CMCをB回(前の例では10000回)で分布を近似する
2.7 Permutation Confidence Interval for
- 違いがないとする帰無仮説について、確率分布をpermutationにて得る方法について述べてきた。ある対立仮説を規定するパラメタについても、同様にpermutationで分布を得ることは可能で、それをすべてのpermutationsについて調べることで、その分布域の上下限域を知ることも可能である
- の近似分布
  - サンプルサイズが大きければ、中心極限定理が適用できて $K(￥delta)=￥frac{￥sum_{i=1}^{n}(X_i-￥delta)}{(￥sum_{i=1}^{n}(X_i-￥delta)^2)^{1/2}}$ 。正規分布にて、上下の確率分布を確認すればCIが得られる
- CMCによるCIの近似
  - Permutation 統計量がモンテカルロ試行分得られている。これが作る分布の $-￥frac{￥alpha }{2} -- ￥frac{￥alpha }{2}$ の間が $￥alpha$ 信頼区間である。P値算出にあたって行った試行について、これを満足するCI下限値、CI上限値を算出する
  - 上下の対称がとれないのでそれぞれ別個に算出する
  - 複数の $￥Delta$ について、permutationによってその分布をとり、それが、求めたい信頼区間の上限・下限のそれぞれに相当するような確率を満たすまで、試してみる。信頼区間の精度は、permutation の試行回数に依存する点に留意



public class signPerm {
	public static void main(String[] args) {

		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[2][perm];

		double[] ya ={19, 22, 18, 18, 24, 30, 26, 28, 15, 30,

				16, 25, 22, 19, 27, 23, 24, 18, 28, 27};

		//double[] yb = {1, 23, 20, 17, 20, 22, 30, 21,20, 29,

				//17, 30, 18, 17, 22, 21, 21, 15, 30, 30};

		//double[] yb = {24, 23, 20, 17, 20, 22, 30, 21,20, 29,

				//17, 30, 18, 17, 22, 21, 21, 15, 30, 30};

		double[] yb = {14, 23, 13, 17, 20, 22, 30, 21, 11, 29,

				17, 20, 18, 17, 22, 21, 21, 15, 24, 22};

		String mval = "";

		String mval2 ="";

		int c = 100000000;

		int seed = (int)(Math.random()*c);

		for(int m=0;m<delta.length;m++){

			System.out.println("\n\n=================\nd value_" + m + "\n=================\n");

			double[] x = new double[ya.length];

			

			

			

			

			//int seed = 10000;

			MersenneTwisterFast mt = new MersenneTwisterFast(seed);

			MersenneTwisterFast mt2 = new MersenneTwisterFast(seed);

			for(int i=0;i<x.length;i++){

				x[i]=ya[i]-yb[i]+delta[m];

				

			}

			double[] to = {sum(x),K(x)};

			//double Ko = K(x);

			//複数の統計量を比較できるように

			permT = SumAndKAndWilcoxonSignPerm(x,perm,mt);

			

			//permutation Pの算出

			double Pvalue = PermTestP(to[0],permT[0]);

			System.out.println("PermP\t" + Pvalue);

			

			//度数分布・信頼区間など

			int[] h=hist(permT[0],20);

			double[][] hdouble = histAccum(permT[0],1000);

//			統計量の信頼区間を算出表示

			double confA = 0.000001;

			double[] confInterval = confInterval(hdouble,confA);

			

			System.out.println("confident interval\t" + confInterval[0] + "\t" + confInterval[1]);

			

//			度数分布を標準出力

//			System.out.println(histSt);

			System.out.println("Histogram\n