カテゴリ:統計 > R

*最終更新日:2020/09/28
*記事を大幅に追加

医療統計に限らず、データを分析する時に集計表を出す場面は多くあります。

特に医療統計では表1(table1)として図表の最初に提示することが多いです。

スクリーンショット 2020-09-12 9.37.17

ポスター発表だとPower Pointでしょうか。
スクリーンショット 2020-09-12 9.38.30

論文であればWordを使うかもしれません。
スクリーンショット 2020-09-12 9.40.28

このtable1を作る時は基本次のような作業が必要になります。
・Excelで集計
・コピー
・WordやPower Pointに張り付け、もしくは手打ち
・体裁を整える

これが意外とめんどくさい作業なのですがRを使えば一度に行えます。
・Excelでデータを集める
・Rで集計→体裁を整える→Word, Power Point, 図として出力する

今回はgtsummaryパッケージを紹介します!


そしてこの記事を書いている最中にgtsumarryパッケージの作者であるDaniel Sjobergさん(@statistishdan)より直接アドバイスを頂きました。この場を借りてお礼申し上げます。



改善点やリクエストがあれば連絡が欲しいとのことでした。自分もいくつかリクエストさせていただきましたが、疑問点や改善点あればご指摘いただければ連絡してみるのはいかがでしょうか。




1.データの準備
今回は【1-11】Rで医療統計で必要なtable1を作るtableoneパッケージについて紹介しますで使用したデータと同じものを使います。


set.seed(1)
年齢 <- floor(rnorm(100,60,10))
性別 <- sample(c("男性","女性"),100,replace = TRUE)
体重 <- floor(rnorm(100,60,7))
MMT <- as.numeric(sample(c(1:5), 100, prob = c(0.1,0.1,0.2,0.3,0.3), replace = TRUE))
術側 <- sample(c("右","左"),100,replace = TRUE)
治療 <- sample(c(0,1),100,replace = TRUE)
data <- data.frame(年齢,性別,体重,MMT,術側,治療)
data
今回は治療(0:対照群、1:治療群)として分けていきます。


もしご自身のExcelデータで取り込み方がわからない場合は【4-0】第4章を進めていく上での準備をご参照ください。プロジェクトの使い方やExcelデータの読み込みについて紹介しています。


2.ライブラリの読み込み

一言でいうとライブラリとはRの機能を拡張してくれる機能の事です。

今回は4つのライブラリを使います
gtsummary(table1や結果の表をいい感じにまとめてくれる)
tidyverse(データの集計などに使用。今回はサブ的に使う)
gt(gtsummaryのデータをもっと細かく修正することができる。gtsummaryだけで完結するので必要ないが、表のタイトルやサブタイトルを付け加えたかったら必要)
flextable(できた表をWordやPower Point形式に貼り付ける際、gtsummaryよりも細かいく修正することができる。gtsummaryだけで完結するので必要ないが、表のタイトルやサブタイトルを付け加えたかったら必要)

rstudioのPackagesをみて自分のパソコンにインストールされているか確認してみてください。
このリストに入っていなかったり、最新のパッケージにしたい場合は以下のコードをコピペして実行してください。この操作は毎回する必要はありません。

install.packages("gtsummary")
install.packages("tidyverse")
install.packages("gt")
install.packages("flextable")

(追加)
記事作成時点(2020/09/23)でのバージョンは1.3.4です。
ただ1.3.4.9008以降で使えるtheme_gtsummary_mean_sd()がとても便利です。
現段階で最新版を使うためにはgithubからインストールする必要があります。

install.packages("remotes")
remotes::install_github("ddsjoberg/gtsummary")

ライブラリがインストールできていれば今度はパッケージを読み込みます。
基本Rstudioを起動すると基本のパッケージしか使えない状態です。
パッケージは毎回呼び出す必要がありlibrary関数を使います。
まずは今回の記事で必要なgtsummaryとtidyverseを読み込みます。

library(gtsummary)
library(tidyverse)

これで準備は完了です。


3.表を作る前に確認しておくこと


table1はそもそも各群のベースラインを比較します。

スクリーンショット 2020-09-12 9.38.30


そのためには各変数(評価)がどんなデータか理解している必要があります


数値かカテゴリーか?

集計を行うにはその変数が数値なのか?カテゴリー変数なのか?の理解が必要です。

数値:年齢・体重

カテゴリー:性別、MMT、術側

MMTは「1が7人、2が5人・・・」と集計を行いたいとすると、値は数値ですがカテゴリー変数として集計を行う必要があります。


数値の場合は正規分布なのか?そうでないのか?

同じ数値でも正規分布かそうでないのか?を決めておくことが必要です。

正規分布であれば平均±標準偏差ですが、正規分布でない場合は中央値(四分位範囲 or 最大値〜最小値)で表現したりします。


4.基本的な使い方

今から実際にgtsummaryを使っていきます。

まずは作ってみる

集計表を作るのはtbl_summary()関数を使います。

dataのデータをtbl_summary関数に入れ込みます。
これを表現するとdata %>% tbl_summary()と表記します。
%>%はパイプ演算子といい、最近Rでデータを扱う時の主流となっています。

data %>% 
  tbl_summary()
スクリーンショット 2020-09-12 13.43.17

とりあえず集計表ができました!ただ群別に表示がされていないので群を分けます。
tbl_summary(by = 群分けする列名)のようにby=○○を付け加えるだけです。

data %>% 
  tbl_summary(by = 治療)
スクリーンショット 2020-09-12 14.01.16
たった2行でそれらしい表ができてきました!


5.テーマを決める

gtsummaryはいくつかのテーマ(設定)があります。
表の細かい修正の前にテーマについて紹介します。


・ジャーナルに合わせたテーマに変える(JAMA / LANCET)
・日本語にする
・デフォルトの中央値(四分位範囲)を平均(標準偏差)に変える
・列幅を狭める
・元の設定に戻す

5−1.ジャーナルに合わせたテーマに変える(JAMA / LANCET)

table1といっても投稿するジャーナルで投稿規定が違います。
gtsummaryでは現在JAMAとLANCETのに合わせたテーマがあります。

#JAMAの投稿規定に合わせる
theme_gtsummary_journal(journal = c("jama"), set_theme = TRUE)
#LANCETの投稿規定に合わせる
theme_gtsummary_journal(journal = c("lancet"), set_theme = TRUE)

#元の設定に戻す
reset_gtsummary_theme()


5−2.日本語にする

英語→日本語にする機能もあります。

#日本語にする
theme_gtsummary_language("ja")
#元の設定に戻す
reset_gtsummary_theme()
スクリーンショット 2020-09-12 14.39.39

注釈などが日本語に変わりました。


5−3.デフォルトの集計・検定を中央値(四分位範囲)を平均(標準偏差)に変える

gtsummaryではデフォルトの集計・検定方法がノンパラメトリックとなっています(中央値・四分位範囲・Wilcoxonの順位和検定など)。
医療統計では平均(標準偏差)で表現することが多いです。個別に設定する方法は後ほど説明しますが、デフォルトでパラメトリックにする設定があります(平均・標準偏差・t検定など)。

theme_gtsummary_mean_sd()

#元の設定に戻す reset_gtsummary_theme()

