山羊の午後

研究関係の備忘録。

一元配置の分散分析

分散分析は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.062.067.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()で読み込む。
  • スタック形式の表について

yaginogogo.hatenablog.jp


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
  • Df:自由度、Sum Sq 平方和、Mean Sq 平均平方、F値、p値
  • Residuals 残差
  • Mean Sq = Sum Sq/Df
  • F値 = 要因の平均平方/残差の平均平方

等分散性のチェック

等分散でないデータに等分散性を仮定した検定を行うと、誤った結果を得る可能性がある。したがって、Oneway.test(var.equal=T)を使う際には、データが等分散であることを確認しておく必要がある。

データのばらつきを見るにはボックスプロット(箱ひげ図)を作ってみる。

boxplot(data=dat, values~ind)

f:id:yk_uminami:20170225005002p:plain

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)

f:id:yk_uminami:20170225003700p:plain

したがって、等分散を仮定した一要因の分散分析ではなく、等分散を仮定しない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)

f:id:yk_uminami:20170225003817p:plain

どちらも大まかにいって正規分布にあてはまっているだろう。

日常的なデータの多くは正規分布すると仮定されている。
しかし、たとえば毎日の欠席者数のように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.61.71.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)

f:id:yk_uminami:20170225003849p:plain

ヒストグラムから、明らかに正規分布していないことが分かる。
このような場合、正規性を仮定しないノンパラメトリックな一元配置分散分析として、クラスカル・ウォリス検定を用いる。

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組の遅刻数には有意な差があると言える。

どれをつかうのかチャート

データは正規分布している 等分散である Oneway.test(var.equal=T)
等分散でない(又はわからない) Oneway.test(var.equal=F)
正規分布していない kruskal.test()

(注)あくまでも目安です。

*1:「差がない」とは言えない。

*2:「等しい」と断言してはいけない。

*3:そのため、最初からWelchの方法を使うべきとする向きもある。

*4:参考:一元配置分散分析 の下部

*5:しかし実際にはこのようなデータは、同じ生徒、同じ曜日、などの繰り返しを考慮すべきだろう。