R統計検定2（一元配置の分散分析）

# Rで行う基本的な統計検定とその選択フローチャート2

このシリーズでは生物学で用いられている基本的な統計検定をRで行う際の手順と実際のコードを特集します。自身のデータセットの特徴と実験目的から適切な統計検定法を選択し、それをRで実行できるようになるための指針となれば幸いです。
第一回：平均値の差の検定（2群） 
**第二回：一元配置の分散分析**
第三回：一元配置の多重比較
第四回：二元配置の分散分析とその後の多重比較 
第五回：比率の差の検定
第六回：サンプルサイズの設計法
第七回：回帰分析

## 第二回　一元配置の分散分析

2.1 適切な検定法の選択
自分の取得した標本が正規分布に従うのか、標本同士の分散が等しいか、を検証して適切な統計検定法を選択できるような手順と必要となる検定のソースコード（R）を紹介します。
　　2.1.1 正規性の検証（Shapiro-Wilk test）
　　2.1.2 等分散性の検証（Bartlett test）
　　2.1.3 等分散性の検証（Levene test）

2.2 パラメトリック検定法のRソースコード
パラメトリック検定のソースコードを紹介します。分散分析の各検定法を紹介します。
　　2.2.1 one-way ANOVA（対応なし）
　　2.2.2 Repeated measures one-way ANOVA（対応あり）

2.3 ノンパラ版分散分析
ノンパラメトリック検定のソースコードを紹介します。ノンパラ版分散分析の各検定法を紹介します。
　　2.3.1 Friedman-test（対応あり）
　　2.3.1 Kruskal-Wallis rank sum test（対応なし）
　
### 2.1 適切な検定法の選択
一元配置の実験デザインにおける平均値の差の多群比較の検定についてです。以下に検定法選択のためのフローチャートをのせておきます。選択の上で重要なのが、測定値について(i)母集団が正規分布に従うと仮定できるか：正規性、と(ii)母集団の分散が標本間で等しいか：等分散性、の2つになります。

検定法選択のためのフローチャート
![スクリーンショット 2022-06-15 15.11.21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/352860/07e84e21-32b4-7376-a8be-d8e026c58df2.png)



#### 2.1.1 正規性の検証（Shapiro-Wilk test）
各群に対して正規性の確認を統計検定によって判断します。1000以下のサンプルに対しては一般にShapiro-Wilk testを用います。
Rで実行してみましょう（第一回で紹介したのと同じ）。
```R
score = c(57, 67, 23, 50, 52, 51, 48, 49, 54, 44, 31, 53, 40, 60)
shapiro.test(score)
```
結果
```
Shapiro-Wilk normality test

data:  score
W = 0.93902, p-value = 0.4059
```
以上から帰無仮説は棄却されず、データは正規分布に従うと仮定してよいと考えられます。
参照：https://data-science.gr.jp/implementation/ist_r_shapiro_wilk_test.html


#### 2.1.2 等分散性の検証（Bartlett test）
データが正規分布に従いそうと分かれば、次は各群の等分散性を検証する。もしデータに正規性がない場合または等分散性がない場合はノンパラメトリック手法を使う（フローチャート参照）。**2標本の等分散性の検定ではF-testを用いたが、多群の等分散性の検定にはBartlett　testを用いる。**

Rで実際に検定を行なってみましょう。帰無仮説は(H0)は各群の分散が等しいことです。
```R
score = c(301, 311, 325, 291, 388, 402, 325, 361, 287, 261, 238, 362, 197, 180, 178, 260, 247, 199, 179, 134, 163, 200, 209, 331, 192, 155, 234, 290, 175, 116, 285, 216, 237, 301, 343, 247, 316, 395, 324, 138, 245, 228, 214, 374, 235)
group=factor(rep(c("A", "B", "C", "D"), c(12, 10, 12, 11)))
bartlett.test(formula=score ~ group)
```
結果
```
Bartlett test of homogeneity of variances

data:  score by group
Bartlett's K-squared = 5.2629, df = 3, p-value = 0.1535
```
以上から帰無仮説は棄却されず、各群の分散は等しいと仮定してよいと考えられます。
すなわちパラメトリックな分散分析あるいは多重比較法を用いてもよいと考えられます。
参照：https://data-science.gr.jp/implementation/ist_r_bartlett_test.html

