カテゴリ: 統計

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


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


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


 


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


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


今回はFisherの正確検定を紹介します



まずFisherの正確検定についてはハルさんのサイトをご参照ください。



また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-05 11.14.15


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

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

デモデータ(Fisherの正確検定)

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


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


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

download.file(url, destfile)

スクリーンショット 2019-11-05 11.18.20


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




4.データの読み込み

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

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

スクリーンショット 2019-11-05 11.19.37



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



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


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

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

スクリーンショット 2019-11-05 11.24.14


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

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


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


5.データテーブルの作成
データの要約は【2-6】Rでgroup_by関数とsummarize関数を使ってグラフ作成に必要な統計量(平均や標準偏差など)を求めるで紹介しました。
group_by関数とsummarize関数を使って要約するならこうなります。
ExcelのSexの列のSが大文字なので注意が必要です。
#データの要約
library(tidyverse)
sex %>% 
  group_by(category, Sex) %>% 
  summarize(n = n())
スクリーンショット 2019-11-05 11.39.34

これでもいいのですが、EZRのようにtable形式にしたい場合はtable関数を使います。

table(sex$Sex, sex$category)
スクリーンショット 2019-11-05 11.39.46




6.Fisherの正確検定を行う

Fishsrの正確検定を行うにはfisher.test関数を使います。そのままでわかりやすいです。
列名を2つ指定するだけです。

fisher.test(1列目, 2列目)

fisher.test(sex$Sex, sex$category)
スクリーンショット 2019-11-05 11.39.59

EZRで0.205なので今回の結果を四捨五入すると同じ結果です。


95%信頼区間も出ているのでグラフを作ってみました。結果ではtrue odds ratio is not equal to 1、つまりオッズ比が1であるかどうかで判断してるので1で線を引きます。

スクリーンショット 2019-11-05 15.30.46

95%信頼区間が1を挟んでいますのでpは0.05以上と判断できます。
今回の信頼区間はかなり広いことも読み取れます。
95%信頼区間はデータ数が増えると幅が狭くなります。


(追記)
7.χ二乗検定を行うには

χ(カイ)二乗検定を行うにはchisq.test関数を使います。
fisherをchisqに変えるだけで中身は同じです。
chisq.test(sex$Sex, sex$category)


8.まとめ
今回はFisherの正確検定を紹介しました。

【4-1】から進めている方は少しずつ慣れてきたでしょうか。
このサイトはそのため第1章から順に読むと徐々に知識が追加され、途中で復習できるよう構成しています。もしわからない箇所が多ければサイトマップを見ていただければ別の発見があるかもしれません。

次回は検定の結果から(p値や信頼区間)のデータを取り出す方法を紹介します。

 


9.今回使ったコード

今回使ったコードをまとめて置いておきます。
95%信頼区間のコードも置いています。

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/04/demo-fishers-exact-test.xlsx"
destfile = "data/demo-fishers-exact-test.xlsx"

download.file(url, destfile)


library(readxl)
sex <- read_excel("data/demo-fishers-exact-test.xlsx", 
                  range = "B2:C42")
View(sex)

#データの要約
library(tidyverse)
sex %>% 
  group_by(category, Sex) %>% 
  summarize(n = n())

#データテーブルの作成
table(sex$Sex, sex$category)

#fihsrの正確検定
fisher.test(sex$Sex, sex$category)


#グラフの作成
res <- fisher.test(sex$Sex, sex$category)

ggplot()+
  geom_errorbar(aes(x = "", ymin = res$conf.int[1], ymax = res$conf.int[2]), width = 0.1) +
  geom_text(aes(x = "", y = res$conf.int[1], label = round(res$conf.int[1], 2)), vjust = -1) +
  geom_text(aes(x = "", y = res$conf.int[2], label = round(res$conf.int[2], 2)), vjust = -1) +
  geom_point(aes(x = "", y = res$estimate)) +
  geom_text(aes(x = "", y = res$estimate, label = round(res$estimate, 2)), vjust = -1) +
  geom_hline(yintercept = 1, color = "red") +
  labs(x = "", y = "") +
  coord_flip()

