トップ > 研究 > R言語で分散分析（と多重比較）

R言語で分散分析（と多重比較）

研究 コンピュータ

例えば、A〜Eの5種類の肥料をそれぞれ別の畑で使ったときのジャガイモの収穫量を表す以下のようなデータがあったとき、肥料による収穫量の明確な違いがあるかどうかを調べようと思います。一元配置の分散分析（対応なし）を行います。（ジャガイモ収穫量の表は、P. Howell "Elementary Statistics"から抜粋）

A B C D E
310 353 366 299 367
284 293 335 264 314
307 306 339 311 377
267 308 312 266 342 
> A <- c(310, 284, 307, 267)
> B <- c(353, 293, 306, 308)
> C <- c(366, 335, 339, 312)
> D <- c(299, 264, 311, 266)
> E <- c(367, 314, 377, 342)
> fertilizer <- factor(c(rep("A",4), rep("B",4), rep("C",4), rep("D",4), rep("E",4)))
> summary(aov(c(A,B,C,D,E)~fertilizer))
            Df  Sum Sq Mean Sq F value   Pr(>F)   
fertilizer   4 12712.0  3178.0   5.406 0.006716 **
Residuals   15  8818.0   587.9                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

肥料の違いによる収穫量にはp<0.01で有意差が出ましたが、このままでは「少なくともどれか1つの肥料が違う」ということしか言えません。そこでTukey's HSDで多重比較をして「どれが違うのか」を調べます。

> TukeyHSD(aov(c(A,B,C,D,E)~fertilizer))
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = c(A, B, C, D, E) ~ fertilizer)

$fertilizer
    diff         lwr         upr     p adj
B-A   23  -29.940872  75.9408721 0.6711716
C-A   46   -6.940872  98.9408721 0.1042010
D-A   -7  -59.940872  45.9408721 0.9935253
E-A   58    5.059128 110.9408721 0.0286488
C-B   23  -29.940872  75.9408721 0.6711716
D-B  -30  -82.940872  22.9408721 0.4355154
E-B   35  -17.940872  87.9408721 0.2941532
D-C  -53 -105.940872  -0.0591279 0.0496798
E-C   12  -40.940872  64.9408721 0.9533486
E-D   65   12.059128 117.9408721 0.0130496

p<0.05レベルで有意差が出たのはE>A、DDの3つです。どうやらEとCが他の肥料よりも優れていそうです。

ここでもし使われた畑が12の別々の場所ではなく4区画だったとしたら、一元配置の分散分析（対応有り）になります。

A B C D E
1 310 353 366 299 367
2 284 293 335 264 314
3 307 306 339 311 377
4 267 308 312 266 342 
> farm<-factor(rep(c(1:4),5))
> summary(aov(c(A,B,C,D,E)~fertilizer+farm))
            Df  Sum Sq Mean Sq F value    Pr(>F)    
fertilizer   4 12712.0  3178.0  15.970 9.466e-05 ***
farm         3  6430.0  2143.3  10.771  0.001013 ** 
Residuals   12  2388.0   199.0                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

となり、肥料と区画両方に有意差が現れます。もし区画による違いを誤差として扱いたければsummary(aov(c(A,B,C,D,E)~fertilizer+Error(farm)))とします。

SPSSならメニューでポン！なのですが、似たような分析を何度も行う必要があるときにはGUIはめんどくさくなってきます（SPSSでも高機能なスクリプト言語が使えるので、それを勉強するのも手ですが）。Rだと基本的にキーボード入力なので試行錯誤に向いているような気がしますし、柔軟な分析が無料で行えるのは素晴らしいです。なにより、必要とあればソースコードが読めるというのが最高です。