(2019年11月24日修正:8の②のfactor関数を修正しています)

第4章は統計を扱います。


今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。


シロート統計学はEZRを使った統計分析をわかりやすく解説されています。


 


第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか?といった視点で進めていきます。


今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください!


今回は繰り返しのある一元配置分散分析を紹介します



まず繰り返しのある一元配置分散分析(ANOVA)についてはハルさんのサイトをご参照ください。



また1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。
もし1〜4まででわからない部分があれば確認してみてください。

 

1.準備

第4章は毎回ExcelデータをダウンロードしてRを使うのでプロジェクトで管理して行うことを勧めています。

 

ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。
これを前提に次から進めていきます。
スクリーンショット 2019-10-20 7.54.14


2.スクリプトファイルの作成

次にRのコードを書くためのスクリプトファイルを作ります。

スクリーンショット 2019-11-23 22.05.34



3.データのダウンロード


【1-13】Rで読み込みやすいExcelファイルの作り方今回もハルさんのサイトのデータを使わせていただきます。
実は分散分析を行う上で覚えておくと便利なExcelデータの作り方(列の並べ方)があります。
詳細は【1-13】Rで読み込みやすいExcelファイルの作り方で紹介しています。

デモデータ(反復測定分散分析)

この章ではRを使ってダウンロードしています。


download.file(url = “ファイルのURL”,
        destfile = “保存したい場所/ファイル名”)
urlはデモデータで右クリック → リンクのアドレスをコピー
destfileは保存場所と保存のファイル名を指定します。


実際のコードは以下になります。
前回のコードのURL(" "の中)とdestfileのdata/以降を変更するだけでOKです。
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-repeated-measures-anova.xlsx"
destfile = "data/demo-repeated-measures-anova.xlsx"
download.file(url, destfile)
スクリーンショット 2019-11-23 22.14.17


dataフォルダにダウンロードできたことを確認します。


4.データの読み込み

データを読み込みます。
今回は【4-0】第4章を進めていく上での準備で行った方法で進めます。

View Fileでデータを確認します。

スクリーンショット 2019-11-23 22.15.57

データが入っているセルを確認します。
B2からC22までデータが入っています(B2:E32と表記)
スクリーンショット 2019-11-23 22.20.35


次にImport Datasetでデータを取り込みます。
スクリーンショット 2019-11-23 22.15.57


Import画面ではName, Sheet,Rangeを指定します。

Name:ハルさんのサイトではgripでしたがgrip_rep_anovaとします(大文字・小文字は別物とされます)
Sheet:このExcelは1つしかデータがないのでDefaultのままでOK
Range:先ほど確認したB2:E32(実は今回のように表以外に余計なものがなければ空欄でもOK)

スクリーンショット 2019-11-23 22.24.40

Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。
library(readxl)
grip_rep_anova <- read_excel("data/demo-repeated-measures-anova.xlsx", 
                             range = "B2:E32")
View(grip_rep_anova)

データが正しく入っていることを確認します。
スクリーンショット 2019-11-23 22.25.10


これでデータの取り込みは完成です。

5.被験者間要因と被験者内要因と分散分析のデザインについて
分散分析を行う上で大切なのが被験者間要因と被験者内要因についてです。
スクリーンショット 2019-11-24 6.33.40
被験者間要因はいわゆる対応のないデータのことで被験者全体をA法、B法の2つに割り振るといった方法です。なので同じ人がA法とB法を行うことはありません。
被験者内要因はいわゆる対応のあるデータのことで同じ人が繰り返し測定します

そして被験者間要因は列の左側に、被験者内要因は列の右側に並べると後々分析しやすくなります。
また被験者間要因は縦に、被験者内要因は横につなげるとEZRで分析を行うには都合がいいです。

仮に下のデータがあったとします。
スクリーンショット 2019-11-24 6.35.12
この並びが分散分析を行う上での基本的な形になります。
要因が3つ(sex, treatment, 繰り返し)あるのでABC
被験者間要因と被験者内要因の間にsを付ける
上記の場合はABsCデザインと言います。

今回はsexの列を使いませんので繰り返しの1要因のみで分散分析を実施します。
(sexの列を使った分析は次の記事で紹介します)

