多重比較の問題を考える前に，p値の分布を見ておく - ケィオスの時系列解析メモランダム

検定を繰り返し何度もやると，「p値の調整をしろ」と言われます．これは多重比較の問題と呼ばれるもので，何度も検定を繰り返すと，p値が小さいものが観測されやすくなります．

この多重比較の問題を理解するためには，

*　p値は絶対的な正しさの目安にはなれないこと

*　p値はきまぐれ (確率的)で，たまたま小さかったり，たまたま大きかったりすること

を知っておく必要があります．

まったく差がない比較では，p値は一様分布に従って決まる

　次の質問をされたらなんと答えますか．

　「２つのデータを比較するときに，両者にまったく差がなければ，p値はどんな値をとる？」

統計的に言えば「帰無仮説が正しいとき」のp値は？という問題です．

答えは，1に近い値です (ウソです)．そうではなく，答えは，

　０から１までの値をとる一様分布に従ってランダムに決まる

ということです．

　この事実は数学的に証明されています．

The p-value follows a uniform distribution under the null hypothesis | The Book of Statistical Proofs

数値実験

　疑似乱数を使って，p値が一様分布になるか調べてみた結果が以下です．今回はRを使いましたが，コンピュータが生成する疑似乱数では，奇妙な現象が見られることがあります．

ここでは，平均０，分散１の標準正規分布に従う乱数をn個ずつ，２セット生成し，t検定，F検定とShapiro-Wilk検定をしてみました．Shapiro-Wilk検定は，２群の比較ではありませんので，１群のみ使いました．

上の図 (n = 10) では，t検定 (0付近が小さい)とShapiro-Wilk検定 (1付近が大きい)で，一様分布からずれている構造が見られますが，これは疑似乱数の性質によるものだと思います．

　以下は nを増やした結果です．

上の図 (n = 20) では，Shapiro-Wilk検定で，1付近の確率密度が高くなっていますが，これは，Rの乱数生成のアルゴリズムの問題だと思います．本当は，下の図のように一様分布に従います．　

上の図 (n = 50) で，やっと，一様分布に近い構造が見られました．　

　以上のように (結果はイマイチでしたが)，差がないとき，帰無仮説が棄却されないとき，p値は一様分布に従います．ですので，まったく差がなくても，検定をN回繰り返せば，p < 0.05になる例は，平均 0.05 N 個くらい観測されます．

　ですので，検定を100回繰り返して，そのうち5個くらい p < 0.05になっていても，差がある根拠にはならないということです．

　では，実際に差がある２群を比較すれば，p値の分布はどのようになるでしょうか？これについては次回，説明します．

Rスクリプト

　Rスクリプトの例を掲載しておきます．

mu1 <- 0
mu2 <- 0
sd1 <- 1
sd2 <- 1
#####################
N <- 10
#####################
n.test <- 10000
p.t <- c()
p.F <- c()
p.S <- c()
for(i in 1:n.test){
  x1 <- rnorm(N,mu1,sd1)
  x2 <- rnorm(N,mu2,sd2)
  #####################
  p.t[i] <- t.test(x1,x2)$p.value
  p.F[i] <- var.test(x1,x2)$p.value
  p.S[i] <- shapiro.test(x2)$p.value
}
#####################################
bin.size <- 0.05
###
par(mar=c(5,6,2,2),mfrow=c(1,3))
hist(p.t,breaks=seq(0,1,bin.size),border=NA,col=4,freq=FALSE,xlab="p value",main="t-test",las=1,cex.lab=2,cex.axis=1.5,ylab="",cex.main=2,ylim=c(0,1.5))
mtext("Density",2,line=4,cex=1.4)
##
hist(p.F,breaks=seq(0,1,bin.size),border=NA,col=4,freq=FALSE,xlab="p value",main="F-test",las=1,cex.lab=2,cex.axis=1.5,ylab="",cex.main=2,ylim=c(0,1.5))
mtext("Density",2,line=4,cex=1.4)
##
hist(p.S,breaks=seq(0,1,bin.size),border=NA,col=4,freq=FALSE,xlab="p value",main="Shapiro-Wilk test",las=1,cex.lab=2,cex.axis=1.5,ylab="",cex.main=2,ylim=c(0,1.5))
mtext("Density",2,line=4,cex=1.4)