#### 2.1.3 等分散性の検証（Levene test）
データが正規分布に従わない場合の、等分散性の検証法としてLevene testを用いる。
Rで実際に検定を行なってみましょう。帰無仮説は(H0)は各群の分散が等しいことです。
```R
score = c(301, 311, 325, 291, 388, 402, 325, 361, 287, 261, 238, 362, 197, 180, 178, 260, 247, 199, 179, 134, 163, 200, 209, 331, 192, 155, 234, 290, 175, 116, 285, 216, 237, 301, 343, 247, 316, 395, 324, 138, 245, 228, 214, 374, 235)
group=factor(rep(c("A", "B", "C", "D"), c(12, 10, 12, 11)))
levene.test(score,group)
```
結果
```
Modified robust Brown-Forsythe Levene-type test based on the absolute deviations from the median

data:  score
Test Statistic = 1.8458, p-value = 0.154
```
以上から帰無仮説は棄却されず、各群の分散は等しいと仮定してよいと考えられます。


### 2.2 one-way ANOVA
母集団が正規分布に従う&各群が等分散と仮定できる場合、各群の平均値に差があるかを調べることができる。またデータに対応があるか（反復測定かどうか）を確認し、対応なしの場合は通常のone-way ANOVAを適用する。
対応ありの場合は、Mauchly’s sphericity testを行なって球面性を確認してRepeated measures one-way ANOVAを適用する。球面性が確認できない場合は、自由度を補正してANOVAを行う必要がある（調査中）。

またANOVAではどの群とどの群に差があるかは言及できない。群間の差を検出するためには、多重比較法を用いる必要がある。F統計量を用いる多重比較法（Scheffe, Games/Howell, Fisher's PLSD）を使う場合は、事前にANOVAを行なって有意差があることを確認しておく必要があるが、F統計量を用いない多重比較法（Dunnet, Tukey-Kramer, Bonferroni）を使う場合は必ずしも事前のANOVAは必要でない。

#### 2.2.1 one-way ANOVA（対応なし）
対応がない（反復測定でない）場合は、通常のone-way ANOVAを用いる。
Rでone-way ANOVAを実行してみる。帰無仮説は(H0)は各群の平均値に差がないことです。
群A	56, 48, 72, 60, 55
群B	60, 62, 76, 84
群C	78, 53, 62, 44, 90, 57
群D	77, 72, 83, 81, 91, 83
のデータを使うとする。
```R
score = c(56, 48, 72, 60, 55, 60, 62, 76, 84, 78, 53, 62, 44, 90, 57, 77, 72, 83, 81, 91, 83)
class = factor(rep(c("A", "B", "C", "D"), c(5, 4, 6, 6)))
anova(aov(score ~ class))
# または
summary(aov(score ~ class))
# または
oneway.test(score ~ class, var.equal=TRUE)
```
結果
```
Analysis of Variance Table

Response: score
          Df Sum Sq Mean Sq F value  Pr(>F)  
class      3 1629.2  543.06  3.9141 0.02708 *
Residuals 17 2358.6  138.74                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
```
以上から各群の測定値の平均値にはなんらかの差があると言える。

#### 2.2.2 Repeated measures one-way ANOVA（対応あり）
対応あり（反復測定をする）場合は、Repeated measures one-way ANOVAを行う。
ただし、球面性を仮定する必要があり、Repeated measures one-way ANOVAを行う前にMauchly’s sphericity testを行う。

参照：https://biostats.w.uib.no/test-for-sphericity-mauchly-test/
参照：https://www1.doshisha.ac.jp/~mjin/R/Chap_13/13.html

6個体に対して時間軸t1~t4で4回ある測定値を得た。このデータセットを使ってRでMauchly’s sphericity testを行なって球面性を確かめてみる。帰無仮説(Ho)は球面性が仮定できる（被験者内条件のすべての可能なペア間の差の分散が等しい）。
	　　t1	t2	t3	t4
個体1	56	60	78	77
個体2	48	62	53	72
個体3	72	76	62	83
個体4	60	84	44	81
個体5	55	72	90	91
個体6	49	65	57	83

