Rで、一元配置分散分析を step by step で計算してみた。
lm() と Anova() を使えばあっという間だが、具体的な一つ一つの計算を自分で組み立ててみるとどうか？
教科書の例題に沿って確認した。
- R で一元配置分散分析を実行するための例題データ
- R で一元配置分散分析 step by step
- R で一元配置分散分析を lm() を使って実行する
- R で一元配置分散分析 EZR で行うとどうか？
- R で一元配置分散分析する際 群のサンプルサイズがアンバランスの場合は？
- まとめ
- 参考書籍
R で一元配置分散分析を実行するための例題データ
サンプルデータ（記事最下段 参考書籍 例題8.1 ）は以下の通り。
A1 <- c(64,72,68,77) A2 <- c(82,78,77,85) A3 <- c(55,64,66,49) # 以下のように変形する。 X.all <- c(A1,A2,A3) A.all <- c(rep("A1",length(A1)),rep("A2",length(A2)),rep("A3",length(A3)))
元のデータセット。
> data.frame(A1,A2,A3) A1 A2 A3 1 64 82 55 2 72 78 64 3 68 77 66 4 77 85 49
以下が解析用データセット。
> data.frame(cbind(X.all, A.all)) X.all A.all 1 64 A1 2 72 A1 3 68 A1 4 77 A1 5 82 A2 6 78 A2 7 77 A2 8 85 A2 9 55 A3 10 64 A3 11 66 A3 12 49 A3
R で一元配置分散分析 step by step
一元配置分散分析は、分散分析とは言うが、要因のグループ（A1, A2, A3）ごとの平均値が、全部同じという帰無仮説を検定する、三群以上の平均値の差の検定だ。
準備として、全体平均（Xbar）、グループの平均（Xjbar）、グループのサンプルサイズ（r）、グループの数（a）、全体のサンプルサイズ（N）を計算する。
次に、三つの平方和を計算する。
全体平方和（SS）、グループ間の平方和（SSA）、グループ内の平方和（SSE）。
SS は、おのおののデータと全体平均の差の二乗和。
SSAは、グループの平均と全体平均の差の二乗和をグループのサンプルサイズ倍したもの。
SSEは、SS から SSA を引いたもの。
さらに、自由度の計算。
全体（nu）は N-1、グループ間平方和の自由度（nuA）は a-1、グループ内平方和の自由度は N-a と計算する。
そして、平方和を自由度で割って、平均平方和を算出する。
グループ間（VA）は SSA/nuA、グループ内（VE）は SSE/nuE で計算できる。
最後に、VAとVEの比を取る。
これが検定統計量 F値になる。
このF値とnuA, nuEを使って、有意確率を計算する。
one.way.anova <- function(X,A){ Xbar <- mean(X) Xjbar <- tapply(X, A, mean) r <- matrix(table(A))[1] a <- length(table(A)) N <- length(X) SS <- sum((X-Xbar)^2) SSA <- r*sum((Xjbar-Xbar)^2) SSE <- SS-SSA nu <- N-1 nuA <- a-1 nuE <- N-a VA <- SSA/nuA VE <- SSE/nuE F <- VA/VE pF <- pf(F, nuA, nuE, lower.tail=FALSE) list(c(SS_A=SSA, nu_A=nuA, V_A=VA, SS_E=SSE, nu_E=nuE, V_E=VE, F=F, p_F=pF, SS=SS, nu=nu)) } one.way.anova1 <- unlist(one.way.anova(X=X.all,A=A.all)) one.way.anova1
結果はこちら。
> one.way.anova1 SS_A nu_A V_A SS_E nu_E 9.695000e+02 2.000000e+00 4.847500e+02 3.227500e+02 9.000000e+00 V_E F p_F SS nu 3.586111e+01 1.351743e+01 1.944638e-03 1.292250e+03 1.100000e+01
つまり、A1, A2, A3の三グループの平均は等しくないと言える。
ちなみに三グループの平均値は、以下の通り。
A3だけ低く見える。
> tapply(X.all, A.all, mean) A1 A2 A3 70.25 80.50 58.50
R で一元配置分散分析を lm() を使って実行する
Anova() を使うため、car パッケージをインストールしていなければ、先にインストールする。
install.packages("car")
インストールは一回のみで OK
X.all を従属変数側に、A.allを独立変数側に書いて チルダ （~）で結ぶ。
分散分析表を Anova() 出力させる。
library(car) # 使うときはこの library で呼び出す lm1 <- lm(X.all ~ A.all) Anova(lm1)
結果はこちら。
> Anova(lm1) Anova Table (Type II tests) Response: X.all Sum Sq Df F value Pr(>F) A.all 969.50 2 13.517 0.001945 ** Residuals 322.75 9 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
A.all の Sum Sq（Sum of Square＝平方和）が、SSAに一致している。
また、Residuals の Sum Sq が SSE に一致しているのがわかる。
F値もp値も同一だ。
R で一元配置分散分析 EZR で行うとどうか？
データセットに書き出す
上記の例題データをCSVファイルに書き出す。
A1 <- c(64,72,68,77) A2 <- c(82,78,77,85) A3 <- c(55,64,66,49) X.all <- c(A1,A2,A3) A.all <- c(rep("A1",length(A1)),rep("A2",length(A2)),rep("A3",length(A3))) tab <- data.frame(cbind(X.all, A.all)) write.csv(tab, "data/anova-data1.csv", row.names=FALSE)
もしくはこれと同じ形にExcelに入力してCSVファイル形式で保存する。
X.all A.all 1 64 A1 2 72 A1 3 68 A1 4 77 A1 5 82 A2 6 78 A2 7 77 A2 8 85 A2 9 55 A3 10 64 A3 11 66 A3 12 49 A3
「3群以上の間の平均値の比較」を用いて解析する
anova-data1.csvをEZRに読み込んだ後、「統計解析」→「連続変数の解析」→「3群以上の間の平均値の比較（一元配置分散分析 one-way ANOVA）」を選択する。
X.allを目的変数にして、A.allを比較する群にして、計算すると以下の結果が出力される。
lm() & Anova() の結果と同じである。
> summary(AnovaModel.1) Df Sum Sq Mean Sq F value Pr(>F) factor(A.all) 2 969.5 484.7 13.52 0.00194 ** Residuals 9 322.7 35.9 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
EZRで分散分析【無料統計ソフトEZRで簡単統計】【動画】
R で一元配置分散分析する際 群のサンプルサイズがアンバランスの場合は？
グループのサンプルサイズがそろっていない場合はどうしたらいいのか？
上記はサンプルサイズがそろっていることが前提だった。
サンプルサイズがそろっていないアンバランスなデータの場合、スクリプトをどう変えるのか？
例題のデータ
同じ教科書の例題8.3を使う。
データは以下の通り。
A1からA4の4つのグループで、サンプルサイズは６，４，６，５とそろっていない。
A1 <- c(205,206,164,190,194,203) A2 <- c(201,221,197,185) A3 <- c(248,265,197,220,212,281) A4 <- c(202,276,237,254,230) X.all <- c(A1,A2,A3,A4) A.all <- c(rep("A1",length(A1)),rep("A2",length(A2)),rep("A3",length(A3)),rep("A4",length(A4))) data.frame(X.all, A.all)
この形にする。
> data.frame(X.all, A.all) X.all A.all 1 205 A1 2 206 A1 3 164 A1 4 190 A1 5 194 A1 6 203 A1 7 201 A2 8 221 A2 9 197 A2 10 185 A2 11 248 A3 12 265 A3 13 197 A3 14 220 A3 15 212 A3 16 281 A3 17 202 A4 18 276 A4 19 237 A4 20 254 A4 21 230 A4
グループのサンプルサイズがそろっていない場合の計算
スクリプトは要因の平方和（SSA）を計算するところを変える。
SSA <- r*sum((Xjbar-Xbar)^2) を SSA <- sum(r*matrix((Xjbar-Xbar)^2)) に変える。
もともとは、グループの平均と全体の平均の差の二乗和を計算してから、グループのサンプルサイズ倍していた。
それを、それぞれのグループの平均と全体平均の差の二乗に、それぞれのグループのサンプルサイズを掛けてから、合計する。
グループのサンプルサイズで重みをつける感じ。
one.way.anova <- function(X,A){ Xbar <- mean(X) Xjbar <- tapply(X, A, mean) r <- matrix(table(A)) a <- length(table(A)) N <- length(X) SS <- sum((X-Xbar)^2) SSA <- sum(r*matrix((Xjbar-Xbar)^2)) SSE <- SS-SSA nu <- N-1 nuA <- a-1 nuE <- N-a VA <- SSA/nuA VE <- SSE/nuE F <- VA/VE pF <- pf(F, nuA, nuE, lower.tail=FALSE) list(c(SS_A=SSA, nu_A=nuA, V_A=VA, SS_E=SSE, nu_E=nuE, V_E=VE, F=F, p_F=pF, SS=SS, nu=nu)) } one.way.anova1 <- unlist(one.way.anova(X=X.all,A=A.all)) one.way.anova1 tapply(X.all, A.all, mean)
結果は以下の通り、F値は5.09で、p値は0.0107で統計学的有意に異なる。
グループごとの平均値を算出するとA1とA2が小さめ、A3とA4が大きめだった。
特にA1は小さい。
> one.way.anova1 SS_A nu_A V_A SS_E nu_E 9.284271e+03 3.000000e+00 3.094757e+03 1.033297e+04 1.700000e+01 V_E F p_F SS nu 6.078216e+02 5.091555e+00 1.072138e-02 1.961724e+04 2.000000e+01 > > tapply(X.all, A.all, mean) A1 A2 A3 A4 193.6667 201.0000 237.1667 239.8000
lm() を使ってみる
lm() を使えば話は早い。
lm() はデータのアンバランスなど気にしない。
まったく同じ書き方で結果が出力される。
lm1 <- lm(X.all ~ A.all) Anova(lm1)
結果はこちら。
> Anova(lm1) Anova Table (Type II tests) Response: X.all Sum Sq Df F value Pr(>F) A.all 9284.3 3 5.0916 0.01072 * Residuals 10333.0 17 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’
A.allのSum SqがSSAと、ResidualsのSum SqがSSEと、そしてF値、p値が、それぞれ一歩一歩計算したものと同じことがわかる。
まとめ
統計ソフトRで、lm() を使えば一瞬で出力できる一元配置分散分析を、一つ一つ順を追ってスクリプトに書いて計算してみた。
こうやると簡単な統計手法も理解が深まる。
今回も思ったのだが、グループ間の平方和で、グループのサンプルサイズ倍するところが、何となくしっくりこない。
こういう計算なのだと頭ではわかっているが、何となく変な感じがする。
グループ間の平均と全体平均の平方和だけのほうがグループ間のばらつきを表しているような気がしてしまうのだ。
しかし、グループのサンプルサイズがそろっていない計算を見て、グループ間の平方和の計算式がしっくり来た。
グループのサンプルサイズを掛ける意味が分かった気がした。
そろっているときは同じ値だけど、そろっていない場合は違う値を掛けるというのはとても筋が通っているし、そろっているときもサンプルサイズ倍していてよかったねという感じだ。
分散分析でグループのサンプルサイズがそろっていないのは計算上も研究上も得策ではないが、予期せぬことは起きるもので、そんなときでも計算できるのは大変助かる。
参考書籍
医学への統計学 8.2.1 一元配置分散分析
最新版はこちら。
