パッケージ読み込み

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6      ✔ purrr   0.3.5 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(readxl)
library(tableone)
library(lubridate)
## 
##  次のパッケージを付け加えます: 'lubridate' 
## 
##  以下のオブジェクトは 'package:base' からマスクされています:
## 
##     date, intersect, setdiff, union

データ読み込み

ds1 <- read_excel("FF1_demo.xlsx")

とりあえずsummaryしてみる。

summary(ds1)

なんとみんなcharacter!これは大変だ。

1.年齢階級別退院患者数のグラフを作成する。

ひとまず年齢階級でテーブルを作る。

CreateTableOne(vars="年齢階級",factorVars="年齢階級",data=ds1)
##               
##                Overall    
##   n            356        
##   年齢階級 (%)            
##      0歳代      29 ( 8.1) 
##      10歳代      5 ( 1.4) 
##      20歳代      7 ( 2.0) 
##      30歳代     10 ( 2.8) 
##      40歳代     15 ( 4.2) 
##      50歳代     34 ( 9.6) 
##      60歳代     95 (26.7) 
##      70歳代     85 (23.9) 
##      80歳代     57 (16.0) 
##      90歳以上   19 ( 5.3)

ggplotで棒グラフ作成。模範解答のように横にならべてみる。

ggplot(data=ds1,aes(x=年齢階級))+
geom_bar() + coord_flip()

2.MDC別手術の有無件数を棒グラフで示す

まずはMDCをOPEでグループ分けしてテーブル作成する。

CreateTableOne(vars="MDC", strata="OPE", factorVars="MDC", data=ds1)
##          Stratified by OPE
##           0           1           p      test
##   n       116         185                    
##   MDC (%)                         <0.001     
##      01    25 (21.6)   16 ( 8.6)             
##      02     0 ( 0.0)    1 ( 0.5)             
##      03     0 ( 0.0)    8 ( 4.3)             
##      04    30 (25.9)   12 ( 6.5)             
##      05    13 (11.2)   52 (28.1)             
##      06    11 ( 9.5)   35 (18.9)             
##      07     4 ( 3.4)    2 ( 1.1)             
##      08     0 ( 0.0)    8 ( 4.3)             
##      09     0 ( 0.0)    1 ( 0.5)             
##      10    13 (11.2)    0 ( 0.0)             
##      11    17 (14.7)    4 ( 2.2)             
##      12     0 ( 0.0)    3 ( 1.6)             
##      13     1 ( 0.9)    0 ( 0.0)             
##      15     0 ( 0.0)    8 ( 4.3)             
##      16     2 ( 1.7)   31 (16.8)             
##      17     0 ( 0.0)    2 ( 1.1)             
##      18     0 ( 0.0)    2 ( 1.1)

ggplotを用いて実際にグラフを作成する

ds2 <- ds1 %>% mutate(f.OPE=as.factor(OPE)) 
ggplot(data=ds2,aes(x=MDC,fill=f.OPE)) +
geom_bar(position = "dodge") + #これでサイドバイサイドに
  scale_x_discrete(na.translate = FALSE)  #これでNAを消す

3.2.を積み上げ棒グラフで示す

ds2 <- ds1 %>% mutate(f.OPE=as.factor(OPE)) 
ggplot(data=ds2,aes(x=MDC,fill=f.OPE)) +
geom_bar(position = "fill") +  #position="fill"で割合の棒グラフに
scale_x_discrete(na.translate = FALSE) #何もしなければ縦積み

4.MDC別の平均在院日数を棒グラフにする

まずはサマリーを作成する

ds3 <- ds2 %>% mutate(nu入院日=as.numeric(入院年月日...17),nu退院日=as.numeric(退院年月日))
ds4 <- ds3 %>% mutate(ymd入院日=ymd(nu入院日),ymd退院日=ymd(nu退院日))
ds5 <- ds4 %>% mutate(LOS=(ymd退院日-ymd入院日+1))
ds5.5<-ds5 %>% mutate(nu.LOS=as.numeric(LOS))

