カテゴリ:統計 > R

RStudioではまずExcelやcsvファイル等を読み込む必要がありますが、そもそもExcelではRに読み込みやすい形とそうでない形があります。Excelの見栄え的には良くても実際にRで読み込もうとすると大変な手間になることもあります。

今回はRで読み込みやすいExcelのポイントを紹介します。



1.上や左に空白はないか?

スクリーンショット 2019-10-20 20.07.06

RでなにもせずにExcelファイルを読み込むとA1を始めとしてデータを読み込みます。よくExcelでかけ線を設定している場合B2から表をつくることがありますが、Rで読み込むときには一手間必要です。



2.一番下に合計の行がないか?

加えてもし一番下に合計や平均の列があったらそれも1データとして読み込んでしまうので注意が必要です。
Excelは集計のために使い分析をRで行う場合はそもそも無いほうが安全です。
もしExcelでも集計を行いたい場合は集計を別のタブで行うかピボットテーブルを使う方法があります。



3.タイトル行が2行に渡っていないか?

タイトル行が2行に渡っているとRではうまく読み込んでくれません。
そしてそういう場合はセル結合されていることも多いです。
こういった場合も対処が必要です。

タイトルが2行に渡った時の対処法に関しては以下の記事で説明しています。




4.無駄な空間を作らない
スクリーンショット 2019-11-24 6.59.39
上記の図のように空欄だとパソコンは読み込んでくれません。
列名をちゃんと書き、データは空欄にせず入力しましょう。

5.桁数に使う , を入れていないか

10,000などカンマ含む場合、Rで読み込むと数字ではなく文字として認識されます。
「どうしても , はつけないとダメ」と上司に言われてもRで , を取る方法はあります(parse_number関数)


6.マイナスの値を▲にしていないか?

-100を▲100にしていたらマイナスに直す必要があります。
文字の置換(ここでは▲ → − )にするにはstr_replaceなどの関数がありますが、手間を考えると−の方が手間が省けます。


7.和暦か西暦か?

和暦だと古いExcelではいつまでたっても令和になりません。
計算や列の順番を考えると西暦のほうが後々対処しやすいです。
加えて2001年を01年としないよう注意しましょう。


8.tidyなデータか?

データ分析を行う上で必要な考え方にtidy(タイディー:整然データ)なデータかどうかという考え方があります。

tidyに関しては以下の記事が具体例も交えながらわかりやすく説明されています。




これからデータ収集を考えている方は参考になると思います。


9.被験者間要因は左側に、被験者内要因は右側に並べる
先程tidyなデータの話がありましたが、実際にデータを取りExcelに打ち込む時に覚えておくと後で便利になる考え方があります。それは被験者間要因被験者内要因です。
スクリーンショット 2019-11-24 6.33.40
被験者間要因はいわゆる対応のないデータのことで被験者全体をA法、B法の2つに割り振るといった方法です。なので同じ人がA法とB法を行うことはありません。

被験者内要因はいわゆる対応のあるデータのことで同じ人が繰り返し測定します

加えて被験者間要因は列の左側に、被験者内要因は列の右側に並べると後々分析しやすくなります。
また被験者間要因は縦に、被験者内要因は横につなげるとEZRで分析を行うには都合がいいです。
もし横に列が長くなりすぎて困る場合は全部縦にしても大丈夫です。あとで横に変換できます。
スクリーンショット 2019-11-24 6.35.12
この方法は特に分散分析を行う場合に便利になります。



10.まとめ

予めExcelのデータをRで読み込みやすい形にするとデータの前処理で心が折れにくくなりますので参考にしてみてください!


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


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


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


 


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


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


今回は一元配置分散分析を紹介します



まず一元配置分散分析(ANOVA)についてはハルさんのサイトをご参照ください。



また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-08 14.32.40


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

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

デモデータ(ANOVA)

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


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


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

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

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


4.データの読み込み

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

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

スクリーンショット 2019-11-08 14.24.20



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


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

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

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

スクリーンショット 2019-11-07 23.40.57



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

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

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

5.正規分布の確認

ハルさんのサイトに沿って正規分布の確認を行います。

正規分布の確認は【4-1】Rでt検定を行う方法で紹介しています。
以下のコードがイメージできない場合はまずこちらの記事をご参照ください。



まずヒストグラムで分布を確認します。
1つずつグラフを作るのは手間がかかるので、まとめてグラフを作ります。

ヒストグラムはgeom_histgram関数を使います。
グラフが重なっていると見にくいのでfacet.gridでカテゴリーごとにグラフを分けます。
facet.gridに関しては【3-4】Rのggplot2でヒストグラムを作るgeom_histogram関数で紹介しています。

library(tidyverse)
ggplot(data = grip_anova) +
  geom_histogram(aes(x = grip, fill = category), bins = 5) +
  facet_grid(category ~ .)
