読者です 読者をやめる 読者になる 読者になる

疫学と医療統計学と遺伝学と時々、大学院生活

疫学を専門とする大学院生の研究に関する備忘録的ページ。

RでSteel-Dwass(スティールドゥワス)検定をする

今日はノンパラメトリックな多重比較の一例として、Steel.Dwass testのRでの実行方法を記載する。
2016年6月13日にコメントされたようにsteel-dwassについてNSM3パッケージにあるコマンドについて追記した。

方法1

Steel.dwass検定は下記のHP(群馬大学・青木先生による)を参考して、自作の関数を作成し、実行することになる。
http://aoki2.si.gunma-u.ac.jp/R/Steel-Dwass.html

早速、上のHPから関数をcopyアンドpasteする。
Steel.Dwass <- function(data,group)
{
        OK <- complete.cases(data, group)
        data <- data[OK]
        group <- group[OK]
        n.i <- table(group)
        ng <- length(n.i)
        t <- combn(ng, 2, function(ij) {
                i <- ij[1]
                j <- ij[2]
                r <- rank(c(data[group == i], data[group == j]))
                R <- sum(r[1:n.i[i]])
                N <- n.i[i]+n.i[j]
                E <- n.i[i]*(N+1)/2
                V <- n.i[i]*n.i[j]/(N*(N-1))*(sum(r^2)-N*(N+1)^2/4)
                return(abs(R-E)/sqrt(V))
        })
        p <- ptukey(t*sqrt(2), ng, Inf, lower.tail=FALSE)
        result <- cbind(t, p)
        rownames(result) <- combn(ng, 2, paste, collapse=":")
        return(result)
}
Steel.Dwass()を実行する。
#個人の中性脂肪値(TG)を目的変数にして、個人の運動習慣(exercise(3群))で検定を行う。
#この関数では指定するグループは数値データであるべきなので、as.numericで変換している。
Steel.Dwass(sample$TG, as.numeric(sample$exercise))

##            t         p
## 1:2 1.631953 0.2321546
## 1:3 0.502737 0.8699921
## 2:3 1.781280 0.1758568

多重比較の結果、特にどの群間にも差は見られなかった。

NSM3パッケージでsteel-dwass検定をする。

コメントを頂いた内容としては、上記の青木先生の方法は漸近検定(p値の求め方)であるとの指摘を受けた。
一方で、NSM3パッケージでは正確検定も実施できるとのコメントであった。
実際に下記の通り、1)青木先生のコマンド、2)pSDCFlig(Exact)、3)pSDCFlig(Asymptotic)を実行した。

非常に簡単な事例をご紹介する(下記のデータを用意する)
A = c(56, 50, 43, 68, 68, 69, 99, 100)
B = c(1, 1, 1, 2, 2, 2, 3, 3) #グループの数値に変換したもの
1)青木先生のコマンドから...
Steel.Dwass(A, B)
##            t         p
## 1:2 1.992633 0.1139883
## 1:3 1.732051 0.1932252
## 2:3 1.777047 0.1773049
2)pSDCFlig(Exact)
pSDCFlig(A, B, method="Exact")

## Ties are present, so p-values are based on conditional null distribution. 
## Group sizes:  3 3 2 
## Using the Exact method: 
##  
## For treatments 1 - 2 , the  Dwass, Steel, Critchlow-Fligner W  Statistic is 2.818 . 
## The smallest experimentwise error rate leading to rejection is 0.0214 .
##   
## For treatments 1 - 3 , the  Dwass, Steel, Critchlow-Fligner W  Statistic is 2.4495 . 
## The smallest experimentwise error rate leading to rejection is 0.325 .
##   
## For treatments 2 - 3 , the  Dwass, Steel, Critchlow-Fligner W  Statistic is 2.5131 . 
## The smallest experimentwise error rate leading to rejection is 0.1321 .
3)pSDCFlig(Asymptotic)
pSDCFlig(A, B, method="Asymptotic")

## Ties are present, so p-values are based on conditional null distribution. 
## Group sizes:  3 3 2 
## Using the Asymptotic method: 
##  
## For treatments 1 - 2 , the  Dwass, Steel, Critchlow-Fligner W  Statistic is 2.818 . 
## The smallest experimentwise error rate leading to rejection is 0.1141 .
##   
## For treatments 1 - 3 , the  Dwass, Steel, Critchlow-Fligner W  Statistic is 2.4495 . 
## The smallest experimentwise error rate leading to rejection is 0.1935 .
##   
## For treatments 2 - 3 , the  Dwass, Steel, Critchlow-Fligner W  Statistic is 2.5131 . 
## The smallest experimentwise error rate leading to rejection is 0.1775 .
最後にこのページの最初に実行している運動習慣3群ごとの中性脂肪値についての比較をpSDCFlig関数で行う。
pSDCFlig(sample$TG, as.numeric(sample$exercise), method="Exact")
## vector("list", count) でエラー:  ベクトルのサイズは∞ではありえません 

あれ、エラー出るの?(試行数の問題?)

Exactではなく、モンテカルロ(Monte Carlo)なら答えが返ってきます。
(何かご存知の方はコメント宜しくお願いします。)

exercise1 <- sample$TG[as.numeric(sample$exercise)==1]
exercise2 <- sample$TG[as.numeric(sample$exercise)==2]
exercise3 <- sample$TG[as.numeric(sample$exercise)==3]
list <- list(dt1=exercise1, dt2=exercise2, dt3=exercise3)
pSDCFlig(list, method="Monte Carlo")

## Ties are present, so p-values are based on conditional null distribution. 
## Group sizes:  297 1294 1392 
## Using the Monte Carlo(with  10000 Iterations)  method: 
##  
## For treatments 1 - 2 , the  Dwass, Steel, Critchlow-Fligner W  Statistic is -2.3079 . 
## The smallest experimentwise error rate leading to rejection is 0.2213 .
##   
## For treatments 1 - 3 , the  Dwass, Steel, Critchlow-Fligner W  Statistic is -0.711 . 
## The smallest experimentwise error rate leading to rejection is 0.8653 .
##   
## For treatments 2 - 3 , the  Dwass, Steel, Critchlow-Fligner W  Statistic is 2.5191 . 
## The smallest experimentwise error rate leading to rejection is 0.1662 .

以上です。

20160529
RF