医師が教えるR言語での医療データ分析入門-2022年改訂 - 215 Rmarkdownでレポートを作成してみる―課題1

課題1:

kadai/dataフォルダに含まれる、data.csvファイルにはある会社の健康診断の（架空の）データが含まれています。この架空データ、列が６個あり、　 * wpid: 職場ID * id :個人のID * yr :健康診断の年度 * bmi :BMI * sbp :収縮期血圧 * dbp :拡張期血圧

です。ここで、個人を特定するにはwpidとidを組み合わせて行う形なっています。なので、id=1の人はwpid毎に一人ずついるので注意してください。　このデータを利用して、Rmrkdownを書いて、会社全体のレポートを上司に向けて作成してみてください。レポートの内容は、上司はパワーポイントでの発表を好むので、　 * スライド1：タイトル * スライド2：全体を集計した表 * スライド3~：職場別に集計した表

としましょう。どのような集計方法とするかは、皆さんにお任せします。なお、会社の標準的なプレゼン資料のテンプレートは、kadai/bunseki_co_template.pptxの内容であるので、それを反映してください。

それでははじめていきましょう。

まずは、実際にRMarkdownを書く前にどんな表とグラフを表示するか、をここで考えておきます。

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

dat <- read_csv("kadai/data/data.csv")

Rows: 2460 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (6): wpid, id, yr, bmi, sbp, dbp

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

全体を集計した表:

tidyverseだけを利用するなら、

hyou1tidyv <- dat %>% 
  group_by(yr) %>% 
  summarise(
    across(c(bmi,sbp,dbp),
           .fn = list(mn = ~{mean(.)},
                      q25 = ~{quantile(.,0.25)},
                      q75 = ~{quantile(.,0.75)}),
           .names = "{.col}_{.fn}")
  ) %>% 
  mutate(across(!yr, ~{signif(.,3)})) %>% 
  mutate(
    bmi_res = str_glue("{bmi_mn}({bmi_q25}-{bmi_q75})"),
    sbp_res = str_glue("{sbp_mn}({sbp_q25}-{sbp_q75})"),
    dbp_res = str_glue("{dbp_mn}({dbp_q25}-{dbp_q75})")
  ) %>% 
  select(`年度`=yr,
         `BMI` = bmi_res,
         `収縮期血圧` = sbp_res,
         `拡張期血圧` = dbp_res)

こんな感じで表を作成してあげます。

across等を駆使して短めに書いてはいますが、ちょっと面倒です

ここで利用している、　 * signifはという有効数字に丸める関数です。 * str_glueは”{列名}“とすると、その部分を置き換えた文字に置き換えてくれる便利な関数です。

ちょっと難しいかもしれません

ただ、とりあえず、このような形の表を作成できれば、knitr::kableという関数に与えてあげると、自動的にmarkdown形式に変換してくれます。

hyou1tidyv

# A tibble: 3 × 4
   年度 BMI             収縮期血圧    拡張期血圧     
  <dbl> <glue>          <glue>        <glue>         
1  2000 22.4(21.3-23.5) 106(82.1-140) 48.4(24-81.7)  
2  2001 22.4(21.3-23.5) 106(82.1-140) 48.4(24-81.6)  
3  2002 22.4(21.3-23.4) 106(82-140)   48.3(23.9-81.9)

knitr::kable(hyou1tidyv)

年度	BMI	収縮期血圧	拡張期血圧
2000	22.4(21.3-23.5)	106(82.1-140)	48.4(24-81.7)
2001	22.4(21.3-23.5)	106(82.1-140)	48.4(24-81.6)
2002	22.4(21.3-23.4)	106(82-140)	48.3(23.9-81.9)

これを利用してもよいですし、

gtsummaryというパッケージを利用するともっと簡単に集計表を作成できて、

tbl <- gtsummary::tbl_summary(dat,by="yr",include = c("bmi","sbp","dbp"))
tbl

Characteristic	2000, N = 820¹	2001, N = 820¹	2002, N = 820¹
bmi	22.34 (21.33, 23.52)	22.35 (21.29, 23.51)	22.42 (21.28, 23.40)
sbp	92 (82, 140)	92 (82, 140)	92 (82, 140)
dbp	34 (24, 82)	34 (24, 82)	33 (24, 82)
¹ Median (IQR)

こんな感じです。英語表記なので、日本語に直す手間は少しかかりますが、こちらの方が、処理の内容としてはだいぶ楽だと思います。

library(gtsummary)

Warning: package 'gtsummary' was built under R version 4.3.2

tbl2 <- tbl_summary(
  data = dat,
  by="yr",
  include = c("bmi","sbp","dbp"),
  label = list(bmi ~ "BMI", 
               sbp ~ "収縮期血圧",
               dbp ~ "拡張期血圧") 
) %>% 
  modify_header(list(label~"**検査値**")) %>%  
  　　　# 検査値を足しています
  
  modify_spanning_header( #この関数で「年度」を足しています
    list(
      stat_1 ~ "**年度**",
      stat_2 ~ "**年度**",
      stat_3 ~ "**年度**"
    )
  ) %>% 
  
  modify_footnote( # この関数でフットノートを日本語にしています
    list(
      stat_1 ~ "中央値(IQR)",
      stat_2 ~ "中央値(IQR)",
      stat_3 ~ "中央値(IQR)"
    )
  )