被験者間要因がありませんので、今回はsAデザインということになります。
繰り返しの群間に差があるかどうかを見たいということです。


6.正規分布かどうかを確認
分散分析を行うにはデータの分布が正規分布であることが必要です。
ハルさんのサイトに沿って正規分布の確認を行います。

正規分布の確認は【4-10】Rで分散分析(一元配置分散分析)を行う方法と基本的には同じですが、今回は繰り返しの部分がwideデータになっています。longデータであれば今までの記事で紹介した方法でまとめて出すことができますが、今回は1つずつ確認します。

今回ヒストグラムはgeom_histgram関数を使いますのでtidyverseパッケージを読み込みます。
EZRで正規性の確認をする時に併せて作成されるヒストグラムは棒が5本なのでbins = 5とすることにします。

library(tidyverse)
ggplot(data = grip_rep_anova) +
  geom_histogram(aes(x = W0), bins = 5) 

ggplot(data = grip_rep_anova) +
  geom_histogram(aes(x = W1), bins = 5) 

ggplot(data = grip_rep_anova) +
  geom_histogram(aes(x = W2), bins = 5) 
スクリーンショット 2019-11-23 22.42.57


7.シャピロ・ウィルク検定
次にハルさんのサイトでは正規性の検定を行っています。
正規性の検討はshapriro.test関数3つのカテゴリーそれぞれに行う必要があります。
ここでもwideデータなので1つずつ行います。

shapiro.test((grip_rep_anova$W0))
shapiro.test((grip_rep_anova$W1))
shapiro.test((grip_rep_anova$W2))
スクリーンショット 2019-11-23 22.44.56


すべて正規分布であることが棄却されませんでした。


8.反復測定分散分析(repeated measures ANOVA)を行う

そして反復測定分散分析を行います。

反復測定分散分析を行うには主に2つの方法があります。

1.carパッケージのAnova関数を使う(EZRもcarパッケージを使ってる)
2.ANOVA君を使う

ANOVA君は井関龍太さんが作成した分散分析に特化したスクリプトです。
使い方に特徴がありますが使えるようになると慣れると簡単でどの方法よりも詳細に分析できます。

ここではまず1の方法を紹介します。

反復測定分散分析を行うには次の方法になります。




①carパッケージを読み込む

まずはcarパッケージを読み込みます。
carパッケージはEZRをインストールする時に一緒にインストールされています。
もしインストールしていなければ先にインストールします。

install.packages("car")

その後library関数で呼び出します。
library(car)

②被験者要因を指定するdata.frameを準備する

今回のデータではsexの列は「MとFの2つのデータがある」とパソコンは読み取れますがW0,W1,W2はwideデータとなっているため「この3列が繰り返しのデータである」ことがパソコンは読み取れません。それがわかるためのものを2つ作ります。

まず繰り返しの順番がわかるtimeというベクトルを作ります。
このデータはfactor型である必要があるためfactor関数を使います。
そしてtimeのデータをdata.frame型に変更したidataを作ります。

 スクリーンショット 2019-11-24 10.11.33
そして大事なのがlevels = です。levels = を指定しないとfacorは五十音順になります。
今回は問題ないのですがもしpre,mid,postのように列名を付けるとlevels = を使わないとmid→post→preの順番になってしまいます。
factor関数に関しては【1-10】Rでよく使われる型について説明します。
 

③lm関数でモデルを作る

次に分析するためのモデルを作ります。
モデルはlm関数を使います。model(名前は何でもいい)という変数名に結果を入れます。
被験者内要因は3列になっているのでcbind関数でつなげます。
被験者間要因は今回無いので1とします。もしあると別の設定も必要ですが次の記事で紹介します。

スクリーンショット 2019-11-24 11.55.22
④Anova関数で分散分析を実行する



反復測定分散分析を行う時はcarパッケージのAnova関数を使います。Aは必ず大文字です。
res(resultの略、何でも良い)という変数名に結果を入れます。

Anova(lm関数で作ったモデル, idata = ◯◯, idesign = ~ ◯◯, type = "Ⅲ")

lm関数で作ったモデル:先程作ったmodel
idata:先程つくったidata
idesign:先程作ったtime(timeの前に を付けます)
type:"Ⅲ"または3とする(反復測定分散分析の場合は必ずⅢ)
スクリーンショット 2019-11-24 12.08.02

