タグ:統計

*最終更新日: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以外の表にも対応しています。
そこに関しては今後の記事で紹介します。



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

第4章ではt検定から始まり色々な種類の統計やp値や信頼性の考え方などを紹介してきました。

そして内容は後半の多変量解析に入り、前回は多変量解析を行う上で大切な考え方を紹介しました。

【4-17】「先にデータを取って後で統計」がなぜ危険なのか?どうすればいいのか?を多変量解析、相関・因果関係、交絡を使って解説してみる


そしてもう1つ重要な考え方があります。それは洞察(因果関係)を見たいのか?それとも予測がしたいのか?です。

これは統計学と機械学習の違いにも通じる話になります。

1.使うデータ

ここでも前回の記事と同じデータを使います。

架空のデータです。

library(tidyverse)
set.seed(25) v1 <- rnorm(100,87,13) v2 <- rnorm(100,63,10.5) v <- c(v1,v2) a1 <- as.integer(rnorm(100,35,5)) a2 <- as.integer(rnorm(100,55,5)) age <- c(a1,a2) sex <- c(sample(c(0,1),size = 100, prob = c(0.7,0.3),replace = TRUE),sample(c(1,0),size = 100, prob = c(0.7,0.3),replace = TRUE)) t1 <- as.integer(rnorm(100,120,20)) t2 <- as.integer(rnorm(100,70,20)) time <- c(t1,t2) df <- data.frame(v,time,age,sex) df

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

v:歩行速度、1分間に歩いた距離(m)とします。

time:片脚立位の時間(s)

age:年齢

sex:性別(0:男性,1:女性)

これをみて片脚立位の時間が長くなれば歩行速度が上がると考えていいかどうかとします。




2.多変量解析はどういう仕組みか?

スクリーンショット 2020-07-17 14.34.46

スクリーンショット 2020-07-17 14.27.48



まず多変量解析(ここでは重回帰分析)は左辺(y)にアウトカム、右辺に要因(x)をおいた式となります。yとxには色々呼び方がありますが同じ意味です。

そして数学的なテクニックを使い(最小二乗法や最尤法)を使い係数を求めます。

コードと解釈の仕方は次の記事で紹介します。

fit <- lm(v ~ time + age + sex)
summary(fit)

スクリーンショット 2020-07-17 14.47.59

色々見るところはあるのですが、とりあえず黄色に注目します。

Estimateというのが係数(β)になります。(Intercept)は切片です。

1分で歩いた距離 = 0.03×片脚立位時間(秒)+ -0.84×年齢(歳) + 0.42(女性であれば) + 108.33

ロジスティック回帰もそうですが、やっていることは係数を求めているということです。


3.重回帰分析とロジスティック回帰の使い分け方

重回帰分析とロジスティック回帰の使い分けはアウトカムを見ます。

重回帰分析:アウトカムが連続した数値になっている
ロジスティック回帰:アウトカムが「自立・非自立」といった2択になっている


ちなみにアウトカムが3択以上のある場合は自分の知りたいことにあわせて2択にしてしまったほうがいいです。
3択以上の解析は複雑になります。


4.xを見たいのか?yを見たいのか?

じつはこの数式には2つの見方があります。それはxを見たいのか?yを見たいのか?です。


洞察(因果関係)に興味がある=xに興味がある

まずxを見るときの視点を説明します。

上の式を再掲します
1分で歩いた距離(m) = 0.03×片脚立位時間(秒)+ -0.84×年齢(歳) + 0.42(女性であれば) + 108.33

この式は右辺に片脚立位時間、年齢、性別という3つの変数があります。

ここで年齢と女性でないという条件で片脚立位の時間が1秒増えるとどうなるでしょうか?

答えは歩行距離が0.03m増えるということになります。

逆に片脚立位・女性でない条件で年齢が1つ増えれば歩行距離は-0.84m減ります

こうやって係数を見ることで他の条件は同じだとして知りたい要因が1つ値が変わるとアウトカムはどう変化するのかを見ることができます。そういう意味ではyの値自体には大して意味がありません。


ここで前回の記事の交絡因子が関係します。

【4-17】「先にデータを取って後で統計」がなぜ危険なのか?どうすればいいのか?を多変量解析、相関・因果関係、交絡を使って解説してみる

