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