人工知能を作りたい!

統計的仮説検定

一元配置分散分析、二元配置分散分析の使い方




一元配置分散分析(対応なし)、Tukey(テューキー)の方法による多重比較

はじめに

t検定も分散分析も、難しい名称が付いていますが、母平均の差の検定を行います。では、t検定と分散分析の違いは何でしょうか。

分散分析は「3つ以上の群それぞれの母平均の差」を検定します。しかし、t検定は「2つの群それぞれの母平均の差」しか検定できません。つまり、まとめると以下のようになります。

2つの群それぞれの母平均の差を検定する場合 → t検定
3つ以上の群それぞれの母平均の差を検定する場合 → 分散分析

分散分析において、母平均に差をもたらす原因を因子といいます。また、因子の中に含まれる個々の条件を水準といいます。そして、因子を1つ配置した場合の分散分析を一元配置分散分析といいます。

全て異なる被験者から得たデータ(被験者間のデータ)は対応がないデータといいます。同じ被験者から異なる水準で得たデータ(被験者内のデータ)は対応があるデータといいます。




例題1:一元配置分散分析(対応なし)

健康な人の血糖値は、空腹時は約70 mg/dL、食後2時間後は約130 mg/dLです。食後のインスリン分泌に遅れのある2型糖尿病患者の血糖値は、空腹時は約70 mg/dLで正常値ですが、食後2時間後は約200 mg/dL以上になります。

新しく開発した4種類のインスリン分泌促進薬の効果を検証するため、インスリン分泌に遅れのある2型糖尿病患者群から40人の被験者を無作為に抽出しました。そして、10人ずつ、4つのグループに分け、グループごとに異なる種類の新薬を割り当てました。

新薬は食前10分前服用してもらい、食後2時間後の血糖値を測定しました。そして、「新薬を服用しない時の食後血糖値」と「新薬を服用した時の食後血糖値」の差を血糖降下度 (mg/dL) として、以下のようにデータ化しました。



<薬A>

 被験者No.  血糖値の差
 1  60
 2  79
 3  127
 4  98
 5  102
 6  79
 7  97
 8  86
 9  109
 10  47


<薬B>

 被験者No.  血糖値の差
 11  58
 12  83
 13  95
 14  102
 15  56
 16  89
 17  69
 18  77
 19  82
 20  64


<薬C>

 被験者No.  血糖値の差
 21  75
 22  89
 23  64
 24  47
 25  49
 26  49
 27  73
 28  102
 29  59
 30  74


<薬D>

 被験者No.  血糖値の差
 31  88
 32  98
 33  132
 34  129
 36  68
 37  95
 37  117
 38  77
 39  89
 40  86


4種類の新薬すべてにおいて、血糖値上昇を抑える効果は等しいと判断していいでしょうか。




帰無仮説と対立仮説の設定

帰無仮説 H0:μA = μB = μC = μD(4群それぞれの母平均は等しい)
対立仮説 H1:少なくとも1つの群の母平均が他の群の母平均と異なる
有意水準 α=0.05

F分布は正の値しか持たないため、常に片側検定となります。




Rによる統計解析

Rには分散分析を実行する関数が用意されています。以下にその2つの関数を示します。

関数 oneway.test()
関数 aov()

まず、下準備を行います。


> グループ1 <-c(60, 79, 127, 98, 102, 79, 97, 86, 109, 47)
> グループ2 <-c(58, 83, 95, 102, 56, 89, 69, 77, 82, 64)
> グループ3 <-c(75, 89, 64, 47, 47, 49, 73, 102 ,59, 74)
> グループ4 <-c(88, 98, 132, 129, 68, 95, 117, 77, 89, 86)


4群のデータをそれぞれ変数「グループ1」、「グループ2」、「グループ3」、「グループ4」に格納します。


> 全グループ<-c(グループ1, グループ2, グループ3, グループ4)


4群のデータすべてを変数「全グループ」に格納します。


> 薬の種類 <-c(rep("A", 10), rep("B", 10), rep("C", 10), rep("D", 10))


変数「薬の種類」を用意します。この変数「薬の種類」に、どの薬を服用したのかを示すAからDという文字型の値を10ずつ繰り返して格納します。

関数 rep(値, 回数) は、指定した値を指定した回数分、繰り返したデータを作成します。

