【R】一元配置分散分析、多重比較検定

Rで統計解析するための備忘録的なやつです。

#使用するデータ
今回使うデータは、これ。

**ChickWeight**

ヒヨコの体重に及ぼすエサの影響を調べたデータらしいです。

#必要なパッケージ
**dplyr**と**ggplot2**を読み込むための**tidyverse**パッケージと、診断プロットを作るための**ggfortify**パッケージを読み込んでおく。

```R
#ライブラリの読み込み
library(tidyverse)
library(ggfortify)
```

#データの読み込みと確認

データを読み込んで、**str()**関数でデータの概要をチェック。

```R
#データの読み込みと確認
DF <- chickwts
str(DF)
```
数値のweight（つまり重さ）、と6水準のfeed（つまり6種類のエサ）のデータであることを把握。

#ひとまずプロット

**ggplot2**でとりあえず、プロット。
箱ひげ図を採用。

```R
#とりあえずプロットしてみる
ggplot(DF, aes(x = feed, y = weight)) +
  geom_boxplot()
```
実行結果がこれ。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/bea2fafd-27bd-c3e4-3761-a52ccb293ede.png)

ここから読み取れること。
エサによって体重が違うのが明らか。
caseinの効果が一番大きそう。
体重の増加はcasein > sunflower > meatmeal > soybean > linseed > horsebeanになってそう。

#分散分析モデルを作る

Rで分散分析モデルを作るときは、**lm()**関数で作ったモデルを**anova()**関数にかけて結果を出すやり方もあれば、**lm()**と**anova()**を組み合わせた分散分析のための**aov()**関数を使うやり方もあります。
ここでは、楽な**aov()**を使ってみます。

```R
#分散分析モデル構築
AOV <- aov(weight ~ feed, data = DF)
```
#前提条件を満たしているかどうか？

**ggfortify**パッケージの**autoplot()**関数で診断プロットを作ります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/e044eac9-38da-d6ba-6759-c99d696f6228.png)

見た感じ、特に問題はなさそうです。

#一元配置分散分析の結果の出力

**summary()**関数を使って結果を出力します。

```R
#一元配置分散分析の結果
summary(AOV)
```
実行結果

```R
> #一元配置分散分析の結果
> summary(AOV)
            Df Sum Sq Mean Sq F value   Pr(>F)    
feed         5 231129   46226   15.37 5.94e-10 ***
Residuals   65 195556    3009                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
```
ここでの帰無仮説は「6種類のエサを与えた群のいずれの平均は同じである」です。
なので、対立仮説は「6種類のエサの群の平均のいずれかのあいだに差がある」です。
分散分析で計算する*F*値は、グループ間の分散とグループ内の分散の比を表していて、グループ間の分散がグループ内の分散より大きいと*F*値が大きくなって、*p*値が小さくなり、「平均値は同じである」という帰無仮説が棄却されます。

結果を見ると、エサの効果は0.1%水準で統計的に優位と言えます。
ただ、分散分析だけでは、どの群とどの群の平均値間に有意差があるかはわかんないです。

#多重比較

多重比較でそれぞれの群の平均値間の差の検定を行います。
テューキーのHSD（honestly significant difference）法で多重比較検定を行ってみました。

```R
#テューキーのHSD法による多重比較
HSD <- TukeyHSD(AOV, "feed")
HSD #検定結果の表示
```
実行結果

```R
> #テューキーのHSD法による多重比較
> HSD <- TukeyHSD(AOV, "feed")
> HSD #検定結果の表示
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = weight ~ feed, data = DF)

$feed
                           diff         lwr       upr     p adj
horsebean-casein    -163.383333 -232.346876 -94.41979 0.0000000
linseed-casein      -104.833333 -170.587491 -39.07918 0.0002100
meatmeal-casein      -46.674242 -113.906207  20.55772 0.3324584
soybean-casein       -77.154762 -140.517054 -13.79247 0.0083653
sunflower-casein       5.333333  -60.420825  71.08749 0.9998902
linseed-horsebean     58.550000  -10.413543 127.51354 0.1413329
meatmeal-horsebean   116.709091   46.335105 187.08308 0.0001062
soybean-horsebean     86.228571   19.541684 152.91546 0.0042167
sunflower-horsebean  168.716667   99.753124 237.68021 0.0000000
meatmeal-linseed      58.159091   -9.072873 125.39106 0.1276965
soybean-linseed       27.678571  -35.683721  91.04086 0.7932853
sunflower-linseed    110.166667   44.412509 175.92082 0.0000884
soybean-meatmeal     -30.480519  -95.375109  34.41407 0.7391356
sunflower-meatmeal    52.007576  -15.224388 119.23954 0.2206962
sunflower-soybean     82.488095   19.125803 145.85039 0.0038845
```

#最後の仕上げの図
こんな感じで。

```R
#最後の仕上げの図
ggplot(DF, aes(x = feed, y = weight, colour = feed)) + #feedによる色分け
  geom_boxplot() +
  geom_point(size = 3, alpha = 0.5) + #点の大きさ（size）、透明度（alpha）
  xlab("Feed type") +
  ylab("Weight (g)") +
  scale_colour_discrete(guide = FALSE) + #凡例を消す
  theme_bw()
```
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/5e0848ff-94f6-0881-301f-c2d38d16f5bd.png)

#最終的なスクリプト

```R
#ニワトリの成長速度に対する一元配置分散分析（one-way ANOVA）

#Rをきれいにする
rm(list = ls())

#ライブラリの読み込み
library(tidyverse)
library(ggfortify)

#データの読み込みと確認
DF <- chickwts
str(DF)

#とりあえずプロットしてみる
ggplot(DF, aes(x = feed, y = weight)) +
  geom_boxplot()

#分散分析モデル構築
AOV <- aov(weight ~ feed, data = DF)

#モデルの診断プロットを確認する
autoplot(AOV)

#一元配置分散分析の結果
summary(AOV)

#テューキーのHSD法による多重比較
HSD <- TukeyHSD(AOV, "feed")
HSD #検定結果の表示

#最後の仕上げの図
ggplot(DF, aes(x = feed, y = weight, colour = feed)) + #feedによる色分け
  geom_boxplot() +
  geom_point(size = 3, alpha = 0.5) + #点の大きさ（size）、透明度（alpha）
  xlab("Feed type") +
  ylab("Weight (g)") +
  scale_colour_discrete(guide = FALSE) + #凡例を消す
  theme_bw()
```