R の分散分析: aov()、分散分析表、Tukey HSD

なぜ t 検定の山ではなく分散分析なのか

t 検定は 2 つの平均を比較します。3 群以上あるとき、すべての組で t 検定をしたくなりますが、それこそが罠です。0.05 の水準で実行する各検定は 5% の偽陽性リスクを負い、そのリスクは蓄積します。4 群なら 6 回の対比較となり、少なくとも 1 つが偽の「有意」になる確率はおよそ 26%、5 群（10 回の検定）ならおよそ 40% です。純粋なノイズから発見をでっち上げてしまうことになります。

分散分析 (ANOVA、ANalysis Of VAriance) は、1 つ の p 値で 1 つ の問いを立てることでこれを解決します。すべての群平均は同じか、それとも少なくとも 1 つは異なるか? 名前に反して比較するのは平均です。ただし分散を分析することでそれを行います。群の平均が、群内のノイズで説明できる以上に散らばっていれば、何か本物のことが起きているのです。

aov() による一元配置分散分析

PlantGrowth はこのために作られたようなデータです。対照条件と 2 つの処理条件のもとでの植物の重量が入っています。aov() で当てはめ、そして重要なことに、summary() で表を出力します。

fit <- aov(weight ~ group, data = PlantGrowth)
summary(fit)

formula は「weight は group に依存するか?」と読みます。グループ化変数は因子でなければなりません。PlantGrowth$group はすでに因子ですが、群が数値（用量、バッチ ID）として符号化されている場合は包みましょう: aov(y ~ factor(dose), ...)。そうしないと aov() は群の平均を比較する代わりに、群のコードを通る 回帰直線 を黙って当てはめます。誤ったモデルなのに、エラーメッセージは出ません。

分散分析表を読む

表は 2 行 — 因子と残差 — と 5 列からなります。セルごとに見ていきましょう。

            Df Sum Sq Mean Sq F value Pr(>F)
group        2  3.766  1.8832   4.846 0.0159
Residuals   27 10.492  0.3886

• Df — 自由度です。因子の行は 群の数 − 1（3 群 → 2）、残差の行は 観測数 − 群の数（30 − 3 = 27）になります。R が意図どおりの計画を認識したかを手早く確認できます。
• Sum Sq — 変動を 2 つの山に分けたものです。group の行は群間の山で、群の平均が全体平均からどれだけ離れているかを表します。Residuals は群内の山で、植物が 自分の 群の平均のまわりでどれだけばらつくかを表します。両者を足すとデータ全体の変動になります。
• Mean Sq — 各 Sum Sq をその Df で割ったもので、変動の山を自由度あたりの比較可能な率に変えます。残差の Mean Sq (0.389) がノイズの水準です。
• F value — その比です。group の Mean Sq を Residuals の Mean Sq で割ります (1.8832 / 0.3886 ≈ 4.85)。すべての群平均が本当に等しければ、この比は 1 のあたりをうろつきます。F が大きくなるほど、群間の差はノイズで説明できる範囲を上回っていきます。
• Pr(>F) — p 値です。3 つの真の平均がすべて同一だったとして、これほど大きな F が生じる確率を表します。ここでは 0.016 で、慣例的な 0.05 の水準では「すべて等しい」を棄却するのに十分小さい値です。いつものように、これは帰無仮説が正しい確率ではありませんし、どの群がどれだけ異なるかについては何も語りません。

最後の点が決定的な限界です。有意な F は「どこかに何らかの差がある」と言っているだけです。 それ以上ではありません。

どの群が異なるのか? TukeyHSD()

続く問いには事後検定が必要です。Tukey の Honest Significant Difference は、すべての比較を通じた ファミリーワイズ の誤り率を 5% に保ちながら、すべての組を検定します。

fit <- aov(weight ~ group, data = PlantGrowth)
TukeyHSD(fit)

組ごとに 1 行、1 行につき 4 つの数値が並びます。

• diff — 推定された平均の差（後に挙げた群から先の群を引いた値）。
• lwr, upr — その差に対する 95% のファミリーワイズ信頼区間。
• p adj — 3 回の比較を行うことをすでに調整済みの p 値。