因果関係(効果があるかどうか)を考える時はどれだけ重要な交絡因子が取り除けているかが説得力に繋がります。そのためにはどの因子が大事な交絡因子なのかを考える必要がありますが、それはその分野の知識によって決まります。
スクリーンショット 2020-07-17 15.10.39


よって大事なことは右辺には興味がある要因と交絡因子をできるだけ入れることが大切です。

そして知りたい要因は介入が可能であるものが基本となります。
例えば知りたい要因を年齢にしてしまうと、もし有意差が出たとしても年齢を変えることはできません。

なぜここを強調するかというと統計の手法には変数選択という自動的に変数を削除することができるオプションがあるからです。

医療統計でいちばん有名なのはステップワイズ法です。
(下図はEZRの重回帰分析の設定画面)
スクリーンショット 2020-07-17 15.20.46

交絡因子の影響を除きたいのに交絡因子自体を計算から除いてしまうとそれは本当にこの研究でいいたかったことなのか?となってしまいます。

そのためxに興味があるときは変数選択は基本行いません(強制投入といいます)。

*ただ本によっては意見が違ったりします。理学療法の研究ではステップワイズが使われていることも多く、もしかしたら査読で「ステップワイズを」と言われる可能性はあるかもしれません。その雑誌に合わせることも必要なのかもしれませんし、生物統計家のご意見も伺ってみたいところです。

以下の2サイトではステップワイズはよくないという記事になります。







予測に興味がある=yに興味がある。

次にy(左辺)を見るときの視点を紹介します。

上の式を再掲します
1分で歩いた距離(m) = 0.03×片脚立位時間(秒)+ -0.84×年齢(歳) + 0.42(女性であれば) + 108.33

これは片脚立位時間、年齢、性別がわかれば左辺の1分間で歩いた距離を予測することができるということになります。

スクリーンショット 2020-07-17 17.02.35

つまり右辺はどれだけいい予測になる因子を集められるかが大切になります。
そのため交絡因子がどうかや介入が可能かどうかは重要な問題ではありません。
ただ、勘違いしやすいのですが因子が多ければいいモデルになるわけではありません
先程と違い精度に貢献しない因子を外すこともします。

そして予測に興味がある場合に大切なのはモデルを作ることより作ったモデルが新しいデータでも通用するかを試すことで信頼性が増すということです。

これは今回のデータのみに当てはまりがよく、未知のデータでは通用しないという過学習という問題が起こっている可能性があるからです。モデルを作っただけでは過学習が起こっているかどうか判断できません。今回のモデルで使わなかったデータで予測がどれくらい的中するか再検証する必要があります。誤差だったり精度・感度・特異度はこういった再検証に用いられます。



加えて医療統計ではあまりみられませんが、取ったデータから新たな変数を作ることも可能です。

例えばアウトカムを歩行距離として、下肢の左右各関節で疼痛(VAS)を測っていたとします。ただ痛いところは人それぞれですし、変数が多くなりすぎるので変数をまとめたいとします。
  • 各関節の最大のVASを採用する
  • 各関節の平均点を採用する
  • VAS5を超えた関節の数を採用する
  • 他の点数とかけ合わせた新たな変数を作る
他にも色々あるかもしれません。そうやって精度の高い予測を作ることになります。
こういった作業は特徴量エンジニアリングといい、機械学習のテクニックになります。

ただし解釈の難しい因子を作ってしまうと解釈(説明可能性)が難しくなります
上の例だと痛い関節の数はイメージが付きますが、疼痛×筋力みたいな数値はその数値にどういった意味があるのか?と聞かれても「予数値で臨床的な意味は・・・だけど予測の精度が上がる数値」となってしまします。

精度を追求するか説明可能性を追求するかはトレードオフな関係ではありますが、医療現場では説明可能性がある因子を使った方が作った方が使ってくれやすいという印象があります。理由を説明しにくいモデルは「ホントにー?」と胡散臭く思われる可能性があります。

まずは説明しやすいモデル、その後に需要があれば難しいモデルに挑戦するような考えがいいかと思います。


5.因果関係でも予測でも入れないほうがいい変数がある

ただ因果関係・予測どちらにも言えることですが入れないほうがいい因子があります。それは多重共線性が関係します。

数ある因子同士の中で相関係数が高すぎるものが混じっていると多重共線性が発生し結果が不安定にあるという特徴があります。正確な表現ではありませんが、計算の過程で◯÷0みたいな変なことが起こってしまうイメージです。

