差分

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

--- メモ:分散分析:1要因分散分析対応あり [2015/12/25 22:44] – 作成 Wiki Editor
+++ メモ:分散分析:1要因分散分析対応あり [2015/12/27 17:59] (現在) – Wiki Editor
@@ 行 113: / 行 113: @@
 Mauchlyの球面性検定を単体で行う関数もあるが、carパッケージのAnova()関数の方がεも同時に計算してくれたり、分散分析も一気にやってくれるので便利。以下でTEST.factorに水準数を設定したベクトルを作ったりしているが、これの意味はまだ調べていない。とりあえず、適当に3要因のベクトルを作成しておく。
+  > library(car);
   > TEST.factor <- factor(c("A","B","C"));
   > res.anova <- Anova(res.mlm, idata=data.frame(TEST.factor), idesign=~TEST.factor, type=3);
@@ 行 147: / 行 148: @@
 出力の下半分で、Mauchlyの球面性検定が行われている。MauchlyのWは「0.14767」、p値は「0.056746」となり、SPSSと一致。イプシロンのGreenhouse-Geisserは「GG eps」の値で「0.53986」、Huynh-Feldtは「HF eps」の値で「0.5818973」。
-Mauchlyの球面性検定が有意なら、球面性が__仮定されない__ので、自由度の補正が必要なようだ。ここでは、分散分析表のSSの「419」がSPSSでいう「タイプIII平方和」で要因の平方和、その右側のError SSが誤差の平方和で、自由度がそれぞれ2,8となっている。球面性が仮定されなければ、自由度に「GG eps」および「HF eps」をかけて自由度を調整し検定をすればよい。オブジェクトから直接値を抜き出す方法がまだわからないので、いまのところ手計算。とりあえずSPSSと同じ結果が得られたので今のところはOKとする。
+Mauchlyの球面性検定が有意なら、球面性が__仮定されない__ので、自由度の補正が必要なようだ。ここでは、分散分析表のSSの「419」がSPSSでいう「タイプIII平方和」で要因の平方和、その右側のError SSが誤差の平方和で、自由度がそれぞれ2,8となっている。球面性が仮定されなければ、自由度に「GG eps」および「HF eps」をかけて自由度を調整し検定をすればよい。<del>オブジェクトから直接値を抜き出す方法がまだわからないので、いまのところ手計算。とりあえずSPSSと同じ結果が得られたので今のところはOKとする</del>。あまり勉強もしないままソースまで追いかけて、半日かかって単純なことに気がつく。値を抜き出す方法が分かった。
+summary関数は結果を表示するものだと思っていたが、これをオブジェクトに入れておくと後で個別の値にアクセスできる。また、GG eps、HF epsで自由度を調整した有意確率についても、実はGG、HFの値の後にすでに計算済みの値が表示されていた。
+  > summary.anova <- summary(res.anova);
+  #要因の平方和
+  > summary.anova$univariate.tests[2,1]
+  #要因の自由度
+  > summary.anova$univariate.tests[2,2]
+  #誤差の平方和
+  > summary.anova$univariate.tests[2,3]
+  #誤差の自由度
+  > summary.anova$univariate.tests[2,4]
+  #Greenhouse-Geisserのイプシロン
+  > summary.anova$pval.adjustments[1]
+  #Greenhouse-Geisserのイプシロンによって自由度を調整した後の有意確率（計算しなくてもここを見れば良かった）
+  > summary.anova$pval.adjustments[2]
+  #Huynh-Feldtのイプシロン
+  > summary.anova$pval.adjustments[3]
+  #Huynh-Feldtのイプシロンによって自由度を調整した後の有意確率（計算しなくてもここを見れば良かった）
+  > summary.anova$pval.adjustments[4]
+遠回りしたが、Rについて幾つかわかったことがあるし、球面性の仮定が成り立たないときの有意確率の計算方法もわかったので、無駄ではなかったということにしよう。
+最後に、多重比較。ここではボンフェローニの方法で検討されているので、以下のようになる。
+  > pairwise.t.test(df$SCORE, df$TEST, paired=T, p.adjust.method="bonferroni");
+          Pairwise comparisons using paired t tests
+  data:  df$SCORE and df$TEST
+    A       B
+  B 5.8e-05 -
+  C 0.3746  0.0055
+  P value adjustment method: bonferroni
+最初の引数は平均を求めるデータ、次がグループ、pairedは対応のある場合、p.adjust.methodに有意確率の調整方法としてボンフェローニの方法。ただし、__対応のある場合の多重比較では、グループごとに個人のデータが同じ順で並んでいる必要__がある。したがって、データフレームをグループ別にソートしたうえで、同グループ内で個人のID順に並べておかなければ、正しい結果は得られない。教科書的には、サンプルデータがそのように並んでいる場合が多いので、対応のあるデータの場合は特に気をつける必要があると思う。