#χ二乗検定
chisq.test(sex$Sex, sex$category)



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


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


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


 


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


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


今回はWilcoxon符号付順位和検定を紹介します



まずWilcoxon符号付順位和検定についてはハルさんのサイトをご参照ください。



また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-04 0.47.12



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

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

デモデータ(wilcoxon符号付順位和検定)


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


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


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

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

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




4.データの読み込み

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

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

スクリーンショット 2019-11-04 7.44.03



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


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

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

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

スクリーンショット 2019-11-04 7.41.10

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

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

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


5.データを要約する

次にデータを要約します。

データの要約は【4-1】Rでt検定を行う方法で紹介しました。
group_by関数とsummarize関数を使って要約しましたが、今回はsummary関数を使います。
summary関数は平均、中央値、最大・最小値、四分位範囲をまとめて出してくれますが、標準偏差はだしてくれません。ただ今回はノンパラメトリックなので標準偏差はいらないだろうという理由です。

そして今回はwideデータになっています。
スクリーンショット 2019-11-04 8.05.00
なのでコードはこうなります。
summary(borg$pre)
summary(borg$post)
スクリーンショット 2019-11-04 19.56.42


6.Wilcoxon符号付順位和検定を行う

次にWilcoxon符号付順位和検定を行います。

【4-3】Rで対応のあるt検定を行う方法で紹介しましたが、対応のあるとなしはt.test関数paired = TRUEをつけるかどうかの違いでした。

実はWilcoxon符号付順位和検定も同じです。

wilcox.testにpaired = FALSEをつける(もしくは何も付けない)とMann-Whitney U 検定
wilcox.testにpaired = TRUEをつけるとWilcoxon符号付順位和検定


ということで、wilcoxテストを行ってみます。
t.test関数もそうでしたが、longデータとwideデータで書き方が違います。

wideデータの場合 → , を使う
wilcox.test(1列目, 2列目, paired = TRUE)

longデータの場合 → を使う
wilcox.test(数値 ~ グループ, paired = TRUE)


今回はwideデータなのでこうなります。
wilcox.test(borg$pre, borg$post, paired = TRUE)
スクリーンショット 2019-11-04 20.38.57

EZRと同じ結果になりましたが、Mann-Whitney U 検定のときに悩ませたたアレが出てきました。
今度はもう1行増えてます。
タイがあるため、正確な p 値を計算することができません 
ゼロ値のため、正確な p 値を計算することができません
ちなみにEZRでも警告が出ています。
スクリーンショット 2019-11-04 20.59.27

cannot compute exact p-value with ties
cannot compute exact p-value with zeroes


7.タイのあるデータの対処法

EZRでも使われているwilcox.test関数はタイ(同順位)があると正確なp値を計算できず、近似値を計算する設定になっていました。





今回も同じ問題が出ています。
ブログではある程度のn数があればEZR(wilcox.test)でもいいのではという話がありました。
ただ警告が気持ち悪い!正確なp値も知りたいというための方法も紹介します。

タイに対しては奥村先生の記事が参考になります。



①coinパッケージのwilcoxsign_test関数

Mann-Whitney U 検定ではタイがあっても正確なp値を計算するcoinパッケージwilcox_test関数がありました。coinパッケージを使ってWilcoxon符号付順位和検定を行う場合はwilcoxsign_test関数を使います。

まだcoinパッケージを1度も使ったことがなければインストールします。

coinパッケージのインストール
install.packages("coin")

wilcoxsign_test関数の書き方はちょっとクセがあります。。。

スクリーンショット 2019-11-04 23.11.32
#パッケージの読み込み
library(coin)

