← 授業トップへ戻る

本日の目標

  1. 多重比較の必要性と Tukey HSD 法の考え方を理解する
  2. Rで多重比較を実行し，どの群間に差があるかを特定できる
  3. 分散分析から多重比較までの結果を統合してレポートにまとめられる
  4. 発展：補正がなぜ全体の誤り率を守れるのかを，式で説明できる

準備：パッケージの読み込み

この回のコードは，すべて tidyverse を読み込んだ状態で動きます。 RStudioを起動したら，まずこれを実行してください。

library(tidyverse)

毎回必要です。読み込まないまま filter() や ggplot() を使うと， 「関数が見つかりません」というエラーになります。

1. 多重比較とは

分散分析の後に必要なこと

分散分析（ANOVA）で「少なくとも1組の群に有意差がある」とわかっても， 「どの群とどの群が異なるのか」はわかりません。

→ 多重比較（post-hoc test） でペアごとの差を検定します。

なぜ専用の方法が必要か？

単純な t検定の繰り返しでは第一種の過誤率が膨らむ（第13回で解説）ため， 検定全体の誤り率を制御する方法を使います。

2. 代表的な多重比較法

方法 特徴 Rの関数
Tukey HSD 全ペア比較，バランスのとれた設計に適す TukeyHSD()
Bonferroni 保守的，検定回数で α を割る pairwise.t.test(..., p.adjust="bonferroni")
Holm Bonferroni より検出力が高い pairwise.t.test(..., p.adjust="holm")
Scheffé 任意の対比を検定，最も保守的 専用パッケージ

本授業では Tukey HSD 法 を中心に扱います。

導出①：ボンフェローニ補正はなぜ効くのか ―― 独立を仮定しない不等式

\(m\) 個の比較を行い，\(i\) 番目の比較で偽陽性を出す事象を \(A_i\) と書きます。守りたいのは「どれか1つでも偽陽性を出す」確率です。

\[\mathrm{FWER} = P(A_1 \cup A_2 \cup \cdots \cup A_m)\]

ここでブールの不等式を使います。和事象の確率は，各事象の確率を単純に足したものを超えません。重なっている部分を二重に数えるぶん，足し算のほうが大きいか等しいからです。

\[P(A_1 \cup \cdots \cup A_m) \le \sum_{i=1}^{m} P(A_i)\]

各比較の水準を \(\alpha/m\) にしておけば，右辺は

\[\sum_{i=1}^{m} \frac{\alpha}{m} = m \cdot \frac{\alpha}{m} = \alpha\]

となり，\(\mathrm{FWER} \le \alpha\) が保証されます。

第13回の式との違いに注目してください。 第13回の \(1-(1-\alpha)^k\) は，\(k\) 回の検定が独立だと仮定した等式でした。一方こちらは不等式で，検定どうしがどう絡み合っていても成り立ちます。総当たり比較のように同じ群のデータを共有していても構いません。何も仮定せずに守れるのが，ボンフェローニの強みです。

代償は保守的になることです。不等式は「最悪でもこれ以下」としか言っていないので，実際の誤り率は \(\alpha\) よりかなり小さくなりがちです。そのぶん，本当にある差を見逃しやすくなります。

数値で確かめる：3群なら比較は \(m = 3\) 組なので，各比較を \(0.05/3 = 0.0167\) で判定することになります。iris の versicolor 対 virginica では，無補正の \(p = 2.77 \times 10^{-9}\) が補正後は \(8.30 \times 10^{-9}\) になります。ちょうど3倍です。

使った数学：ブールの不等式（和事象の確率は各確率の和以下）。独立性はどこにも使っていません。

導出②：Tukey HSD は「いちばん離れた2つ」を基準にする

ボンフェローニは，比較どうしの関係を一切使わずに安全側で抑える方法でした。Tukey HSD は別の道を通ります。

\(k\) 個の群があり，すべてのペアを比べるとします。ここで考えてみてください。どれか1つでも偽陽性を出してしまうのは，どういう場合でしょうか。それは，ペアの中でいちばん差が大きく出たものが，しきい値を超えてしまう場合です。

つまり守るべき相手は，個々のペアではなく，最大の平均と最小の平均の差ただ1つです。

\[\max_j \bar{x}_j - \min_j \bar{x}_j\]

