統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

サンプルサイズ計算ファイル&サンプルデータファイル
全品半額セール中!5月5日まで。HHA SHOPへGO!

https://happyhappygk.base.ec/

統計ソフトRで一元配置分散分析一歩一歩

ブログランキングに参加しています。
まずはぽちぽちっとお願いします。
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
にほんブログ村 科学ブログ 数学へ

統計ソフトRで、一元配置分散分析を one by one でじっくり計算してみた。

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

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

教科書

医学への統計学 8.2.1 一元配置分散分析

最新版はこちら。

例題のデータ

例題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

一元配置分散分析一歩一歩

一元配置分散分析は、分散分析とは言うが、要因のグループ(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

結果はこちら。F値は13.52、p値は0.001945で、統計学的有意。つまり、A1, A2, A3の三グループの平均は等しくないと言える。

> 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 

ちなみに三グループの平均値は、以下の通り。A3だけ低く見える。

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

lm() を使ってみるとどうか?

Anova() を使うため、car パッケージをインストールしていなければ、先にインストールする。インストールの方法は以下の過去記事参照。

toukeier.hatenablog.com

X.all を従属変数側に、A.allを独立変数側に書いて チルダ (~)で結ぶ。分散分析表を Anova() 出力させる。

library(car)
lm1 <- lm(X.all ~ A.all)
Anova(lm1)

結果はこちら。A.all の Sum Sq(Sum of Square=平方和)が、SSAに一致している。また、Residuals の Sum Sq が SSE に一致しているのがわかる。F値もp値も同一だ。

> 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

EZRで行うとどうか?(2020年9月13日追記)

データセットに書き出す

上記の例題データを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で簡単統計】【動画】

youtu.be

まとめ

統計ソフトRで、lm() を使えば一瞬で出力できる一元配置分散分析を、一つ一つ順を追ってスクリプトに書いて計算してみた。こうやると簡単な統計手法も理解が深まる。

今回も思ったのだが、グループ間の平方和で、グループのサンプルサイズ倍するところが、何となくしっくりこない。こういう計算なのだと頭ではわかっているが、何となく変な感じがする。グループ間の平均と全体平均の平方和だけのほうがグループ間のばらつきを表しているような気がしてしまうのだ。

ちなみに、もしあなたが、分散分析の結果だけでは気が済まずに、どのグループ間に差があるのか知りたいと思うなら、最初から多重比較検定をすべきであると付け加えて、今回は終わる。

toukeier.hatenablog.com

2019年1月10日追記:アンバランスの場合は?

グループのサンプルサイズがそろっていない場合はどうしたらいいのか?上記はサンプルサイズがそろっていることが前提だった。サンプルサイズがそろっていないアンバランスなデータの場合、スクリプトをどう変えるのか?

例題のデータ

同じ教科書の例題8.3を使う。データは以下の通り。A1からA4の4つのグループで、サンプルサイズは6,4,6,5とそろっていない。

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)

結果はこちら。A.allのSum SqがSSAと、ResidualsのSum SqがSSEと、そしてF値、p値が、それぞれ一歩一歩計算したものと同じことがわかる。

> 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 ‘ ’ 

まとめ その2

グループのサンプルサイズがそろっていない計算を見て、グループ間の平方和の計算式がしっくり来た。

グループのサンプルサイズを掛ける意味が分かった気がした。

そろっているときは同じ値だけど、そろっていない場合は違う値を掛けるというのはとても筋が通っているし、そろっているときもサンプルサイズ倍していてよかったねという感じだ。

分散分析でグループのサンプルサイズがそろっていないのは計算上も研究上も得策ではないが、予期せぬことは起きるもので、そんなときでも計算できるのは大変助かる。