*この機能はパッケージ開発者のgtsumarryパッケージの作者であるDaniel Sjobergさん(@statistishdan)に自分がリクエストしたところ機能追加していただきました。バージョン1.3.4.9008以降で使えます。ありがとうございます!



5−4.幅を狭くする

できた表はゆとりがありますがコンパクトにまとめることもできます。
theme_gtsummary_compact()を実行するだけです。

#幅をコンパクトにまとめる
theme_gtsummary_compact()

#もとの設定に戻す
reset_gtsummary_theme()
スクリーンショット 2020-09-12 14.45.41

5−5.設定を戻す

今まで説明したtheme_gtsummary_〇〇関数ですが、全てrecet_gtsummary_theme()を使えばデフォルトの設定に戻ります。

reset_gtsummary_theme()


6.p値の列を作る

次にp値の列を作ります。
作り方は%>%を使い、先程のコードにadd_p()を加えるだけです。

data %>% 
  tbl_summary(by = 治療) %>% 
  add_p()
スクリーンショット 2020-09-12 14.52.49
注釈までつけてくれました。
年齢と性別にはWilcoxonの順位和検定(マン・ホットニーのU検定)
性別と術側にはカイ二乗検定
MMTにFisherの正確確率検定が使われています。
「Wilcoxonの順位和検定じゃなくてt検定にして欲しい」に関しては先程説明した5−3.デフォルトの集計・検定を中央値(四分位範囲)を平均(標準偏差)に変えるをご参照ください。中には平均±標準偏差にしたり、これは平均だけどこっちは中央値など細かい調整をしたい場合もあるかもしれません。細かい設定は下にある14−3.集計方法を変更するをご参照ください。

*検定自体に関しては第4章に解説があります。サイトマップに一覧がありますので検定に関して知りたい場合は確認してください。



7.N数を加える
人数を加える場合はadd_n()を加えます。

data %>% 
  tbl_summary(by = 治療) %>% 
  add_p() %>% 
  add_n() 
スクリーンショット 2020-09-12 15.01.43
欠損値がないので全て100ですが、欠損値があるデータだと変わることはあります。


8.全体の集計の列を作る

もし全体の集計列を加えたければadd_overall()を加えます。
今回はadd_n()を外してadd_overall()を加えていますが、もちろんどちらも加えることも可能です。


data %>% 
  tbl_summary(by = 治療) %>% 
  add_p() %>% 
  add_overall() 
スクリーンショット 2020-09-12 15.05.19


9.変数名を太字にする
このままでもいいですが、変数名を太字にしてわかりやすくします。
bold_labels()を使います。

data %>% 
  tbl_summary(by = 治療) %>% 
  add_p() %>% 
  add_overall() %>%   
  bold_labels()
スクリーンショット 2020-09-12 16.04.08




10.出力する

まだまだ調整は必要ですが、とりあえずここで出力してみます。
スクリーンショット 2020-09-12 15.50.09

出力には出力のファイル形式(docやpdf,図など)にあわせてgtsummary内の型を変える必要があります。

スクリーンショット 2020-09-12 15.27.50
公式サイトより引用

出力の形式はgt,kable,flextableなどがあります。
デフォルトはgtとなっています。
gtとはRだけでいい感じの表を作るためのgtパッケージの形式です。
今まで見てきた表です。確かにとても綺麗です。




HTMLの場合

HTML形式の場合はこのまま出力できます。
右下のViwerからSave as Web Pageを選択します。
Save as ImageやCopy to Clipboardは日本語が入るとうまく機能しません。
画像として残したいならパソコンでスクリーンショットを撮ったほうが早いです。
スクリーンショット 2020-09-12 16.07.28


Office(Word, Power Point)の場合

Officeの場合はflextable形式に変換します。
先程のコードにas_flex_table()を加えます。
そしてtable1(他の名前でもいい)という変数を付けます。
今までずっとtheme_gtsummary_compact()を設定していたので小さなサイズの表になっています。
今回は一度リセットし再度日本語にした状態で行ってみます。

reset_gtsummary_theme()
theme_gtsummary_language("ja")
table1 <- data %>% tbl_summary(by = 治療) %>% add_p() %>% add_overall() %>% bold_labels() %>% as_flex_table() table1
スクリーンショット 2020-09-12 16.18.20
gt形式と雰囲気が変わりました。こちらのほうが論文で見た形に近いかもしれません。

次にこれをWordやPower Pointに変換しますが1行のコードでできてしまいます。
手っ取り早く変換するにはprint(今作った変数名, preview="○○")を使います。
preview="○○"はpreview="docx"preview="pptx"の2択です。

#Wordで開く場合
print(table1, preview = "docx")
スクリーンショット 2020-09-12 16.22.24


#Power Pointで開く場合
print(table1, preview = "pptx")
スクリーンショット 2020-09-12 16.26.25


どうでしょうか!!!
この後の文字サイズや修正もWordやPower Pointで直接できてしまいます!
これだけでもかなり時短になるかと思います。

ただ集計方法を平均(標準偏差)やt検定を使うためにはもう少しコードを変える必要があります。


11.様々な箇所を修正する
ここからは応用編になりますが、gtsummaryは色々な箇所を変更することができます。

スクリーンショット 2020-09-27 20.54.57


変数の順番を変えたい、いらない変数を減らしたい:gtsummaryの前に修正しておく
列内の文字の修正(0→対照群)|:gtsummaryの前に修正しておく
集計の方法や表示を修正したい:tbl_summary()
検定の方法を変えたい:add_p()
左上のCharacteristicを変えたい:modify_header(label=)
ヘッダーを編集したい:modify_header(stat_by=)
overallの列を編集したい:add_overall()
ヘッダーの上にもう一つ追加したい:modify_spanning_eader()|
Nの列を編集したい:add_n()
脚注を編集したい:modify_footnote()
タイトル、サブタイトルを加えたい:gtsummaryで直接できない。(gtパッケージやflextableパッケージで編集する)

上記の()内にコードを追加することで細かい修正ができます。


12.変数を並べ替えたい

変数はdataの変数の並んでいる順になっています。
そのため並べ替えるにはgtsummaryを使う前にdataの並び順を変える必要があります。
dataを並べ替えるには【2-4】Rで指定した列や行だけを取り出すselect関数、slice関数、filter関数を紹介しますで紹介したtidyverseパッケージselect関数を使います。
data %>% select関数で並べ替える %>% tbl_summary() %>% ・・・とtbl_summary()の前に入れます。

data %>% 
  select(治療,年齢,体重,性別,術側,MMT) %>% 
  tbl_summary(by = 治療) %>% 
  add_p() %>% 
  add_overall() %>%   
  bold_labels()
スクリーンショット 2020-09-12 16.38.50


13.列内の文字の修正(0→対照群など)

今は治療が0,1となっています。
修正するにはgtsummaryの前にdataのデータ自体を修正します。
列を修正するには【2-2】Rのmutate関数を使って列の追加や修正を行うで紹介したmutate関数を使います。先ほどと同じでtbl_sumamry()の前に入れます。
mutate()の中はfactor関数を使います。factor(列名, labels=c(○○))に関しては【1-10】Rでよく使われる型について説明します。のfactor関数の箇所をご参照ください。

data %>% 
  select(治療,年齢,体重,性別,術側,MMT) %>% 
  mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
  tbl_summary(by = 治療) %>% 
  add_p() %>% 
  add_overall() %>%   
  bold_labels()