PlantGrowth では、trt2-trt1 の差が約 0.87、p adj ≈ 0.012、区間はゼロを含みません。処理 2 は処理 1 より成長がよいということです。処理対対照の 2 行 (trt1-ctrl、trt2-ctrl) はどちらも区間がゼロをまたぎ、p adj は 0.05 をかなり上回ります。この標本では、どちらの処理も対照と区別できません。経験則はどこでの信頼区間とも同じです。区間がゼロを含まない ⇔ p adj が 0.05 未満。

全体の F は「差が存在する」と言い、Tukey がそれをちょうど 1 つの組に特定していることに注目してください。この 2 段階の構造こそが設計そのものです。まず 1 つの誠実な全体検定、次に適切に補正された対ごとの探索です。

二元配置分散分析: 2 つの因子とその交互作用

グループ化変数が 2 つあるとき、1 回の呼び出しで両方を検定し、さらに両者が 交互作用 するかどうかも検定できます。ToothGrowth はサプリメントの種類 (supp) と用量（0.5、1、2 mg — 数値なので factor() が必要）を交差させています。

fit2 <- aov(len ~ supp * factor(dose), data = ToothGrowth)
summary(fit2)

supp * factor(dose) は 3 つの効果に展開され、それぞれが表の 1 行になります。

• supp — 用量を通じて平均したとき、サプリメントの種類は効くか?（効きます: p ≈ 0.0002。）
• factor(dose) — サプリメントを通じて平均したとき、用量は効くか?（明白に効きます: p はごく小さい値です。）
• supp:factor(dose) — 交互作用です。サプリメントの効果 は 用量に応じて 変化するか? ここでは p ≈ 0.022 で、変化します。低用量ではオレンジジュースがアスコルビン酸を上回りますが、2 mg の用量では差が縮まります。

有意な交互作用は主効果に対する注意書きです。「サプリメントの種類が効く」は 平均的には 真であるにすぎず、その平均は用量依存の物語を隠しています。交互作用が有意なときは、主効果だけを報告するのではなく、組み合わせを記述しましょう（TukeyHSD(fit2, "supp:factor(dose)") や、グループごとの集約による群平均の表など）。意図的に交互作用なしのモデルを使いたいときにだけ、* の代わりに + を使います。

前提 — そして順位ベースの代替手段

分散分析の数学は 3 つの前提に依拠します。妥協の余地が小さい順に並べます。

• 独立性 — 観測値が互いに影響しないこと。違反を後から修正する手立てはありません。研究計画の性質です。
• 群間で等しい分散 — 残差の Mean Sq は 1 つにプールされたノイズの推定値なので、群のノイズはおおむね等しいべきです。1 行での確認: bartlett.test(weight ~ group, data = PlantGrowth)（p 値が大きければ、分散が等しくないという証拠はないということです）。
• 残差がおおむね正規 — 均衡した計画とそれなりの群サイズなら最も影響が小さい前提です。residuals(fit) を使って残差のヒストグラムや箱ひげ図を目視しましょう。

データが大きく歪んでいる、順序尺度である、外れ値だらけであるといった場合、一元配置分散分析に対する順位ベースの代替手段は kruskal.test(weight ~ group, data = PlantGrowth) です。wilcox.test() の多群版にあたります。1 行で書け、いくらかの検出力と引き換えに頑健性を得られます。

この章のまとめ

• 分散分析は 3 群以上の平均について 1 つの p 値で 1 つの問いを立てます。t 検定の多重性に対する解決策です。
• fit <- aov(y ~ group, data = df); summary(fit) を使い、グループ化変数は因子でなければなりません。
• F 値は群間のシグナルを群内のノイズで割った値です。Pr(>F) が小さいということは「どこかに何らかの差がある」という意味であり、それ以上ではありません。
• TukeyHSD(fit) が どの 組が異なるかを、多重比較の補正込みで見つけます。
• a * b は 2 つの因子とその交互作用を当てはめます。交互作用が有意なら、主効果だけでは物語を語れません。
• 等分散 (bartlett.test) と残差の正規性を確認しましょう。kruskal.test() が順位ベースの代替手段です。

次は、群の平均の比較から関係のモデル化へ — lm() による線形回帰です。