Fisher's combined probability test - ryamadaの遺伝学・遺伝統計学メモ

昨日の記事でp値を併せる話しをした
$\sum_{i=1}^k -2\log(p_i)$ というフィッシャーの方法についても書いた
そのとき $-2\log(p)$ というのは、p値(が一様分布であるとして)に対応する自由度２のカイ二乗統計量にしているんだ、ということも書いた
それでは、「どうして自由度２」なの？となる
自由度はなんでもよいのか(最終的なp値の分布を一様分布にするだけならばなんでもよいが、 $p_1 \times p_2$ の大小の順序を守るには自由度２である必要がある)

n.pt <- 10000
# 統計量空間の次元
df <- 3
# df次元正規乱数
X <- matrix(rnorm(df*n.pt),ncol=df)
# この空間に複数の自由度１テストを表す単位ベクトルを置く
# テストは相互に非独立
v1 <- rnorm(df)
v1 <- v1/sqrt(sum(v1^2))
v2 <- v1 + rnorm(df)*0.1
v2 <- v2/sqrt(sum(v2^2))
v3 <- v1 + rnorm(df)*0.4
v3 <- v3/sqrt(sum(v3^2))
# この空間での乱点から、各テストベクトルへの射影の２乗はカイ二乗値(自由度１)
chi2.1 <- (X%*%v1)^2
chi2.2 <- (X%*%v2)^2
chi2.3 <- (X%*%v3)^2

# 自由度１でｐ値化する(自由度１で検定しているから〜ベクトルへの射影で考えることは自由度１ということだから)
p.1 <- pchisq(chi2.1,df=1,lower.tail=FALSE)
p.2 <- pchisq(chi2.2,df=1,lower.tail=FALSE)
p.3 <- pchisq(chi2.3,df=1,lower.tail=FALSE)

# Fisher'sの方法では自由度２で統計量に戻したが
# df.2 として自由度2で統計量に戻し、その上で足し合わせよう
df.2 <- 2
x <- qchisq(p.1,df.2,lower.tail=FALSE)+qchisq(p.2,df.2,lower.tail=FALSE)+qchisq(p.3,df.2,lower.tail=FALSE)
#x <- -2*log(p.1)-2*log(p.2)-2*log(p.3)

# 統計量の標本平均と標本分散を出す
m.x <- mean(x)
v.x <- var(x)
# ガンマ分布の平均はtheta k、分散はtheta^2 k なので
# 標本平均・標本分散とから、theta,kの推定値を出す
theta <- v.x/m.x
k <- m.x/theta
# 推定したガンマ分布パラメタを使って、ガンマ分布乱数を作ってみよう
y <- rgamma(n.pt,scale=theta,shape=k)

# 確かに、標本統計量との関係がとれている
xlim <- ylim <- range(c(x,y))
plot(sort(x),sort(y),xlim=xlim,ylim=ylim)
abline(0,1,col=2)


plot(data.frame(p.1,p.2,p.3))