タグ:mutate

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

*最終更新日:2020/09/30
*gtsummaryのバージョンが上がったことでインストールのコードを一部修正



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

特に医療統計では表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からインストールする必要があります。
2020/09/30にバージョン1.3.5がリリースされました。
上記のinstall.packages("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以外の表にも対応しています。
そこに関しては今後の記事で紹介します。



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

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で連絡いただければ幸いです。

(2019.09.17 更新)

前回はmutate関数を使って列の追加や修正を行いました。
そして最後に以下のコメントで終了しました。
mutate関数を使うと「FIMの合計点」とかだけでなく、年齢を「年代ごとに分類する」だったり、ある評価の「カットオフ値以上とそれ以下」に分けた列をExcelを使わずに作成することができます。
【2-2】Rのmutate関数を使って列の追加や修正を行うより)

今回は連続変数を条件に応じて複数のカテゴリーに変換していきます。

Rではいくつか方法があるのですが、今回はifelse関数、case_when関数、cut関数を紹介します。


今回使う架空のデータは年齢、性別、病名が入っています。
下のコードをまとめてスクリプトにコピペして実行してください。
library(tidyverse)
set.seed(1) 年齢 <- floor(rnorm(100,60,20)) 性別 <- sample(c("男性","女性"),100,replace = TRUE) 病名 <- sample(c("脳梗塞", "脳出血", "骨折", "靭帯損傷","心筋梗塞"), 100, replace = TRUE) data <- data.frame(年齢,性別,病名) head(data)
スクリーンショット 2019-02-10 22.48.49



今回は%>%やmutate関数を使って作業します。
これらの使い方に不安がある場合は前回の記事をご参照ください



1.ifelse関数をつかって分類する

条件から2つに分類する場合

ifelse関数は条件からTRUEとFALSEの2つに分けることができます。

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



まず例としてdata_20という変数名に以下のように分類していきます。
・20歳未満を「未成年」、20歳以上を「成人」
・性別の男性を「male」、女性を「female」

data_20 <- data %>% 
  mutate(成人 = ifelse(.$年齢 < 20 , "未成年", "成人")) %>% 
  mutate(性別 = ifelse(.$性別 == "男性", "male", "female"))
head(data_成人)

スクリーンショット 2019-02-11 0.44.03

スクリーンショット 2019-02-11 0.46.48



条件から3つ以上分類する場合

3つ以上分類するはifelse関数のFALSEにifelseを入れ込みます。

ifelse(条件, TRUE, ifelse(条件2, TRUE, ifelse(条件3, TRUE, FALSE)))

このような形を「入れ子」と呼ぶそうです。繰り返せばいくらでも条件を増やすことができます。


今度は「年代」という変数名で20歳ごとに区切ってみます。
20歳未満なら「0〜」、
そうでなければ40歳未満なら「20〜」
そうでなければ60歳未満なら「40〜」・・・

上のイメージをコードにするとこうなります。

data_age <- data %>% 
  mutate(年代 = ifelse(.$年齢 < 20,"0〜", 
                      		ifelse(.$年齢 <40 ,"20〜", 
                        		ifelse(.$年齢 < 60,"40〜",
                        			ifelse(.$年齢 < 80, "60〜",
                        				ifelse(.$年齢 < 100, "80〜", "100〜")))))) %>% 
  mutate(年代 = factor(.$年代, levels = c("0〜", "20〜", "40〜", "60〜", "80〜", "100〜"), ordered = TRUE))
head(data_age)
スクリーンショット 2019-02-12 21.22.24

 mutate(年代 = factor(.$年代, levels = c("0〜", "20〜", "40〜", "60〜", "80〜", "100〜"), ordered = TRUE)の部分はfactorの順番を並べ替えています。

factorは50音順に並ぶので、もしこれをしなかったらグラフなど作る時に「0〜, 100〜, 20〜・・・」と100歳以上の場所がずれてしまいます。それを修正するためにlevels = でfactorの並べ替えをしています。
参照:【1-10】Rでよく使われる型について説明します。


個人的に感じる「ifelse関数+入れ子」のメリットは「他の関数を覚えなくていい」で、デメリットは終わりの))))))))の数が多すぎて間違えやすく、エラーに悩まされることです(笑)



2.case_when関数を使って分類する

3つ以上の分類で入れ子構造対策で慣れるとわかりやすいのがcase_when関数かもしれません。