⑤summary関数で結果を表示する

summary(Anova関数の結果, multivariate = FALSE)

multivariate = FALSEと入れるとEZRと同じ結果が表示されます。

①〜⑤をまとめると以下になります
library(car)
time<- factor(c("W0", "W1", "W2"), levels = c("W0","W1", "W2")) idata <- as.data.frame(time) model <- lm(cbind(W0,W1,W2) ~ 1, data = grip_rep_anova) res <- Anova(model, idata = idata, idesign = ~time, type = 3 ) summary(res, multivariate = FALSE)
スクリーンショット 2019-11-24 21.02.57
行が長くて大変ですが実は色の部分以外はコピペ可能です。
(被験者間要因があればまた変わります。次回記事で紹介します)

スクリーンショット 2019-11-24 12.51.33
このあたりはハルさんのサイトと同じです。
EZRと同じ結果になります。

9.多重比較を行う

EZRの多重比較はEZR独自の関数を使っています。
そのため今回はpairwise.t.test関数を使います。
(被験者間要因があると使えないので注意!)
ただpairwise.t.test関数はlongデータで行う必要があります。

wideデータをlongデータに変えるにはgather関数、またはpivot_longer関数を使います。
wideデータとlongデータの切り替えは【2-5】Rでデータを集計するのに便利なlongデータとgather関数で紹介しています。

pairwise.t.test(数値, グループ, p.adjust.method = "◯◯")
〇〇はここでは3つの方法を紹介します。

bonferroni
holm
hochberg


grip_rep_anova_long <- 
  grip_rep_anova %>% 
  gather(key = "key", value = "value", W0:W2)

pairwise.t.test(grip_rep_anova_long$value, grip_rep_anova_long$key, paired = TRUE, p.adjust.method = "bonferroni")
pairwise.t.test(grip_rep_anova_long$value, grip_rep_anova_long$key, paired = TRUE, p.adjust.method = "holm")
pairwise.t.test(grip_rep_anova_long$value, grip_rep_anova_long$key, paired = TRUE, p.adjust.method = "hochberg")

スクリーンショット 2019-11-24 13.17.18


10.まとめ

今回は反復測定分散分析を紹介しました。
今までの手法に比べると手続きが必要ですが、慣れるとコピペでもできるようになります。
次回は被験者間要因も含めた混合計画反復測定分散分析を紹介します。


11.今回使ったコード

今回使ったコードは以下になります。
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-repeated-measures-anova.xlsx"
destfile <- "data/demo-repeated-measures-anova.xlsx"

download.file(url, destfile)

#データの読み込み
library(readxl)
grip_rep_anova <- read_excel("data/demo-repeated-measures-anova.xlsx", 
                             range = "B2:E32")
View(grip_rep_anova)

#正規分布の確認
library(tidyverse)
ggplot(data = grip_rep_anova) +
  geom_histogram(aes(x = W0), bins = 5) 

ggplot(data = grip_rep_anova) +
  geom_histogram(aes(x = W1), bins = 5) 

ggplot(data = grip_rep_anova) +
  geom_histogram(aes(x = W2), bins = 5) 

shapiro.test((grip_rep_anova$W0))
shapiro.test((grip_rep_anova$W1))
shapiro.test((grip_rep_anova$W2))


#carパッケージのインストールしたことなければインストール
install.packages("car")

#反復測定分散分析を行う
library(car)
time<- factor(c("W0", "W1", "W2"), levels = c("W0","W1", "W2"))
idata <- as.data.frame(time)
model <- lm(cbind(W0,W1,W2) ~ 1, data = grip_rep_anova)
res <- Anova(model, idata = idata, idesign = ~time, type = 3)
summary(res, multivariate = FALSE)

#多重比較
grip_rep_anova_long <- 
  grip_rep_anova %>% 
  gather(key = "key", value = "value", W0:W2)

pairwise.t.test(grip_rep_anova_long$value, grip_rep_anova_long$key, paired = TRUE, p.adjust.method = "bonferroni")
pairwise.t.test(grip_rep_anova_long$value, grip_rep_anova_long$key, paired = TRUE, p.adjust.method = "holm")
pairwise.t.test(grip_rep_anova_long$value, grip_rep_anova_long$key, paired = TRUE, p.adjust.method = "hochberg")



