目的

比率の差の多重比較（Ryan 法，Tukey 法）を行う

使用法

p.multi.comp(n, r, alpha=0.05, method=c("ryan", "tukey"))

引数

n       各群の例数
r       各群の陽性数
alpha   有意水準
method  方法（省略時は Ryan 法）

ソース

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

# ライアンの方法とチューキーの方法による比率の対比較
p.multi.comp <- function(    n,                              # 標本サイズ
                                r,                              # 陽性サイズ
                                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(r), k == length(n), n > 0, r >= 0, r <= n, floor(n) == n, floor(r) == r)
        method <- match.arg(method)                          # 引数の補完
        o <- order(r/n)                                              # 標本比率の大きさの順位
        sr <- r[o]
        sn <- n[o]
        sr.sn <- sr/sn                                               # 割合
        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)) {
                                prop <- sum(sr[small:big]) / sum(sn[small:big])              # 推定母比率
                                se <- sqrt(prop*(1-prop)*(1/sn[small]+1/sn[big]))    # 標準誤差
                                diff <- sr.sn[big]-sr.sn[small]                              # 比率の差
                                if (method == "ryan") { # Ryan の方法
                                        nominal.alpha <- 2*alpha/(k*(m-1))           # 名義的有意水準
                                        rd <- se*qnorm(nominal.alpha/2, lower.tail=FALSE)    # 差があると見なせる差の大きさ
                                        z <- if (se == 0) 0 else diff/se
                                        p <- pnorm(z, lower.tail=FALSE)*2
                                        result <- p <= nominal.alpha
                                }
                                else { # Tukey の方法
                                        if (m == k) {           # 最大群と最小群の比較のとき
                                                qq <- q <- qtukey(alpha, k, Inf, lower.tail=FALSE)
                                        }
                                        else {                  # 上記以外の群の比較のとき
                                                qq <- (q+qtukey(alpha, m, Inf, lower.tail=FALSE))/2
                                        }
                                        WSD <- se*qq/sqrt(2)
                                        result <- diff != 0 && diff >= WSD
                                }
                                if (result) { # 有意であるとき
                                        num.significant <- 1
                                        printf("p[%2i]=%7.5f vs. p[%2i]=%7.5f : diff.= %7.5f, ",
                                                o[small], sr.sn[small], o[big], sr.sn[big], diff)
                                        if (method == "ryan") {
                                                printf("RD=%7.5f : P=%7.5f, alpha'=%7.5f\n", rd, p, nominal.alpha)
                                        }
                                        else {
                                                printf("WSD=%7.5f\n", WSD)
                                        }
                                }
                                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.")
        }
}


使用例

p.multi.comp(c(30, 35, 47, 21, 45), c(2, 4, 14, 13, 39), method="ryan")
p.multi.comp(c(30, 35, 47, 21, 45), c(2, 4, 14, 13, 39), method="tukey")

出力結果例

> p.multi.comp(c(30, 35, 47, 21, 45), c(2, 4, 14, 13, 39), method="ryan")
p[ 1]=0.06667 vs. p[ 5]=0.86667 : diff.= 0.80000, RD=0.32472 : P=0.00000, alpha'=0.00500
p[ 1]=0.06667 vs. p[ 4]=0.61905 : diff.= 0.55238, RD=0.33341 : P=0.00001, alpha'=0.00667
p[ 2]=0.11429 vs. p[ 5]=0.86667 : diff.= 0.75238, RD=0.30528 : P=0.00000, alpha'=0.00667
p[ 1]=0.06667 vs. p[ 3]=0.29787 : diff.= 0.23121, RD=0.23054 : P=0.00979, alpha'=0.01000
p[ 2]=0.11429 vs. p[ 4]=0.61905 : diff.= 0.50476, RD=0.32612 : P=0.00007, alpha'=0.01000
p[ 3]=0.29787 vs. p[ 5]=0.86667 : diff.= 0.56879, RD=0.26479 : P=0.00000, alpha'=0.01000
p[ 3]=0.29787 vs. p[ 4]=0.61905 : diff.= 0.32118, RD=0.29877 : P=0.01239, alpha'=0.02000

> p.multi.comp(c(30, 35, 47, 21, 45), c(2, 4, 14, 13, 39), method="tukey")
p[ 1]=0.06667 vs. p[ 5]=0.86667 : diff.= 0.80000, WSD=0.31555
p[ 1]=0.06667 vs. p[ 4]=0.61905 : diff.= 0.55238, WSD=0.32546
p[ 2]=0.11429 vs. p[ 5]=0.86667 : diff.= 0.75238, WSD=0.29800
p[ 1]=0.06667 vs. p[ 3]=0.29787 : diff.= 0.23121, WSD=0.22695
p[ 2]=0.11429 vs. p[ 4]=0.61905 : diff.= 0.50476, WSD=0.32104
p[ 3]=0.29787 vs. p[ 5]=0.86667 : diff.= 0.56879, WSD=0.26067
p[ 3]=0.29787 vs. p[ 4]=0.61905 : diff.= 0.32118, WSD=0.30102

 解説ページ