この「最大と最小の幅」が，帰無仮説のもとでどれくらいの大きさになるか。その分布をスチューデント化された範囲（studentized range，記号 \(q\)）と呼びます。\(q\) の上位 \(5\%\) 点をしきい値にすれば，最大の差が守られ，したがってすべてのペアが同時に守られます。

\[HSD = q_{\alpha}(k,\ df) \times \sqrt{\frac{MS_{群内}}{n}}\]

数値で確かめる：iris は \(k=3\)，\(df=147\)，\(MS_{群内}=0.2650\)，各群 \(n=50\) です。\(q_{0.05}(3,\ 147) = 3.3484\) なので

\[HSD = 3.3484 \times \sqrt{0.2650/50} = 0.2438\]

平均差がこれを超えれば有意です。実際の3つの平均差は \(0.930\)，\(1.582\)，\(0.652\) で，いずれも大きく上回っています。

同じ条件でボンフェローニの最小有意差を計算すると \(0.2493\) となり，Tukey の \(0.2438\) よりわずかに大きい。つまり Tukey のほうが，小さい差でも有意と判定できます。臨界値を \(t\) の尺度に揃えて比べると，ボンフェローニが \(2.4217\)，Tukey が \(2.3677\) です。

Tukey のほうが検出力が高いのは，比較どうしが無関係ではないことを利用しているからです。 ボンフェローニは「どんな絡み方をしていても大丈夫なように」と身構えて余裕を取りますが，Tukey は総当たり比較という具体的な状況での \(q\) の分布を使うので，その余裕を削れます。

群の数が増えるほど，この差は広がります。

群数 \(k\) 比較数 \(m\) ボンフェローニ Tukey
3 3 2.4217 2.3677
5 10 2.8500 2.7619
10 45 3.3265 3.2129

上の表の臨界値は，いずれも自由度 \(147\)（iris と同じ）で計算したものです。

この \(HSD\) の式は，各群の \(n\) が等しい（バランス型）ことを使っています。 分母の \(\sqrt{MS_{群内}/n}\) に \(n\) が1つしか出てこないのは，そのためです。 群のサイズが違う場合，Rの TukeyHSD() は自動的に \(\sqrt{\dfrac{MS_{群内}}{2}\left(\dfrac{1}{n_i}+\dfrac{1}{n_j}\right)}\) に置き換えて計算します （Tukey–Kramer法）。演習3の mtcars は各群 11・7・14 名でサイズが違うので，こちらが使われます。

使った数学：最大値と最小値の差の分布という考え方。なお \(q\) 分布そのものを導くのは，この科目の範囲外です。

導出③：ホルム法 ―― 段階的に条件を緩める

表に挙げたホルム法は，ボンフェローニを改良したものです。考え方は単純で，\(p\) 値を小さい順に並べ，しきい値を段階的に緩めていきます。

\(m\) 個の \(p\) 値を小さい順に \(p_{(1)} \le p_{(2)} \le \cdots \le p_{(m)}\) と並べ，

\[p_{(1)} \text{ は } \frac{\alpha}{m}，\quad p_{(2)} \text{ は } \frac{\alpha}{m-1}，\quad \dots，\quad p_{(i)} \text{ は } \frac{\alpha}{m-i+1}\]

と比べます。順に判定していき，初めて基準を満たさなかったところで止めて，それ以降はすべて有意ではないとします。

いちばん小さい \(p\) 値に対する基準 \(\alpha/m\) はボンフェローニと同じですが，2番目以降は分母が小さくなるぶん基準が緩みます。したがって，ボンフェローニで有意になるものはホルムでも必ず有意になり，その逆は成り立ちません。同じ保証のもとで検出力だけが高いので，2つのうちならホルムを選んで困ることはありません。

数値で確かめる：iris の3つの比較のうち，versicolor 対 virginica は無補正の \(p = 2.77 \times 10^{-9}\) が最大の \(p\) 値です。ボンフェローニでは一律に3倍されて \(8.30 \times 10^{-9}\) になりますが，ホルムでは分母が \(m-i+1 = 1\) となるため \(2.77 \times 10^{-9}\) のまま据え置かれます。

使った数学：\(p\) 値を大きさの順に並べることと，段階的な判定。

3. Tukey HSD 法の実行

例：iris データ（3種のがく片の長さ）

