差分

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

--- メモ:分散分析 [2017/09/07 09:56] – Wiki Editor
+++ メモ:分散分析 [2017/09/08 19:44] (現在) – Wiki Editor
@@ 行 11: / 行 11: @@
 ===== 備忘録 =====
-ネストした分散分析について
+==== ネストした分散分析について ====
 ネスト構造にあるデータの分散分析では、親要素のみで分散分析をした場合と、子要素もモデルに含めて分析した場合では、親要素の平方和は変わらない。子要素を投入すると、誤差の平方和が子要素で説明できる分だけ減少し、その分が子要素の平方和となる。
@@ 行 67: / 行 68: @@
 こうすることで、子は親にネストされているものとして平方和が計算される。具体的には、子要素の平方和は親要素との交互作用で計算される。そのため、親要素をまたがって子要素のコードが重なっていても、親要素が異なれば計算上は無関係となるので、期待通りの計算結果が得られる。
+==== 単純主効果が有意な場合の多重比較について ====
+要因の分散分析において交互作用が有意になった場合、単純主効果の検定を行う。
+ここでは竹内理・水本篤『外国語教育研究ハンドブック 【改訂版】 研究手法のより良い理解のために 』2014年8月の第7章を例に、筆者が提供してくれる[[http://mizumot.com/handbook/?page_id=246|データ]]（繰り返しなしの二元配置分散分析）を用いてRで同様の分析する方法をメモしておく。ポイントは単純主効果が有意な場合の多重比較の方法であるが、ついでなので7章の分析方法を全て再現しておく。
+SPSSでは一般線形モデルで分析を行うようである。分散分析表で平方和のTypeがIIIとなっているので、RでもType IIIの平方和を計算するためにcarパッケージを読込、Anova関数を利用する。
+  > library(car) # パッケージの読込
+  # Type III平方和を計算する際には、factorのダミーコードの方法を「contr.sum」に変更しないと計算結果がおかしくなる
+  > options(contrasts=c("contr.sum","contr.sum"))
+  # 筆者のサイトからダウンロードしてデータを読み込む
+  > DAT <- read.csv("データ")
+  # 数値でコーディングされているのでfactorとして認識させる
+  > DAT$教室 <- factor(DAT$教室);
+  > DAT$クラスサイズ <- factor(DAT$クラスサイズ);
+  # carパッケージのAnova関数を利用して分散分析
+  > Anova(lm(点数 ~ 教室 * クラスサイズ, data=DAT), type=3);
+  Anova Table (Type III tests)
+  Response: 点数
+                    Sum Sq  Df   F value    Pr(>F)
+  (Intercept)       664301   1 6575.8307 < 2.2e-16 ***
+  教室                  26   1    0.2618    0.6095
+  クラスサイズ       11238   2   55.6219 < 2.2e-16 ***
+  教室:クラスサイズ   4057   2   20.0823 1.421e-08 ***
+  Residuals          17578 174
+  # 偏イータ2乗の計算
+  # 対象の変数の平方和 / （対象の変数の平方和 + 誤差平方和）
+  # Anovaの結果を変数に投入
+  > res.anova <- Anova(lm(点数 ~ 教室 * クラスサイズ, data=DAT), type=3);
+  # res.anovaから平方和を取り出して計算
+  # 教室の偏イータ2乗
+  > res.anova$"Sum Sq"[2] / (res.anova$"Sum Sq"[2] + res.anova$"Sum Sq"[5])
+  [1] 0.001502481
+  # クラスサイズの偏イータ2乗
+  > res.anova$"Sum Sq"[3] / (res.anova$"Sum Sq"[3] + res.anova$"Sum Sq"[5])
+  [1] 0.3899955
+  # 交互作用の偏イータ2乗
+  > res.anova$"Sum Sq"[4] / (res.anova$"Sum Sq"[4] + res.anova$"Sum Sq"[5])
+  [1] 0.187541
+  # 単純主効果の検定
+  # データを一方の要因の水準毎に分割して他方の要因で1要因分散分析を実施
+  # ただし、検定の際の誤差項は全体の分散分析の値を用いる必要がある
+  # クラスサイズが10人のデータに限定して、教室要因で分散分析
+  # ここでのF値や検定結果は使わない
+  > Anova(lm(点数 ~ 教室, data=DAT[DAT$クラスサイズ == "1",]), type=3);
+  Anova Table (Type III tests)
+  Response: 点数
+              Sum Sq Df  F value    Pr(>F)
+  (Intercept) 283731  1 3725.137 < 2.2e-16 ***
+  教室          1421  1   18.657 6.199e-05 ***
+  Residuals     4418 58
+  # 結果を仮に「res.クラスサイズ1」に代入しておく
+  > res.クラスサイズ1 <- Anova(lm(点数 ~ 教室, data=DAT[DAT$クラスサイズ == "1",]), type=3);
+  # 全体の分散分析結果「res.anova」と水準毎の結果「res.クラスサイズ1」からF値を計算
+  # F値は要因の平均平方 / 誤差の平均平方で計算する。
+  # Anova関数では平均平方が出力されないので、平均平方を自由度で除して計算する
+  # クラスサイズの水準毎に分割したデータを用いた分散分析の「教室」要因の平均平方を、
+  #  最初に行った全体の分散分析の誤差平均平方で除した値が、検定に用いるF値となる
+  > (res.クラスサイズ1$"Sum Sq"[2]/res.クラスサイズ1$Df[2]) / (res.anova$"Sum Sq"[5]/res.anova$Df[5])
+    [1] 14.06695
+  # 変数に保持しておく
+  F.クラスサイズ1 <- (res.クラスサイズ1$"Sum Sq"[2]/res.クラスサイズ1$Df[2]) / (res.anova$"Sum Sq"[5]/res.anova$Df[5])
+  # 求めたF値からp値を求める。
+  # pf関数に、F値と自由度を与える。lower.tailは上限からの累積確率を求めるオプション
+  > pf(F.クラスサイズ1, res.クラスサイズ1$Df[2], res.anova$Df[5], lower.tail=FALSE)
+  [1] 0.0002400081
+  # p値の値が小さいので書籍の数値と一致しているかわからないので、クラスサイズが20人の場合も計算して確認
+  > res.クラスサイズ2 <- Anova(lm(点数 ~ 教室, data=DAT[DAT$クラスサイズ == "2",]), type=3);
+  > res.クラスサイズ2
+  Anova Table (Type III tests)
+  Response: 点数
+              Sum Sq Df   F value Pr(>F)
+  (Intercept) 241808  1 2913.2487 <2e-16 ***
+  教室            23  1    0.2749 0.6021
+  Residuals     4814 58
+  > F.クラスサイズ2 <- (res.クラスサイズ2$"Sum Sq"[2]/res.クラスサイズ2$Df[2]) / (res.anova$"Sum Sq"[5]/res.anova$Df[5])
+  > F.クラスサイズ2
+  [1] 0.2258592
+  > pf(F.クラスサイズ2, res.クラスサイズ1$Df[2], res.anova$Df[5], lower.tail=FALSE)
+  [1] 0.6352074
+単純主効果が有意な場合の多重比較。ものによっては、一方の水準毎に分割したデータで分散分析を行ってTukeyHSDなどで多重比較するだけでよいと説明しているところもあるが、SPSSと同様の結果を得るにはやや面倒な手続きが必要。
+  # 要因1「教室」と要因2「クラスサイズ」を掛け合わせて新しい変数を作成する
+  # つまり、教室の各水準とクラスサイズの各水準を掛け合わせてデータを分割する
+  # 教室の各水準とクラスサイズの各水準の値を組み合わせて変数を作成
+  DAT$CROSS <- paste("教", DAT$教室, "ク", DAT$クラスサイズ, sep="")
+  # pasteで作成したCROSSという変数の型は文字列なのでfactorに変換する
+  DAT$CROSS <- factor(DAT$CROSS)
+  # このように2つの変数を掛け合わせて作成した新しい変数でデータを分割して、それぞれの多重比較を行う
+  # ここでは教室が2水準で、それぞれクラスサイズが3水準あるので、カテゴリは6つになる。
+  # pairwise.t.test関数を利用してDAT$CROSSについて1対毎にT検定を行う
+  # ただし、ここでp値の調整オプションをつけてはいけない。Bonferroniの方法で調整を行う際の帰無仮説の数は
+  #  6つの全てのカテゴリではなく、「教室」要因毎であれば「クラスサイズ」の3水準の多重比較なので、
+  #  3(3-1)/2=3回の検定の繰り返しを調整しなければならない。pairewise.t.test関数でp.adjust.method="bonferroni"を
+  #  指定してしまうと、6(6-1)/2=15で調整してしまう。
+    > pairwise.t.test(DAT$点数, DAT$CROSS, p.adjust.method="none")
+          Pairwise comparisons using t tests with pooled SD
+  data:  DAT$点数 and DAT$CROSS
+         教1ク1    教1ク2    教1ク3    教2ク1    教2ク2
+  教1ク2 0.00032   -         -         -         -
+  教1ク3 < 2e-16   1.8e-13   -         -         -
+  教2ク1 0.00024   0.93866   2.8e-13   -         -
+  教2ク2 5.2e-05   0.63521   2.9e-12   0.69099   -
+  教2ク3 6.2e-10   0.00452   8.3e-07   0.00569   0.01736
+  P value adjustment method: none
+  # ここで、単純に先ほど作成した「CROSS」の水準毎にt検定を行ってもよいように思うが、
+  # t値を求める際の標準偏差としてプールされた標準偏差を利用する必要があるため、1対毎のt検定ではだめ
+  # 交互作用がある場合の多重比較では全ての2要因をクロスして作成した全てのグループのプールされた標準偏差を
+  # 利用するようだ。SPSSの出力を再現するには、この手順を踏む必要がある。
+  # 結果を変数に保持したうえでBonferrniの方法でp値を調整する
+  > 多重比較p値 <- pairwise.t.test(DAT$点数, DAT$CROSS, p.adjust.method="none")
+  # 例えば、教室が「CALL教室」の場合のクラスサイズの各水準の多重比較ならば
+  #  「教2ク1」「教2ク2」「教2ク3」の検定結果を取り出して、p値を調整する
+  # 必要なp値を取り出してベクトルに格納
+  > 対象p値 <- c(多重比較p値$p.value["教2ク2","教2ク1"], 多重比較p値$p.value["教2ク3","教2ク1"], 多重比較p値$p.value["教2
+ク3","教2ク2"])
+  # もしくは
+  > 対象p値 <- c(多重比較p値$p.value[4:5,4], 多重比較p値$p.value[5,5])
+  # Bonferroniの方法でp値を調整
+  > p.adjust(対象p値, method="bonferroni")
+  [1] 1.00000000 0.01705889 0.05208608
+  # ついでのholmの方法も。Rの分散分析でよく使われるAnova君ではこの値が出力されるようだ。
+  > p.adjust(対象p値, method="holm")
+  [1] 0.69098604 0.01705889 0.03472406
+　以上が、外国語教育研究ハンドブックの対応のない2元配置の分散分析のRでの再現方法となる。
+　素直に「ANOVA君」を使えばよいのだが、計算方法（というか計算はRが勝手にやるが）の理解のためにRで手順を追って調べてみた。