wilcoxsign_test(borg$pre ~ borg$post, distribution = "exact", zero.method="Wilcoxon")
スクリーンショット 2019-11-04 23.12.57


exactRankTestsパッケージのwilcox.exact関数

もう1つexactRankTestsパッケージがあります。
このパッケージは開発が終わっており、インストールするとcoinパッケージ使ってねと警告が出ます。それでもcoinパッケージのwilcoxsign_test関数と同じ結果になります。
スクリーンショット 2019-11-05 0.06.33

まずexactRankTestsパッケージをインストールします。
#exactRankTestsパッケージのインストール
install.packages("exactRankTests")
wilcox.exact関数はwilcox.testと似たような書き方ができるのでわかりやすいのが特徴です。
スクリーンショット 2019-11-04 23.53.16
library(exactRankTests)
wilcox.exact(borg$pre, borg$post, paired = TRUE, exact=TRUE)
スクリーンショット 2019-11-04 23.55.13

どちらも同じになりました。


8.まとめ
今回はWilcoxon符号付順位和検定を紹介しました。
Mann-Whitney U 検定との共通点や相違点を比較するとイメージが深まると思います。

4章を順に見ていくと重複する箇所も出てきますので検索で来られた方はサイトマップを見ていただければ別の発見があるかもしれません。

 


9.今回使ったコード

今回使ったコードをまとめて置いておきます。

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/04/demo-wilcoxon-rank-sum-test.xlsx"
destfile = "data/demo-wilcoxon-rank-sum-test.xlsx"

download.file(url, destfile)


library(readxl)
borg <- read_excel("data/demo-wilcoxon-rank-sum-test.xlsx", 
                   range = "B2:C22")
View(borg)

summary(borg$pre)
summary(borg$post)

summary_pre <- summary(borg$pre)
summary_post <- summary(borg$post)


#Wilcoxon符号付順位和検定
wilcox.test(borg$pre, borg$post, paired = TRUE)


#coinパッケージのインストール
install.packages("coin")

library(coin)
wilcoxsign_test(borg$pre ~ borg$post, distribution = "exact", zero.method="Wilcoxon")

#exactRankTestsパッケージのインストール install.packages("exactRankTests") library(exactRankTests) wilcox.exact(borg$pre, borg$post, exact=TRUE, paired = TRUE)

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


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


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


 


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


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


今回は対応のあるt検定を紹介します



まず対応のあるt検定についてはハルさんのサイトをご参照ください。



また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-03 8.11.46


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

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

デモデータ(対応のあるt検定)


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


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


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

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



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



4.データの読み込み

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

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

スクリーンショット 2019-11-03 18.53.56


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

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

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

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

スクリーンショット 2019-11-03 19.01.56


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

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


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


5.正規性の確認

正規性の確認は【4-1】Rでt検定を行う方法でも紹介しました

今回ハルさんのサイトではdifferenceの列の正規性をQQプロットを使って確認しています。

QQプロットの出し方はいくつかありますが、ここでは2つ紹介します。

まずはqqline関数です。
qqline(y = gait$difference)
スクリーンショット 2019-11-03 19.36.51
もしEZRと同じグラフが出したい場合はcarパッケージqqPlot関数を使います(Pが大文字なので注意!)
Packagesにcarパッケージが入っていれなければinstall.packages関数でインストール後library関数で呼び出します。

ちなみに先程のqqline関数y = でしたが、qqPlot関数x = となっています。 

install.packages("car")

library(car) qqPlot(x = gait$difference)
スクリーンショット 2019-11-03 19.43.26

こうすることでEZRを開かなくても同様のことが行なえます。
確認するのが目的であればEZRでいいと思いますが、1年後に再現しようと思ったらスクリプトに残しておくと再現性が高まります。

ヒストグラムはgeom_histgramもしくはhist関数を使います。
library(tidyverse)
ggplot(data = gait)+
  geom_histogram(aes(x = difference), bins = 5)