RStudioではまずExcelやcsvファイル等を読み込む必要がありますが、そもそもExcelではRに読み込みやすい形とそうでない形があります。Excelの見栄え的には良くても実際にRで読み込もうとすると大変な手間になることもあります。

今回はRで読み込みやすいExcelのポイントを紹介します。



1.上や左に空白はないか?

スクリーンショット 2019-10-20 20.07.06

RでなにもせずにExcelファイルを読み込むとA1を始めとしてデータを読み込みます。よくExcelでかけ線を設定している場合B2から表をつくることがありますが、Rで読み込むときには一手間必要です。



2.一番下に合計の行がないか?

加えてもし一番下に合計や平均の列があったらそれも1データとして読み込んでしまうので注意が必要です。
Excelは集計のために使い分析をRで行う場合はそもそも無いほうが安全です。
もしExcelでも集計を行いたい場合は集計を別のタブで行うかピボットテーブルを使う方法があります。



3.タイトル行が2行に渡っていないか?

タイトル行が2行に渡っているとRではうまく読み込んでくれません。
そしてそういう場合はセル結合されていることも多いです。
こういった場合も対処が必要です。

タイトルが2行に渡った時の対処法に関しては以下の記事で説明しています。




4.無駄な空間を作らない
スクリーンショット 2019-11-24 6.59.39
上記の図のように空欄だとパソコンは読み込んでくれません。
列名をちゃんと書き、データは空欄にせず入力しましょう。

5.桁数に使う , を入れていないか

10,000などカンマ含む場合、Rで読み込むと数字ではなく文字として認識されます。
「どうしても , はつけないとダメ」と上司に言われてもRで , を取る方法はあります(parse_number関数)


6.マイナスの値を▲にしていないか?

-100を▲100にしていたらマイナスに直す必要があります。
文字の置換(ここでは▲ → − )にするにはstr_replaceなどの関数がありますが、手間を考えると−の方が手間が省けます。


7.和暦か西暦か?

和暦だと古いExcelではいつまでたっても令和になりません。
計算や列の順番を考えると西暦のほうが後々対処しやすいです。
加えて2001年を01年としないよう注意しましょう。


8.tidyなデータか?

データ分析を行う上で必要な考え方にtidy(タイディー:整然データ)なデータかどうかという考え方があります。

tidyに関しては以下の記事が具体例も交えながらわかりやすく説明されています。




これからデータ収集を考えている方は参考になると思います。


9.被験者間要因は左側に、被験者内要因は右側に並べる
先程tidyなデータの話がありましたが、実際にデータを取りExcelに打ち込む時に覚えておくと後で便利になる考え方があります。それは被験者間要因被験者内要因です。
スクリーンショット 2019-11-24 6.33.40
被験者間要因はいわゆる対応のないデータのことで被験者全体をA法、B法の2つに割り振るといった方法です。なので同じ人がA法とB法を行うことはありません。

被験者内要因はいわゆる対応のあるデータのことで同じ人が繰り返し測定します

加えて被験者間要因は列の左側に、被験者内要因は列の右側に並べると後々分析しやすくなります。
また被験者間要因は縦に、被験者内要因は横につなげるとEZRで分析を行うには都合がいいです。
もし横に列が長くなりすぎて困る場合は全部縦にしても大丈夫です。あとで横に変換できます。
スクリーンショット 2019-11-24 6.35.12
この方法は特に分散分析を行う場合に便利になります。



10.まとめ

予めExcelのデータをRで読み込みやすい形にするとデータの前処理で心が折れにくくなりますので参考にしてみてください!


第4章は統計を扱います。


今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。


シロート統計学はEZRを使った統計分析をわかりやすく解説されています。


 


第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか?といった視点で進めていきます。


今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください!


今回はKruskal-Wallis(クラスカル・ウォリス)検定を紹介します



まずKruskal-Wallis(クラスカル・ウォリス)検定についてはハルさんのサイトをご参照ください。



また1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。
もし1〜4まででわからない部分があれば確認してみてください。

 

1.準備

第4章は毎回ExcelデータをダウンロードしてRを使うのでプロジェクトで管理して行うことを勧めています。

 

ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。
これを前提に次から進めていきます。
スクリーンショット 2019-10-20 7.54.14


