MENU
| 統計ER
  1. ホーム
  2. EZR
  3. R で一元配置分散分析を行う方法

R で一元配置分散分析を行う方法

EZR R 一元配置分散分析 分散分析

Rで、一元配置分散分析を step by step で計算してみた。

lm() と Anova() を使えばあっという間だが、具体的な一つ一つの計算を自分で組み立ててみるとどうか？

教科書の例題に沿って確認した。

＞＞もう統計で悩むのは終わりにしませんか？ 

↑期間・数量限定で無料プレゼント中！

目次

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

F値は13.52、p値は0.001945で、統計学的有意。

つまり、A1, A2, A3の三グループの平均は等しくないと言える。

ちなみに三グループの平均値は、以下の通り。

A3だけ低く見える。

> tapply(X.all, A.all, mean)
A1    A2    A3
70.25 80.50 58.50

＞＞もう統計で悩むのは終わりにしませんか？ 

↑1万人以上の医療従事者が購読中

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値も同一だ。

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 一元配置分散分析

医学への統計学 (統計ライブラリー)

医学への統計学 (統計ライブラリー)

Amazon

最新版はこちら。

医学への統計学 (統計ライブラリー)

医学への統計学 (統計ライブラリー)

Amazon

EZR R 一元配置分散分析 分散分析
よかったらシェアしてね！
  • URLをコピーしました！
  • URLをコピーしました！

リサーチクエスチョン探し？データ分析？論文投稿？、、、で、もう悩まない！

第1章：臨床研究ではなぜ統計が必要なのか？計画することの重要性
  • 推定ってどんなことをしているの？
  • 臨床研究を計画するってどういうこと？
  • どうにかして標本平均を母平均に近づけられないか？
第2章：研究目的をどれだけ明確にできるのかが重要
  • データさえあれば解析でどうにかなる、という考え方は間違い
  • 何を明らかにしたいのか？ という研究目的が重要
  • 研究目的は4種類に分けられる
  • 統計専門家に相談する上でも研究目的とPICOを明確化しておく
第3章：p値で結果が左右される時代は終わりました
  • アメリカ統計協会（ASA）のp値に関する声明で指摘されていること
  • そうは言っても、本当に有意差がなくてもいいの…？
  • なぜ統計専門家はp値を重要視していないのか
  • 有意差がない時に「有意な傾向があった」といってもいい？
  • 統計を放置してしまうと非常にまずい
第4章：多くの人が統計を苦手にする理由
  • 残念ながら、セミナー受講だけで統計は使えません。
  • インプットだけで統計が使えない理由
  • どうやったら統計の判断力が鍛えられるか？
  • 統計は手段なので正解がないため、最適解を判断する力が必要
第5章：統計を使えるようになるために今日から何をすれば良いか？
  • 論文を読んで統計が使えるようになるための5ステップ
第6章：統計を学ぶために重要な環境
  • 統計の3つの力をバランスよく構築する環境

『統計を身につける5つのステップ』を
通常1,650円 → 無料でプレゼント

以下のボタンをクリックして、画面に出てくる指示に従って、必要事項を記入してください。

通常価格：1,650円　→　無料でお届け

【無料】統計を使いこなすための秘策を学ぶ

この記事を書いた人

toukei-erのアバター toukei-er

統計 ER ブログ執筆者

元疫学研究者

統計解析が趣味

関連記事

コメント

コメント一覧 （2件）

コメントする

目次