1．分散分析
実験、観測、調査などでは同じの条件であっても、計測の誤差やノイズなどが混入され、得られているデータにはずれが多かれ少なかれ生じる。また同様な実験、観測、調査を、条件を変えて行ったとき、計測の誤差やノイズなど以外に条件の影響で違いが生み出される可能性がある。実験、観測、調査の結果に影響をおよぼすと考えられる要因はいろいろあるが、その実験、観測、調査で取り上げている要因を因子(factor)と呼び、因子を細分類したグループ(群)を水準(levels)と呼ぶ。
分散分析(analysis of
variance; ANOVA)は、得られた各水準の平均値が因子の影響により変化されていると言えるかどうかに関するデータ分析の方法である。本稿では、一元(one-way )分散分析と二元(two-way)分散分析の簡単な例を用いてRによる分散分析について説明する。
2．一元分散分析
実験、観測、調査で取り上げている要因・因子(factor)が1つである分散分析を一元分散分析と呼ぶ。そのデータの形式(一元配置)を表1に示す。表1の中の
(
)は因子Aの水準である。
表1 一元配置データ表
|
標本 |
因子 |
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
平均 |
|
|
|
|
|
|
一元分散分析は対応なし、対応ありに大きく分けることができる。対応とは、
は同一の対象
(
)についてそれぞれの水準の下で得られたデータであるか否かである。具体的な例を用いて説明する。
2.1 対応なしの場合
例えば、ある機器の性能について、異なる温度の下でそれぞれ10回のテストをランダムに行い、表2のようなデータが得られたとする。この場合、各水準(列)内のデータの行の順序には意味がない。例えば、表2のセル(1，1)のデータ63を第1列内であればどの行と入れ換えてもよい。かつそのとき他の列の第1行のデータをセル(1，1)に連動させなくてもデータが持っている情報には変わりがない。問題を簡単にするためここでは機器の疲労の影響などについては無視する。
表1 繰り返しの機器のテスト結果(架空)
(対応なし)
|
回数 |
要因・因子A |
|||
|
－ |
|
|
|
|
|
1 2 3 4 5 6 7 8 9 10 |
63 58 64 58 77 66 52 64 49 66 |
64 64 68 61 56 71 64 65 85 75 |
59 87 79 71 65 65 65 71 74 58 |
83 79 65 67 80 72 80 75 72 84 |
|
平均 |
61.7 |
67.3 |
69.4 |
75.7 |
ここでは、機器の平均性能が条件の影響を受けて変化していると言えるかどうかについて興味を持っている。要因・因子の影響による観測値の変化を因子効果と呼ぶ。一元分散分析はこのような多群の平均値の検定に関するデータ分析方法である。
Rでは、分散分析を行う関数aovが用意されている。関数aovに用いるデータの形式はデータフレーム型で、データ
を１列に並べ、それに関連する因子及び水準の情報をそれぞれ異なる列に記録する必要がある。例えば、水準
はa1、水準
はa2、水準
はa3、水準
はa4のように表記する。ここのa1、a2、a3、a4は水準の群(groups)を識別するラベルである。よって、文字A、B、C、D、あ、い、う、え、1、2、3、4などを用いてもよい。
次にデータフレームの作成に関するコマンドを示す。データフレームbunsan1のA列が因子の水準に関する情報で、y列がテストで得られた数値である。関数data.frame中のfactorはAに入力するのは文字であること、rep(“a1”,10)は文字列a1を10回生成するコマンドである。
>a1<-c(63,58,64,58,77,66,52,64,49,66)
>a2<-c(64,64,68,61,56,71,64,65,85,75)
>a3<-c(59,87,79,71,65,65,65,71,74,58)
>a4<-c(83,79,65,67,80,72,80,75,72,84)
>bunsan1<-data.frame(A=factor(c(rep("a1",10),rep("a2",10),rep("a3",10),rep("a4",10))),y=c(a1,a2,a3,a4))
> bunsan1
A y
1 a1 63
2 a1 58
＜中略＞
39 a4 72
40 a4 84
まず視覚的にデータの大まかな状況を把握するためにデータセットbunsan1の箱ひげ図を示す。次のコマンドで図1に示す箱ひげ図が作成される。
>boxplot(y~A,data=bunsan1,col="lightblue")
図１から水準a1、a2、a3、a4の中心値が右上がりになっていることがわかる。この差の有意性について分散分析を行ってみよう。
図1 表1のデータの箱ひげ図
関数aovにデータセットbunsan1を適したコマンドを次に示す。関数aovを用いた分散分析の書式はaov(y~A,data=bunsan1)となるが、分散分析の結果の要約を、次のように関数summaryを用いて返すことができる。
>
summary(aov(y~A,data=bunsan1))
Df Sum Sq Mean Sq F
value Pr(>F)
A
3 1003.27 334.42 5.302 0.003941 **
Residuals 36 2270.70 63.08
Signif.
codes: 0 `***' 0.001 `**' 0.01 `*'
0.05 `.' 0.1 ` ' 1
返された結果を分散表と呼ぶ。各データを記号で対応付けた表を次に示す。
表3 分散分析表
|
|
自由度 Df |
変動 Sum
Sq |
不偏分散 Mean
Sq |
分散比 F
value |
P値 Pr(>F) |
|
因子 |
|
|
|
|
pf(F) |
|
誤差 |
|
|
|
表3の中の記号列について表1の一元配置データ表に基づいてその定義を説明する。データ
は次のように分解することができる。
ただし、記号
、
は次のように定義されている。
総平均：
群(水準)内平均：
分解された
は群内の平均と総平均とのずれであり、
は観測値と群内の平均とのずれである。すべての観測値についてこのようなずれを分析するため次のような統計量を導入する。
次の式により定義される水準の違いにより生じるばらつきを群間の変動(variation between groups)と呼び、
(sums of squares between)で示す。
水準・群の数から1を引いた値
を群間の変動の自由度と言う。群間の変動を自由度で割った値は群間の不偏分散で、
(mean squares between)で示す。
水準内のデータのばらつきを次の式で定義したものを群内の変動(variation within
groups)と呼び、
(sums of squares within)で示す。群内の変動を変動誤差・残差(Residuals)とも呼ぶ。
群内の変動
の自由度はデータの総数
から群数
を引いた値
である。群内の不偏分散を
(mean squares within)で示す。
データの群間に差があるかどうかを比較するには、群間の不偏分散と群内の不偏分散について等分散検定を行うことが考えられる。等分散検定が採択されると群の平均値には有意な差がないと言える。分散の比は
分布にしたがう。群間と群内の不偏分散の比を次に示す。
関数aovの結果には検定統計量
に対応する
値[Pr(>F)]が返される。このPr(>F)を限界水準とも呼ぶ。上記の例ではPr(>F)値が約0.004である。この値は、通常仮説検定に用いる有意水準0.01より小さい。よって、有意水準0.01で「各水準の平均値が等しい」という仮説が棄却され、群間には差があると統計的に判断される。この結果は機器の性能は安定せず、温度の変化に伴い性能が変化していることを意味する。
各水準のデータの数が一致しない場合でも同様な分析の方法が適応できる。
2.2対応ありの場合
前項の例のデータを10台の機器を異なる条件の下で行ったテストの結果として書き直し、表4に再掲する。表4のデータは見かけ上では表1のデータと似ているが、本質的には異なる。表1では同水準においてデータの行の順序には意味がないが、表4では1行が1台の機器に対するテストの結果であるので行のデータがセットになっている。このようなデータを対応ありのデータと言う。
表4 10台の機器のテスト結果
(繰り返しなし、対応あり)
|
機器 の 番号 |
要因・因子A |
|||
|
－ |
|
|
|
|
|
No.1 No.2 No.3 No.4 No.5 No.6 No.7 No.8 No.9 No.10 |
63 58 64 58 77 66 52 64 49 66 |
64 64 68 61 56 71 64 65 85 75 |
59 87 79 71 65 65 65 71 74 58 |
83 79 65 67 80 72 80 75 72 84 |
|
平均 |
61.7 |
67.3 |
69.4 |
75.7 |
3．二元分散分析
表5のように何らかの2つの要因を考慮したデータを二元(Two-way)配置と言う。表5には2つ形式で示しているが本質的には同じである。二元分散分析では、表5のような二元配置データにおける因子Aの効果、因子Bの効果、因子AとBとの交互作用効果などについて分析を行う。ここで言う交互作用効果(interaction)とは、異なる因子の効果が絡んで生じる効果である。
表5 二元配置のデータ(対応なし)
表 (a)
|
|
因 子 |
|||
|
|
|
|
||
|
因 子 |
|
3 3 5 4 |
6 4 5 6 |
6 7 8 7 |
|
|
3 5 2 4 |
3 4 5 3 |
3 4 3 2 |
|
表 (b)
|
因子 |
|
|
||||
|
因子 |
|
|
|
|
|
|
|
|
3 3 5 4 |
6 4 5 6 |
6 7 8 7 |
3 5 2 4 |
3 4 5 3 |
3 4 3 2 |
まず表5の二次元配置データを次のようにR用のデータセットbunsan3を作成する。
>a1<-c(3,3,5,4,6,4,5,6,6,7,8,7)
>a2<-c(3,5,2,4,3,4,5,3,3,4,3,2)
>bunsan3<-data.frame(A=factor(c(rep("a1",12),rep("a2",12))),B=factor(rep(c(rep("b1",4),
rep("b2",4), rep("b3",4)),2)),y= c(a1,a2))
> bunsan3
A B y
1
a1 b1 3
2
a1 b1 3
23
a1 b3 3
24
a1 b3 2
交互作用効果を無視し、両因子の効果のみについて分析を行う場合は、関数aovを次のように用いる。
>
summary(aov(y~A+B,data=bunsan3))
Df Sum Sq Mean Sq F value
Pr(>F)
A
1 22.042 22.042 13.8482
0.001349 **
B
2 7.750 3.875 2.4346 0.113161
Residuals 20 31.833
1.592
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
返された結果のAの行が因子Aの分散分析の結果である。そのPr(>F)値は約0.0014で、有意水準0.01より小さい。よって、A因子の2つの水準
は有意水準0.01でその平均値には差があると判断される。Bの行が因子Bの分散分析の結果であり、そのPr(>F)値は0.113である。よって、B因子の3つの水準
は有意水準0.1でも有意の差があると言えない。
表5の因子Aの水準
のみを見た限りでは水準
には異なりがありそうにも関わらずこのような大きなPr(>F)値が得られている。そこで、交互作用効果を考慮した分散分析の統計量を求めてみよう。交互作用効果を考慮した場合は、次のように因子の間に記号*を用いる。
>
summary(aov(y~A*B,data=bunsan3))
Df Sum Sq Mean Sq F
value Pr(>F)
A
1 22.0417 22.0417 23.0000 0.0001447 ***
B
2 7.7500 3.8750 4.0435 0.0354521 *
A:B 2
14.5833 7.2917 7.6087 0.0040287 **
Residuals 18 17.2500 0.9583
Signif.
codes: 0 `***' 0.001 `**' 0.01 `*'
0.05 `.' 0.1 ` ' 1
返された結果では、因子BのPr(>F)値は0.035で有意水準0.05より小さくなっており有意水準0.05で
に有意な差があると判断される。A:Bの行が交互作用効果に関する分散分析の結果である。そのPr(>F)値が0.004で、0.01より小さいので有意水準0.01で交互作用効果があると判断される。
交互作用効果に関しては折れ線グラフでも視覚的に考察することができる。Rには交互作用グラフを作成するための関数interaction.plotが用意されている。次に関数interaction.plotを用いた交互作用図の作成コマンドと結果を示す。
>attach(bunsan3)
>interaction.plot(A,B,y)
図2 交互作用図
交互作用効果が存在しない場合は、折れ線が平行になり、交互作用効果が存在する場合は平行にならず、交互作用効果が大きい場合は交差現象が伴う。図2では交互作用効果があることがわかる。交互作用に関しては、本誌の連載「統計数理はいま･･･」で連載中である(第56回(2004年2月)～)分散分析の専門家広津千尋教授の文章が詳しい。
二元分散分析は、データの対応ありとなしで表6に示すように4つのタイプに分けられる。紙面上の都合により、本稿では二元配置についてタイプⅠの例のみを示している。分散分析は各タイプのデータ構造に基づいて分析を行わなければならない。分散分析のソフトを用いる際には、データ構造のタイプごとに分析方法を選択し、引数などの設定を行わなければならない。初心者にとっては多少煩雑に感じるであろう。
表6 二元配置のタイプ
因子A |
因子B |
|
タイプⅠ |
対応なし |
対応なし |
タイプⅡ |
対応あり |
対応あり |
タイプⅢ |
対応なし |
対応あり |
タイプⅣ |
対応あり |
対応なし |
中部大学人文学部心理学科松井孝雄助教授は「データ解析テクニカルブック、森・吉田編著、北大道書房」の中の十数タイプのデータを用いて「言語Rによる分散分析」というタイトルでR用のデータセットの作成及び分散分析のコマンドを公開している。
http://mat.isc.chubu.ac.jp/
3．分散分析の関連関数