Tukey 検定

テューキーの検定は、各群の母集団が正規分布に従い、かつ、各群の母分散が互いに等しいに、その母平均が帰無仮説に属するかどうかを検定する方法である。テューキーの検定を行う上で、データが次のように配置できるとする。実験群は 1, 2, ..., a まであるとする。

実験群データ個数平均分散
1x11, x12, .. ,x1n1n1μ1σ21
2x21, x22, .. ,x2n2n2μ2σ22
3x31, x32, .. ,x3n3n3μ3σ23
axa1, xa2, .. ,xananaμaσ2a

テューキー検定

表中のデータをもとに、次のようにして誤差自由度 φE、誤差分散 VE を求める。

\[ \phi_{E} = N - a = \sum_{j=1}^{a}n_{j} - a \] \[ V_{E} = \frac{\sum_{j=1}^{a}(n_{j} - 1)\sigma^{2}_{j}}{\phi_{E}} \]

誤差自由度と誤差分散を利用して i 群と j 群の t 統計量を求める。

\[ t_{ij} = \frac{|\mu_{i} - \mu_{j}|}{\sqrt{V_{E}\left( \frac{1}{n_{i}} + \frac{1}{n_{j}} \right)}} \]

すべての組み合わせについて t 統計量を求め、次の表のようにまとめる。

t 値表12...a
1
2t21
:::
ata1ta2...

最後に、すべての t 統計値を境界値と比較する。境界値はステューデント化した分布(Q 分布)を利用して算出する。具体的には、有意水準を α とすると、tij が次式を満たすならば、i 群と j 群の母平均に有意な差がある、と判定を下す。

\[ t_{ij} \ge \frac{Q(a, \phi_{E}, \alpha)}{\sqrt{2}} \]

R を利用したテューキー検定

テューキー検定を行う関数は R の SimComp パッケージに実装されている。ここでサンプルデータを用いてテューキー検定を行う方法を示す。サンプルデータは、19 匹のラットの体重を測定したものである。実験は C 群、G 群、L 群および PF 群の計 4 つの実験群からなる。各実験群には 5 匹のラットを用意したが、L 群では 1 個体が測定できなかったため欠損値が存在する。

library(SimComp)
library(reshape)

d <- read.table("https://stat.biopapyrus.net/data/ratweight.csv", header = TRUE, sep = ",")
d

# データを整形して m に代入
m <- melt(d)
m <- m[(!is.na(m[,2])),]   # 欠損値(NA)を取り除く

head(m)
##   variable value
## 1        C   266
## 2        C   273
## 3        C   260
## 4        C   290
## 5        C   285
## 6        G   208

# Tukey 検定
# grp には m 中のグループの列名(variable)、resp には値の列名(value)を指定
SimTestDiff(data = m, grp = "variable", resp = "value", type = "Tukey", covar.equal = T)
## Test for differences of means of multiple endpoints 
## Assumption: Homogeneous 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      G - C    value      0   -64.40    -9.903      0.0000      0.0000
## 2      L - C    value      0   -21.05    -3.052      0.0081      0.0360
## 3     PF - C    value      0   -93.40   -14.362      0.0000      0.0000
## 4      L - G    value      0    43.35     6.285      0.0000      0.0001
## 5     PF - G    value      0   -29.00    -4.459      0.0005      0.0021
## 6     PF - L    value      0   -72.35   -10.489      0.0000      0.0000

テューキー検定を行うスクリプトを自作するさいに以下のようにする。

tukey.test <- function(data) {
  # 0. この関数で必要な変数を準備
  data.num  <- rep(0, length = ncol(data))  # 各群の個体数
  data.mean <- rep(0, length = ncol(data))  # 各群の平均
  data.var  <- rep(0, length = ncol(data))  # 各群の不偏分散
  phi.E <- 0  # 誤差自由度
  V.E   <- 0  # 誤差分散
  t.values <- p.values <- matrix(0, ncol = ncol(data), nrow = ncol(data))  # t 統計量, p 値
  colnames(t.values) <- rownames(p.values) <- colnames(data)   # 列名をつける
  rownames(t.values) <- rownames(p.values) <- ccolnames(data)   # 行名をつける

  # 1. 各群の個体数、平均、不偏分散を求める。
  data.num  <- apply(!is.na(data), 2, sum)
  data.mean <- apply(data, 2, mean, na.rm = TRUE)
  data.var  <- apply(data, 2, var, na.rm = TRUE)
  
  # 2. 誤差自由度と誤差分散を求める。
  phi.E <- sum(data.num) - ncol(data)
  V.E   <- sum((data.num -1) * data.var) / phi.E

  # 3. t 統計量と p 値を計算する。
  for (i in 1:ncol(data)) {
    for (j in i:ncol(data)) {
      # t 値を計算する
      t.values[i, j] <- abs(data.mean[i] - data.mean[j]) / sqrt(V.E * ((1 / data.num[i]) + (1 / data.num[j])))
      # Q 分布から p 値を計算する
      p.values[i, j] <- ptukey(t.values[i, j] * sqrt(2), ncol(data), phi.E, lower.tail=FALSE)
    }
  }

  # t 値と p 値を返す
  res <- list(t.values = t.values, p.values = p.values)
  return(res)
}


# データを取り込みむ
d <- read.table("ratweight.csv", header=TRUE, sep=",")
 
# テューキー検定を行う
res <- tukey.test(d)


# t 値を確認
res$t.values
##    C        G        L        PF
## C  0 9.902771 3.051734 14.362094
## G  0 0.000000 6.284688  4.459323
## L  0 0.000000 0.000000 10.488978
## PF 0 0.000000 0.000000  0.000000

# p 値を確認
res$p.values
##    C            G            L           PF
## C  1 3.119912e-07 3.626499e-02 1.980294e-09
## G  0 1.000000e+00 7.775854e-05 2.311751e-03
## L  0 0.000000e+00 1.000000e+00 1.462386e-07
## PF 0 0.000000e+00 0.000000e+00 1.000000e+00

t 統計量にルート 2 をかけた数値が、群数(ncol(data) = 4)、自由度(phi.E = 19 - 4 = 15)のステューデント化した範囲の分布(Q 分布)に従う。すなわち、t 統計量にルート 2 をかけた値が、Q(4, 15, α) 以上ならば、i 群と j 群の母平均に危険率 α で差があると判定する。Q 分布表から Q(4, 15, 5%) = 4.1 なので、上記の例では、すべての群の母平均に有意差が認められる。もちろん、p 値もすべて0.05 よりも小さいことが書くにできる。