# 分散分析の実行（第13回の復習）
result_anova <- aov(Sepal.Length ~ Species, data = iris)
summary(result_anova)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## Species       2  63.21  31.606   119.3 <2e-16 ***
## Residuals   147  38.96   0.265                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Tukey HSD 多重比較
result_tukey <- TukeyHSD(result_anova)
result_tukey
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Sepal.Length ~ Species, data = iris)
## 
## $Species
##                       diff       lwr       upr p adj
## versicolor-setosa    0.930 0.6862273 1.1737727     0
## virginica-setosa     1.582 1.3382273 1.8257727     0
## virginica-versicolor 0.652 0.4082273 0.8957727     0

結果の読み方

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Sepal.Length ~ Species, data = iris)

$Species
                      diff       lwr       upr   p adj
versicolor-setosa    0.930 0.6862273 1.1737727 0.0e+00
virginica-setosa     1.582 1.3382273 1.8257727 0.0e+00
virginica-versicolor 0.652 0.4082273 0.8957727 8.3e-09
列 意味
diff 平均差（後ろの群 − 前の群）
lwr / upr 平均差の95%信頼区間の下限・上限（ファミリー全体で95%）
p adj 調整済みp値（< .05 で有意）

p adj の表示に注意。 上2つは 0.0e+00 と出ていますが， 「ちょうど0」ではなく表示できないほど小さいという意味です。 3つ目の 8.3e-09 は \(8.3 \times 10^{-9}\) で，これも十分小さい値です。 レポートには p adj = 0 ではなく p < .001 と書いてください。

diff の符号は行の名前が決めます。 versicolor-setosa は「versicolor の平均 − setosa の平均」です。 \(0.930 > 0\) なので versicolor のほうが長い，と読みます。

信頼区間が「95%」である意味も，第10回とは少し違います。 ここでの95%は3つの区間をまとめての信頼度（family-wise）です。 3本の区間が同時にすべて正しく真の差を含む確率が95%という意味で， そのぶん1本ずつの区間は，個別に作るより広くなっています。

4. Tukey HSD 結果の可視化

# confint() でデータフレームに変換して描く
tukey_df <- as.data.frame(result_tukey$Species)
tukey_df$comparison <- rownames(tukey_df)

ggplot(tukey_df, aes(x = comparison, y = diff)) +
  geom_point(size = 3, color = "steelblue") +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1, color = "steelblue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Tukey HSD：群間の平均差と95%CI",
       x = "比較ペア", y = "平均差（cm）") +
  coord_flip() +
  theme_bw()

信頼区間が 0 をまたがっていなければ有意差あり。

5. 全体の結果をグラフで示す

有意差ブラケット付きの箱ひげ図

# ggplot2で箱ひげ図を作成
p_box <- ggplot(iris, aes(x = Species, y = Sepal.Length, fill = Species)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.7) +
  geom_jitter(width = 0.1, alpha = 0.3, size = 1.5) +
  scale_fill_brewer(palette = "Set2") +
  labs(title = "種ごとのがく片の長さ（Tukey HSD: すべての群間に p < .001）",
       x = "種", y = "がく片の長さ (cm)") +
  theme_bw() +
  theme(legend.position = "none")
p_box

6. 完全な分析フローのまとめ

Step 1：記述統計

iris |>
  group_by(Species) |>
  summarise(n = n(), M = mean(Sepal.Length), SD = sd(Sepal.Length)) |>
  mutate(across(where(is.numeric), ~ round(., 2)))

Step 2：前提条件の確認

bartlett.test(Sepal.Length ~ Species, data = iris)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Sepal.Length by Species
## Bartlett's K-squared = 16.006, df = 2, p-value = 0.0003345

第13回で見たとおり，このデータは p = .0003 で等分散が成り立っていません。 それでも aov() → TukeyHSD() で進めるのは， 群ごとの \(n\) が 50 ずつ等しく，分散のズレの影響が比較的小さいためです。 \(n\) が群によって大きく違い，かつ分散も違う場合は，この流れは安全ではありません。 その場合は oneway.test() と，等分散を仮定しない多重比較 （pairwise.t.test(..., pool.sd = FALSE, p.adjust.method = "holm")）を使ってください。

# 等分散を仮定しない多重比較（参考）
pairwise.t.test(iris$Sepal.Length, iris$Species,
                pool.sd = FALSE, p.adjust.method = "holm")