スクリーンショット 2019-11-08 15.27.34

6.シャピロ・ウィルク検定
次にハルさんのサイトでは正規性の検定を行っています。
正規性の検討はshapriro.test関数3つのカテゴリーそれぞれに行う必要があります。

今回はsplit関数map関数を使いまとめて行ってみます。
下の図は【4-1】Rでt検定を行う方法で使った方法ですが今回のコードと見比べてみてください。
コードがほとんど同じです。コードを使うと他でも役に立つことがわかります。

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

grip_anova %>% 
  split(.$category) %>% 
  map(~shapiro.test(.$grip))
スクリーンショット 2019-11-08 16.19.14

まとめて結果が出ました。すべて正規分布であることが棄却されませんでした。


7.等分散性の確認

次に等分散性の確認を行います。

3群以上で等分散性を検定するにはBartlett検定を使います。

bartlett.test(数値の列 ~ グループ)

data = を使えば$は必要ありません。どちらでも大丈夫です。
bartlett.test(grip_anova$grip ~ grip_anova$category)

bartlett.test(grip ~ category, data = grip_anova)



8.分散分析(ANOVA)を行う

それでは本題の一元配置分散分析と多重比較を行います。

実は一元配置分散分析と重回帰分析(t検定も)は同じ分析方法です(詳しく知りたい方は一般線形モデルで検索してみてください)。出てきた結果を分散分析の形式で表示するか重回帰分析の結果で表示するかの違いだけです。
そのため一元配置分散分析を行うには分散分析のaov関数重回帰分析のlm関数で行います。
ただその後に行う多重比較ではlm関数が使えないものもあるので基本はaov関数で大丈夫です。

どちらの関数もまずモデルを作りサマリーを表示するというのが一般的な流れです。

全体的な流れは以下の通りです。
sumamry関数だけ結果の表示が変わるので注意が必要です。
スクリーンショット 2019-11-11 20.05.02


aov関数の場合

モデル名:model_grip_anova(好きな名前でOK)
従属変数:grip
独立変数:category
データ名:grip_anova
model_grip_anova <- aov(grip ~ category, data = grip_anova)
summary(model_grip_anova)
スクリーンショット 2019-11-11 20.09.17


lm関数の場合

モデル名:model_grip_lm(好きな名前でOK)
従属変数:grip
目的変数:category
データ名:grip_anova
model_grip_lm <- lm(grip ~ category, data = grip_anova)
summary.aov(model_grip_lm)
スクリーンショット 2019-11-11 20.10.45



結果はどちらもEZRと同じになりました。


9.多重比較を行う

ハルさんのサイトでは多重比較でTukeyを紹介していますが他も紹介します。
Tukeyだけ関数が違います。
スクリーンショット 2019-11-11 22.18.55


Tukeyの場合
Tukeyの場合は先程のaov関数の結果を使います。
lm関数の結果は使えませんので注意が必要です。
TukeyHSD(model_grip_anova)
スクリーンショット 2019-11-11 22.26.49


Bonferroniの場合
pairwise.t.test関数はモデルの結果を使いません。
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "bonferroni")
スクリーンショット 2019-11-11 22.24.31



Holmの場合
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "holm")
スクリーンショット 2019-11-11 22.24.40


Hochbergの場合
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "hochberg")
スクリーンショット 2019-11-11 22.24.46



10.まとめ
今回は一元配置分散分析を紹介しました。
今回は割愛しましたが、一元配置分散分析を行うときはlongデータである必要があります。

Excelでデータ収集するときに注意するか、wideデータの場合はlongデータに直してから分析を行ってください。

longデータとwideデータについてはで【2-5】Rでデータを集計するのに便利なlongデータとgather関数
紹介しています。

次回はKruskal-Wallis検定を行います。

(追記)グラフの作成
EZRでは検定と同時にグラフを作成しますが、Rでは直接作成します。
【4-1】t検定でグラフを使ったときのコードがそのまま使えます。


グラフを作るために必要な要約のコードもそのまま使います。
#データの要約
grip_anova_summary <- 
  grip_anova %>% 
  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_anova_summary
#グラフ化 ggplot(data = grip_anova_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-11-15 20.13.23
コードの違いはわかったでしょうか?
t検定で使ったデータであるgripを今回のデータであるgrip_anovaに変更しただけです。
それ以外のコードの中身は全く一緒です!

11.今回のコード
今回使ったコードのまとめになります。
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-anova.xlsx"
destfile = "data/demo-anova.xlsx"

download.file(url, destfile)

library(readxl)
grip_anova <- read_excel("data/demo-anova.xlsx", 
                         range = "B2:C62")
View(grip_anova)


#正規分布の確認
library(tidyverse)
ggplot(data = grip_anova) +
  geom_histogram(aes(x = grip, fill = category), bins = 5) +
  facet_grid(category ~ .)


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