確認のため、出力してみます。


> 薬の種類
 [1] "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "B" "B" "B" "B" "B" "B" "B" "B"
[19] "B" "B" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "D" "D" "D" "D" "D" "D"
[37] "D" "D" "D" "D"


「薬の種類」を出力すると"A"のように " が付いています。これは「薬の種類」に含まれているデータが文字型であることを示しています。

しかし、分散分析を実行する関数において、グループ分けに使う関数は文字型から因子型に直す必要があります。因子型とは「それぞれの文字型の値に対してRが自動的に対応する数値を割り当てたデータ型」です。


> 薬の種類2 <-factor(薬の種類)


関数 factor() を使って文字型データを因子型データに変換します。


試しに「薬の種類2」を出力してみます。


> 薬の種類2
 [1] A A A A A A A A A A B B B B B B B B B B C C C C C C C C C C D D D D D D
[37] D D D D
Levels: A B C D


「薬の種類2」を出力すると「薬の種類」を出力した時と異なり、Aを挟んでいた " が消え、「Levels: A B C D」という行が追加されています。「因子の水準:A、B、C、D」を意味します。

これは、Rの内部で自動的に文字型データA、B、C、Dに「A=1、B=2、C=3、D=4」のように数値を割り当てて「因子の水準」を作り、因子型データに変換したことを示しています。


一元分散分析を実行する準備が整いました。それでは、まず、関数 oneway.test() を利用して統計解析を実行することにします。

関数 oneway.test() の書式は
oneway.test(分散分析の対象となる変数~群分けのための変数)
です。

関数 oneway.test() はデフォルトでは、「全ての群の母分散の等質性は満たされない」と仮定して統計解析が実行されます。ですので、var.equal=TRUEというオプションを指定して、「全ての群の母分散は等しい」と指定します。


> oneway.test(全グループ~薬の種類2, var.equal=TRUE)

        One-way analysis of means

data: 全グループ and 薬の種類2
F = 4.2585, num df = 3, denom df = 36, p-value = 0.0113


Fの実現値、分子の自由度、分母の自由度、p値が示されています。p 値が0.0113となり、有意水準0.05より小さいので帰無仮説は棄却されました。つまり、「対立仮説 H1:少なくとも1つの群の母平均が他の群の母平均と異なる」が採択されました。


続いて、関数aov() を利用してみましょう。

関数aov() の書式は
aov(分散分析の対象となる変数~群分けのための変数)
です。


> aov(全グループ~薬の種類2)
Call:
   aov(formula = 全グループ ~ 薬の種類2)

Terms:
                薬の種類2 Residuals
Sum of Squares 5094.075 14354.700
Deg. of Freedom 3 36

Residual standard error: 19.96852
Estimated effects may be unbalanced


関数 aov() を使って出力すると、結果は見やすいものになっていません。Fの実現値も出力されていません。関数 summary() を利用して結果を見やすくしてみましょう。

関数 summary() と関数 aov() を合わせた書式は
summary(aov(分散分析の対象となる変数~群分けのための変数))
です。


> summary(aov(全グループ~薬の種類2))
            Df Sum Sq Mean Sq F value Pr(>F)
薬の種類2 3 5094 1698.0 4.258 0.0113 *
Residuals 36 14355 398.7
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


上記の表を分散分析表と呼びます。p値が0.0113となり、関数 oneway.test() の結果と同じになりました。




Tukey(テューキー)の方法による多重比較とは

上記の例題において、Rによる統計計解析の結果、「帰無仮説 H0:μA = μB = μC = μD(4群それぞれの母平均は等しい)は棄却され、「対立仮説 H1:少なくとも1つの群の母平均が他の群の母平均と異なる」が採択されました。

しかし、この結果では「どの薬とどの薬の母平均に差があるのか」は具体的にわかりません。差のある「水準と水準の組み合わせ」を判明させるためには、複数の「水準と水準の組み合わせ」において、それぞれ比較する必要があります。

このように、複数の「水準と水準の組み合わせ」において、それぞれ比較することを多重比較と呼びます。

