目的

スティール（Steel）の方法による多重比較を行う。

使用法

Steel(data, group)

引数

data     データベクトル
group    群変数ベクトル

ソース

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

# スティール（Steel）の方法による多重比較
Steel <- function(   data,                                   # データベクトル
                        group)                                  # 群変数ベクトル
{
        get.rho <- function(ni)                                      # ρを計算する
        {
                k <- length(ni)
                rho <- outer(ni, ni, function(x, y) { sqrt(x/(x+ni[1])*y/(y+ni[1])) })
                diag(rho) <- 0
                sum(rho[-1, -1])/(k-2)/(k-1)
        }

        OK <- complete.cases(data, group)                    # 欠損値を持つケースを除く
        data <- data[OK]
        group <- factor(group[OK])                           # 群変数データを factor に変関る
        ni <- table(group)                                   # 各群のデータ数
        a <- length(ni)                                              # 群の数
        control <- data[group == 1]                          # 対照群のデータを取り出す
        n1 <- length(control)                                        # 対照群のデータ数
        t <- numeric(a)
        rho <- ifelse(sum(n1 == ni) == a, 0.5, get.rho(ni))  # ρを決める
        for (i in 2:a) {
                r <- rank(c(control, data[group == i]))              # 対照群と対照群をまとめて順位をつける
                R <- sum(r[1:n1])                            # 検定統計量
                N <- n1+ni[i]                                        # 対照群と対照群のデータ数の合計
                E <- n1*(N+1)/2                                      # 検定統計量の期待値
                V <- n1*ni[i]/N/(N-1)*(sum(r^2)-N*(N+1)^2/4) # 検定統計量の分散
                t[i] <- abs(R-E)/sqrt(V)                     # 検定統計量（t 分布による漸近近似）
        }
        result <- cbind(t, rho)[-1,]                         # 結果をまとめる
        rownames(result) <- paste(1, 2:a, sep=":")
        return(result)
}


使用例

data <- c(
    50, 55, 65, 63, 60, 68, 69, 60, 52, 49,  # 第 1 群（対照群）のデータ，10 例
    80, 86, 74, 66, 79, 81, 70, 62, 60, 72,  # 第 2 群（処理群）のデータ，10 例
    42, 48, 58, 63, 62, 55, 63, 60, 53, 45   # 第 3 群（処理群）のデータ，10 例
)
group <- rep(1:3, each=10)                   # 群を表す数値
Steel(data, group)

出力結果例

> Steel(data, group)
           t rho
1:2 2.952566 0.5
1:3 1.175674 0.5

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