Rで統計学を学ぶ(6)

分散分析(ANOVA, Analysis Of Variance) とは、３つ以上の群の平均値差の検定や、条件(の組み合わせ)の違いによって母平均が異なるかどうかを検定する方法です (注意:「分散」という名前に惑わされないように)

注: ここでは7.1.2節の説明を先に行います.

データにおいて変数がとる値は、以下の『要因』によるものと「誤差」によるものがあるとみなせます:

• 要因: 母平均に違いをもたらす原因　　　　　
• 水準: ある要因の中に含まれる個々の条件

# A, B, C, Dは「対応のない一元配置分散分析」の例からとった。1要因4水準
A <- c(15,9,18,14,18)
B <- c(13,8,8,12,7)
C <- c(10,6,11,7,12)
D <- c(10,7,3,5,7)

全データ <- cbind(A,B,C,D)	
群平均 <- colMeans(全データ)	# 水準ごとの平均のベクトル
# 水準ごとの平均のベクトルを行数だけ並べた行列を作る
群平均行列 <- matrix(rep(群平均, nrow(全データ)), nrow=nrow(全データ), 
                     ncol=ncol(全データ),byrow=TRUE)

群内 <- 全データ - 群平均行列	# 各データと、その群内の平均との差
群内平方和 <- sum(群内^2)		# 　群内平方和

分母 <- 群内平方和 / ( (nrow(全データ)-1)*ncol(全データ) )

一方、「要因によるばらつき」は、群間のばらつきで表されます。それは「全データの平均と各群の平均の差」の2乗の和を、群間の自由度(群の数-1、つまり水準数-1)で割った値と計算されます:

# 全データの和の平均
全体平均 <- sum(colMeans(全データ))/ncol(全データ)	

# 各群ごと「群平均-全平均」からなる行列
群間 <- matrix(rep(群平均 - 全体平均, nrow(全データ)),nrow(全データ), byrow=TRUE) 
群間平方和 <- sum(群間^2)	# 　群内平方和

分子 <- 群間平方和 / (ncol(全データ) - 1)

> 分子
[1] 61.33333
> 分母
[1] 8.625
> ( F <- 分子/分母 )
[1] 7.111111

検定の結果として表示される分散分析表の構成:

要因	自由度(Df)	平方和(Sum Sq)	平均平方(Mean Sum Sq)	F値(F value)	P値	評価
要因名	水準数-1	群間平方和	群間平方平均	F値	p値	印(表の下に説明)
残差(Residuals)	水準数×(１群あたりのデータ数-1)	群内平方和	群内平方平均
合計	データ数-1

評価の印: 星印は有意(***は0.1%、**は1%、*は5%水準)、.は10%(「有意傾向」という)を表す

Rを用いた計算方法の紹介

1. 帰無仮説と対立仮説の設定
   　 ・帰無仮説 H₀：4群の母平均は等しい
   　 ・対立仮説 H₁：4群の母平均は等しくない
   注意: 4群のうち1つだけ母平均が異なるような場合でも対立仮説は成立
2. 検定統計量の選択 --- 次のFを検定統計量とする:
   これは帰無仮説（すべての群の母平均が等しいとき）のもとで、F分布に従う
   F分布の自由度は２種類：分子と分母の自由度
   F分布の描画: df(x,分子の自由度,分母の自由度) 関数を使用

   curve(df(x, 3, 16), 0, 5)
   

   
   
3. 有意水準αの決定: ここでは有意水準は5%、つまりα=0.05とする
   F分布のグラフから明らかなように、F分布は正値を取る⇒　分散分析は常に片側検定
4. 検定統計量の実現値を求める: 以下の３種類を紹介
   1. oneway.test関数 : 一元配置分散分析のみ実行可能
   2. aov関数: 最も一般的
   3. anova関数: 複数のモデルの比較など高度な分析に対応
   Rを用いた計算:

   > A <- c(15,9,18,14,18)
   > B <- c(13,8,8,12,7)
   > C <- c(10,6,11,7,12)
   > D <- c(10,7,3,5,7)
   
   > 統計テスト2 <- c(A,B,C,D)		# 4群のデータをまとめる
   > 指導法 <- c(rep("A",5),rep("B",5),rep("C",5),rep("D",5))
   > 指導法2 <- factor(指導法)		# 要因型ベクトルに変換
   > 指導法2 
    [1] A A A A A B B B B B C C C C C D D D D D
   Levels: A B C D
   

   以上で準備ができたので、「帰無仮説の棄却/採択の決定」も含めて検定統計量の実現値を求める:

   oneway関数

   # var.equal=TRUEにより母分散の等質性を仮定
   > oneway.test(統計テスト2 ~ 指導法2, var.equal=TRUE )
           One-way analysis of means
   data:  統計テスト2 and 指導法2 
   F = 7.1111,   num df = 3,   denom df = 16, p-value = 0.002988	
   # F値は7.11, p値は0.002 < 0.05 
   

   母分散が等質と仮定したoneway.test関数の結果、F値は7.11、p値が0.003であるから、帰無仮説が棄却される

   aov関数

   > aov(統計テスト2 ~ 指導法2)
   Call:   aov(formula = 統計テスト2 ~ 指導法2)
   Terms:
                        指導法2     Residuals
   Sum of Squares       184             138
   Deg. of Freedom        3              16
   Residual standard error: 2.936835 
   Estimated effects may be unbalanced
   > summary(aov(統計テスト2~指導法2))
               Df    Sum Sq    Mean Sq       F value     Pr(>F)   
   指導法2      3   184.000    61.333        7.1111     0.002988 **
   Residuals   16   138.000      8.625                    
   ---
   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
   

   aov関数により分散分析表が表示される。この表では、F値が7.11、p値が 0.003 であるから、帰無仮説が棄却される。

   anova関数

   > anova(lm(統計テスト2 ~ 指導法2))
   Analysis of Variance Table
   
   Response: 統計テスト2
             Df Sum Sq Mean Sq  F value   Pr(>F)   
   指導法2    3    184  61.333  7.1111 0.002988 **
   Residuals 16    138   8.625                    
   ---
   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
   

   aov関数と同様に、分散分析表が出力される。F値は7.11、p値は 0.003 であるから、帰無仮説が棄却される。

   実際に使うのは、上の3つの方法のどれを使ってもよい。すべて帰無仮説が棄却されたことから、検定の結果、4つの指導方法の間に有意な平均値差が見られたことが導ける。

