タグ:RStudio

第4章では統計の中でも検定を扱っています。

ここまでいろいろな検定を行ってきましたが、グラフを作る時はどうするでしょうか?

【4-6】Rのt検定の結果からp値や信頼区間の数値を取り出す方法では検定した結果を取り出す方法を紹介しました。

今回は取り出したデータで95%信頼区間のグラフを作る方法を紹介します。

データとt検定のコードは以下になります。
前回の記事でこのコードの解説をしていますのでわからない場合は先に確認をお願いします。


set.seed(2019)
male <- rnorm(50, 170, 10)
set.seed(2019)
female <- rnorm(50, 165, 10)
height <- c(male, female)
sex <- c(rep("male", length(male)), rep("female", length(female)))
dat_height <- tibble(height, sex)
head(dat_height)

res_height <- t.test(height ~ sex, var = TRUE, data = dat_height)


1.データを取り出す
まずはp値、95%信頼区間のデータを取り出してみます。
これも前回の記事で紹介しています。

p:p値
conf.low:95%信頼区間の小さい方
conf.high:95%信頼区間の大きい方
p <- res_height$p.value
conf.low <- res_height$conf.int[1]
conf.high <- res_height$conf.int[2]
p conf.low conf.high
スクリーンショット 2019-11-07 18.28.31
この数字をグラフに乗せると小数点が多すぎるのでround関数で丸めます。
今回はp値は小数点第3位、95%信頼区間は小数点第2位までにしてみます。
p <- round(res_height$p.value, 3)
conf.low <- round(res_height$conf.int[1], 2)
conf.high <- round(res_height$conf.int[2], 2)
p conf.low conf.high
スクリーンショット 2019-11-07 18.32.21
これでグラフに貼り付けやすくなりました。

2.グラフを作成する

今回はうまくいかないグラフも載せることでグラフを作る流れもお見せできればと思います。
目標は以下のようなグラフを作ることです。
スクリーンショット 2019-11-07 22.35.17


第3章で使ってきたtidyverseパッケージのggplot関数を使います。
もしggplotの使い方がわからない場合はこちらをご参照ください。




最初にtidyverseパッケージを呼び出します。
library(tidyverse)
まずグラフを作ってみます。
今回使うのはgeom_errorbar関数を使います。
geom_errorbar関数はaes関数の中に3つの要素が必要です。

x
:x軸
ymin:エラーバーの最小値
ymax:エラーバーの最大値

x軸に値するものは今回ないので空欄にしておきます。
yminとymaxはconf.lowとconf.highになります。

x:""
ymin:conf.low
ymax:conf.high

ggplot() +
  geom_errorbar(aes(x = "", ymin = conf.low, ymax = conf.high))
スクリーンショット 2019-11-07 20.35.07
なんだかすごいグラフができました。
まずは横向きにします。
グラフを横向きにするのはcoord_flip関数です。()には何も入れません。
ggplot() +
  geom_errorbar(aes(x = "", ymin = conf.low, ymax = conf.high)) +
  coord_flip()
スクリーンショット 2019-11-07 20.36.15
これで横向きになりましたが線が幅長いので短くします。
width = で指定します。今回は0.1にしました。widthはaes関数の外に配置します。
ggplot() +
  geom_errorbar(aes(x = "", ymin = conf.low, ymax = conf.high), width = 0.1) +
  coord_flip()
スクリーンショット 2019-11-07 20.39.18
これでエラーバーっぽくなりました。
今度は数値を入れます。
文字を打ち込むのでgeom_text関数を使います。
geom_text関数のaes()では3つ指定します。
coord_frip()を使っているのでx軸が縦、y軸が横になっていることに注意します。

x:文字を置くx軸の位置→今回は空欄
y:文字を置くy軸の位置→conf.lowとconf.highの2つ
label:実際の文字→conf.lowとconf.highの2つ


conf.lowとconf.highの2つのデータを入れるのでc関数でつなげます。

x:""
y:c(conf.low, conf.high)
label:c(conf.low, conf.high)