## 
##  Pairwise comparisons using t tests with non-pooled SD 
## 
## data:  iris$Sepal.Length and iris$Species 
## 
##            setosa  versicolor
## versicolor < 2e-16 -         
## virginica  < 2e-16 1.9e-07   
## 
## P value adjustment method: holm

このデータでは，どちらで計算しても3組すべて有意という結論は変わりません。

Step 3：分散分析

result_anova <- aov(Sepal.Length ~ Species, data = iris)
summary(result_anova)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## Species       2  63.21  31.606   119.3 <2e-16 ***
## Residuals   147  38.96   0.265                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Step 4：多重比較

TukeyHSD(result_anova)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Sepal.Length ~ Species, data = iris)
## 
## $Species
##                       diff       lwr       upr p adj
## versicolor-setosa    0.930 0.6862273 1.1737727     0
## virginica-setosa     1.582 1.3382273 1.8257727     0
## virginica-versicolor 0.652 0.4082273 0.8957727     0

Step 5：効果量

ss <- anova(result_anova)[["Sum Sq"]]   # 1番目が群間，2番目が群内
round(ss[1] / sum(ss), 3)               # η²
## [1] 0.619

7. 結果の報告（論文・レポート形式）

方法：3種のアヤメ（setosa, versicolor, virginica）のがく片の長さを 一元配置分散分析によって比較した。多重比較には Tukey HSD 法を用いた。

結果：種の主効果が有意であった，F(2, 147) = 119.3, p < .001, η² = .62。 Tukey HSD 法による多重比較の結果， setosa と versicolor（p < .001），setosa と virginica（p < .001）， versicolor と virginica（p < .001）の間でそれぞれ有意差が認められた。

演習

演習1：多重比較の実行

iris データで Petal.Length の一元配置分散分析を実行し， Tukey HSD で多重比較を行ってください。

どの群間に有意差がありましたか？

result <- aov(___ ~ ___, data = iris)
TukeyHSD(result)

演習2：視覚化

演習1の Tukey HSD の結果を，信頼区間のプロットで表示してください。

tukey_df は5節で作った変数です。5節のコードを実行していなければ存在しません。 演習1から始めた場合は，次の2行で作り直してください。

tukey_df <- as.data.frame(TukeyHSD(result)$Species)
tukey_df$comparison <- rownames(tukey_df)

演習3：総合レポート（最終課題の練習）

以下の問いに対して，Step 1〜5 の分析フローを実行し， 結果を論文の「方法・結果」の節として1段落にまとめてください。

問い：mtcars データで，シリンダー数（cyl：4・6・8）によって 馬力（hp）に違いがあるかを検討してください。

  1. 記述統計（平均・SD・n）
  2. 箱ひげ図
  3. 一元配置分散分析（aov()）
  4. Tukey HSD
  5. 効果量 η²
  6. 結果の記述

⚠ cyl は数値のままでは分散分析になりません

mtcars$cyl は 4・6・8 という数値の列です。そのまま aov(hp ~ cyl, data = mtcars) と書くと，Rは「4・6・8 という量」として扱い， 自由度1の回帰を当てはめてしまいます（群の比較になりません）。 続けて TukeyHSD() を実行すると 「当てはめモデル中には因子がありません」というエラーになります。

先に因子（factor）に変換してください。

mtcars2 <- mtcars |> mutate(cyl_f = factor(cyl))
result  <- aov(hp ~ cyl_f, data = mtcars2)
summary(result)
TukeyHSD(result)

summary(result) の Df が 2（\(= k-1\)）になっていれば，正しく群比較になっています。 1 になっていたら，factor() を忘れています。

ヒント：3組のうち1組だけ有意になりません。 「すべて有意」ではない結果をどう記述するかが，この演習の狙いです。 また各群の \(n\) が 11・7・14 と違うので，導出②の Tukey–Kramer法 が使われます。

本日のまとめ

  • 分散分析で有意差が出た後，Tukey HSD でどの群間かを特定する
  • TukeyHSD() で調整済みp値と信頼区間を確認
  • 分析フローは「記述統計 → 前提確認 → ANOVA → 多重比較 → 効果量 → 報告」
  • 結果の記述には F値・df・p値・η²・多重比較の結果を含める

宿題・予習

  • 演習3（総合問題）を完成させ，スクリプトを保存する
  • 次回（第15回）は最終まとめと習熟度確認テスト（第3回）です
  • 第1〜14回の内容全体を見直しておいてください

← 授業トップへ戻る