「統計学関連なんでもあり」の過去ログ--- 046

No.21073　Rでのt検定とU検定で例数8の時だけ現れる謎の現象、あるいは私のミス　　【shimp】　2014/05/30(Fri) 19:00

実験計画における例数の重要性をできるだけ生々しく同僚に示すために，t検定とマン・ホイットニーのU検定それぞれで，p<0.05の「勝ち取りやすさ」の例数依存性をシミュレーションしてグラフ化する関数を書いてみました(末尾記載)。
引数の例示を兼ねて簡単に説明します。

-------------
n=3から10までのそれぞれの場合について，
A群平均 1，B群平均 5
A群SD 2，B群SD 4
の条件下にて，rnorm(1000)でそれら2群を作ってt検定とU検定を行い，p値が0.05を下回った割合を記録し，最後に横軸をn数でプロットします。

sim.norm.p_n(1000, 3, 10, 1, 5, 2, 4)
-------------

一見，期待した通りの結果を示すのですが，どういう訳かn=8の時にだけ，U検定とt検定それぞれでp<0.05となる割合が接近するのです。
稚拙とはいえシンプルな関数なので，誤りを探しても見つけられずにいます。
t検定やU検定の仕組み，あるいはRでの実装に由来する現象という事は考えられるでしょうか?

宜しくお願い致します。

#---------------------

sim.norm.p_n <- function(trial, #試行回数
size.min, #試行する最小例数
size.max, #試行する最大例数
A.mean, #A群の平均値
B.mean, #B群の平均値
A.sd, #A群の標準偏差
B.sd #B群の標準偏差
)
{
n <- trial
result <- data.frame(t.result=1:n, U.result=1:n) #n回試行時の一時データフレーム
result2 <- data.frame(size=size.min:size.max, t.result2=size.min:size.max, U.result2=size.min:size.max)
#最終データ格納用のデータフレーム
result2$size <- NA #初期化
result2$t.result2 <- NA
result2$U.result2 <- NA
y = size.min:size.max #例数範囲の設定
for (j in y) #例数カウンタ
{
result$t.result <- NA #一時データフレーム初期化
result$U.result <- NA
x = 1:n #試行回数範囲の設定
for (i in x) #試行カウンタ
{
A.sample <- rnorm(j, mean=A.mean, sd=A.sd) #A群データ生成
B.sample <- rnorm(j, mean=B.mean, sd=B.sd) #B群データ生成
t.result <- t.test(A.sample,B.sample) #t.test
U.result <- wilcox.test(A.sample,B.sample) #U.test
result$t.result[i] <- t.result$p.value #一時データフレームにp値収納
result$U.result[i] <- U.result$p.value
}
result2$size[j-size.min+1] <- j #最終データに例数を格納
result2$t.result2[j-size.min+1] <- sum(result$t.result < 0.05)/n #最終データにp<0.05となる割合を格納
result2$U.result2[j-size.min+1] <- sum(result$U.result < 0.05)/n
}
x.axis.lim <- c(size.min, size.max) #以下，グラフ化
y.axis.lim <- c(min(result2$t.result2,result2$U.result2),max(result2$t.result2,result2$U.result2))
title.main <- paste("試算条件:各群平均=",A.mean,",",B.mean,"，各群SD=",A.sd,",", B.sd)
par(mfrow=c(2,1))
plot(result2$size, result2$t.result2,pch=20,col="red", xlim=x.axis.lim, ylim=y.axis.lim, xlab="", ylab="")
par(new=TRUE)
plot(result2$size, result2$U.result2,pch=20,col="blue",xlim=x.axis.lim, ylim=y.axis.lim, xlab="例数", ylab="p<0.05となる割合",main=title.main)
legend("bottomright",c("t検定の場合","U検定の場合"),pch=19,col=c("red","blue"),cex=0.7,bty="n")

#----- ここから下はおまけ -----
norm.x.min <- min(A.mean, B.mean)-4*(max(A.sd,B.sd))
norm.x.max <- max(A.mean, B.mean)+4*(max(A.sd,B.sd))
curve(dnorm(x,mean=A.mean, sd=A.sd), from=norm.x.min, to=norm.x.max, ylab="")
curve(dnorm(x,mean=B.mean, sd=B.sd), from=norm.x.min, to=norm.x.max, add=TRUE)
}

No.21074　Re: Rでのt検定とU検定で例数8の時だけ現れる謎の現象，あるいは私のミス　　【青木繁伸】　2014/05/31(Sat) 05:01

原因は，ウィルコクソンの符号順位和検定の検定統計量が不連続なためです。

