まずは、単純な群間比較。

データの準備

このままSPSSのように入力するとうまくいかないので、クラスを分ける変数を作成し以下のようなデータを作成する。IDは特に必要ないが、一応。

> id <- 1:20;
> class <- c("A","A","A","A","A","B","B","B","B","B","C","C","C","C","C","D","D","D","D","D");
> score <- c(8,7,9,12,10,5,6,6,7,3,12,10,9,14,14,11,10,15,12,13);
> df <- data.frame(ID=id, CLASS=class, SCORE=score);

各列のベクトルを作成し、データフレームにまとめる。CSVファイルから読む込む方が面倒はないが、ちょっと練習。

lm()関数で推定。CLASS変数はfactorなので、自動的にダミー変数に展開される。

> res.anova <- lm(SCORE ~ CLASS, df);
> summary(res.anova);

Call:
lm(formula = SCORE ~ CLASS, data = df)

Residuals:
   Min     1Q Median     3Q    Max
 -2.80  -1.35   0.00   1.00   2.80

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   9.2000     0.8631  10.659 1.12e-08 ***
CLASSB       -3.8000     1.2207  -3.113  0.00669 **
CLASSC        2.6000     1.2207   2.130  0.04904 *
CLASSD        3.0000     1.2207   2.458  0.02577 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.93 on 16 degrees of freedom
Multiple R-squared:  0.7115,    Adjusted R-squared:  0.6573
F-statistic: 13.15 on 3 and 16 DF,  p-value: 0.0001383

回帰係数表を、記述統計出力に代用。Estimateに平均値が出力される。各郡の平均は、（Intercept）に各郡の推定値を加える。標準誤差はやや数値が異なる。一般線形モデルの場合と従来型の分散分析とは計算方法が異なるので、ここではあまり気にしない。

次に分散分析表

> anova(res.anova);
Analysis of Variance Table

Response: SCORE
          Df Sum Sq Mean Sq F value    Pr(>F)
CLASS      3 146.95  48.983   13.15 0.0001383 ***
Residuals 16  59.60   3.725

平方和の合計は出力されないが、足し算すれば良いだけ。単純な分散分析は、簡単。

続いて、Tukey HSDによる多重比較。

> res.aov <- aov(SCORE ~ CLASS, df);
> TukeyHSD(res.aov);
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = SCORE ~ CLASS, data = df)

$CLASS
    diff        lwr        upr     p adj
B-A -3.8 -7.2923198 -0.3076802 0.0306239
C-A  2.6 -0.8923198  6.0923198 0.1857885
D-A  3.0 -0.4923198  6.4923198 0.1057960
C-B  6.4  2.9076802  9.8923198 0.0004223
D-B  6.8  3.3076802 10.2923198 0.0002236
D-C  0.4 -3.0923198  3.8923198 0.9873918

lm()関数が返すオブジェクトではダメなので、aov()関数の結果をTukeyHSD関数に渡す。