人工知能を作りたい!

統計的仮説検定

統計的仮説検定の使い方




母平均に関する t 検定(1群、母分散が未知)

はじめに

母集団分布が正規分布、母分散σ2が既知の場合、検定統計量 Z を求める数式は以下のようになります。

Z = Z 分子/ Z 分母
Z 分子:標本平均- μ
Z 分母: ( σ2 /サンプルデータ数 n )の平方根
( σ2 の値は既知の値を使う)

しかし、母集団分布が正規分布、母分散σ2 が未知の場合、検定統計量 Z の値を計算できません。そこで、検定統計量 t を使います。tを求める数式は以下のようになります。

t = t 分子/ t 分母
t 分子:標本平均-μ
t 分母: ( σ2 /サンプルデータ数 n )の平方根
( σ2 の値は不偏分散の値を使う)

不偏分散は「サンプルデータと標本平均の偏差の二乗和」を「サンプルデータ数 n -1」で割った値です。

不偏分散 = 不偏分散の分子 / 不偏分散の分母
不偏分散の分子:(データ1-標本平均)2+(データ2-標本平均)2+ . . . +(データ n -標本平均)2
不偏分散の分母:サンプルデータ数 n -1

ちなみに、標本分散は「サンプルデータと標本平均の偏差の二乗和」を「サンプルデータ数 n 」で割った値です。

標本分散 = 標本分散の分子/標本分散の分母
標本分散の分子:(データ1-標本平均)2+(データ2-標本平均)2+ . . . +(データ n -標本平均)2
標本分散の分母:サンプルデータ数 n

母数の推定量である不偏推定量を優先するなら不偏分散を使います。最尤(さいゆう)推定量を優先するなら標本分散を使います。

検定統計量 t は自由度 df (degrees of freedom) = n -1の t 分布に従います。それでは、 t 分布を用いた検定( t 検定と言います)により一つの平均値の検定を行います。




例題(両側検定の場合)

新しい数学習得法のもとで数学を学んだ学生群から無作為に20人を選び、数学テストを受けてもらいました。その得点データは以下のようになりました。

41, 52, 60, 73, 32, 53, 62, 70, 81, 48, 55, 62, 45, 85, 27, 72, 66, 77, 58, 40

過去の数年のデータの蓄積から、従来の数学習得法を受けた学生の数学テスト得点は、母平均 μ =40の正規分布に従うことが知られています。つまり、母集団分布は X~N(40, σ2) です。

新しい数学習得法を受けた20人の数学テスト得点の母平均は40から変化しているでしょうか。




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

帰無仮説 H0: μ =40(数学テストの母平均は40である)
対立仮説 H1: μ ≠40(数学テストの母平均は40ではない)
有意水準 α =0.05
両側検定である




Rによる統計解析

Rに用意されている関数 t.test() を使って t 検定を行います。

関数 t.test() の書式は
t.test(検定対象となるデータ, mu=帰無仮説に設定された平均値)
です。


> x <- c(41, 52, 60, 73, 32, 53, 62, 70, 81, 48, 55, 62, 45, 85, 27, 72, 66, 77, 58, 40)


数学テストの得点を変数「x」に格納します。


> t.test(x,mu=40)

        One Sample t-test

data: x
t = 5.0126, df = 19, p-value = 7.73e-05
alternative hypothesis: true mean is not equal to 40
95 percent confidence interval:
 50.45488 65.44512
sample estimates:
mean of x
    57.95


「One Sample t-test」は「一つのサンプルの t 検定」というタイトルです。

「data:」の横に検定の対象となる変数が示されます。下の行に行くと、検定統計量 t の実現値、t分布の自由度 ( df = n -1=20-1=19) 、 p 値(デフォルトでは両側検定を行う)が示されます。

p 値が7.73e-05、つまり、0.0000773となり、有意水準0.05より小さいので「帰無仮説 H0: μ =40(数学テストの母平均は40である)」は棄却され、検定の結果は有意となります。

「alternative hypothesis: true mean is not equal to 40」は「対立仮説 H1: μ ≠40(数学テストの母平均は40ではない)」を示しています。これより、両側検定であることが分かります。

