タグ:boxplot

第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)


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

 


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


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


ハルさん、ありがとうございます!


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


 


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


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


今回はMann-Whitney U 検定を紹介します



まずMann-Whitney U 検定についてはハルさんのサイトをご参照ください。




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



1.準備

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



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


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

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

スクリーンショット 2019-10-29 19.25.29


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

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

スクリーンショット 2019-11-01 1.21.02


デモデータ(Mann-Whitney U検定)

前回Rを使ってダウンロードしましたが、まだ慣れないと思いますので再度説明します。


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


実際のコードは以下になります。
前回のコードのURL(" "の中)とdestfileのdata/以降を変更するだけでOKです。
url <- "https://haru-reha.com/wp-content/uploads/2018/04/demo-mann-whitney-u-test.xlsx"
destfile = "data/demo-mann-whitney-u-test.xlsx"

download.file(url, destfile)
スクリーンショット 2019-10-29 19.43.04


dataフォルダにダウンロードできました!
スクリーンショット 2019-10-29 19.36.48


4.データの読み込み

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

View Fileでデータを確認します。
スクリーンショット 2019-10-29 19.46.19

データが入っているセルを確認します。
B2からC62までデータが入っています(B2:C62と表記)

スクリーンショット 2019-10-29 19.47.48

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


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

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




Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。
library(readxl)
mmt <- read_excel("data/demo-mann-whitney-u-test.xlsx", 
                  range = "C2:D62")
View(mmt)
スクリーンショット 2019-10-29 21.03.46

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

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


5.Mann-Whitney U 検定を行う

早速Mann-Whitney U 検定を行ってみます。

Mann-Whitney U 検定(マン・ホイットニーのU検定)は別な言い方もあります。

Wilcoxon-Mann-Whitney検定
Wilcoxon の順位和検定


EZRで使われているMann-Whitney U 検定はwilcox.test関数を使っています。

wilcox.test(MMT ~ category, data = mmt)

スクリーンショット 2019-11-01 1.02.29
スクリーンショット 2019-11-01 0.40.49


EZRの結果はこうです。
スクリーンショット 2019-11-01 1.07.20


EZRと同じ結果がでました。警告が気になりますがひとまず置きます。

もしp値だけがほしいときは変数名を付けて変数名$p.valueとします。
mmt_wilcox <- wilcox.test(MMT ~ category, data = mmt)
mmt_wilcox$p.value


6.グラフの作成
EZRでは同時に箱ひげ図を作成しますがRでは自作する必要があります。

箱ひげ図の作成に関しては第3章で紹介しています。



ggplot(data = mmt) +
  geom_boxplot(aes(x = category, y = MMT))
スクリーンショット 2019-11-01 1.16.42

もしEZRと同じグラフにしたい場合はboxplot関数を使うという選択肢もあります。
boxplot(MMT ~ category, data = mmt)
スクリーンショット 2019-11-01 1.19.01



7.一旦まとめ

RでMann-Whitney U 検定を行う方法を紹介しました。
wlcox.testを使えばEZRと同じ結果になります。

ただ警告が出てしまい心配な方はいると思います。

以下からは少し応用編となりますがこの警告について進めていきます。




8.タイがあるため、正確な p 値を計算することができません

スクリーンショット 2019-11-01 0.40.49

Rでノンパラメトリックの検定を行っているとこういった画面が出ることがあります。

cannot compute exact p-value with ties
タイがあるため、正確な p 値を計算することができません
正確なp値が計算できないと言われるととても心配になります。

実は先程のEZRでも警告は出ていました。
一番下の行です。
スクリーンショット 2019-11-01 1.35.45

cannot compute exact p-value with ties


8.タイ(ties)って何?

タイというのは同順(同じ数値)のことです。

タイがないデータ:1,2,4,5,7,10,11,12
タイがあるデータ:1,2,2,4,6,6,6,7,9,11


詳細は成書に譲りますがMann-Whitney U 検定は平均ではなくデータを少ない順に並べてデータに順位を付け計算するため、同順位があるとそれをどう計算するかといった問題や様々な対処法が出てきます。


9.タイに対してどう考えるか?

正確なp値を求めるには膨大な計算が必要になるらしく、そのため正規近似というものを使って計算するようです。EZR(wilcox.test)はタイがなければ正確なp値を求めるオプションがありますが、タイがあるとオプションを付けても正確なp値を求めることができないようです。近似値ではなく正確なp値を求めるにはcoinパッケージwilcox_test関数を使う必要があります。coinパッケージに関しては奥村先生の記事や裏 RjpWikiの記事に詳しく説明があります。