CreateTableOne(vars="nu.LOS", strata="MDC", data=ds5.5)
##                     Stratified by MDC
##                      01            02        03            04           
##   n                     41            1          8            42        
##   nu.LOS (mean (SD)) 20.90 (17.57) 6.00 (NA) 20.50 (27.82) 15.60 (12.45)
##                     Stratified by MDC
##                      05            06            07            08           
##   n                     65            46             6             8        
##   nu.LOS (mean (SD)) 11.25 (13.93) 16.02 (19.06) 15.33 (23.91) 24.00 (24.68)
##                     Stratified by MDC
##                      09        10           11            12         
##   n                     1         13           21            3       
##   nu.LOS (mean (SD)) 3.00 (NA) 15.00 (1.41) 15.71 (13.69) 8.67 (0.58)
##                     Stratified by MDC
##                      13         15            16            17         
##   n                      1          8            33            2       
##   nu.LOS (mean (SD)) 12.00 (NA) 18.50 (21.99) 15.21 (15.20) 2.00 (1.41)
##                     Stratified by MDC
##                      18            p   test
##   n                      2                 
##   nu.LOS (mean (SD)) 13.00 (11.31)  NA
ds3 <- ds2 %>% mutate(nu入院日=as.numeric(入院年月日...17),nu退院日=as.numeric(退院年月日))
ds4 <- ds3 %>% mutate(ymd入院日=ymd(nu入院日),ymd退院日=ymd(nu退院日))
ds5 <- ds4 %>% mutate(LOS=(ymd退院日-ymd入院日+1))
ds5 %>% group_by(MDC) %>% summarize(meanLOS=mean(LOS))->ds6
ggplot(data=ds6, aes(x=MDC,y=meanLOS))+
 geom_bar(stat = "identity") #x軸が離散型、y軸が連続変数の棒グラフ(goem_bar)の時にはstat="identity"をいれる

5.退院経路を集計

CreateTableOne(vars = "退院先",factorVars = "退院先",data=ds5)
##             
##              Overall    
##   n          356        
##   退院先 (%)            
##      1       220 (61.8) 
##      2        42 (11.8) 
##      3         6 ( 1.7) 
##      4        45 (12.6) 
##      5         7 ( 2.0) 
##      6         2 ( 0.6) 
##      7         8 ( 2.2) 
##      8        26 ( 7.3)

6.65歳未満および、65歳以上の患者における入院中の大腿骨頚部骨折の発生率を集計

大腿骨頚部骨折のICD10コードはS7200  入院後発症ICD10コードは10まである。これをor=|でつなぐ

ds7 <- ds5 %>% mutate(入院後頚部骨折=case_when(入院後発症ICD10コード1=="S7200"~"yes",入院後発症ICD10コード2=="S7200"~"yes",入院後発症ICD10コード3=="S7200"~"yes",入院後発症ICD10コード4=="S7200"~"yes",入院後発症ICD10コード5=="S7200"~"yes",入院後発症ICD10コード6=="S7200"~"yes",入院後発症ICD10コード7=="S7200"~"yes",入院後発症ICD10コード8=="S7200"~"yes",入院後発症ICD10コード9=="S7200"~"yes",入院後発症ICD10コード10=="S7200"~"yes",TRUE~"no"))
CreateTableOne(vars = "入院後頚部骨折" ,strata = "65歳以上高齢者",factorVars = "入院後頚部骨折",data=ds7)
##                           Stratified by 65歳以上高齢者
##                            0          1          p      test
##   n                        164        192                   
##   入院後頚部骨折 = yes (%)   1 (0.6)    5 (2.6)   0.296

7.発症3日以内に入院した急性脳梗塞患者の入院死亡を集計

脳卒中の発症時期が1の症例のうち主病名が脳梗塞症例、すなわちI630-639を抽出する。さらに、退院転帰が生存or死亡でデータ確認する。退院時転帰として死亡となった場合は6(最も医療資源投入したもので死亡) or 7(6以外のもので死亡)でデータ取得する。

ds5 %>% filter(脳卒中の発症時期==1) %>%
  filter(主傷病ICD10コード=="I630"|主傷病ICD10コード=="I631"|主傷病ICD10コード=="I632"|主傷病ICD10コード=="I633"|主傷病ICD10コード=="I634"|主傷病ICD10コード=="I635"|主傷病ICD10コード=="I636"|主傷病ICD10コード=="I637"|主傷病ICD10コード=="I638"|主傷病ICD10コード=="I639") %>%
  mutate(入院時死亡=case_when(退院時転帰=="6"~"yes",退院時転帰=="7"~"yes",TRUE~"no")) ->ds8
CreateTableOne(vars = "入院時死亡", factorVars = "入院時死亡",data=ds8)
##                       
##                        Overall   
##   n                    18        
##   入院時死亡 = yes (%)  4 (22.2)