「95 percent confidence interval:」は、95%信頼区間の下限と上限の値を示しています。「sample estimates:」は、標本から計算された標本平均の実現値を示しています。




例題(片側検定の場合)

新しい数学習得法のもとで数学を学んだ学生群から無作為に20人を選び、数学テストを受けてもらいました。その得点データは以下のようになりました。

41, 52, 60, 73, 32, 53, 62, 70, 81, 48, 55, 62, 45, 85, 27, 72, 66, 77, 58, 40

過去の数年のデータの蓄積から、従来の数学習得法を受けた学生の数学テスト得点は、母平均 μ =40の正規分布に従うことが知られています。つまり、母集団分布は X~N(40, σ2) です。

新しい数学習得法を受けた20人の数学テスト得点の母平均は40より大きくなっているでしょうか。




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

帰無仮説 H0: μ =40(数学テストの母平均は40である)
対立仮説 H1: μ >40(数学テストの母平均は40より大きい)
有意水準 α =0.05
片側検定である




Rによる統計解析

Rに用意されている関数 t.test() を使って t 検定を行います。

関数 t.test() の書式は
t.test(検定対象となるデータ, mu=帰無仮説に設定されたμの値, alternative="greater")
です。

デフォルトでは両側検定が行われますが、以下のようにオプションにて対立仮説の方向を指定することができます。
上片側検定: alternative="greater"
下片側検定: alternative="less"


> x <- c(41, 52, 60, 73, 32, 53, 62, 70, 81, 48, 55, 62, 45, 85, 27, 72, 66, 77, 58, 40)


数学テストの得点を変数「x」に格納します。


> t.test(x, mu=40, alternative="greater")

        One Sample t-test

data: x
t = 5.0126, df = 19, p-value = 3.865e-05
alternative hypothesis: true mean is greater than 40
95 percent confidence interval:
 51.75798 Inf
sample estimates:
mean of x
    57.95


「One Sample t-test」は「一つのサンプルの t 検定」というタイトルです。

「data:」の横に検定の対象となる変数が示されます。下の行に行くと、検定統計量 t の実現値、 t 分布の自由度 ( df = n -1=20-1=19) 、 p 値(上片側検定を設定しました)が示されます。

p 値が3.865e-05、つまり、0.00003865となり、有意水準0.05より小さいので「帰無仮説 H0: μ =40(数学テストの母平均は40である)」は棄却され、検定の結果は有意となります。

「alternative hypothesis: true mean is greater than 40」は「対立仮説 H1: μ >40(数学テストの母平均は40より大きい)」を示しています。これより、上片側検定であることが分かります。

「95 percent confidence interval:」は、95%信頼区間の下限と上限の値を示しています。「sample estimates:」は、標本から計算された標本平均の実現値を示しています。






母分散の等質性に関するF検定、母平均の差に関する t 検定(2群、対応なし)、Welch(ウェルチ)の検定

はじめに

t 検定(2群、対応なし)を実行するためには、以下に示す3つの前提条件を満たす必要があります。

●標本抽出が無作為に行われている。(無作為抽出)
●母集団分布が正規分布に従っている。(正規性)
●2つの群の母分散が等質である。(母分散の等質性)

2つの群の母分散が等質でない場合、Welch(ウェルチ)の検定という特殊な t 検定を行う必要があります。




例題(母分散が等質である場合)

英語テストを作成し、全国から無作為に抽出した高校生1年生男子10人と高校一年生女子10人にその英語テストを受けてもらいました。そして、以下のような英語テスト得点のデータを得ました。

<男子の英語テスト得点>
53, 37, 49, 72, 43, 62, 41, 52, 76, 65

<女子の英語テスト得点>
45, 57, 68, 35, 54, 48, 85, 40, 55

2つの母集団(全国の高校1年生男子と全国の高校1年生女子)それぞれの英語テスト平均点に有意な差があるといえるでしょうか。




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

「2つの群ぞれぞれの平均点」とは、「2つの群ぞれぞれの母平均」のことです。

帰無仮説 H0: μ1 = μ2 (2つの群それぞれの母平均は等しい)
対立仮説 H1: μ1 ≠ μ2 (2つの群それぞれの母平均は異なる)
有意水準 α =0.05




Rによる統計解析