スクリーンショット 2020-09-12 17.46.55



14.集計を変更する(tbl_summaryのいろいろな設定)

集計方法を指定するのはtbl_summary()です。
公式ドキュメントはこちら



既にby=は説明していますが、ここではtype(変数のタイプ)、label(一番左の変数名)、statistics(集計方法)、digits(小数点第何位まで示すか)を紹介します。

どれも○○ =変えたい変数名を選択 ~ 何に変えるかといった表記方法になります。

14−1.集計のタイプを指定する(type=)

今回のデータでいうと、MMTはカテゴリーとして表にする予定で特に修正は必要ありませんでした。もしMMTがカテゴリーでなく数値として認識されるとこうなります。
スクリーンショット 2020-09-12 17.05.44
これは列ごとに集計のタイプが決まっているからです。
タイプというのは以下のの3種類です。

"continuous(数値)"
"categorical(カテゴリー)"
"dichotomous(2択)" 

gtsummaryでは自動的にタイプを決めてくれるのですが、もし変更したい場合場合はtbl_summary()内でtype = 変えたい変数名 ~ "○○"を加えます。
○○は先程の3択です。

例えばMMTをカテゴリー変数として扱いたければ以下のようになります。

data %>% 
  tbl_summary(by = 治療,
              type = MMT ~ "categorical")


14−2.一番左の変数名を変更する(label=)

左の列は変数名が載っていますが、単位など付け加えたい場合もあります。
変数名を直接編集したい場合はtbl_summary()のlabel=""で編集します。
もし変えたい変数が2つ以上の場合はlist()を使う必要があります。
下の2つのコードはlabel=しか変わっていないので見比べてください。上のツイートも参考になります。

#変えたい列が1つだけのとき
data %>% mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% select(治療,年齢,体重,性別,術側,MMT) %>% tbl_summary(by = 治療, label = 年齢 ~ "年齢(歳)") %>% add_p() %>% add_overall() %>% bold_labels()

#変えたい列が2つ以上の時はlist()の()内に入れる
data %>% mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% select(治療,年齢,体重,性別,術側,MMT) %>% tbl_summary(by = 治療, label = list(年齢 ~ "年齢(歳)", 体重 ~ "体重(kg)")) %>% add_p() %>% add_overall() %>% bold_labels()
スクリーンショット 2020-09-12 18.06.04


14−3.集計方法を変える(statistic=)
集計方法のデフォルトが中央値(四分位範囲)となっています。
最新バージョン(バージョン1.3.4.9008以降)ではtheme_gtsummary_mean_sd()関数を使えば平均(標準偏差)になりますが、平均±標準偏差平均と中央値を同時使いしたい場合はtbl_summary()内でstatistic=○○を指定します。
tbl_summary(statistic = 変えたい変数名 ~ "集計方法")という書き方をします。

data %>% 
  select(治療,年齢,体重,性別,術側,MMT) %>% 
  mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
  tbl_summary(by = 治療,
              statistic = c(年齢, 体重) ~ "{mean}({sd})") %>% 
  add_p() %>% 
  add_overall() %>%   
  bold_labels()
スクリーンショット 2020-09-12 18.11.46
脚注を見ると平均(標準偏差)に変わっています。


今回はc(年齢, 体重)と列名を指定しました。
もし全ての数値の列を選びたい時はall_continuous()を使うこともできます。
data %>% 
  select(治療,年齢,体重,性別,術側,MMT) %>% 
  mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
  tbl_summary(by = 治療,
              statistic = all_continuous() ~ "{mean}({sd})") %>% 
  add_p() %>% 
  add_overall() %>%   
  bold_labels()

もし年齢だけ指定すると、年齢は平均(標準偏差)、体重はデフォルトの中央値(四分位範囲)となります。

data %>% 
  select(治療,年齢,体重,性別,術側,MMT) %>% 
  mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
  tbl_summary(by = 治療,
              statistic = 年齢 ~ "{mean}({sd})") %>% 
  add_p() %>% 
  add_overall() %>%   
  bold_labels()

スクリーンショット 2020-09-12 18.19.20


もし中央値[四分位範囲]など2箇所変更するならlabel=の時と同様にlist()を使います。
data %>% 
  select(治療,年齢,体重,性別,術側,MMT) %>% 
  mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
  tbl_summary(by = 治療,
              statistic = list(年齢 ~ "{mean}({sd})",
                                 体重 ~ "{median}[{p25},{p75}]")) %>% 
  add_p() %>% 
  add_overall() %>%   
  bold_labels()
スクリーンショット 2020-09-12 18.21.16

ここで右側の説明を行います。
{集計方法}とすることで集計ができます。{}に入っていない箇所はその文字がそのまま使われます(ここでは( )や[]やカンマなど黒文字の箇所。%や人など文字を入れることも可能)。
スクリーンショット 2020-09-12 18.32.59

他にも集計方法はあります。
{n}:集計した人数
{p}:集計した人数の%
{N}:母数
{p○○}:○○%時の値(0から100の数値を指定。"{p10},{p90}"など)
{var}:分散
集計方法を変えると自動的に脚注が変わるのもありがたいところです。

もし平均±標準偏差にしたい場合はどうしたらいいでしょうか?
"{mean}±{sd}"とすればOKです。±の間にスペースを入れることだって可能です。

14−4.小数点をあわせる(digits=)
デフォルトでは整数になっていますが、小数点第何位まで指定することができます。
人数は整数(0)がいいですし、小数点第何位で指定したい場合もあります。

data %>% 
  select(治療,年齢,体重,性別,術側,MMT) %>% 
  mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
  tbl_summary(by = 治療,
              statistic = list(年齢 ~ "{mean}({sd})",
                                 体重 ~ "{median}[{p25},{p75}]"),
              digits = list(c(年齢,体重) ~ c(1,2), 
                            c(性別,術側,MMT)  ~ c(0,1))) %>% 
  add_p() %>% 
  add_overall() %>%   
  bold_labels()
スクリーンショット 2020-09-12 18.48.13


digits = list(c(年齢,体重) ~ c(1,2), 
                  c(性別,術側,MMT)  ~ c(0,1)))
~の左側は今までと同様に列の選択を行っています。
今回~の右側は2つ数値を指定しています。
これは{平均}({標準偏差})という2つの数値を使っているからです。
年齢と体重の1つめの1は平均の小数点は第1位、2つめの2は小数点第2位という意味になります。

ここで体重の75%を確認してみてください。小数点第1位になっています。これは体重は{中央値}、{25%}、{75%}の3つの数値を指定する必要があるのに3つ目を指定していないためです。
そのため正しくは以下になります。


data %>% 
  select(治療,年齢,体重,性別,術側,MMT) %>% 
  mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
  tbl_summary(by = 治療,
              statistic = list(年齢 ~ "{mean}({sd})",
                                 体重 ~ "{median}[{p25},{p75}]"),
              digits = list(年齢 ~ c(1,2), 
                            体重 ~ c(1,2,2),
                            c(性別,術側,MMT)  ~ c(0,1))) %>% 
  add_p() %>% 
  add_overall() %>%   
  bold_labels()



15.add_p()のいろいろな設定

14ではtbl_summary()のいろいろな設定について説明しましたが、今度はadd_p()の設定について説明します。

