ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【Rで統計】R で ANOVA(分散分析)への最短ルート

ANOVA(分散分析)は、複数グループの平均に差があるかどうかを一度に判定するための手法です。本来は前提条件や理論的背景を理解した上で使うべき解析ですが、ここではそれを一切扱いません。
 この解説の目的はただ一つです。「最短手順で、解析結果を出すこと
 私は理解を大切にする研究者ですが、最近の学生はその重要性をなかなか理解してくれません。ということで、「理解しなくても、迷わず進めば結果が出る」説明をしてみます。

1. データの準備

 ANOVA では、最終的に条件の違いなどを区別する「グループ」と、観測結果の「数値」が対応していれば適用できます。まず、RでANOVAを試すために、R で、そのようなデータを人工的に作ってみます。

1.1 ベクトルで作る方法:1因子の分散分析(1元配置 ANOVA)

 ここでは、因子が 1 つだけの ANOVA(1元配置分散分析)を想定してデータを作ります。つまり、調べたい要因は グループ(group)ただ 1 つであり、その水準(A、B、C など)によって数値データの平均が異なるかを調べます。
 この場合、データの作り方は非常に単純で、「すべての測定値を 1 本の数値ベクトルにまとめ、その値がどのグループに属するかを別の変数で指定する」だけです。

# 数値データ(1元配列)
y <- c(
  rnorm(10, 0.0, 1),
  rnorm(10, 0.5, 1),
  rnorm(10, 1.0, 1)
)

# 因子(1因子:グループ)
group <- factor(rep(c("A", "B", "C"), each = 10))

# data.frame にまとめる
dat <- data.frame(group, y)

ここで用いている因子は group の 1 つだけです。y はすべての測定値を並べた数値ベクトルであり、group は各値がどの水準(A、B、C)に属するかを示す 1因子の factor です。1因子 ANOVA に必要なのは、この 2 つだけです。この形になっていれば、1元配置分散分析はそのまま実行できます。

aov(y ~ group, data = dat)

余分な変数や構造を加える必要はありません。因子が 1 つである限り、この形式が最短で確実な準備方法です。

1.2 ベクトルで作る方法:2因子の分散分析(2元配置 ANOVA)

 次に、因子が 2 つある ANOVA(2元配置分散分析)を想定してデータを作ります。つまり、数値データに影響すると考える要因が 2 種類あり、それぞれの因子と、その組み合わせによって平均が異なるかを調べます。ただし、データの作り方自体は 1 因子の場合とほとんど変わりません。違いは、factor が 2 本になるという点だけです。

# 各条件のサンプル数
n_each <- 10

# 因子1(例:グループ)
group <- factor(rep(c("A", "B"), each = n_each * 2))

# 因子2(例:条件)
condition <- factor(rep(c("X", "Y"), each = n_each, times = 2))

# 数値データ
y <- c(
  rnorm(n_each, 0.0, 1),  # A-X
  rnorm(n_each, 0.5, 1),  # A-Y
  rnorm(n_each, 0.3, 1),  # B-X
  rnorm(n_each, 0.8, 1)   # B-Y
)

# data.frame にまとめる
dat <- data.frame(group, condition, y)

ここで用いている因子は、

  • group
  • condition

2 つです。y はすべての測定値を 1 本の数値ベクトルとして並べ、各値がどの 因子の組み合わせ(例:A–X、B–Y)に対応するかを、2 本の factor で指定しています。2因子 ANOVA に必要なのは、数値ベクトル 1 本と factor 2 本だけです。 この形になっていれば、2元配置分散分析はそのまま実行できます。

aov(y ~ group * condition, data = dat)

1因子の場合と同様に、余分な構造や複雑なデータ形式は必要ありません。 因子が増えても、「数値 + 因子」の基本構造は変わらない、ということを覚えておいてください。

2. 最短とはいえデータの構造を必ず目で見ろ

 最短手順で ANOVA を実行する場合であっても、解析の前に一度はデータを可視化してください。ANOVA は数値だけを入力すれば結果が出ますが、データの形そのものが破綻している場合でも計算自体は止まらないからです。
 ここでは、各グループの分布を一目で確認できる箱ひげ図を用います。