<母分散の等質性の検定>
まず、母分散の等質性の検定を、関数 ver.test() を使って行います。

関数 ver.test() の書式は
ver.test(群1のデータ, 群2のデータ)
です。

母分散の等質性の検定において、帰無仮説と対立仮説と有意水準は以下のようになります。

帰無仮説 H0:2つの群それぞれの母分散は等しい
対立仮説 H1:2つの群それぞれの母分散は異なる
有意水準 α =0.05


> 男 <- c(53, 37, 49, 72, 43, 62, 41, 52, 76, 65)
> 女 <- c(45, 57, 68, 35, 54, 48, 85, 40, 55)


男子の得点データを変数「男」に、女子の得点データを変数「女」に格納します。


> var.test(男, 女)
        F test to compare two variances

data: 男 and 女
F = 0.76871, num df = 9, denom df = 8, p-value = 0.6996
alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval:
 0.1764218 3.1532175 sample estimates:
ratio of variances
         0.7687108


「F test to compare two variances」は「2つの分散の比較のためのF検定」というタイトルです。

「data:」の横に検定の対象となる変数が「男 and 女」と示されます。下の行に行くと、検定統計量Fの実現値、F分布の自由度 (「num df」が分子の自由度、「denom df」が分母の自由度) 、p 値が示されます。

p 値が0.6996となり、有意水準0.05より大きいので帰無仮説は棄却されません。つまり、「帰無仮説 H0:2つの群それぞれの母分散は等しい」が採択されました。

「alternative hypothesis: true ratio of variances is not equal to 1」は、「対立仮説:母分散の比率が1(比が1対1)ではない」を示しています。つまり、「対立仮説 H1:2つの群それぞれの母分散は異なる」です。

「95 percent confidence interval:」は、95%信頼区間の下限と上限の値を示しています。「sample estimates:」は、標本から計算された不偏分散の比を示しています。



<t 検定(2群、対応なし)>
それでは、 t 検定の前提である「母分散の等質性」は満たされているということなので、関数 t.test() を使って、t 検定(2群、対応なし)を実行します。

関数 t.test() の書式は
t.test(群1のデータ, 群2のデータ, var.equal=TRUE)
です。

var.equal=TRUE がないデフォルトでは、Welch(ウェルチ)の検定が実行されます。ですので、 t 検定(2群、対応なし)を実行するには var.equal=TRUE というオプションを指定する必要があります。


> t.test(男, 女, var.equal=TRUE)

        Two Sample t-test

data: 男 and 女
t = 0.13614, df = 17, p-value = 0.8933
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -12.88655 14.66433
sample estimates:
mean of x mean of y
 55.00000 54.11111


「Two Sample t-test」は、「2つのサンプルの t 検定」というタイトルです。

「data:」の横に検定の対象となる変数が「男 and 女」と示されます。下の行に行くと、検定統計量 t の実現値、 t 分布の自由度、 p 値が示されます。

p 値が0.8933となり、有意水準0.05より大きいので、「帰無仮説 H0: μ1 = μ2 (2つの群それぞれの母平均は等しい)」が採択されました。

「alternative hypothesis: true difference in means is not equal to 0」は、「対立仮説: μ1 ≠ μ2 」、つまり「対立仮説: μ1 - μ2 ≠ 0」を示しています。これより両側検定であることが分かります。

「95 percent confidence interval:」は、95%信頼区間の下限と上限の値を示しています。「sample estimates:」は、標本から計算された標本平均の実現値を示しています。




例題(母分散が等質でない場合)

英語テストを作成し、全国から無作為に抽出した高校生1年生男子10人と高校一年生女子10人にその英語テストを受けてもらいました。そして、以下のような英語テスト得点のデータを得ました。

<男子の英語テスト得点>
53, 37, 49, 72, 43, 62, 41, 52, 76, 65

<女子の英語テスト得点>
25, 57, 68, 15, 54, 88, 85, 20, 75

2つの母集団(全国の高校1年生男子と全国の高校1年生女子)それぞれの英語テスト平均点に有意な差があるといえるでしょうか。




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

「2つの群ぞれぞれの平均点」とは、「2つの群ぞれぞれの母平均」のことです。