スクリーンショット 2020-09-27 20.54.57


15−1.検定方法を変える

集計のデフォルトは中央値(四分位範囲)でした。
統計を見るとデフォルトはWilcoxonの順位和検定(マン・ホットニーのU検定)やクラスカル・ウォリス検定といったノンパラメトリック検定となっています。

これをt検定や分散分析で行うには修正が必要です。
手っ取り早いのは最新バージョンで使えるtheme_gtsummary_mean_sd()です。
(上にある5−3をご参照ください)
これを使えば集計も検定も自動的にパラメトリックに変わります。

theme_gtsummary_mean_sd()


もし個別に修正するにはadd_p(列名 ~ "検定方法")と記載します。

で指定します。2つ移乗あればlist()を使います。このあたりの記載は14集計を変更するで説明していますのでさかのぼって確認してみてください。
"t.test" for a t-test,(t検定)
"aov" for a one-way ANOVA test,(1元配置分散分析)
"wilcox.test" for a Wilcoxon rank-sum test,(マン・ホットニーのU検定)
"kruskal.test" for a Kruskal-Wallis rank-sum test,(クラスカル・ウォリス検定)
"chisq.test" for a chi-squared test of independence,(χ二乗検定)
"chisq.test.no.correct" for a chi-squared test of independence without continuity correction,
"fisher.test" for a Fisher's exact test,(フィッシャーの正確確率検定)
"lme4" for a random intercept logistic regression model to account for clustered data, lme4::glmer(by ~ variable + (1 | group), family = binomial). The by argument must be binary for this option.(階層モデルのロジスティック回帰?)

Tests default to "kruskal.test" for continuous variables, "chisq.test" for categorical variables with all expected cell counts >=5, and "fisher.test" for categorical variables with any expected cell count <5. A custom test function can be added for all or some variables. See below for an example.
公式サイトより引用


例えば全ての数値データをWilcoxonの順位和検定でなくt検定にするには add_p(all_continuous() ~ "t.test")となります。

またt検定とWilcoxonの順位和検定を同時使いしたければlistを使って以下のようになります。

data %>% 
  select(治療,年齢,体重,性別,術側,MMT) %>% 
  mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
  tbl_summary(by = 治療,
              statistic = list(年齢 ~ "{mean}({sd})",
                                 体重 ~ "{median}[{p25},{p75}]"),
              digits = list(年齢 ~ c(1,2), 
                              体重 ~ c(1,2,2),
                              c(性別,術側,MMT)  ~ c(0,1))) %>% 
  add_p(list(年齢 ~ "t.test",
                  体重 ~ "wilcox.test")) %>% 
  add_overall() %>%   
  bold_labels()

スクリーンショット 2020-09-27 12.00.58


16.左上の変数(Characteristic)を修正したい
スクリーンショット 2020-09-27 20.54.57


左上の「変数(Characteristic)」と書かれている所を編集するにはmodify_header(label="○○")を追加します。

data %>% 
  tbl_summary(by = 治療) %>% 
  add_p() %>% 
  add_overall() %>%   
  modify_header(label = "○○") %>% 
  bold_labels()
スクリーンショット 2020-09-27 21.52.30


17.ヘッダーを編集したい

上の表では対照群, N = 53となっています。
スクリーンショット 2020-09-27 20.54.57

オレンジの箇所をヘッダーといい、変えたり改行するにはmodify_header(stat_by=)で編集します。
<br>は改行記号で、{style_percent(p)}は人数を%で表記することができます。


スクリーンショット 2020-09-27 23.24.04

このような修正も可能です。
スクリーンショット 2020-09-27 23.30.18

一つ前の左上の「変数(Characteristic)」を変えるのもmodify_header(label=)でした。
同じ関数なので両方とも修正したい場合はまとめることができます。

  data %>% 
    mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
    select(治療,年齢,体重,性別,術側,MMT) %>% 
    tbl_summary(by = 治療) %>% 
    add_p() %>% 
    add_overall() %>%   
    modify_header(label = "○○",
                  stat_by = "**{level}**<br>N =  {n}") %>% 
    bold_labels()
スクリーンショット 2020-09-27 23.25.37
よくみると全体(Overall)は改行されていません。これはOverallを修正するのはadd_overall()内になるからです。

18.全体(Overall)の列を修正したい

スクリーンショット 2020-09-27 20.54.57

全体(Overall)に関してはadd_overall()で修正を行います。
先程modify_header()でヘッダーを2段にしました。
ここではOverallのヘッダーも2段にします。
Overallのヘッダーはadd_overall(col_label = "")を使います。

スクリーンショット 2020-09-27 23.47.33

  data %>% 
    mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
    select(治療,年齢,体重,性別,術側,MMT) %>% 
    tbl_summary(by = 治療) %>% 
    add_p() %>% 
    add_overall(col_label = "**全体**<br>N = {N}") %>%   
    modify_header(label = "○○",
                  stat_by = "**{level}**<br>N =  {n}") %>% 
    bold_labels()
スクリーンショット 2020-09-27 23.49.35

19.ヘッダーの上に列を追加したい

スクリーンショット 2020-09-27 20.54.57

ピンクの箇所になりますが、ヘッダーの上に「治療」などもう一列つけることもできます。
modify_spanning_header(starts_with("stat_") ~ "**○○**") とします。
starts_with("stat_")は集計をしている列という意味になります。
**で挟まれたところは太字になります。

スクリーンショット 2020-09-27 23.55.35


20.脚注を編集したい

スクリーンショット 2020-09-27 20.54.57

脚注(黄色)の箇所はmodify_footnote()で編集します。

脚注自体をなくす
modify_footnote(update = everything() ~ NA)を追加。

集計方法(1となっている列)を編集
modify_footnote(update = starts_with("stat_") ~ "○○")を追加。

検定(2となっている列)を編集
modify_footnote(update = starts_with("p") ~ "○○")を追加。

data %>% 
  mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
  select(治療,年齢,体重,性別,術側,MMT) %>% 
  tbl_summary(by = 治療) %>% 
  add_p() %>% 
  add_overall() %>% 
  modify_footnote(update = starts_with("stat_") ~ "11111") %>% 
  modify_footnote(update = starts_with("p") ~ "22222") %>% 
  bold_labels()
スクリーンショット 2020-09-28 0.33.37

21.タイトル・サブタイトルを加える
表のタイトル・サブタイトルは実はgtsummaryで作成できません。

もう一度出力のイメージを確認します。
今までは全て「集計」と書かれたところで作業していました。
ただタイトル・サブタイトルはas_○○と出力に合わせた形式に変換してから作成します。

スクリーンショット 2020-09-12 15.50.09

21−1.HTML形式でタイトル・サブタイトルを加える

今回はas_tg()形式で行います。
もともとgtsummaryのデフォルトはas_gt()となっているので修正はいりませんがわかりやすく加えています。

table1 <-
  data %>% 
  mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
  select(治療,年齢,体重,性別,術側,MMT) %>% 
  tbl_summary(by = 治療) %>% 
  add_p() %>% 
  add_overall() %>%   
  bold_labels() %>% 
  as_gt()
次にgtパッケージを使って編集します。gtパッケージはgtsummaryの元になっているパッケージです。
まずはgtパッケージを呼び出し、tab_header()を使います。
 library(gt)