多重比較の方法としてよく使われているのが、Tukey(テューキー)の方法です。Tukey(テューキー)の方法は、各郡の拇分散は等しいと仮定して、すべての「水準と水準の組み合わせ」を一つ残らず比較します。
全ての「水準と水準の組み合わせ」は以下の6通りです。

薬Aと薬B
薬Aと薬C
薬Aと薬D
薬Bと薬C
薬Bと薬D
薬Cと薬D




例題2:Tukey(テューキー)の方法による多重比較

例題1において「対立仮説 H1:少なくとも1つの群の母平均が他の群の母平均と異なる」が採択されたとしたら、どの薬とどの薬の母平均に差があるのでしょうか。




Rによる統計解析

Tukey(テューキー)の方法による多重比較を行うために、Rに用意された関数 TukeyHSD() を利用します。

関数 TukeyHSD() の引数(ひきすう)として、一元配置分散分析(対応なし)を実行した関数 aov(全グループ~薬の種類2) を指定します。


> TukeyHSD(aov(全グループ~薬の種類2))
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = 全グループ ~ 薬の種類2)

$`薬の種類2`
     diff lwr upr p adj
B-A -10.9 -34.951036 13.151036 0.6181037
C-A -20.5 -44.551036 3.551036 0.1180038
D-A 9.5 -14.551036 33.551036 0.7133637
C-B -9.6 -33.651036 14.451036 0.7067349
D-B 20.4 -3.651036 44.451036 0.1206944
D-C 30.0 5.948964 54.051036 0.0096084


D-Cのところを見ると、p値が0.0096084です。有意水準0.05より小さいことから、「C群とD群の間に5%水準で有意な差がある」ということになります。つまり、「薬Cと薬Dの母平均は異なる」ということです。






一元配置分散分析(対応あり)

はじめに

分散分析において、母平均に差をもたらす原因を因子といいます。また、因子の中に含まれる個々の条件を水準といいます。そして、因子を1つ配置した場合の分散分析を一元配置分散分析といいます。

全て異なる被験者から得たデータ(被験者間のデータ)は対応がないデータといいます。同じ被験者から異なる水準で得たデータ(被験者内のデータ)は対応があるデータといいます。

一元配置分散分析において、「各データの値は全データの平均値からどれだけ離れているのか」、つまり、「各データの値のばらつき」は、対応なしの場合と対応ありの場合では以下のように分解され、説明されます。

<一元配置分散分析(対応なし)>
各データの値のばらつき = 群分けによる効果 + 残差

この残差は、「群の違いによる効果では説明できない何か」です。

<一元配置分散分析(対応あり)>
各データの値のばらつき = 群分けによる効果 + 個人差の効果 + 残差

一元配置分散分析(対応あり)では、同じ人が複数の条件で測定されるため、データに個人の特徴が現れます。この残差は、「群分けによる効果でも個人差の効果でも説明できない何か」です。




例題

健康な人の血糖値は、空腹時は約70 mg/dL、食後2時間後は約130 mg/dLです。食後のインスリン分泌に遅れのある2型糖尿病患者の血糖値は、空腹時は約70 mg/dLで正常値ですが、食後2時間後は約200 mg/dL以上になります。

新しく開発した3種類のインスリン分泌促進薬それぞれの効果を検証するため、インスリン分泌に遅れのある2型糖尿病患者群から5人の被験者を無作為に抽出しました。

3種類の新薬それぞれを、被験者一人一人に食前10分前服用してもらい、食後2時間後の血糖値を測定しました。そして、「新薬を服用しない時の食後血糖値」と「新薬を服用した時の食後血糖値」の差を血糖降下度 (mg/dL) として、以下のようにデータ化しました。


   薬A  薬B  薬C
 田中  60  75  88
 福山  79  89  98
 平松  127  64  132
 松田  98  47  129
 高野  102  47  68

3種類の薬の効果に有意な差があるといえるでしょうか。




帰無仮説と対立仮説の設定

帰無仮説 H0: μA = μB = μC(3群それぞれの母平均は等しい)
対立仮説 H1:少なくとも1つの母平均が他の母平均と異なる
有意水準 α =0.05




Rによる統計解析

分散分析(対応あり)を実行するため、関数 aov() を利用します。書式は2通りあります。どちらも同じ結果を表示します。

aov(分散分析の対象となる変数~群分けの変数+個人差の変数)
aov(分散分析の対象となる変数~群分けの変数+Error(個人差の変数))