帰無仮説 H0: μ1 = μ2 (2つの群それぞれの母平均は等しい)
対立仮説 H1:μ1 ≠ μ2(2つの群それぞれの母平均は異なる)
有意水準 α=0.05




Rによる統計解析

<母分散の等質性の検定>
まず、母分散の等質性の検定を、関数 ver.test() を使って行います。

関数 ver.test() の書式は
ver.test(群1のデータ, 群2のデータ)
です。

母分散の等質性の検定において、帰無仮説と対立仮説と有意水準は以下のようになります。

帰無仮説 H0:2つの群それぞれの母分散は等しい
対立仮説 H1:2つの群それぞれの母分散は異なる
有意水準 α =0.05


> 男 <- c(53, 37, 49, 72, 43, 62, 41, 52, 76, 65)
> 女 <- c(25, 57, 68, 15, 54, 88, 85, 20, 75)


男子の得点データを変数「男」に、女子の得点データを変数「女」に格納します。


> var.test(男, 女)

        F test to compare two variances

data: 男 and 女
F = 0.2253, num df = 9, denom df = 8, p-value = 0.03922
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.05170811 0.92418815
sample estimates:
ratio of variances
         0.2253043


p 値が0.03922となり、有意水準0.05より小さいので帰無仮説は棄却されます。つまり、「対立仮説 H1:2つの群それぞれの母分散は異なる」が採択されました。



<Welch(ウェルチ)の検定>
t 検定の前提である「母分散の等質性」は満たされていないということなので、t 検定(2群、対応なし)に代えて、Welch(ウェルチ)の検定を実行します。

Welch(ウェルチ)の検定を行う関数 t.test の書式は
t.test(群1のデータ, 群2のデータ)
です。

var.equal=TRUE がないデフォルトでは、Welch(ウェルチ)の検定が実行されます。


> t.test(男, 女)

        Welch Two Sample t-test

data: 男 and 女
t = 0.086778, df = 11.165, p-value = 0.9324
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -21.61563 23.39341
sample estimates:
mean of x mean of y
 55.00000 54.11111


「Welch Two Sample t-test」は「ウェルチの2つのサンプルの t 検定」というタイトルです。

p 値が0.9324となり、有意水準0.05より大きいので帰無仮説は棄却されません。つまり、「帰無仮説 H0: μ1 = μ2 (2つの群それぞれの母平均は等しい)」が採択されました。




まとめ

t 検定(2群、対応なし)を行う場合、まず、関数 var.test(群1のデータ, 群2のデータ) を実行して「母分散の等質性」をチェックします。

母分散が等質である場合、関数 t.test(群1のデータ, 群2のデータ, var.equal=TRUE) を使って、t 検定(2群、対応なし)を実行します。

母分散が等質でない場合、関数 t.test(群1のデータ, 群2のデータ) を使って、Welchの検定を実行します。






母平均の差に関する t 検定(2群、対応あり)

はじめに

母平均に差をもたらす原因を因子といいます。また、因子の中に含まれる個々の条件を水準といいます。

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

ここでちょっと混乱が起こります。 t 検定(2群)も一元配置分散分析も、母平均差の検定を行います。そして、両方とも因子は1つです。では、t 検定(2群)と一元配置分散分析の違いは何でしょうか。

一元配置分散分析は「3つ以上の水準(標本群)の母平均差」を検定します。しかし、t 検定(2群)は「2つの水準(標本群)の母平均差」しか検定できません。つまり、まとめると以下のようになります。

●2つの標本群の母平均差を検定する場合 → t 検定(2群)
●3つ以上の標本群の母平均差を検定する場合 → 一元配置分散分析




例題

飲料会社Aは体脂肪を減らす効果があると思われるお茶Aを開発しました。17人の従業員を無作為に抽出し、一か月間、毎朝飲んでもらいました。その結果、実験開始前と一か月後の体脂肪率 (%) は以下のように変化しました。


 実験開始前  一か月後
 23.5  22.5
 27.5  27.3
 34.2  32.4
 14.6  15.0
 30.45  27.5
 27.2  25.4
 31.3  29.3
 23.7  24.4
 21.5  20.6
 16.5  17.6
 19.5  20.2
 29.4  26.5
 17.4  16.4
 20.6  21.4
 23.6  21.3
 32.3  30.1
 13.4  14.3