#タイトル・サブタイトルを加えるにはtab_header()を用いる
table1 <- table1 %>% tab_header(title = "タイトル", subtitle = "サブタイトル")
スクリーンショット 2020-09-28 0.50.42


21−2.office形式でタイトル・サブタイトルを加える

今回はflextable形式で行います。
flextable形式にするにはas_flex_table()を使います。
table1 <-
  data %>% 
  mutate(治療 = factor(治療, labels = c("対照群", "治療群"))) %>% 
  select(治療,年齢,体重,性別,術側,MMT) %>% 
  tbl_summary(by = 治療) %>% 
  add_p() %>% 
  add_overall() %>%   
  bold_labels() %>% 
  as_flex_table()
table1
次にflextableパッケージを使って編集します。
まずはflextableパッケージを呼び出すのですが、このパッケージにはタイトル・サブタイトルを指定する関数がありません。代わりにadd_header_lines()を使って1番上に1行増やします。
ポイントはサブタイトル→タイトルの順に作ることです。
フォントサイズを変更することもできます。
 library(flextable)
table1 <- table1 %>% add_header_lines("サブタイトル") %>% add_header_lines("タイトル") %>% fontsize(i = c(1,2), size = c(20,14), part = "header")
table1
スクリーンショット 2020-09-28 0.58.52
この後はprint関数でWordやPower Pointに変換します。

#Wordの場合
print(table1, preview = "docx")
#Power Pointの場合
print(table1, preview = "pptx")


22.まとめ

かなり長くなりましたが、gtsummaryでtable1を作成する方法を紹介しました。

ただなれないうちは全てgtsummaryで作成しなくても、ある程度作成したら残りはdocxやpptxに変換して修正できることもできます。。

まずは触ってみながらどこまでをRでするか、どこからをOfficeでするかなど自分なりの方法を探してみるのもいかがでしょうか?

そしてgtsummaryは回帰分析の結果など、table1以外の表にも対応しています。
そこに関しては今後の記事で紹介します。



前回の記事では相関係数と相関係数の検定を行う方法を紹介しました。


ただ複数の列がある場合1つずつ相関係数を調べるのは大変な作業です。
今回は相関行列について紹介します。


1.使うデータ

今回はこちらのダミーデータを使います。

歩行が自立しているか(0:非自立、1:自立)と各種データが入っています。
スクリーンショット 2020-07-15 5.42.46



性別が男性・女性になっています。

url <- "https://raw.githubusercontent.com/mitti1210/myblog/master/data01.csv"
dat <- read.csv(url)
dat$氏名 <- NULL  #氏名の列を削除


2.GGallyパッケージのggpairs関数

相関係数を手早く確認する方法としてGGallyパッケージggpairs関数があります。
この方法のメリットは相関係数と散布図が同時に見れサクッとわかるです。

library(tidyverse)
library(GGally)
ggpairs(dat)
スクリーンショット 2020-07-15 9.30.36


左下に散布図、右上には相関係数とまとめて表示されます。
数値と数値の場合は散布図、カテゴリーと数値の場合はヒストグラムと箱ひげ図が表示されます。

もしmacで日本語が□□と文字化けしていたらもう1行追加します。

ggpairs(dat)+
  theme_gray(base_family = "HiraKakuPro-W3")
スクリーンショット 2020-07-15 9.30.47



そしてggpairsの便利な機能としてはあるカテゴリーごとに色を付けることができます。
aes_string(colour="カテゴリーの列", alpha=0.5)を追加します。
カテゴリーの列はcharactor型、もしくはfactor型である必要があります(数値型だとダメ)
alphaは色の透過の程度なので0.7あたりでも大丈夫です。
ggpairs(dat,aes_string(colour="性別", alpha=0.5))
スクリーンショット 2020-07-15 9.29.32





3.cor関数とcorrplotパッケージを使う

次にcor関数corrplotパッケージで相関行列を作成します。
この方法のメリットはp値が計算できる、pythonっぽい相関行列のグラフが作れるです。

散布図行列を求めるにはcor関数を使います。

ただ男性・女性のような文字は対応できないので0,1に置き換えます。
dat$性別 <- ifelse(dat$性別 == "男性", 0, 1)

c <- cor(dat)
c
スクリーンショット 2020-07-15 9.42.04


四捨五入する時はround関数を使います。
round(c,2)
スクリーンショット 2020-07-15 9.42.38


もしスピアマンの相関係数を求めるのであればmethod="spearman"を追加します。
c <- cor(dat, method = "spearman")

ただ数値の羅列は中々見にくいものです。そこでcorrplotパッケージのcorrplot関数を使いグラフにします。
corrplot(c, method = "color", addCoef.col = TRUE)
method="color"はグラフの見た目、addCoef.col = TRUEはグラフに数値を入れるオプションなのでそのままコピペで大丈夫です。変更するのはcのところで、ここには先程もとめた相関行列を指定します。

c <- cor(dat, method = "spearman")
cの部分です。もちろん違う名前でも大丈夫です。

もしmacで日本語文字化けの問題があれば1行追加します。
par()関数はフォントの指定ができます。

library(corrplot)
par(family="HiraKakuPro-W3") #macで日本語が文字化けするときに追加。ここでフォントの指定ができる
corrplot(c, method = "color", addCoef.col = TRUE)
スクリーンショット 2020-07-15 9.47.46
これでかなり見やすくなりました。pythonっぽいです。


検定をしてp値を求める

p値を求めるのであればcorrplotパッケージのcor.mtest()関数を使います。
cor.mtestはp値、95%信頼区間の下限、上限の3つが表示されます。スピアマンの相関係数の場合はmethod = "spearman"を追加するだけなのですが、同値(タイ)があると95%信頼区間が出ないようです。

p <- cor.mtest(dat)
p
スクリーンショット 2020-07-15 7.20.20

p$p:p値
p$lowCI:95%信頼区間の下限
p$uppCI:95%信頼区間の上限

先程の相関行列にp<0.05かそうでないかを出すこともできます。
p.mat:p値の行列(ここではp$p。上の図をよく見ると$pと書いてあります。)
sig.level = ◯◯(0.05 にしたり、0.2に設定することもあります)

corrplot(c, method = "color", addCoef.col = TRUE, p.mat = p$p, sig.level = 0.05)
スクリーンショット 2020-07-15 10.33.07

これで有意差があるものがどれか一目瞭然です。


4.corrrパッケージを使う

もう一つcorrrパッケージを紹介します。
この方法のメリットは相関が高い組み合わせ順が知りたい!というときに便利です。

相関行列を出す時はcor関数でなくcorrelate関数を使います。
使い方はcor関数と同じです。cor→correlateに変わるだけです。
スピアマンの相関係数にしたい時はmethod="spearman"を追加すればOKです。

c <- correlate(dat, method = "spearman")
c
スクリーンショット 2020-07-15 10.41.13


マイナスの値が赤字になっただけでここまでは今までと違いはありませんが、corrrの便利な所はここからになります。

まずcorrr関数は第2章で扱ったtidyverseと同じ%>%が使えます。

右上と左下は同じものなので片方を削除することができます。削除するにはshave関数を使います。
c %>% shave()
スクリーンショット 2020-07-15 10.45.16



これで左下だけになりました。