tbl2

検査値	年度
検査値	2000, N = 820¹	2001, N = 820¹	2002, N = 820¹
BMI	22.34 (21.33, 23.52)	22.35 (21.29, 23.51)	22.42 (21.28, 23.40)
収縮期血圧	92 (82, 140)	92 (82, 140)	92 (82, 140)
拡張期血圧	34 (24, 82)	34 (24, 82)	33 (24, 82)
¹ 中央値(IQR)

modify~という関数で設定を色々変えています。

stat_1など、表の列名をしらべるには、

gtsummary::show_header_names(tbl)

ℹ As a usage guide, the code below re-creates the current column headers.

modify_header(
  label = "**Characteristic**",
  stat_1 = "**2000**, N = 820",
  stat_2 = "**2001**, N = 820",
  stat_3 = "**2002**, N = 820"
)



Column Name   Column Header      
------------  -------------------
label         **Characteristic** 
stat_1        **2000**, N = 820  
stat_2        **2001**, N = 820  
stat_3        **2002**, N = 820

としてあげます。gtsummaryの詳しい使い方は本コースの対象範囲からはずれるので、ここでは上のようなパッケージや関数を使うと表が作れるんだというくらいの認識でよいと思います。

今回は、この表を利用しましょう。

尚、職場別に集計した結果は、このgtsummary::tbl_summaryに与えるデータをfilterを利用して職場単位に絞り込めばOKです。なので、

dat %>% count(wpid)

# A tibble: 5 × 2
   wpid     n
  <dbl> <int>
1     1   300
2     2   240
3     3   600
4     4   120
5     5  1200

このように、wpidが１から５まで、５つあるので、５回、処理を繰り返す必要があります。

dat %>% 
  filter(wpid==1) %>% 
tbl_summary(
  data = .,
  by="yr",
  include = c("bmi","sbp","dbp"),
  label = list(bmi ~ "BMI", 
               sbp ~ "収縮期血圧",
               dbp ~ "拡張期血圧") 
) %>% 
  modify_header(list(label~"**検査値**")) %>%  
  # 検査値を足しています
  
  modify_spanning_header( #この関数で「年度」を足しています
    list(
      stat_1 ~ "**年度**",
      stat_2 ~ "**年度**",
      stat_3 ~ "**年度**"
    )
  ) %>% 
  
  modify_footnote( # この関数でフットノートを日本語にしています
    list(
      stat_1 ~ "中央値(IQR)",
      stat_2 ~ "中央値(IQR)",
      stat_3 ~ "中央値(IQR)"
    )
  )

検査値	年度
検査値	2000, N = 100¹	2001, N = 100¹	2002, N = 100¹
BMI	22.52 (21.36, 23.70)	22.84 (21.94, 23.87)	22.79 (21.68, 23.62)
収縮期血圧	124.61 (123.63, 125.46)	124.66 (123.83, 125.69)	124.74 (123.83, 125.54)
拡張期血圧	66.76 (65.58, 67.56)	66.83 (65.79, 67.58)	66.86 (65.83, 67.65)
¹ 中央値(IQR)

さすがに面倒なので、関数化しておきましょう。

make_hyou <- function(.data){
  .data %>% 
    tbl_summary(
      data = .,
      by="yr",
      include = c("bmi","sbp","dbp"),
      label = list(bmi ~ "BMI", 
                   sbp ~ "収縮期血圧",
                   dbp ~ "拡張期血圧") 
    ) %>% 
    modify_header(list(label~"**検査値**")) %>%  
    # 検査値を足しています
    
    modify_spanning_header( #この関数で「年度」を足しています
      list(
        stat_1 ~ "**年度**",
        stat_2 ~ "**年度**",
        stat_3 ~ "**年度**"
      )
    ) %>% 
    
    modify_footnote( # この関数でフットノートを日本語にしています
      list(
        stat_1 ~ "中央値(IQR)",
        stat_2 ~ "中央値(IQR)",
        stat_3 ~ "中央値(IQR)"
      )
    )
}

.data を与えてあげて、表ができあがります。

dat %>% filter(wpid==2) %>% make_hyou()

検査値	年度
検査値	2000, N = 80¹	2001, N = 80¹	2002, N = 80¹
BMI	22.44 (21.32, 23.26)	22.19 (21.12, 23.00)	22.13 (21.17, 23.19)
収縮期血圧	93.45 (92.22, 94.52)	93.26 (92.26, 94.11)	93.36 (92.37, 94.31)
拡張期血圧	35.41 (34.44, 36.01)	35.42 (34.32, 36.35)	35.80 (34.19, 36.44)
¹ 中央値(IQR)

この関数を用いて、kadai/kadai1_1.Rmdでレポートが出力できるか見ていきましょう。