【分散分析】3群以上の平均の差を検定しよう

こんにちは、たくろー(@takuro_data )です！

前々回、前回で１群、２群の仮説検定について学びました。

【仮説検定】仮説検定の用語を確認して実施してみよう

【仮説検定】2群の平均の差を検定をしよう

今回は、３群以上の平均の差を検定できる分散分析について、学びます。

要点

一元配置分散分析は、実験の目的となる要因（因子）が１つの場合の、もっとも基本的な分散分析。
二元配置分散分析は、それぞれの要因の主効果のほかに、交互作用の存在についても検定。
Rコマンドはどちらもsummary(anova())。

分散分析とは

３群以上の平均値の差を検定します。目的となる要因効果が、誤差効果よりも大きいとき、その分散比であるF値が大きくなることを利用します。

グループ（群）をつくる処理条件や水準が３つ以上の場合でも、平均の差を検定できる（t検定は２グループの場合だけ）
複数要因の交互作用（相乗効果や相殺効果）を検定できる

一元配置分散分析

実験の目的となる要因（因子）が１つの場合の、もっとも基本的な分散分析です。人為的に処理条件（水準）を変化させた要因が、実験の結果に影響を与えているかどうか確かめます。

データ全体のバラツキ（総変動）は、目的となる要因の効果によるバラツキ（群間変動）と、目的以外の要因である誤差の効果によるバラツキ（郡内変動）から構成されます。

目的となる要因によって生じる不偏分散を分子、誤差によって生じる不偏分散を分母とした検定統計量F値を求め、F分布を用いて、上側確率の片側検定をします。

RのサンプルデータInsectSprays（昆虫スプレーの有効性）で実行してみます。異なる殺虫剤で処理した農業実験単位での昆虫の数です。

data("InsectSprays")
str(InsectSprays)
plot(count ~ spray, data = InsectSprays)
#boxplot(InsectSprays[,1]~InsectSprays[,2])

> str(InsectSprays)
'data.frame':	72 obs. of  2 variables:
 $ count: num  10 7 20 14 14 12 10 23 17 20 ...
 $ spray: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...

#一元配置分散分析
summary(aov(count~spray,data=InsectSprays))

> summary(aov(count~spray,data=InsectSprays))
Analysis of Variance Table

Response: count
          Df Sum Sq Mean Sq F value    Pr(>F)    
spray      5 2668.8  533.77  34.702 < 2.2e-16 ***
Residuals 66 1015.2   15.38                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

sprayの要因は6で、母平均の差があるか検定してます。検定結果をみると、sprayの自由度が5（群間平均数6-総平均1）で、群間変動は2668.8で、不偏分散は533.77(群間変動÷自由度)です。郡内変動は1015.2で、自由度は66（データ数72-群間平均の数6）で、不偏分散は15.38です。検定統計量F値は34.702（533.77÷15.38）で、p値が0.00%なので、有意水準5%で、群間の母平均には差がない（要因効果がない）という帰無仮説を棄却し、群間の母平均には差がある（要因効果がある）という対立仮説を採択します。

多群の等分散の検定

分散分析も、等分散が前提なので、等分散を帰無仮説とした検定（Bartlett検定）を実施するのが望ましいです。InsectSpraysデータで実行してみます。

#bartlett.test(InsectSprays$count, InsectSprays$spray)
bartlett.test(count ~ spray, data = InsectSprays)

> bartlett.test(count ~ spray, data = InsectSprays)

	Bartlett test of homogeneity of variances

data:  count by spray
Bartlett's K-squared = 25.96, df = 5, p-value = 9.085e-05

自由度5（要因数6-1）のカイ二乗分布に従う統計量が25.96で、p値が0.00%なので、各群の母分散には差がない（等分散）という帰無仮説は棄却され、各群の母分散には１対以上で差がある（異分散）という対立仮説を採択します。

箱ひげ図から、C,D,Eが等分散に見えるので、抽出して検定してみます。

InsectSprays_cde <- rbind(InsectSprays[InsectSprays$spray == "C",],InsectSprays[InsectSprays$spray == "D",],InsectSprays[InsectSprays$spray == "E",])
bartlett.test(count ~ spray, data = InsectSprays_cde)

> bartlett.test(count ~ spray, data = InsectSprays_cde)

	Bartlett test of homogeneity of variances

