4.10 分散分析

ここまでの内容では、２グループ間の平均の差に関する分析を捉えた。しかしながら、三つ以上のグループ間の平均の差に関心があることもある。例えば、異なる地域における売上高の差を比較したい場合が挙げられる。このような目的を持つ場合によく用いられるのが分散分析（Analysis of Variance; ANOVA）である。

ANOVAの構造については、要因と水準という二つの要素から説明が行われる。要因とは、観測値に影響を与えていると考えうるカテゴリ変数のことを指し、水準とは、要因を構成するいくつかの条件やグループを指す。例えばある小売企業における各店舗の一定期間内の売上高が出店エリア特性によって差があるのかという問いに関心があるとする。その際各店舗を、都市エリア、郊外エリア、農村エリアという三つのグループに分類し、それぞれのグループにおける標本平均を求めれば、エリアごとの差を分析できるだろう。この場合、「地域」が売上高に影響を与えうる要因であり、「都市、郊外、農村」という三つのグループが水準だと言える。

分析において取り上げる要因が一つである分散分析を一元配置分散分析と呼ぶ。二元配置分散分析については以下の節でモデル表現と共に紹介している。

4.10.1 分散分析についてのモデル表現*{#ANOAtheory}

本節では、分散分析の構造をモデルによって表現することでより明確な理解を促すことを目指す。ここでは、要因をA、グループ（水準）の数をJとする。第\(j\)グループで\(n_j\)個の観測値が得られる。第\(j\)水準における\(i\)番目の個体の観測値を\(y_{ij}\)と表し、以下が成り立つとする。 \[ y_{ij}=\mu_j+e_{ij},~~~~e_{ij} \sim N(0,\sigma^2), \]

ただし、\(i=1,2,...,n_j\); \(j=1,2,...,J\)であるとする。すべての\(y_{ij}\)は独立であり、 \[ y_{ij}\sim N(\mu_j, \sigma^2) \] が成り立つ。したがって、\(E(y_{ij})=\mu_j\)と、\(V(y_{ij}=\sigma^2)\)も成り立つ。

上記モデルの関心は各グループで平均に差が存在するか否かである。そこで、以下の\(\mu\)を全体平均とする。 \[ \mu=\frac{n_1}{n}\mu_1+\frac{n_2}{n}\mu_2+...+\frac{n_J}{n}\mu_J=\frac{1}{n}\sum^J_{j=1}n_j\mu_j. \]

また、 \(\alpha_j=\mu_j-\mu\)を第 \(j\)グループの効果と呼び、以下のように表現できる。 \[ y_{ij}=\mu+\alpha_j+e_{ij} , \] このとき定義より、 \[ \frac{1}{n}\sum^J_{j=1}n_j\alpha_j=0, \]が成り立つ。もし、\(\alpha_1=\alpha_2=...=\alpha_J=0\)ならば、以下のような水準間の差が存在しないモデルとなる \[ y_{ij}=\mu+e_{ij}, \]

分散分析の第一の目的は、要因が観測値に対して影響を与えているかどうか、すなわち、\(\mu_1=\mu_2=...=\mu_J\)が成立しているかどうかを調べることである。これは全てのグループで主効果\(\alpha_j\)がゼロであるかどうかを調べることに等しい。そのため、以下を検定することと等しい。

\[ \begin{cases} |t|>t_{\alpha/2}(n-2) & \Rightarrow \text{H0を棄却する。}\\ |t|\leq t_{\alpha/2}(n-2)& \Rightarrow \text{H0を採択する。} \end{cases} \]

\[ \begin{cases} H_0:~ &\alpha_1=\alpha_2=...=\alpha_J=0\\ H_1:~ &\text{少なくとも１つの}\alpha_j \text{は0ではない} \end{cases} \]

ここで、総変動: \(S_T=\sum^J_{j=1}\sum^{n_j}_{i=1}(y_{ij}-\bar{y})^2\)を、群間変動: \(S_B=\sum^J_{j=1}n_j(\bar{y_j}-\bar{y})^2\)と、郡内変動: \(S_W=\sum^J_{j=1}\sum^{n_j}_{i=1}(y_{ij}-\bar{y_j})^2\)とに分ける(ただし、\(S_T=S_B+S_W\))。このとき、群間変動が大きいほど各\(\alpha\)の推定値\(\hat{\alpha_j}=\bar{y_j}-\bar{y}\)のばらつき、すなわちグループ間での平均の違いが大きく、対立仮説\(H_1\)が支持されやすい。

しかしながら、この値自体は単位に依存し意味が無いため、以下の群間変動の郡内変動に対する比に基づいて検定を行う。 \[ F=\frac{S_B/(J-1)}{S_W/(n-J)}=\frac{n-J}{J-1}\times \frac{S_B}{S_W} \]

Fは帰無仮説が正しい時に自由度（\(J-1\), \(n-J\)）のF分布に従うことが知られている。そのため、有意水準 \(\alpha\) に基づく\(H_0\)に関する検定は、以下のように表すことができる。

\[ \begin{cases} F>F_\alpha(J-1, n-J)\Rightarrow H_0\text{を棄却}\\ F\leq F_\alpha(J-1, n-J)\Rightarrow H_0\text{を採択} \end{cases} \]

4.10.2 二元配置分散分析*{#twowayANOA}

一元配置分散分析観測値に影響を与える要因が2個あるモデルを考える。この場合の最も特徴的な違いが交互作用効果である。2つの要因をAとBとし、Aのグループの数をJ、Bのグループの数をKとする。（\(A_j\), \(B_k\)）につきr個の観測値\(y_{ijk}\)（\(i=1,2,...,r\)）が得られるとし、以下のモデルを考える。 \[ y_{ijk}=\mu_{jk}+e_{ijk}, ~e_{ijk}\sim N(0,\sigma^2) \] ただし、\(i=1,2,...,r\); \(j=1,2,...,J\); \(k=1,2,...,K\)。

前節の内容と同様に全体平均と要因効果とに分解し、全体平均と呼ぶ。 \[ \mu=\frac{1}{JK}\sum^J_{j=1}\sum^K_{k=1}=\mu_{jk}, \] そして、全体平均を以下のように分解し、前者を行平均、後者を列平均と呼ぶ。 \[\begin{eqnarray} \mu_{j.}=\frac{1}{K}\sum^K_{k=1}=\mu_{jk}, ~~ \mu_{.k}=\frac{1}{J}\sum^J_{j=1}=\mu_{jk}, \end{eqnarray}\]

このとき、要因Aの第\(j\)グループの効果（\(\alpha_j\)）、要因Bの第\(k\)グループの効果（\(\beta_k\)）をそれぞれ、以下のように定義する。 \[ \alpha_j=\mu_{j.}-\mu, ~~\beta_k=\mu_{.k}-\mu \] この時、以下の恒等式が成り立つ。 \[ \mu_jk=\mu+(\mu_{j.}-\mu)+(\mu_{.k}-\mu)+(\mu_{jk}-\mu_{j.}-\mu_{.k}-\mu)\nonumber\\ =\mu+\alpha_j+\beta_k+(\alpha \beta)_{jk}, \] ただし、効果\(\alpha_j\), \(\beta_k\)を要因AとBの主効果、\((\alpha \beta)_{jk}\)をABの交互作用効果と呼ぶ。そして、定義より以下を示すことができる。 \[ \sum^J_{j=1}\alpha_j=0, ~~\sum^K_{k=1}\beta_k=0, ~~\sum^J_{j=1}(\alpha \beta)_{jk}=0, ~~ \sum^K_{k=1}(\alpha \beta)_{jk}=0 \] なお、交互作用効果が無いモデルは、以下のように示すことができる。 \[\begin{eqnarray} \mu_{jk}=\mu+\alpha_j+\beta_k \end{eqnarray}\]

4.10.3 二元配置モデルの効果検定*{#twowaytest}

二元配置モデルにおいても、以下のような要因Aの主効果に関する仮説に関心がある。 \[ \begin{cases} H_0:~ \alpha_1=\alpha_2=...=\alpha_J=0\\ H_1:~ \text{少なくとも１つの}\alpha_jは0\text{ではない} \end{cases} \]

また、要因Bの主効果に関する仮説は以下のように定義できる。 \[ \begin{cases} H_0:~ \beta_1=\beta_2=...=\beta_K=0\\ H_1:~ \text{少なくとも１つの}\beta_kは0\text{ではない} \end{cases} \]

これらに加え、交互作用効果に関する仮説は以下のように示すことができる。

\[ \begin{cases} H_0:~ (\alpha \beta)_{11}=(\alpha \beta)_{12}=...=(\alpha \beta)_{JK}=0\\ H_1:~ \text{少なくとも１つの}(\alpha \beta)_{jk}\text{ではない} \end{cases} \]

これらの仮説の検定ではまず、以下のように「総変動」、「A間変動」、「B間変動」、「交互作用変動」、「グループ内変動」を定義する。ただし以下では、\(S_T=S_A+S_B+S_{AB}+S_W\) とする。 \[ \text{総変動:} S_T=\sum^r_{i=1}\sum^J_{j=1}\sum^K_{k=1}=(y_{ijk}-\bar{y})^2 \] \[ \text{A間変動:} ~S_A=Kr\sum^J_{j=1}(\bar{y_{j.}}-\bar{y})^2 \] \[ \text{B間変動:}~S_B=Jr\sum^K_{k=1}(\bar{y_{.k}}-\bar{y})^2 \] \[ \text{交互作用変動:}~S_{AB}=r\sum^K_{k=1}\sum^J_{j=1}(\bar{y_{jk}}-\bar{y_{j.}}-\bar{y_{.k}}-\bar{y})^2 \] \[ \text{グループ内変動:}~S_W=\sum^r_{i=1}\sum^K_{k=1}\sum^J_{j=1}(y_{ijk}-\bar{y_{jk}})^2 \]

これらを用いて、以下のような \(F_A\)、\(F_B\)、\(F_{AB}\) を検定統計量として検定を行う。 \[ F_A=\frac{S_A/(J-1)}{S_W/JK(r-1)} \] \[ F_B=\frac{S_B/(K-1)}{S_W/JK(r-1)} \] \[ F_{AB}=\frac{S_{AB}/(J-1)(K-1)}{S_W/JK(r-1)} \]

4.10.4 多重比較とそれに伴う統計的問題

ANOVAにおける帰無仮説の棄却は、少なくとも1つの群は全体と異なる平均値を持っているという結論につながる。しかし、どれが異なる値を持つのかがわからないため、それを特定化するために多重比較と呼ばれる分析を行う。この多重比較には統計的問題が伴うと言われている。本節では多重比較に伴う統計的問題について考える。

多重比較に伴う問題は、第一種の誤り（帰無仮説が正しいときに帰無仮説を棄却する誤り）を起こす確率（有意水準）に着目し議論される。ANOVAにおける帰無仮説は以下のとおりである。 \[ H_0:~ \mu_1=\mu_2=...=\mu_J \] この帰無仮説を有意水準5%を設定し、 (8) 式の統計量を用いて検定した場合の\(H_0\)に関する有意水準はもちろん5%になる。しかしながら、各グループごとの対となる検定を5%水準で行った場合、\(H_0:~ \mu_1=\mu_2=...=\mu_J\)に対する有意水準は5%よりも大きくなってしまうと言われている。

この問題では、先述の\(H_0\)に基づき、分散分析が捉えている全体としての第一種の誤りと多重比較で捉えている個別の第一種の誤りという構造を把握する必要がある。言い換えると、多重比較により複数の検定を繰り返し行う場合、その複数の検定を「検定のセット」と考え、セット全体での第一種の誤りについて考えるということになる。

この統計的問題についてより明確に説明するために、1要因3グループの分散分析という具体例を考える。この時、分散分析としての全体の帰無仮説は\(H_0:~ \mu_1=\mu_2=\mu_3\)となる。ここで対となる2群環の平均の差の検定を行うと、\({}_3C_2=3\) 通りのt検定を行うことになる。例えば、グループ1対グループ2の検定においては \(H_0^1:~\mu_1=\mu_2\) という帰無仮説を採用して検定を行うことになる。これら各対の検定における帰無仮説が棄却されれば、分散分析全体として捉えていた帰無仮説（\(H_0:~ \mu_1=\mu_2=\mu_3\)）も当然棄却されることになる。そのため、各対に関する帰無仮説が少なくともひとつ棄却されることで、全体としての帰無仮説が棄却されるという構造を捉える必要が出てくる。この構造が、第一種の誤りについての問題につながる。

各対に関する帰無仮説に対して第一種の誤りが起こる確率は、「ひとつ目の対の検定において誤りが起きる、またはふたつ目の対の検定において誤りが起きる、またはみっつ目の対の検定において誤りが起きる」という事象の確率を捉えることになる。仮に各対の検定において有意水準5%を設定した場合、3回に1回以上第一種の誤りを起こす確率（全体としての有意水準）は、それぞれの検定において第一種の誤りを犯すという事象が互いに独立だと仮定すると、以下の通り0.05よりも多くなる。 \[ 1-P(\text{3回中1回も第一種の誤りを犯さない})=1-(0.95)^3\approx 0.143 \]

そのため、分散分析に対応する複数の検定を繰り返し行う場合には、複数の検定を検定のセットと捉え、セット全体での第1種の誤りを犯す確率が事前に設定した有意水準（例：0.05等）以下に抑えられるような方法が求められる。この問題に対応した方法が多重比較法であり、代表的な例としてTukeyの多重比較法がある (Tukey’sのPairwise法については例えば南風原（2010)を参照してほしい)。

4.10.5 Rによる分散分析の実行

Rで分散分析を実行すること自体は難しくない。ここでは reshape2 データを用いて、以下の通りチップ額が曜日によって異なるか否かを分析する。ANOVAの実行においてはaov()関数によって分析するモデルとデータを指定し、anova() によって分析結果を出力する（summary() を用いることも可能である）。なお、上記のデータサマリーより、day という要因には四水準 (levels) 含まれていることがうかがえる。

library(reshape2)

s <- aov(tip ~ day, data = tips)
anova(s)

## Analysis of Variance Table
## 
## Response: tip
##            Df Sum Sq Mean Sq F value Pr(>F)
## day         3   9.53  3.1753  1.6724 0.1736
## Residuals 240 455.69  1.8987

分析結果における Sum Sq は、水準間変動和（または、平方和）と呼ばれ、特定要因の水準間によって説明される観測値の変動を表している。 Mean Sq は平均平方と呼ばれ、Sum SqをDF(自由度)で割ったものである。F value は F値という検定統計量の実現値であり、Pr(>F) は本検定の p値を表している。また、Residuals の行で示されているのは、残差平方和と呼ばれ、郡内変動、つまり同グループ（水準）内での値のばらつきの程度を表している。

本結果の場合、曜日間でのチップ額についての差についての推定値（F value）は1.6724 であった。そして、これが統計的に有意なのかについて検討するためにp値を参照する。このときp値は 0.1736 であり、帰無仮説が正しいと仮定したときに，手元のデータから計算した検定統計量以上に極端な値をとる確率は17％以上あるといえる。このことから、有意水準 5% で仮説検定を行っても帰無仮説を棄却できないと解釈でき、曜日による統計的に有意な差があるとは言えない。ANOVA実行している検定では、「すべての水準間で平均値は同じ」という帰無仮説を検定しており、仮に帰無仮説が棄却された場合、「少なくとも一つの水準では値が異なる」という対立仮説を支持する。そのため、ANOVAにおける帰無仮説の棄却は、少なくとも1つの群は全体と異なる平均値を持っているという結論につながる。今回の結果では、このような帰無仮説も棄却するに至らなかったということである。

もし仮にANOVAにおける帰無仮説が棄却されたとしても、それだけでは具体的にどの水準間に差があるのかはわからない。そこで、ANOVAを用いた研究では事後分析として、多重比較と呼ばれる分析を行うことが多い。しかし、この多重比較には統計的な問題が伴うと言われている（詳細は4.10.4節参照）。その問題に対応した手法として広く用いられているのが、Tukeyのpair-wise 比較である。Rではこの分析を、TukeyHSD() という関数で実行できる。具体的には、aov() でストアしたANOVAの分析結果を用いて以下のように実施する。また、plot()を用いて、pair-wise比較の結果を95%信頼区間とともに図示化することもできる。

tukey_result <- TukeyHSD(s)
tukey_result

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = tip ~ day, data = tips)
## 
## $day
##                 diff        lwr       upr     p adj
## Sat-Fri   0.25836661 -0.6443694 1.1611026 0.8806455
## Sun-Fri   0.52039474 -0.3939763 1.4347658 0.4558054
## Thur-Fri  0.03671477 -0.8980753 0.9715049 0.9996235
## Sun-Sat   0.26202813 -0.2976929 0.8217492 0.6203822
## Thur-Sat -0.22165184 -0.8141430 0.3708394 0.7678581
## Thur-Sun -0.48367997 -1.0937520 0.1263921 0.1724212

plot(tukey_result)

分析結果における diff 列は平均値の差を表している。 lwer と upr は信頼区間の下限と上限を表しており、一番右の列はp値を示している。分析の結果、P-valueが10%水準よりも低い結果がないため、どのペアに関する検定でも有意な差は確認できなかった。したがって、分析結果は必ずしもチップ額が曜日によって変化するとは言えないことを示した。この結果は、アメリカにおけるチップ額が会計額に対する割合や提供されたサービス品質によって決まるという慣習から考えると妥当な結果である。しかしながら、前節で注意した通り、このような統計的に非有意な結果をもって「曜日はチップ額に影響を与えない」と結論づけるのは不適切である。

なお、今回の分析においてはANOVAもTukeyのpair-wise比較も有意な結論を得ることができなかったという点で、両者の結果に一貫性があった。しかしながら、ANOVAでは有意だが、Tukeyの分析ではどの組み合わせも有意ではないという一見整合でない結果を得ることもある。その場合には、慣習としてTukeyの多重比較結果を優先して解釈を提示することが多い。しかしながら、研究者が自身の実施した検定の「意味」を理解し、解釈や議論を提示することも重要である。例えば、ANOVAとTukeyでは用いている帰無仮説が異なるため、異なる比較対象を用いた検定を実施している。そのため、自身が実行した分析がどのような帰無仮説を採用しており、何と何の比較を行っているのかを正確に把握し、実施した検定の意味に適した解釈や議論を展開する事が重要になる。

本章では、基礎的な統計学の復習として、主に区間推定と統計的仮説検定について説明した。区間推定では、主に信頼区間の計算に着目し、信頼区間の意味についてきちんと理解、解釈することの重要性を強調した。また、統計的仮説検定では、母平均の検定を起点とし統計的検定の基礎的な構造と考え方について説明した。検定においては、母集団の統計的特徴に関する予測である帰無仮説と対立仮説を設計することが重要である。また、グループ間の差異に着目した検定を行う場合には、関心のある未知パラメータについての差や比に着目し検定統計量を作成する事が多いが、基本的な統計的検定の考え方と手順は変わらないという点についても説明した。

4.10.6 回帰分析と分散分析の関係

分散分析は回帰分析の特別な場合（回帰分析における独立変数をダミー変数とした場合）に対応する。ダミー変数とは以下の様な2値変数である。 \[ x_i=\left\{ \begin{array}{l} 1~~~\text{if 被験者iは1群に属する}\\ 0~~~\text{otherwise} \end{array} \right. \] そのため、2標本問題は、\(y=\beta_0 +\beta_1x+e,\) において、以下を検定することと等しい。 \[ \begin{array}{l} H_0:~ \beta_1=0\\ H_1:~ \beta_1\not=0 \end{array} \]

また、ダミー変数の効果は以下のように定数項（切片）の違いとして表すことができる。 \[ E(y|x=1)=\beta_0+\beta_1 \] \[ E(y|x=0)=\beta_0 \]

二元配置に関しては、以下のように、交差項 (interaction term) を含む回帰分析として分析を行えば良い。 \[ y=\beta_0 +\beta_1x+\beta_2z+\beta_3(x\times z)+e, \]

この場合、xがyに与える影響は以下のように表される。 \[ \frac{\partial y}{\partial x}=\beta_1+\beta_3z \] 一方で、zがyに与える影響は以下のように表される。 \[ \frac{\partial y}{\partial z}=\beta_2+\beta_3x \] そのため、xとzがyへ与える影響は互いがそれぞれに依存しており、\(\beta_3\)によってそれぞれの影響に対する調整効果が表されている。