要因が2つで、1つ1つの観測値がすべて別人の場合の分散分析。教科書の例では、国別の男女で3カ国の男女についてそれぞれ観測値がある。各国の男女に対応はないし、国の間でも対応がない。つまり、観測値はすべて別人のもの。

まずはデータ。見やすい表で整理した場合。

分析用データは、ID、性別、国、スコア変数を作成して入力していく。

> ID <- rep(1:30);
> SEX <- c(rep("男", 15), rep("女", 15));
> CN <- rep(c(rep("日本", 5), rep("韓国", 5), rep("米国", 5)), 2);
> SCORE <- c(3,3,1,3,5,4,3,4,5,7,6,6,6,4,8,6,6,5,4,6,2,6,3,6,4,3,2,3,2,1);
> df <- data.frame(ID=ID, SEX=SEX, CN=CN, SCORE=SCORE);
> df
   ID SEX   CN SCORE
1   1  男 日本     3
2   2  男 日本     3
3   3  男 日本     1
4   4  男 日本     3
5   5  男 日本     5
6   6  男 韓国     4
7   7  男 韓国     3
8   8  男 韓国     4
9   9  男 韓国     5
10 10  男 韓国     7
11 11  男 米国     6
12 12  男 米国     6
13 13  男 米国     6
14 14  男 米国     4
15 15  男 米国     8
16 16  女 日本     6
17 17  女 日本     6
18 18  女 日本     5
19 19  女 日本     4
20 20  女 日本     6
21 21  女 韓国     2
22 22  女 韓国     6
23 23  女 韓国     3
24 24  女 韓国     6
25 25  女 韓国     4
26 26  女 米国     3
27 27  女 米国     2
28 28  女 米国     3
29 29  女 米国     2
30 30  女 米国     1

2要因以上の分散分析の場合、タイプ1平方和では問題がでる。タイプ1平方和だと、要因を検討する順序によって、平方和が変化する。それに対応するため、SPSSではタイプ3平方和が使われる。Rでタイプ3平方和を計算するにはcarパッケージのAnova()関数が使えるが、その際、対比の方法を適切に設定してやる必要がある。Rでは、カテゴリ変数は適切にダミー変数に変換されるが、その変換方式によって平方和の計算が変わってくる。タイプ3平方和を適切に計算するには、contrastsをcontr.sumに変更しておく必要がある（1行目）。contrastsの2番目に指定した要素は、順序カテゴリカル変数の場合の指定。今回は無関係。

> res.lm <- lm(SCORE ~ SEX*CN, df);
> Anova(res.lm, type=3);
Anova Table (Type III tests)

Response: SCORE
            Sum Sq Df  F value    Pr(>F)
(Intercept) 537.63  1 293.2545 5.871e-15 ***
SEX           2.70  1   1.4727 0.2367243
CN            0.47  2   0.1273 0.8810838
SEX:CN       48.20  2  13.1455 0.0001395 ***
Residuals    44.00 24

まず、lm()関数で推定し、その結果を分散分析する。lm()関数では、性別と国に加えて、両者の交互作用も含めて検討する。

結果は、主効果は有意ではないが、交互作用が有意になっている。なので、単純主効果の検討に入る。単純主効果の検討にはいろいろ方式があるようだ。SPSSでは、第1要因ごとに群分けしたデータで第2要因について分散分析を行うが、誤差項のとりかたは群ごとの分散分析の誤差平方和および自由度の合計として検討している。なので、Rで同様の分析を行うにはやや複雑になる。単純に、第1要因ごとに群分けして第2要因について分散分析を行うだけで済ましている場合も見かける。ここでは、SPSSと同じ分析を行うことが主眼なので、面倒でもやってみよう。

追記：せっかく計算したが、実は上記の分散分析表の誤差平方和（Residuals）の値であった。

#全体データから国ごとのサブセットを作成
> dfsub1 <- subset(df, CN=="日本");
> dfsub2 <- subset(df, CN=="韓国");
> dfsub3 <- subset(df, CN=="米国");

#サブセット（国）ごとに分散分析
> res.anova1 <- anova(lm(SCORE ~ SEX, dfsub1));
> res.anova2 <- anova(lm(SCORE ~ SEX, dfsub2));
> res.anova3 <- anova(lm(SCORE ~ SEX, dfsub3));

#各分散分析の誤差平方和を足し合わせる
> errSS <- res.anova1$"Sum Sq"[2]+res.anova2$"Sum Sq"[2]+res.anova3$"Sum Sq"[2];
#同様に、各分散分析の自由度を足し合わせる
> errDF <- res.anova1$Df[2]+res.anova2$Df[2]+res.anova3$Df[2];

#F値を求める
#前半は平方和を自由度で除して平均平方を求める
#後半は先ほど求めた全体の誤差方法和を誤差自由度で除して誤差の平均平方を求める
#前者を後者で除してF値を算出
> F1 <- (res.anova1$"Sum Sq"[1] / res.anova1$Df[1]) / (errSS / errDF);
> F2 <- (res.anova2$"Sum Sq"[1] / res.anova2$Df[1]) / (errSS / errDF);
> F3 <- (res.anova3$"Sum Sq"[1] / res.anova3$Df[1]) / (errSS / errDF);

#F値から有意確率を求めて検定する
> round(1-pf(F1, res.anova1$Df[1], errDF), 3);
[1] 0.01
> round(1-pf(F2, res.anova2$Df[1], errDF), 3);
[1] 0.645
> round(1-pf(F3, res.anova3$Df[1], errDF), 3);
[1] 0

要因数が多いと面倒だが、手順は単純なのでうまく関数にまとめておけば、便利そうだ。