data:  count by spray
Bartlett's K-squared = 1.504, df = 2, p-value = 0.4714

p値が47.1%で有意水準5%で等分散の帰無仮説は棄却されませんでした。InsectSpraysについてグループ分けして一元配置分散分析をしたほうが望ましそうです。

対応のある一元配置分散分析

各条件下で、同一個体を測定している場合には、対応のある一元配置分散分析を実施することで、個体差が考慮された検定が可能になります。

対応のない場合の群内変動（誤差によるバラツキ）は、個体差による変動（被験者間変動）を含んでいますが、対応のある場合にはそれぞ取り除くことができます。

対応のない一元配置分散分析で用いたデータを加工して、対応のある一元配置分散分析を実行してみます。

#各水準のデータ数確認
library(dplyr)
df <- group_by(InsectSprays_cde, spray)
summarize(df, n=n())
#対応のある情報を付加
InsectSprays_cde["NO"] <- factor(rep(1:12,3))
head(InsectSprays_cde)

#対応のある分散分析実行
#summary(aov(count ~ A + spray, InsectSprays_cde))
anova(aov(count ~ spray + NO , InsectSprays_cde))

#対応のない分散分析実行
anova(aov(count ~ spray, InsectSprays_cde))

> summarize(df, n=n())
# A tibble: 3 x 2
  spray     n
  <fct> <int>
1 C        12
2 D        12
3 E        12
> head(InsectSprays_cde)
   count spray NO
25     0     C  1
26     1     C  2
27     7     C  3
28     2     C  4
29     3     C  5
30     1     C  6

> anova(aov(count ~ spray + NO , InsectSprays_cde))
Analysis of Variance Table

Response: count
          Df Sum Sq Mean Sq F value   Pr(>F)   
spray      2 48.167 24.0833  6.6367 0.005555 **
NO        11 65.000  5.9091  1.6284 0.159021   
Residuals 22 79.833  3.6288                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> anova(aov(count ~ spray, InsectSprays_cde))
Analysis of Variance Table

Response: count
          Df  Sum Sq Mean Sq F value   Pr(>F)   
spray      2  48.167 24.0833  5.4873 0.008763 **
Residuals 33 144.833  4.3889                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

対応のある一元配置分散分析結果を、対応のない場合と比べると、Residuals（群内変動）の一部がNO（被験者間変動）に分離されている。sprayの効果の差がないという帰無仮説は棄却され、被験者間の母平均は等しいという帰無仮説は棄却されない結果となりました。

二元配置分散分析

目的となる要因が２つ（二元配置）以上の分散分析では、それぞれの要因の主効果のほかに、交互作用の存在についても検定することができます。

総変動は、要因効果による変動（群間変動）と、誤差効果による変動（郡内変動）とに分けることができます。また、要因効果による変動は、それぞれの主効果による変動と交互作用による変動とに分けることができます。

２因子でそれぞれ２水準の場合のサンプルデータを作成して実行してみます。

#サンプルデータ作成
a1 <- c(0,2,8,10)
a2 <- c(6,8,9,13)
A <- factor(c(1,1,1,1,2,2,2,2))
B <- factor(c(1,1,2,2,1,1,2,2))
y <- c(a1,a2)
sample_anova <- data.frame(y,A,B)
sample_anova

> sample_anova
   y A B
1  0 1 1
2  2 1 1
3  8 1 2
4 10 1 2
5  6 2 1
6  8 2 1
7  9 2 2
8 13 2 2

#交互作用なし二元配置分散分析実行
summary(aov(y~A+B,data=sample_anova))
#交互作用あり二元配置分散分析実行
summary(aov(y~A*B,data=sample_anova))

> summary(aov(y~A+B,data=sample_anova))
            Df Sum Sq Mean Sq F value  Pr(>F)   
A            1     32    32.0   7.273 0.04295 * 
B            1     72    72.0  16.364 0.00987 **
Residuals    5     22     4.4                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> summary(aov(y~A*B,data=sample_anova))
            Df Sum Sq Mean Sq F value Pr(>F)  
A            1     32    32.0   9.143 0.0390 *
B            1     72    72.0  20.571 0.0105 *
A:B          1      8     8.0   2.286 0.2051  
Residuals    4     14     3.5                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

要因AとBの主効果がないという帰無仮説は5%水準で棄却されましたが、A×Bの交互作用効果がないという帰無仮説は棄却されませんでした。