目的

ボンフェローニの方法，ホルムの方法，シェイファーの方法，ホランド・コペンハーバーの方法による多重比較を行う。

使用法

Bonferroni(data, group, method="Bonferroni", "Holm", "Shaffer", "Holland-Copenhaver"), alpha=0.05)

引数

data     観察値ベクトル
group    群変数ベクトル
method   有意水準の調整方法。標準値は Bonferroni。
alpha    有意水準（デフォールトは 0.05）

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/Bonferroni.R", encoding="euc-jp")

# ボンフェローニの方法および関連する手法による多重比較
Bonferroni <- function(      data,                                           # 観察値ベクトル
                        group,                                          # 群変数ベクトル
                        method=c("Bonferroni", "Holm", "Shaffer", "Holland-Copenhaver"),
                        # ボンフェローニ，ホルム，シェーファー，ホランド・コペンハーバー
                        alpha=0.05)                                     # 有意水準
{
        OK <- complete.cases(data, group)                            # データと群変数がともに欠損値ではないデータを抽出する
        data <- data[OK]
        group <- factor(group[OK])
        method <- match.arg(method)                                  # method に与えられた略語による指定を補完する
# ボンフェローニ法　Bonferroni
        n <- table(group)                                            # 各群のデータ数
        a <- length(n)                                                       # 群の数
        k <- a*(a-1)/2                                                       # 2 群の総当たり数
        m <- tapply(data, group, mean)                                       # 各群の平均値
        v <- tapply(data, group, var)                                        # 各群の不偏分散
        phi <- length(data)-a                                                # 自由度
        V <- sum((n-1)*tapply(data, group, var))/phi                 # 誤差分散
        t <- combn(a, 2, function(ij) abs(diff(m[ij]))/sqrt(sum(V/n[ij])))   # abs(m[i]-m[j])/sqrt(V*(1/n[i]+1/n[j])))
        P <- pt(t, phi, lower.tail=FALSE)*2                          # P 値を計算
        result1 <- cbind("n"=n, "mean"=m, "variance"=v)                      # 各群の，標本サイズ，平均値，不偏分散
        rownames(result1) <- paste("Group", 1:a, sep="")             # 表側を作る
        result2 <- cbind("t value"=t, "P value"=P)                   # 二群の比較結果の，t 値，P 値
        rownames(result2) <- combn(a, 2, paste, collapse=":")                # 表側を作る
        if (method == "Bonferroni") {                                   # ボンフェローニ法 Bonferroni のとき，
                P.threshold <- alpha/k                                       # P 値の限界値と，
                judge <- result2[,2] <= P.threshold                       # 判定（有意のときには 1，有意でないときには 0 が表示される）
        }
        else {                                                          # ボンフェローニ法以外のとき
                result2 <- result2[order(result2[,2]),]
                if (method != "Holm" && a > 9) {                  # ホルムの方法以外で，群数が 9 より大きいときには対応できないのでホルムの方法にする
                        warning("Too many groups. Use Holm.")
                        method = "Holm"
                }
                if (method == "Holm") {                                 # ホルムの方法 Holm のとき
                        P.threshold <- alpha/(k:1)                   # P 値の限界値を設定
                }
                else {
                        h <- c(3,1,1, 6,rep(3,3),2,1, 10,rep(6,4),4,4:1, 15,rep(10,5),rep(7,3),6,4,4:1,
                               21,rep(15,6),rep(11,4),10,9,7,7:1, 28,rep(21,7),rep(16,5),15,13,13:1,
                               36,rep(28,8),rep(22,6),21,rep(18,3),16,16,15,13,13:1)[(a*(a-1)*(a-2)/6):((a^3-a)/6-1)]
                        P.threshold <- if (method == "Shaffer") alpha/h else 1-(1-alpha)^(1/h)       # P 値の限界値を設定
                }
                judge <- cumprod(result2[,2] <= P.threshold) != 0 # 判定（有意のときには 1，有意でないときには 0 が表示される）
                if (method == "Holm" || method == "Shaffer" || method == "Holland-Copenhaver") {
                        min.pos <- which.min(judge)                  # 途中で有意でないものが出てきたら，それ以降は全て結果を保留
                        if (judge[min.pos] == 0) {                      # 初めて出てくる 0 の位置
                                judge[min.pos:k] <- NA                       # 保留を表すという意味で NA を代入
                        }
                }
        }
        result2 <- cbind(result2,                                    # 結果の追加
                    "P threshold"=P.threshold, "judge"=judge)
        return(list(method=method, alpha=alpha, n.of.tests=k, result2=result2,
                remark="judge=1 significant, judge=0 not significant, judge=NA suspend", result1=result1, phi=phi, V=V))
}


