[連載]フリーソフトによるデータ解析・マイニング 12回

Rと分散分析


 

 

1.分散分析

 実験、観測、調査などでは同じの条件であっても、計測の誤差やノイズなどが混入され、得られているデータにはずれが多かれ少なかれ生じる。また同様な実験、観測、調査を、条件を変えて行ったとき、計測の誤差やノイズなど以外に条件の影響で違いが生み出される可能性がある。実験、観測、調査の結果に影響をおよぼすと考えられる要因はいろいろあるが、その実験、観測、調査で取り上げている要因を因子(factor)と呼び、因子を細分類したグループ()を水準(levels)と呼ぶ。

分散分析(analysis of variance; ANOVA)は、得られた各水準の平均値が因子の影響により変化されていると言えるかどうかに関するデータ分析の方法である。本稿では、一元(one-way )分散分析と二元(two-way)分散分析の簡単な例を用いてRによる分散分析について説明する。

 

2.一元分散分析

実験、観測、調査で取り上げている要因・因子(factor)1つである分散分析を一元分散分析と呼ぶ。そのデータの形式(一元配置)を表1に示す。表1の中()は因子Aの水準である。

 

1 一元配置データ表

標本

 

因子

*

 

 

 

平均

 

 

 

一元分散分析は対応なし、対応ありに大きく分けることができる。対応とは、は同一の対象()についてそれぞれの水準の下で得られたデータであるか否かである。具体的な例を用いて説明する。

2.1 対応なしの場合

例えば、ある機器の性能について、異なる温度の下でそれぞれ10回のテストをランダムに行い、表2のようなデータが得られたとする。この場合、各水準()内のデータの行の順序には意味がない。例えば、表2のセル(11)のデータ63を第1列内であればどの行と入れ換えてもよい。かつそのとき他の列の第1行のデータをセル(11)に連動させなくてもデータが持っている情報には変わりがない。問題を簡単にするためここでは機器の疲労の影響などについては無視する。

 

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に用いるデータの形式はデータフレーム型で、データを1列に並べ、それに関連する因子及び水準の情報をそれぞれ異なる列に記録する必要がある。例えば、水準a1、水準a2、水準a3、水準a4のように表記する。ここのa1a2a3a4は水準の群(groups)を識別するラベルである。よって、文字ABCD、あ、い、う、え、1234などを用いてもよい。

次にデータフレームの作成に関するコマンドを示す。データフレームbunsan1A列が因子の水準に関する情報で、y列がテストで得られた数値である。関数data.frame中のfactorAに入力するのは文字であること、rep(“a1”,10)は文字列a110回生成するコマンドである。

 

>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")

 

図1から水準a1a2a3a4の中心値が右上がりになっていることがわかる。この差の有意性について分散分析を行ってみよう。

 

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

 

Rでこのデータを、関数aovを用いて分析するためには、前項のデータセットbunsan1をそのまま用いることはできない。それはデータセットbunsan1には対応ありに関する情報が欠けているからである。対応ありに関する情報を付加したデータセットを作成するコマンドを次に示す。ここでは前項で作成したデータベクトルa1a2a3a4を用いる。コマンドの中のAは水準に関する情報であり、Noはデータの「対応あり」に関する情報である。

 

>bunsan2<-data.frame(A= factor(c(rep("a1",10), rep("a2",10), rep("a3",10), rep("a4",10))),No= factor(rep(1:10, 4)),y=c(a1,a2,a3,a4))

> bunsan2

     A  No  y

1   a1   1  63

2   a1   2  58

<中略>

39  a4   9  72

40  a4  10  84

 

 関数aovを用いた書式と返された結果を次に示す。

 

> summary(aov(y ~ A+No, bunsan2))

            Df  Sum Sq Mean Sq F value  Pr(>F) 

A            3 1003.27  334.42  4.3695 0.01243 *

No           9  204.22   22.69  0.2965 0.96952 

Residuals   27 2066.48   76.54                 

Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

 

返されたAの行が水準・群に関する分散分析の検定統計量であり、そのPr(>F)値・限界水準は約0.012である。このPr(>F)値は通常の仮説検定に用いる有意水準0.05より小さい。よって因子Aの水準・群の平均値は有意水準0.05で差があると言える。Noの行が機器に関する分散分析の検定統計量で、そのPr(>F)値は約0.97 である。この値は通常用いられている有意水準よりはるかに大きく、10台の機器には差があると判断できない。つまり、因子Aの水準の平均値には差が認められるが、機器の間には有意の差が認められない結果となっている。

このように、見かけ上は同じデータであってもデータを収集する際の前提条件が異なると分散分析を行う方法や結果の解釈も異なる。分散分析は、基本的な考え方は同じであるがデータセットのタイプによってその解析方法が微妙に異なる。

 

3.二元分散分析

5のように何らかの2つの要因を考慮したデータを二元(Two-way)配置と言う。表5には2つ形式で示しているが本質的には同じである。二元分散分析では、表5のような二元配置データにおける因子Aの効果、因子Bの効果、因子ABとの交互作用効果などについて分析を行う。ここで言う交互作用効果(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

 

 返された結果では、因子BPr(>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(20042))分散分析の専門家広津千尋教授の文章が詳しい。

 二元分散分析は、データの対応ありとなしで表6に示すように4つのタイプに分けられる。紙面上の都合により、本稿では二元配置についてタイプTの例のみを示している。分散分析は各タイプのデータ構造に基づいて分析を行わなければならない。分散分析のソフトを用いる際には、データ構造のタイプごとに分析方法を選択し、引数などの設定を行わなければならない。初心者にとっては多少煩雑に感じるであろう。

 

6 二元配置のタイプ

 
因子A
因子B
タイプT
対応なし
対応なし  
タイプU
対応あり
対応あり 
タイプV
対応なし
対応あり
タイプW
対応あり
対応なし

 

中部大学人文学部心理学科松井孝雄助教授は「データ解析テクニカルブック、森・吉田編著、北大道書房」の中の十数タイプのデータを用いて「言語Rによる分散分析」というタイトルでR用のデータセットの作成及び分散分析のコマンドを公開している。

http://mat.isc.chubu.ac.jp/

 

3.分散分析の関連関数

本稿では、分散分析についてRに用意された関数aovを用いた分散分析の例を示している。Rでは線形モデルの関数lmと関数anovaを用いて同様な分散分析を行うことも可能である。これらの関数による分散分析は、誤差が標準正規分布にしたがうと言う仮定のもとで行っている。しかし、本稿では分散の正規性などについてはまったく触れる余裕がなかった。

一元分散分析はこれらの関数以外に、Rに用意されている関数pairwise.t.testoneway.testなどを用いて平均検定を行うことも可能である。また関数kruskal.testによる分散の正規性の仮定を置かないノンパラメトリックな分散分析を行うことも可能である。

二元(多重)配置のノンパラメトリックな分散分析方法としてはフリードマン(Firdman)の関数friedman.testRに用意されている。

また群馬大学社会情報学部青木教授は自作のプログラムを次のホームページで公開している。

http://aoki2.si.gunma-u.ac.jp/R/