下準備として、まず、「分散分析の対象となる変数」「群分の変数」「個人差の変数」を作ります。



分散分析の対象となる変数

> 薬Aによる血糖値差 <- c(60, 79, 87, 98, 52)
> 薬Bによる血糖降下度 <- c(55, 89, 64, 97, 47)
> 薬Cによる血糖降下度 <- c(88, 98, 112, 129, 98)
> 全ての血糖降下度 <- c(薬Aによる血糖降下度, 薬Bによる血糖降下度, 薬Cによる血糖降下度)


薬Aによる血糖値差を変数「薬Aによる血糖降下度」へ、 薬Bによる血糖降下度を変数「薬Bによる血糖降下度」へ、 薬Cによる血糖降下度を変数「薬Cによる血糖降下度」へ格納します。

そして次に、変数「薬Aによる血糖降下度」「薬Bによる血糖降下度」「薬Cによる血糖降下度差」を変数「全ての血糖降下度」へ格納します。



群分けの変数

> 薬の種類 <- c(rep("薬A", 5), rep("薬B", 5), rep("薬C", 5))


変数「薬の種類」に、どの薬による血糖降下度なのかを示す「薬A」、「薬B」、「薬C」という文字型データを5回ずつ繰り返して格納しておきます。


> 薬の種類2 <- factor(薬の種類)


関数 factor() を利用して、変数「薬の種類」に格納された文字型データを要因型データに変換し、変数「薬の種類2」に格納しておきます。



個人差の変数

> 人 <- c(rep(c("田中", "福山", "平松", "松田", "高野"), 3))


変数「人」に、誰のデータなのかを示す「田中」、「福山」、「平松」、「松田」、「高野」という文字型データを3回ずつ繰り返して格納しておきます。


> 人2 <- factor(人)

関数 factor() を利用して、変数「人」に格納された文字型データを要因型データに変換し、変数「人2」に格納しておきます。



分散分析表の出力

関数 aov() は分散分析表を出力しないので、関数 summary() を合わせて利用し、分散分析表を出力します。


> summary(aov(全ての血糖降下度~薬の種類2+人2))

            Df Sum Sq Mean Sq F value Pr(>F)
薬の種類2 2 3514 1756.9 20.23 0.000743 ***
人2 4 3639 909.8 10.47 0.002881 **
Residuals 8 695 86.9
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


「薬の種類2」の行のp値は0.000743ですので、群分けによる効果に有意差があることが分かります。つまり、帰無仮説は棄却され、「対立仮説 H1:少なくとも1つの母平均が他の母平均と異なる」が採択されました。ちなみに、「人2」の行のp値は0.002881ですので、個人差による効果にも有意差があることが分かります。


> summary(aov(全ての血糖降下度~薬の種類2+Error(人2)))

Error: 人2
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 4 3639 909.8

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
薬の種類2 2 3514 1756.9 20.23 0.000743 ***
Residuals 8 695 86.9
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘


「薬の種類2」の行のp値は0.000743ですので、群分けによる効果に有意差があることが分かります。つまり、帰無仮説は棄却され、「対立仮説 H1:少なくとも1つの母平均が他の母平均と異なる」が採択されました。summary(aov(全ての血糖降下度~薬の種類2+人2)) の場合と結果は同じです。






くり返し数が等しい二元配置分散分析(対応なし)

はじめに

分散分析において、母平均に差をもたらす原因を因子といいます。また、因子の中に含まれる個々の条件を水準といいます。そして、因子を2つ配置した場合の分散分析を二元配置分散分析といいます。

因子Aの水準と因子Bの水準のどの組み合わせでも測定回数が等しい場合、くり返し数が等しいといいます。因子Aの水準と因子Bの水準のどの組み合わせでも測定回数が1回のみの場合、くり返しがないといいます。

くり返しのない二元配置分酸分析の場合、数式の構造上、交互作用効果を計算することはできません。ですので、もし、2つの因子に交互作用がありそうなら、くり返し数が等しいデータを得た後、くり返し数が等しい二元配置分散分析を行った方が良いです。
全て異なる被験者から得たデータ(被験者間のデータ)は対応がないデータといいます。同じ被験者から異なる水準で得たデータ(被験者内のデータ)は対応があるデータといいます。




