p値の分布

仮説検定を繰り返してp値の分布を調べる.

最初に正規分布に従う乱数を作る.以下のコードは平均0，標準偏差1の正規分布に従う乱数を10個作る

rnorm(10, 0, 1)

##  [1]  0.6049668 -1.3782510  1.3977854  0.8731403 -2.5349832 -0.2550702
##  [7] -1.7642630 -0.1718184 -0.8647703 -0.9769723

次の関数coinはn個の正規分布(平均0，標準偏差sd)に従う乱数を二列作り，その平均値に差があるかどうかをt検定する関数である

coin内の二つの分布の平均は同じだから，帰無仮説は正しい．

coin<-function(n,sd){
x <- rnorm(n, 0, sd)#n個の正規乱数をxに代入
y <- rnorm(n, 0, sd)#n個の正規乱数をyに代入
z<-t.test(x,y,var.equal=T)#
z$p.value}#t検定の結果からp値だけを返す

coin(1000,1)

## [1] 0.2551716

次の関数testは関数coinをk回反復する

test<-function(k,n,sd){
x <- c()   # 空の（要素数ゼロの）リストを作る
for (i in 1:k) {
  y <- coin(n,sd) #　n個のサンプルを使って，t検定の結果からp値のみ返す
  x <- c(x,  y)   }#結果ベクトルに追加していく
x}

p値のヒストグラムをプロットする

hist(test(k=1000,n=10000,sd=5), breaks = 100)

シミュレーションにより，帰無仮説が正しいときは，p値の分布が一様分布に近づく様子が分かる

次に関数を一般化して，帰無仮説が正しくない場合を表現しよう

coin2<-function(n,m1,m2,sd){
x <- rnorm(n, m1, sd)#n個の正規乱数をxに代入
y <- rnorm(n, m2, sd)#n個の正規乱数をyに代入
z<-t.test(x,y,var.equal=T)#
z$p.value}#t検定の結果からp値だけを返す

今度は帰無仮説が正しくない場合， m1,m2を明示的に指定して，母集団の平均値に差がある状態を作る．

coin2(n=100,m1=0,m2=1,sd=1)

## [1] 7.2594e-15

test2<-function(k,n,m1,m2,sd){
x <- c()   # 空の（要素数ゼロの）リストを作る
for (i in 1:k) {
  y <- coin2(n,m1,m2,sd) #　n個のサンプルを使って，t検定の結果からp値のみ返す
  x <- c(x,  y)   }#結果ベクトルに追加していく
x}

p値のヒストグラムをプロットする

hist(test2(k=1000,n=10000,m1=0,m2=0.5,sd=5), breaks = 100)

対立仮説が正しい場合は，p値の分布が0付近に集中することが分かる