TukeyHSD法で多重比較を行い、a, b, cなどの文字でグループ間の有意差を表す。
> res<-aov(value~fac,data) > summary(res)
もし、fac因子の水準間に有意差があれば
> TukeyHSD(res, "fac")
とする。牛の乳量のデータでは
> data<-read.table("cow.csv", sep=",", header=T)
> res<-aov(Yield~Food+Cow+Season,data)
> summary(res)
# Foodの水準間でのみ有意差がある
> TukeyHSD(res, "Food")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Yield ~ Food + Cow + Season, data = data)
$Food
diff lwr upr p adj
B-A 0.50 -1.230856 2.23085569 0.7553967
C-A -1.75 -3.480856 -0.01914431 0.0478429
D-A -1.25 -2.980856 0.48085569 0.1571372
C-B -2.25 -3.980856 -0.51914431 0.0160634
D-B -1.75 -3.480856 -0.01914431 0.0478429
D-C 0.50 -1.230856 2.23085569 0.7553967
Foodの水準ごとの平均は
> by(data$Yield, data$Food, mean) data$Food: A [1] 9.25 ------------------------------------------------------ data$Food: B [1] 9.75 ------------------------------------------------------ data$Food: C [1] 7.5 ------------------------------------------------------ data$Food: D [1] 8
これを表にまとめる。
| 餌 | 乳量 (kg) | |
|---|---|---|
| A | 9.25 | ab |
| B | 9.75 | a |
| C | 7.50 | c |
| D | 8.00 | bc |
まず、一番値が大きいBにアルファベットaを振り、それと有意差のない水準に同じaをつけていく。TukeyHSD関数の結果から、AとBの間には5%水準で有意な差がないので、Aの数字の横にもBと同じaを振る。次に、２番目に大きいAにbを振り、すでに比較が終わっているB以外の水準を見ていくと、AとCの間には有意差があり、AとDの間にはないので、Dの数字の横にbを振る。３番目のDにcをつけ、B、Aとの比較は終わっているので、Cとの比較をすると、有意差はないのでCにcをつける。 もし、最後のD-C間に有意差があれば、Cにdを振るが、この場合、Dにつけたcが無意味になってしまうので、Dにはcをつけず、Cにだけcをつける。
初めて分析するときはパッケージインストーラでmultcompパッケージをダウンロードする。この時、「依存パッケージも含める」にチェックを入れておくと、multcompが内部で使っている別のパッケージも同時にインストールされる。
multcompが使っているパッケージ：MASS, mvtnorm, sandwich, survival, TH.data, zoo (MASSとsurvivalはおそらく最初からインストールされている）
> library(multcomp) # 一度実行すればRを終了するまで有効 > res<-aov(value~fac,data) > summary(res) > tuk<-glht(res,linfct=mcp(fac="Tukey")) # mcp(の後のfacはaovでグループ分けに使った因子 > cld(tuk)
とする。牛の乳量のデータでは
> data<-read.table("cow.csv", sep=",", header=T)
> res<-aov(Yield~Food+Cow+Season,data)
> summary(res)
# Foodの水準間でのみ有意差がある
> tuk<-glht(res,linfct=mcp(Food="Tukey"))
> cld(tuk)
A B C D
"bc" "c" "a" "ab"
上の例では大きい順にアルファベットを振ったが、cld関数は何も指定しなければ小さい順に振るのでaとcが逆になっている。cld(tuk, decreasing=T)と降順を指定すれば上の結果と同じになる。