ryamada の弟子日記 （統計遺伝学 ・ 数学 ・ R）

■[ryamada本][多重検定][R]17.1.2 多重検定時のp値の期待値

• 師匠のブログのこちらを参照。

• あるサンプルセットに対し、独立な個の検定を行った。
• 一番小さい値はどのくらいかを考える。
  • 、つまり、あるサンプルセットに対し、１つの検定を行った。
  • その時の、値は0〜1の均一分布で、その期待値は0.5となる。

• あるサンプルセットに対し、独立な個の検定を行ったときの最小の値の期待値は？

# 「独立なk個の検定」を行ったときの最小のp値の期待値を求める。
# kは１回に行う独立な検定の個数。k=5個。
# Nは「独立なk個の検定」を施行する回数。N=1000回。
# つまり、独立な5個の検定を10000回行ったときの、最小のp値の期待値をシミュレーションする。

N <- 10000
k < -5
Xs <- matrix(runif(N*k),N,k)
Mins <- apply(Xs,1,min)
mean(Mins)
1/(k+1)

N<-10000
k<-5

N
[1] 10000
k
[1] 5

# 一様分布(min=0,max=1)に従う乱数をN*k個生成させる。
head(runif(n=N*k,min=0,max=1))
[1] 0.12189329 0.03397693 0.10567576 0.61410688 0.48816964 0.64391405
N*k
[1] 50000
length(runif(n=N*k,min=0,max=1))
[1] 50000

# 一様分布(min=0,max=1)に従う乱数をN*k個生成させ、
# N行k列の行列にする。（要素は列順に並べる）
# 行が１回の施行結果。
# 列が１回の施行の独立な５個の検定のp値。
head(matrix(runif(n=N*k,min=0,max=1),nrow=N,ncol=k,byrow=F))
          [,1]       [,2]      [,3]      [,4]      [,5]
[1,] 0.7458530 0.25808329 0.3861962 0.1314322 0.4038830
[2,] 0.8115991 0.09255577 0.2124800 0.8695457 0.2540804
[3,] 0.7692341 0.03943485 0.6950468 0.9832358 0.1309257
[4,] 0.3467765 0.59617752 0.9089028 0.1529305 0.6858364
[5,] 0.5660528 0.58521075 0.8597751 0.7954644 0.6009608
[6,] 0.5002232 0.47567372 0.6173216 0.5150203 0.6271794

# 10000行5列の行列が生成されることが分かる。
dim(matrix(runif(n=N*k,min=0,max=1),nrow=N,ncol=k,byrow=F))
[1] 10000     5

# 行列をXsに付値する。
Xs <- matrix(runif(n=N*k,min=0,max=1),nrow=N,ncol=k,byrow=F)

# 行列Xsの行ごとの最小値を求める。
# apply(行列,1,関数)　で、引数１は行を示す。
head(apply(Xs,1,min))
[1] 0.06932860 0.01760936 0.03746944 0.11688327 0.08593789 0.05195531

length(apply(Xs,1,min))
[1] 10000

# 各行の最小値(p値)がMinsに付値される。
Mins <- apply(Xs,1,min)

# 各行の最小値(p値)の平均が算出される。
mean(Mins)
[1] 0.1651497

# 「独立なk個の検定」を行ったときの最小のp値の期待値が1/(k+1)になるかを確認する。
1/(k+1)
[1] 0.1666667

# 「独立なk個の検定」を行ったときの最小のp値の期待値がほぼ1/(k+1)になることが分かった。

• あるサンプルセットに「独立な個の検定」を行ったときに、番目に小さい値の期待値は？

# あるサンプルセットに「独立なk個(k=10)の検定」の施行をN(N=1000)回行った場合のシミュレーション。
N <- 1000
k <- 10
Xs <- matrix(runif(N*k),N,k)
Sorted <- matrix(apply(Xs,1,sort),N,k,byrow=TRUE)
plot(apply(Sorted,2,mean))

N <- 1000
k <- 10

N
[1] 1000
k
[1] 10

# 「独立なk個(k=10)の検定」をN(N=1000)回行った場合のp値を
# N行k列の行列として生成する。（ここまでは同じ）
# 今回は独立な検定がk=10個。
# 施行回数がN=1000回。
 head(matrix(runif(n=N*k,min=0,max=1),nrow=N,ncol=k,byrow=F))
           [,1]      [,2]      [,3]       [,4]       [,5]       [,6]      [,7]
[1,] 0.06023224 0.8509787 0.2893429 0.16933603 0.81986046 0.14548850 0.7211218
[2,] 0.64031295 0.4219804 0.2813534 0.94702960 0.61250452 0.02528587 0.6811834
[3,] 0.69209141 0.8233987 0.3419297 0.03508989 0.08794964 0.33928631 0.4850175
[4,] 0.04609405 0.1299671 0.4189534 0.55838081 0.04271290 0.25727952 0.1604458
[5,] 0.36383763 0.5724046 0.4750929 0.70507964 0.56467813 0.43842702 0.9518281
[6,] 0.60831131 0.9462469 0.3352221 0.64257990 0.11072018 0.95654425 0.8776467
          [,8]      [,9]     [,10]
[1,] 0.3937247 0.7042586 0.9109493
[2,] 0.3479346 0.5568777 0.2359857
[3,] 0.9395950 0.3455220 0.5115202
[4,] 0.4364787 0.2851158 0.7837359
[5,] 0.6311239 0.6492596 0.9332059
[6,] 0.7413012 0.7415595 0.4949100

# 行列は1000行10列になる。
dim(matrix(runif(n=N*k,min=0,max=1),nrow=N,ncol=k,byrow=F))
[1] 1000   10

Xs <- matrix(runif(n=N*k,min=0,max=1),nrow=N,ncol=k,byrow=F)