そのため独立変数同士で相関行列を求め、相関係数が高いものがあるかを調べたりVIFというものを確認したりします。

もし多重共線性が見つかれば相関が高い因子のうちのどちらかを省くか、合成した新しい変数を採用するかになります。

相関行列の求め方は【4-16】Rで相関行列を作成するGGally, cor, clrrplot, corrr関数を紹介しますで紹介しています。



6.まとめ

重回帰分析やロジスティック回帰という回帰分析は計算方法やコードの書き方は同じでも医療統計として使う場面と機械学習として使う場面があります。知りたい要因の因果関係に興味があるのか?アウトカムの予測に興味があるのかで集めるデータや解析や検証の行い方に違いが出てきます。


知りたい要因の因果関係に興味がある→医療統計的な考え方
スクリーンショット 2020-07-17 15.10.39
→データを取り始める前に重要な交絡因子が無いかを十分に検討する
→選択バイアス、情報バイアスもできるだけ入らないよう注意する
→多重共線性がないか確認し、あれば対処する
→分析の結果から考察する


アウトカムの予測に興味がある→機械学習的な考え方
スクリーンショット 2020-07-17 17.02.35
→データを取り始める前に予測に重要な要因は何かを検討する
→多重共線性を調べる
→データを予測作成用のデータと検証用のデータに分ける
→予測作成用のデータでモデルを作る。できたモデルの精度を確認する
→検証用のデータで精度を再検証する
(モデルを作成しただけだと実は本当の精度はわからない)


医療統計ではほとんど前者になることが多い、もしくはどっちつかずになっていることが多い印象ですが、どちらに興味があるかが意識できると多変量解析の勘所が上がるのではないかなと思います。


前回・今回と長くなってしまいましたが、次は重回帰分析の実際について紹介します!


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


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


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


 


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


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


今回はFriedman(フリードマン)検定を紹介します



まずFriedman(フリードマン)検定についてはハルさんのサイトをご参照ください。



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

 

1.準備

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

 

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


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

次にRのコードを書くためのスクリプトファイルを作ります。

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

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


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

デモデータ(Friedman検定)

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


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


実際のコードは以下になります。
前回のコードのURL(" "の中)とdestfileのdata/以降を変更するだけでOKです。
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-friedman-rank-sum-test.xlsx"
destfile = "data/demo-friedman-rank-sum-test.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からC22までデータが入っています(B2:D32と表記)
スクリーンショット 2019-12-08 20.31.58



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



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

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

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


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

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


これでデータの取り込みは完成です。
Friedman検定は繰り返しのデータ(被験者内要因)なのでデータを横につなげています。


5.Friedman検定を行う(wideデータの場合)

ハルさんのサイトでは早速Friedman検定を行っていますのでこちらも早速行います。
今回はwideデータなのでまずwideデータの場合を紹介します。

Friedman検定はfriedman.test関数を使います。

friedman.test(cbind(繰り返しデータのある列))
スクリーンショット 2019-12-09 0.52.59

friedman.test(cbind(grip_friedman$W0,grip_friedman$W1,grip_friedman$W2))
スクリーンショット 2019-12-09 0.54.04

EZRと同じ結果が出ました。

6.Friedman検定を行う(longデータの場合)
次にlongデータでFriedman検定を行う方法も紹介します。

friedman.test(データ ~ グループ | 個人を識別するid, data = ○○)

スクリーンショット 2019-12-09 0.58.27

Friedman検定は対応のあるデータになります。

今回は30人がW0,W1,W2の3回データを取っているので3*30=90データあります。
そのため作るlongデータは上記のイメージになります。

今回はid, time, gripとしてますが自分がわかれば好きな名前を付けてOKです。
id:30人の内誰のデータかわかるためのid(氏名でも可)
time:W0, W1, W2の3つ
grip:今回測定した握力のデータ
wideデータをlongデータに変えるにはgather関数、もしくはpivot_longer関数を使います。
gather関数に関しては【2-5】Rでデータを集計するのに便利なlongデータとgather関数で紹介しています。



データの扱いはtidyverseパッケージに含まれていますので。まずlibrary関数でtidyverseパッケージを読み込みます。今回はgrip_friedman_longという名前にします。

library(tidyverse)

