第4章は統計を扱います。
今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。
シロート統計学はEZRを使った統計分析をわかりやすく解説されています。
第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか？といった視点で進めていきます。
今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください！
今回は一元配置分散分析を紹介します
まず一元配置分散分析（ANOVA）についてはハルさんのサイトをご参照ください。
また1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。
もし1〜4まででわからない部分があれば確認してみてください。
1.準備
第4章は毎回ExcelデータをダウンロードしてRを使うのでプロジェクトで管理して行うことを勧めています。
ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。
これを前提に次から進めていきます。
2.スクリプトファイルの作成
次にRのコードを書くためのスクリプトファイルを作ります。
3.データのダウンロード
今回もハルさんのサイトのデータを使わせていただきます。
デモデータ（ANOVA）
この章ではRを使ってダウンロードしています。
download.file(url = “ファイルのURL”,destfile = “保存したい場所/ファイル名”)
urlはデモデータで右クリック → リンクのアドレスをコピー
destfileは保存場所と保存のファイル名を指定します。
実際のコードは以下になります。
前回のコードのURL（" "の中）とdestfileのdata/以降を変更するだけでOKです。
#データのダウンロード url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-anova.xlsx" destfile = "data/demo-anova.xlsx" download.file(url, destfile)
dataフォルダにダウンロードできたことを確認します。
4.データの読み込み
データを読み込みます。
今回は【4-0】第4章を進めていく上での準備で行った方法で進めます。
View Fileでデータを確認します。
データが入っているセルを確認します。
B2からC22までデータが入っています（B2:C62と表記）
次にImport Datasetでデータを取り込みます。
Import画面ではName, Sheet，Rangeを指定します。
Name：ハルさんのサイトではgripでしたがgrip_anovaとします（大文字・小文字は別物とされます）
Sheet：このExcelは1つしかデータがないのでDefaultのままでOK
Range：先ほど確認したB2:C62
Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。
library(readxl)
grip_anova <- read_excel("data/demo-anova.xlsx",
range = "B2:C62")
View(grip_anova)データが正しく入っていることを確認します。
これでデータの取り込みは完成です。
5.正規分布の確認
ハルさんのサイトに沿って正規分布の確認を行います。
正規分布の確認は【4-1】Rでt検定を行う方法で紹介しています。
以下のコードがイメージできない場合はまずこちらの記事をご参照ください。
まずヒストグラムで分布を確認します。
1つずつグラフを作るのは手間がかかるので、まとめてグラフを作ります。
ヒストグラムはgeom_histgram関数を使います。
グラフが重なっていると見にくいのでfacet.gridでカテゴリーごとにグラフを分けます。
facet.gridに関しては【3-4】Rのggplot2でヒストグラムを作るgeom_histogram関数で紹介しています。
library(tidyverse) ggplot(data = grip_anova) + geom_histogram(aes(x = grip, fill = category), bins = 5) + facet_grid(category ~ .)
6.シャピロ・ウィルク検定
次にハルさんのサイトでは正規性の検定を行っています。
正規性の検討はshapriro.test関数を3つのカテゴリーそれぞれに行う必要があります。
今回はsplit関数とmap関数を使いまとめて行ってみます。
下の図は【4-1】Rでt検定を行う方法で使った方法ですが今回のコードと見比べてみてください。
コードがほとんど同じです。コードを使うと他でも役に立つことがわかります。
grip_anova %>% split(.$category) %>% map(~shapiro.test(.$grip))
まとめて結果が出ました。すべて正規分布であることが棄却されませんでした。
7.等分散性の確認
次に等分散性の確認を行います。
3群以上で等分散性を検定するにはBartlett検定を使います。
bartlett.test(数値の列 ~ グループ)
data = を使えば$は必要ありません。どちらでも大丈夫です。
bartlett.test(grip_anova$grip ~ grip_anova$category)
bartlett.test(grip ~ category, data = grip_anova)
8.分散分析（ANOVA）を行う
それでは本題の一元配置分散分析と多重比較を行います。
実は一元配置分散分析と重回帰分析（t検定も）は同じ分析方法です（詳しく知りたい方は一般線形モデルで検索してみてください）。出てきた結果を分散分析の形式で表示するか重回帰分析の結果で表示するかの違いだけです。
そのため一元配置分散分析を行うには分散分析のaov関数か重回帰分析のlm関数で行います。
ただその後に行う多重比較ではlm関数が使えないものもあるので基本はaov関数で大丈夫です。
ただその後に行う多重比較ではlm関数が使えないものもあるので基本はaov関数で大丈夫です。
どちらの関数もまずモデルを作りサマリーを表示するというのが一般的な流れです。
全体的な流れは以下の通りです。
sumamry関数だけ結果の表示が変わるので注意が必要です。
aov関数の場合
モデル名：model_grip_anova（好きな名前でOK）
従属変数：grip
独立変数：category
データ名：grip_anova
model_grip_anova <- aov(grip ~ category, data = grip_anova) summary(model_grip_anova)
lm関数の場合
モデル名：model_grip_lm（好きな名前でOK）
従属変数：grip
目的変数：category
データ名：grip_anova
model_grip_lm <- lm(grip ~ category, data = grip_anova) summary.aov(model_grip_lm)
9.多重比較を行う
ハルさんのサイトでは多重比較でTukeyを紹介していますが他も紹介します。
Tukeyだけ関数が違います。
Tukeyだけ関数が違います。
Tukeyの場合
Tukeyの場合は先程のaov関数の結果を使います。
lm関数の結果は使えませんので注意が必要です。
Tukeyの場合は先程のaov関数の結果を使います。
lm関数の結果は使えませんので注意が必要です。
TukeyHSD(model_grip_anova)
Bonferroniの場合
pairwise.t.test関数はモデルの結果を使いません。
pairwise.t.test関数はモデルの結果を使いません。
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "bonferroni")
Holmの場合
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "holm")
Hochbergの場合
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "hochberg")
10.まとめ
今回は一元配置分散分析を紹介しました。
今回は割愛しましたが、一元配置分散分析を行うときはlongデータである必要があります。
今回は割愛しましたが、一元配置分散分析を行うときはlongデータである必要があります。
Excelでデータ収集するときに注意するか、wideデータの場合はlongデータに直してから分析を行ってください。
longデータとwideデータについてはで【2-5】Rでデータを集計するのに便利なlongデータとgather関数で
紹介しています。
次回はKruskal-Wallis検定を行います。
（追記）グラフの作成
EZRでは検定と同時にグラフを作成しますが、Rでは直接作成します。
【4-1】t検定でグラフを使ったときのコードがそのまま使えます。
グラフを作るために必要な要約のコードもそのまま使います。
コードの違いはわかったでしょうか？
t検定で使ったデータであるgripを今回のデータであるgrip_anovaに変更しただけです。
それ以外のコードの中身は全く一緒です！
11.今回のコード
今回使ったコードのまとめになります。
全記事の目次はこちらをご参照ください
紹介しています。
次回はKruskal-Wallis検定を行います。
（追記）グラフの作成
EZRでは検定と同時にグラフを作成しますが、Rでは直接作成します。
【4-1】t検定でグラフを使ったときのコードがそのまま使えます。
グラフを作るために必要な要約のコードもそのまま使います。
#データの要約
grip_anova_summary <-
grip_anova %>%
group_by(category) %>%
summarize(平均 = mean(grip),
標準偏差 = sd(grip),
'0%' = quantile(grip, 0),
'25%' = quantile(grip, 0.25),
'50%' = quantile(grip, 0.5),
'75%' = quantile(grip, 0.75),
'100%' = quantile(grip, 1),
n = n())
grip_anova_summary
#グラフ化
ggplot(data = grip_anova_summary)+
geom_bar(aes(x = category, y = 平均), color = "black", fill = "gray", stat = "identity") +
geom_errorbar(aes(x = category, ymin = 平均 - 標準偏差, ymax = 平均 + 標準偏差), width = 0.1) +
labs(y = "grip") +
theme_classic(base_size = 15) +
scale_y_continuous(expand = c(0,0), limits = c(0,50))
コードの違いはわかったでしょうか？
t検定で使ったデータであるgripを今回のデータであるgrip_anovaに変更しただけです。
それ以外のコードの中身は全く一緒です！
11.今回のコード
今回使ったコードのまとめになります。
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-anova.xlsx"
destfile = "data/demo-anova.xlsx"
download.file(url, destfile)
library(readxl)
grip_anova <- read_excel("data/demo-anova.xlsx",
range = "B2:C62")
View(grip_anova)
#正規分布の確認
library(tidyverse)
ggplot(data = grip_anova) +
geom_histogram(aes(x = grip, fill = category), bins = 5) +
facet_grid(category ~ .)
grip_anova %>%
split(.$category) %>%
map(~shapiro.test(.$grip))
#Bartlett検定
bartlett.test(grip_anova$grip ~ grip_anova$category)
bartlett.test(grip ~ category, data = grip_anova)
#一元配置分散分析
#aov関数の場合
model_grip_anova <- aov(grip ~ category, data = grip_anova)
summary(model_grip_anova)
#lm関数の場合
model_grip_lm <- lm(grip ~ category, data = grip_anova)
summary.aov(model_grip_lm)
#多重比較
#Tukey
TukeyHSD(model_grip_lm)
#Bonferroni
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "bonferroni")
#Holm
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "holm")
#Hochberg
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "hochberg")
#データの要約
grip_anova_summary <-
grip_anova %>%
group_by(category) %>%
summarize(平均 = mean(grip),
標準偏差 = sd(grip),
'0%' = quantile(grip, 0),
'25%' = quantile(grip, 0.25),
'50%' = quantile(grip, 0.5),
'75%' = quantile(grip, 0.75),
'100%' = quantile(grip, 1),
n = n())
grip_anova_summary
#グラフ化
ggplot(data = grip_anova_summary)+
geom_bar(aes(x = category, y = 平均), color = "black", fill = "gray", stat = "identity") +
geom_errorbar(aes(x = category, ymin = 平均 - 標準偏差, ymax = 平均 + 標準偏差), width = 0.1) +
labs(y = "grip") +
theme_classic(base_size = 15) +
scale_y_continuous(expand = c(0,0), limits = c(0,50))
全記事の目次はこちらをご参照ください