2.スクリプトファイルの作成

次にRのコードを書くためのスクリプトファイルを作ります。


スクリーンショット 2019-11-12 7.48.34


3.データのダウンロード

今回もハルさんのサイトのデータを使わせていただきます。

デモデータ(Kraskal-Wallis検定)

この章ではRを使ってダウンロードしています。


download.file(url = “ファイルのURL”,
        destfile = “保存したい場所/ファイル名”)
urlはデモデータで右クリック → リンクのアドレスをコピー
destfileは保存場所と保存のファイル名を指定します。


実際のコードは以下になります。
前回のコードのURL(" "の中)とdestfileのdata/以降を変更するだけでOKです。
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-kruskal-wallis-test.xlsx"
destfile = "data/demo-kruskal-wallis-test.xlsx"

download.file(url, destfile)
スクリーンショット 2019-11-12 20.43.14


dataフォルダにダウンロードできたことを確認します。


4.データの読み込み

データを読み込みます。
今回は【4-0】第4章を進めていく上での準備で行った方法で進めます。

View Fileでデータを確認します。

スクリーンショット 2019-11-12 7.50.38




データが入っているセルを確認します。
B2からC92までデータが入っています(B2:C92と表記)
スクリーンショット 2019-11-12 7.51.15



次にImport Datasetでデータを取り込みます。
スクリーンショット 2019-11-12 7.50.38


Import画面ではName, Sheet,Rangeを指定します。

Name:ハルさんのサイトではBWでしたがbw_kwとします(大文字・小文字は別物とされます)
Sheet:このExcelは1つしかデータがないのでDefaultのままでOK
Range:先ほど確認したB2:C92

スクリーンショット 2019-11-12 20.43.55


Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。
#データの読み込み
library(readxl)
bw_kw <- read_excel("data/demo-kruskal-wallis-test.xlsx", 
                    range = "B2:C92")
View(bw_kw) 

データが正しく入っていることを確認します。
スクリーンショット 2019-11-12 20.46.07


これでデータの取り込みは完成です。

5.正規分布の確認

ハルさんのサイトに沿って正規分布の確認を行います。

正規分布の確認は【4-1】Rでt検定を行う方法で紹介しています。
以下のコードがイメージできない場合はまずこちらの記事をご参照ください。

 

まずヒストグラムで分布を確認します。
1つずつグラフを作るのは手間がかかるので、まとめてグラフを作ります。

ヒストグラムはgeom_histgram関数を使います。
グラフが重なっていると見にくいのでfacet.gridでカテゴリーごとにグラフを分けます。
ちなみに【4-10】Rで分散分析(一元配置分散分析)を行う方法で作ったヒストグラムのコピペでできてしまいます。変数名を変えるだけです。




library(tidyverse)
ggplot(data = bw_kw) +
  geom_histogram(aes(x = BW, fill = category), bins = 5) +
  facet_grid(category ~ .)
スクリーンショット 2019-11-08 15.27.34

6.シャピロ・ウィルク検定
次にハルさんのサイトでは正規性の検定を行っています。
正規性の検討はshapriro.test関数3つのカテゴリーそれぞれに行う必要があります。

今回はsplit関数map関数を使いまとめて行ってみます。
これも【4-10】Rで分散分析(一元配置分散分析)を行う方法で作ったコードと全く同じです。
変数名を変えるだけです。

コードを使うと他でも役に立つことがわかります。

bw_kw %>% 
  split(.$category) %>% 
  map(~shapiro.test(.$BW))
スクリーンショット 2019-11-12 20.54.55


まとめて結果が出ました。全てp < 0.05ですので正規分布であるという仮説が棄却されました。

7.グラフを作る
今回はノンパラメトリックの検定を行うので箱ひげ図を作成します。
箱ひげ図の作り方は【3-6】Rのggplot2で箱ひげ図を作るgeom_boxplot関数で紹介しています。

ggplot(data = bw_kw) +
  geom_boxplot(aes(x = category, y = BW))
スクリーンショット 2019-11-15 8.01.00


もしEZRと同じようなグラフであればboxplot関数も使えるので併せて消化しします。

boxplot(数値 ~ グループ, data = ◯◯)
boxplot(数値 ~ グループ)


