目的

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

R の multcomp ライブラリで，glht を使うという手もある（aov も使う）。

使用法

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 値の性質上，むやみに精度を上げても意味はないし，計算時間もかかる。

# aov, glht を使う

> library(multcomp)
> factor.group <- factor(group)
> ans <- glht(aov(data~factor.group), linfct=mcp(factor.group="Dunnett"))
> summary(ans)

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = data ~ factor.group)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)    
2 - 1 == 0   1.1964     0.6453   1.854  0.21560    
3 - 1 == 0   2.9048     0.6937   4.188  < 0.001 ***
4 - 1 == 0   4.7143     0.6665   7.074  < 0.001 ***
5 - 1 == 0   2.7143     0.6665   4.073  0.00115 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported -- single-step method)

 解説ページ