例題

健康な人の血糖値は、空腹時は約70 mg/dL、食後2時間後は約130 mg/dLです。食後のインスリン分泌に遅れのある2型糖尿病患者の血糖値は、空腹時は約70 mg/dLで正常値ですが、食後2時間後は約200 mg/dL以上になります。

新しく開発したインスリン分泌促進薬の効果を確かめるため、インスリン分泌に遅れのある2型糖尿病患者群から18人の被験者を無作為に抽出しました。

それぞれの被験者に、新薬を食前10分前服用してもらい、食後2時間後に血糖値を測定しました。そして、「新薬を服用しない時の食後血糖値」と「新薬を服用した時の食後血糖値」の差を血糖降下度 (mg/dL) としてデータ化しました。

因子とその水準を以下の通り設定しました。

因子A:種類の違い(水準:新薬A、新薬B)
因子B:服用量の違い(水準:100μg、500μg、1000μg)

6通りある因子A水準と因子B水準の組み合わせそれぞれに、18人の被験者から3人ずつを割り当てました。つまり、以下のデータは、繰り返し数3回で対応のないデータです。

   100μg  500μg  1000μg
 薬A  30
 25
 28
 49
 59
 68
 80
 76
 82
 薬B  80
 65
 59
 79
 89
 78
 70
 86
 65

種類と服用量の違いによって、血糖降下度の母平均は異なるといえるでしょうか。




主効果と交互作用

くり返し数が等しい二元配置分酸分析では、二つの因子それぞれの主効果と、因子同士がお互いに影響を与え合って生じる交互作用を考える必要があります。

因子Aと因子Bに交互作用がない場合のデータを以下に示します。


   水準B1  水準B2  水準B3
 水準A1  3
 4
 5
 7
 8
 9
 4
 5
 6
 水準A2  5
 6
 7
 9
 10
 11
 6
 7
 8

「因子Aの水準と因子Bの水準の組み合わせ」ごとの平均値を求めると以下のようになります。


   水準B1  水準B2  水準B3
 水準A1  4  8  5
 水準A2  6  10  7

関数 interaction.plot() を使って、平均値のプロットを描いてみましょう。

関数 interaction.plot の書式は
interaction.plot(横軸に置く一つ目の因子, 軸外に置く2つ目の因子, 平均値を求めるデータ)
です。

ます、関数 interaction.plot() に入力する変数を作成します。


> データ <- c(3, 4, 5, 5, 6, 7, 7, 8, 9, 9, 10, 11, 4, 5,6, 6, 7, 8)


平均していない数値を、変数「データ」に格納します。


> 文字型の因子A <- rep(c(rep("水準A1", 3), rep("水準A2", 3)), 3) > 因子A <- factor(文字型の因子A)
> 因子A  [1] 水準A1 水準A1 水準A1 水準A2 水準A2 水準A2 水準A1 水準A1 水準A1 水準A2 [11] 水準A2 水準A2 水準A1 水準A1 水準A1 水準A2 水準A2 水準A2 Levels: 水準A1 水準A2


関数 rep(値, 回数) は、指定した値を指定した回数分、繰り返したデータを作成します。関数 rep() を利用して、数値の順番に対応するように水準A1と水準A2を並べ、変数「文字型の因子A」に格納します。次に、変数「文字型の因子A」を文字型データから因子型データに変換してから、変数「因子A」に格納します。最後に確認のため出力しました。


> 文字型の因子B <- c(rep("水準B1", 6), rep("水準B2", 6), rep("水準B3", 6)) > 因子B <- factor(文字型の因子B)
> 因子B  [1] 水準B1 水準B1 水準B1 水準B1 水準B1 水準B1 水準B2 水準B2 水準B2 水準B2 [11] 水準B2 水準B2 水準B3 水準B3 水準B3 水準B3 水準B3 水準B3 Levels: 水準B1 水準B2 水準B3


数値の順番に対応するように水準B1と水準B2と水準B3を並べてから、変数「文字型の因子B」に格納します。次に変数「文字型の因子B」を文字型データから因子型データに変換してから、変数「因子B」に格納します。最後に確認のため出力しました。


変数の下準備ができたので、関数 interaction.plot() を実行します。