boxplot(BW ~ category, data = bw_kw)
スクリーンショット 2019-11-15 8.01.06



8.Kruskal-Wallis(クラスカル・ウォリス)検定を行う

正規性が棄却されましたので今回はノンパラメトリックのKruskal-Wallis(クラスカル・ウォリス)検定を行います。


Kruskal-Wallis検定はEZRではkruskal.test関数を使います。

kruskal.test(目的変数 ~ factor(グループ), data = データ)

目的変数:BW
グループ:factor(category)
data = bw_kw

kruskal.test(BW ~ factor(category), data=bw_kw)
スクリーンショット 2019-11-13 16.53.51

EZRと同じ関数を使っているので同じ結果になりました。


またMann-Whitney U 検定のときに紹介したcoinパッケージでも行うことができます。
もしcoinパッケージを使ったことがなければインストールします。

#coinパッケージのインストール
#1度でもしていればしなくてOK
install.packages("coin")

coinパッケージのkruskal_test関数はkruskal.test関数と同じ使い方です。
使う前にlibrary関数で呼び出します。
そしてグループにはfactor関数もしくはas.factor関数でfactor型に変更します。

kruskal_test(目的変数 ~ factor(グループ), data = データ)

library(coin)

kruskal_test(BW ~ factor(category), data=bw_kw)
スクリーンショット 2019-11-13 22.38.20

Asymptotic(漸近的:正規近似)とあるので正確検定ではないようです。
そのためEZRと同じ結果になりました。
実は3以上のグループではなく2変数だと正確検定を行うことができます。



9.多重比較(Steel-Dwass法)を行う

EZRで使っていたSteel-Dwass法はEZR専用のパッケージなので、Rで行うにはNSM3パッケージpSDCFlig関数を使う必要があります。
ちなみにEZRのSteel-Dwass法は正確検定ではなく正規近似を行います。
NSM3パッケージではデータ数が少なければ正確近似ができますが、データ数が大きくなると計算量が膨大になりすぎて実施できません。モンテカルロ法もしくは正規近似を行うことになります。今回の90例では正確近似はできませんでした。
正規近似と正確検定、Steel-Dwass法の使い方に関しては多重比較 Steel-Dwass 正規近似と正確検定: U 検定が基準に詳しく紹介されています。
#EZRと同じ正規近似で行う場合
pSDCFlig(目的変数, factor(グループ), method = "Asymptotic")

#正確検定を行いたい場合(データが多いと自動的にモンテカルロ法になる)
pSDCFlig(目的変数, factor(グループ))

目的変数:bw_kw$BW
グループ:factor(bw_kw$category)
*グループはfactor型にする必要があるためfactor関数as.factor関数を使います。
*今回もdata = ◯◯ は使えないので を使います。

#Steel-Dwass法を行うためにNSM3パッケージをインストール
#1度でもインストールしてたらこの1行は必要なし
install.packages("NSM3")

#NSM3パッケージの呼び出し
library(NSM3)

#EZRと同じ正規近似で行う場合
pSDCFlig(bw_kw$BW, factor(bw_kw$category), method = "Asymptotic")

#正確検定を行いたい場合(データが多いと自動的にモンテカルロ法になる)
#結果に時間がかかるので注意! pSDCFlig(bw_kw$BW, factor(bw_kw$category))
スクリーンショット 2019-11-15 2.39.39

それぞれのp値は以下になりました。
大きな差はありません。
A-B:正規近似0.014、モンテカルロ0.0142
A-C:正規近似0.0255、モンテカルロ0.0249
B-C:正規近似0.9006、モンテカルロ0.902


10.ひとまずまとめ

今回はRでKruskal-Wallis(クラスカル・ウォリス)検定と多重比較であるSteel-Dwass法の紹介をしました。ハルさんのサイトではここまでですがBonferroniやHolm法でp値の調整をするにはどうしたら良いのか?Mann-Whitney U検定でさんざん問題になったタイと正確検定に関してはどうなったんだ?という方は続きがありますので興味がある方は続きをご覧ください。

次回は反復測定分散分析を紹介します。



11.他の多重比較の方法

一元配置分散分析でも紹介しましたが、他の多重比較にはBonferroniHolmHochbergなどがあります。