仮定: 各群のデータ数（サンプルサイズ)nが等しく、各群の母分散も等しい(分散の等質性)
検定統計量q:
Rを用いて説明：指導法AとDの比較

> mean(A)
[1] 14.8
> mean(D)
[1] 6.4
> nrow(全データ)
[1] 5
> 群内平均平方
[1] 8.625
> q <- abs(mean(A)-mean(D))/sqrt(群内平均平方/nrow(全データ))
> q			# 検定統計量q
[1] 6.395651

# qtukey(下側確率、平均値の数、群内の自由度）
> qtukey(0.95,4,16)
[1] 4.046093
> qtukey(0.05,4,16,lower.tail=FALSE)		# 上側確率
[1] 4.046093

> TukeyHSD(aov(統計テスト2 ~ 指導法2))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = 統計テスト2 ~ 指導法2)

$指導法2
         diff         lwr            upr           p adj
B-A  	 -5.2  	-10.514108 	  0.1141085  	 0.0562227
C-A   	 -5.6  	-10.914108	 -0.2858915  	 0.0371222
D-A 	 -8.4 	-13.714108	 -3.0858915  	 0.0017736
C-B  	 -0.4  	 -5.714108  	  4.9141085 	 0.9963241
D-B   	 -3.2  	 -8.514108  	  2.1141085 	 0.3446966
D-C  	 -2.8  	 -8.114108 	  2.5141085 	 0.4561325

> 好意度 <- c(7,8,9,5,6,5,4,7,1,3,8,6,7,2,5)
> 科目 <- factor(c(rep("線形代数",5),rep("微分積分",5),rep("確率統計",5)))
> summary(aov(好意度~科目))
            Df Sum Sq Mean Sq F value Pr(>F)
科目         2  22.53  11.267   2.641  0.112
Residuals   12  51.20   4.267     

例題: 以下の表は5名の学生に、線形代数、微分積分、確率統計に対する好意度の評定データである。 これら3つの科目に対する好意度に差があるだろうか？（念のため：これは１要因３水準)

学生	線形代数	微分積分	確率統計
田中	7	5	8
岸	8	4	6
大引	9	7	7
吉川	5	1	2
荻野	6	3	5

Rを用いた計算方法の紹介

1. 帰無仮説と対立仮説の設定
   　
   • 帰無仮説 H₀：3科目の好意度の母平均は等しい 　
   • 対立仮説 H₁：3科目の好意度の母平均は等しくない
2. 検定統計量の選択 --- 次のFを検定統計量とする:
   これは帰無仮説（すべての群の母平均が等しいとき）のもとで、F分布に従う
3. 有意水準αの決定: ここでは有意水準は5%、つまりα=0.05とする
   F分布のグラフから明らかなように、F分布は正値を取る⇒　分散分析は常に片側検定
4. 検定統計量の実現値を求める: ここでは、最も一般的なaov関数を用いる:

   > 好意度 <- c(7,8,9,5,6,5,4,7,1,3,8,6,7,2,5)
   > 科目 <- factor(c(rep("線形代数",5),rep("微分積分",5),rep("確率統計",5)))
   > 人 <- factor(rep(c("田中","岸","大引","吉川","萩野"),3))
   > summary(aov(好意度~科目+人))			# 好意度を科目と人の観点から分析する
              　Df  Sum Sq   Mean Sq    F value    Pr(>F)    
   科目         2   22.533    11.267     14.696     0.002095 ** 		# 検定統計量F = 14.7
   人           4   45.067    11.267     14.696     0.000931 ***
   Residuals    8    6.133     0.767                     
   ---
   Signif. Codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
   

   [『人』がなぜ必要？もしくは「好意度~科目+人」の意味]
   まず『人』を含めないaovの結果を示しましょう---これはデータに「対応がない」ときの分析と同じものです:

   > summary(aov(好意度~科目))
               Df Sum Sq Mean Sq F value Pr(>F)
   科目         2  22.53  11.267   2.641  0.112
   Residuals   12  51.20   4.267     
   

   これが示すように、科目の効果だけでは好意度の平均に有意な差が認められません。 そもそもこのデータにおいて『対応がある』というのは、好意度のデータが個人別だったからです。 言い換えれば、「好意度」データを説明する要因として『科目』と『人』があり、『人』を考慮することでデータに「対応がある」と認められたのです。それをRでは 好意度 ~ 科目+人と表したのです。
5. 帰無仮説の棄却か採択かの決定: 結果の出力の「科目」の行において、p値が 0.0021 と得られているので、5%の有意水準で、 科目間に有意な差があることが認められた。