```R
install.packages("car")
library(car)

t1 <- c(56, 48, 72, 60, 55, 49)
t2<-c(60, 62, 76, 84,72, 65)
t3<-c(78, 53, 62, 44, 90, 57)
t4<-c(77, 72, 83, 81, 91, 83)

neo <- matrix(c(t1,t2,t3,t4), nrow = 6, ncol = 4)
mauchly.test(lm(neo ~ 1), X = ~ 1)
```
結果
```
Mauchly's test of sphericity
	Contrasts orthogonal to
	~1


data:  SSD matrix from lm(formula = neo ~ 1)
W = 0.060575, p-value = 0.07235
```
以上から帰無仮説は棄却されず、球面性を仮定できると考えて良い。

正規性、等分散性、球面性の仮定を全てクリアできたらやっと分散分析（Repeated measures one-way ANOVA）を行うことができる。
Rで実行してみる。帰無仮説(H0)は時間軸で平均値に差がない。
```R
install.packages("car")　
library(car)

t1 <- c(56, 48, 72, 60, 55, 49)
t2<-c(60, 62, 76, 84,72, 65)
t3<-c(78, 53, 62, 44, 90, 57)
t4<-c(77, 72, 83, 81, 91, 83)
dat<-data.frame(A=factor(c(rep("t1", 6), rep("t2",6), rep("t3", 6),rep("t4", 6))), No= factor(rep(1:6, 4)),y=c(t1, t2, t3,  t4)) 
#Aに因子を格納、Noで対応ありであることを指定、yに各群のスコアを格納
summary(aov(y ~ A + No, dat))
```
結果
```
Df Sum Sq Mean Sq F value  Pr(>F)   
A            3 1926.8   642.3   6.024 0.00667 **
No           5  859.8   172.0   1.613 0.21669   
Residuals   15 1599.2   106.6                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
```
以上から帰無仮説は棄却され、個体間に差はなく、各群の測定値の平均値にはなんらかの差があると言える。

### 2.3 ノンパラ版分散分析
母集団が正規分布に従わない、または各群が等分散でない場合、ノンパラメトリックな検定手法を使って、多群の中で平均値に差があるか（ノンパラ版分散分析）、どの群に差があるか（ノンパラ多重比較）を検定することになります。ノンパラ版分散分析はFriedman-testとKruskal-Wallis rank sum testの2つが主に使われている。

#### 2.3.1 Friedman-test（対応あり）
3群以上の独立したサンプルの比較を行うノンパラメトリック検定の手法です。正規性を仮定できない場合や外れ値を多く含む場合、サンプルサイズが小さい場合などに用いられます。Friedman-testはデータ間に対応のある場合に用いられる。
参照：https://data-science.gr.jp/implementation/ist_r_friedman_test.html

Rで実行してみる。帰無仮説(H0)は各群の値の分布はそれぞれの群間で等しいこと。
```R
A=c(56, 48, 72, 60)
B=c(60, 62, 76, 84)
C=c(78, 53, 62, 44)
D=c(77, 72, 83, 81)
friedman.test(y=matrix(c(A, B, C, D), ncol=4))
```
結果
```
Friedman rank sum test

data:  matrix(c(A, B, C, D), ncol = 4)
Friedman chi-squared = 6, df = 3, p-value = 0.1116
```
以上から帰無仮説は棄却されず、各群の値の分布に差があるとはいえない。

#### 2.3.2 Kruskal-Wallis rank sum test（対応なし）
Kruskal-Wallis testは、3群以上の独立したサンプルの比較を行うノンパラメトリック検定の手法です。正規性を仮定できない場合や外れ値を多く含む場合、サンプルサイズが小さい場合などに用いられます。対応がないときに適用する。
参照：https://www.stats-guild.com/analytics/13251

Rで実行してみる。帰無仮説(H0)は各群の値の分布はそれぞれの群間で等しいこと。
```R
A=c(56, 48, 72, 60)
B=c(60, 62, 76, 84)
C=c(78, 53, 62, 44, 90, 57)
D=c(77, 72, 83, 81, 91, 83)
kruskal.test(x=list(A, B, C, D))
```
結果
```
Kruskal-Wallis rank sum test

data:  list(A, B, C, D)
Kruskal-Wallis chi-squared = 7.2766, df = 3, p-value = 0.06358
```
以上から帰無仮説は棄却されず、各群の値の分布に差があるとはいえない。