> interaction.plot(因子A, 因子B, データ)


因子Aを横軸に取ります。平均値を結んだ3つの直線はすべて平行になっています。


> interaction.plot(因子B, 因子A, データ)



因子Bを横軸に取ります。平均値を結んだ2つの折れ線はすべて平行になっています。



因子Aと因子Bに交互作用がある場合のデータを以下に示します。


   水準B1  水準B2  水準B3
 水準A1  6
 7
 8
 7
 8
 9
 8
 9
 10
 水準A2  5
 6
 7
 10
 11
 12
 13
 14
 15


「因子Aの水準と因子Bの水準の組み合わせ」ごとの平均値を求めると以下のようになります。


   水準B1  水準B2  水準B3
 水準A1  7  8  9
 水準A2  6  11  14

関数 interaction.plot() を使って、平均値のプロットを描いてみましょう。まず、変数の下準備をします。


> データ <- c(6, 7, 8, 5, 6, 7, 7, 8, 9, 10, 11, 12, 8, 9,10, 13, 14, 15)


> 因子A <- factor(rep(c(rep("水準A1", 3), rep("水準A2", 3)), 3))
> 因子A
 [1] 水準A1 水準A1 水準A1 水準A2 水準A2 水準A2 水準A1 水準A1 水準A1 水準A2
[11] 水準A2 水準A2 水準A1 水準A1 水準A1 水準A2 水準A2 水準A2
Levels: 水準A1 水準A2


> 因子B <- factor(c(rep("水準B1", 6), rep("水準B2", 6), rep("水準B3", 6)))
> 因子B
 [1] 水準B1 水準B1 水準B1 水準B1 水準B1 水準B1 水準B2 水準B2 水準B2 水準B2
[11] 水準B2 水準B2 水準B3 水準B3 水準B3 水準B3 水準B3 水準B3
Levels: 水準B1 水準B2 水準B3


変数の下準備ができたので、関数 interaction.plot() を実行します。


> interaction.plot(因子A, 因子B, データ)



因子Aを横軸に取ります。平均値を結んだ3つの直線は平行になっていません。



> interaction.plot(因子B, 因子A, データ)






因子Bを横軸に取ります。平均値を結んだ2つの折れ線はすべて平行になっていません。


平均値を結んだ直線や折れ線が平行である場合、グラフで見た通りに交互作用効果はないと言えます。

しかし、平行でない場合、交互作用があるように見えますが、分散分析表を出力してみると、交互作用の効果はそれほどではなく、「有意ではない」と判定されることもあります。

ですので、交互作用効果があるかどうかは、きちんと分散分析表を出力して、最終確認する必要があります。




帰無仮説と対立仮説の設定

因子Aの主効果
帰無仮説 H0:因子Aの水準が違っても母平均は等しい(因子Aに主効果はない)
対立仮説 H1:因子Aの水準の違いにより母平均に差がある(因子Aには主効果がある)

因子Bの主効果
帰無仮説 H0:因子Bの水準が違っても母平均は等しい(因子Bに主効果はない)
対立仮説 H1:因子Bの水準の違いにより母平均に差がある(因子Bには主効果がある)

因子Aと因子Bの交互作用効果
帰無仮説 H0:因子Aの水準と因子Bの水準は互いに影響を及ぼしあっていない(因子Aと因子Bに交互作用はない)
対立仮説 H1:因子Aの水準と因子Bの水準は互いに影響を及ぼしあっている(因子Aと因子Bには交互作用がある)




Rによる統計解析

関数 aov() を使って、くり返し数の等しい二元配置分散分析を実行してみましょう。

関数 aov() の書式は
aov(データ~一つ目の因子*二つ目の因子)
です。「一つ目の因子*二つ目の因子」は、「一つ目の因子の主効果」と「二つ目の因子の主効果」と「一つ目の因子と二つ目の因子の交互作用」の3つすべてを含めるという意味です。

まず、関数 aov() に入力する3つの変数を作成しましょう。


> 血糖降下度 <- c(30, 25, 28,49, 59, 68, 80, 76, 82, 80, 65, 59, 79, 89, 78, 70, 86, 65)
> 血糖降下度
  [1] 30 25 28 49 59 68 80 76 82 80 65 59 79 89 78 70 86 65