10.coinパッケージを使ってみる

wilcox_test関数を使うためにはcoinパッケージをインストールする必要があります。
パッケージのインストールは今までinstall.packages関数を紹介していましたが、今回はRStudioのPackagesタブからインストールしてみます。
PackagesタブのInstallをクリックしてパッケージ名を入れます。

スクリーンショット 2019-10-29 21.37.10

これでインストールは完了です。
これでlibrary関数を使ってcoinパッケージを呼び出せます。
wilcox.testとの違いは主に3点あります。
①グループの箇所は必ずfactor型に変換(as.factor関数またはfactor関数を使う)
(わからないときは【4-1】Rでt検定を行う方法を参照)
②正確なp値を求めるためにはdistribution = "exact"を追加
③実はwilcox.testはlongデータでもwideデータデータでも使えるが、wilcox_testはlongデータしか扱えない
(わからないときは【1-10】Rでよく使われる型について説明しますを参照)
スクリーンショット 2019-11-02 4.52.27

library(coin)
mmt_coin <- wilcox_test(mmt$MMT ~ as.factor(mmt$category), distribution="exact") mmt_coin
スクリーンショット 2019-11-02 4.35.30


もしp値だけを出したいときはpvalue関数、95%信頼区間を出したいときはpvalue_interval関数を使います。
pvalue(mmt_coin)
pvalue_interval(mmt_coin)
スクリーンショット 2019-11-02 4.35.45


これで正確なp値を求めることができました。


11.じゃあEZRはダメなのか?

ただ必ず正確なp値を求めないといけないかというとそうではないようです。

ちゅんたろさんからアドバイスを頂きました。







ある程度n数があればEZRでも利用可能と考えることもできます。
このあたりは統一見解があるわけではなさそうです。

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


12.今度こそまとめ

今回はMann-Whitney U 検定を紹介しました。
実はwilcox.test(EZR)でタイが無いときの正確なp値の求め方などwilcox.testのオプションなど話は割愛させてもらっています。

qiitaの記事ではwilcox.testとwilcox_testのオプションやタイの有無で結果がどう変わるかについての解説をしていますので、興味があればぜひこちらも併せてご覧ください。



13.今回使ったコード

今回使ったコードをまとめて置いておきます。
coinパッケージのインストールに関してもコードにしています。

#Mann-Whitney U 検定を行うためにcoinパッケージのインストール
#1度でもしていればしなくてOK
install.packages("coin") #もしtidyverseパッケージをインストールしていなければインストール。 #1度でもしていればしなくてOK install.packages("tidyverse") #パッケージの読み込み library(coin) library(readxl) library(tidyverse) #データのダウンロード url <- "https://haru-reha.com/wp-content/uploads/2018/04/demo-mann-whitney-u-test.xlsx" destfile = "data/demo-mann-whitney-u-test.xlsx" download.file(url, destfile) #データの読み込み mmt <- read_excel("data/demo-mann-whitney-u-test.xlsx", range = "B2:C62") View(mmt) #wilcox.test wilcox.test(MMT ~ category, data = mmt)
boxplot(MMT ~ category, data = mmt) #wilcox.testでp値を求めたい時(EZRと同じ) mmt_wilcox <- wilcox.test(MMT ~ category, data = mmt) mmt_wilcox$p.value #coinパッケージのwilcox_testを使う mmt_coin <- wilcox_test(mmt$MMT ~ as.factor(mmt$category), distribution="exact") mmt_coin pvalue(mmt_coin) pvalue_interval(mmt_coin)

第3章ではggplot2を使ったグラフの作り方を紹介しています。


【3-1】ExcelにはないRでグラフを作るメリットと特徴

【3-2】ggplot2でグラフを作る流れを説明します

【3-3】Rのggplot2で散布図を作るgeom_point関数

【3-4】Rのggplot2でヒストグラムを作るgeom_histogram関数

【3-5】Rのggplot2で密度曲線を作るgeom_density関数

今回は箱ひげ図(boxplot)を紹介します。

箱ひげ図は正規分布でないデータやマンホイットニーのU検定などノンパラメトリックの検定で使われます。RコマンダーやEZRでも検定と一緒に箱ひげ図が出てきますが、変数名が日本語だと□□と文字化けしてしまいます。Rstudioでggplot関数を使えば文字化け対策やグループ毎に色を付けるなどきれいなグラフが作れます。

1.データの読み込み

今回もggplotパッケージが含まれているtidyverseパッケージを読み込みます。
#tidyverseパッケージをインストールしていなければインストール。していれば次へ