hist(gait$difference)
スクリーンショット 2019-11-03 19.52.12

シャピロウィルク検定は【4-1】Rでt検定を行う方法でも紹介しましたがshapiro.testでした。
shapiro.test(gait$difference)
スクリーンショット 2019-11-03 19.55.47


6.対応のあるt検定を行う

対応のあるt検定は実は【4-1】Rでt検定を行う方法でも紹介したt.test関数paired = TRUEを足すだけです。ちなみに対応のないt検定はpaired = FALSEをしていることになります。何も書かなければpaired = FALSEになります。

t検定は2種類の書き方があります。その違いはデータの形にあります。
ハルさんのサイトでEZRでは右の形しかできませんでしたがRならどちらも可能です
ちなみに左をlongデータ(縦に長い)、右をwideデータ(横に広い)といいます。
スクリーンショット 2019-11-03 20.03.44

下の図【4-1】Rでt検定を行う方法で紹介しましたが、これにpaired = TRUEを付け加えます。ただlongデータで行う場合はAの並び順とBの並び順が違うと計算結果が異なりますので注意が必要です。
スクリーンショット 2019-10-27 8.09.41

今回はwideデータなので以下のようになります。
var = は外しています。

t.test(gait$pre, gait$post, paired = TRUE)

スクリーンショット 2019-11-03 20.12.08

EZRと同じ結果になりました。

更に95%信頼区間は41.46340 - 86.86993となっており、差の平均(期待値)は64.16667となっています。

どういうことかと言うと、対応のあるt検定は2つの差(ここではpost - pre)が0であると仮定しています。そして95%信頼区間に0が含まれているとpは0.05以上と判断されます。
例えば差の平均(期待値)が20でも95%信頼区間が-30 〜50といった具合です。

スクリーンショット 2019-11-03 20.37.08
上のsampleの場合、差の平均(期待値)は20とプラスなのでpostの方がより高い値と言えそうですが、95%信頼区間に0やマイナスが含まれています。言い換えると20だったけど-10にだってなり得るし、0にだってなり得るとなります。そうなればpostの方が高い値といえません。

ちなみに今回は差(期待値)が0を仮定していたので0でしたが、例えばロジスティック回帰分析だとodd比が1であると仮定するので、そのときは1を挟むかどうかを確認することになります。

p値だけでなく信頼区間を確認するとまた発見があるかもしれません。


7.まとめ

今回は対応のあるt検定を行いました。
実際には対応のないt検定にpaired = TRUEを付け加えるだけでした。

4章を順に見ていくと重複する箇所も出てきますので検索で来られた方はサイトマップを見ていただければ別の発見があるかもしれません。




8.今回使ったコード

今回使ったコードをまとめて置いておきます。
95%信頼区間のコードも載せています。


#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/04/demo-paired-t-test.xlsx"
destfile = "data/demo-paired-t-test.xlsx"

download.file(url, destfile)

#データの読み込み
library(readxl)
gait <- read_excel("data/demo-paired-t-test.xlsx", 
                   range = "B2:D32")
View(gait)

#正規性の検定

#qqplot
qqline(y = gait$difference)

#carパッケージの
install.packages("car")
library(car)
qqPlot(x = gait$difference)

#ヒストグラム
library(tidyverse)
ggplot(data = gait)+
  geom_histogram(aes(x = difference), bins = 5)

hist(gait$difference)

#シャピロ・ウィルク検定
shapiro.test(gait$difference)

#t検定
t.test(gait$pre, gait$post, paired = TRUE)

#信頼区間
gait_ttest <- t.test(gait$pre, gait$post, paired = TRUE)