case_when(条件A ~ 結果A, 条件B ~ 結果B, 条件C ~ 結果C)

今回は0歳以上 かつ 20歳未満といったように「かつ」が入る時は & を使います

年齢が0歳以上   かつ 年齢が20歳未満 を「0〜」
年齢が20歳以上 かつ 年齢が40歳未満 を「20〜」
年齢が40歳以上 かつ 年齢が60歳未満 を「60〜」・・・

上のイメージをコードにするとこうなります。

data_age <- data %>% 
  mutate(年代 = case_when(
    .$年齢 >= 0  & .$年齢 <20 ~ "0〜",
    .$年齢 >= 20 & .$年齢 <40 ~ "20〜",
    .$年齢 >= 40 & .$年齢 <60 ~ "40〜",
    .$年齢 >= 60 & .$年齢 <80 ~ "60〜",
    .$年齢 >= 80 & .$年齢 <100 ~ "80〜",
    .$年齢 >= 100 ~ "100〜")) %>% 
  mutate(年代 = factor(.$年代, levels = c("0〜", "20〜", "40〜", "60〜", "80〜", "100〜"), ordered = TRUE))

head(data_age)
スクリーンショット 2019-02-12 23.58.42


入れ子よりも( )の数が少ないのがありがたいところです。


数値でなくカテゴリーを条件で更に分類する場合

上記の場合は数値でしたが次のような場合はどうでしょう?

最初に行ったdata_20という変数名に以下のように分類していきます。
・成人という変数名で20歳未満を「未成年」、20歳以上を「成人」とする
・性別の男性を「male」、女性を「female」に変更
・疾患分類という変数名で、骨折・靭帯損傷を「運動器疾患」、脳梗塞・脳出血を「脳血管疾患」、心筋梗塞を「心疾患」とする

data_age <- data %>% 
  mutate(年代 = case_when(
    .$年齢 < 20 ~ "未成年",
    .$年齢 <= 20 ~ "成人")) %>% 
  mutate(性別 = case_when(
    .$性別 == "男性" ~ "male",
    .$性別 == "女性" ~ "female")) %>%
  mutate(疾患分類 = case_when(
    .$病名 %in% c("骨折", "靭帯損傷") ~ "運動器疾患",
    .$病名 %in% c("脳梗塞", "脳出血") ~ "脳血管疾患",
    .$病名 %in% "心筋梗塞" ~ "心疾患"))
head(data_age)
スクリーンショット 2019-02-13 0.58.43


ここで新たな ==%in% が出てきました。

スクリーンショット 2019-02-11 2.04.40
A == Bと A %in% Bは似ていますが、Bにあたる部分で==であれば1つしか入りません。
%in%は1つでも複数でも大丈夫です。

上記のコードで言うと .$性別 == "男性" ~ "male" の==は %in%に変えても大丈夫ですし、.$病名 %in% "心筋梗塞" ~ "心疾患" の %in% は == に変えても大丈夫です。

しかし.$病名 %in% c("骨折", "靭帯損傷") ~ "運動器疾患" に関しては ==に変えるとエラーが出ます。

上記図は今回の記事だけでなく条件式を作る時に必要になる考え方です。


case_when関数のメリットは入れ子構造にならないことで、デメリットとしては2つに分類するならifelseの方がコードが短くなります。



その他を作る(2019.09.17追記)

Exploratoryというソフトを開発しているKan Nishidaさんのツイートでcase_whenの解説を見つけました。ExploratoryはRをベースにしたデータ分析のソフトです。Rはプログラムを勉強してからでしか使えないところがありますが、ExploratoryはRコマンダーのようにプログラムを使わずにデータ分析ができます。

すべての人にデータサイエンスを

最新のデータサイエンスでは、さまざまな種類のデータを扱うことができ、機械学習や統計アルゴリズムを使用して、深い統計の知識がなくてもデータに隠れているパターンや傾向を見つけることができます。しかし、それはデータサイエンティストとプログラマーに限られています。