#gather関数の場合 grip_friedman_long <- grip_friedman %>% rowid_to_column("id") %>% gather(key = time, value = grip, W0:W2, factor_key = TRUE)

#pivot_longer関数の場合 grip_friedman_long <- grip_friedman %>% rowid_to_column("id") %>% pivot_longer(cols = c(W0,W1,W2), names_to = "time", values_to = "grip", names_ptypes = list(time = factor()))

ポイントはもともとのExcelデータにidにあたる列がないことです。
そこでrowid_to_column関数を使います。""の中身が列名になります。
スクリーンショット 2019-12-09 1.14.03

これで準備が整いました。

longデータでFriedan検定を行う場合は改めてこうなります。

friedman.test(データ ~ グループ | 個人を識別するid, data = ○○)
friedman.test(grip ~ time | id, data = grip_friedman_long)
スクリーンショット 2019-12-09 0.54.04
もちろん結果は同じです。

7.多重比較を行う
Friedman検定で多重比較を行いたいのですがwideデータだとEZR専用のパッケージでしかできません。
そのためRで行うにはlongデータで行う必要があります。

longデータで多重比較を行うにはpairwise.wilcox.test関数を使います。簡単に言うと全ての組み合わせでMann-Whitney U検定を行いその結果を元にp値の調整を行います
どの方法を使うかはp.adjust.method = で指定します。
またタイのエラーを回避するためにexact = FALSEを入れています。
data = 

pairwise.wilcox.test(目的変数factor(グループ),  p.adjust.method = "◯◯", paired = TRUE, exact = FALSE)

#Bonferroniの方法
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "bonferroni", paired = TRUE, exact = FALSE)
#Holmの方法
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "holm", paired = TRUE, exact = FALSE)
スクリーンショット 2019-12-09 1.39.44
EZRとほぼ同じ結果になりました。

8.どうしてもwideデータのままなら手計算が必要

wideデータのまま多重比較を行う方法がEZRのパッケージしかありません。

そのためどうしてもwideデータのまま行うにはwilcox.test関数、もしくはcoinパッケージwilcox_test関数を使い、1つずつ組み合わせを行い、BonferroniHolmHochbergそれぞれの方法で自分でp値を直接計算することになります。

この方法は統計ERさんの記事に詳しく紹介されています。
更にBonferroniHolmHochbergって何をしてるの?という疑問も解決するかもしれません。
仕組みがわかればp値に2とか3とかかけるだけなので意外と単純です。
 


8.まとめ

今回はFriedman検定を紹介しました。

もしRで行う場合はwideデータで分析するかlongで分析するかをあらかじめイメージしていた方がその後の解析もスムーズかもしれません。

次は相関係数について紹介します!

9.今回使ったコード
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-friedman-rank-sum-test.xlsx"
destfile = "data/demo-friedman-rank-sum-test.xlsx"

download.file(url, destfile)

library(readxl)
grip_friedman <- read_excel("data/demo-friedman-rank-sum-test.xlsx")
View(grip_friedman)

friedman.test(cbind(grip_friedman$W0,grip_friedman$W1,grip_friedman$W2))

#pivot_longer関数の場合
grip_friedman_long <-
  grip_friedman %>% 
  rowid_to_column("id") %>% 
  pivot_longer(cols = c(W0,W1,W2), names_to = "time", values_to = "grip", names_ptypes = list(time = factor())) 

#gather関数の場合
grip_friedman_long <-
  grip_friedman %>% 
  rowid_to_column("id") %>% 
  gather(key = time, value = grip, W0:W2, factor_key = TRUE)

friedman.test(grip ~ time | id, data = grip_friedman_long)

pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "bonferroni", paired = TRUE, exact = FALSE)
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "holm", paired = TRUE, exact = FALSE)
 

記事一覧はこちらになります。


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


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


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


 


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


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


今回はWilcoxon符号付順位和検定を紹介します



まずWilcoxon符号付順位和検定についてはハルさんのサイトをご参照ください。



また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-04 0.47.12



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

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

デモデータ(wilcoxon符号付順位和検定)


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


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


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

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

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




4.データの読み込み

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

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

スクリーンショット 2019-11-04 7.44.03



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


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

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

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

スクリーンショット 2019-11-04 7.41.10

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

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

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


5.データを要約する

次にデータを要約します。