boxplot(y ~ group, data = dat,
        col  = c("skyblue", "lightgreen", "salmon"),
        xlab = "Group",
        ylab = "y",
        main = "Group-wise distribution")

この図を見る目的は、統計的な検定条件を厳密に確認することではありません。 「このデータで平均を比べること自体が不自然ではないか」を直感的に判断するためです。
 具体的には、グループ間で平均がある程度ずれていそうかどうか、ばらつきが極端に異なっていないか、あるいは一部の外れ値が全体を支配していないか、といった点を確認します。これらに明らかな違和感がある場合、ANOVA の結果は数値としては出ても、解釈に意味がなくなります。ANOVA は平均の差を検定する手法です。そのため、計算を始める前に、「そもそも平均を比較する状況として妥当かどうか」を図で確認することが、最短手順であっても省略してはいけない唯一の工程になります。

可視化の例

 1 因子の場合、プロットの目的は単純で、グループごとの分布を並べて表示し、平均がずれていそうかを確認することだけでした。そのため、横軸に因子を置き、箱ひげ図を並べるだけで十分でした。一方、2 因子の場合は、確認すべき点が 1 つ増えます。平均の違いを生み出している可能性のある要因が 2 種類あるため、「それぞれの因子の違い」と「因子どうしの組み合わせによる違い」を同時に意識する必要があります。ただし、やることが本質的に難しくなるわけではありません。
 実務的には、1 因子の箱ひげ図をそのまま拡張すると考えれば十分です。片方の因子を横軸に使い、もう片方の因子を色分けやパネル分けとして表現します。これにより、1 因子では見えなかった「同じグループ内で条件が変わったときの違い」や、「条件によってグループ差の出方が変わっていないか」といった点を視覚的に確認できます。
 つまり、1 因子では「平均との差がありそうか」を見るだけでしたが、2 因子ではそれに加えて、「その差がもう一つの因子によって変わっていないか」を見ることになります。プロットの作り方は少し増えますが、考え方は 1 因子の延長線上にあり、特別な図を新しく覚える必要はありません。

3. ANOVA を実行するコマンド

 ここでは、何も難しいことは考えません。理屈や前提条件は後回しにして、とりあえず ANOVA を実行してみます
 R では、ANOVA は 1 行のコマンドで実行できます。すでに用意したデータ dat に対して、次のように書くだけです。

fit <- aov(y ~ group, data = dat)
summary(fit)

このコードの意味を細かく理解する必要はありません。やっていることは単純で、y という数値データが group という因子によって違うかどうかを、R に判定させているだけです。    

 aov() は分散分析を行うための関数で、その結果を fit という名前で保存しています。続く summary() は、その結果を 人が読める形で表示するためのコマンドです。この 2 行は常にセットで使うと考えて構いません。
 これで、あなたも「ANOVA使い」になれました。正しいかどうかを深く考える前に、まずコマンドを打って出力を確認します。ANOVA は、この 2 行を書けるようになれば、とりあえず回せる解析です。
 因子が 2 つある場合でも、やることはほとんど変わりません。違うのは、式の右辺に因子が 2 つ並ぶ、という点だけです。  すでに用意した 2 因子データ dat に対しては、次のように書きます。

fit <- aov(y ~ group * condition, data = dat)
summary(fit)

ここでは、y という数値データが、groupcondition という 2 つの因子、およびその組み合わせによって違うかどうかを、R にまとめて判定させているだけです。
 group * condition と書くことで、

  • group の効果
  • condition の効果
  • 両者の組み合わせの効果

一度に調べる形になりますが、現時点では「2 因子の場合はこう書く」ことが多いということです。

 1 因子のときと同様に、

fit <- aov(...)
summary(fit)

という 2 行セットは変わりません。

4. 出力結果の見方

■ 1因子 ANOVA

 上で紹介したANOVAのコマンドを実行すると、次のような表が表示されます。

> summary(fit)
            Df Sum Sq Mean Sq F value Pr(>F)
group        2  2.547   1.274   1.176  0.324
Residuals   27 29.243   1.083

1 因子 ANOVA の場合、注目するのは 1 行、1 つの数値だけです。それは、group 行の Pr(>F) です。
 この例では、

Pr(>F) = 0.324

