タグ:tidyverse

第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の理解も深まりやすいと思います。



前回は今までの総復習を行いました。
【演習】R初心者が統計をかけるための前準備の流れを復習します : 独学で始める統計×データサイエンス



上記の記事を進めると以下のように統計量をまとめて出すことができますが、summary関数では標準偏差が出ません。
スクリーンショット 2019-02-04 21.49.17


今回は標準偏差を出すいろいろな方法を復習も含め紹介します。

繰り返しのfor関数やapply関数、パイプ演算子(%>%)もはじめて出てきますが、使えるようになるとFIMのような項目が多いものでも1度にまとめて解析できるようになり、いよいよExcelではできない体験をすることができるようになります。

まだ今回の内容であればExcelでも可能ですが、1つずつ進めていきましょう。

1.sd関数で1つずつ出す
2.for関数で繰り返す
3.apply関数を使う
4.tableoneパッケージを使う


今回は前回の記事のデータを使います。

data01.xlsx 

ダウンロードした後、プロジェクトの指定フォルダにファイルを移動させておけば、以下のコマンドで上記の画面まで進みます。

library(readxl)
data01 <- read_excel("data01.xlsx", sheet = "日本語")
data01$性別 <- factor(data01$性別, levels = c("男性","女性"))
data01$歩行 <- as.factor(data01$歩行)
summary(data01)

*注意
もしプロジェクトファイルを作成していない場合はgetwd()と打ち込むと作業フォルダが表示されます、作業フォルダにExcelデータを入れると後は同じです。

getwd()




1.sd関数で1つずつ出す

標準偏差を出す関数はsd関数です。

sd(ベクトル)
sd(data01[[3]])
sd(data01[[4]])
sd(data01[[5]])
sd(data01[[6]])
sd(data01[[7]])
sd(data01[[8]])
スクリーンショット 2019-02-04 21.53.36


data.frameからベクトルを表示するには$または[[ ]]を使います。
[[ ]]の中の数字は列番号です。
$でもいいのですが、このようなケースでは列番号のほうが手っ取り早いです。

列番号を確認するにはnames関数t関数(もしくはdata.frame関数)を使うと便利です。

t(names(data01))
スクリーンショット 2019-02-04 23.20.52

data.frame(names(data01))
スクリーンショット 2019-02-04 23.22.06


2.for関数で繰り返す

1個ずつ出すのもいいですが、繰り返しになる作業はfor関数で繰り返すことができます。
for(i in 3:8){
  x <- sd(data01[[i]])
  print(x)
  }

for関数は(◯ in △)と{}の2つに別れています。
スクリーンショット 2019-02-04 23.51.05
スクリーンショット 2019-02-04 23.59.46


スクリーンショット 2019-02-05 0.03.59

ただ結果は出るのですが結果の数字しか出ません。


3-1.apply関数を使う

library(tidyverse)
apply(data01, 2, sd)

スクリーンショット 2019-02-05 0.17.43

スクリーンショット 2019-02-05 0.08.54


下にエラーが出ます。これは「氏名と性別で集計できないからNAにしてますよ!」と怒られています。なのでselect関数で集計する列だけ抜き出します。


library(tidyverse) #既に実行していたらなくても可
data01_select <- select(data01, 3:8) data01_select apply(data01_select, 2, sd)
スクリーンショット 2019-02-05 0.32.27

今回はdata01_selectという変数名にselect関数で3〜8列目だけ抜き出しapply関数を使いました。


スクリーンショット 2019-02-05 0.33.14


ただこの方法はdata01_selectのように間に変数を挟む必要があります。
tidyverseパッケージはこれをスッキリさせるパイプ演算子(%>%)というものがあります。


3−2.パイプ演算子を使ったapply関数

パイプ演算子(%>%)は解析の流れをスムーズにしてくれます。

データ① %>%
    プログラム② %>%
    プログラム③ %>%
    プログラム④

とすると先程のdata01_selectの用な変数を介することもなく、データ①をプログラム②に送り、その結果をプログラム③に送り、その結果をプログラム④に送るといったことが可能になります

library(tidyverse) #既に実行していたらなくても可
data01 %>% #data01に対して
  select(3:8) %>% #3列目(年齢)〜8列目(MMSE)までを選択し
  apply(., 2, sd) #それぞれの列で標準偏差を求めて!

スクリーンショット 2019-02-05 0.55.55
スクリーンショット 2019-02-05 0.57.27


もし四捨五入するならround関数を使います。

library(tidyverse)#既に実行していたらなくても可
data01 %>% #data01に対して
  select(3:8) %>% #3列目(年齢)〜8列目(MMSE)までを選択し
  apply(., 2, sd) %>% #それぞれの列で標準偏差を求めて!
  round(.,2) #小数点2位より下を四捨五入する
スクリーンショット 2019-02-05 1.13.08


4.tableoneパッケージを使う

tableoneパッケージを使うのも1つの方法です。

tableoneパッケージについてはこの記事で紹介しています。
Rで医療統計で必要なtable1を作るtableoneパッケージについて紹介します : 独学で始める統計×データサイエンス

library(tableone)
CreateTableOne(data = data01,
               vars = c("年齢", "身長", "体重", "SIAS", "BBS", "MMSE"))
スクリーンショット 2019-02-05 1.20.05


まとめ

今回は標準偏差を出すいろいろな方法を通じてfor関数やapply関数、パイプ演算子の使い方を紹介しました。

今回は標準偏差(sd)でしたが、平均(mean)や中央値(median)など他の統計量も求めることができます。

for関数やtidyverse(dplyrなど)に関してはまとまった記事もありますが、これからもちょくちょく出てくると思います。その都度復習しながら使い方のイメージを付けていただければいいかと思います。








↑このページのトップヘ