文書の過去の版を表示しています。
分散分析のメモ
Rで各種分散分析を実行する際のメモ。
竹原卓真（2007）『SPSSのススメ1 2要因の分散分析をすべてカバー』北大路書房を参照して、Rで同じ分析結果を再現しながら勉強。
備忘録
ネストした分散分析について
ネスト構造にあるデータの分散分析では、親要素のみで分散分析をした場合と、子要素もモデルに含めて分析した場合では、親要素の平方和は変わらない。子要素を投入すると、誤差の平方和が子要素で説明できる分だけ減少し、その分が子要素の平方和となる。
ネスト構造にあるデータでは、共通の親要素の範囲内に限定して平方和を計算するようだ。Rではネストしたデータの分散分析を以下のように行う。
summary(aov(SCORE ~ A / b, dat=DATA))
この場合、Aが親要素で、Aの中にbがネストされている。この式ではRは自動的に
summary(aov(SCORE ~ A + A:b, dat=DATA))
と展開して分散分析をを実施する。つまり親要素の主効果と、親要素と子要素の交互作用の2つの変数を要因とした分散分析になる。子要素は親要素ごとに分析されるということを考えれば、親要素と子要素の交互作用がネストされたデータの子要素を示す要因であることが理解できそうである。
ところで、ネストされたデータは
親 子 SCORE A a 10 A b 20 A c 30 B d 20 B e 25 B f 30 C g 15 C h 11 C i 23
というデータ構造になる。このようにコーディングされたデータでは、以下の分析結果は同じになる。
summary(aov(SCORE ~ 親 + 子, dat=DATA)) #親・子を独立した変数として分析 summary(aov(SCORE ~ 親 / 子, dat=DATA)) #子が親にネストされたデータとして分析 summary(aov(SCORE ~ 親 + 親:子, dat=DATA)) #ネストされたデータの別表現
このデータでは子要素は親要素に関係なく、すべてのサンプルで異なったコードが割り当てられている。そのため平方和の計算では、親子が独立した変数として投入されても、ネストされていても同じになる。前者では親要素が同時にモデルに存在するので親要素が統制され、親要素ごとに平方和が計算される。ネストされていても同様の計算がされるので、考え方は異なってもやっている計算結果は結局同じ。
ここでデータ構造を少し変えてみよう。
親 子 SCORE A a 10 A b 20 A c 30 B a 20 B b 25 B c 30 C a 15 C b 11 C c 23
先のデータから少しコーディングを変更して、子要素のコードが親要素ごとに再利用するようにした。この時、子要素は親要素にネストされているものなので、子要素のコードは単独では意味がなく親要素と組み合わせて初めて一意なデータを示す。そうすると、親子を独立した変数としてモデルに投入する場合とは結果が異なってくる。
summary(aov(SCORE ~ 親 + 子, dat=DATA)) #親・子を独立した変数として分析
この場合、親子が独立した変数として投入されているので、各親要素の中にある子要素（例えばaやb）は子要素変数の同じ水準を示すことになる。しかし、ネストされているデータでは親Aの子aと親Bの子aは無関係である。親要素ごとにデータが収集された際に、親ごとに連番を与えている場合などで、このようなことが生じる。そこで、きちんとネスト関係を指示してやる。
summary(aov(SCORE ~ 親 / 子, dat=DATA)) #子が親にネストされたデータとして分析
こうすることで、子は親にネストされているものとして平方和が計算される。具体的には、子要素の平方和は親要素との交互作用で計算される。そのため、親要素をまたがって子要素のコードが重なっていても、親要素が異なれば計算上は無関係となるので、期待通りの計算結果が得られる。
単純主効果が有意な場合の多重比較について
2要因の分散分析において交互作用が有意になった場合、単純主効果の検定を行う。
ここでは竹内理・水本篤『外国語教育研究ハンドブック 【改訂版】 研究手法のより良い理解のために 』2014年8月の第7章を例に、筆者が提供してくれるデータ（繰り返しなしの二元配置分散分析）を用いて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