データの要約は【4-1】Rでt検定を行う方法で紹介しました。
group_by関数とsummarize関数を使って要約しましたが、今回はsummary関数を使います。
summary関数は平均、中央値、最大・最小値、四分位範囲をまとめて出してくれますが、標準偏差はだしてくれません。ただ今回はノンパラメトリックなので標準偏差はいらないだろうという理由です。

そして今回はwideデータになっています。
スクリーンショット 2019-11-04 8.05.00
なのでコードはこうなります。
summary(borg$pre)
summary(borg$post)
スクリーンショット 2019-11-04 19.56.42


6.Wilcoxon符号付順位和検定を行う

次にWilcoxon符号付順位和検定を行います。

【4-3】Rで対応のあるt検定を行う方法で紹介しましたが、対応のあるとなしはt.test関数paired = TRUEをつけるかどうかの違いでした。

実はWilcoxon符号付順位和検定も同じです。

wilcox.testにpaired = FALSEをつける(もしくは何も付けない)とMann-Whitney U 検定
wilcox.testにpaired = TRUEをつけるとWilcoxon符号付順位和検定


ということで、wilcoxテストを行ってみます。
t.test関数もそうでしたが、longデータとwideデータで書き方が違います。

wideデータの場合 → , を使う
wilcox.test(1列目, 2列目, paired = TRUE)

longデータの場合 → を使う
wilcox.test(数値 ~ グループ, paired = TRUE)


今回はwideデータなのでこうなります。
wilcox.test(borg$pre, borg$post, paired = TRUE)
スクリーンショット 2019-11-04 20.38.57

EZRと同じ結果になりましたが、Mann-Whitney U 検定のときに悩ませたたアレが出てきました。
今度はもう1行増えてます。
タイがあるため、正確な p 値を計算することができません 
ゼロ値のため、正確な p 値を計算することができません
ちなみにEZRでも警告が出ています。
スクリーンショット 2019-11-04 20.59.27

cannot compute exact p-value with ties
cannot compute exact p-value with zeroes


7.タイのあるデータの対処法

EZRでも使われているwilcox.test関数はタイ(同順位)があると正確なp値を計算できず、近似値を計算する設定になっていました。





今回も同じ問題が出ています。
ブログではある程度のn数があればEZR(wilcox.test)でもいいのではという話がありました。
ただ警告が気持ち悪い!正確なp値も知りたいというための方法も紹介します。

タイに対しては奥村先生の記事が参考になります。



①coinパッケージのwilcoxsign_test関数

Mann-Whitney U 検定ではタイがあっても正確なp値を計算するcoinパッケージwilcox_test関数がありました。coinパッケージを使ってWilcoxon符号付順位和検定を行う場合はwilcoxsign_test関数を使います。

まだcoinパッケージを1度も使ったことがなければインストールします。

coinパッケージのインストール
install.packages("coin")

wilcoxsign_test関数の書き方はちょっとクセがあります。。。

スクリーンショット 2019-11-04 23.11.32
#パッケージの読み込み
library(coin)

wilcoxsign_test(borg$pre ~ borg$post, distribution = "exact", zero.method="Wilcoxon")
スクリーンショット 2019-11-04 23.12.57


exactRankTestsパッケージのwilcox.exact関数

もう1つexactRankTestsパッケージがあります。
このパッケージは開発が終わっており、インストールするとcoinパッケージ使ってねと警告が出ます。それでもcoinパッケージのwilcoxsign_test関数と同じ結果になります。
スクリーンショット 2019-11-05 0.06.33

まずexactRankTestsパッケージをインストールします。
#exactRankTestsパッケージのインストール
install.packages("exactRankTests")
wilcox.exact関数はwilcox.testと似たような書き方ができるのでわかりやすいのが特徴です。
スクリーンショット 2019-11-04 23.53.16
library(exactRankTests)
wilcox.exact(borg$pre, borg$post, paired = TRUE, exact=TRUE)
スクリーンショット 2019-11-04 23.55.13

どちらも同じになりました。


8.まとめ
今回はWilcoxon符号付順位和検定を紹介しました。
Mann-Whitney U 検定との共通点や相違点を比較するとイメージが深まると思います。

4章を順に見ていくと重複する箇所も出てきますので検索で来られた方はサイトマップを見ていただければ別の発見があるかもしれません。

 


9.今回使ったコード

今回使ったコードをまとめて置いておきます。

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

download.file(url, destfile)