ggplot() +
  geom_errorbar(aes(x = "", ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(x = "", y = c(conf.low, conf.high), label = c(conf.low, conf.high))) +
  coord_flip()
スクリーンショット 2019-11-07 20.51.35
数字が出てきましたがグラフとぶつかってしまいます。
vjust =で調整します。
上に上げるときはマイナスの値を入れます。今回は試してみて-1.5にしました。
ggplot() +
  geom_errorbar(aes(x = "", ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(x = "", y = c(conf.low, conf.high), label = c(conf.low, conf.high)), vjust = -1.5) +
  coord_flip()
スクリーンショット 2019-11-07 22.24.29
これでグラフができましたが、y軸(グラフを横にしたから下がy軸になっている)の名前を変えてみます。labs関数を使って「95%信頼区間」と入れてみます。xも消します(空欄にする)。
ggplot() +
  geom_errorbar(aes(x = "", ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(x = "", y = c(conf.low, conf.high), label = c(conf.low, conf.high)), vjust = -1.5) +
  labs(x = "", y = "95%信頼区間") + 
  coord_flip()
スクリーンショット 2019-11-07 21.00.00
Windowの場合はこれで完成かもしれませんが、Macだと日本語が□□□と豆腐になってしまいます。
Macの方はtheme_◯◯関数base_family = を指定します。
今回はヒラギノ角ゴproW3を指定します。
ggplot() +
  theme_gray(base_family = "HiraKakuPro-W3") +
  geom_errorbar(aes(x = "", ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(x = "", y = c(conf.low, conf.high), label = c(conf.low, conf.high)), vjust = -1.5) +
  labs(x = "", y = "95%信頼区間") + 
  coord_flip()
スクリーンショット 2019-11-07 21.03.37
今回のt検定は2群の差が0と仮定していました。0のラインで赤線を足してみます。
y軸に垂直な線を引くにはgeom_hline関数を使います。
かならずいるのはyintercept = です。color = は付けなければ黒になります。
aesは付けなくて大丈夫です。

geom_hline(yintercept = ◯, color = "色名")

ggplot() +
  theme_gray(base_family = "HiraKakuPro-W3") +
  geom_hline(yintercept = 0, color = "red") +
  geom_errorbar(aes(x = "", ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(x = "", y = c(conf.low, conf.high), label = c(conf.low, conf.high)), vjust = -1.5) +
  labs(x = "", y = "95%信頼区間") + 
  coord_flip()
スクリーンショット 2019-11-07 22.35.17

3.まとめ
今回は検定の結果から値を取り出してグラフを作成しました。

グラフの要素を1つずつ追加しました。ブログに書かれているコードは長くて読みにくいかもしれませんが、「今度はこの要素を付け加えたい」という順番にコードを書けばその通りにグラフができるのもRの特徴の1つです。

今回の順番は1例ですので順番を変えて1つずつ作ってみていただければggplotの理解も深まりやすいと思います。



第4章では統計の中でも検定を扱っています。

ここまでいろいろな検定を行ってきましたが、グラフを作る時はどうするでしょうか?

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

これは【4-1】Rでt検定を行う方法で行ったt検定の結果です。
p値があり、下には95%信頼区間、その下には各グループの平均が記載されています。

p = 0.00975
95%信頼区間:-12.154877 〜 -1.745123
A:28.12
B:35.07

この数値をグラフなどにするには手で打ち直すかコピペをすることが多いと思います。
今回はこの数値がどこに入っているのか?どうやって取り出せばいいのかという話です。

今回は架空のデータを使いますがもし第4章を実践された方はそのデータでも大丈夫です。


1.使うデータ

こんなデータを作りました。
男性:50人、平均170cm、標準偏差10cm
女性:50人、平均165cm、標準偏差10cm

set.seed(2019)
male <- rnorm(50, 170, 10)
set.seed(2019)
female <- rnorm(50, 165, 10)
height <- c(male, female)
sex <- c(rep("male", length(male)), rep("female", length(female)))
dat_height <- tibble(height, sex)
head(dat)

rnorm関数(個数、平均、標準偏差)を指定すると疑似データを作ってくれます。
heightの列に男女のデータをc関数でくっつけ100人のデータにします。
sexの列は(男性,男性,男性,,,,女性,女性,女性,,,,,)としたいのでrep関数を使いました。
rep(繰り返したいもの, 回数)と使います。
tibble関数でheightとsexをつなげました。

2.t検定をする

これをt検定してみます。
t検定は【4-1】Rでt検定を行う方法で紹介しています。


等分散かどうかはそもそも等分散になるようデータを作りましたのでvar = TRUEを付けます。

t.test(height ~ sex, var = TRUE, data = dat_height)
スクリーンショット 2019-11-06 20.37.39

ここまでできました。

このままだとデータは取り出せません。結果を変数に入れます。
今回はres_heightとします
res_height <- t.test(height ~ sex, var = TRUE, data = dat_height)

3.結果からデータを取り出す

すると右上のEnvironmentタブに結果が格納されます。
スクリーンショット 2019-11-06 20.59.45


スクリーンショット 2019-11-06 21.05.31

データをみるとList of 10とあります。
リストというのはExcelのシートが10個あって、それぞれにデータが入っているというイメージです。
一番左がExcelでいうシート名です。
スクリーンショット 2019-11-07 0.13.47


リストの各データを呼び出すのは$もしくは[[ ]]を使います。
p値はp.valueに入っていますのでこうなります。
スクリーンショット 2019-11-07 0.30.11
res_height$p.value

res_height[["p.value"]]
スクリーンショット 2019-11-07 1.06.18

もし小数点以下を四捨五入したい場合はround関数を使います。

round(res_height$p.value, 3)

round(res_height[["p.value"]], 3)
スクリーンショット 2019-11-07 1.06.07


95%信頼区間はconf.intに入っています。
ただそのまま出すと2つのデータが出てきます(95%信頼区間の最小と最大)。
それぞれの値を取り出すには[ ] を使います。[[ ]]ではないので注意してください。
[ ]は【1-12】Rで特定の条件にあう要素を抜き出す方法で紹介しましたが、条件式ではなくただの数値をいれると◯番目のデータを取り出してくれます。

[1]だと1つ目のデータ(つまり最小)
[2]だと2つ目のデータ(つまり最大)が取り出せます。

res_height$conf.int[1]
res_height[["conf.int"]][2]
スクリーンショット 2019-11-07 1.03.41

そして点推定である各平均はestimateに入っています。
res_height$estimate
res_height$estimate[1]
res_height[["estimate"]][2]
スクリーンショット 2019-11-07 2.04.43



4.データを取り出すメリット
結果をコピペしても問題はないのですが、こうすることでのメリットがあります。

それはもし後でデータが変わったとしてもコード自体は何も変えなくていいことです。
考えたくはありませんが、後でデータが増えた、実は個々の数値が間違ってた等あると結果も変わります。しかしデータをコードで取り出しておくと、データが変わっても自動的に数値が更新されるのでミスする可能性が減ります。元データを修正するだけでいいというのはExcelやEZRでは簡単にはできません。

5.まとめ
今回は検定の結果からデータを取り出す方法を紹介しました。
検定の種類によっては取り出し方が違うこともありますが、まずは結果が格納されたデータをみるとヒントがあるかもしれません。

次回は取り出した値を使って95%信頼区間のグラフを作る方法を紹介します。



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

↑このページのトップヘ