目的

平均値の差の多重比較（Ryan 法，Tukey 法）を行う

使用法

m.multi.comp(n, me, s, alpha=0.05, method=c("ryan", "tukey"))

引数

n       各群の例数
me      各群の平均値
s       各群の標準偏差
alpha   有意水準
method  方法（省略時は Ryan 法）

ソース

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

# ライアンの方法とチューキーの方法による平均値の対比較
m.multi.comp <- function(    n,                              # 標本サイズベクトル
                                me,                             # 平均値ベクトル
                                s,                              # 標準偏差ベクトル
                                alpha=0.05,                     # 有意水準
                                method=c("ryan", "tukey"))      # 方法
{
        printf <- function(fmt, ...)
        {
                cat(sprintf(fmt, ...))
        }

        check <- function(s, b)                                      # 検定しようとしている二群が，それまでに有意でないとされた二群に挟まれているか
        {
                if (ns.n > 1) {
                        for (i in 1:ns.n) {
                                if (ns.s[i] <= s && s <= ns.b[i] && ns.s[i] <= b && b <= ns.b[i]) {
                                        return(FALSE)           # 検定するまでもなく有意でないとする
                                }
                        }
                }
                return(TRUE)                                    # 検定しなくてはならない
        }

        k <- length(n)                                               # 群の数
        stopifnot(k == length(me), k == length(s), n > 0, floor(n) == n, s > 0)
        method <- match.arg(method)                          # 引数の補完
        o <- order(me)                                               # 平均値の大きさの順位
        sn <- n[o]                                           # 並べ替えた標本サイズ
        sm <- me[o]                                          # 並べ替えた平均値
        ss <- s[o]                                           # 並べ替えた標準偏差
        nt <- sum(sn)                                                # 全体の標本サイズ
        mt <- sum(sn*sm)/nt                                  # 全体の平均値
        dfw <- nt-k                                          # 群内平方和の自由度
        vw <- sum((sn-1)*ss^2)/dfw                           # 群内分散
        num.significant <- ns.n <- 0
        ns.s <- ns.b <- numeric(k*(k-1)/2)                        # 有意でない群の記録用
        for (m in k:2) {                                        # 検定対象の選定
                for (small in 1:(k-m+1)) {
                        big <- small+m-1
                        if (check(small, big)) {
                                t0 <- (sm[big]-sm[small])/sqrt(vw*(1/sn[big]+1/sn[small]))   # 検定統計量
                                if (method == "ryan") {                                         # Ryan の方法
                                        P <- pt(t0, dfw, lower.tail=FALSE)*2                 # 有意確率
                                        nominal.alpha <- 2*alpha/(k*(m-1))                   # 名義的有意水準
                                        result <- P <= nominal.alpha                              # 検定結果
                                }
                                else { # Tukey の方法
                                        t0 <- t0*sqrt(2)
                                        q <- (qtukey(0.05, k, dfw, lower.tail=FALSE) + qtukey(0.05, m, dfw, lower.tail=FALSE))/2
                                        WSD <- q*sqrt(vw*(1/sn[big]+1/sn[small])/2)          # Wholly Significant Difference
#                                       P <- ptukey(t0, m, dfw, lower.tail=FALSE)            # 有意確率　以下で置き換え
                                        P <- uniroot(function(p)((qtukey(p, k, dfw, lower.tail=FALSE)+
                                                                  qtukey(p, m, dfw, lower.tail=FALSE))/2 -t0), c(0,1),tol=1e-7)$root
                                        result <- sm[big]-sm[small] >= WSD                   # 検定結果
                                }
                                if (result) {                   # 有意であるとき
                                        num.significant <- 1
                                        o.small <- o[small]  # 群番号は，入力時（実験時）の順序番号を表示する
                                        o.big <- o[big]
                                        printf("mean[%2i]=%7.5f vs. mean[%2i]=%7.5f : diff.= %7.5f, ",
                                                o.small, me[o.small], o.big, me[o.big], me[o.big]-me[o.small])
                                        if (method == "ryan") {
                                                printf("t=%7.5f : P=%7.5f, alpha'=%7.5f\n", t0, P, nominal.alpha)
                                        }
                                        else {
                                                printf("WSD=%7.5f : t=%7.5f : P=%7.5f\n", WSD, t0, P)
                                        }
                                }
                                else {                          # 有意でないとき
                                        ns.n <- ns.n+1
                                        ns.s[ns.n] <- small
                                        ns.b[ns.n] <- big
                                }
                        }
                }
        }
        if (num.significant == 0) {                             # 有意差のある群は一つもなかった
                print("Not significant at all.")
        }
}


使用例

m.multi.comp(c(8, 11, 22, 6), c(135.83, 160.49, 178.35, 188.06), c(19.59, 12.28, 15.01, 9.81), method="ryan")
m.multi.comp(c(8, 11, 22, 6), c(135.83, 160.49, 178.35, 188.06), c(19.59, 12.28, 15.01, 9.81), method="tukey")

出力結果例

> m.multi.comp(c(8, 11, 22, 6), c(135.83, 160.49, 178.35, 188.06), c(19.59, 12.28, 15.01, 9.81), method="ryan")
mean[ 1]=135.83000 vs. mean[ 4]=188.06000 : diff.= 52.23000, t=6.53866 : P=0.00000, alpha'=0.00833
mean[ 1]=135.83000 vs. mean[ 3]=178.35000 : diff.= 42.52000, t=6.96308 : P=0.00000, alpha'=0.01250
mean[ 2]=160.49000 vs. mean[ 4]=188.06000 : diff.= 27.57000, t=3.67279 : P=0.00066, alpha'=0.01250
mean[ 1]=135.83000 vs. mean[ 2]=160.49000 : diff.= 24.66000, t=3.58814 : P=0.00085, alpha'=0.02500
mean[ 2]=160.49000 vs. mean[ 3]=178.35000 : diff.= 17.86000, t=3.26998 : P=0.00212, alpha'=0.02500
> m.multi.comp(c(8, 11, 22, 6), c(135.83, 160.49, 178.35, 188.06), c(19.59, 12.28, 15.01, 9.81), method="tukey")
mean[ 1]=135.83000 vs. mean[ 4]=188.06000 : diff.= 52.23000, WSD=21.34697 : t=9.24706 : P=0.00000
mean[ 1]=135.83000 vs. mean[ 3]=178.35000 : diff.= 42.52000, WSD=15.57114 : t=9.84728 : P=0.00000
mean[ 2]=160.49000 vs. mean[ 4]=188.06000 : diff.= 27.57000, WSD=19.14118 : t=5.19411 : P=0.00258
mean[ 1]=135.83000 vs. mean[ 2]=160.49000 : diff.= 24.66000, WSD=16.11328 : t=5.07440 : P=0.00195
mean[ 2]=160.49000 vs. mean[ 3]=178.35000 : diff.= 17.86000, WSD=12.80554 : t=4.62444 : P=0.00482
平均値の差（diff.）が WSD より大きければ有意（P値は参考）

 解説ページ