使用例

data <- c(
    10.7, 9.7, 8.5, 9.4, 8.8, 8.4, 10.6, # 第 1 群のデータ，7 例
    8.1, 8.3, 8.7, 6.9, 5.7, 9.5, 6.7,   # 第 2 群のデータ，7 例
    7.9, 7.5, 7.4, 9.2, 5.7, 8.3, 9.7,   # 第 3 群のデータ，7 例
    6.2, 7.1, 5.5, 4.7, 6.3, 6.9, 7.5    # 第 4 群のデータ，7 例
)
group <- rep(1:4, each=7)                # 群を表す数値

Bonferroni(data, group)                     # ボンフェローニの方法
Bonferroni(data, group, method="Holm")      # ホルムの方法
Bonferroni(data, group, method="Shaffer")   # シェイファーの方法
Bonferroni(data, group, method="Holland-Copenhaver")	# ホランド・コペンハーバーの方法

出力結果例

> Bonferroni(data, group)	# ボンフェローニの方法
$method
[1] "Bonferroni"

$alpha
[1] 0.05

$n.of.tests
[1] 6

$result2
      t value      P value P threshold judge
1:2 2.8351656 9.148816e-03 0.008333333     0
1:3 2.4168625 2.362029e-02 0.008333333     0
1:4 5.0893546 3.315126e-05 0.008333333     1
2:3 0.4183031 6.794449e-01 0.008333333     0
2:4 2.2541890 3.358738e-02 0.008333333     0
3:4 2.6724921 1.331927e-02 0.008333333     0

$remark
[1] "judge=1 significant, judge=0 not significant, judge=NA suspend"

$result1
       n     mean  variance
Group1 7 9.442857 0.8961905
Group2 7 7.700000 1.7333333
Group3 7 7.957143 1.7195238
Group4 7 6.314286 0.9414286

$phi
[1] 24

$V
[1] 1.322619

> Bonferroni(data, group, method="Holm")	# ホルムの方法
$method
[1] "Holm"
　　：　中略
$result2
      t value      P value P threshold judge
1:4 5.0893546 3.315126e-05 0.008333333     1
1:2 2.8351656 9.148816e-03 0.010000000     1
3:4 2.6724921 1.331927e-02 0.012500000    NA
1:3 2.4168625 2.362029e-02 0.016666667    NA
2:4 2.2541890 3.358738e-02 0.025000000    NA
2:3 0.4183031 6.794449e-01 0.050000000    NA
　　：　後略

> Bonferroni(data,group, method="Shaffer")	# シェイファーの方法
$method
[1] "Shaffer"
　　：　中略
$result2
      t value      P value P threshold judge
1:4 5.0893546 3.315126e-05 0.008333333     1
1:2 2.8351656 9.148816e-03 0.016666667     1
3:4 2.6724921 1.331927e-02 0.016666667     1
1:3 2.4168625 2.362029e-02 0.016666667    NA
2:4 2.2541890 3.358738e-02 0.025000000    NA
2:3 0.4183031 6.794449e-01 0.050000000    NA
　　：　後略

> Bonferroni(data,group, method="Holland-Copenhaver")	# ホランド・コペンハーバーの方法
$method
[1] "Holland-Copenhaver"
　　：　中略
$result2
      t value      P value P threshold judge
1:4 5.0893546 3.315126e-05 0.008512445     1
1:2 2.8351656 9.148816e-03 0.016952428     1
3:4 2.6724921 1.331927e-02 0.016952428     1
1:3 2.4168625 2.362029e-02 0.016952428    NA
2:4 2.2541890 3.358738e-02 0.025320566    NA
2:3 0.4183031 6.794449e-01 0.050000000    NA
　　：　後略

 参考文献
永田靖・吉田道弘(1997)「統計的多重比較法の基礎」サイエンティスト社