このデータから、お茶Aは体脂肪を減らす効果があると判断してよいでしょうか。




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

この検定は、 t 検定(2群、対応あり)ですが、一元分散分析(対応あり)と同じように、「因子と水準は何か」を考えると以下のようになります。

繰り返しますが、一元分散分析(対応あり)は3つ以上の水準(あるいは標本群)を持ちます。

1つの因子:
お茶Aの摂取

2つの水準(標本群)
水準(群)1:「お茶Aを飲まない」
水準(群)2:「お茶Aを飲んだ」

実験開始前の母平均を μA 、一か月後の母平均を μB とします。

帰無仮説 H0: μA = μB (実験開始前の母平均と一か月後の母平均に差はない)
対立仮説 H1: μA > μB (実験開始前の母平均は一か月後の母平均より大きい)
有意水準 α =0.05
片側検定となります。




Rによる統計解析

関数 t.test() を使って、母平均の差に関する t 検定(2群、対応あり)を実行します。

関数 t.test() の書式は
t.test(群1のデータ, 群2のデータ, paired=TRUE, alternative="greater") です。

2つの変数が対応のある場合、 paired=TRUE とオプションを指定します。

「対立仮説 H1: μA > μB 」は
「対立仮説 H1: μA - μB >0」つまり、
「対立仮説 H1: μA - μB greater than 0」と言い換えられます。ですので、 alternative="greater" とオプション指定をします。片側検定となります。


> 実験開始前 <-c(23.5, 27.5, 34.2, 14.6, 30.4, 27.2, 31.3, 23.7, 21.5, 16.5, 19.5, 29.4, 17.4, 20.6, 23.6, 32.3, 13.4)



実験開始前の体脂肪率のデータを変数「実験開始前」に格納します。


> 一か月後 <-c(22.5, 27.3, 32.4, 15.0, 27.5, 25.4, 29.3, 24.4, 20.6, 17.6, 20.2, 26.5, 16.4, 21.4, 21.3, 30.1, 14.3)


一か月後の体脂肪率のデータを変数「一か月後」に格納します。


> t.test(実験開始前, 一か月後, paired=TRUE, alternative="greater")

        Paired t-test

data: 実験開始前 and 一か月後
t = 2.4741, df = 16, p-value = 0.01247
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
 0.2493275 Inf
sample estimates:
mean of the differences
              0.8470588


「Paired t-test」は、「対応のある t 検定」というタイトルです。

「data:」の横に検定の対象となる変数が「実験開始前 and 一か月後」と示されます。下の行に行くと、検定統計量 t の実現値、 t 分布の自由度 df ( n -1=17-1=16)、p値が示されます。

p 値が0.01247となり、有意水準0.05より大きいので帰無仮説は棄却されます。つまり、「対立仮説 H1: μA > μB (実験開始前の母平均は一か月後の母平均より大きい)」が採択されました。

「alternative hypothesis: true difference in means is greater than 0」は、
「対立仮説 H1: μA > μB 」、つまり、「対立仮説 H1: μA - μB >0」を示しています。これより片側検定であることが分かります。

「95 percent confidence interval:」は、95%信頼区間の下限と上限の値を示しています。「sample estimates:」は、標本から計算された差の平均値を示しています。






独立性の検定(カイ二乗検定)

はじめに

数値という量的変数どうしの関係を調べる場合、散布図を描き、相関係数の検定無相関検定)を行います。

一方、「好き嫌い」という質的変数どうしの関係を調べる場合、クロス集計表を描き、独立性の検定カイ二乗検定)を行います。

質的変数どうしに関係がある時、「2つの質的変数に連関がある」、あるいは「2つの質的変数は独立ではない」と表現します。

質的変数どうしに関係がない時、「2つの質的変数に連関はない」、あるいは「2つの質的変数は独立である」と表現します。




例題

高校生20人にチョコレートとコーヒーの好き嫌いを聞き、以下の結果を得ました。(高校生一人一人に通し番号を付け、以下のように番号順に並べました。)

<チョコレートの好き嫌い>
"好き", "好き", "嫌い", "好き", "好き", "嫌い", "嫌い", "嫌い", "好き", "嫌い", "好き", "嫌い", "好き", "好き", "好き", "嫌い", "嫌い", "好き", "嫌い", "好き"

