First  upload: Feb 02, 2007
　　　　　　　　　　　　　　　　　　　　Last modified: Feb 10, 2007

1. はじめに

2. 方法

シミュレーションプログラムのソース

===== begin =====
BF.sim <- function(type, n1=20, n2=20, sd2=1, loop=5000)
{
        calc <- function(x1, x2, n1, n2)
        {
                # median test                             メディアン検定
                md <- median(c(x1, x2))
                r1 <- sum(x1 < md)
                r2 <- sum(x2 < md)
                tbl <- matrix(c(r1, r2, n1-r1, n2-r2), 2, 2)
                chisq.n <- chisq.test(tbl, correct=FALSE)$p.value
                chisq.c <- chisq.test(tbl, correct=TRUE)$p.value
                fisher <- fisher.test(tbl)$p.value
                # U test                          マン・ホイットニーの U 検定
                U.n <- wilcox.test(x1, x2, correct=FALSE, exact=FALSE)$p.value
                U.c <- wilcox.test(x1, x2, correct=TRUE, exact=FALSE)$p.value
                # t test                          t 検定
                t.w <- t.test(x1, x2, var.equal=FALSE)$p.value
                t.n <- t.test(x1, x2, var.equal=TRUE)$p.value
                # equality of valiance                    等分散性の検定
                vp <- var.test(x1, x2)$p.value          # パラメトリック
                an <- ansari.test(x1, x2)$p.value       # ノンパラメトリック
                return(c(
                        chisq.n=chisq.n,                # 1 メディアン検定を　カイ二乗検定で
                        chisq.c=chisq.c,                # 2 メディアン検定を　連続性の補正をするカイ二乗検定検定で
                        fisher=fisher,                  # 3 メディアン検定を　フィッシャーの正確検定で
                        U.n=U.n,                        # 4 マン・ホイットニーの U 検定
                        U.c=U.c,                        # 5 連続性の補正をするマン・ホイットニーの U 検定
                        t.w=t.w,                        # 6 ウエルチの方法による t 検定
                        t.n=t.n,                        # 7 普通の t 検定
                        vp=vp,                          # 8 等分散性の検定
                        an=an))                         # 9 アンサリ検定
        }
        if (type == "normal") {                         # 正規分布
                g1 <- matrix(rnorm(n1*loop, sd=1), loop, n1)
                g2 <- matrix(rnorm(n2*loop, sd=sd2), loop, n2)
        }
        else if (type == "lognormal") {                 # 対数正規分布みたいな外形をした分布（平均値，標準偏差を調整する）
                g1 <- matrix(scale(rlnorm(n1*loop, sdlog=0.4)), loop, n1)
                g2 <- matrix(scale(rlnorm(n2*loop, sdlog=0.4))*sd2, loop, n2)
        }
        else if (type == "uniform") {                   # 一様分布
                g1 <- matrix(scale(runif(n1*loop)), loop, n1)
                g2 <- matrix(scale(runif(n2*loop))*sd2, loop, n2)
        }
        result <- t(sapply(1:loop, function(i) calc(g1[i,], g2[i,], n1, n2)))
        return(result)
}
draw <- function(type, n1, n2, loop, result)                    # 結果を図にする
{
        color <- c("", "red","magenta","orange","blue","purple","brown","grey1")
        linetype <- c(4, 4, 1, 2, 3, 4, 1, 2, 3, 4)
        rng <- range(as.vector(result[,2:8]))
        pdf(paste(type, n1, n2, "pdf", sep="."), onefile=FALSE, width=800/72, height=600/72)
        old <- par(mar=c(4,4,1,1), cex=2)
        plot(c(1,31), rng, type="n", ylim=c(0, 0.25), xaxt="n", xlab="sd2 / sd1",
                ylab="Type I error", main=paste(type, ", n1=", n1, ", n2=", n2, ", trial=", loop, sep=""))
        sapply(2:8, function(i) lines(1:31, result[,i], col=color[i], lwd=4, lty=linetype[i]))
        abline(h=0.05, v=16, lty=3, lwd=1)
        legend(18, 0.25, c("chisq.test", "chisq.test(c)", "fisher.test", "wilcox.test",
                "wilcox.test(c)", "t.test(w)", "t.test"), col=color[2:8], lwd=6)
        axis(1, at=1:31, labels=c("1/4", "1/3.8", "1/3.6", "1/3.4", "1/3.2", "1/3",
                "1/2.8", "1/2.6", "1/2.4", "1/2.2", "1/2", "1/1.8", "1/1.6", "1/1.4",
                "1/1.2", "1", "1.2", "1.4", "1.6", "1.8", "2", "2.2", "2.4", "2.6",
                "2.8", "3", "3.2", "3.4", "3.6", "3.8", "4"))
        par(old)
        dev.off()
}
driver2 <- function(type, n1, n2, sd2, loop)
{
        result <- BF.sim(type, n1, n2, sd2, loop)
        resultp <- result < 0.05
        return(colSums(resultp)/loop)
}
driver <- function(type, n1, n2, sd2, loop=5000)        # 分布の形，各群のサンプルサイズ，標準偏差の比
{
        sink(paste(type, n1, n2, "log", sep="."))
        cat("\n\nn1=", n1, "   n2=", n2, "   trial=", loop, "   distribution=", type, "\n\n")
        result <- t(sapply(sd2, function(sd2) c(sd2=sd2, driver2(type, n1, n2, sd2, loop))))
        print(result)
        draw(type, n1, n2, loop, result)
        sink()
}
# 第二群の標準偏差のベクトル　1〜4 を 0.2 刻みで，逆方向も
sd2 <- seq(1, 4, by=0.2)                                # 大きくなる方
sd2 <- c(1/rev(sd2), sd2[-1])                           # 小さくなる方
loop <- 50                                              # 試行回数
driver("normal", 5, 5, sd2, loop)                       # シミュレーション実行
===== end =====

図1. 左右対称な分布

図2. 歪みのある分布

図3. 一様な分布

3. 結果

参考：三群以上の比較の場合の結果

3.1. 左右対称な分布の場合

3.1.1. サンプルサイズが等しい場合（n1=20, n2=20）

3.1.2. サンプルサイズが約 1:2 の場合（n1=15, n2=25)

3.1.3. サンプルサイズが 1:3 の場合（n1=10, n2=30）

3.2. 歪みのある分布の場合

3.2.1. サンプルサイズが等しい場合（n1=20, n2=20）

3.2.2. サンプルサイズが約 1:2 の場合（n1=15, n2=25)

3.2.3. サンプルサイズが 1:3 の場合（n1=10, n2=30）

3.3. 一様な分布の場合

3.3.1. サンプルサイズが等しい場合（n1=20, n2=20）

3.3.2. サンプルサイズが約 1:2 の場合（n1=15, n2=25)

3.3.3. サンプルサイズが 1:3 の場合（n1=10, n2=30）

3.4. その他のサンプルサイズの場合

4. 結論

5. まとめ

追記

参考文献

ご意見・ご希望，間違いの指摘等がありましたら，
下記宛へ E-mail いただくと，たいへんうれしく思います。