conf <- data.frame(conf.low = c(-30, round(gait_ttest$conf.int[[1]], 3)), conf.high = c(50,round(gait_ttest$conf.int[[2]],3)), mean = c(20,round(gait_ttest$estimate, 3)), x = c("sample", "post - pre"))
ggplot(data = conf) +
  geom_errorbar(aes(x = x, ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_hline(yintercept = 0, color = "red") +
  geom_text(aes(label = conf.low, x = x, y = conf.low), vjust = -1) +
  geom_text(aes(label = conf.high, x = x, y = conf.high), vjust = -1) +
  geom_text(aes(label = mean, x = x, y = mean), vjust = -1) +
  geom_point(aes(x = x, y = mean)) +
  labs(x = "", y = "")+
  coord_flip()

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

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


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


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


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




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


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


今回はt検定を紹介します



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

 



1.準備

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



ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。
これを前提に次から進めていきます。


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

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

スクリーンショット 2019-10-25 12.15.42

完成です。
スクリーンショット 2019-10-25 12.18.51


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

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

デモデータ(t検定)

これをダウンロードしてdataフォルダに入れればいいのですが実はRでできてしまいます
download.file関数を使います。" "を忘れないようにしてください。

url <- "https://haru-reha.com/wp-content/uploads/2018/03/demo-t-test.xlsx"
destfile = "data/demo-t-test.xlsx"

download.file(url, destfile)

以下説明します。

download.file(url = “ファイルのURL”,
        destfile = “保存したい場所/ファイル名”)


urlはデモデータで右クリック → リンクのアドレスをコピー

destfileは保存場所と保存のファイル名を指定します。
保存場所は今回プロジェクトを使っているのでR練習フォルダになります。加えてdata/を付け足すことでR練習フォルダ内にあるdataフォルダという意味になります。
ファイル名は自由に決められますが今回は元のファイルと同じにしました。拡張子も忘れないようにしましょう。

もしプロジェクトを使っていなければ保存場所はgetwd関数で出てきたフォルダになります。
getwd関数の()には何も入れません。

getwd()


この方法を使う最大のメリットは、次回使うExcelデータはurlの部分を変えるだけでできてしまうことです。毎回右クリックでアドレスを保存 → 保存したファイルを指定したところに移動させて・・・といった手作業必要ありません。こういった作業もRで行っていくことでRにも早く慣れてくると思います。



4.データの読み込み

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

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

今回は握力のデータです。A群とB群を比較します。
ただA1にデータが入っていません。実際にはB2からC62までデータが入っていることを確認します。
スクリーンショット 2019-10-25 14.11.22

次はImport Datasetを選びます。
スクリーンショット 2019-10-25 14.07.18

ポイントは2つです。
①データの名前(変数名)を付ける
何でもいいのですが今回はハルさんのサイトと同じgripにしました。

②読み込む範囲を指定する
今回A1からのデータではないので先程確認したB2からC62を指定します。
B2:C62のように左上と右下を:でつなげます
スクリーンショット 2019-10-25 14.17.12


そして右下にコードが自動的に作られます。Importを押せば完了なのですが、このコードをコピーしスクリプトに貼り付けておけば1年たった後でも同じことができます。EZRでもスクリプトを保存することができないわけではないのですが、再現性(後でしても、他の人がしても同じ事ができる)を保つためにもこういったコードを残しておく習慣をつけるようにしましょう。
スクリーンショット 2019-10-25 14.19.26


コードの一番下にあるView関数を使うことでRStudio内でもデータの確認ができます。このタブを消してもデータに影響はありません。もう1回View関数を使えばまた表示できます。
スクリーンショット 2019-10-25 14.28.00

View関数はEZRで言う表示と同じです(下図はハルさんのサイトより。比べてみてください)



これでデータの取り込みは完了です!



5.データの要約

ハルさんは次にデータの要約をしています。
EZRでのデータの要約と全く同じ機能はないですが、第2章で紹介したtidyverseパッケージのgroup_by関数とsummarize関数が使えます。group_by関数とsummarize関数に関してはこちらで紹介しています。

%>%やgroup_by関数、summarize関数はtidyverseパッケージに含まれていますのでtidyverseパッケージを呼び出します。もしtidyverseパッケージを全く使ったことが無い方はパッケージをインストールします。1度でも使ったことがあれば次の1行は必要ありません。

install.packages("readxl")

実際のコードは以下になります。イメージ図も添付します。

library(tidyverse)
grip %>% 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())

スクリーンショット 2019-10-25 20.52.13

スクリーンショット 2019-10-26 7.55.30


これでA群とB群のデータのばらつきを確認することができます。
ちなみにこのコードをコピーして色がついた箇所を変更すれば他の場面でも使えます!

このままでもいいのですが、データ要約は後でグラフ作成に使うのでgrip_summaryという名前をつけます。

grip_summary <- 
grip %>% 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())