<コーヒーの好き嫌い>
"嫌い", "嫌い", "好き", "好き", "好き", "嫌い", "嫌い", "嫌い", "好き", "嫌い", "好き", "嫌い", "好き", "嫌い", "嫌い", "好き", "嫌い", "好き", "嫌い", "嫌い"

「チョコレートの好き嫌い」と「コーヒーの好き嫌い」の間に有意な連関があるといえるでしょうか。有意水準5%で検定を行ってください。




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

帰無仮説 H0:2つの変数「チョコレートの好き嫌い」と「コーヒーの好き嫌い」は独立である(連関がない)
対立仮説 H1:2つの変数「チョコレートの好き嫌い」と「コーヒーの好き嫌い」は独立ではない(連関がある)




Rによる統計解析

> チョコレート <- c("好き", "好き", "嫌い", "好き", "好き", "嫌い", "嫌い", "嫌い", "好き", "嫌い", "好き", "嫌い", "好き", "好き", "好き", "嫌い", "嫌い", "好き", "嫌い", "好き")
> table(チョコレート)
チョコレート
嫌い 好き
     9 11


見やすいように修正して以下に表示します。


 チョコレート嫌い  チョコレート好き
 9  11


変数「チョコレート」にデータを格納します。関数 table() を使って度数分布を確認しました。


> コーヒー <- c("嫌い", "嫌い", "好き", "好き", "好き", "嫌い", "嫌い", "嫌い", "好き", "嫌い", "好き", "嫌い", "好き", "嫌い", "嫌い", "好き", "嫌い", "好き", "嫌い", "嫌い")
> table(コーヒー)
コーヒー
嫌い 好き
   12 8


見やすいように修正して以下に表示します。


 コーヒー嫌い  コーヒー好き
 12  8


変数「コーヒー」にデータを格納します。関数 table() を使って度数分布を確認しました。


次に、変数「チョコレート」と変数「コーヒー」を組み合わせて、関数 table() を使ってみましょう。


> table(チョコレート, コーヒー)
                       コーヒー
チョコレート 嫌い 好き
               嫌い 7 2
               好き 5 6


見やすいように修正して以下に表示します。


   コーヒー嫌い  コーヒー好き
 チョコレート嫌い  7  2
 チョコレート好き  5  6




関数 chisq.test() を使ってカイ二乗検定を実行してみましょう。

関数 chisq.test() の書式は
chisq.test(クロス集計表, correct=FALSE)
です。

関数 table で作成したクロス集計表の結果を引数として指定します。correct=FALSE というオプションを指定すると、デフォルトで実行される連続性の補正を行いません。


> クロス集計表 <- table(チョコレート, コーヒー)


まず最初に、変数「クロス集計表」に、上記の関数 table(チョコレート, コーヒー) の結果を格納します。


> chisq.test(クロス集計表, correct=FALSE)
   Pearson's Chi-squared test

data: クロス集計表
X-squared = 2.1549, df = 1, p-value = 0.1421

 警告メッセージ:
 chisq.test(クロス集計表, correct = FALSE) で:
   カイ自乗近似は不正確かもしれません


「Pearson's Chi-squared test」は「ピアソンのカイ二乗検定」というタイトルです。

「data:」の横に検定の対象となる変数である「クロス集計表」が示されます。下の行に行くと、検定統計量 χ2 の実現値、カイ二乗分布の自由度、 p 値が示されます。

自由度は
(行の数-1)×(列の数-1)
で求められます。2×2のクロス集計表では
(2-1)×(2-1)=1
となります。

p 値が0.1421となり、有意水準0.05より大きいので、「帰無仮説 H0:2つの変数「チョコレートの好き嫌い」と「コーヒーの好き嫌い」は独立である(連関がない)」は棄却されません。

「chisq.test(クロス集計表, correct = FALSE) で: カイ自乗近似は不正確かもしれません」という警告メッセージは、データから計算される期待度数が小さい時に発せられます。カイ二乗検定の結果が不正確である可能性を教えてくれます。






一元配置分散分析(対応なし)、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に交互作用がある」が採択されました。