パッケージ読み込み
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!これは大変だ。
ひとまず年齢階級でテーブルを作る。
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()
まずは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を消す
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) #何もしなければ縦積み
まずはサマリーを作成する
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"をいれる
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)
大腿骨頚部骨折の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
脳卒中の発症時期が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)