次に行列形式ではなくlong形式に変換します。streach関数を使います。
streach(na.rm=TRUE)とするとNAとなっているいらない列を消すことができます。

c %>% shave() %>% stretch(na.rm = TRUE)
スクリーンショット 2020-07-15 10.55.02



次に並べ替えます。arrange関数を使います。何も指定しなかければ相関係数が高い順から並びます。
c %>% shave() %>% stretch(na.rm = TRUE) %>% arrange()
スクリーンショット 2020-07-15 10.55.14
スクリーンショット 2020-07-15 10.55.31


ただ注意しないといけないのが相関係数は-1〜1の値を取ります。
相関がないのは0付近で、-1は1と同様に強い相関があります。
なので相関が強い順に並べるには相関係数の絶対値で並び替える事が必要になります。
相関係数の列名はrになっています。abs関数が絶対値を求める関数なのでabs(r)となります。

c %>% shave() %>% stretch(na.rm = TRUE) %>% arrange(abs(r))
スクリーンショット 2020-07-15 10.55.31


そしてこれだと相関係数が0に近い方から並んでしまったので順番を逆にします。
順番を逆にするにはdesc(並べ替えたい列)を使います。
並べ替えたい列はabs(r)でした。desc(abs(r))とします。
()が多くなってエラーが出やすくなるので注意してください。

c %>% shave() %>% stretch(na.rm = TRUE) %>% arrange(desc(abs(r)))
スクリーンショット 2020-07-15 10.55.45

相関行列を1個ずつ見てると見落としがちですが、これなら見落とすことはありません。


5.相関係数の有意差がないってどういう意味?

スクリーンショット 2020-07-15 10.33.07

ところで×がついた項目(p値が有意水準に満たさなかった)というのはどういう意味でしょうか?
これは相関係数が低いから有意差が出ないわけではありません。

相関係数の検定は「相関係数が0である」と仮説を立てます。
そしてデータから信頼区間を求め、信頼区間に0が含まれているかを調べます。
その時もし信頼区間に0が含まれていたら有意差なしとなります。

例えば上の図でいうと性別×BBSを見ると相関係数が-0.13となっています。
ただ95%信頼区間を見ると-0.27〜0.007となっています。
つまり今回は計算上-0.13となったけど実は信頼区間にマイナスの値も0もプラスの値も入っていて、符号が変わる可能性がある→正の相関があるとも負の相関があるとも言えないよね。といったニュアンスになります。もし95%信頼区間が-0.27〜-0.001であれば符号が変わらないため有意差あり(p<0.05)となります。ただ相関係数0.13というのはほとんど相関がないという意味ではありますが...

つまり相関係数の検定のp値は出した相関係数の符号が変わる可能性が高いのか低いのか?を見る検定と考え、p値が低いから相関係数が強いんだという意味ではないことに注意が必要です。相関の強さはあくまでもp値ではなく相関係数の値で評価します。



6.まとめ

今回は相関行列を紹介しました。

相関行列を出すだけならEZRでもできますが、検定やグラフの作成はEZRだけではできません。

他にもいろいろな方法はあるのですが、今回は3パターン紹介しました。

ただGGallyやcorrrはもっといろいろなことができます。

もっと詳しく知りたい方は以下のサイトがおすすめです!








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


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


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


 


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


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


今回は相関係数(ピアソン・スピアマン)を紹介します



まず相関係数についてはハルさんのサイトをご参照ください。



また1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。
もし1〜4まででわからない部分があれば確認してみてください。

 

1.準備

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

 

ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。
これを前提に次から進めていきます。
スクリーンショット 2019-10-20 7.54.14


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

次にRのコードを書くためのスクリプトファイルを作ります。
ここでは【4-15】相関係数.Rとしておきます。

スクリーンショット 2020-07-14 23.41.44


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


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

デモデータ(pearsonの相関係数)

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


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


実際のコードは以下になります。
前回のコードのURL(" "の中)とdestfileのdata/以降を変更するだけでOKです。
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/06/demo-pearson.xlsx"
destfile = "data/demo-pearson.xlsx"
download.file(url, destfile)
スクリーンショット 2019-11-23 22.14.17


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


4.データの読み込み

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

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

スクリーンショット 2019-12-08 20.31.39


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




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



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

Name:ハルさんのサイトと同じでborgとします(大文字・小文字は別物とされます)
Sheet:このExcelは1つしかデータがないのでDefaultのままでOK
Range:先ほど確認したB2:D32(実は今回のように表以外に余計なものがなければ空欄でもOK)

スクリーンショット 2020-07-14 23.53.34



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

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



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

5.散布図で確認

相関係数を見るときはいきなり相関係数を求めるのではなく、先に散布図を書くほうが事故が減ります。
【3-3】Rのggplot2で散布図を作るgeom_point関数で散布図を書く方法を紹介しています。
散布図を書くのはgeom_point関数、回帰直線を書くのはgeom_smooth関数です。
geom_smooth関数では(method="lm")を必ず付ける必要があります。
これは直線ですよという意味です。直線以外にも色々あるんです。
ちなみに色を付けるときはcolor=を、se=FALSEは今回はあってもなくても構いません。

library(tidyverse)
ggplot(data=borg,aes(x=`6MWD`,y=VC))+
  geom_point()+
  geom_smooth(method = "lm", color="red", se=FALSE)
注意!
列名や変数名の最初の文字が数字の時(今回でいうと6MWD)の場合、x=6MWDのようにするとエラーが出ます。こういった場合は` `(shiftキーを押しながら@マークを押す)でエラーを回避できます。でもなるべくなら最初の文字は数字じゃない方が後々便利です。


スクリーンショット 2020-07-15 0.06.04


他にも「もっとかんたんに散布図だけみたい」というときはplot関数があります。

plot(borg$`6MWD`, borg$VC)
ここでも`6MWD`と` `を忘れないようにします。
スクリーンショット 2020-07-15 0.05.51

こうみると6MWDとVCは直線上の関係にあるのが見て取れます。
こういった場合はピアソンの相関係数が使えます。
逆に外れ値が多いデータや直線でないデータでは相関係数は使えなくなります。

外れ値がない場合
スクリーンショット 2020-07-15 0.34.07

外れ値が5個あった場合
スクリーンショット 2020-07-15 0.33.58


そして二次関数のようなデータはうまく行きません。他にも「10代〜30代は良好、40〜50代は不良、60代以降は良好」みたいなデータもうまくいきません。
(相関係数は「一方が上がる」するともう一方が「上がるor下がる」ということだけで、2時間数のように途中で傾向が変わるデータには対処できない)
スクリーンショット 2020-07-15 0.36.06
上記のような例は相関係数を出すだけではわかりません。散布図で確認して初めて分かることなので、相関係数は「まずは散布図!」とイメージしておいてください。


6.相関係数の検定の実施

相関係数の検定は1行で終わりますが主に2種類あることを知る必要があります。

パラメトリック:ピアソンの相関係数
ノンパラメトリック、外れ値が多い:スピアマンの相関係数



どちらも相関係数の検定はcor.test(1つ目の列, 2つ目の列)です。

cor.test(borg$`6MWD`,borg$VC)
スクリーンショット 2020-07-15 1.02.00

