varix main analysis
市田 親正
2023/4/20
必要なlibraryの読み込み
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorslibrary(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 を持ちました##
## 次のパッケージを付け加えます: 'summarytools'
##
## 以下のオブジェクトは 'package:tibble' からマスクされています:
##
## viewlibrary(naniar)##
## 次のパッケージを付け加えます: 'naniar'
##
## 以下のオブジェクトは 'package:skimr' からマスクされています:
##
## n_completelibrary(devtools)## 要求されたパッケージ usethis をロード中ですlibrary(reader)## 要求されたパッケージ NCmisc をロード中です
##
## 次のパッケージを付け加えます: 'reader'
##
## 以下のオブジェクトは 'package:NCmisc' からマスクされています:
##
## cat.path, get.ext, rmv.extlibrary(stringr)
library(missForest)
library(mice)##
## 次のパッケージを付け加えます: 'mice'
##
## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## filter
##
## 以下のオブジェクトは 'package:base' からマスクされています:
##
## cbind, rbindlibrary(geepack)
library(ggplot2)
library(cobalt)## cobalt (Version 4.4.1, Build Date: 2022-11-03)library(tableone)
library(sandwich)
library(survey)## 要求されたパッケージ grid をロード中です
## 要求されたパッケージ Matrix をロード中です
##
## 次のパッケージを付け加えます: 'Matrix'
##
## 以下のオブジェクトは 'package:tidyr' からマスクされています:
##
## expand, pack, unpack
##
## 要求されたパッケージ survival をロード中です
##
## 次のパッケージを付け加えます: 'survey'
##
## 以下のオブジェクトは 'package:graphics' からマスクされています:
##
## dotchartlibrary(WeightIt)
library(MatchIt)##
## 次のパッケージを付け加えます: 'MatchIt'
##
## 以下のオブジェクトは 'package:cobalt' からマスクされています:
##
## lalondelibrary(twang)## To reproduce results from prior versions of the twang package, please see the version="legacy" option described in the documentation.library(naniar)
pacman::p_load(norm2, miceadds, lmtest, car, ROCR, pROC)データ読み込み
data <- read_excel("varix_icu_ok.xlsx")データの確認
str(data)## tibble [790 × 50] (S3: tbl_df/tbl/data.frame)
## $ hosp_id : num [1:790] 1004 1022 1003 1001 1004 ...
## $ pt_id : num [1:790] 192 500 129 16 217 618 336 63 76 206 ...
## $ ad_num : num [1:790] 1 1 1 1 1 1 1 1 1 1 ...
## $ year : num [1:790] 2014 2013 2014 2016 2015 ...
## $ antibio_exposure : num [1:790] 0 0 0 1 0 1 0 0 0 0 ...
## $ antibio_type_a : chr [1:790] "99" "99" "99" "d" ...
## $ antibio_type_b : chr [1:790] "99" "99" "99" "99" ...
## $ antibio_term : num [1:790] 0 0 0 5 0 1 0 0 0 0 ...
## $ age_num : num [1:790] 86 63 61 38 51 40 69 66 75 66 ...
## $ age_cate : num [1:790] 2 0 0 0 0 0 1 1 2 1 ...
## $ sex : chr [1:790] "F" "M" "M" "M" ...
## $ bmi_num : num [1:790] 17.1 25.7 22.5 27.7 30.2 ...
## $ smoking : num [1:790] 0 860 0 0 380 300 0 0 0 400 ...
## $ brthel_cate : num [1:790] 2 NA 2 1 1 1 0 NA NA 2 ...
## $ child_numl : num [1:790] 8 7 7 NA 8 8 7 9 7 6 ...
## $ child_score : chr [1:790] "b" "b" "b" NA ...
## $ jcs_all : num [1:790] 1 1 0 0 1 0 0 0 0 0 ...
## $ cci_num : num [1:790] 4 3 5 4 4 3 3 4 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] 1 0 0 0 0 0 0 1 0 0 ...
## $ alocohol : num [1:790] 0 0 0 1 0 1 0 0 1 0 ...
## $ past_varix_rup : num [1:790] 0 0 0 0 0 0 0 0 1 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 ...
## $ vaso : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ map_day0 : num [1:790] 0 4 4 0 4 4 4 0 4 4 ...
## $ shock : num [1:790] 0 0 1 0 0 1 1 0 1 0 ...
## $ t_bil : num [1:790] 0.4 0.3 0.69 1.8 0.6 1.3 0.4 0.3 1.19 0.9 ...
## $ ast : num [1:790] 31 17 25 169 37 24 20 30 26 37 ...
## $ alt : num [1:790] 17 19 23 36 21 14 14 33 13 18 ...
## $ wbc : num [1:790] 3500 8400 13700 9300 5000 8100 18000 6000 26300 11600 ...
## $ hb : num [1:790] 9.4 6.1 7.1 10.9 9.1 5 3.5 7.1 10.6 9.4 ...
## $ plt : num [1:790] 93 119 166 90 134 102 291 218 204 165 ...
## $ alb : num [1:790] 2.7 3.3 3.2 3.2 3.1 3 2.5 2.7 3.2 3.7 ...
## $ eGFR : num [1:790] 53.9 52.7 95.5 167.7 52.7 ...
## $ eGFR30 : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ crp : num [1:790] 0.09 0.01 0.144 0.09 0.95 ...
## $ pt : num [1:790] 72.9 77.9 57 66.7 62.7 36.9 85 74 58 60.9 ...
## $ aptt : num [1:790] 17.6 18.9 20.3 20.8 21 21 21.2 21.2 21.4 21.4 ...
## $ 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 0 0 0 0 0 ...
## $ infection_all : num [1:790] 0 0 0 0 0 0 0 0 0 0 ...
## $ los : num [1:790] 16 9 19 5 20 5 9 10 8 11 ...データの整理
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),
vaso=factor(vaso),
map_day0= as.integer(map_day0),
shock=factor(shock),
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)
)整理後の確認
str(df)## tibble [790 × 50] (S3: tbl_df/tbl/data.frame)
## $ hosp_id : int [1:790] 1004 1022 1003 1001 1004 1024 1008 1003 1003 1004 ...
## $ pt_id : int [1:790] 192 500 129 16 217 618 336 63 76 206 ...
## $ ad_num : int [1:790] 1 1 1 1 1 1 1 1 1 1 ...
## $ year : int [1:790] 2014 2013 2014 2016 2015 2014 2020 2014 2017 2015 ...
## $ antibio_exposure : int [1:790] 0 0 0 1 0 1 0 0 0 0 ...
## $ antibio_type_a : Factor w/ 10 levels "a","b","c","d",..: 10 10 10 4 10 4 10 10 10 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 5 0 1 0 0 0 0 ...
## $ age_num : int [1:790] 86 63 61 38 51 40 69 66 75 66 ...
## $ age_cate : Factor w/ 4 levels "0","1","2","3": 3 1 1 1 1 1 2 2 3 2 ...
## $ sex : Factor w/ 2 levels "M","F": 2 1 1 1 2 1 1 1 1 1 ...
## $ bmi_num : num [1:790] 17.1 25.7 22.5 27.7 30.2 ...
## $ smoking : int [1:790] 0 860 0 0 380 300 0 0 0 400 ...
## $ brthel_cate : Factor w/ 3 levels "0","1","2": 3 NA 3 2 2 2 1 NA NA 3 ...
## $ child_numl : int [1:790] 8 7 7 NA 8 8 7 9 7 6 ...
## $ child_score : Factor w/ 3 levels "a","b","c": 2 2 2 NA 2 2 2 2 2 1 ...
## $ jcs_all : Factor w/ 10 levels "0","1","2","3",..: 2 2 1 1 2 1 1 1 1 1 ...
## $ cci_num : int [1:790] 4 3 5 4 4 3 3 4 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": 2 1 1 1 1 1 1 2 1 1 ...
## $ alocohol : Factor w/ 2 levels "0","1": 1 1 1 2 1 2 1 1 2 1 ...
## $ past_varix_rup : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 2 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 ...
## $ vaso : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ map_day0 : int [1:790] 0 4 4 0 4 4 4 0 4 4 ...
## $ shock : Factor w/ 2 levels "0","1": 1 1 2 1 1 2 2 1 2 1 ...
## $ t_bil : num [1:790] 0.4 0.3 0.69 1.8 0.6 1.3 0.4 0.3 1.19 0.9 ...
## $ ast : num [1:790] 31 17 25 169 37 24 20 30 26 37 ...
## $ alt : num [1:790] 17 19 23 36 21 14 14 33 13 18 ...
## $ wbc : num [1:790] 3500 8400 13700 9300 5000 8100 18000 6000 26300 11600 ...
## $ hb : num [1:790] 9.4 6.1 7.1 10.9 9.1 5 3.5 7.1 10.6 9.4 ...
## $ plt : num [1:790] 93 119 166 90 134 102 291 218 204 165 ...
## $ alb : num [1:790] 2.7 3.3 3.2 3.2 3.1 3 2.5 2.7 3.2 3.7 ...
## $ eGFR : num [1:790] 53.9 52.7 95.5 167.7 52.7 ...
## $ eGFR30 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ crp : num [1:790] 0.09 0.01 0.144 0.09 0.95 ...
## $ pt : num [1:790] 72.9 77.9 57 66.7 62.7 36.9 85 74 58 60.9 ...
## $ aptt : num [1:790] 17.6 18.9 20.3 20.8 21 21 21.2 21.2 21.4 21.4 ...
## $ 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 1 1 1 1 1 ...
## $ infection_all : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ los : int [1:790] 16 9 19 5 20 5 9 10 8 11 ...確認その2
dfSummary(df)## Data Frame Summary
## df
## Dimensions: 790 x 50
## 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 133 (16.8%) III
## 4. 3 2 ( 0.3%)
##
## 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 vaso 1. 0 764 (96.7%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 26 ( 3.3%) (100.0%) (0.0%)
##
## 30 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) : : .
## : : : .
##
## 31 shock 1. 0 475 (62.0%) IIIIIIIIIIII 766 24
## [factor] 2. 1 291 (38.0%) IIIIIII (97.0%) (3.0%)
##
## 32 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) : :
## : : : .
##
## 33 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) :
## : .
##
## 34 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) :
## : .
##
## 35 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) : :
## : : :
##
## 36 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) : : : : : .
## . : : : : : : . .
##
## 37 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) : :
## : : .
##
## 38 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) : : : :
## : : : : : .
##
## 39 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) : : : : .
## : : : : : : . .
##
## 40 eGFR30 1. 0 716 (92.4%) IIIIIIIIIIIIIIIIII 775 15
## [factor] 2. 1 59 ( 7.6%) I (98.1%) (1.9%)
##
## 41 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) :
## : .
##
## 42 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) . : : : : : .
## . : : : : : : : . .
##
## 43 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) :
## : :
##
## 44 outcome 1. 0 724 (91.6%) IIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 66 ( 8.4%) I (100.0%) (0.0%)
##
## 45 death_day2 1. 0 742 (93.9%) IIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 48 ( 6.1%) I (100.0%) (0.0%)
##
## 46 rebleed_end_hemo 1. 0 775 (98.1%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 15 ( 1.9%) (100.0%) (0.0%)
##
## 47 sbp 1. 0 775 (98.1%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 15 ( 1.9%) (100.0%) (0.0%)
##
## 48 infection 1. 0 753 (95.3%) IIIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 37 ( 4.7%) (100.0%) (0.0%)
##
## 49 infection_all 1. 0 744 (94.2%) IIIIIIIIIIIIIIIIII 790 0
## [factor] 2. 1 46 ( 5.8%) I (100.0%) (0.0%)
##
## 50 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) :
## : .
## ------------------------------------------------------------------------------------------------------------------------確認その3
dfSummary(df) %>% view()## Switching method to 'browser'## Output file written: /var/folders/n9/tf_wmwpn3gl2t7l1cz4tqk000000gn/T//RtmpTjqcTJ/file3f1214cd8eb0.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", "shock", "eGFR30", "outcome", "death_day2", "rebleed_end_hemo", "sbp", "infection", "infection_all")
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.235
## 0 322 ( 57.7) 143 (61.6)
## 1 144 ( 25.8) 46 (19.8)
## 2 90 ( 16.1) 43 (18.5)
## 3 2 ( 0.4) 0 ( 0.0)
## 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
## shock = 1 (%) 197 ( 36.5) 94 (41.4) 0.236
## 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
## 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.170 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
## shock = 1 (%) 0.100 3.0
## 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連続量をアウトカムも含め可視化
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", "pt", "aptt")
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","vaso", "shock", "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
## 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
## age_cate (%) 0.235
## 0 322 (57.7) 143 (61.6)
## 1 144 (25.8) 46 (19.8)
## 2 90 (16.1) 43 (18.5)
## 3 2 ( 0.4) 0 ( 0.0)
## 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
## vaso = 1 (%) 19 ( 3.4) 7 ( 3.0) 0.953
## shock = 1 (%) 197 (36.5) 94 (41.4) 0.236
## 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
## pt (mean (SD)) 0.018 5.8
## aptt (mean (SD)) 0.163 11.4
## age_cate (%) 0.170 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
## vaso = 1 (%) 0.022 0.0
## shock = 1 (%) 0.100 3.0
## 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 411 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=3, #mは代入の結果として作成するデータセットの数
maxit =3, #maxitは反復回数
method="pmm",
printFlag = FALSE,
seed = 10
)## Warning: Number of logged events: 162plot(imp1)
主要評価解析
# 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:3) {
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 + 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
stabilized_weights <- with(imputed_data, antibio_exposure/(ps) + (1 - antibio_exposure)/(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:3) {
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")
results <- lapply(outcomes, perform_iptw_analysis)## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Waiting for profiling to be done...## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Waiting for profiling to be done...## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Waiting for profiling to be done...## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Waiting for profiling to be done...## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Waiting for profiling to be done...## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Waiting for profiling to be done...## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Waiting for profiling to be done...## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Waiting for profiling to be done...## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Waiting for profiling to be done...## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Waiting for profiling to be done...## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Waiting for profiling to be done...## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Waiting for profiling to be done...## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!# 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.2790394 0.9054478 1.806777 0.1625864
## 2 death_day2 0.9145928 0.6042601 1.384305 0.6728956
## 3 rebleed_end_hemo 1.7454318 0.8805209 3.459921 0.1105962
## 4 sbp 1.6072185 0.8128773 3.177788 0.1724669# 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"))
# 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# Print c-statistics
print(c_statistics)## [[1]]
## Area under the curve: 0.6275
##
## [[2]]
## Area under the curve: 0.6349
##
## [[3]]
## Area under the curve: 0.6227# Plot the propensity score overlap for each imputed dataset
for (i in 1:3) {
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"
) +
theme_bw() +
scale_fill_discrete(labels = c("Non-prophylaxis group", "Prophylaxis group")) +
xlim(0, 1) +
guides(fill = guide_legend(reverse = TRUE, title = NULL)) -> p
print(p)
}## Warning: Removed 4 rows containing missing values (`geom_bar()`).
## Warning: Removed 4 rows containing missing values (`geom_bar()`).
## Warning: Removed 4 rows containing missing values (`geom_bar()`).
forest plotの名称変更
# Install and load the plyr package
if (!requireNamespace("plyr", quietly = TRUE)) {
install.packages("plyr")
}
library(plyr)## ------------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## ------------------------------------------------------------------------------##
## 次のパッケージを付け加えます: 'plyr'## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## 以下のオブジェクトは 'package:purrr' からマスクされています:
##
## compact# Modify labels in the data frame
forest_plot_data$Outcome <- revalue(forest_plot_data$Outcome,
c("outcome" = "Composite outcome",
"death_day2" = "Death in hospitalization",
"rebleed_end_hemo" = "Rebleeding",
"sbp" = "SBP"))
# Reorder forest_plot_data by the desired order of outcomes
forest_plot_data$Outcome <- factor(forest_plot_data$Outcome,
levels = c("Composite outcome",
"Death in hospitalization",
"Rebleeding",
"SBP"))
# Display results as a table
print(forest_plot_data)## Outcome OR Lower_CI Upper_CI P_value
## 1 Composite outcome 1.2790394 0.9054478 1.806777 0.1625864
## 2 Death in hospitalization 0.9145928 0.6042601 1.384305 0.6728956
## 3 Rebleeding 1.7454318 0.8805209 3.459921 0.1105962
## 4 SBP 1.6072185 0.8128773 3.177788 0.1724669# 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) + # Adjust vjust value to change p-value position
scale_x_discrete(limits = c("SBP", "Rebleeding", "Death in hospitalization", "Composite outcome"))
# Display forest plot
forest_plot_ordered
# Save the forest plot as a PNG
ggsave(filename = "IPTW_main_forest.png", plot = forest_plot_ordered, width = 8, height = 6, units = "in", dpi = 300)負の二項分布で解析してみる
# 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:3) {
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 + 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
stabilized_weights <- with(imputed_data, antibio_exposure/(ps) + (1 - antibio_exposure)/(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...# Save the results from each imputed dataset
params <- list()
for (i in 1:3) {
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.03239cat("\n95% Confidence Interval: (", pooled_lower, ", ", pooled_upper, ")\n")##
## 95% Confidence Interval: ( 0.9554436 , 1.115534 )cat("\nP-value:", pooled_p, "\n")##
## P-value: 0.419881サマリー表の作成 補完後IPTW前
# 補完データの1つを選択 (例: 最初の補完データ)
selected_data <- complete(imp1, 1)
# IPTW前のデータのTable 1
unweighted_tableone <- CreateTableOne(vars = c(col_continuous, col_factors), strata = "antibio_exposure",data = selected_data, factorVars = col_factors)
print(unweighted_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.33 (1.93) 8.44 (2.04) 0.476
## 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.25 (2.03) 0.110
## ast (mean (SD)) 74.76 (82.28) 82.66 (106.22) 0.262
## alt (mean (SD)) 38.93 (42.28) 42.34 (48.85) 0.325
## wbc (mean (SD)) 8153.75 (4236.98) 8921.64 (5438.22) 0.034
## hb (mean (SD)) 8.73 (2.49) 8.98 (2.38) 0.180
## plt (mean (SD)) 116.53 (61.66) 116.91 (84.92) 0.943
## alb (mean (SD)) 2.93 (0.59) 2.93 (0.64) 0.903
## crp (mean (SD)) 0.75 (1.36) 0.71 (1.56) 0.726
## pt (mean (SD)) 57.58 (16.65) 57.69 (17.28) 0.931
## aptt (mean (SD)) 31.83 (12.15) 34.20 (16.64) 0.026
## age_cate (%) 0.235
## 0 322 (57.7) 143 (61.6)
## 1 144 (25.8) 46 (19.8)
## 2 90 (16.1) 43 (18.5)
## 3 2 ( 0.4) 0 ( 0.0)
## sex = F (%) 141 (25.3) 51 (22.0) 0.374
## brthel_cate (%) 0.364
## 0 225 (40.3) 92 (39.7)
## 1 190 (34.1) 70 (30.2)
## 2 143 (25.6) 70 (30.2)
## malignancy = 1 (%) 65 (11.6) 29 (12.5) 0.829
## child_score (%) 0.798
## a 102 (18.3) 45 (19.4)
## b 308 (55.2) 122 (52.6)
## c 148 (26.5) 65 (28.0)
## 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
## vaso = 1 (%) 19 ( 3.4) 7 ( 3.0) 0.953
## shock = 1 (%) 203 (36.4) 95 (40.9) 0.260
## eGFR30 = 1 (%) 41 ( 7.3) 19 ( 8.2) 0.795
## Stratified by antibio_exposure
## SMD Missing
## n
## child_numl (mean (SD)) 0.055 0.0
## cci_num (mean (SD)) 0.008 0.0
## map_day0 (mean (SD)) 0.093 0.0
## t_bil (mean (SD)) 0.123 0.0
## ast (mean (SD)) 0.083 0.0
## alt (mean (SD)) 0.075 0.0
## wbc (mean (SD)) 0.158 0.0
## hb (mean (SD)) 0.106 0.0
## plt (mean (SD)) 0.005 0.0
## alb (mean (SD)) 0.009 0.0
## crp (mean (SD)) 0.027 0.0
## pt (mean (SD)) 0.007 0.0
## aptt (mean (SD)) 0.163 0.0
## age_cate (%) 0.170 0.0
## 0
## 1
## 2
## 3
## sex = F (%) 0.077 0.0
## brthel_cate (%) 0.111 0.0
## 0
## 1
## 2
## malignancy = 1 (%) 0.026 0.0
## child_score (%) 0.052 0.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
## vaso = 1 (%) 0.022 0.0
## shock = 1 (%) 0.094 0.0
## eGFR30 = 1 (%) 0.031 0.0補完後IPTW後の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)
cat("Mean SMD for all covariates (rounded to 4 decimal places):\n")## Mean SMD for all covariates (rounded to 4 decimal places):print(mean_smd_df)## Covariate Mean_SMD
## 26 steroid_user01 0.0712
## 1 child_numl 0.0408
## 13 aptt 0.0397
## 11 crp 0.0369
## 10 alb -0.0329
## 12 pt -0.0251
## 2 cci_num 0.0241
## 4 t_bil 0.0230
## 7 wbc 0.0171
## 3 map_day0 0.0157
## 5 ast 0.0098
## 8 hb -0.0081
## 9 plt -0.0073
## 14 age_cate 0.0044
## 18 child_score 0.0041
## 16 brthel_cate 0.0015
## 6 alt 0.0013
## 15 sex 0.0000
## 17 malignancy 0.0000
## 19 hd 0.0000
## 20 hcc 0.0000
## 21 alocohol 0.0000
## 22 past_varix_rup 0.0000
## 23 antiplate_user_01 0.0000
## 24 anticoag_user01 0.0000
## 25 nsaids_user01 0.0000
## 27 ppi_user01 0.0000
## 28 vaso 0.0000
## 29 shock 0.0000
## 30 eGFR30 0.0000love plotで比較
# Calculate pre-imputation SMD
pre_smd_continuous <- sapply(col_continuous, function(x) {
cov_unexposed <- data[data$antibio_exposure == 0, x, drop = FALSE]
cov_exposed <- data[data$antibio_exposure == 1, x, drop = FALSE]
mean_unexposed <- mean(cov_unexposed[[x]], na.rm = TRUE)
mean_exposed <- mean(cov_exposed[[x]], na.rm = TRUE)
sd_unexposed <- sd(cov_unexposed[[x]], na.rm = TRUE)
sd_exposed <- sd(cov_exposed[[x]], na.rm = TRUE)
smd <- (mean_exposed - mean_unexposed) / sqrt((sd_exposed^2 + sd_unexposed^2) / 2)
return(smd)
})
pre_smd_factors <- calculate_smd(data, rep(1, nrow(data)), c(), col_factors)
pre_smd <- c(pre_smd_continuous, pre_smd_factors)
# Calculate the absolute SMD
abs_pre_smd_continuous <- abs(pre_smd_continuous)
abs_pre_smd_factors <- abs(unlist(pre_smd_factors))
abs_mean_smd <- abs(mean_smd_df$Mean_SMD)
# Combine pre- and post-imputation SMDs
smd_df_pre_post <- data.frame(
Covariate = rep(c(col_continuous, col_factors), 2),
SMD = c(abs_pre_smd_continuous, abs_pre_smd_factors, abs_mean_smd),
Imputation = rep(c("Pre IPTW", "Post IPTW"), each = length(col_continuous) + length(col_factors))
)plot作成のためのpreデータ作成
# パッケージの読み込みは不要
# データの読み込み
pre_table <- read.csv("before.csv", col.names = c("covariate", "non_user", "user", "smd", "missing"))
# NAを削除
pre_table <- pre_table[!is.na(pre_table$smd),]
# Imputation列を追加
pre_table$Imputation <- "Before IPTW"
# 新しい行名を指定
new_row_names <- c("child_numl", "cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "alb", "crp", "pt", "aptt", "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", "vaso", "shock", "eGFR30")
# 新しい行名のマッピングデータフレームを作成
mapping <- data.frame(covariate = pre_table$covariate, new_covariate = new_row_names)
# covariate列の名前を変更
pre_table <- merge(pre_table, mapping, by = "covariate", all.x = TRUE)
pre_table$covariate <- NULL
colnames(pre_table)[colnames(pre_table) == "new_covariate"] <- "covariate"
# covariate, smd, Imputation列だけを選択
pre_table <- pre_table[, c("covariate", "smd", "Imputation")]
pre_table## covariate smd Imputation
## 1 age_cate 0.170 Before IPTW
## 2 alb 0.021 Before IPTW
## 3 alocohol 0.214 Before IPTW
## 4 alt 0.082 Before IPTW
## 5 anticoag_user01 0.054 Before IPTW
## 6 antiplate_user_01 0.012 Before IPTW
## 7 aptt 0.163 Before IPTW
## 8 ast 0.086 Before IPTW
## 9 brthel_cate 0.101 Before IPTW
## 10 cci_num 0.008 Before IPTW
## 11 child_numl 0.054 Before IPTW
## 12 child_score 0.058 Before IPTW
## 13 crp 0.030 Before IPTW
## 14 eGFR30 0.034 Before IPTW
## 15 hb 0.119 Before IPTW
## 16 hcc 0.096 Before IPTW
## 17 hd 0.012 Before IPTW
## 18 malignancy 0.026 Before IPTW
## 19 map_day0 0.093 Before IPTW
## 20 nsaids_user01 0.012 Before IPTW
## 21 past_varix_rup 0.102 Before IPTW
## 22 plt 0.001 Before IPTW
## 23 ppi_user01 0.154 Before IPTW
## 24 pt 0.018 Before IPTW
## 25 sex 0.077 Before IPTW
## 26 shock 0.100 Before IPTW
## 27 steroid_user01 0.085 Before IPTW
## 28 t_bil 0.129 Before IPTW
## 29 vaso 0.022 Before IPTW
## 30 wbc 0.154 Before IPTWpostSMDも計算する。
mean_smd_df %>% mutate(Imputation = "After IPTW") ->mean_smd_df
mean_smd_df <- mean_smd_df[, c("Covariate", "Mean_SMD", "Imputation")]
colnames(mean_smd_df) <- c("covariate", "smd", "Imputation")
# smd列の値を絶対値に変換
mean_smd_df$smd <- abs(mean_smd_df$smd)
mean_smd_df## covariate smd Imputation
## 26 steroid_user01 0.0712 After IPTW
## 1 child_numl 0.0408 After IPTW
## 13 aptt 0.0397 After IPTW
## 11 crp 0.0369 After IPTW
## 10 alb 0.0329 After IPTW
## 12 pt 0.0251 After IPTW
## 2 cci_num 0.0241 After IPTW
## 4 t_bil 0.0230 After IPTW
## 7 wbc 0.0171 After IPTW
## 3 map_day0 0.0157 After IPTW
## 5 ast 0.0098 After IPTW
## 8 hb 0.0081 After IPTW
## 9 plt 0.0073 After IPTW
## 14 age_cate 0.0044 After IPTW
## 18 child_score 0.0041 After IPTW
## 16 brthel_cate 0.0015 After IPTW
## 6 alt 0.0013 After IPTW
## 15 sex 0.0000 After IPTW
## 17 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
## 27 ppi_user01 0.0000 After IPTW
## 28 vaso 0.0000 After IPTW
## 29 shock 0.0000 After IPTW
## 30 eGFR30 0.0000 After IPTWpre,postの結合
smd_df_pre_post_main <- rbind(pre_table, mean_smd_df)
smd_df_pre_post_main %>% arrange(smd)## covariate smd Imputation
## 1 sex 0.0000 After IPTW
## 2 malignancy 0.0000 After IPTW
## 3 hd 0.0000 After IPTW
## 4 hcc 0.0000 After IPTW
## 5 alocohol 0.0000 After IPTW
## 6 past_varix_rup 0.0000 After IPTW
## 7 antiplate_user_01 0.0000 After IPTW
## 8 anticoag_user01 0.0000 After IPTW
## 9 nsaids_user01 0.0000 After IPTW
## 10 ppi_user01 0.0000 After IPTW
## 11 vaso 0.0000 After IPTW
## 12 shock 0.0000 After IPTW
## 13 eGFR30 0.0000 After IPTW
## 14 plt 0.0010 Before IPTW
## 15 alt 0.0013 After IPTW
## 16 brthel_cate 0.0015 After IPTW
## 17 child_score 0.0041 After IPTW
## 18 age_cate 0.0044 After IPTW
## 19 plt 0.0073 After IPTW
## 20 cci_num 0.0080 Before IPTW
## 21 hb 0.0081 After IPTW
## 22 ast 0.0098 After IPTW
## 23 antiplate_user_01 0.0120 Before IPTW
## 24 hd 0.0120 Before IPTW
## 25 nsaids_user01 0.0120 Before IPTW
## 26 map_day0 0.0157 After IPTW
## 27 wbc 0.0171 After IPTW
## 28 pt 0.0180 Before IPTW
## 29 alb 0.0210 Before IPTW
## 30 vaso 0.0220 Before IPTW
## 31 t_bil 0.0230 After IPTW
## 32 cci_num 0.0241 After IPTW
## 33 pt 0.0251 After IPTW
## 34 malignancy 0.0260 Before IPTW
## 35 crp 0.0300 Before IPTW
## 36 alb 0.0329 After IPTW
## 37 eGFR30 0.0340 Before IPTW
## 38 crp 0.0369 After IPTW
## 39 aptt 0.0397 After IPTW
## 40 child_numl 0.0408 After IPTW
## 41 anticoag_user01 0.0540 Before IPTW
## 42 child_numl 0.0540 Before IPTW
## 43 child_score 0.0580 Before IPTW
## 44 steroid_user01 0.0712 After IPTW
## 45 sex 0.0770 Before IPTW
## 46 alt 0.0820 Before IPTW
## 47 steroid_user01 0.0850 Before IPTW
## 48 ast 0.0860 Before IPTW
## 49 map_day0 0.0930 Before IPTW
## 50 hcc 0.0960 Before IPTW
## 51 shock 0.1000 Before IPTW
## 52 brthel_cate 0.1010 Before IPTW
## 53 past_varix_rup 0.1020 Before IPTW
## 54 hb 0.1190 Before IPTW
## 55 t_bil 0.1290 Before IPTW
## 56 ppi_user01 0.1540 Before IPTW
## 57 wbc 0.1540 Before IPTW
## 58 aptt 0.1630 Before IPTW
## 59 age_cate 0.1700 Before IPTW
## 60 alocohol 0.2140 Before IPTWlove plot作成
# Assuming smd_df_pre_post_main is your dataframe with columns:
# covariate, smd, Imputation (with values "Pre IPTW" and "Post IPTW")
# 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" = "Prothrombin time",
"aptt" = "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",
"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]
# 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(abs(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))
# Save the plot
ggsave("love_plot_smd_before_after.png", width = 10, height = 8, dpi = 300)
# Show the love plot
print(love_plot_main)
素データでの背景サマリーの作成
# Rename columns
col_continuous <- c("Child-Pugh score", "Charlson Comorbidity Index", "RBC transfusion", "Total bilirubin", "Aspartate aminotransferase", "Alanine aminotransferase", "White blood cell", "Hemoglobin", "Platelet", "Albumin", "C-reactive protein", "Prothrombin time", "Activated partial thromboplastin time")
col_factors <- c("Age", "Sex", "Barthel index", "Child-Pugh classification", "Maintenance hemodialysis", "Hepatic cancer", "Malignant tomor history", "Alcohol related disease", "Past varix rupture history", "Antiplatelet use", "Anticoagulant use", "NSAIDs use", "Corticosteroid use", "Acid blocker use", "Vasopressor","Shock index > 1", "eGFR < 30")
# Update column names in the dataframe
colnames(df)[colnames(df) %in% c("child_numl", "cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "alb", "crp", "pt", "aptt")] <- col_continuous
colnames(df)[colnames(df) %in% c("age_cate", "sex", "brthel_cate", "child_score", "hd", "hcc", "malignancy", "alocohol", "past_varix_rup", "antiplate_user_01", "anticoag_user01", "nsaids_user01", "steroid_user01", "ppi_user01", "vaso", "shock", "eGFR30")] <- col_factors
# Change order of columns
ordered_colnames <- c("Age", "Sex", "Barthel index", "Child-Pugh classification", "Child-Pugh score", "Charlson Comorbidity Index", "Maintenance hemodialysis", "Hepatic cancer", "Malignant tomor history", "Alcohol related disease", "Past varix rupture history", "Antiplatelet use", "Anticoagulant use", "NSAIDs use", "Corticosteroid use", "Acid blocker use","Vasopressor", "RBC transfusion", "Total bilirubin", "Aspartate aminotransferase", "Alanine aminotransferase", "Albumin", "White blood cell", "Hemoglobin", "Platelet", "eGFR < 30", "C-reactive protein", "Prothrombin time", "Activated partial thromboplastin time", "Shock index > 1")
# Create table
tableone <- CreateTableOne(vars = c(ordered_colnames, "antibio_exposure"), strata = "antibio_exposure", data = df, factorVars = col_factors)
# Print table
print(tableone, smd = T, missing = T, test = T, explain = T)## Stratified by antibio_exposure
## 0
## n 558
## Age (%)
## 0 322 (57.7)
## 1 144 (25.8)
## 2 90 (16.1)
## 3 2 ( 0.4)
## Sex = F (%) 141 (25.3)
## Barthel index (%)
## 0 186 (41.5)
## 1 152 (33.9)
## 2 110 (24.6)
## Child-Pugh classification (%)
## a 93 (18.8)
## b 266 (53.7)
## c 136 (27.5)
## Child-Pugh score (mean (SD)) 8.30 (1.96)
## Charlson Comorbidity Index (mean (SD)) 4.51 (1.32)
## Maintenance hemodialysis = 1 (%) 65 (11.6)
## Hepatic cancer = 1 (%) 8 ( 1.4)
## Malignant tomor history = 1 (%) 112 (20.1)
## Alcohol related disease = 1 (%) 246 (44.1)
## Past varix rupture history = 1 (%) 127 (22.8)
## Antiplatelet use = 1 (%) 8 ( 1.4)
## Anticoagulant use = 1 (%) 8 ( 1.4)
## NSAIDs use = 1 (%) 13 ( 2.3)
## Corticosteroid use = 1 (%) 2 ( 0.4)
## Acid blocker use = 1 (%) 486 (87.1)
## Vasopressor = 1 (%) 19 ( 3.4)
## RBC transfusion (mean (SD)) 2.86 (2.92)
## Total bilirubin (mean (SD)) 2.01 (1.88)
## Aspartate aminotransferase (mean (SD)) 74.97 (82.14)
## Alanine aminotransferase (mean (SD)) 38.83 (42.49)
## Albumin (mean (SD)) 2.91 (0.57)
## White blood cell (mean (SD)) 8194.22 (4242.15)
## Hemoglobin (mean (SD)) 8.69 (2.47)
## Platelet (mean (SD)) 116.49 (61.73)
## eGFR < 30 = 1 (%) 40 ( 7.3)
## C-reactive protein (mean (SD)) 0.76 (1.39)
## Prothrombin time (mean (SD)) 57.00 (16.56)
## Activated partial thromboplastin time (mean (SD)) 31.95 (12.69)
## Shock index > 1 = 1 (%) 197 (36.5)
## antibio_exposure (mean (SD)) 0.00 (0.00)
## Stratified by antibio_exposure
## 1 p
## n 232
## Age (%) 0.235
## 0 143 (61.6)
## 1 46 (19.8)
## 2 43 (18.5)
## 3 0 ( 0.0)
## Sex = F (%) 51 (22.0) 0.374
## Barthel index (%) 0.486
## 0 83 (40.5)
## 1 63 (30.7)
## 2 59 (28.8)
## Child-Pugh classification (%) 0.776
## a 42 (19.4)
## b 110 (50.9)
## c 64 (29.6)
## Child-Pugh score (mean (SD)) 8.41 (2.03) 0.508
## Charlson Comorbidity Index (mean (SD)) 4.52 (1.33) 0.917
## Maintenance hemodialysis = 1 (%) 29 (12.5) 0.829
## Hepatic cancer = 1 (%) 3 ( 1.3) 1.000
## Malignant tomor history = 1 (%) 38 (16.4) 0.269
## Alcohol related disease = 1 (%) 127 (54.7) 0.008
## Past varix rupture history = 1 (%) 63 (27.2) 0.221
## Antiplatelet use = 1 (%) 3 ( 1.3) 1.000
## Anticoagulant use = 1 (%) 5 ( 2.2) 0.675
## NSAIDs use = 1 (%) 5 ( 2.2) 1.000
## Corticosteroid use = 1 (%) 0 ( 0.0) 0.892
## Acid blocker use = 1 (%) 213 (91.8) 0.077
## Vasopressor = 1 (%) 7 ( 3.0) 0.953
## RBC transfusion (mean (SD)) 3.13 (2.98) 0.231
## Total bilirubin (mean (SD)) 2.26 (2.03) 0.098
## Aspartate aminotransferase (mean (SD)) 83.16 (106.55) 0.248
## Alanine aminotransferase (mean (SD)) 42.57 (49.00) 0.285
## Albumin (mean (SD)) 2.92 (0.64) 0.789
## White blood cell (mean (SD)) 8945.32 (5438.02) 0.039
## Hemoglobin (mean (SD)) 8.98 (2.39) 0.133
## Platelet (mean (SD)) 116.46 (84.84) 0.996
## eGFR < 30 = 1 (%) 19 ( 8.3) 0.769
## C-reactive protein (mean (SD)) 0.72 (1.57) 0.699
## Prothrombin time (mean (SD)) 57.31 (17.01) 0.816
## Activated partial thromboplastin time (mean (SD)) 34.43 (17.33) 0.035
## Shock index > 1 = 1 (%) 94 (41.4) 0.236
## antibio_exposure (mean (SD)) 1.00 (0.00) <0.001
## Stratified by antibio_exposure
## test SMD Missing
## n
## Age (%) 0.170 0.0
## 0
## 1
## 2
## 3
## Sex = F (%) 0.077 0.0
## Barthel index (%) 0.101 17.3
## 0
## 1
## 2
## Child-Pugh classification (%) 0.058 10.0
## a
## b
## c
## Child-Pugh score (mean (SD)) 0.054 12.8
## Charlson Comorbidity Index (mean (SD)) 0.008 0.0
## Maintenance hemodialysis = 1 (%) 0.026 0.0
## Hepatic cancer = 1 (%) 0.012 0.0
## Malignant tomor history = 1 (%) 0.096 0.0
## Alcohol related disease = 1 (%) 0.214 0.0
## Past varix rupture history = 1 (%) 0.102 0.0
## Antiplatelet use = 1 (%) 0.012 0.0
## Anticoagulant use = 1 (%) 0.054 0.0
## NSAIDs use = 1 (%) 0.012 0.0
## Corticosteroid use = 1 (%) 0.085 0.0
## Acid blocker use = 1 (%) 0.154 0.0
## Vasopressor = 1 (%) 0.022 0.0
## RBC transfusion (mean (SD)) 0.093 0.0
## Total bilirubin (mean (SD)) 0.129 3.2
## Aspartate aminotransferase (mean (SD)) 0.086 2.2
## Alanine aminotransferase (mean (SD)) 0.082 2.2
## Albumin (mean (SD)) 0.021 4.6
## White blood cell (mean (SD)) 0.154 1.8
## Hemoglobin (mean (SD)) 0.119 1.8
## Platelet (mean (SD)) <0.001 1.8
## eGFR < 30 = 1 (%) 0.034 1.9
## C-reactive protein (mean (SD)) 0.030 4.7
## Prothrombin time (mean (SD)) 0.018 5.8
## Activated partial thromboplastin time (mean (SD)) 0.163 11.4
## Shock index > 1 = 1 (%) 0.100 3.0
## antibio_exposure (mean (SD)) Inf 0.0tableobe作成の詳細
library(dplyr)
library(tidyr)
detach("package:plyr", unload = TRUE)## Warning: 'plyr' 名前空間はアンロード出来ません
## 名前空間 'plyr' は 'pROC', 'reshape2', 'rapportools' によってインポートされているのでアンロード出来ません# Continuous variables for which we want median and IQR
continuous_vars <- c("Child-Pugh score", "Charlson Comorbidity Index", "Total bilirubin", "Aspartate aminotransferase", "Alanine aminotransferase", "Albumin", "White blood cell", "Hemoglobin", "Platelet", "C-reactive protein", "Prothrombin time", "Activated partial thromboplastin time")
# 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")## 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(continuous_vars)
##
## # Now:
## data %>% select(all_of(continuous_vars))
##
## 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.# 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
print(tableone)## # A tibble: 29 × 3
## Variable `0` `1`
## <chr> <chr> <chr>
## 1 Acid blocker use 0% 0%
## 2 Activated partial thromboplastin time (29.6, 26.4 - 33.8) (31.1, 28.4 - 35.4)
## 3 Age 0% 0%
## 4 Alanine aminotransferase (27, 19 - 42) (30.5, 20 - 47)
## 5 Albumin (2.9, 2.5 - 3.3) (2.9, 2.5 - 3.3)
## 6 Alcohol related disease 0% 0%
## 7 Anticoagulant use 0% 0%
## 8 Antiplatelet use 0% 0%
## 9 Aspartate aminotransferase (47, 31 - 83) (54.5, 32.2 - 94.8)
## 10 Barthel index 0% 0%
## # ℹ 19 more rowstableone作成の詳細
# Continuous variables for which we want median and IQR
continuous_vars <- c("Child-Pugh score", "Charlson Comorbidity Index", "Total bilirubin", "Aspartate aminotransferase", "Alanine aminotransferase", "Albumin", "White blood cell", "Hemoglobin", "Platelet", "C-reactive protein", "Prothrombin time", "Activated partial thromboplastin time", "RBC transfusion")
# 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)) %>%
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))) %>%
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
print(tableone)## # A tibble: 30 × 3
## Variable `0` `1`
## <chr> <chr> <chr>
## 1 Acid blocker use 0% 0%
## 2 Activated partial thromboplastin time (29.6, 26.4 - 33.8) (31.1, 28.4 - 35.4)
## 3 Age 0% 0%
## 4 Alanine aminotransferase (27, 19 - 42) (30.5, 20 - 47)
## 5 Albumin (2.9, 2.5 - 3.3) (2.9, 2.5 - 3.3)
## 6 Alcohol related disease 0% 0%
## 7 Anticoagulant use 0% 0%
## 8 Antiplatelet use 0% 0%
## 9 Aspartate aminotransferase (47, 31 - 83) (54.5, 32.2 - 94.8)
## 10 Barthel index 0% 0%
## # ℹ 20 more rows