スティール・ドゥワス(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

関連する記事

  • 平均的に分類する方法の考察(3)平均的に分類する方法の考察(3) 前回は、定量データをスコア順に並び替えたのち、この順番でグループに割り振っていく方法を見た。今回は、定量データをスコア順に並び替えるところは同じだが、割り振り方を変更することにより、より平均的に分類できないかを見ていく。 前回と同じく、100人の学生を3つのクラスA、B、Cに分ける方法を考えてみる。 まず、学生をスコア順にA、B、Cに一人ずつ割り振る。 次に […]
  • Ubuntu14.04でPython3に対応したmatplotlibを使用するための手順Ubuntu14.04でPython3に対応したmatplotlibを使用するための手順 Ubuntu14.04でPython3に対応したmatplotlibを使用するための手順をご紹介する。 1. […]
  • アンケートの自由記述からニーズを抽出する方法アンケートの自由記述からニーズを抽出する方法 アンケートは様々な場面で有効活用されていると思うが、特に、「年齢」などの数量回答や、「男性・女性」、「満足度の5段階評価」などの単一回答、「購入動機を3つまで回答してください」などの複数回答は、集計や分析がしやすいため、重点的に活用していることだろう。 回答者の立場で考えると、数量回答・単一回答・複数回答は、設問者からの問いかけに対し、返答があらかじめ用意されていること […]
  • R言語 CRAN Task View:系統学、特に比較方法R言語 CRAN Task View:系統学、特に比較方法 CRAN Task View: Phylogenetics, Especially Comparative Methodsの英語での説明文をGoogle翻訳を使用させていただき機械的に翻訳したものを掲載しました。 Maintainer: Brian O'Meara Contact: omeara.brian at […]
  • R言語 CRAN Task View:RグラフィカルモデルR言語 CRAN Task View:Rグラフィカルモデル CRAN Task View: gRaphical Models in Rの英語での説明文をGoogle翻訳を使用させていただき機械的に翻訳したものを掲載しました。 Maintainer: Soren Hojsgaard Contact: sorenh at […]
R スティール・ドゥワス(Steel-Dwass)法