Exploratoryは誰もが、プログラミングを必要とせずに、モダンで最先端のデータサイエンスのアルゴリズムに簡単にアクセスできるようにします。私たちは誰もが最新のテクノロジーにアクセスし、データを通じて世界をより深く理解していくべきだと強く信じています。
(https://exploratory.io/より)





Exploratoryのサイトでもcase_when紹介がありました。

 


もし運動器疾患脳血管疾患以外をその他にしたければ以下のようになります。
TRUEはここでは残り全部といった意味になります。
data_age <- data %>% 
  mutate(年代 = case_when(
    .$年齢 < 20 ~ "未成年",
    .$年齢 <= 20 ~ "成人")) %>% 
  mutate(性別 = case_when(
    .$性別 == "男性" ~ "male",
    .$性別 == "女性" ~ "female")) %>%
  mutate(疾患分類 = case_when(
    .$病名 %in% c("骨折", "靭帯損傷") ~ "運動器疾患",
    .$病名 %in% c("脳梗塞", "脳出血") ~ "脳血管疾患",
    TRUE ~ "その他"

勉強になりました!


3.cut関数を使って分類する


数値を分類するに限って言うとcut関数も使えます。
cut関数はコードはスッキリしますが、少しクセがあり、おまじないが必要です。

cut(目的の列,
      breaks = (下限, カットする数値, 上限),
      right = FALSE,
      include.lowest = TRUE,
      labels =c("カテゴリー毎の名前"))

スクリーンショット 2019-02-13 1.49.04


「10代・20代」などでは0〜19のような分け方をするので、right = FALSEが必要になります。
加えてinclude.lowest = TRUEを加えないと端の値を読み込まず<NA>とデータなしとみなされてしまいます。


これらをふまえると以下のコードになります。

 data_age <- data %>% 
  mutate(年代 = cut(.$年齢, 
                  breaks = c(0, 20, 40, 60, 80, 100, 120),
                  right = FALSE, 
                  include.lowest = TRUE,
                  labels = c("0〜", "20〜", "40〜", "60〜", "80〜", "100〜")))
head(data_age)
スクリーンショット 2019-02-13 1.51.09


cut関数を使えば細かい条件が必要にならない分、right = や include.lowest = を忘れないように注意が必要です。


補足.mutate関数の組み合わせについて

mutate関数を複数回行うにはいくつかの方法があります。

スクリーンショット 2019-02-16 19.31.02
(追記)2019/02/16更新
すみません。mutate関数では②の方法は使えませんでした。
②の方法は次で紹介するfilter関数で使うことができます。


本来は,でつなぐと1つの行で済むのですが、初心者だと、()や,の数がうまくあわなかったりします。プログラミングの上級者からは意見があるかもしれませんが、プログラミング未経験・初心者はまず自分の覚えやすいものから始めていいと思っています。



まとめ

今回は条件によって分類する方法についてifelse関数、case_when関数、cut関数を紹介しました。

それぞれに特徴がありますので、必要に応じて使い分けてください。

また今回は演算子(<, ==, %in%など)の紹介も行いました。

演算子はデータを集計するときに必要になりますのでまた紹介していきたいと思います。

第2章ではデータハンドリングの基礎について紹介してきました。

【2-1】Rのfor関数、apply関数を使ってまとめて標準偏差などの統計量を求める方法

【2-2】Rのmutate関数を使って列の追加や修正を行う



今回は上記の復習として実際にデータを前処理し、集計をかけるところまで行います。

スクリーンショット 2019-02-08 10.33.40


1.データの準備

前回同様FIMのデータを使います。

FIM.xlsx 

今回は「入院時, 1ヶ月後, 退院時」全てのタブの使います。

スクリーンショット 2019-02-19 22.25.41


ダウンロードした後、プロジェクトの指定フォルダにファイルを移動させます。

プロジェクトの使い方がわからない場合は【1-6】Rstudioのプロジェクトについて解説しますをご参照ください。


まだ複数のデータを結合する方法を紹介していません。以下のコードを実行してください。

#必要なパッケージを読み込む
library(tidyverse)
library(readxl)

#それぞれのタブのデータフレームを作る
fim_in <- read_excel("FIM.xlsx", sheet = "入院時")
fim_1 <- read_excel("FIM.xlsx", sheet = "1ヶ月")
fim_out <- read_excel("FIM.xlsx", sheet = "退院時")

#bind_rows関数で3つのデータフレームを縦につなげる
fim <- bind_rows(fim_in, fim_1, fim_out)

bind_rows関数は複数のデータフレームを縦につなげる関数です。
fimという変数に入院時、1ヶ月、退院時全てのデータを縦につなげました。

これで準備完了です。


2.今回の目標

データにはFIM(18項目 + 運動合計 + 認知合計 + 全体の合計)のデータがあります。

更に今回同じ患者に3回反復測定を行っています。


<目標>
  • FIM各項目の平均点が入院時→1ヶ月→退院時でどう変化しているのかを表にする



課題①

まずデータの確認を行います。
  1. head関数、str関数を使ってデータを確認します。

回答は下にスクロールするとあります↓















head(fim)
スクリーンショット 2019-02-25 3.26.39

str(fim)
スクリーンショット 2019-02-25 3.27.15




課題②

データをlongデータに変えます。

まず列名と列番号を取得します。

t関数(またはdata.frame関数)とnames関数を組み合わせて列名と列番号を取得します。














t(names(fim))
スクリーンショット 2019-02-25 3.27.40

または
data.frame(names(fim))
スクリーンショット 2019-02-25 3.28.03



課題③

次はlongデータに変更します

  • パイプ演算子(%>%)を使います
  • 時期はbind_rows関数を使って縦に結合したので既にlongデータになっています
  • fimの「食事」〜「FIM全体」の列をlongデータに変えます
  • 今回はkeyの列名を「項目」、valueの項目を「点数」とします
  • 「項目」のfactorの要素を五十音順にせず、列で並んだ順で表示するようにします
















fim %>% 
  gather(5:25, key = 項目, value = 点数, factor_key = TRUE)
スクリーンショット 2019-02-25 3.55.14







課題④

次はsummarize関数を使い集計します。

  • fim_summarizeという変数名に作ります
  • 各項目が時期によってどう変化するのかが見たいのでした
  • 結果が項目→時期→点数と並ぶようにします。
  • 出た結果を見ると1箇所望まない結果になっている箇所があります。どこでしょう。














fim_summarize <- fim %>% 
  gather(5:25, key = 項目, value = 点数, factor_key = TRUE) %>% 
  group_by(項目,時期) %>% 
  summarize(平均 = mean(点数), 標準偏差 = sd(点数))
fim_summarize
スクリーンショット 2019-02-25 3.35.13



課題⑤
時期を見ると1ヶ月→退院時→入院時となっています。

時期のclassを確認します。















class(fim_summarize$時期)




課題④のコードに1行足して時期をfactor型に変え、入院時→1ヶ月→退院時の順に並ぶようにしてください。



















fim_summarize <- fim %>% 
  mutate(時期 = factor(時期, levels = c("入院時", "1ヶ月", "退院時"))) %>% 
  gather(5:25, key = 項目, value = 点数, factor_key = TRUE) %>% 
  group_by(項目,時期) %>% 
  summarize(平均 = mean(点数),
  標準偏差 = sd(点数))
fim_summarize
スクリーンショット 2019-02-25 3.29.54




課題⑥

今回の結果から「運動合計」「認知合計」「FIM合計」を取り除くにはどうすればいいでしょう。
fim_summarizeをつかって求めてください。
  • fim_summarize2 という変数名に作ります。
  • filter関数を使います
  • もし困ったら下の図を参考にしてください。
スクリーンショット 2019-02-11 2.04.40




















fim_summarize2 <- fim_summarize %>% 
  filter(!項目 %in% c("運動合計", "認知合計", "FIM合計"))
スクリーンショット 2019-02-25 3.31.41


まとめ

今回は第2章の一部の復習を行いました。

第2章では他にもifelse関数やcase_when関数を使って新たな変数名を作ったり、select関数で列を抽出しています。

今回のデータでも年代を作ったり、年代ごとにグループ分けすることも可能です。



fim %>% 
  mutate(年代 = cut(.$年齢, 
                  breaks = c(50, 60, 70, 80, 90),
                  right = FALSE, 
                  include.lowest = TRUE,
                  labels = c("50代", "60代", "70代", "80代")),
         年代 = as_factor(年代),
         時期 = factor(時期, levels = c("入院時", "1ヶ月", "退院時"))) %>% 
  gather(食事:FIM合計, key = 項目, value = 点数, factor_key = TRUE) %>% 
  group_by(時期, 項目, 年代) %>% 
  summarize(平均 = mean(点数), 
            標準偏差 = sd(点数),
            人数 = n())

スクリーンショット 2019-02-25 3.44.47


n関数はまだ紹介していませんでしたが、グループの数を表示することができます。
それ以外は過去の記事で紹介したものなので、まだしっくりこない方はサイトマップから過去の記事を探してみてください。



また色々な切り口があると思いますので、ぜひ色々試してみてください。


今までRでデータの変数名を変更したり条件でグループ化したりしました。







そして前回の記事では集計して平均や標準偏差などの要約を出したり、グラフを作るためにgather関数を使ってlongデータを作りました。

【2-5】Rでデータを集計するのに便利なtidyデータとgather関数



今回は棒グラフや折れ線グラフ作成に必要な平均や標準偏差などの統計量を求めます。



データは前回と同じデータを使います。

FIM.xlsx

スクリーンショット 2019-02-19 22.25.41


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

library(tidyverse)
library(readxl)
fim <- read_excel("FIM.xlsx", sheet = "入院時")
fim_long <- fim %>% 
  gather(食事:FIM合計, key = 項目, value = 点数, factor_key = TRUE) 
head(fim_long)


1,groop_by関数でグループ化したい項目を指定する

ますグループ化するためにはgroop_by関数を使います。

group_by(データ,列名)

%>%を使うとデータの部分は省略できます。

fim_long %>% 
  group_by(項目) 

スクリーンショット 2019-02-21 21.42.39



何も変わってないように見えますが、薄い文字のところに# Groups:項目 [21]とあります。


ちなみにグループの種類が複数でも可能です。

fim_long %>% 
  group_by(項目,性別) 

スクリーンショット 2019-02-21 21.45.43



2.列名の順番について

下の2つのコードはどういった違いがあるのでしょうか?

スクリーンショット 2019-02-22 0.15.27

表の並び順が違うだけでその後グラフを作るときには影響はないのですがイメージとしては上記のようになります。





3.統計量を出すsummarize関数


平均などの統計量などを出すにはsummarize関数を使います。

summarize(名前1 = 関数1, 名前2 = 関数2)

summarize関数の前にgroup_by関数を使っていると、グループごとの集計が出てきます。

fim_summarize <-  fim_long %>% 
  group_by(項目,性別) %>% 
  summarize(平均 = mean(点数), 
            標準偏差 = sd(点数), 
            最小値 = min(点数), 
            最大値 = max(点数))
fim_summarize

スクリーンショット 2019-02-22 0.49.09




もし保存したい場合はwrite.csv関数を使います。

wite.csv(保存する変数名, "ファイル名.csv")
ファイル名には" "と.csvを入れます

write.csv(fim_summarize,"FIM集計")
スクリーンショット 2019-02-22 0.57.54
右のfilesビューにFIM集計.csvができました。

csvファイルをExcelで読み込む時は【1-11】Rで医療統計で必要なtable1を作るtableoneパッケージについて紹介しますをご参照ください。

ポイントとしてはファイルの出力先を$A$1ではなく$B$1にしてください。
$A$1だとなぜかエラーが出ます。
スクリーンショット 2019-02-22 1.04.57

スクリーンショット 2019-02-22 1.07.22

Excelにするとわかるのですが、小数点10桁まで表示されます。
もし小数点第一位までの表示にしたい時はround関数を使います。

round(数値, x)
たとえばxが1だと小数点第二位を四捨五入して小数点第一位まで表示します。

fim_summarize <- fim_long %>%
group_by(項目,性別) %>%
summarize(平均 = round(mean(点数), 1),
   標準偏差 = round(sd(点数), 1),
   最小値 = min(点数),
   最大値 = max(点数))
fim_summarize
スクリーンショット 2019-02-22 1.18.28



4.summarize関数を使うときの注意点。

group_by関数の順番に影響する

2でも紹介しましたが、groop_by関数で指定したグループが複数の場合、summarize関数で表示される順番はgroop_by関数の影響を受けます。

スクリーンショット 2019-02-22 0.15.27

あとでExcelで並べ替えるのはただ面倒です。もしcsvで保存をする時は何を示したいかをあらかじめイメージしておくことが重要になります。



2.欠損地があるとNAとなる

平均を求めるmean関数などはどこか1つでも欠損値(空欄)があると結果はNAとなります。

氏名 <- c("A", "B", "C", "D")
年齢 <- c(55, 63, 67, 71)
test_1回目 <- c(1,2,3,NA)
test_2回目 <- c(5:8) test_3回目 <- c(9:12)
data <- data_frame(氏名, 年齢, test_1回目, test_2回目, test_3回目)
data data %>%
gather(3:5, key = 回数, value = 点数) %>%
group_by(回数) %>%
summarize(平均 = mean(点数))
スクリーンショット 2019-02-22 1.53.38


もし結果にNAが出た時はまずはそもそものデータの入れ損ねがないか確認をし対応します。

それでも欠損値がある時はNAを省いて計算する欠損値を統計の技術を使って代入するといった方法があります。

もし欠損値を省いて平均を出す場合はna.rm = TRUEを付け加えます。

data %>%
gather(3:5, key = 回数, value = 点数) %>%
group_by(回数) %>%
summarize(平均 = mean(点数, na.rm = TRUE))
スクリーンショット 2019-02-22 2.06.51


ただ実は欠損値を省いた方がいいのかどうかという問題があります。ただ業務で傾向を確認したいなどであれば欠損値を省いてもいいと思いますが、きちんと出さないと行けない場面ではこれはこれできちんと勉強する必要があります。欠測データに関しては医療統計の本では紹介されていない事が多く専門書が必要かもしれません。


欠測データ処理: Rによる単一代入法と多重代入法 (統計学One Point)



summarize関数に入れられる関数は単一の値が出るものに限る。

summarize関数で使える統計量はmean関数,sd関数,median関数など単一の値になります。

ただ、関数の中には最小値と最大値を一度に出してくれるrange関数など複数の値を出すものがあります。

test_2回目 <- c(5:8)
range(test_2回目)
スクリーンショット 2019-02-22 2.27.00


range関数をsummarize関数に入れようとするとエラーが出ます。

data %>% gather(3:5, key = 回数, value = 点数) %>% group_by(回数) %>% summarize(平均 = mean(点数, na.rm = TRUE), 範囲 = range(点数))

スクリーンショット 2019-02-22 2.36.38

「1つの値しか入れられないのに2つ入ってるよ!」と怒られています。

スクリーンショット 2019-02-22 2.41.49

特にあるのが四分位範囲です。
四分位範囲はquantile関数を使って以下のように一度に値を出すことができます。
quantile(fim$年齢, c(0.1, 0.25, 0.5, 0.75, 0.9))
スクリーンショット 2019-02-22 2.45.59


しかしsummarize関数ではまとめて入れられないので1つずつ入れる必要があります。
fim_summarize <- fim_long %>% group_by(項目,性別) %>% summarize(平均 = mean(点数), 標準偏差 = sd(点数), 最小値 = min(点数), percent_25 = quantile(点数, 0.25), 中央値 = median(点数), percent_75 = quantile(点数, 0.75), 最大値 = max(点数)) fim_summarize
スクリーンショット 2019-02-22 3.01.32



そもそもの列のfactorの順番があってるのか?

gather関数のkey列に関してはfactor_key=TRUEで五十音順ではなく、元の順番に戻すことができます。

しかし他の列でfactorの順番が合っていない可能性もあります。

もしほかの列でfactorの順番を合わせるにはmutate関数とfactor関数を組み合わせて使うことができます。

今回は女性と男性を入れ替えてみます。
fim_summarize <- fim_long %>%
mutate(性別 = factor(性別, levels = c("男性","女性"))) %>%
group_by(項目,性別) %>%
summarize(平均 = mean(点数),
     標準偏差 = sd(点数),
      最小値 = min(点数),
   percent_25 = quantile(点数, 0.25),
   中央値 = median(点数),
   percent_75 = quantile(点数, 0.75),
      最大値 = max(点数))
fim_summarize
スクリーンショット 2019-02-22 3.01.32


factor関数の使い方は【1-10】Rでよく使われる型について説明しますをご参照ください


まとめ

今回はgroup_by関数とsummarize関数と、実際に集計を出し保存するところまで紹介しました。

これでひとまず第2章で行う予定だった「データを集計しやすいように形を整え集計する」が終わりです。

スクリーンショット 2019-02-08 10.33.40


ただ今回は「できるだけExcelの時点でデータは編集しやすい形式で保存している」ことを前提に話を進めています。

もっと形が整っていないデータの扱いやここでは説明できなかった項目も多くあります。

もし「もっと知りたい」「ここの情報では足りない!」ということであれば下記のサイトなどもご参照ください。


データハンドリング入門
https://kazutan.github.io/kazutanR/hands_on_170730/index.html




↑このページのトップヘ