6.正規性の確認

次に正規性の確認を行います。ハルさんのサイトではヒストグラムを作成しました。
ヒストグラムの作り方はこちらで紹介しています。


ハルさんのヒストグラムは棒が5本だったので同じ形にするようbins = 5とします。
A群とB群の棒を横に並べるときはposition = "dodge"を使います。
ggplot()で行を変える時は %>% ではなく + を使うので注意してください。

ggplot(data = grip) + 
  geom_histogram(aes(x = grip, fill = category), position = "dodge", bins = 5) 

スクリーンショット 2019-10-25 20.55.06



加えてシャピロ・ウィルク検定も紹介されています。シャピロ・ウィルク検定はshapriro.test関数を使います。
ただA群とB群それぞれで行いますのでgrip$gripの列をA群とB群に分ける必要があります。

shapiro.test(grip$grip[category == "A"])
shapiro.test(grip$grip[category == "B"])



データの中で特定の条件だけを抜き出すには[ ]を使います。



ハルさんのサイトのEZRで行う場合と見比べてみてください。
category == "A" と書いてあるところの意味が見えてきます。
ちなみに変数(1つ選択)のgripはgrip$gripの色のついた部分です。

スクリーンショット 2019-10-25 21.13.57


結果は以下のとおりです。
スクリーンショット 2019-10-25 19.36.44

どちらも0.5を超えているので正規分布であるという仮説は棄却されませんでした。
このあたりの解釈はハルさんのサイトをご参照ください。


7.1度に計算するsplit + map関数

shapiro.test(grip$grip[category == "A"])
shapiro.test(grip$grip[category == "B"])

上記のように1つずつ計算する方法もいいのですが群の数だけ繰り返します。プログラミングであるRは繰り返しに強いという特徴があります。


群ごとにデータを分割し、まとめて計算する方法として今回はsplit関数map関数を使います。


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

スクリーンショット 2019-10-26 22.42.55

まずgripのデータを
split関数を使ってA群とB群の2つのデータに分割し、map(~シャピロウィルク検定)でシャピロウィルク検定を繰り返します。

スクリーンショット 2019-10-26 23.01.43

EZRだと1つ1つ検定を繰り返す必要があります。
Rを活用すると1度に作業が終わるため、繰り返すことでのやり間違えの予防や時間の短縮にも繋がります。


8.等分散性を確認する

先程は正規性を調べましたが、今回は等分散性を調べます。
F検定と言いますがRではvar.test関数を使います。



Rで行うとこうなります。
var.test(目的変数 
~ グループ)

var.test(grip$grip ~ grip$category)

もしくはdata = を使えば$が要らなくなります。

var.test(grip ~ category, data = grip)
スクリーンショット 2019-10-27 5.57.05

p値 > 0.05だったので等分散性は棄却されませんでした。



9.t検定を行う

いままで正規分布かどうか、等分散性かどうかを確認しました。

スクリーンショット 2019-10-27 6.17.29

これでt検定が使えます。t検定はt.test関数を使います。
t検定は2種類の書き方があります。その違いはデータの形にあります。
ハルさんのサイトでEZRでは左の形しかできませんでしたがRならどちらも可能です
ちなみに左をlongデータ(縦に長い)、右をwideデータ(横に広い)といいます。
スクリーンショット 2019-10-27 21.34.00