もしスピアマンの相関係数を求めるときはmethod = "spearman"を追記します。
cor.test(borg$`6MWD`,borg$VC, method = "spearman")
スクリーンショット 2020-07-15 4.00.34


【4-2】RでMann-Whitney U 検定を行う方法でも紹介しましたが、ノンパラメトリックの手法では列の中に同じ数値があると警告がでます。

スクリーンショット 2020-07-15 4.07.07
このままでも大丈夫ではありますが(EZRも同じ字問題が起こっているが警告が出ないようにしている)、もし「警告をどうにかしたい!」とあれば、coinパッケージのspearman_test関数を使うという方法もあります(そこまでしている人はいるのでしょうか?)

spearman_test(1つ目の列 ~ 2つ目の列, data=データ)という書き方になります。

library(coin)
coin::spearman_test(`6MWD`~VC, data=borg)

cor(borg$`6MWD`,borg$VC, method = "spearman")
スクリーンショット 2020-07-15 4.17.43


ただspearman_test関数は検定はしてくれますが、相関係数自体を出してくれません。
そのためcor関数で相関係数を求める必要があります。


7.まとめ

今回は相関係数の求め方を紹介しました。

今回紹介していませんが、ハルさんのサイトにもあるように相関関係と因果関係は別物という考えはとても大切です。

そして、相関係数を出す前に散布図を作るという習慣をつけることが大切です。

散布図を出した上で直線的な関係性がありそうであれば相関係数を求めてください。

【1-3】Rのソフトはどれを使えばいいか?目的別チャートを作りましたではRを使うためのソフトを紹介しました。


その中で統計を勉強してみたいという方はRstudio + EZR(Rコマンダー)の組み合わせをおすすめしました。

Rを使っている方は「EZR(Rコマンダー)しか使っていない」や「Rstudioしか使っていない」という方が多いのですが、実はRstudioからRコマンダーやEZRを簡単に呼び出すことができます。

今回はRstudioからEZRを呼び出す方法を紹介します。


1.EZRのインストールをRstudioだけで行う。

EZRのダウンロード、インストールは自治医科大学附属さいたま医療センターのページで紹介されていますが、Rstudioからでも行えます。

①RstudioのInstallボタンを押す
Rstudioの真ん中あたりにInstallボタンがあります。これはRの様々なパッケージをプログラムを書かずにダウンロード、インストールすることができます。
スクリーンショット 2020-04-24 19.34.48

②RコマンダーとEZRをインストール
次にRコマンダーとEZRをダウンロードします。

Packagesという欄にrcmdrと入力します。
そうすると候補にRcmdrとRcmdrPlugin.EZRという2つが見つかります。
スクリーンショット 2020-04-24 19.34.58

まずはRcmdrをインストールし、次にRcmdrPlugin.EZRをインストールします。

この作業は今回やってしまえば次回以降は行う必要はありません。


2.EZRを呼び出す

EZRを呼び出すには2つ手順が必要ですが慣れるとすぐにできるようになります。

①Rコマンダーを呼び出す(2種類の方法があるので好きな方でいいです)
②RコマンダーからEZRを呼び出す


①Rコマンダーを呼び出す

1つめの方法はPackagesからRcmdrをクリックして呼び出す方法です。
スクリーンショット 2020-04-24 19.35.03
2つめの方法は下のコードを入力する方法です
library(Rcmdr)

慣れるとコードを書いた方が早いですが、最初はクリックする方法の方が早いかもしれません。
(自分がそうでした)


②RコマンダーからEZRを呼び出す

EZRはRコマンダーから呼び出します。

ツール→Rcmdrのプラグインをロード…を選択
スクリーンショット 2020-04-24 19.35.13


RcmdrPlguin.EZRが選択されていることを確認してOK
スクリーンショット 2020-04-24 19.35.21

OKを押して再起動
スクリーンショット 2020-04-24 19.35.28

これで完了です。
見た目は似ていますが「標準メニュー」となっていればEZRになっています。
スクリーンショット 2020-04-24 19.35.36



3.注意点
・RstudioのPackageにRcmdrPlugin.EZRがあるのですが、クリックしても何も反応しません。上記の手順が必要です。

・EZRを使うと結果はRstudioのコンソール画面に出てきます

・1度開いたEZRを×で閉じてしまうと、同じ方法を行っても開かないことがあります。一旦Rstudioを再起動するかプロジェクトを開き直すとまたできるようになります。

・グラフはRstudioのPlotではなく別ウィンドウで開くことがあります。

この辺りはパソコンのOSでも挙動が違うかもしれません。私はmacを使っていますが、もしWindowsで違う挙動があればこの記事でも紹介しますので連絡いただけると嬉しいです。


4.まとめ
今回はRstudioとEZRを同時使いするための方法を紹介しました。

Rstudioだとどう使えばいいかわからない。でもいずれ使えるようになってみたい。という方は同時使いをがおすすめです。EZRの手軽さ、それがどのようなコードで書かれているのかを同時に体験することができると思います。

またEZRの使い方はハルさんのシロート統計学講座がおすすめです!
無料ですが、これだけでEZRの使い方はほとんど理解できると思います。






データを扱う中で「外れ値を外す方法」に関して検討してみます。

apply関数の記事で以下のコメントがありました。




大変勉強になりました。
列ごとに、ファイルの値(NAがある)が2SDより高い値をNAにしたいです。
Data01 <-read_excel("raw.xlsx",1)
Data02 <-apply(Data01, 2, sd, na.rm = TRUE)
Data01[Data01 > 2*Data02]<- NA
以上のコードをして最後は列ごとにうまく出来なかったです。
ご提言があれば大変嬉しいです。

自分もそうですが、実際にやってみるとうまくいかないことってたくさんあります。
今回は質問を一緒に考えながらどう考えるか?どんな解決方法があるかを考えて、今までの理解を深めるきっかけになればと思います。

1.データの作成

コードを確認してみます。
Data01 <-read_excel("raw.xlsx",1)
これはraw.xlsxの1シート目を読み込むという意味なので、今回は仮のデータを作ります。

set.seed(2020)
a <- c(rnorm(998,10,1),NA,30)
b <- c(rnorm(998,20,1),NA,40)
c <- c(rnorm(998,30,1),NA,50)
Data01 <- data.frame(a,b,c)
head(Data01)
簡単に言うと以下のとおりです。
a : 平均10, 標準偏差1のデータ998個にNAと30を加えた1000個のデータ
b : 平均20, 標準偏差1のデータ998個にNAと40を加えた1000個のデータ
c : 平均30, 標準偏差1のデータ998個にNAと50を加えた1000個のデータ

スクリーンショット 2020-03-28 2.34.45
この赤色が2SDを超えたデータです。やりたいことはこの赤い点のデータを削除することです(これの是非については今回問いません)

2.apply関数のおさらい

次の行を確認します。
Data02 <-apply(Data01, 2, sd, na.rm = TRUE)
apply関数【2-1】Rのfor関数、apply関数を使ってまとめて標準偏差などの統計量を求める方法で紹介しました。下はイメージ図です。
スクリーンショット 2020-03-27 13.46.36
上図ではsd関数(標準偏差)を求めました。
また欠損値(NA)があるのでrm.na = TRUEをつけることでNAの値を外して計算してくれます。