install.packages("tidyverse")

#既にtidyverseパッケージをインストールしている方は以下でもOK

library(tidyverse)

#データ取り込みます。今回はdatという変数にデータを入れます

url <- "https://github.com/mitti1210/myblog/blob/master/heights.csv?raw=true"
dat <- read.csv(url)

このデータはヒストグラムを作る時に作ったデータと同じものを使用しています。



2.箱ひげ図(boxplot)とは?


箱ひげ図は以下のようなグラフです。

スクリーンショット 2019-07-21 23.28.32


棒グラフは「平均値」だけですが、箱ひげ図は「最小値, 第一四分位数, 中央値, 第三四分位数, 最大値」を示します。
また「ひげ」の外にあるデータは「外れ値」として点で表示され、大まかな分布を知ることができます。



3.四分位とは

四分位範囲も含め箱ひげ図で表示されるデータは全て「順位」で決められています。スクリーンショット 2019-07-22 1.00.16

これらはqantile関数で直接求めることができます。



4.箱ひげ図の特徴

(平均ではなく)中央値を比較するのに向いている

箱ひげ図は平均ではなく中央値を表示します。
そのため正規分布でないデータやノンパラメトリックの検定でよく使われます。

スクリーンショット 2019-07-22 1.18.14


簡易的な分布がわかる

棒グラフ±標準偏差ではある程度の分布は示してくれますが、実際今回のデータの最大値・最小値がどうだったのか等は教えてくれません。それに比べ箱ひげ図はある程度の分布を表してくれます。

ただ注意したいのが、データに二峰性の分布があると箱ひげ図でも分布がうまく反映できません。


スクリーンショット 2019-07-22 1.10.06



グラフを使って何を示したいのか?が大切になってきます。


5.箱ひげ図の基本的な作り方

ggplot2でヒストグラムを作る時にはgeom_boxplot関数を使います。
ggplot()+
  theme_gray(base_family = "HiraKakuPro-W3")+
  geom_boxplot(aes(x = 性別, y = 身長), data = dat)
スクリーンショット 2019-07-17 1.12.03

6.線の色を変える

線の色を変える時はcolor = ○○を使います。

ggplot()+
  theme_gray(base_family = "HiraKakuPro-W3")+
  geom_boxplot(aes(x = 性別, y = 身長), color = "red", data = dat)


スクリーンショット 2019-07-17 1.12.10


7.中に色をぬる

中の色はfill = ○○を使います。
ggplot()+
  theme_gray(base_family = "HiraKakuPro-W3")+
  geom_boxplot(aes(x = 性別, y = 身長), fill = "red", data = dat)

スクリーンショット 2019-07-17 1.12.19


色を薄くする時はalpha = ○○を使います。
alphaは0〜1の値を選択します。

ggplot()+
  theme_gray(base_family = "HiraKakuPro-W3")+
  geom_boxplot(aes(x = 性別, y = 身長), fill = "red", alpha = 0.5, data = dat)


スクリーンショット 2019-07-17 1.12.26


8.グループ毎に色を変える

グループ毎に色を指定するにはaes関数の中にcolorやfillを入れます。
ggplot()+
  theme_gray(base_family = "HiraKakuPro-W3")+
  geom_boxplot(aes(x = 性別, y = 身長, fill = 性別), alpha = 0.5, data = dat)


スクリーンショット 2019-07-17 1.12.32

9.応用編:グラフを横にする

グラフの向きを横にするにはcoord_flip()を使います。()の中は何も入れません。

ggplot()+
  theme_gray(base_family = "HiraKakuPro-W3")+
  geom_boxplot(aes(x = 性別, y = 身長, fill = 性別), alpha = 0.5, data = dat)+
  coord_flip()



スクリーンショット 2019-07-17 1.12.38




10.応用編:ひげ(wisker)の求め方


箱ひげ図は「最小値, 第一四分位数, 中央値, 第三四分位数, 最大値」を示します。
ただ上記のグラフのようにひげは外れ値がなければひげの上端・下端が最大・最小値になるのですが、男性のように外れ値があると髭の長さは上記のどれにも当てはまりません。

もしひげの上端(下端)の値が実際何なんなのか?をこちらの記事で紹介しています。



ひげの上端・下端を求める関数を自作しました。もう少し勉強してみたい方はこちらもご参照ください。


11.まとめ

今回は箱ひげ図を紹介しました。

RコマンダーやEZRで納得の行く箱ひげ図にならなかった場合に是非試してみてください。

↑このページのトップヘ