スクリーンショット 2019-10-27 8.09.41

今回はlongデータなのでこのようになります。

t.test(grip$grip ~ grip$category, var = TRUE)

data = gripを使うとgrip$は外せます。どちらでも好きな方で大丈夫です。

t.test(grip ~ category, var = TRUE, data = grip)

スクリーンショット 2019-10-27 21.42.47
EZRと同じ結果になりました!


10.グラフの作成
EZRでは丁寧にグラフも作ってくれますがRでは自作する必要があります。
EZRでは平均と標準偏差の棒グラフを作成していました。
先程のgrip_summaryに平均と標準偏差を求めていたのでそれを使います。

棒グラフの作り方はここで紹介しています。


①最低限のグラフ(見た目は気にしない)

グラフを作る上で最低限必要な要素としては
棒グラフ:geom_bar(aes(x軸の指定、y軸の指定), stat = "identity")
エラーバー:geom_errobar(aes(x軸の指定, エラーバーの下端, エラーバーの上端))
y軸の名前が「平均」になるので「grip」に変更する

ggplot(data = grip_summary)+
  geom_bar(aes(x = category, y = 平均), stat = "identity") +
  geom_errorbar(aes(x = category, ymin = 平均 - 標準偏差, ymax = 平均 + 標準偏差)) + 
  labs(y = "grip") 
スクリーンショット 2019-10-27 22.59.56
ただこれでは見るに耐えません・・・



②見栄えを変更
次は以下を修正します。
棒グラフ:周りの線をblack、中の色をgray
エラーバー:幅を補足する(今回は0.1)

ggplot(data = grip_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")
スクリーンショット 2019-10-27 23.00.05
だいぶ見た目が良くなりました。



③EZRのグラフにできるだけ近づける
EZRのグラフに似せるためにもうひと工夫します。
背景を白にする(グラフのテーマをclassicに変える)
文字を大きくする(今回はsizeを15に変更)
棒グラフが浮いているように見えるのを修正する

ggplot(data = grip_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))
スクリーンショット 2019-10-27 23.00.12
これでEZRにだいぶ似ました。



11.まとめ

かなり長くなりましたが今回はt検定を紹介しました。

最初だったので1つ1つ説明しました。ボリュームが多すぎて慣れないところもあると思いますが、これから何度も出てくるものもありますので少しずつ慣れて貰えればと思います。


12.コードをまとめました
今回のコードを1つにまとめました。
パッケージはlibrary()を1度だけ行えばいいので最初に出しています。

#パッケージの読み込み
library(readxl)
library(tidyverse)

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/03/demo-t-test.xlsx"
destfile = "data/demo-t-test.xlsx"

download.file(url, destfile)

#データの読み込み
grip <- read_excel("data/demo-t-test.xlsx", 
                   range = "B2:C62")
View(grip)  

#データの要約
grip_summary <- 
  grip %>% 
  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_summary

#正規性を調べる
#ヒストグラム
ggplot(data = grip) + 
  geom_histogram(aes(x = grip, fill = category), position = "dodge", bins = 5) 

#シャピロ・ウィルク検定
shapiro.test(grip$grip[category == "A"])
shapiro.test(grip$grip[category == "B"])

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


#等分散性を確認
var.test(grip ~ category, data = grip)

#t検定
t.test(grip$grip ~ grip$category, var = TRUE)

t.test(grip ~ category, var = TRUE, data = grip)

#グラフ化
ggplot(data = grip_summary)+
  geom_bar(aes(x = category, y = 平均), stat = "identity") +
  geom_errorbar(aes(x = category, ymin = 平均 - 標準偏差, ymax = 平均 + 標準偏差)) + 
  labs(y = "grip") 

ggplot(data = grip_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")

ggplot(data = grip_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))

                


↑このページのトップヘ