library(readxl)
borg <- read_excel("data/demo-wilcoxon-rank-sum-test.xlsx", 
                   range = "B2:C22")
View(borg)

summary(borg$pre)
summary(borg$post)

summary_pre <- summary(borg$pre)
summary_post <- summary(borg$post)


#Wilcoxon符号付順位和検定
wilcox.test(borg$pre, borg$post, paired = TRUE)


#coinパッケージのインストール
install.packages("coin")

library(coin)
wilcoxsign_test(borg$pre ~ borg$post, distribution = "exact", zero.method="Wilcoxon")

#exactRankTestsパッケージのインストール install.packages("exactRankTests") library(exactRankTests) wilcox.exact(borg$pre, borg$post, exact=TRUE, paired = TRUE)

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


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


ハルさん、ありがとうございます!


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


 


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


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


今回はMann-Whitney U 検定を紹介します



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




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



1.準備

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



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


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

次にRのコードを書くためのスクリプトファイルを作ります。

スクリーンショット 2019-10-29 19.25.29


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

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

スクリーンショット 2019-11-01 1.21.02


デモデータ(Mann-Whitney U検定)

前回Rを使ってダウンロードしましたが、まだ慣れないと思いますので再度説明します。


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


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

download.file(url, destfile)
スクリーンショット 2019-10-29 19.43.04


dataフォルダにダウンロードできました!
スクリーンショット 2019-10-29 19.36.48


4.データの読み込み

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

View Fileでデータを確認します。
スクリーンショット 2019-10-29 19.46.19

データが入っているセルを確認します。
B2からC62までデータが入っています(B2:C62と表記)

スクリーンショット 2019-10-29 19.47.48

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


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

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




Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。
library(readxl)
mmt <- read_excel("data/demo-mann-whitney-u-test.xlsx", 
                  range = "C2:D62")
View(mmt)
スクリーンショット 2019-10-29 21.03.46

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

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


5.Mann-Whitney U 検定を行う

早速Mann-Whitney U 検定を行ってみます。

Mann-Whitney U 検定(マン・ホイットニーのU検定)は別な言い方もあります。

Wilcoxon-Mann-Whitney検定
Wilcoxon の順位和検定


EZRで使われているMann-Whitney U 検定はwilcox.test関数を使っています。

wilcox.test(MMT ~ category, data = mmt)

スクリーンショット 2019-11-01 1.02.29
スクリーンショット 2019-11-01 0.40.49


EZRの結果はこうです。
スクリーンショット 2019-11-01 1.07.20


EZRと同じ結果がでました。警告が気になりますがひとまず置きます。

もしp値だけがほしいときは変数名を付けて変数名$p.valueとします。
mmt_wilcox <- wilcox.test(MMT ~ category, data = mmt)
mmt_wilcox$p.value


6.グラフの作成
EZRでは同時に箱ひげ図を作成しますがRでは自作する必要があります。

箱ひげ図の作成に関しては第3章で紹介しています。



ggplot(data = mmt) +
  geom_boxplot(aes(x = category, y = MMT))
スクリーンショット 2019-11-01 1.16.42

もしEZRと同じグラフにしたい場合はboxplot関数を使うという選択肢もあります。
boxplot(MMT ~ category, data = mmt)
スクリーンショット 2019-11-01 1.19.01



7.一旦まとめ

RでMann-Whitney U 検定を行う方法を紹介しました。
wlcox.testを使えばEZRと同じ結果になります。

ただ警告が出てしまい心配な方はいると思います。

以下からは少し応用編となりますがこの警告について進めていきます。




8.タイがあるため、正確な p 値を計算することができません

スクリーンショット 2019-11-01 0.40.49

Rでノンパラメトリックの検定を行っているとこういった画面が出ることがあります。

cannot compute exact p-value with ties
タイがあるため、正確な p 値を計算することができません
正確なp値が計算できないと言われるととても心配になります。

実は先程のEZRでも警告は出ていました。
一番下の行です。
スクリーンショット 2019-11-01 1.35.45

cannot compute exact p-value with ties


8.タイ(ties)って何?

タイというのは同順(同じ数値)のことです。

タイがないデータ:1,2,4,5,7,10,11,12
タイがあるデータ:1,2,2,4,6,6,6,7,9,11


詳細は成書に譲りますがMann-Whitney U 検定は平均ではなくデータを少ない順に並べてデータに順位を付け計算するため、同順位があるとそれをどう計算するかといった問題や様々な対処法が出てきます。


