スティール・ドゥワス(Steel-Dwass)法とは、テューキー(Tukey)法の多重比較に対応するノンパラメトリックな多重比較である。
スティール・ドゥワス法を簡単に言うと、正規分布を仮定しない各群間を順位を用いて多重比較で調べる方法である。

Rで、スティール・ドゥワス法を使う場合は、「スティール・ドゥワス(Steel-Dwass)の方法による多重比較」のページにソースコードが紹介されている。

ここでは、このソースコードを少しばかり改変したものをご紹介する。主な対応は次のとおりである。

  • 引数をformula形式に対応(response ~ group)
  • 群を表すベクトルを1からの数値ベクトル以外のものに対応
  • 返り値をデータフレームに変更

次が改変したソースコードである。標準パッケージstats内のソースコードも参考にさせていただいた。


steel.dwass.test <- function(x, ...) UseMethod("steel.dwass.test")

steel.dwass.test.default <- 
function(x, g, ...)
{
    if (is.list(x)) {
        if (length(x) < 2L)
            stop("'x' must be a list with at least 2 elements")
        if (!missing(g))
            warning("'x' is a list, so ignoring argument 'g'")
        DNAME <- deparse(substitute(x))
        x <- lapply(x, function(u) u <- u[complete.cases(u)])
        if (!all(sapply(x, is.numeric)))
            warning("some elements of 'x' are not numeric and will be coerced to numeric")
        k <- length(x)
        l <- sapply(x, "length")
        if (any(l == 0L))
            stop("all groups must contain data")
        g <- factor(rep.int(seq_len(k), l))
        x <- unlist(x)
    }
    else {
        if (length(x) != length(g))
            stop("'x' and 'g' must have the same length")
        DNAME <- paste(deparse(substitute(x)), "and",
                       deparse(substitute(g)))
        OK <- complete.cases(x, g)
        x <- x[OK]
        g <- g[OK]
        if (!all(is.finite(g)))
            stop("all group levels must be finite")
        g <- factor(g)
        k <- nlevels(g)
        if (k < 2L)
            stop("all observations are in the same group")
    }

    ## number of data in each group
    n.i <- table(g)

    ## calculate t-value
    t <- combn(levels(g), 2, function(ij) {
        i <- ij[1]
        j <- ij[2]
        ## ranking group i,j
        r <- rank(c(x[g == i], x[g == j]))                              
        ## test statistic
        R <- sum(r[1 : length(x[g == i])])                              
        ## total number of the 2 group data
        N <- n.i[i] + n.i[j]                                            
        ## expectation of test statistic
        E <- n.i[i] * (N + 1) / 2                                       
        ## variance of test statistic
        V <- n.i[i] * n.i[j] / (N * (N - 1)) * (sum(r^2) - N * (N + 1)^2 / 4) 
        ## return t-value
        return(abs(R - E) / sqrt(V))                                    
        })

    ## calculate p-value
    p <- ptukey(t * sqrt(2), length(n.i), Inf, lower.tail = FALSE)      

    df <- data.frame(
            comparison = combn(levels(g), 2, paste, collapse = ":"),
            t.value = t,
            p.value = p)
    return(df)
}

steel.dwass.test.formula <-
function(formula, data, subset, na.action, ...)
{
    if(missing(formula) || (length(formula) != 3L))
        stop("'formula' missing or incorrect")
    m <- match.call(expand.dots = FALSE)
    if(is.matrix(eval(m$data, parent.frame())))
        m$data <- as.data.frame(data)
    ## need stats:: for non-standard evaluation
    m[[1L]] <- quote(stats::model.frame)
    mf <- eval(m, parent.frame())
    if(length(mf) > 2L)
        stop("'formula' should be of the form response ~ group")
    DNAME <- paste(names(mf), collapse = " by ")
    names(mf) <- NULL
    y <- do.call("steel.dwass.test", as.list(mf))
    y
}

元のソースコードと結果が一致するかを確認する。次が元のソースコードの結果である。


> data <- c(
    6.9, 7.5, 8.5, 8.4, 8.1, 8.7, 8.9, 8.2, 7.8, 7.3, 6.8, # 第 1 群のデータ,11 例
    9.6, 9.4, 9.5, 8.5, 9.4, 9.9, 8.7, 8.1, 7.8, 8.8,      # 第 2 群のデータ,10 例
    5.7, 6.4, 6.8, 7.8, 7.6, 7.0, 7.7, 7.5, 6.8, 5.9,      # 第 3 群のデータ,10 例
    7.6, 8.7, 8.5, 8.5, 9.0, 9.2, 9.3, 8.0, 7.2, 7.9, 7.8  # 第 4 群のデータ,11 例
    )
> group <- rep(1:4, c(11, 10, 10, 11))                     # 群の識別変数
> Steel.Dwass(data, group)

           t           p
1:2 2.680234 0.036960431
1:3 2.539997 0.053980573
1:4 1.282642 0.574011771
2:3 3.746076 0.001031145
2:4 2.046776 0.170965537
3:4 3.384456 0.003976894

次が改変したソースコードの結果である。意図的に群を表すベクトルを文字列にして、formulaを用いている。


> df <- data.frame(response = data, group = rep(letters[1:4], c(11, 10, 10, 11)))
> steel.dwass.test(response ~ group, data = df)

  comparison  t.value     p.value
1        a:b 2.680234 0.036960431
2        a:c 2.539997 0.053980573
3        a:d 1.282642 0.574011771
4        b:c 3.746076 0.001031145
5        b:d 2.046776 0.170965537
6        c:d 3.384456 0.003976894

結果が同様になることが確認できた。

introducing-steel-dwass-in-r-1

関連する記事

  • R スティール(Steel)法R スティール(Steel)法 スティール(Steel)法とは、ダネット(Dunnett)法の多重比較に対応するノンパラメトリックな多重比較である。 スティール法を簡単に言うと、正規分布を仮定しない1つの対照群と2つ以上の処理群間を順位を用いて多重比較で調べる方法である。 Rで、スティール法を使う場合は、「スティール(Steel)の方法による多重比較」のページにソースコードが紹介されている。 […]
  • R MASSパッケージcorresp関数のエラー対処方法 RのMASSパッケージ内のcorresp関数を用いたとき、エラーに悩まされたので、皆様と共有しておく。 環境 Ubuntuのバージョン $ cat /etc/lsb-release […]
  • Journal of Statistical Software: 記事一覧 Journal of Statistical Software の記事一覧をご紹介する。英語での説明文をgoogle翻訳を使用させていただき機械的に翻訳したものを掲載した。 確認日:2017/03/24 論文数:1089 Introduction to stream: An Extensible Framework for Data Stream […]
  • UbuntuにRをインストールするための手順UbuntuにRをインストールするための手順 UbuntuにRをインストールするための手順をお伝えする。 Ubuntuのバージョン確認 Ubuntuのバージョンを確認するために、端末を起動し(ショートカットキー:Ctrl+Alt+t)、以下のコマンドを実行する。 DISTRIB_CODENAMEの行を控えておこう。下の内容では「trusty」の部分を控えておく。 $ cat […]
  • R 関数に時間制限を設ける方法 ある処理にとても時間が掛かるため、一定時間経過後はその処理を途中で打ち切りたいときがある。 例えば、for文で、あるループだけが重いため全体として時間が掛かってしまう場合、その処理を一旦スキップしてfor文の先に処理を進めたい、などである。 ここでは、そのひとつの解決策として、関数に時間制限を設けて、一定時間経過後はその関数を強制終了するコードをご紹介する。 ただし […]
R スティール・ドゥワス(Steel-Dwass)法