スティール・ドゥワス(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 (!is.factor(mf[[2]])) {
      mf[[2]] <- factor(mf[[2]])
    }
    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のバージョンが3.6.xでは動作していたが、4.1.0では次のエラーが出て、動作しない事例をご報告いただきました。
上記のコードは以下のエラーが出ないように変更したコードになります。
Rのバージョン3.6.3と4.1.0で調べてみたところ、eval()関数またはparent.frame()関数の挙動がバージョン間で異なっておりました。
ご報告いただいた方へ、この場を借りてお礼申し上げます。


teel.dwass.test.default(c(6.9, 7.5, 8.5, 8.4, 8.1, 8.7, 8.9,  でエラー: 
  all group levels must be finite 
R スティール・ドゥワス(Steel-Dwass)法