9.タイに対してどう考えるか?

正確なp値を求めるには膨大な計算が必要になるらしく、そのため正規近似というものを使って計算するようです。EZR(wilcox.test)はタイがなければ正確なp値を求めるオプションがありますが、タイがあるとオプションを付けても正確なp値を求めることができないようです。近似値ではなく正確なp値を求めるにはcoinパッケージwilcox_test関数を使う必要があります。coinパッケージに関しては奥村先生の記事や裏 RjpWikiの記事に詳しく説明があります。






10.coinパッケージを使ってみる

wilcox_test関数を使うためにはcoinパッケージをインストールする必要があります。
パッケージのインストールは今までinstall.packages関数を紹介していましたが、今回はRStudioのPackagesタブからインストールしてみます。
PackagesタブのInstallをクリックしてパッケージ名を入れます。

スクリーンショット 2019-10-29 21.37.10

これでインストールは完了です。
これでlibrary関数を使ってcoinパッケージを呼び出せます。
wilcox.testとの違いは主に3点あります。
①グループの箇所は必ずfactor型に変換(as.factor関数またはfactor関数を使う)
(わからないときは【4-1】Rでt検定を行う方法を参照)
②正確なp値を求めるためにはdistribution = "exact"を追加
③実はwilcox.testはlongデータでもwideデータデータでも使えるが、wilcox_testはlongデータしか扱えない
(わからないときは【1-10】Rでよく使われる型について説明しますを参照)
スクリーンショット 2019-11-02 4.52.27

library(coin)
mmt_coin <- wilcox_test(mmt$MMT ~ as.factor(mmt$category), distribution="exact") mmt_coin
スクリーンショット 2019-11-02 4.35.30


もしp値だけを出したいときはpvalue関数、95%信頼区間を出したいときはpvalue_interval関数を使います。
pvalue(mmt_coin)
pvalue_interval(mmt_coin)
スクリーンショット 2019-11-02 4.35.45


これで正確なp値を求めることができました。


11.じゃあEZRはダメなのか?

ただ必ず正確なp値を求めないといけないかというとそうではないようです。

ちゅんたろさんからアドバイスを頂きました。







ある程度n数があればEZRでも利用可能と考えることもできます。
このあたりは統一見解があるわけではなさそうです。

そしてそもそもp値がより小さければより差が大きいわけではありませんので、ここのp値をどうするかよりも研究デザインが適切かどうかの方が大切になるかもと言うのが個人的なイメージです。


12.今度こそまとめ

今回はMann-Whitney U 検定を紹介しました。
実はwilcox.test(EZR)でタイが無いときの正確なp値の求め方などwilcox.testのオプションなど話は割愛させてもらっています。

qiitaの記事ではwilcox.testとwilcox_testのオプションやタイの有無で結果がどう変わるかについての解説をしていますので、興味があればぜひこちらも併せてご覧ください。



13.今回使ったコード

今回使ったコードをまとめて置いておきます。
coinパッケージのインストールに関してもコードにしています。

#Mann-Whitney U 検定を行うためにcoinパッケージのインストール
#1度でもしていればしなくてOK
install.packages("coin") #もしtidyverseパッケージをインストールしていなければインストール。 #1度でもしていればしなくてOK install.packages("tidyverse") #パッケージの読み込み library(coin) library(readxl) library(tidyverse) #データのダウンロード url <- "https://haru-reha.com/wp-content/uploads/2018/04/demo-mann-whitney-u-test.xlsx" destfile = "data/demo-mann-whitney-u-test.xlsx" download.file(url, destfile) #データの読み込み mmt <- read_excel("data/demo-mann-whitney-u-test.xlsx", range = "B2:C62") View(mmt) #wilcox.test wilcox.test(MMT ~ category, data = mmt)
boxplot(MMT ~ category, data = mmt) #wilcox.testでp値を求めたい時(EZRと同じ) mmt_wilcox <- wilcox.test(MMT ~ category, data = mmt) mmt_wilcox$p.value #coinパッケージのwilcox_testを使う mmt_coin <- wilcox_test(mmt$MMT ~ as.factor(mmt$category), distribution="exact") mmt_coin pvalue(mmt_coin) pvalue_interval(mmt_coin)

↑このページのトップヘ