varix main analysis
市田 親正
2023/5/13
必要なlibraryの読み込み
library(tidyverse)
library(readxl)
library(tableone)
library(lubridate)
library(skimr)
library(svglite)
library(haven)
library(summarytools)## Warning in system2("/usr/bin/otool", c("-L", shQuote(DSO)), stdout = TRUE):
## 命令 ''/usr/bin/otool' -L
## '/Library/Frameworks/R.framework/Resources/library/tcltk/libs//tcltk.so''
## の実行は状態 1 を持ちましたlibrary(naniar)
library(devtools)
library(reader)
library(stringr)
library(missForest)
library(mice)
library(geepack)
library(ggplot2)
library(cobalt)
library(tableone)
library(sandwich)
library(survey)
library(WeightIt)
library(MatchIt)
library(twang)
library(naniar)
pacman::p_load(norm2, miceadds, lmtest, car, ROCR, pROC, knitr)データ読み込み
data <- read_excel("varix_icu_ok_2.xlsx")データの確認
str(data)## tibble [790 × 54] (S3: tbl_df/tbl/data.frame)
## $ hosp_id : num [1:790] 1001 1001 1001 1001 1001 ...
## $ pt_id : num [1:790] 1 2 4 5 6 9 12 14 16 17 ...
## $ ad_num : num [1:790] 1 1 1 1 1 1 1 1 1 1 ...
## $ year : num [1:790] 2012 2011 2011 2010 2017 ...
## $ antibio_exposure : num [1:790] 0 0 0 0 0 0 0 0 1 0 ...
## $ antibio_type_a : chr [1:790] "99" "99" "99" "99" ...
## $ antibio_type_b : chr [1:790] "99" "99" "99" "99" ...
## $ antibio_term : num [1:790] 0 0 0 0 0 0 0 0 5 0 ...
## $ age_num : num [1:790] 50 80 44 67 47 73 65 41 38 62 ...
## $ age_cate : num [1:790] 0 2 0 1 0 1 1 0 0 0 ...
## $ sex : chr [1:790] "M" "M" "M" "F" ...
## $ bmi_num : num [1:790] NA 25.3 14.5 NA 21 ...
## $ smoking : num [1:790] 0 0 240 0 270 1000 0 210 0 0 ...
## $ brthel_cate : num [1:790] NA 2 1 NA NA 0 1 1 1 2 ...
## $ child_numl : num [1:790] 11 6 8 NA 11 9 6 8 NA 9 ...
## $ child_score : chr [1:790] "c" "a" "b" NA ...
## $ jcs_all : num [1:790] 0 0 0 0 0 0 0 0 0 3 ...
## $ cci_num : num [1:790] 4 4 4 4 4 4 5 5 4 4 ...
## $ malignancy : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ hd : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ hcc : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ alocohol : num [1:790] 1 0 0 0 1 1 0 1 1 0 ...
## $ past_varix_rup : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ antiplate_user_01: num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ anticoag_user01 : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ nsaids_user01 : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ steroid_user01 : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ ppi_user01 : num [1:790] 1 1 1 1 1 1 1 1 1 1 ...
## $ beta_user01 : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ vaso : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ map_day0 : num [1:790] 0 2 6 2 0 4 0 6 0 8 ...
## $ shock : num [1:790] 1 0 1 1 0 1 0 1 0 NA ...
## $ t_bil : num [1:790] 2.2 1.2 3.4 1.2 7.5 2.2 0.4 1.5 1.8 0.5 ...
## $ ast : num [1:790] 217 31 129 52 112 55 55 56 169 56 ...
## $ alt : num [1:790] 63 22 46 36 53 20 67 12 36 30 ...
## $ wbc : num [1:790] 7400 5000 9100 3900 9800 11100 2900 8700 9300 14600 ...
## $ hb : num [1:790] 6.9 10.8 10.7 6.3 12.2 5.6 5.7 4.9 10.9 3.7 ...
## $ plt : num [1:790] 115 77 162 63 69 124 61 129 90 437 ...
## $ alb : num [1:790] 2.2 3.2 2.9 2.8 2.6 2.3 4 2.2 3.2 2.1 ...
## $ eGFR : num [1:790] 58 57.7 112.4 123.6 89.3 ...
## $ eGFR30 : num [1:790] 0 0 0 0 0 0 1 0 0 0 ...
## $ crp : num [1:790] 0.92 0.29 0.29 0.29 0.1 NA 0.04 0.44 0.09 0.17 ...
## $ pt : num [1:790] 37.8 55 37.8 74.6 28.5 45.9 71.5 43.6 66.7 51.3 ...
## $ pt_cate : num [1:790] 2 1 2 0 2 1 0 1 1 1 ...
## $ aptt : num [1:790] 29.4 29 35.3 30.1 37.8 27.5 25 24.8 20.8 30.9 ...
## $ aptt_cate : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ outcome : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ death_day2 : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ rebleed_end_hemo : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ sbp : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ infection : num [1:790] 0 0 0 0 0 1 0 0 0 0 ...
## $ infection_all : num [1:790] 0 0 0 0 0 1 0 0 0 0 ...
## $ los : num [1:790] 12 7 10 3 6 8 5 8 5 2 ...
## $ cdi : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...データの整理
df <-
data|> #miceパッケージが順序ありの因子型でないと読めないためこれを指定。
mutate(
hosp_id=as.integer(hosp_id),
year=as.integer(year),
pt_id=as.integer(pt_id),
ad_num=as.integer(ad_num),
antibio_exposure=as.integer(antibio_exposure), #geeglmを回すにはfactorでは回らない
antibio_type_a= factor(antibio_type_a, levels = c("a", "b","c","d","e","f","g","h","i","99")),
antibio_type_b= factor(antibio_type_b, levels = c("c","d","f","h","i","99")),
antibio_term= as.integer(antibio_term ),
age_num=as.integer(age_num),
age_cate=factor(age_cate,levels = c("0", "1", "2", "3")),
sex= factor(sex, levels = c("M", "F")),
smoking= as.integer(smoking),
brthel_cate= factor(brthel_cate, levels = c("0", "1", "2")),
child_numl= as.integer(child_numl),
child_score=factor(child_score, levels = c("a", "b", "c")),
jcs_all=factor(jcs_all, levels = c("0", "1", "2","3","10","20","30","100","200","300")),
cci_num=as.integer(cci_num),
malignancy=factor(malignancy),
hd=factor(hd),
hcc=factor(hcc),
alocohol=factor(alocohol),
past_varix_rup=factor(past_varix_rup),
antiplate_user_01=factor(antiplate_user_01),
anticoag_user01=factor(anticoag_user01),
nsaids_user01=factor(nsaids_user01),
steroid_user01=factor(steroid_user01),
ppi_user01=factor(ppi_user01),
beta_user01=factor(beta_user01),
vaso=factor(vaso),
map_day0= as.integer(map_day0),
shock=factor(shock),
pt_cate= factor(pt_cate, levels = c("0", "1", "2")),
aptt_cate=factor(aptt_cate,levels = c("0", "1", "2", "3")),
eGFR30=factor(eGFR30),
outcome=factor(outcome),
death_day2=factor(death_day2),
rebleed_end_hemo=factor(rebleed_end_hemo),
sbp=factor(sbp),
infection=factor(infection),
infection_all=factor(infection_all),
los=as.integer(los),
cdi=factor(cdi)
)整理後の確認
str(df)## tibble [790 × 54] (S3: tbl_df/tbl/data.frame)
## $ hosp_id : int [1:790] 1001 1001 1001 1001 1001 1001 1001 1001 1001 1001 ...
## $ pt_id : int [1:790] 1 2 4 5 6 9 12 14 16 17 ...
## $ ad_num : int [1:790] 1 1 1 1 1 1 1 1 1 1 ...
## $ year : int [1:790] 2012 2011 2011 2010 2017 2010 2013 2014 2016 2019 ...
## $ antibio_exposure : int [1:790] 0 0 0 0 0 0 0 0 1 0 ...
## $ antibio_type_a : Factor w/ 10 levels "a","b","c","d",..: 10 10 10 10 10 10 10 10 4 10 ...
## $ antibio_type_b : Factor w/ 6 levels "c","d","f","h",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ antibio_term : int [1:790] 0 0 0 0 0 0 0 0 5 0 ...
## $ age_num : int [1:790] 50 80 44 67 47 73 65 41 38 62 ...
## $ age_cate : Factor w/ 4 levels "0","1","2","3": 1 3 1 2 1 2 2 1 1 1 ...
## $ sex : Factor w/ 2 levels "M","F": 1 1 1 2 1 1 1 2 1 2 ...
## $ bmi_num : num [1:790] NA 25.3 14.5 NA 21 ...
## $ smoking : int [1:790] 0 0 240 0 270 1000 0 210 0 0 ...
## $ brthel_cate : Factor w/ 3 levels "0","1","2": NA 3 2 NA NA 1 2 2 2 3 ...
## $ child_numl : int [1:790] 11 6 8 NA 11 9 6 8 NA 9 ...
## $ child_score : Factor w/ 3 levels "a","b","c": 3 1 2 NA 3 2 1 2 NA 2 ...
## $ jcs_all : Factor w/ 10 levels "0","1","2","3",..: 1 1 1 1 1 1 1 1 1 4 ...
## $ cci_num : int [1:790] 4 4 4 4 4 4 5 5 4 4 ...
## $ malignancy : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ hd : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ hcc : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ alocohol : Factor w/ 2 levels "0","1": 2 1 1 1 2 2 1 2 2 1 ...
## $ past_varix_rup : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ antiplate_user_01: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ anticoag_user01 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ nsaids_user01 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ steroid_user01 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ ppi_user01 : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
## $ beta_user01 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ vaso : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ map_day0 : int [1:790] 0 2 6 2 0 4 0 6 0 8 ...
## $ shock : Factor w/ 2 levels "0","1": 2 1 2 2 1 2 1 2 1 NA ...
## $ t_bil : num [1:790] 2.2 1.2 3.4 1.2 7.5 2.2 0.4 1.5 1.8 0.5 ...
## $ ast : num [1:790] 217 31 129 52 112 55 55 56 169 56 ...
## $ alt : num [1:790] 63 22 46 36 53 20 67 12 36 30 ...
## $ wbc : num [1:790] 7400 5000 9100 3900 9800 11100 2900 8700 9300 14600 ...
## $ hb : num [1:790] 6.9 10.8 10.7 6.3 12.2 5.6 5.7 4.9 10.9 3.7 ...
## $ plt : num [1:790] 115 77 162 63 69 124 61 129 90 437 ...
## $ alb : num [1:790] 2.2 3.2 2.9 2.8 2.6 2.3 4 2.2 3.2 2.1 ...
## $ eGFR : num [1:790] 58 57.7 112.4 123.6 89.3 ...
## $ eGFR30 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 2 1 1 1 ...
## $ crp : num [1:790] 0.92 0.29 0.29 0.29 0.1 NA 0.04 0.44 0.09 0.17 ...
## $ pt : num [1:790] 37.8 55 37.8 74.6 28.5 45.9 71.5 43.6 66.7 51.3 ...
## $ pt_cate : Factor w/ 3 levels "0","1","2": 3 2 3 1 3 2 1 2 2 2 ...
## $ aptt : num [1:790] 29.4 29 35.3 30.1 37.8 27.5 25 24.8 20.8 30.9 ...
## $ aptt_cate : Factor w/ 4 levels "0","1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ outcome : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ death_day2 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ rebleed_end_hemo : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ sbp : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ infection : Factor w/ 2 levels "0","1": 1 1 1 1 1 2 1 1 1 1 ...
## $ infection_all : Factor w/ 2 levels "0","1": 1 1 1 1 1 2 1 1 1 1 ...
## $ los : int [1:790] 12 7 10 3 6 8 5 8 5 2 ...
## $ cdi : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...確認その2
dfSummary(df)## Data Frame Summary
## df
## Dimensions: 790 x 54
## Duplicates: 0
##
## ------------------------------------------------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ------------------- ----------------------------- --------------------- ---------------------- ---------- ---------
## 1 hosp_id Mean (sd) : 1021.6 (21.9) 44 distinct values : 790 0
## [integer] min < med < max: : (100.0%) (0.0%)
## 1001 < 1017 < 1075 :
## IQR (CV) : 20 (0) : : :
## : : : : : :
##
## 2 pt_id Mean (sd) : 416.6 (243.3) 679 distinct values : . 790 0
## [integer] min < med < max: : : : . : : : : (100.0%) (0.0%)
## 1 < 425.5 < 839 : : : : : : : :
## IQR (CV) : 431.5 (0.6) : : : : : : : : .
## : : : : : : : : :
##
## 3 ad_num Mean (sd) : 1.2 (0.6) 1 : 672 (85.1%) IIIIIIIIIIIIIIIII 790 0
## [integer] min < med < max: 2 : 79 (10.0%) II (100.0%) (0.0%)
## 1 < 1 < 6 3 : 26 ( 3.3%)
## IQR (CV) : 0 (0.5) 4 : 9 ( 1.1%)
## 5 : 2 ( 0.3%)
## 6 : 2 ( 0.3%)
##
## 4 year Mean (sd) : 2015.9 (3.8) 13 distinct values : . : 790 0
## [integer] min < med < max: : : : (100.0%) (0.0%)
## 2010 < 2016 < 2022 : : . : . :
## IQR (CV) : 7 (0) : : : : : : : : : :
## : : : : : : : : : :
##
## 5 antibio_exposure Min : 0 0 : 558 (70.6%) IIIIIIIIIIIIII 790 0
## [integer] Mean : 0.3 1 : 232 (29.4%) IIIII (100.0%) (0.0%)
## Max : 1
##
## 6 antibio_type_a 1. a 4 ( 0.5%) 790 0
## [factor] 2. b 32 ( 4.1%) (100.0%) (0.0%)
## 3. c 51 ( 6.5%) I
## 4. d 104 (13.2%) II
## 5. e 2 ( 0.3%)
## 6. f 12 ( 1.5%)
## 7. g 2 ( 0.3%)
## 8. h 22 ( 2.8%)
## 9. i 3 ( 0.4%)
## 10. 99 558 (70.6%) IIIIIIIIIIIIII
##
## 7 antibio_type_b 1. c 1 ( 0.1%) 790 0
## [factor] 2. d 5 ( 0.6%) (100.0%) (0.0%)
## 3. f 2 ( 0.3%)
## 4. h 4 ( 0.5%)
## 5. i 2 ( 0.3%)
## 6. 99 776 (98.2%) IIIIIIIIIIIIIIIIIII
##
## 8 antibio_term Mean (sd) : 1.3 (3.3) 23 distinct values : 790 0
## [integer] min < med < max: : (100.0%) (0.0%)
## 0 < 0 < 35 :
## IQR (CV) : 1 (2.5) :
## : .
##
## 9 age_num Mean (sd) : 61.1 (13) 65 distinct values : 790 0
## [integer] min < med < max: . : : : (100.0%) (0.0%)
## 24 < 62 < 93 : : : : .
## IQR (CV) : 19 (0.2) : : : : : : :
## . : : : : : : : .
##
## 10 age_cate 1. 0 465 (58.9%) IIIIIIIIIII 790 0
## [factor] 2. 1 190 (24.1%) IIII (100.0%) (0.0%)
## 3. 2 114 (14.4%) II
## 4. 3 21 ( 2.7%)
##
## 11 sex 1. M 598 (75.7%) IIIIIIIIIIIIIII 790 0
## [factor] 2. F 192 (24.3%) IIII (100.0%) (0.0%)
##
## 12 bmi_num Mean (sd) : 23.8 (12.6) 554 distinct values : 706 84
## [numeric] min < med < max: : (89.4%) (10.6%)
## 13.7 < 22.8 < 339.8 :
## IQR (CV) : 5 (0.5) :
## :
##
## 13 smoking Mean (sd) : 257.5 (390) 116 distinct values : 694 96
## [integer] min < med < max: : (87.8%) (12.2%)
## 0 < 0 < 2200 :
## IQR (CV) : 410 (1.5) :
## : : . .
##
## 14 brthel_cate 1. 0 269 (41.2%) IIIIIIII 653 137
## [factor] 2. 1 215 (32.9%) IIIIII (82.7%) (17.3%)
## 3. 2 169 (25.9%) IIIII
##
## 15 child_numl Mean (sd) : 8.3 (2) 11 distinct values : . : 689 101
## [integer] min < med < max: : : : : (87.2%) (12.8%)
## 5 < 8 < 15 : : : : :
## IQR (CV) : 3 (0.2) : : : : : . .
## : : : : : : : .
##
## 16 child_score 1. a 135 (19.0%) III 711 79
## [factor] 2. b 376 (52.9%) IIIIIIIIII (90.0%) (10.0%)
## 3. c 200 (28.1%) IIIII
##
## 17 jcs_all 1. 0 628 (79.5%) IIIIIIIIIIIIIII 790 0
## [factor] 2. 1 93 (11.8%) II (100.0%) (0.0%)
## 3. 2 20 ( 2.5%)
## 4. 3 17 ( 2.2%)
## 5. 10 14 ( 1.8%)
## 6. 20 2 ( 0.3%)
## 7. 30 6 ( 0.8%)
## 8. 100 6 ( 0.8%)
## 9. 200 2 ( 0.3%)
## 10. 300 2 ( 0.3%)
##
## 18 cci_num Mean (sd) : 4.5 (1.3) 11 distinct values : 790 0
## [integer] min < med < max: : (100.0%) (0.0%)
## 3 < 4 < 13 :
## IQR (CV) : 1 (0.3) : :
## : : : .
##
## 19 malignancy 1. 0 696 (88.1%) IIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 94 (11.9%) II (100.0%) (0.0%)
##
## 20 hd 1. 0 779 (98.6%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 11 ( 1.4%) (100.0%) (0.0%)
##
## 21 hcc 1. 0 640 (81.0%) IIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 150 (19.0%) III (100.0%) (0.0%)
##
## 22 alocohol 1. 0 417 (52.8%) IIIIIIIIII 790 0
## [factor] 2. 1 373 (47.2%) IIIIIIIII (100.0%) (0.0%)
##
## 23 past_varix_rup 1. 0 600 (75.9%) IIIIIIIIIIIIIII 790 0
## [factor] 2. 1 190 (24.1%) IIII (100.0%) (0.0%)
##
## 24 antiplate_user_01 1. 0 779 (98.6%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 11 ( 1.4%) (100.0%) (0.0%)
##
## 25 anticoag_user01 1. 0 777 (98.4%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 13 ( 1.6%) (100.0%) (0.0%)
##
## 26 nsaids_user01 1. 0 772 (97.7%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 18 ( 2.3%) (100.0%) (0.0%)
##
## 27 steroid_user01 1. 0 788 (99.7%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 2 ( 0.3%) (100.0%) (0.0%)
##
## 28 ppi_user01 1. 0 91 (11.5%) II 790 0
## [factor] 2. 1 699 (88.5%) IIIIIIIIIIIIIIIII (100.0%) (0.0%)
##
## 29 beta_user01 1. 0 738 (93.4%) IIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 52 ( 6.6%) I (100.0%) (0.0%)
##
## 30 vaso 1. 0 764 (96.7%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 26 ( 3.3%) (100.0%) (0.0%)
##
## 31 map_day0 Mean (sd) : 2.9 (2.9) 13 distinct values : 790 0
## [integer] min < med < max: : (100.0%) (0.0%)
## 0 < 4 < 26 : :
## IQR (CV) : 4 (1) : : .
## : : : .
##
## 32 shock 1. 0 475 (62.0%) IIIIIIIIIIII 766 24
## [factor] 2. 1 291 (38.0%) IIIIIII (97.0%) (3.0%)
##
## 33 t_bil Mean (sd) : 2.1 (1.9) 254 distinct values : 765 25
## [numeric] min < med < max: : (96.8%) (3.2%)
## 0.2 < 1.4 < 13.7 :
## IQR (CV) : 1.6 (0.9) : :
## : : : .
##
## 34 ast Mean (sd) : 77.4 (90.1) 195 distinct values : 773 17
## [numeric] min < med < max: : (97.8%) (2.2%)
## 13 < 49 < 984 :
## IQR (CV) : 56 (1.2) :
## : .
##
## 35 alt Mean (sd) : 39.9 (44.5) 124 distinct values : 773 17
## [numeric] min < med < max: : (97.8%) (2.2%)
## 7 < 28 < 562 :
## IQR (CV) : 25 (1.1) :
## : .
##
## 36 wbc Mean (sd) : 8417.8 (4639.6) 345 distinct values : 776 14
## [numeric] min < med < max: : : (98.2%) (1.8%)
## 1400 < 7400 < 58400 : :
## IQR (CV) : 5027.5 (0.6) : :
## : : :
##
## 37 hb Mean (sd) : 8.8 (2.4) 119 distinct values . . : 776 14
## [numeric] min < med < max: : : : (98.2%) (1.8%)
## 2.1 < 8.6 < 16.5 : : : :
## IQR (CV) : 3.3 (0.3) : : : : : .
## . : : : : : : . .
##
## 38 plt Mean (sd) : 116.5 (69.4) 217 distinct values : 776 14
## [numeric] min < med < max: : (98.2%) (1.8%)
## 18 < 102 < 1073 : .
## IQR (CV) : 68 (0.6) : :
## : : .
##
## 39 alb Mean (sd) : 2.9 (0.6) 35 distinct values : 754 36
## [numeric] min < med < max: : : (95.4%) (4.6%)
## 1.1 < 2.9 < 4.9 : : :
## IQR (CV) : 0.8 (0.2) : : : :
## : : : : : .
##
## 40 eGFR Mean (sd) : 71.6 (30.8) 737 distinct values : : 775 15
## [numeric] min < med < max: : : . (98.1%) (1.9%)
## 4.3 < 69.2 < 196 : : :
## IQR (CV) : 40.6 (0.4) : : : : .
## : : : : : : . .
##
## 41 eGFR30 1. 0 716 (92.4%) IIIIIIIIIIIIIIIIII 775 15
## [factor] 2. 1 59 ( 7.6%) I (98.1%) (1.9%)
##
## 42 crp Mean (sd) : 0.8 (1.4) 360 distinct values : 753 37
## [numeric] min < med < max: : (95.3%) (4.7%)
## 0 < 0.3 < 18.6 :
## IQR (CV) : 0.7 (1.9) :
## : .
##
## 43 pt Mean (sd) : 57.1 (16.7) 383 distinct values . . : 744 46
## [numeric] min < med < max: : : : (94.2%) (5.8%)
## 13.4 < 57.9 < 107.9 : : : : :
## IQR (CV) : 23.8 (0.3) . : : : : : .
## . : : : : : : : . .
##
## 44 pt_cate 1. 0 161 (21.7%) IIII 743 47
## [factor] 2. 1 462 (62.2%) IIIIIIIIIIII (94.1%) (5.9%)
## 3. 2 120 (16.2%) III
##
## 45 aptt Mean (sd) : 32.7 (14.3) 234 distinct values : 700 90
## [numeric] min < med < max: : (88.6%) (11.4%)
## 17.6 < 30.2 < 200 :
## IQR (CV) : 7.6 (0.4) :
## : :
##
## 46 aptt_cate 1. 0 621 (88.7%) IIIIIIIIIIIIIIIII 700 90
## [factor] 2. 1 70 (10.0%) II (88.6%) (11.4%)
## 3. 2 2 ( 0.3%)
## 4. 3 7 ( 1.0%)
##
## 47 outcome 1. 0 724 (91.6%) IIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 66 ( 8.4%) I (100.0%) (0.0%)
##
## 48 death_day2 1. 0 742 (93.9%) IIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 48 ( 6.1%) I (100.0%) (0.0%)
##
## 49 rebleed_end_hemo 1. 0 775 (98.1%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 15 ( 1.9%) (100.0%) (0.0%)
##
## 50 sbp 1. 0 775 (98.1%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 15 ( 1.9%) (100.0%) (0.0%)
##
## 51 infection 1. 0 753 (95.3%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 37 ( 4.7%) (100.0%) (0.0%)
##
## 52 infection_all 1. 0 744 (94.2%) IIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 46 ( 5.8%) I (100.0%) (0.0%)
##
## 53 los Mean (sd) : 13 (15) 58 distinct values : 790 0
## [integer] min < med < max: : (100.0%) (0.0%)
## 2 < 9 < 217 :
## IQR (CV) : 9 (1.2) :
## : .
##
## 54 cdi 1. 0 788 (99.7%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 2 ( 0.3%) (100.0%) (0.0%)
## ------------------------------------------------------------------------------------------------------------------------確認その3
dfSummary(df) %>% view()## Switching method to 'browser'## Output file written: /var/folders/n9/tf_wmwpn3gl2t7l1cz4tqk000000gn/T//RtmpNt1TFV/file5ac766c1c82a.html欠測評価
gg_miss_var(df,show_pct = T)## Warning: The `guide` argument in `scale_*()` cannot be `FALSE`. This was deprecated in
## ggplot2 3.3.4.
## ℹ Please use "none" instead.
## ℹ The deprecated feature was likely used in the naniar package.
## Please report the issue at <]8;;https://github.com/njtierney/naniar/issueshttps://github.com/njtierney/naniar/issues]8;;>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
欠測評価2
vis_miss(df,cluster = T)## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## ℹ Please use `gather()` instead.
## ℹ The deprecated feature was likely used in the visdat package.
## Please report the issue at <]8;;https://github.com/ropensci/visdat/issueshttps://github.com/ropensci/visdat/issues]8;;>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
tableoneアウトカムも含めて作成
col_cont <-
c("antibio_term","child_numl","cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "alb", "crp", "pt", "aptt", "los")
col_fact <-
c("antibio_type_a", "antibio_type_b", "age_cate","sex", "brthel_cate", "malignancy", "child_score", "hd", "hcc", "alocohol", "past_varix_rup", "antiplate_user_01", "anticoag_user01", "nsaids_user01" ,"steroid_user01", "ppi_user01", "beta_user01","shock","pt_cate", "aptt_cate","eGFR30", "outcome", "death_day2", "rebleed_end_hemo", "sbp", "infection", "infection_all", "cdi")
df|> select(c(col_cont, col_fact, "antibio_exposure"))|>
CreateTableOne(vars = c(col_cont, col_fact), strata = "antibio_exposure", factorVars = col_fact) -> tableone2## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(col_cont)
##
## # Now:
## data %>% select(all_of(col_cont))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(col_fact)
##
## # Now:
## data %>% select(all_of(col_fact))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.print(tableone2, smd = T, missing = T, test = T, explain = T) ## Stratified by antibio_exposure
## 0 1 p test
## n 558 232
## antibio_term (mean (SD)) 0.00 (0.00) 4.59 (4.83) <0.001
## child_numl (mean (SD)) 8.30 (1.96) 8.41 (2.03) 0.508
## cci_num (mean (SD)) 4.51 (1.32) 4.52 (1.33) 0.917
## map_day0 (mean (SD)) 2.86 (2.92) 3.13 (2.98) 0.231
## t_bil (mean (SD)) 2.01 (1.88) 2.26 (2.03) 0.098
## ast (mean (SD)) 74.97 (82.14) 83.16 (106.55) 0.248
## alt (mean (SD)) 38.83 (42.49) 42.57 (49.00) 0.285
## wbc (mean (SD)) 8194.22 (4242.15) 8945.32 (5438.02) 0.039
## hb (mean (SD)) 8.69 (2.47) 8.98 (2.39) 0.133
## plt (mean (SD)) 116.49 (61.73) 116.46 (84.84) 0.996
## alb (mean (SD)) 2.91 (0.57) 2.92 (0.64) 0.789
## crp (mean (SD)) 0.76 (1.39) 0.72 (1.57) 0.699
## pt (mean (SD)) 57.00 (16.56) 57.31 (17.01) 0.816
## aptt (mean (SD)) 31.95 (12.69) 34.43 (17.33) 0.035
## los (mean (SD)) 12.92 (13.73) 13.06 (17.84) 0.905
## antibio_type_a (%) <0.001
## a 0 ( 0.0) 4 ( 1.7)
## b 0 ( 0.0) 32 (13.8)
## c 0 ( 0.0) 51 (22.0)
## d 0 ( 0.0) 104 (44.8)
## e 0 ( 0.0) 2 ( 0.9)
## f 0 ( 0.0) 12 ( 5.2)
## g 0 ( 0.0) 2 ( 0.9)
## h 0 ( 0.0) 22 ( 9.5)
## i 0 ( 0.0) 3 ( 1.3)
## 99 558 (100.0) 0 ( 0.0)
## antibio_type_b (%) <0.001
## c 0 ( 0.0) 1 ( 0.4)
## d 0 ( 0.0) 5 ( 2.2)
## f 0 ( 0.0) 2 ( 0.9)
## h 0 ( 0.0) 4 ( 1.7)
## i 0 ( 0.0) 2 ( 0.9)
## 99 558 (100.0) 218 (94.0)
## age_cate (%) 0.155
## 0 322 ( 57.7) 143 (61.6)
## 1 144 ( 25.8) 46 (19.8)
## 2 75 ( 13.4) 39 (16.8)
## 3 17 ( 3.0) 4 ( 1.7)
## sex = F (%) 141 ( 25.3) 51 (22.0) 0.374
## brthel_cate (%) 0.486
## 0 186 ( 41.5) 83 (40.5)
## 1 152 ( 33.9) 63 (30.7)
## 2 110 ( 24.6) 59 (28.8)
## malignancy = 1 (%) 65 ( 11.6) 29 (12.5) 0.829
## child_score (%) 0.776
## a 93 ( 18.8) 42 (19.4)
## b 266 ( 53.7) 110 (50.9)
## c 136 ( 27.5) 64 (29.6)
## hd = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## hcc = 1 (%) 112 ( 20.1) 38 (16.4) 0.269
## alocohol = 1 (%) 246 ( 44.1) 127 (54.7) 0.008
## past_varix_rup = 1 (%) 127 ( 22.8) 63 (27.2) 0.221
## antiplate_user_01 = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## anticoag_user01 = 1 (%) 8 ( 1.4) 5 ( 2.2) 0.675
## nsaids_user01 = 1 (%) 13 ( 2.3) 5 ( 2.2) 1.000
## steroid_user01 = 1 (%) 2 ( 0.4) 0 ( 0.0) 0.892
## ppi_user01 = 1 (%) 486 ( 87.1) 213 (91.8) 0.077
## beta_user01 = 1 (%) 26 ( 4.7) 26 (11.2) 0.001
## shock = 1 (%) 197 ( 36.5) 94 (41.4) 0.236
## pt_cate (%) 0.977
## 0 113 ( 21.7) 48 (21.5)
## 1 324 ( 62.3) 138 (61.9)
## 2 83 ( 16.0) 37 (16.6)
## aptt_cate (%) 0.757
## 0 436 ( 89.3) 185 (87.3)
## 1 47 ( 9.6) 23 (10.8)
## 2 1 ( 0.2) 1 ( 0.5)
## 3 4 ( 0.8) 3 ( 1.4)
## eGFR30 = 1 (%) 40 ( 7.3) 19 ( 8.3) 0.769
## outcome = 1 (%) 44 ( 7.9) 22 ( 9.5) 0.550
## death_day2 = 1 (%) 34 ( 6.1) 14 ( 6.0) 1.000
## rebleed_end_hemo = 1 (%) 9 ( 1.6) 6 ( 2.6) 0.531
## sbp = 1 (%) 10 ( 1.8) 5 ( 2.2) 0.957
## infection = 1 (%) 33 ( 5.9) 4 ( 1.7) 0.019
## infection_all = 1 (%) 38 ( 6.8) 8 ( 3.4) 0.095
## cdi = 1 (%) 1 ( 0.2) 1 ( 0.4) 1.000
## Stratified by antibio_exposure
## SMD Missing
## n
## antibio_term (mean (SD)) 1.345 0.0
## child_numl (mean (SD)) 0.054 12.8
## cci_num (mean (SD)) 0.008 0.0
## map_day0 (mean (SD)) 0.093 0.0
## t_bil (mean (SD)) 0.129 3.2
## ast (mean (SD)) 0.086 2.2
## alt (mean (SD)) 0.082 2.2
## wbc (mean (SD)) 0.154 1.8
## hb (mean (SD)) 0.119 1.8
## plt (mean (SD)) <0.001 1.8
## alb (mean (SD)) 0.021 4.6
## crp (mean (SD)) 0.030 4.7
## pt (mean (SD)) 0.018 5.8
## aptt (mean (SD)) 0.163 11.4
## los (mean (SD)) 0.009 0.0
## antibio_type_a (%) 10.677 0.0
## a
## b
## c
## d
## e
## f
## g
## h
## i
## 99
## antibio_type_b (%) 0.358 0.0
## c
## d
## f
## h
## i
## 99
## age_cate (%) 0.183 0.0
## 0
## 1
## 2
## 3
## sex = F (%) 0.077 0.0
## brthel_cate (%) 0.101 17.3
## 0
## 1
## 2
## malignancy = 1 (%) 0.026 0.0
## child_score (%) 0.058 10.0
## a
## b
## c
## hd = 1 (%) 0.012 0.0
## hcc = 1 (%) 0.096 0.0
## alocohol = 1 (%) 0.214 0.0
## past_varix_rup = 1 (%) 0.102 0.0
## antiplate_user_01 = 1 (%) 0.012 0.0
## anticoag_user01 = 1 (%) 0.054 0.0
## nsaids_user01 = 1 (%) 0.012 0.0
## steroid_user01 = 1 (%) 0.085 0.0
## ppi_user01 = 1 (%) 0.154 0.0
## beta_user01 = 1 (%) 0.244 0.0
## shock = 1 (%) 0.100 3.0
## pt_cate (%) 0.017 5.9
## 0
## 1
## 2
## aptt_cate (%) 0.085 11.4
## 0
## 1
## 2
## 3
## eGFR30 = 1 (%) 0.034 1.9
## outcome = 1 (%) 0.057 0.0
## death_day2 = 1 (%) 0.002 0.0
## rebleed_end_hemo = 1 (%) 0.068 0.0
## sbp = 1 (%) 0.026 0.0
## infection = 1 (%) 0.220 0.0
## infection_all = 1 (%) 0.153 0.0
## cdi = 1 (%) 0.046 0.0連続量をアウトカムも含め可視化
df |> #全体
select(col_cont) |>
pivot_longer(cols = col_cont, names_to = "name", values_to = "value") |>
ggplot()+
geom_histogram(aes(x = value), color = "black")+
facet_wrap(~ name, scales = "free", ncol = 5) +
theme_bw()+
theme(text = element_text(size = 12))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 411 rows containing non-finite values (`stat_bin()`).
素データでのtable1の作成
col_continuous <-
c("child_numl","cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "alb", "crp")
col_factors <-
c("age_cate","sex", "brthel_cate", "malignancy", "child_score", "hd", "hcc", "alocohol", "past_varix_rup", "antiplate_user_01", "anticoag_user01", "nsaids_user01" ,"steroid_user01", "ppi_user01","beta_user01", "vaso", "shock", "pt_cate", "aptt_cate", "eGFR30")
df|> select(c(col_continuous, col_factors, "antibio_exposure"))|>
CreateTableOne(vars = c(col_continuous, col_factors), strata = "antibio_exposure", factorVars = col_factors) -> tableone## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(col_continuous)
##
## # Now:
## data %>% select(all_of(col_continuous))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(col_factors)
##
## # Now:
## data %>% select(all_of(col_factors))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.print(tableone, smd = T, missing = T, test = T, explain = T) ## Stratified by antibio_exposure
## 0 1 p test
## n 558 232
## child_numl (mean (SD)) 8.30 (1.96) 8.41 (2.03) 0.508
## cci_num (mean (SD)) 4.51 (1.32) 4.52 (1.33) 0.917
## map_day0 (mean (SD)) 2.86 (2.92) 3.13 (2.98) 0.231
## t_bil (mean (SD)) 2.01 (1.88) 2.26 (2.03) 0.098
## ast (mean (SD)) 74.97 (82.14) 83.16 (106.55) 0.248
## alt (mean (SD)) 38.83 (42.49) 42.57 (49.00) 0.285
## wbc (mean (SD)) 8194.22 (4242.15) 8945.32 (5438.02) 0.039
## hb (mean (SD)) 8.69 (2.47) 8.98 (2.39) 0.133
## plt (mean (SD)) 116.49 (61.73) 116.46 (84.84) 0.996
## alb (mean (SD)) 2.91 (0.57) 2.92 (0.64) 0.789
## crp (mean (SD)) 0.76 (1.39) 0.72 (1.57) 0.699
## age_cate (%) 0.155
## 0 322 (57.7) 143 (61.6)
## 1 144 (25.8) 46 (19.8)
## 2 75 (13.4) 39 (16.8)
## 3 17 ( 3.0) 4 ( 1.7)
## sex = F (%) 141 (25.3) 51 (22.0) 0.374
## brthel_cate (%) 0.486
## 0 186 (41.5) 83 (40.5)
## 1 152 (33.9) 63 (30.7)
## 2 110 (24.6) 59 (28.8)
## malignancy = 1 (%) 65 (11.6) 29 (12.5) 0.829
## child_score (%) 0.776
## a 93 (18.8) 42 (19.4)
## b 266 (53.7) 110 (50.9)
## c 136 (27.5) 64 (29.6)
## hd = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## hcc = 1 (%) 112 (20.1) 38 (16.4) 0.269
## alocohol = 1 (%) 246 (44.1) 127 (54.7) 0.008
## past_varix_rup = 1 (%) 127 (22.8) 63 (27.2) 0.221
## antiplate_user_01 = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## anticoag_user01 = 1 (%) 8 ( 1.4) 5 ( 2.2) 0.675
## nsaids_user01 = 1 (%) 13 ( 2.3) 5 ( 2.2) 1.000
## steroid_user01 = 1 (%) 2 ( 0.4) 0 ( 0.0) 0.892
## ppi_user01 = 1 (%) 486 (87.1) 213 (91.8) 0.077
## beta_user01 = 1 (%) 26 ( 4.7) 26 (11.2) 0.001
## vaso = 1 (%) 19 ( 3.4) 7 ( 3.0) 0.953
## shock = 1 (%) 197 (36.5) 94 (41.4) 0.236
## pt_cate (%) 0.977
## 0 113 (21.7) 48 (21.5)
## 1 324 (62.3) 138 (61.9)
## 2 83 (16.0) 37 (16.6)
## aptt_cate (%) 0.757
## 0 436 (89.3) 185 (87.3)
## 1 47 ( 9.6) 23 (10.8)
## 2 1 ( 0.2) 1 ( 0.5)
## 3 4 ( 0.8) 3 ( 1.4)
## eGFR30 = 1 (%) 40 ( 7.3) 19 ( 8.3) 0.769
## Stratified by antibio_exposure
## SMD Missing
## n
## child_numl (mean (SD)) 0.054 12.8
## cci_num (mean (SD)) 0.008 0.0
## map_day0 (mean (SD)) 0.093 0.0
## t_bil (mean (SD)) 0.129 3.2
## ast (mean (SD)) 0.086 2.2
## alt (mean (SD)) 0.082 2.2
## wbc (mean (SD)) 0.154 1.8
## hb (mean (SD)) 0.119 1.8
## plt (mean (SD)) <0.001 1.8
## alb (mean (SD)) 0.021 4.6
## crp (mean (SD)) 0.030 4.7
## age_cate (%) 0.183 0.0
## 0
## 1
## 2
## 3
## sex = F (%) 0.077 0.0
## brthel_cate (%) 0.101 17.3
## 0
## 1
## 2
## malignancy = 1 (%) 0.026 0.0
## child_score (%) 0.058 10.0
## a
## b
## c
## hd = 1 (%) 0.012 0.0
## hcc = 1 (%) 0.096 0.0
## alocohol = 1 (%) 0.214 0.0
## past_varix_rup = 1 (%) 0.102 0.0
## antiplate_user_01 = 1 (%) 0.012 0.0
## anticoag_user01 = 1 (%) 0.054 0.0
## nsaids_user01 = 1 (%) 0.012 0.0
## steroid_user01 = 1 (%) 0.085 0.0
## ppi_user01 = 1 (%) 0.154 0.0
## beta_user01 = 1 (%) 0.244 0.0
## vaso = 1 (%) 0.022 0.0
## shock = 1 (%) 0.100 3.0
## pt_cate (%) 0.017 5.9
## 0
## 1
## 2
## aptt_cate (%) 0.085 11.4
## 0
## 1
## 2
## 3
## eGFR30 = 1 (%) 0.034 1.9連続量の可視化
df |> #全体
select(col_continuous) |>
pivot_longer(cols = col_continuous, names_to = "name", values_to = "value") |>
ggplot()+
geom_histogram(aes(x = value), color = "black")+
facet_wrap(~ name, scales = "free", ncol = 5) +
theme_bw()+
theme(text = element_text(size = 12))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 275 rows containing non-finite values (`stat_bin()`).
losのmedianとIQRを求める
library(dplyr)
# antibio_exposureの0.1を基準にデータフレームを群別
grouped_df <- df %>%
mutate(group = ifelse(antibio_exposure >= 0.1, "Group1", "Group0"))
# 各群のlosのmedianとIQRを計算
results <- grouped_df %>%
group_by(group) %>%
summarise(median_los = median(los),
lower_quartile = quantile(los, 0.25),
upper_quartile = quantile(los, 0.75)) %>%
mutate(IQR_los = paste0(lower_quartile, " - ", upper_quartile))
# 不要な列を削除
results <- results[, c("group", "median_los", "IQR_los")]
print(results)## # A tibble: 2 × 3
## group median_los IQR_los
## <chr> <dbl> <chr>
## 1 Group0 9 6 - 15
## 2 Group1 8 5 - 15単変量の解析をみてみる
# Load necessary libraries
library(tidyverse)
library(MASS)##
## 次のパッケージを付け加えます: 'MASS'## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## selectlibrary(broom)
# Univariate logistic regression for outcome, death_day2, rebleed_end_hemo, and sepsis
outcome_model <- glm(outcome ~ antibio_exposure, data = df, family = binomial(link = "logit"))
death_day2_model <- glm(death_day2 ~ antibio_exposure, data = df, family = binomial(link = "logit"))
rebleed_end_hemo_model <- glm(rebleed_end_hemo ~ antibio_exposure, data = df, family = binomial(link = "logit"))
sbp_model <- glm(sbp ~ antibio_exposure, data = df, family = binomial(link = "logit"))
# Univariate negative binomial regression for los
los_model <- glm.nb(los ~ antibio_exposure, data = df)
# Calculate odds ratios
outcome_or <- coef(outcome_model) %>% exp()
death_day2_or <- coef(death_day2_model) %>% exp()
rebleed_end_hemo_or <- coef(rebleed_end_hemo_model) %>% exp()
sbp_or <- coef(sbp_model) %>% exp()
# Calculate rate ratio for LOS
los_rr <- coef(los_model) %>% exp()
# Get odds ratios and 95% confidence intervals for logistic regression models
outcome_ci <- confint(outcome_model, method = "Wald") %>% as.data.frame() %>%
rownames_to_column("term") %>%
mutate(estimate = exp(`2.5 %`), conf.low = exp(`2.5 %`), conf.high = exp(`97.5 %`))## Waiting for profiling to be done...death_day2_ci <- confint(death_day2_model, method = "Wald") %>% as.data.frame() %>%
rownames_to_column("term") %>%
mutate(estimate = exp(`2.5 %`), conf.low = exp(`2.5 %`), conf.high = exp(`97.5 %`))## Waiting for profiling to be done...rebleed_end_hemo_ci <- confint(rebleed_end_hemo_model, method = "Wald") %>% as.data.frame() %>%
rownames_to_column("term") %>%
mutate(estimate = exp(`2.5 %`), conf.low = exp(`2.5 %`), conf.high = exp(`97.5 %`))## Waiting for profiling to be done...sbp_ci <- confint(sbp_model, method = "Wald") %>% as.data.frame() %>%
rownames_to_column("term") %>%
mutate(estimate = exp(`2.5 %`), conf.low = exp(`2.5 %`), conf.high = exp(`97.5 %`))## Waiting for profiling to be done...# Get rate ratios and 95% confidence intervals for negative binomial regression model
los_ci <- confint(los_model, method = "Wald") %>% as.data.frame() %>%
rownames_to_column("term") %>%
mutate(estimate = exp(`2.5 %`), conf.low = exp(`2.5 %`), conf.high = exp(`97.5 %`))## Waiting for profiling to be done...# Print results with rounding
print("Outcome:")## [1] "Outcome:"print(paste("OR:", round(outcome_or["antibio_exposure"], 4), "95% CI:", round(outcome_ci$conf.low[2], 4), "to", round(outcome_ci$conf.high[2], 4)))## [1] "OR: 1.2238 95% CI: 0.7045 to 2.0704"print("Death day 2:")## [1] "Death day 2:"print(paste("OR:", round(death_day2_or["antibio_exposure"], 4), "95% CI:", round(death_day2_ci$conf.low[2], 4), "to", round(death_day2_ci$conf.high[2], 4)))## [1] "OR: 0.9897 95% CI: 0.5052 to 1.8429"print("Rebleed end hemo:")## [1] "Rebleed end hemo:"print(paste("OR:", round(rebleed_end_hemo_or["antibio_exposure"], 4), "95% CI:", round(rebleed_end_hemo_ci$conf.low[2], 4), "to", round(rebleed_end_hemo_ci$conf.high[2], 4)))## [1] "OR: 1.6195 95% CI: 0.5374 to 4.5435"print("SBP:")## [1] "SBP:"print(paste("OR:", round(sbp_or["antibio_exposure"], 4), "95% CI:", round(sbp_ci$conf.low[2], 4), "to", round(sbp_ci$conf.high[2], 4)))## [1] "OR: 1.207 95% CI: 0.3725 to 3.4367"print("LOS:")## [1] "LOS:"print(paste("RR:", round(los_rr["antibio_exposure"], 4), "95% CI:", round(los_ci$conf.low[2], 4), "to", round(los_ci$conf.high[2], 4)))## [1] "RR: 1.0109 95% CI: 0.9 to 1.137"多重代入
#imp1 <- mice(df,
# m=100, #mは代入の結果として作成するデータセットの数
# maxit =200, #maxitは反復回数
# method="pmm",
# printFlag = FALSE,
# seed = 10
# )
#
#save(imp1, file = "imp.Rdata")load("imp.Rdata")plot(imp1)
densityplot(imp1)
ATTでIPTW
# Function to perform IPTW analysis on different outcomes
perform_iptw_analysis <- function(outcome_name) {
gee_results <- list()
summary_tables <- list()
conf_ints <- list()
weighted_data <- list()
ps_data <- list()
for (i in 1:100) {
imputed_data <- complete(imp1, i)
# Fit GEE propensity score model
propensity_model <- geeglm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hd + hcc + alocohol + past_varix_rup + antiplate_user_01 + anticoag_user01 + nsaids_user01 + steroid_user01 + ppi_user01 + beta_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt + aptt,
data = imputed_data,
id = hosp_id,
family = binomial("logit"),
corstr = "exchangeable")
# Calculate propensity scores
ps <- predict(propensity_model, type = "response", newdata = imputed_data)
imputed_data$ps <- ps
# Calculate IPTW and fit the logistic regression model by ATT
stabilized_weights <- with(imputed_data, antibio_exposure + (1 - antibio_exposure) * ps / (1 - ps))
iptw_model <- glm(as.formula(paste0(outcome_name, " ~ antibio_exposure")), data = imputed_data, weights = stabilized_weights, family = binomial("logit"))
# Store the results
gee_results[[i]] <- iptw_model
summary_tables[[i]] <- summary(iptw_model)
conf_ints[[i]] <- confint(iptw_model, level = 0.95, method = "Wald")
# Store weighted data
imputed_data$weight <- stabilized_weights
weighted_data[[i]] <- imputed_data
ps_data[[i]] <- imputed_data[, c("ps", "antibio_exposure")]
}
# Save the results from each imputed dataset
params <- list()
for (i in 1:100) {
params[[i]] <- data.frame(Estimate = coef(gee_results[[i]])[2],
Std.Error = coef(summary_tables[[i]])[2, "Std. Error"],
df = Inf)
}
#Calculate the pooled estimates
Qbar <- mean(sapply(params, function(x) x$Estimate))
se_Qbar <- sqrt(sum(sapply(params, function(x) x$Std.Error^2))/length(params))
#Calculate the pooled odds ratio and confidence intervals
pooled_or <- exp(Qbar)
pooled_lower <- exp(Qbar - 1.96 * se_Qbar)
pooled_upper <- exp(Qbar + 1.96 * se_Qbar)
#Calculate the pooled p-value (two-tailed)
pooled_p <- 2 * (1 - pnorm(abs(Qbar / se_Qbar)))
return(list(OR = pooled_or, Lower_CI = pooled_lower, Upper_CI = pooled_upper, P_value = pooled_p, PS_data = ps_data))
}
#Perform IPTW analysis for each outcome
outcomes <- c("outcome", "death_day2", "rebleed_end_hemo", "sbp", "cdi")
results <- lapply(outcomes, perform_iptw_analysis)## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...#Combine results into a data frame for plotting
forest_plot_data <- data.frame(
Outcome = outcomes,
OR = sapply(results, function(x) x$OR),
Lower_CI = sapply(results, function(x) x$Lower_CI),
Upper_CI = sapply(results, function(x) x$Upper_CI),
P_value = sapply(results, function(x) x$P_value)
)
#Display results as a table
print(forest_plot_data)## Outcome OR Lower_CI Upper_CI P_value
## 1 outcome 1.0990299 0.55574323 2.173426 0.7860638
## 2 death_day2 0.8650817 0.39164495 1.910829 0.7200018
## 3 rebleed_end_hemo 1.3607936 0.35446958 5.224029 0.6535291
## 4 sbp 1.2618382 0.29962549 5.314086 0.7512138
## 5 cdi 3.4103519 0.03452265 336.894775 0.6006042#Create forest plot with the specified order and adjusted p-value position
forest_plot_ordered <- ggplot(forest_plot_data, aes(x = Outcome, y = OR, ymin = Lower_CI, ymax = Upper_CI)) +
geom_pointrange(size = 0.5) +
geom_hline(yintercept = 1, linetype = "dashed") +
scale_y_log10() +
coord_flip() +
theme_bw() +
labs(
x = "Outcome",
y = "Odds Ratio (95% Confidence Interval)"
) +
geom_text(aes(label = sprintf("p = %.3f", P_value)), vjust = -1.0, size = 3) +
scale_x_discrete(limits = c("sbp", "rebleed_end_hemo", "death_day2","outcome"), labels = c( "SBP", "Rebleeding", "In hospital mortality", "Composite outcome"))
#Display forest plot
forest_plot_ordered
#Calculate and display c-statistics for each imputed dataset
c_statistics <- lapply(results[[1]]$PS_data, function(data) {
roc_obj <- roc(data$antibio_exposure, data$ps)
auc(roc_obj)
})## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases#Print c-statistics
print(c_statistics)## [[1]]
## Area under the curve: 0.6381
##
## [[2]]
## Area under the curve: 0.6399
##
## [[3]]
## Area under the curve: 0.6377
##
## [[4]]
## Area under the curve: 0.6351
##
## [[5]]
## Area under the curve: 0.6436
##
## [[6]]
## Area under the curve: 0.6372
##
## [[7]]
## Area under the curve: 0.643
##
## [[8]]
## Area under the curve: 0.6363
##
## [[9]]
## Area under the curve: 0.6358
##
## [[10]]
## Area under the curve: 0.6354
##
## [[11]]
## Area under the curve: 0.6415
##
## [[12]]
## Area under the curve: 0.6414
##
## [[13]]
## Area under the curve: 0.6478
##
## [[14]]
## Area under the curve: 0.6361
##
## [[15]]
## Area under the curve: 0.6459
##
## [[16]]
## Area under the curve: 0.6336
##
## [[17]]
## Area under the curve: 0.6396
##
## [[18]]
## Area under the curve: 0.6388
##
## [[19]]
## Area under the curve: 0.6374
##
## [[20]]
## Area under the curve: 0.6448
##
## [[21]]
## Area under the curve: 0.6401
##
## [[22]]
## Area under the curve: 0.6418
##
## [[23]]
## Area under the curve: 0.6448
##
## [[24]]
## Area under the curve: 0.6404
##
## [[25]]
## Area under the curve: 0.6436
##
## [[26]]
## Area under the curve: 0.6408
##
## [[27]]
## Area under the curve: 0.6404
##
## [[28]]
## Area under the curve: 0.6391
##
## [[29]]
## Area under the curve: 0.6388
##
## [[30]]
## Area under the curve: 0.6421
##
## [[31]]
## Area under the curve: 0.6432
##
## [[32]]
## Area under the curve: 0.637
##
## [[33]]
## Area under the curve: 0.6404
##
## [[34]]
## Area under the curve: 0.6386
##
## [[35]]
## Area under the curve: 0.6378
##
## [[36]]
## Area under the curve: 0.6423
##
## [[37]]
## Area under the curve: 0.6372
##
## [[38]]
## Area under the curve: 0.6422
##
## [[39]]
## Area under the curve: 0.6404
##
## [[40]]
## Area under the curve: 0.6412
##
## [[41]]
## Area under the curve: 0.6373
##
## [[42]]
## Area under the curve: 0.6394
##
## [[43]]
## Area under the curve: 0.6421
##
## [[44]]
## Area under the curve: 0.6437
##
## [[45]]
## Area under the curve: 0.6445
##
## [[46]]
## Area under the curve: 0.6365
##
## [[47]]
## Area under the curve: 0.643
##
## [[48]]
## Area under the curve: 0.6443
##
## [[49]]
## Area under the curve: 0.6414
##
## [[50]]
## Area under the curve: 0.6398
##
## [[51]]
## Area under the curve: 0.6387
##
## [[52]]
## Area under the curve: 0.6421
##
## [[53]]
## Area under the curve: 0.6348
##
## [[54]]
## Area under the curve: 0.6416
##
## [[55]]
## Area under the curve: 0.6363
##
## [[56]]
## Area under the curve: 0.6454
##
## [[57]]
## Area under the curve: 0.6395
##
## [[58]]
## Area under the curve: 0.6421
##
## [[59]]
## Area under the curve: 0.6333
##
## [[60]]
## Area under the curve: 0.6415
##
## [[61]]
## Area under the curve: 0.6348
##
## [[62]]
## Area under the curve: 0.6387
##
## [[63]]
## Area under the curve: 0.6444
##
## [[64]]
## Area under the curve: 0.6348
##
## [[65]]
## Area under the curve: 0.6386
##
## [[66]]
## Area under the curve: 0.642
##
## [[67]]
## Area under the curve: 0.6398
##
## [[68]]
## Area under the curve: 0.6384
##
## [[69]]
## Area under the curve: 0.6461
##
## [[70]]
## Area under the curve: 0.6401
##
## [[71]]
## Area under the curve: 0.6434
##
## [[72]]
## Area under the curve: 0.6403
##
## [[73]]
## Area under the curve: 0.6392
##
## [[74]]
## Area under the curve: 0.6375
##
## [[75]]
## Area under the curve: 0.6444
##
## [[76]]
## Area under the curve: 0.6367
##
## [[77]]
## Area under the curve: 0.6417
##
## [[78]]
## Area under the curve: 0.6414
##
## [[79]]
## Area under the curve: 0.6391
##
## [[80]]
## Area under the curve: 0.6373
##
## [[81]]
## Area under the curve: 0.6404
##
## [[82]]
## Area under the curve: 0.6371
##
## [[83]]
## Area under the curve: 0.6398
##
## [[84]]
## Area under the curve: 0.6353
##
## [[85]]
## Area under the curve: 0.6344
##
## [[86]]
## Area under the curve: 0.6468
##
## [[87]]
## Area under the curve: 0.6416
##
## [[88]]
## Area under the curve: 0.642
##
## [[89]]
## Area under the curve: 0.6399
##
## [[90]]
## Area under the curve: 0.6352
##
## [[91]]
## Area under the curve: 0.6399
##
## [[92]]
## Area under the curve: 0.6369
##
## [[93]]
## Area under the curve: 0.6363
##
## [[94]]
## Area under the curve: 0.6385
##
## [[95]]
## Area under the curve: 0.6366
##
## [[96]]
## Area under the curve: 0.6409
##
## [[97]]
## Area under the curve: 0.6352
##
## [[98]]
## Area under the curve: 0.6405
##
## [[99]]
## Area under the curve: 0.6437
##
## [[100]]
## Area under the curve: 0.641#Plot the propensity score overlap for each imputed dataset
for (i in 1:5) {
ggplot(results[[1]]$PS_data[[i]], aes(x = ps, fill = as.factor(antibio_exposure))) +
geom_histogram(position = "identity", alpha = 0.5, bins = 30) +
labs(
title = paste0("Imputed dataset ", i),
x = "Propensity Score",
y = "Count",
fill = "Antibio Exposure"
) +
theme_bw() +
scale_fill_discrete(labels = c("No Exposure", "Exposure")) +
xlim(0, 1) +
guides(fill = guide_legend(reverse = TRUE)) -> p
print(p)
}
ATTで負の二項分布で解析してみる
# Install and load the MASS package
if (!requireNamespace("MASS", quietly = TRUE)) {
install.packages("MASS")
}
library(MASS)
# LOS
# 本解析: IPTW models on imputed data with Negative Binomial distribution
gee_results <- list()
summary_tables <- list()
conf_ints <- list()
weighted_data <- list()
for (i in 1:100) {
imputed_data <- complete(imp1, i)
# Fit GEE propensity score model
propensity_model <- geeglm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hd + hcc + alocohol + past_varix_rup + antiplate_user_01 + anticoag_user01 + nsaids_user01 + steroid_user01 + ppi_user01 + beta_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt + aptt,
data = imputed_data,
id = hosp_id,
family = binomial("logit"),
corstr = "exchangeable")
# Calculate propensity scores
ps <- predict(propensity_model, type = "response", newdata = imputed_data)
imputed_data$ps <- ps
# Calculate IPTW and fit the negative binomial regression model by ATT
stabilized_weights <- with(imputed_data, antibio_exposure + (1 - antibio_exposure) * ps / (1 - ps))
iptw_model <- glm.nb(los ~ antibio_exposure, data = imputed_data, weights = stabilized_weights)
# Store the results
gee_results[[i]] <- iptw_model
summary_tables[[i]] <- summary(iptw_model)
conf_ints[[i]] <- confint(iptw_model, level = 0.95, method = "Wald")
# Store weighted data
imputed_data$weight <- stabilized_weights
weighted_data[[i]] <- imputed_data
}## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...# Save the results from each imputed dataset
params <- list()
for (i in 1:100) {
params[[i]] <- data.frame(estimate = coef(gee_results[[i]])[2],
std.error = coef(summary_tables[[i]])[2, "Std. Error"],
df = Inf) # Add infinite degrees of freedom
}
# Manually pool the results
Qbar <- mean(sapply(params, function(x) x$estimate))
se_Qbar <- sqrt(sum(sapply(params, function(x) x$std.error^2))/length(params))
# Calculate the pooled rate ratio and confidence intervals
pooled_rr <- exp(Qbar)
pooled_lower <- exp(Qbar - 1.96 * se_Qbar)
pooled_upper <- exp(Qbar + 1.96 * se_Qbar)
# Calculate the pooled z-value
pooled_z <- Qbar / se_Qbar
# Calculate the pooled p-value (two-tailed)
pooled_p <- 2 * (1 - pnorm(abs(pooled_z)))
cat("\nPooled results (IPTW with Negative Binomial distribution)\n")##
## Pooled results (IPTW with Negative Binomial distribution)cat("\nRate Ratio:", pooled_rr, "\n")##
## Rate Ratio: 1.032026cat("\n95% Confidence Interval: (", pooled_lower, ", ", pooled_upper, ")\n")##
## 95% Confidence Interval: ( 0.8854464 , 1.202872 )cat("\nP-value:", pooled_p, "\n")##
## P-value: 0.6866938IPTW後のSMD(post-SMD)
# 1. Calculate SMD for all covariates
calculate_smd <- function(data, weights, col_continuous, col_factors) {
smd_continuous <- sapply(col_continuous, function(x) {
cov_unexposed <- data[data$antibio_exposure == 0, x]
cov_exposed <- data[data$antibio_exposure == 1, x]
w_unexposed <- weights[data$antibio_exposure == 0]
w_exposed <- weights[data$antibio_exposure == 1]
mean_unexposed <- sum(w_unexposed * cov_unexposed) / sum(w_unexposed)
mean_exposed <- sum(w_exposed * cov_exposed) / sum(w_exposed)
sd_unexposed <- sqrt(sum(w_unexposed * (cov_unexposed - mean_unexposed)^2) / sum(w_unexposed))
sd_exposed <- sqrt(sum(w_exposed * (cov_exposed - mean_exposed)^2) / sum(w_exposed))
smd <- (mean_exposed - mean_unexposed) / sqrt((sd_exposed^2 + sd_unexposed^2) / 2)
return(smd)
})
smd_factors <- sapply(col_factors, function(x) {
prop_unexposed <- tapply(weights[data$antibio_exposure == 0], data[data$antibio_exposure == 0, x], sum) / sum(weights[data$antibio_exposure == 0])
prop_exposed <- tapply(weights[data$antibio_exposure == 1], data[data$antibio_exposure == 1, x], sum) / sum(weights[data$antibio_exposure == 1])
common_levels <- intersect(names(prop_unexposed), names(prop_exposed))
prop_unexposed <- prop_unexposed[common_levels]
prop_exposed <- prop_exposed[common_levels]
smd <- sum((prop_exposed - prop_unexposed) / sqrt((prop_exposed * (1 - prop_exposed) + prop_unexposed * (1 - prop_unexposed)) / 2), na.rm = TRUE)
return(smd)
})
return(c(smd_continuous, smd_factors))
}
# 2. Calculate SMD for each imputed dataset
smd_results <- lapply(weighted_data, function(data) {
calculate_smd(data, data$weight, col_continuous, col_factors)
})
# 3. Create a data frame with the mean SMD for each covariate
mean_smd_dataframe <- function(smd_results, col_continuous, col_factors) {
mean_smd_continuous <- sapply(1:length(col_continuous), function(i) {
mean(sapply(smd_results, function(x) x[i]))
})
mean_smd_factors <- sapply((length(col_continuous) + 1):(length(col_continuous) + length(col_factors)), function(i) {
mean(sapply(smd_results, function(x) x[i]))
})
mean_smd_continuous_rounded <- round(mean_smd_continuous, 4)
mean_smd_factors_rounded <- round(mean_smd_factors, 4)
smd_df_continuous <- data.frame(Covariate = col_continuous, Mean_SMD = mean_smd_continuous_rounded)
smd_df_factors <- data.frame(Covariate = col_factors, Mean_SMD = mean_smd_factors_rounded)
smd_df <- rbind(smd_df_continuous, smd_df_factors)
smd_df <- smd_df[order(abs(smd_df$Mean_SMD), decreasing = TRUE),]
return(smd_df)
}
# 4. Display the mean SMD dataframe
mean_smd_df <- mean_smd_dataframe(smd_results, col_continuous, col_factors)
mean_smd_df$Mean_SMD<- round(abs(mean_smd_df$Mean_SMD), 4)
mean_smd_df %>% mutate(Imputation ="After IPTW")->mean_smd_df
mean_smd_df %>% rename(SMD=Mean_SMD)->mean_smd_df
rownames(mean_smd_df) <- NULL
print(mean_smd_df)## Covariate SMD Imputation
## 1 alb 0.1170 After IPTW
## 2 aptt_cate 0.1037 After IPTW
## 3 crp 0.0920 After IPTW
## 4 hb 0.0707 After IPTW
## 5 age_cate 0.0584 After IPTW
## 6 plt 0.0429 After IPTW
## 7 t_bil 0.0421 After IPTW
## 8 map_day0 0.0126 After IPTW
## 9 child_score 0.0098 After IPTW
## 10 child_numl 0.0097 After IPTW
## 11 alt 0.0078 After IPTW
## 12 ast 0.0062 After IPTW
## 13 pt_cate 0.0049 After IPTW
## 14 brthel_cate 0.0034 After IPTW
## 15 wbc 0.0011 After IPTW
## 16 cci_num 0.0007 After IPTW
## 17 sex 0.0000 After IPTW
## 18 malignancy 0.0000 After IPTW
## 19 hd 0.0000 After IPTW
## 20 hcc 0.0000 After IPTW
## 21 alocohol 0.0000 After IPTW
## 22 past_varix_rup 0.0000 After IPTW
## 23 antiplate_user_01 0.0000 After IPTW
## 24 anticoag_user01 0.0000 After IPTW
## 25 nsaids_user01 0.0000 After IPTW
## 26 steroid_user01 0.0000 After IPTW
## 27 ppi_user01 0.0000 After IPTW
## 28 beta_user01 0.0000 After IPTW
## 29 vaso 0.0000 After IPTW
## 30 shock 0.0000 After IPTW
## 31 eGFR30 0.0000 After IPTWpre-smdの算出 tableoneから引っ張る
print(tableone, smd = T, missing = T, test = F, explain = F) |> as.data.frame()|>
mutate(SMD = ifelse(SMD == "<0.001", 0.001, SMD)) |>
write.csv("before2.csv")## Stratified by antibio_exposure
## 0 1 SMD Missing
## n 558 232
## child_numl 8.30 (1.96) 8.41 (2.03) 0.054 12.8
## cci_num 4.51 (1.32) 4.52 (1.33) 0.008 0.0
## map_day0 2.86 (2.92) 3.13 (2.98) 0.093 0.0
## t_bil 2.01 (1.88) 2.26 (2.03) 0.129 3.2
## ast 74.97 (82.14) 83.16 (106.55) 0.086 2.2
## alt 38.83 (42.49) 42.57 (49.00) 0.082 2.2
## wbc 8194.22 (4242.15) 8945.32 (5438.02) 0.154 1.8
## hb 8.69 (2.47) 8.98 (2.39) 0.119 1.8
## plt 116.49 (61.73) 116.46 (84.84) <0.001 1.8
## alb 2.91 (0.57) 2.92 (0.64) 0.021 4.6
## crp 0.76 (1.39) 0.72 (1.57) 0.030 4.7
## age_cate 0.183 0.0
## 0 322 (57.7) 143 (61.6)
## 1 144 (25.8) 46 (19.8)
## 2 75 (13.4) 39 (16.8)
## 3 17 ( 3.0) 4 ( 1.7)
## sex = F 141 (25.3) 51 (22.0) 0.077 0.0
## brthel_cate 0.101 17.3
## 0 186 (41.5) 83 (40.5)
## 1 152 (33.9) 63 (30.7)
## 2 110 (24.6) 59 (28.8)
## malignancy = 1 65 (11.6) 29 (12.5) 0.026 0.0
## child_score 0.058 10.0
## a 93 (18.8) 42 (19.4)
## b 266 (53.7) 110 (50.9)
## c 136 (27.5) 64 (29.6)
## hd = 1 8 ( 1.4) 3 ( 1.3) 0.012 0.0
## hcc = 1 112 (20.1) 38 (16.4) 0.096 0.0
## alocohol = 1 246 (44.1) 127 (54.7) 0.214 0.0
## past_varix_rup = 1 127 (22.8) 63 (27.2) 0.102 0.0
## antiplate_user_01 = 1 8 ( 1.4) 3 ( 1.3) 0.012 0.0
## anticoag_user01 = 1 8 ( 1.4) 5 ( 2.2) 0.054 0.0
## nsaids_user01 = 1 13 ( 2.3) 5 ( 2.2) 0.012 0.0
## steroid_user01 = 1 2 ( 0.4) 0 ( 0.0) 0.085 0.0
## ppi_user01 = 1 486 (87.1) 213 (91.8) 0.154 0.0
## beta_user01 = 1 26 ( 4.7) 26 (11.2) 0.244 0.0
## vaso = 1 19 ( 3.4) 7 ( 3.0) 0.022 0.0
## shock = 1 197 (36.5) 94 (41.4) 0.100 3.0
## pt_cate 0.017 5.9
## 0 113 (21.7) 48 (21.5)
## 1 324 (62.3) 138 (61.9)
## 2 83 (16.0) 37 (16.6)
## aptt_cate 0.085 11.4
## 0 436 (89.3) 185 (87.3)
## 1 47 ( 9.6) 23 (10.8)
## 2 1 ( 0.2) 1 ( 0.5)
## 3 4 ( 0.8) 3 ( 1.4)
## eGFR30 = 1 40 ( 7.3) 19 ( 8.3) 0.034 1.9 before_table <-read.csv("before2.csv", col.names = c("covariate", "non_user", "user", "smd", "missing"))
before_table %>% filter(!is.na(smd))-> before_table
before_table %>% mutate(Imputation ="Before IPTW")->before_table
before_table %>% dplyr::select(covariate,SMD=smd,Imputation)->before_table
# 新しい行名を指定
new_row_names <- c("child_numl", "cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "alb", "crp", "age_cate", "sex", "brthel_cate", "malignancy", "child_score", "hd", "hcc", "alocohol", "past_varix_rup", "antiplate_user_01", "anticoag_user01", "nsaids_user01", "steroid_user01", "ppi_user01", "beta_user01","vaso", "shock", "pt_cate", "aptt_cate", "eGFR30")
# 行名を変更
before_table %>% mutate(covariate, new_row_names)->before_table
before_table %>% dplyr::select(Covariate=new_row_names,SMD,Imputation)->before_table
before_table## Covariate SMD Imputation
## 1 child_numl 0.054 Before IPTW
## 2 cci_num 0.008 Before IPTW
## 3 map_day0 0.093 Before IPTW
## 4 t_bil 0.129 Before IPTW
## 5 ast 0.086 Before IPTW
## 6 alt 0.082 Before IPTW
## 7 wbc 0.154 Before IPTW
## 8 hb 0.119 Before IPTW
## 9 plt 0.001 Before IPTW
## 10 alb 0.021 Before IPTW
## 11 crp 0.030 Before IPTW
## 12 age_cate 0.183 Before IPTW
## 13 sex 0.077 Before IPTW
## 14 brthel_cate 0.101 Before IPTW
## 15 malignancy 0.026 Before IPTW
## 16 child_score 0.058 Before IPTW
## 17 hd 0.012 Before IPTW
## 18 hcc 0.096 Before IPTW
## 19 alocohol 0.214 Before IPTW
## 20 past_varix_rup 0.102 Before IPTW
## 21 antiplate_user_01 0.012 Before IPTW
## 22 anticoag_user01 0.054 Before IPTW
## 23 nsaids_user01 0.012 Before IPTW
## 24 steroid_user01 0.085 Before IPTW
## 25 ppi_user01 0.154 Before IPTW
## 26 beta_user01 0.244 Before IPTW
## 27 vaso 0.022 Before IPTW
## 28 shock 0.100 Before IPTW
## 29 pt_cate 0.017 Before IPTW
## 30 aptt_cate 0.085 Before IPTW
## 31 eGFR30 0.034 Before IPTWpre,postの結合
smd_df_pre_post_main <- rbind(before_table, mean_smd_df)
smd_df_pre_post_main## Covariate SMD Imputation
## 1 child_numl 0.0540 Before IPTW
## 2 cci_num 0.0080 Before IPTW
## 3 map_day0 0.0930 Before IPTW
## 4 t_bil 0.1290 Before IPTW
## 5 ast 0.0860 Before IPTW
## 6 alt 0.0820 Before IPTW
## 7 wbc 0.1540 Before IPTW
## 8 hb 0.1190 Before IPTW
## 9 plt 0.0010 Before IPTW
## 10 alb 0.0210 Before IPTW
## 11 crp 0.0300 Before IPTW
## 12 age_cate 0.1830 Before IPTW
## 13 sex 0.0770 Before IPTW
## 14 brthel_cate 0.1010 Before IPTW
## 15 malignancy 0.0260 Before IPTW
## 16 child_score 0.0580 Before IPTW
## 17 hd 0.0120 Before IPTW
## 18 hcc 0.0960 Before IPTW
## 19 alocohol 0.2140 Before IPTW
## 20 past_varix_rup 0.1020 Before IPTW
## 21 antiplate_user_01 0.0120 Before IPTW
## 22 anticoag_user01 0.0540 Before IPTW
## 23 nsaids_user01 0.0120 Before IPTW
## 24 steroid_user01 0.0850 Before IPTW
## 25 ppi_user01 0.1540 Before IPTW
## 26 beta_user01 0.2440 Before IPTW
## 27 vaso 0.0220 Before IPTW
## 28 shock 0.1000 Before IPTW
## 29 pt_cate 0.0170 Before IPTW
## 30 aptt_cate 0.0850 Before IPTW
## 31 eGFR30 0.0340 Before IPTW
## 32 alb 0.1170 After IPTW
## 33 aptt_cate 0.1037 After IPTW
## 34 crp 0.0920 After IPTW
## 35 hb 0.0707 After IPTW
## 36 age_cate 0.0584 After IPTW
## 37 plt 0.0429 After IPTW
## 38 t_bil 0.0421 After IPTW
## 39 map_day0 0.0126 After IPTW
## 40 child_score 0.0098 After IPTW
## 41 child_numl 0.0097 After IPTW
## 42 alt 0.0078 After IPTW
## 43 ast 0.0062 After IPTW
## 44 pt_cate 0.0049 After IPTW
## 45 brthel_cate 0.0034 After IPTW
## 46 wbc 0.0011 After IPTW
## 47 cci_num 0.0007 After IPTW
## 48 sex 0.0000 After IPTW
## 49 malignancy 0.0000 After IPTW
## 50 hd 0.0000 After IPTW
## 51 hcc 0.0000 After IPTW
## 52 alocohol 0.0000 After IPTW
## 53 past_varix_rup 0.0000 After IPTW
## 54 antiplate_user_01 0.0000 After IPTW
## 55 anticoag_user01 0.0000 After IPTW
## 56 nsaids_user01 0.0000 After IPTW
## 57 steroid_user01 0.0000 After IPTW
## 58 ppi_user01 0.0000 After IPTW
## 59 beta_user01 0.0000 After IPTW
## 60 vaso 0.0000 After IPTW
## 61 shock 0.0000 After IPTW
## 62 eGFR30 0.0000 After IPTWpre,postの結合と名称変更
smd_df_pre_post_main <- rbind(before_table, mean_smd_df)
# Create a mapping of covariate names
covariate_name_mapping <- c(
"age_cate"="Age",
"child_numl" = "Child-Pugh score",
"cci_num" = "Charlson Comorbidity Index",
"map_day0" = "RBC transfusion",
"t_bil" = "Total bilirubin",
"ast" = "Aspartate aminotransferase",
"alt" = "Alanine aminotransferase",
"wbc" = "White blood cell",
"hb" = "Hemoglobin",
"plt" = "Platelet",
"alb" = "Albumin",
"crp" = "C-reactive protein",
"pt_cate" = "Prothrombin time",
"aptt_cate" = "Activated partial thromboplastin time",
"sex" = "Sex",
"brthel_cate" = "Barthel index",
"malignancy" = "Malignant tomor history",
"child_score" = "Child-Pugh classification",
"hd" = "Maintenance hemodialysis",
"hcc" = "Hepatic cancer",
"alocohol" = "Alcohol related disease",
"past_varix_rup" = "Past varix rupture history",
"antiplate_user_01" = "Antiplatelet use",
"anticoag_user01" = "Anticoagulant use",
"nsaids_user01" = "NSAIDs use",
"steroid_user01" = "Corticosteroid use",
"ppi_user01" = "Acid blocker use",
"beta_user01" = "β blocker use",
"vaso" = "Vasopressor",
"shock" = "Shock index > 1",
"eGFR30" = "eGFR < 30"
)
# Replace covariate names with the new names
smd_df_pre_post_main$Covariate <- covariate_name_mapping[smd_df_pre_post_main$Covariate]
smd_df_pre_post_main## Covariate SMD Imputation
## 1 Child-Pugh score 0.0540 Before IPTW
## 2 Charlson Comorbidity Index 0.0080 Before IPTW
## 3 RBC transfusion 0.0930 Before IPTW
## 4 Total bilirubin 0.1290 Before IPTW
## 5 Aspartate aminotransferase 0.0860 Before IPTW
## 6 Alanine aminotransferase 0.0820 Before IPTW
## 7 White blood cell 0.1540 Before IPTW
## 8 Hemoglobin 0.1190 Before IPTW
## 9 Platelet 0.0010 Before IPTW
## 10 Albumin 0.0210 Before IPTW
## 11 C-reactive protein 0.0300 Before IPTW
## 12 Age 0.1830 Before IPTW
## 13 Sex 0.0770 Before IPTW
## 14 Barthel index 0.1010 Before IPTW
## 15 Malignant tomor history 0.0260 Before IPTW
## 16 Child-Pugh classification 0.0580 Before IPTW
## 17 Maintenance hemodialysis 0.0120 Before IPTW
## 18 Hepatic cancer 0.0960 Before IPTW
## 19 Alcohol related disease 0.2140 Before IPTW
## 20 Past varix rupture history 0.1020 Before IPTW
## 21 Antiplatelet use 0.0120 Before IPTW
## 22 Anticoagulant use 0.0540 Before IPTW
## 23 NSAIDs use 0.0120 Before IPTW
## 24 Corticosteroid use 0.0850 Before IPTW
## 25 Acid blocker use 0.1540 Before IPTW
## 26 β blocker use 0.2440 Before IPTW
## 27 Vasopressor 0.0220 Before IPTW
## 28 Shock index > 1 0.1000 Before IPTW
## 29 Prothrombin time 0.0170 Before IPTW
## 30 Activated partial thromboplastin time 0.0850 Before IPTW
## 31 eGFR < 30 0.0340 Before IPTW
## 32 Albumin 0.1170 After IPTW
## 33 Activated partial thromboplastin time 0.1037 After IPTW
## 34 C-reactive protein 0.0920 After IPTW
## 35 Hemoglobin 0.0707 After IPTW
## 36 Age 0.0584 After IPTW
## 37 Platelet 0.0429 After IPTW
## 38 Total bilirubin 0.0421 After IPTW
## 39 RBC transfusion 0.0126 After IPTW
## 40 Child-Pugh classification 0.0098 After IPTW
## 41 Child-Pugh score 0.0097 After IPTW
## 42 Alanine aminotransferase 0.0078 After IPTW
## 43 Aspartate aminotransferase 0.0062 After IPTW
## 44 Prothrombin time 0.0049 After IPTW
## 45 Barthel index 0.0034 After IPTW
## 46 White blood cell 0.0011 After IPTW
## 47 Charlson Comorbidity Index 0.0007 After IPTW
## 48 Sex 0.0000 After IPTW
## 49 Malignant tomor history 0.0000 After IPTW
## 50 Maintenance hemodialysis 0.0000 After IPTW
## 51 Hepatic cancer 0.0000 After IPTW
## 52 Alcohol related disease 0.0000 After IPTW
## 53 Past varix rupture history 0.0000 After IPTW
## 54 Antiplatelet use 0.0000 After IPTW
## 55 Anticoagulant use 0.0000 After IPTW
## 56 NSAIDs use 0.0000 After IPTW
## 57 Corticosteroid use 0.0000 After IPTW
## 58 Acid blocker use 0.0000 After IPTW
## 59 β blocker use 0.0000 After IPTW
## 60 Vasopressor 0.0000 After IPTW
## 61 Shock index > 1 0.0000 After IPTW
## 62 eGFR < 30 0.0000 After IPTWlove plot作成
# Sort the dataframe by absolute SMD of Pre IPTW in descending order
smd_df_pre_post_main_sorted <- smd_df_pre_post_main %>%
filter(Imputation == "Before IPTW") %>%
arrange(SMD) %>%
mutate(Covariate = factor(Covariate, levels = Covariate))
# Merge sorted Pre IPTW data with Post IPTW data
smd_df_pre_post_main_sorted <- smd_df_pre_post_main %>%
filter(Imputation == "After IPTW") %>%
mutate(Covariate = factor(Covariate, levels = smd_df_pre_post_main_sorted$Covariate)) %>%
bind_rows(smd_df_pre_post_main_sorted)
# Create Love plot
love_plot_main <- ggplot(smd_df_pre_post_main_sorted, aes(x = Covariate, y = SMD, color = Imputation)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
coord_flip() +
theme_bw() +
labs(
x = "Covariate",
y = "Absolute Standardized Mean Difference"
) +
guides(color = guide_legend(title = NULL))
# Show the love plot
print(love_plot_main)
tableobe作成の詳細
library(dplyr)
library(tidyr)
# Continuous variables for which we want median and IQR
continuous_vars <- c("child_numl","cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "alb", "crp")
# Function to calculate median and IQR
calc_median_iqr <- function(x) {
q <- quantile(x, c(0.25, 0.5, 0.75), na.rm = TRUE)
paste0("(", round(q[2], 1), ", ", round(q[1], 1), " - ", round(q[3], 1), ")")
}
# Calculate median and IQR for continuous variables
continuous_summary <- df %>%
group_by(antibio_exposure) %>%
summarise(across(all_of(continuous_vars), calc_median_iqr), .groups = "drop") %>%
pivot_longer(cols = continuous_vars, names_to = "Variable", values_to = "Value")
# Calculate percentages for categorical variables
categorical_summary <- df %>%
group_by(antibio_exposure) %>%
summarise(across(all_of(col_factors), ~mean(. == "M", na.rm = TRUE)), .groups = "drop") %>%
pivot_longer(cols = col_factors, names_to = "Variable", values_to = "Value") %>%
mutate(Value = paste0(round(Value * 100, 1), "%"))
# Combine continuous and categorical summaries
tableone <- bind_rows(continuous_summary, categorical_summary) %>%
spread(key = antibio_exposure, value = Value)
# Print table
kable(tableone, row.names = FALSE)| Variable | 0 | 1 |
|---|---|---|
| age_cate | 0% | 0% |
| alb | (2.9, 2.5 - 3.3) | (2.9, 2.5 - 3.3) |
| alocohol | 0% | 0% |
| alt | (27, 19 - 42) | (30.5, 20 - 47) |
| anticoag_user01 | 0% | 0% |
| antiplate_user_01 | 0% | 0% |
| aptt_cate | 0% | 0% |
| ast | (47, 31 - 83) | (54.5, 32.2 - 94.8) |
| beta_user01 | 0% | 0% |
| brthel_cate | 0% | 0% |
| cci_num | (4, 4 - 5) | (4, 4 - 5) |
| child_numl | (8, 7 - 10) | (8, 7 - 10) |
| child_score | 0% | 0% |
| crp | (0.3, 0.1 - 0.8) | (0.3, 0.1 - 0.7) |
| eGFR30 | 0% | 0% |
| hb | (8.5, 6.9 - 10.2) | (9, 7.3 - 10.5) |
| hcc | 0% | 0% |
| hd | 0% | 0% |
| malignancy | 0% | 0% |
| map_day0 | (2.5, 0 - 4) | (4, 0 - 4) |
| nsaids_user01 | 0% | 0% |
| past_varix_rup | 0% | 0% |
| plt | (103, 75 - 144) | (99, 72 - 139) |
| ppi_user01 | 0% | 0% |
| pt_cate | 0% | 0% |
| sex | 74.7% | 78% |
| shock | 0% | 0% |
| steroid_user01 | 0% | 0% |
| t_bil | (1.4, 0.9 - 2.4) | (1.6, 1 - 2.9) |
| vaso | 0% | 0% |
| wbc | (7300, 5200 - 10400) | (7720, 5900 - 10700) |