となっています。これは、有意水準 5%(0.05)で判断すると、グループ間の平均に有意な差は認められないことを意味します。つまり、このデータからは、「A・B・C のいずれかの平均が他と異なる」とは言えません。
 他の列に並んでいる DfSum SqMean SqF value は、ANOVA の計算途中で使われる量ですが、これについては以前の解説を参考にしてください。今回は、最短、最低限です。
 まとめると、この結果から、最短、最低限の説明をするとすれば、

1 因子 ANOVA の結果、グループによる平均の差は統計的に有意ではなかった(p = 0.324)。

ということです。

■ 2因子 ANOVA

 2 因子 ANOVA を実行すると、次のような表が表示されます。

> summary(fit)
                Df Sum Sq Mean Sq F value Pr(>F)
group            1   0.00   0.000   0.000  0.996
condition        1   2.07   2.074   1.427  0.240
group:condition  1   1.72   1.716   1.181  0.284
Residuals       36  52.31   1.453

2 因子の場合でも、見るべき場所は決まっています。それぞれの因子とその組み合わせに対応する Pr(>F) の列です。
 この表には、判断対象が 3 つ並んでいます。
 まず、group 行は、group という因子そのものの効果を表します。この例では p 値が 0.996 であり、有意水準 5% では group による平均の差は認められません。
 次に、condition 行は、condition という因子単独の効果を表します。ここでは p 値が 0.240 となっており、condition による平均の差も有意とは言えません。  最後に、group:condition 行は、2 つの因子の組み合わせによる効果、いわゆる交互作用を表します。この例では p 値が 0.284 であり、group と condition の組み合わせによって平均の差が変わる、とは言えない結果になっています。
 このように、2 因子 ANOVA では、

  • 因子 1 の効果
  • 因子 2 の効果
  • 因子どうしの組み合わせの効果

それぞれ独立に判断します。この例では、いずれについても有意な差は認められませんでした。  最短、最低限の説明では、

2 因子 ANOVA の結果、group、condition、および両者の交互作用のいずれについても有意な効果は認められなかった。

ということです。

5. 事後検定:Tukey の HSD(TukeyHSD)

ANOVA の結果が有意だった場合、次に気になるのは

「では、どのグループ同士が違うのか?」

という点です。ANOVA 自体は、「少なくともどこかは違うか」を判定するだけで、どの組み合わせが違うかまでは教えてくれません。そこで用いるのが 事後検定です。
 ここでは、R で最もよく使われる Tukey の HSD(Honestly Significant Difference)だけを扱います。理由は単純で、良く使われているからです。

■ とりあえず使うコマンド

 ANOVA の結果を fit に保存していれば、事後検定は次の 1 行で実行できます。

TukeyHSD(fit)

これ以上の準備は必要ありません。ANOVA が有意だったら、とりあえずこの行を打ちます。

■ 1因子 ANOVA が有意だった場合と事後検定の解釈

 1 因子 ANOVA の結果が次のようになったとします。

> summary(fit)
            Df Sum Sq Mean Sq F value Pr(>F)  
group        2  9.024   4.512   4.598 0.0191 *
Residuals   27 26.494   0.981
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

この結果では、group 行の Pr(>F) = 0.0191 となっており、有意水準 5% で グループ間の平均に有意な差があると判断できます。つまり、「A、B、C のどこかには平均の異なるグループが存在する」ことまでは、この ANOVA から分かりました。ただし、この段階では どのグループ同士が違うのか は分かりません。そこで次に行うのが、事後検定です。
 ANOVA が有意だったため、Tukey の HSD を実行します。結果は次のようになります。

> TukeyHSD(fit)
$group
         diff        lwr      upr     p adj
B-A 0.9347179 -0.1636787 2.033115 0.1066370
C-A 1.3029917  0.2045951 2.401388 0.0176371
C-B 0.3682738 -0.7301229 1.466670 0.6871725

ここでも、すべての列を理解する必要はありません。見るべきなのは、各行の p adj(調整済み p 値)だけです。
 この出力は、グループどうしを 2 つずつ比較した結果を示しています。

  • B − A
    p adj = 0.1066 → 有意水準 5% では有意差なし

  • C − A
    p adj = 0.0176 → 有意差あり

  • C − B
    p adj = 0.6872 → 有意差なし

したがって、このデータでは、

  • C 群と A 群の間にのみ有意な平均差がある
  • B 群は A 群とも C 群とも、有意には異ならない