一元配置分散分析の多重比較ではpairwise.t.test関数を使いましたが、Kruskal-Wallisではpairwise.wilcox.test関数を使います。簡単に言うと全ての組み合わせでMann-Whitney U検定を行いその結果を元にp値の調整を行います
どの方法を使うかはp.adjust.method = で指定します。

pairwise.wilcox.test(目的変数factor(グループ),  p.adjust.method = "◯◯")

pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "bonferroni")
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "holm")
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "hochberg")
スクリーンショット 2019-11-15 19.11.42

ただ、この方法だと警告が出ます。タイ(同順位)があるためです。
タイに関してはMann-Whitney U検定で紹介しました。




今回はずっと正規近似で行っているので正確近似で出すのであれば exact = FASEを追記します。
そうすればエラーは出ません。

pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "bonferroni", exact = FALSE)
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "holm", exact = FALSE)
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "hochberg", exact = FALSE)
スクリーンショット 2019-11-15 19.21.16


11.どうしても正確検定を行いたい!!!という方は手計算が必要

多重比較は各組み合わせでMann-Whitney U検定を行いその結果を元にp値の調整を行います。
ただpairwise_wilcox_test関数が存在せず自分で手計算するしかありません。

そのためどうしても正確検定を行いたい場合はcoinパッケージwilcox_test関数を使い、1つずつ組み合わせを行い、BonferroniHolmHochbergそれぞれの方法で自分でp値を直接計算することになります。

この方法は統計ERさんの記事に詳しく紹介されています。
更にBonferroniHolmHochbergって何をしてるの?という疑問も解決するかもしれません。
仕組みがわかればp値に2とか3とかかけるだけなので意外と単純です。



12.本当にまとめ
今回はRでKruskal-Wallis(クラスカル・ウォリス)検定と多重比較であるSteel-Dwass法の紹介をしました。ノンパラメトリック検定はタイと正確検定をするか正規近似を行うか?の問題がありますが、Mann-Whitney U検定のちゅんたろさんのツイートがここでも参考になると思います。


そしてそもそもp値がより小さければより差が大きいわけではありませんので、正確検定か正規近似か?よりも研究デザインが適切かどうかの方が大切になるかもと言うのが個人的なイメージです。

次回は反復測定分散分析を紹介します!

13.今回使ったコード

今回使ったコードを紹介します。
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-kruskal-wallis-test.xlsx"
destfile = "data/demo-kruskal-wallis-test.xlsx"

download.file(url, destfile)

#データの読み込み
library(readxl)
bw_kw <- read_excel("data/demo-kruskal-wallis-test.xlsx", 
                    range = "B2:C92")
View(bw_kw) 

#正規性の確認
library(tidyverse)
ggplot(data = bw_kw) +
  geom_histogram(aes(x = BW, fill = category), bins = 5) +
  facet_grid(category ~ .)


bw_kw %>% 
  split(.$category) %>% 
  map(~shapiro.test(.$BW))

#グラフの作成
ggplot(data = bw_kw) +
  geom_boxplot(aes(x = category, y = BW))

boxplot(BW ~ category, data = bw_kw)

#Kruskal-Wallis検定
kruskal.test(BW ~ factor(category), data = bw_kw)

#coinパッケージでKruskal-Wallis検定を行う場合
library(coin)
kruskal_test(BW ~ factor(category), data = bw_kw)


#Steel-Dwass法を行うためにNSM3パッケージをインストール
#1度でもインストールしてたらこの1行は必要なし
install.packages("NSM3")

#NSM3パッケージの呼び出し
library(NSM3)

#EZRと同じ正規近似で行う場合
pSDCFlig(bw_kw$BW, factor(bw_kw$category), method = "Asymptotic")

#正確検定を行いたい場合(データが多いと自動的にモンテカルロ法になる)
#結果に時間がかかるので注意!
pSDCFlig(bw_kw$BW, factor(bw_kw$category))


#他の多重比較

pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "bonferroni")
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "holm")
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "hochberg")

#正規近似を使えば警告は出ない。exact = FALSEを追加すれば良い
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "bonferroni", exact = FALSE)
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "holm", exact = FALSE)
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "hochberg", exact = FALSE)


全記事の目次はこちらをご参照ください

 


↑このページのトップヘ