スクリーンショット 2020-03-27 13.54.01
それぞれの列の標準偏差が表示されました。
ただ元データが2標準偏差を超えているかどうかを確認するにはこれだけでは足りません。


3.どう考えるか?
スクリーンショット 2020-03-28 4.00.30

それぞれのデータが2SDを超えているかどうかを判定する必要があります。
具体的には(平均-2SD) < それぞれの値 < (平均 + 2SD)の条件を満たしているかどうかとなります。

そのため以下の2つを行えばなんとかできそうです。
①a,b,cそれぞれの平均±2標準偏差の値を求める
②条件に合えばそのまま、条件にあわなければNAと変換する


 4.平均±2SDを求める

まずa,b,cの平均と標準偏差を求めます。
これには前述のapply関数かsummarize関数を使います。
まずはapply関数でやってみます。

Data01_mean <-apply(Data01, 2, mean, na.rm = TRUE)
Data01_sd <-apply(Data01, 2, sd, na.rm = TRUE)
Data01_mean
Data01_sd

low <- Data01_mean - 2*Data01_sd
high <- Data01_mean + 2*Data01_sd
low
high
スクリーンショット 2020-03-28 3.17.41

ここからlow[2]、もしくはlow["b"]とすればbの平均-2SDが取り出せます。
high[1]、もしくはhigh["a"]であればaの平均+2SDとなります。

5.条件に合わないものをNAにする

次に条件に合わないものを抽出しNAにします。
データを抽出するにはデータ[条件]を使います

(aのデータ < 平均-2sd) または (平均+2SD < aのデータ)だとこうなります
Data01$a[(Data01$a < low["a"]) | (high["a"] < Data01$a)]
スクリーンショット 2020-03-28 4.08.45

条件式の記号については以下をご参照ください。今回はまたはの意味の を使いました。
スクリーンショット 2019-02-11 2.04.40
抽出した条件に <- を使うことでデータを置き換えることができます。
Data02 <- Data01
Data02$a[(Data02$a < low["a"]) | (high["a"] < Data02$a)] <- NA
Data02$b[(Data02$b < low["b"]) | (high["b"] < Data02$b)] <- NA
Data02$c[(Data02$c < low["c"]) | (high["c"] < Data02$c)] <- NA
これで完成です!
ただ他にも方法があります。どれもよく使う方法なので、もし余力があればぜひ比べてみてください。


6.tidyverseパッケージを使う(dplyrパッケージでも同じ)

第2章で紹介しているtidyverseパッケージのsummarize関数とmutate関数を使うこともできます。
tidyverseパッケージを使うときには先にlibrary関数で呼び出します。なければ先にインストールします。
#もし初めての場合はインストール
install.packages("tidyverse")
#インストール済ならlibrary関数で呼び出す
library(tidyverse)
low_high <-
  summarize(Data01,
          a_low = mean(a, na.rm = TRUE)-2*sd(a, na.rm = TRUE),
          a_high = mean(a, na.rm = TRUE)+2*sd(a, na.rm = TRUE),
          b_low = mean(b, na.rm = TRUE)-2*sd(b, na.rm = TRUE),
          b_high = mean(b, na.rm = TRUE)+2*sd(b, na.rm = TRUE),
          c_low = mean(c, na.rm = TRUE)-2*sd(c, na.rm = TRUE),
          c_high = mean(c, na.rm = TRUE)+2*sd(c, na.rm = TRUE)) 
class(low_high) #型を確認するとdata.frameとなっている
low_high <- unlist(low_high)
class(low_high) #型を確認するとnumeric(数値型)となっている
low_high
スクリーンショット 2020-03-28 3.37.43
文字は多いですが、コピペを使えば作業自体は楽になります。
ただsummarize関数の結果は表形式(data.frame型)で、このままでは数値として計算できません。
そのためunlist関数を使って数値型に戻します。
low_high[5]、もしくはlow_high["c_low"]とすればcの平均-2SDとなります。


既に第2章を読まれた方は%>%を使う事もできます。
#%>%使う場合
low_high <- summarize(Data01, a_low = mean(a, na.rm = TRUE)-2*sd(a, na.rm = TRUE), a_high = mean(a, na.rm = TRUE)+2*sd(a, na.rm = TRUE), b_low = mean(b, na.rm = TRUE)-2*sd(b, na.rm = TRUE), b_high = mean(b, na.rm = TRUE)+2*sd(b, na.rm = TRUE), c_low = mean(c, na.rm = TRUE)-2*sd(c, na.rm = TRUE), c_high = mean(c, na.rm = TRUE)+2*sd(c, na.rm = TRUE)) %>%
unlist()
low_high
スクリーンショット 2020-03-28 3.32.12

次に条件に応じて分類するifelse関数を使うこともできます。

ifelse(条件TRUEの場合FALSEの場合)

列の追加や修正を行うmutate関数を使うことでこうなります。
Data02 <- 
  Data01 %>% 
  mutate(a = ifelse((low_high["a_low"] < a) & (a < low_high["a_high"]), a, NA),
         b = ifelse((low_high["b_low"] < b) & (b < low_high["b_high"]), b, NA),
         c = ifelse((low_high["c_low"] < c) & (c < low_high["c_high"]), c, NA))
スクリーンショット 2020-03-28 4.57.19

条件式を少し変えているので確認してみてください。今回は&を使っています。

tidyverseパッケージは現在のRを使う上でとても重要になります。
ifelse,mutate関数に関してはこちらで紹介していますので難しかったという方はご参照ください。




7.関数を自作しapply関数でまとめて処理する
他の方法としては関数を自作するという方法があります。

関数名 <- function(○) {中身のコード}

今回はdrop_2sdという名前の関数を作ります(名前は何でもいいけど文字の最初に数字を使うのは×)。


a = ifelse((aの平均-2SD < a) & (a < aの平均+2SD), a, NA),
b = ifelse((bの平均-2SD < b) & (b < bの平均+2SD), b, NA),
c = ifelse((cの平均-2SD < c) & (c < cの平均+2SD), c, NA)
 ↓
function(x){
 ifelse((xの平均-2SD < x) & (x < xの平均+2SD), x, NA)
}
と共通部分をxにします。ちなみにxでなくてもcolとか別の名前をつけても大丈夫です。

drop_2sd <- function(x){
  ifelse((mean(x, na.rm = TRUE) - 2 * sd(x, na.rm = TRUE) < x) & (x < mean(x, na.rm = TRUE) + 2 * sd(x, na.rm = TRUE)), x, NA)
}
そしてapply関数で各列にdrop_2sd関数を使います。
Data02 <- apply(Data01, 2, drop_2sd)
スクリーンショット 2020-03-28 4.50.32
2行のコードで終わりました。
関数を作るのは最初は慣れないですが、使えるようになると強力です。


8.まとめ

今回は標準偏差±2SDの除去を3つの方法で紹介しました。

本やサイトで独学すると1つの方法でしか紹介してないことが多いです。

本当は自分の知っている知識でもできるはずなのに、知らない方法を偶然読んでしまいドツボにはまることもあります。そういった紆余曲折は勉強する上で必要な苦労でもありますが、いくつかの方法を見比べることでショートカットできることもあるのではと思っています。

また自分が躓いたところや他の記事でもコメントがあれば検討していきますので、ブログのコメントやtwitterで連絡いただければ幸いです。

↑このページのトップヘ