目的

ダネットの方法による多重比較を行う。

R の multcomp ライブラリで，simtest(formula, type="Dunnett") を使うという手もある。
同時に，simtest(formula, type="Williams") もチェックのこと。
注：multcomp の中には simtest はなくなった（少なくとも R 2.10.0 では）

使用法

dunnett(data, group)

引数

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

ソース

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

# ダネットの方法による多重比較
dunnett <- 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)
        }

        pdunnett <- function(x, a, df, r)                       # P 値を計算する
        {
                corr <- diag(a-1)
                corr[lower.tri(corr)] <- r
                1-pmvt(lower=-x, upper=x, delta=numeric(a-1), df=df, corr=corr, abseps=0.0001)
        }

        OK <- complete.cases(data, group)                       # 欠損値を持つケースを除く
        data <- data[OK]
        group <- factor(group[OK])                              # 群変数は factor に変換
        ni <- table(group)                                      # 各群のデータ数
        a <- length(ni)                                         # 群の数
        n <- length(data)                                       # 全体のデータ数
        mi <- tapply(data, group, mean)                         # 各群の平均値
        vi <- tapply(data, group, var)                          # 各群の不偏分散
        Vw <- sum(vi*(ni-1))/(n-a)                              # 群内分散
        rho <- get.rho(ni)                                      # ρ
        t <- (abs(mi-mi[1])/sqrt(Vw*(1/ni+1/ni[1])))[2:a]       # 対照群と各群の比較における t 値
        p <- sapply(t, pdunnett, a, n-a, rho)                   # P 値
        result <- cbind(t, p)
        rownames(result) <- paste(1, 2:a, sep=":")
        return(result)
}


使用例

data <- c(
    14, 15, 14, 16, 15, 17, 17,        # 第 1 群（対照群）のデータ，7 例
    17, 16, 17, 16, 15, 18, 19, 15,    # 第 2 群（処理群）のデータ，8 例
    18, 19, 20, 19, 17, 17,            # 第 3 群（処理群）のデータ，6 例
    20, 21, 19, 20, 19, 22, 20,        # 第 4 群（処理群）のデータ，7 例
    19, 20, 19, 17, 17, 17, 18         # 第 5 群（処理群）のデータ，7 例
)
group <- rep(1:5, c(7, 8, 6, 7, 7))     # 群を表す数値

library(mvtnorm)
dunnett(data, group)

出力結果例

> library(mvtnorm)	# mvtnorm を利用するために，R を起動するごとに必要
> dunnett(data, group)
           t            p
1:2 1.854090 2.156676e-01
1:3 4.187543 8.023633e-04
1:4 7.073685 5.706414e-08
1:5 4.072727 1.169213e-03

注：この dunnett 関数は，P 値を求めるために，pdnunnett 関数の中で pmvt 関数を呼んでいる。
pmvt 関数は乱数を発生させて P 値を求めているため，毎回の P 値は変わる。
大幅に変わるわけではなく，その精度も制御できる（現在は bseps=0.0001）としているが，
この数値を小さくすれば精度は高くなる。
しかし，P 値の性質上，むやみに精度を上げても意味はないし，計算時間もかかる。

multcomp の中には simtest はなくなった（少なくとも R 2.10.0 では）

SimComp の中の SimTestDiff が該当するのかな（シミュレーションによる検定なので，上の dunnett 関数の結果とはずいぶん違う）

> data <- c(
+     14, 15, 14, 16, 15, 17, 17,        # 第 1 群（対照群）のデータ，7 例
+     17, 16, 17, 16, 15, 18, 19, 15,    # 第 2 群（処理群）のデータ，8 例
+     18, 19, 20, 19, 17, 17,            # 第 3 群（処理群）のデータ，6 例
+     20, 21, 19, 20, 19, 22, 20,        # 第 4 群（処理群）のデータ，7 例
+     19, 20, 19, 17, 17, 17, 18         # 第 5 群（処理群）のデータ，7 例
+ )
> group <- rep(1:5, c(7, 8, 6, 7, 7))     # 群を表す数値
> d <- data.frame(data=data, group=factor(group))
> library(SimComp)
> SimTestDiff(d, "group", "data", type="Dunnett")
 
Test for differences of means of multiple endpoints 
Assumption: Heterogeneous covariance matrices for the groups 
Alternative hypotheses: True differences not equal to the margins 
 
  comparison endpoint margin estimate statistic p.value.raw p.value.adj
1      2 - 1     data      0    1.196     1.729      0.1076      0.2927
2      3 - 1     data      0    2.905     4.211      0.0015      0.0057
3      4 - 1     data      0    4.714     7.505      0.0000      0.0000
4      5 - 1     data      0    2.714     4.082      0.0015      0.0061
 
> SimTestDiff(d, "group", "data", type="Williams")
 
Test for differences of means of multiple endpoints 
Assumption: Heterogeneous covariance matrices for the groups 
Alternative hypotheses: True differences not equal to the margins 
 
  comparison endpoint margin estimate statistic p.value.raw p.value.adj
1        C 1     data      0    2.714     4.082      0.0015      0.0040
2        C 2     data      0    3.714     6.517      0.0000      0.0001
3        C 3     data      0    3.471     6.347      0.0001      0.0001
4        C 4     data      0    2.821     5.275      0.0005      0.0011

 解説ページ