差分

このページの2つのバージョン間の差分を表示します。

--- メモ:分散分析:2要因分散分析対応なし [2015/12/27 17:38] – 作成 Wiki Editor
+++ メモ:分散分析:2要因分散分析対応なし [2015/12/28 12:44] (現在) – Wiki Editor
@@ 行 25: / 行 25: @@
 |  3|  男|  日本|  1|
 |  ...  ||||
+  > 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  男 日本     3
+   2  男 日本     3
+   3  男 日本     1
+   4  男 日本     3
+   5  男 日本     5
+   6  男 韓国     4
+   7  男 韓国     3
+   8  男 韓国     4
+   9  男 韓国     5
+10  男 韓国     7
+11  男 米国     6
+12  男 米国     6
+13  男 米国     6
+14  男 米国     4
+15  男 米国     8
+16  女 日本     6
+17  女 日本     6
+18  女 日本     5
+19  女 日本     4
+20  女 日本     6
+21  女 韓国     2
+22  女 韓国     6
+23  女 韓国     3
+24  女 韓国     6
+25  女 韓国     4
+26  女 米国     3
+27  女 米国     2
+28  女 米国     3
+29  女 米国     2
+30  女 米国     1
+要因以上の分散分析の場合、タイプ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
+要因数が多いと面倒だが、手順は単純なのでうまく関数にまとめておけば、便利そうだ。
+単純主効果について検討できたところで、SPSSでは各水準における多重比較が行なわれるが、この推定値の求め方がよくわからない。いきなり有意確率が表示されているが、おそらくペアごとにt検定を行った結果であると思われる。しかし、通常通りt検定を実施しても同一の結果とならない。ここは、やはり群ごとにサブデータを作って分散分析を行う方法の方が、簡便のようだ。。（しかし気になる）。