# 行を昇順に並べ替える。
# apply(Xs,1,sort) この表現式がどのようになるかを試す。
x < -matrix(sample(1:12,12,replace=F),nrow=3,ncol=4,byrow=F)
x
     [,1] [,2] [,3] [,4]
[1,]   11    8    2    3
[2,]    7    6    5   12
[3,]    4    9    1   10

apply(x,1,sort)
     [,1] [,2] [,3]
[1,]    2    5    1
[2,]    3    6    4
[3,]    8    7    9
[4,]   11   12   10
# 結果は、n行目は昇順ソートされ、n列目に配置される。
# m行n列の行列は、n行m列の行列になる。
# 昇順ソートされたm行は、n行目となる。
# apply(Xs,1,sort)は、1000行10列だったので、10行1000列の行列となる。

# Sorted <- matrix(apply(Xs,1,sort),nrow=N,ncol=k,byrow=TRUE) この表現式を考える。
x
     [,1] [,2] [,3] [,4]
[1,]   11    8    2    3
[2,]    7    6    5   12
[3,]    4    9    1   10

s <- apply(x,1,sort)
s
     [,1] [,2] [,3]
[1,]    2    5    1
[2,]    3    6    4
[3,]    8    7    9
[4,]   11   12   10

matrix(s,nrow=3,ncol=4,byrow=T)
     [,1] [,2] [,3] [,4]
[1,]    2    3    8   11
[2,]    5    6    7   12
[3,]    1    4    9   10

t(s) # 行列を転置しても同じ結果になる。
     [,1] [,2] [,3] [,4]
[1,]    2    3    8   11
[2,]    5    6    7   12
[3,]    1    4    9   10
# 元の行列を各行について昇順ソートした行列となった。

Sorted <- matrix(apply(Xs,1,sort),nrow=N,ncol=k,byrow=TRUE)

dim(Sorted)
[1] 1000   10

head(Sorted)
           [,1]       [,2]      [,3]      [,4]      [,5]       [,6]      [,7]
[1,] 0.06932860 0.07322612 0.4970635 0.5384743 0.9186192 0.01760936 0.2616903
[2,] 0.03746944 0.31653477 0.4647877 0.6278671 0.8971589 0.11688327 0.2341168
[3,] 0.08593789 0.36115840 0.5654091 0.8559646 0.8679085 0.05195531 0.6148712
[4,] 0.35854093 0.43637956 0.7411988 0.9542136 0.9700300 0.40958418 0.4423390
[5,] 0.14023083 0.28714450 0.3068490 0.4612815 0.5631951 0.14943989 0.4773012
[6,] 0.14223002 0.26894026 0.3975098 0.4886170 0.8488136 0.42328664 0.5770654
          [,8]      [,9]     [,10]
[1,] 0.6161067 0.8977828 0.9293514
[2,] 0.3167031 0.6411231 0.8007600
[3,] 0.6615881 0.6640558 0.9775101
[4,] 0.4683319 0.8252976 0.9379207
[5,] 0.5132673 0.6821402 0.8988834
[6,] 0.8143040 0.8861237 0.9529390

# 列ごとの平均を算出し、ベクトルにする。
head(apply(Sorted,2,mean))
[1] 0.1635236 0.3337900 0.5015983 0.6748455 0.8375406 0.1705790

length(apply(Sorted,2,mean))
[1] 10

# 独立した10個の検定を1000回行ったときの、各検定ごとの平均値(p値)がベクトルとなる。
apply(Sorted,2,mean)
 [1] 0.1635236 0.3337900 0.5015983 0.6748455 0.8375406 0.1705790 0.3409559
 [8] 0.5062973 0.6722367 0.8368279

# 独立した10個の検定を1000回行ったときの、各検定ごとの平均値(p値)を
# 昇順ソートしてみる。
sort(apply(Sorted,2,mean))
 [1] 0.1635236 0.1705790 0.3337900 0.3409559 0.5015983 0.5062973 0.6722367
 [8] 0.6748455 0.8368279 0.8375406

# 上の結果が、独立なk(k=10)個の検定を行ったとき、i番目に小さい
# p値の期待値がi/(k+1)になるか確認してみる。
k <- 10
pith <- NULL
for(i in 1:10){
 pith <- c(pith,i/(k+1))
}

pith
 [1] 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 0.54545455
 [7] 0.63636364 0.72727273 0.81818182 0.90909091

sort(apply(Sorted,2,mean))
 [1] 0.1635236 0.1705790 0.3337900 0.3409559 0.5015983 0.5062973 0.6722367
 [8] 0.6748455 0.8368279 0.8375406
# ん〜、pithとsort(apply(Sorted,2,mean))の結果が同じような違うような・・・

# 「独立なk(k=10)個の検定」のp値の期待値をプロットしてみる。
plot(apply(Sorted,2,mean))

• 「独立な個の検定」を回行ったときの、各回ごとの昇順ソートした値の分布を調べる。

# 独立なk(k=10)個の検定。
# N(N=1000)回行う。
N <- 1000
k <- 10
Xs <- matrix(runif(n=N*k,min=0,max=1),nrow=N,ncol=k,byrow=F)
Sorted <- matrix(apply(Xs,1,sort),nrow=N,ncol=k,byrow=TRUE)
# Sortedは、行にシミュレーションの回数。
# 列は、各シミュレーションごとの検定のp値が昇順に並んでいる。

# 列ごとに昇順ソートし、その昇順ソートした結果は列のままとなる。
# 1000行10列の行列は、1000行10列の行列で行列の次元は変わらない。
doubleSorted <- apply(Sorted,2,sort)
dim(doubleSorted)
[1] 1000   10

matplot(doubleSorted,type="l")

Permalink | コメント(0) | トラックバック(0) | 03:51