#Bartlett検定

bartlett.test(grip_anova$grip ~ grip_anova$category)

bartlett.test(grip ~ category, data = grip_anova)

#一元配置分散分析
#aov関数の場合
model_grip_anova <- aov(grip ~ category, data = grip_anova)
summary(model_grip_anova)

#lm関数の場合
model_grip_lm <- lm(grip ~ category, data = grip_anova)
summary.aov(model_grip_lm)

#多重比較
#Tukey
TukeyHSD(model_grip_lm)

#Bonferroni
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "bonferroni")

#Holm
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "holm")

#Hochberg
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "hochberg")

#データの要約 grip_anova_summary <- grip_anova %>% 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_anova_summary
#グラフ化 ggplot(data = grip_anova_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))


全記事の目次はこちらをご参照ください






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


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


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


 


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


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


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



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



また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-08 14.29.27




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

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

デモデータ(McNemar検定)

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


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


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

download.file(url, destfile)

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


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




4.データの読み込み

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

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

スクリーンショット 2019-11-07 23.37.11


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




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



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

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

スクリーンショット 2019-11-07 23.40.57



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

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



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


5.データテーブルの作成

Fisherの正確検定と同様にまずはデータテーブルを作成してみます。
データテーブルを作成するのはtable関数でした。
table(bw$Before, bw$After)
スクリーンショット 2019-11-08 6.14.43

これでデータテーブルができました。後で使うのでtable_bwという名前を付けておきます。
(必須ではありません)
table_bw <- table(bw$Before, bw$After)

6.McNemar検定を行う

次にMcNemar検定を行います。McNemar検定はmcnemar.test関数をを使います。

ncnemar.test関数は2つの記述方法があります。

①集計前の表(Excelの形式)を使う
今回のExcelの形でmcnemar.test関数を使う事ができます。

mcnemar.test(1列目, 2列目)
mcnemar.test(bw$Before, bw$After)

②集計後の表(テーブル形式)を使う
先程のテーブル形式(table_bw)をそのまま使うこともできます。
他にも既にExcelで集計してしまった場合です。

mcnemar.test(tableデータ)
mcnemar.test(table_bw)
スクリーンショット 2019-11-08 12.39.55

どちらもEZRと同じ結果になります。


7.まとめ

今回はRでMcNemar検定を行う方法を紹介しました。
Fisherの正確検定、カイ二乗検定、McNemar検定は似た使い方ですので続けて記事を見られた方は行いやすかったかもしれません。

またこのページの最後に連続補正について紹介しています。

次回は分散分析に移ります。



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

download.file(url, destfile)

library(readxl)
bw <- read_excel("data/demo-mcnemar-test.xlsx", 
                 range = "B2:C22")
View(bw)

#テーブルの作成
table_bw <- table(bw$Before, bw$After)
table_bw

#McNemar検定
mcnemar.test(bw$Before, bw$After)

mcnemar.test(table_bw)

9.(追記)連続修正による正規近似
McNemar検定はカイ二乗検定を元に近似値を求めています。


mcnemar.test関数はデフォルトで連続修正を用いた正規近似を用いています。
もし連続修正を用いない場合はcorrect = FALSEを付ける必要があります。
mcnemar.test(bw$Before, bw$After, correct = FALSE)

サイトによって連続補正を付けるか付けないか意見はあるようです。

連続補正あり(厳しめな検定、第1種の過誤を生みやすい)
直接検定
連続補正なし(甘めな検定、第1種の過誤が減る)


総じて言うと上記の傾向のようですが、気になる方は下記をご参照ください。



こちらでは連続補正を施したほうが近似値が良くなるとありました。
この近似はかなり正確なことと、p値がだいたい0.01くらいまでは連続修正を施した方が近似が良くなることがわかると思います。 そのためウィルコクソンの1標本検定と違って、通常はzoの値によらず全て連続修正を施します。


Rコマンダーのサイトでは理由は書いてありませんでしたが、correct = FALSEと言っています。
mcnemar.test(.Table, correct=FALSE)を追記して下さい。

マ ッ チドペ アにおける 2 × 2表の検定の第 1種の 過誤 ... - J-Stage


https://www.jstage.jst.go.jp › article › jscswabun › _pdf › -char

こちらの記事では連続補正をすると第1種の過誤が生まれやすいので無条件に使うには連続補正を行わないほうが良いといった記載がありありました。



こちらの記事ではnデータ数によって連続補正あり・なしを指定しています。



また、もともとはパソコンの性能上厳密な検定ができずカイ二乗検定をつかって近似値を求めていたようですが、二項検定やパッケージを使ってp値を厳密に計算することもできるようです。それに関しては奥村先生の記事をご参照ください。






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



↑このページのトップヘ