という結果になります。最短、最低限のアプローチでも、

1 因子 ANOVA の結果、グループ間に有意差が認められた(p = 0.019)。Tukey の HSD による事後検定の結果、C 群と A 群の間にのみ有意な差が認められた(p < 0.05)。

ということができます。
 重要なポイントは、

  • ANOVA が有意でなければ、原則として事後検定は行わない
  • 事後検定の結果だけを単独で報告しない

ということです。「ANOVAで差があるかも」となった場合に限り、事後検定を試すという手順を守ることが、最低限のルールです。

■ 2因子 ANOVA が有意だった場合と事後検定の解釈

 次に、2 因子 ANOVA の結果が次のようになったとします。

> summary(fit)
                Df Sum Sq Mean Sq F value  Pr(>F)   
group            1   0.07   0.066   0.065 0.79968   
condition        1  10.85  10.854  10.757 0.00231 **
group:condition  1   2.07   2.071   2.053 0.16055   
Residuals       36  36.33   1.009                   

2 因子 ANOVA では、見るべき行は 3 つあります。

  • group:因子1の主効果
  • condition:因子2の主効果
  • group:condition:交互作用

この例では、結果は次のように整理できます。
 まず、group の p 値は 0.7997 であり、有意水準 5% では group による平均差は認められません。 次に、condition の p 値は 0.00231 であり、condition による平均差は有意です。 一方、group:condition の p 値は 0.1606 であり、交互作用は有意ではありません
 つまり、この ANOVA から分かることは、

  • group の違いによる平均差はない
  • condition の違いによる平均差はある
  • group によって condition の効果が変わるとは言えない

ということです。この時点で、「どの条件同士が違うのか」を具体的に確認するために、事後検定を行います。

> TukeyHSD(fit)

出力は、因子ごとに分かれて表示されます。

> TukeyHSD(fit)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = y ~ group * condition, data = dat)

$group
          diff        lwr       upr     p adj
B-A 0.08120663 -0.5630273 0.7254405 0.7996806

$condition
        diff       lwr     upr     p adj
Y-X 1.041836 0.3976018 1.68607 0.0023104

$`group:condition`
              diff         lwr       upr     p adj
B:X-A:X -0.3739217 -1.58380496 0.8359616 0.8387368
A:Y-A:X  0.5867074 -0.62317590 1.7965906 0.5651872
B:Y-A:X  1.1230423 -0.08684095 2.3329256 0.0769453
A:Y-B:X  0.9606291 -0.24925420 2.1705123 0.1604911
B:Y-B:X  1.4969640  0.28708075 2.7068473 0.0103190
B:Y-A:Y  0.5363349 -0.67354832 1.7462182 0.6345845

 ここで、group の比較

$group
B-A 0.08120663  p adj = 0.79968

については、ANOVA の結果と同様に、A と B の間に有意差はありません

 つぎに、condition の比較

$condition
Y-X 1.041836  p adj = 0.00231

については、Y と X の間に有意差があります。これは、ANOVA で condition の主効果が有意だったことと一致しています。

 さらに、group × condition(交互作用)の比較

$group:condition
B:Y-B:X  p adj = 0.0103
(その他は p adj ≥ 0.05)

については、複数の組み合わせ比較が並びますが、結果は、

  • ほとんどの組み合わせでは有意差なし
  • B 群における Y と X の差(B:Y − B:X)のみが有意

となっています。ただし、ANOVA 本体では交互作用が有意ではなかったため、これらの結果は補足的な情報として扱うのが妥当です。交互作用が有意でない場合、組み合わせごとの差を過度に強調するべきではありません。
 この結果を、最短・最低限でまとめると次のようになります。

2 因子 ANOVA の結果、condition の主効果のみが有意であり(p = 0.0023)、group の主効果および交互作用は有意ではなかった。Tukey の HSD による事後検定の結果、condition の水準 X と Y の間に有意な差が認められた。

6. まとめ

 「すぐに結果を知りたい」という気持ちはわかります。ですので、まずは、お手軽にRで試してみれば良いでしょう。しかし、「なぜ」という気持ちを大切にして、理解を深めることも忘れないでください。過去の記事も参考にしてください。

chaosmemo.com

chaosmemo.com

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。