第4章では様々な検定を紹介しています。

検定を行うとp値に目が向きがちですが、期待値と95%信頼区間を出す必要性も言われています。

今回はp値と95%信頼区間の関係性についてです。

1. p < 0.05 と95%信頼区間の関係

統計ソフトで検定をかけるとp値が出てきます。
ただ実はp値を見なくても95%信頼区間を見ればp < 0.05かどうかはグラフを見れば一目でわかります。
p < 0.05は95%信頼区間が基準を跨ぐかどうかと全く同じ意味です

スクリーンショット 2019-11-05 21.58.19
例えばt検定など2つの基準の差を見る検定では差が0という仮説を立てます。
95%信頼区間に0を含むかどうかとp < 0.05 かどうかは全く同じ意味です。

またFisherの正確検定やロジスティック回帰分析のようにオッズ比をみる検定ではオッズ比が1という仮説を立てます。
95%信頼区間に1を含むかどうかとp < 0.05 かどうかは全く同じ意味です。



1.p値だけ vs 信頼区間


たとえば架空のダイエットのメソッドがあったとします。数字は完全に適当です。

3つの架空のメソッドで3ヶ月トレーニングをした研究があったとします。ちなみに1つの研究ではなく、3つの別々研究だったとします。それぞれN数も違うとします。交絡因子などの共変量もここでは考えないこととします。

①p値のみで判断

メソッドA:p < 0.00001
メソッドB:p = 0.35
メソッドC:p = 0.11

この場合どれが効果ありと感じるでしょうか?

②期待値と信頼区間を見る

次に以下の期待値と95%信頼区間をみるといかがでしょうか?
マイナスになっているというのはそれだけ減量したという意味とします。
p < 0.05は95%信頼区間が0を跨ぐかどうかと全く同じ意味です

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


③N数も追加

更にN数も追加してみます

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


④95%信頼区間の特徴

95%信頼区間の幅はN数が増えると短くなるという特徴があります。
そのためどんな意味のない些細な効果でもN数を増やすと有意差が出やすくなります。
(倫理的にどうかはここでは言及しません)

Aは有意差はあるのですが、逆に言うと3ヶ月トレーニングして1kgの効果しかないメソッドとも言えます。更にN数を増やすと区間が狭くなるので有意差はより出るかもしれませんが、減量効果が上がる(点の位置が左に動く)訳ではありません。

Cはもしかしたらサンプルサイズの設計(取るべきN数)が足りなかっただけという研究の設定段階の問題なのかもしれません。もちろんN数を増やすことで信頼区間の幅だけでなく期待値もどう変わるかわかりませんが、有意差がなかったので効果が無いと断言してしまっていいのかどうかは気になります。


このようにp < 0.05の有無だけではどのくらい(意味のある)効果があったのかを判断できないこともありますし、p >= 0.05 だけで効果が無いと断言できるものでもありません(本当に意味のある効果がなかったのか、本来あったはずなんだけどN数が足りなくて効果がなかったのか判断つかない)。

そこに95%信頼区間があればp値のみよりも判断する材料が増えるかもしれません。


2.まとめ

今回はp値と信頼区間の関係性について紹介しました。
もちろん研究デザイン(共変量の調整など)やバイアスを評価することが大切なのですが、加えてp値だけにとらわれないように自分も気をつけたいと思います。





第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%信頼区間のグラフを作る方法を紹介します。



↑このページのトップヘ