wil <- function(n) {
	q <- qwilcox(0.025, n, n)
	W <- 0:q
	d <- dwilcox(W, n, n)
	p <- cumsum(d)
	return(round(cbind(W, d, p), 5))
}

というようなプログラムを書いて，n=7, 8, 9 の場合の結果を見ると，

> wil(7) # n = 7

      W       d       p
 [1,] 0 0.00029 0.00029
 [2,] 1 0.00029 0.00058
 [3,] 2 0.00058 0.00117
 [4,] 3 0.00087 0.00204
 [5,] 4 0.00146 0.00350
 [6,] 5 0.00204 0.00554
 [7,] 6 0.00321 0.00874
 [8,] 7 0.00437 0.01311
 [9,] 8 0.00583 0.01894 ここまでが棄却域
[10,] 9 0.00758 0.02652

> wil(8) # n = 8

       W       d       p
 [1,]  0 0.00008 0.00008
 [2,]  1 0.00008 0.00016
 [3,]  2 0.00016 0.00031
 [4,]  3 0.00023 0.00054
 [5,]  4 0.00039 0.00093
 [6,]  5 0.00054 0.00148
 [7,]  6 0.00085 0.00233
 [8,]  7 0.00117 0.00350
 [9,]  8 0.00171 0.00521
[10,]  9 0.00218 0.00738
[11,] 10 0.00295 0.01033
[12,] 11 0.00373 0.01406
[13,] 12 0.00490 0.01896
[14,] 13 0.00598 0.02494 ここまでが棄却域
[15,] 14 0.00754 0.03248

> wil(9) # n = 9

       W       d       p
 [1,]  0 0.00002 0.00002
 [2,]  1 0.00002 0.00004
 [3,]  2 0.00004 0.00008
 [4,]  3 0.00006 0.00014
 [5,]  4 0.00010 0.00025
 [6,]  5 0.00014 0.00039
 [7,]  6 0.00023 0.00062
 [8,]  7 0.00031 0.00093
 [9,]  8 0.00045 0.00138
[10,]  9 0.00062 0.00200
[11,] 10 0.00082 0.00282
[12,] 11 0.00107 0.00389
[13,] 12 0.00142 0.00531
[14,] 13 0.00179 0.00710
[15,] 14 0.00228 0.00938
[16,] 15 0.00284 0.01222
[17,] 16 0.00352 0.01573
[18,] 17 0.00426 0.01999 ここまでが棄却域
[19,] 18 0.00516 0.02515

まあ，たまたまといえばたまたま，n = 8 のときは，以下のように，実際の p/2 値と有意水準（α/2 = 0.025）の差が他に比べると小さいということ。

> diff <- numeric(20)
> for (n in 3:20) {
+ 	diff[n] <- round(pwilcox(qwilcox(0.025, n, n)-1, n, n)-0.025, 5)
+ }
> data.frame(n=3:20, diff=diff[3:20])
    n     diff
1   3 -0.02500
2   4 -0.01071
3   5 -0.00913
4   6 -0.00444
5   7 -0.00606
6   8 -0.00006 この差が小さい
7   9 -0.00501
8  10 -0.00337
9  11 -0.00135
10 12 -0.00255
11 13 -0.00294
12 14 -0.00013 他にも差の小さい場合がある
13 15 -0.00233
14 16 -0.00159
15 17 -0.00066
16 18 -0.00146
17 19 -0.00017 他にも差の小さい場合がある
18 20 -0.00045

No.21075　Re: Rでのt検定とU検定で例数8の時だけ現れる謎の現象，あるいは私のミス　　【shimp】　2014/05/31(Sat) 21:32

有難う御座いました。
なるほど，ところどころに絶妙の値が巡ってくるのですね。

> wil(8) # n = 8
[14,] 13 0.00598 0.02494

> wil(14) # n = 14
[56,] 55 0.00268 0.02487

> wil(19) # n = 19
[114,] 113 0.00171 0.02483

実のところ，最初に回答を拝見した時はどこを見れば良いのやら途方に暮れたのですが，
判ってしまえば当たり前のことにすら思えてきました。
ちょっと感動です。

重ねてお礼申し上げます。
有難う御座いました。

No.21076　Re: Rでのt検定とU検定で例数8の時だけ現れる謎の現象，あるいは私のミス　　【shimp】　2014/05/31(Sat) 21:45

あと，プログラムの書き方が参考になりました。

こういう基本的なところはもちろん，
> diff <- numeric(20)

こういう使い方が大変参考になりました。
W <- 0:q
d <- dwilcox(W, n, n)

せっかくRを使っているのだから，すこし自分のコードを見直してみます。
併せてお礼申し上げます。