> 種類 <- factor(c(rep("A", 9), rep("B", 9)))
> 種類
 [1] A A A A A A A A A B B B B B B B B B
Levels: A B
> 服用量 <- factor(rep(c(rep("100", 3), rep("500", 3), rep("1000", 3)), 2))
> 服用量
 [1] 100 100 100 500 500 500 1000 1000 1000 100 100 100 500 500
[15] 500 1000 1000 1000
Levels: 100 1000 500


関数 aov() を実行する準備が整ったので、くり返し数の等しい二元分散分析を実行してみましょう。


> summary(aov(血糖降下度~種類*服用量))

            Df Sum Sq Mean Sq F value Pr(>F)
種類 1 1682.0 1682.0 26.54 0.000241 ***
服用量 2 2732.1 1366.1 21.55 0.000107 ***
種類:服用量 2 1623.0 811.5 12.80 0.001056 **
Residuals 12 760.7 63.4
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


分散分析表のp値を見ると以下の結果が導かれます。

●種類:5%水準で有意な効果がある (p = 0.000241)
つまり、「対立仮説 H1:因子Aの水準の違いにより母平均に差がある(因子Aには主効果がある)」が採択されました。

●服用量:5%水準で有意な効果がある (p = 0.000107)
つまり、「対立仮説 H1:因子Bの水準の違いにより母平均に差がある(因子Bには主効果がある)」が採択されました。

●種類と服用量の交互作用:5%水準で有意な効果がある (p = 0.001056)
つまり、「対立仮説 H1:因子Aの水準と因子Bの水準は互いに影響を及ぼしあっている(因子Aと因子Bには交互作用がある)」が採択されました。


関数 interaction.plot() を使って、平均値のプロットを描いてみましょう。


> interaction.plot(種類, 服用量, 血糖降下度)



> interaction.plot(服用量, 種類, 血糖降下度)




上記の「主効果と交互作用」の章で載せた、交互作用がないデータと交互作用があるデータは、「因子Aの水準と因子Bの水準の組み合わせごとの平均値」を折れ線グラフにしただけで、本当のところ交互作用があるのかないのかはわからないままです。

それぞれのデータの分散分析表を出力してみましょう。



交互作用がないデータの分散分析表

> データ <- c(3, 4, 5, 5, 6, 7, 7, 8, 9, 9, 10, 11, 4, 5,6, 6, 7, 8)
> 因子A <- factor(rep(c(rep("水準A1", 3), rep("水準A2", 3)), 3))
> 因子B <- factor(c(rep("水準B1", 6), rep("水準B2", 6), rep("水準B3", 6)))
> summary(aov(データ~因子A*因子B))

            Df Sum Sq Mean Sq F value Pr(>F)
因子A 1 18 18 18 0.00114 **
因子B 2 52 26 26 4.35e-05 ***
因子A:因子B 2 0 0 0 1.00000
Residuals 12 12 1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


因子Aと因子Bの交互作用を表す「因子A:因子B 」の行のp値は1ですので、5%水準で有意な効果はありません。つまり、「帰無仮説 H0:因子Aと因子Bに交互作用はない」が採択されました。



交互作用があるデータの分散分析表

> データ <- c(6, 7, 8, 5, 6, 7, 7, 8, 9, 10, 11, 12, 8, 9,10, 13, 14, 15)
> 因子A <- factor(rep(c(rep("水準A1", 3), rep("水準A2", 3)), 3))
> 因子B <- factor(c(rep("水準B1", 6), rep("水準B2", 6), rep("水準B3", 6)))
> summary(aov(データ~因子A*因子B))

            Df Sum Sq Mean Sq F value Pr(>F)
因子A 1 24.5 24.5 24.5 0.000337 ***
因子B 2 76.0 38.0 38.0 6.43e-06 ***
因子A:因子B 2 28.0 14.0 14.0 0.000729 ***
Residuals 12 12.0 1.0
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’


因子Aと因子Bの交互作用を表す「因子A:因子B 」の行のp値は0.000729ですので、5%水準で有意な効果があります。つまり、「帰無仮説 H0:因子Aと因子Bに交互作用はない」は棄却され、「対立仮説 H1:因子Aと因子Bに交互作用がある」が採択されました。