一元配置の分散分析
分散分析は3群以上に分けられるデータに差があるかどうかを比較するときに行います。
データを分ける要因が一つである場合、一要因の分散分析を行います。
たとえば、「三年生の国語のテスト結果」をクラスごとに比較する場合、生徒のデータは「クラス」だけで区切られています。よって一要因の分散分析をもちいます。
要因:クラス | |||
---|---|---|---|
A組 | B組 | C組 | |
国語の点数 | 生徒1 : 67点 生徒2 : 63点 生徒3 : 77点 生徒4 : 70点 生徒5 : 80点 生徒6 : 97点 生徒7 : 85点 |
生徒8 : 45点 生徒9 : 65点 生徒10 : 87点 生徒11 : 46点 生徒12 : 88点 生徒13 : 48点 生徒14 : 55点 |
生徒15 : 64点 生徒16 : 67点 生徒17 : 62点 生徒18 : 58点 生徒19 : 60点 生徒20 : 87点 生徒21 : 75点 |
平均点 | 77.0 | 62.0 | 67.6 |
- 国語の点数は各生徒から一度だけ測定されているので「データにくり返しはない」といいます。
- データ数 n = 21
等分散性を仮定した一要因の分散分析
class.A <- c(67, 63, 77, 70, 80, 97, 85) class.B <- c(45, 65, 87, 46, 88, 48, 55) class.C <- c(64, 67, 62, 58, 60, 87, 75) #データの入力 dat <- data.frame(class.A,class.B,class.C) # 上記表とおなじく行に生徒の点数、列にクラスをとったwide形式の表 dat <- stack(dat) # Rが扱いやすいstack形式の表に変換する # ind 列にクラス名(要因の違い)、values 列に点数(従属変数)が並ぶ # このとき各行は生徒個人を表しています。 oneway.test(data=dat, values~ind, var.equal=T) # 繰り返しのない一要因の分散分析(等分散性を仮定)を実行する # データ dat の values を ind の違いによって比較
- Rでは ~ (チルダ)を使って Y~Xとかくと「YをXによって」という意味を示す
- oneway.test()
- 一要因の分散分析を実行する。var.equal=T で各群の分散は等しいと仮定。Fで等分散性を仮定しないWelchの方法を実行する。
- 大量のデータを入力する場合、エクセル、テキストファイルで表をつくっておいてread.table()で読み込む。
- スタック形式の表について
Oneway.test(var.equal=T)の結果の見方
One-way analysis of means
data: values and ind
F = 2.0441, num df = 2, denom df = 18, p-value = 0.1585
- F: F value、F値 分散分析の統計量
- num df: numenator df、分子の自由度、要因の自由度
- denom df: denominator df、分母の自由度、多くの場合「残差の自由度」
- p-value:p値、num df, denom dfのF分布表でのF値の有意確率
分散分析では、要因の自由度、残差の自由度から決定されるF分布表をもちいて、算出されたF値の有意確立pを求める。したがって、p値を示すときには F(要因の自由度, 残差の自由度)= F値 もあわせて示す。
p値が有意水準α(通常0.05)より小さいとき、values(点数)はind(クラス)によって異なる、といえる。
ここでは p > 0.05 であるため有意ではない。
したがって「国語の点数はクラスごとに差があるとは言えない」*1。
繰り返しのない一要因の分散分析(等分散性を仮定)の結果は、関数 aov()、lm()でも同じものが得られる。
comp.test <- aov(data=dat, values~ind) summary(comp.test) comp.test.2 <- lm(data=dat, values~ind) anova(comp.test.2)
結果は下記のように分散分析表で表示される。
Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
ind | 2 | 48.1 | 24.07 | 0.18 | 0.837 |
Residuals | 12 | 1601.6 | 133.47 |
等分散性のチェック
等分散でないデータに等分散性を仮定した検定を行うと、誤った結果を得る可能性がある。したがって、Oneway.test(var.equal=T)を使う際には、データが等分散であることを確認しておく必要がある。
データのばらつきを見るにはボックスプロット(箱ひげ図)を作ってみる。
boxplot(data=dat, values~ind)
B組のばらつきがやや大きいように見える。
このデータが統計的に等分散といえるかどうかはバートレット検定を用いて確認できる。
bartlett.test(data=dat, values~ind)
結果はこのように表示される。
Bartlett test of homogeneity of variances
data: values by ind
Bartlett's K-squared = 2.3881, df = 2, p-value = 0.303
バートレット検定では帰無仮説H0は「各群の分散は等しい」であり、p > 有意水準α(0.05)のとき、帰無仮説は棄却できない。したがって、「等しくない、とはいえない」と結論付けられる。*2
このデータでは p = 0.30 > α = 0.05 であるので、等しくないとは言えない。つまり等分散性を仮定した分散分析を用いても悪くない。
いっぽう、p < 0.05のとき帰無仮説は棄却され、「等しくない」という対立仮説H1を採択する。
このとき、等分散性を仮定しないWelchの方法を使うべきである。
一要因の分散分析 Welchの方法
Welchの方法は等分散であってもなくても利用できる。*3
たとえば、下記のようなデータではクラスごとの点数のばらつきは等分散とはいえない。
class.D <- c(67, 70, 77, 70, 80, 70, 85) #平均 77.0 class.E <- c(45, 40, 55, 80, 92, 48, 55) #平均 62.0 class.F <- c(64, 67, 62, 58, 60, 87, 75) #平均 67.6 dat2 <- stack(data.frame(class.D, class.E, class.F)) bartlett.test(data=dat2, values~ind)
Bartlett test of homogeneity of variances
data: values by ind
Bartlett's K-squared = 6.2796, df = 2, p-value = 0.04329
ボックスプロット(箱ひげ図)をみるとE組の点数のばらつきがD、F組より大きいことが分かる。
boxplot(data=dat2, values~ind)
したがって、等分散を仮定した一要因の分散分析ではなく、等分散を仮定しないWelchの方法で一要因の分散分析をおこなうべきである。
oneway.test(data=dat2, values~ind, var.equal=F)
One-way analysis of means (not assuming equal variances)
data: values and ind
F = 2.3196, num df = 2.000, denom df = 10.642, p-value = 0.1458
等分散性を仮定した方法とは、denom dfおよびF値の算出方法が異なる。*4
のでlm()およびaov()では実行できないし、分散分析表も作らない。しかしF値はnum dfとdenom dfに従うのでF(num df, denom df)=F値 と書く。
ノンパラメトリック法による一元配置分散分析:クラスカル・ウォリス法
データの分布が正規分布に当てはまらないときは、ノンパラメトリックな手法を用いることが適切である。
データの分布を確認するには、ヒストグラムと確率密度曲線を描いて判断することが簡便。
hist(dat$values) par(new=T) #グラフへの上書き許可 plot(density(dat$values),xlim=c(40,100),axes=F,ann=F) # xlim には表示されたヒストグラムのx軸下限、上限を指定する。 # axes=F 軸の非表示、ann=F ラベルの非表示 hist(dat2$values) par(new=T) plot(density(dat2$values),xlim=c(40,100),axes=F,ann=F)
どちらも大まかにいって正規分布にあてはまっているだろう。
日常的なデータの多くは正規分布すると仮定されている。
しかし、たとえば毎日の欠席者数のように0からカウントアップされ、頻度が多くないデータでは、その分布は正規分布に従わないことが多い。
要因:クラス | |||
---|---|---|---|
A組 | B組 | C組 | |
十日間の遅刻者数 | 1日目 : 1人 2日目 : 0人 3日目 : 1人 4日目 : 0人 5日目 : 0人 6日目 : 2人 7日目 : 1人 8日目 : 1人 9日目 : 0人 10日目 : 0人 |
1日目 : 3人 2日目 : 1人 3日目 : 2人 4日目 : 3人 5日目 : 1人 6日目 : 2人 7日目 : 1人 8日目 : 2人 9日目 : 1人 10日目 : 1人 |
1日目 : 5人 2日目 : 2人 3日目 : 2人 4日目 : 1人 5日目 : 1人 6日目 : 3人 7日目 : 2人 8日目 : 1人 9日目 : 1人 10日目 : 0人 |
平均人数 | 0.6 | 1.7 | 1.8 |
- n = 30 ここでは1~10日目を独立のデータとして扱う。*5
class.A <- c(1, 0, 1, 0, 0, 2, 1, 1, 0, 0) class.B <- c(3, 1, 2, 3, 1, 2, 1, 2, 1, 1) class.C <- c(5, 2, 2, 1, 1, 3, 2, 1, 1, 0) dat3 <- stack(data.frame(class.A, class.B, class.C)) hist(dat3$values)
ヒストグラムから、明らかに正規分布していないことが分かる。
このような場合、正規性を仮定しないノンパラメトリックな一元配置分散分析として、クラスカル・ウォリス検定を用いる。
kruskal.test(data=dat3, values~ind)
Kruskal-Wallis rank sum test
data: values by ind
Kruskal-Wallis chi-squared = 8.4607, df = 2, p-value = 0.01455
p < 0.05 であり、3組の遅刻数には有意な差があると言える。