analysis
市田 親正
2023-04-13
必要なlibraryの読み込み
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.1 ✔ 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(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_analysis_0412.xlsx")データの確認
str(data)## tibble [791 × 48] (S3: tbl_df/tbl/data.frame)
## $ hosp_id : num [1:791] 1001 1001 1001 1001 1001 ...
## $ pt_id : num [1:791] 1 2 4 5 6 9 12 14 16 17 ...
## $ ad_num : num [1:791] 1 1 1 1 1 1 1 1 1 1 ...
## $ year : num [1:791] 2012 2011 2011 2010 2017 ...
## $ antibio_exposure : num [1:791] 0 0 0 0 0 0 0 0 1 0 ...
## $ antibio_type_a : chr [1:791] "99" "99" "99" "99" ...
## $ antibio_type_b : chr [1:791] "99" "99" "99" "99" ...
## $ antibio_term : num [1:791] 0 0 0 0 0 0 0 0 5 0 ...
## $ age_num : num [1:791] 50 80 44 67 47 73 65 41 38 62 ...
## $ sex : chr [1:791] "M" "M" "M" "F" ...
## $ bmi_num : num [1:791] NA 25.3 14.5 NA 21 ...
## $ smoking : num [1:791] 0 0 240 0 270 1000 0 210 0 0 ...
## $ brthel_cate : num [1:791] NA 2 1 NA NA 0 1 1 1 2 ...
## $ child_numl : num [1:791] 11 6 8 NA 11 9 6 8 NA 9 ...
## $ child_score : chr [1:791] "c" "a" "b" NA ...
## $ cci_num : num [1:791] 4 4 4 4 4 4 5 5 4 4 ...
## $ malignancy : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ hd : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ hcc : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ alocohol : num [1:791] 1 0 0 0 1 1 0 1 1 0 ...
## $ past_varix_rup : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ antiplate_user_01: num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ anticoag_user01 : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ nsaids_user01 : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ steroid_user01 : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ ppi_user01 : num [1:791] 1 1 1 1 1 1 1 1 1 1 ...
## $ map_day0 : num [1:791] 0 2 6 2 0 4 0 6 0 8 ...
## $ shock : num [1:791] 1 0 1 1 0 1 0 1 0 NA ...
## $ t_bil : num [1:791] 2.2 1.2 3.4 1.2 7.5 2.2 0.4 1.5 1.8 0.5 ...
## $ ast : num [1:791] 217 31 129 52 112 55 55 56 169 56 ...
## $ alt : num [1:791] 63 22 46 36 53 20 67 12 36 30 ...
## $ wbc : num [1:791] 7400 5000 9100 3900 9800 11100 2900 8700 9300 14600 ...
## $ hb : num [1:791] 6.9 10.8 10.7 6.3 12.2 5.6 5.7 4.9 10.9 3.7 ...
## $ plt : num [1:791] 115 77 162 63 69 124 61 129 90 437 ...
## $ alb : num [1:791] 2.2 3.2 2.9 2.8 2.6 2.3 4 2.2 3.2 2.1 ...
## $ eGFR : num [1:791] 58 57.7 112.4 123.6 89.3 ...
## $ eGFR30 : num [1:791] 0 0 0 0 0 0 1 0 0 0 ...
## $ crp : num [1:791] 0.92 0.29 0.29 0.29 0.1 NA 0.04 0.44 0.09 0.17 ...
## $ pt : num [1:791] 37.8 55 37.8 74.6 28.5 45.9 71.5 43.6 66.7 51.3 ...
## $ aptt : num [1:791] 29.4 29 35.3 30.1 37.8 27.5 25 24.8 20.8 30.9 ...
## $ death_day2 : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ ivr_day2 : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ ope_day2 : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ ivr_ope_day2 : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ rebleed_end_hemo : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ sepsis : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ outcome : num [1:791] 0 0 0 0 0 0 0 0 0 0 ...
## $ los : num [1:791] 12 7 10 3 6 8 5 8 5 2 ...データの整理
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),
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")),
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),
map_day0= as.integer(map_day0),
shock=factor(shock),
eGFR30=factor(eGFR30),
death_day2=factor(death_day2),
ivr_day2=factor(ivr_day2),
ope_day2=factor(ope_day2),
ivr_ope_day2=factor(ivr_ope_day2),
rebleed_end_hemo=factor(rebleed_end_hemo),
sepsis=factor(sepsis),
outcome=factor(outcome),
los=as.integer(los)
)整理後の確認
str(df)## tibble [791 × 48] (S3: tbl_df/tbl/data.frame)
## $ hosp_id : int [1:791] 1001 1001 1001 1001 1001 1001 1001 1001 1001 1001 ...
## $ pt_id : int [1:791] 1 2 4 5 6 9 12 14 16 17 ...
## $ ad_num : int [1:791] 1 1 1 1 1 1 1 1 1 1 ...
## $ year : int [1:791] 2012 2011 2011 2010 2017 2010 2013 2014 2016 2019 ...
## $ antibio_exposure : int [1:791] 0 0 0 0 0 0 0 0 1 0 ...
## $ antibio_type_a : Factor w/ 10 levels "a","b","c","d",..: 10 10 10 10 10 10 10 10 4 10 ...
## $ antibio_type_b : Factor w/ 6 levels "c","d","f","h",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ antibio_term : int [1:791] 0 0 0 0 0 0 0 0 5 0 ...
## $ age_num : int [1:791] 50 80 44 67 47 73 65 41 38 62 ...
## $ sex : Factor w/ 2 levels "M","F": 1 1 1 2 1 1 1 2 1 2 ...
## $ bmi_num : num [1:791] NA 25.3 14.5 NA 21 ...
## $ smoking : int [1:791] 0 0 240 0 270 1000 0 210 0 0 ...
## $ brthel_cate : Factor w/ 3 levels "0","1","2": NA 3 2 NA NA 1 2 2 2 3 ...
## $ child_numl : int [1:791] 11 6 8 NA 11 9 6 8 NA 9 ...
## $ child_score : Factor w/ 3 levels "a","b","c": 3 1 2 NA 3 2 1 2 NA 2 ...
## $ cci_num : int [1:791] 4 4 4 4 4 4 5 5 4 4 ...
## $ malignancy : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ hd : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ hcc : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ alocohol : Factor w/ 2 levels "0","1": 2 1 1 1 2 2 1 2 2 1 ...
## $ past_varix_rup : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ antiplate_user_01: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ anticoag_user01 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ nsaids_user01 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ steroid_user01 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ ppi_user01 : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
## $ map_day0 : int [1:791] 0 2 6 2 0 4 0 6 0 8 ...
## $ shock : Factor w/ 2 levels "0","1": 2 1 2 2 1 2 1 2 1 NA ...
## $ t_bil : num [1:791] 2.2 1.2 3.4 1.2 7.5 2.2 0.4 1.5 1.8 0.5 ...
## $ ast : num [1:791] 217 31 129 52 112 55 55 56 169 56 ...
## $ alt : num [1:791] 63 22 46 36 53 20 67 12 36 30 ...
## $ wbc : num [1:791] 7400 5000 9100 3900 9800 11100 2900 8700 9300 14600 ...
## $ hb : num [1:791] 6.9 10.8 10.7 6.3 12.2 5.6 5.7 4.9 10.9 3.7 ...
## $ plt : num [1:791] 115 77 162 63 69 124 61 129 90 437 ...
## $ alb : num [1:791] 2.2 3.2 2.9 2.8 2.6 2.3 4 2.2 3.2 2.1 ...
## $ eGFR : num [1:791] 58 57.7 112.4 123.6 89.3 ...
## $ eGFR30 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 2 1 1 1 ...
## $ crp : num [1:791] 0.92 0.29 0.29 0.29 0.1 NA 0.04 0.44 0.09 0.17 ...
## $ pt : num [1:791] 37.8 55 37.8 74.6 28.5 45.9 71.5 43.6 66.7 51.3 ...
## $ aptt : num [1:791] 29.4 29 35.3 30.1 37.8 27.5 25 24.8 20.8 30.9 ...
## $ death_day2 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ ivr_day2 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ ope_day2 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ ivr_ope_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 ...
## $ sepsis : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ outcome : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ los : int [1:791] 12 7 10 3 6 8 5 8 5 2 ...確認その2
dfSummary(df)## Data Frame Summary
## df
## Dimensions: 791 x 48
## Duplicates: 0
##
## ------------------------------------------------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ------------------- ----------------------------- --------------------- ---------------------- ---------- ---------
## 1 hosp_id Mean (sd) : 1021.6 (21.9) 44 distinct values : 791 0
## [integer] min < med < max: : (100.0%) (0.0%)
## 1001 < 1017 < 1075 :
## IQR (CV) : 20 (0) : : :
## : : : : : :
##
## 2 pt_id Mean (sd) : 416.6 (243.1) 680 distinct values : . 791 0
## [integer] min < med < max: : : : . : : : : (100.0%) (0.0%)
## 1 < 425 < 839 : : : : : : : :
## IQR (CV) : 431 (0.6) : : : : : : : : .
## : : : : : : : : :
##
## 3 ad_num Mean (sd) : 1.2 (0.6) 1 : 673 (85.1%) IIIIIIIIIIIIIIIII 791 0
## [integer] min < med < max: 2 : 79 (10.0%) I (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 : . : 791 0
## [integer] min < med < max: : : : (100.0%) (0.0%)
## 2010 < 2016 < 2022 : : . : . :
## IQR (CV) : 7 (0) : : : : : : : : : :
## : : : : : : : : : :
##
## 5 antibio_exposure Min : 0 0 : 558 (70.5%) IIIIIIIIIIIIII 791 0
## [integer] Mean : 0.3 1 : 233 (29.5%) IIIII (100.0%) (0.0%)
## Max : 1
##
## 6 antibio_type_a 1. a 4 ( 0.5%) 791 0
## [factor] 2. b 32 ( 4.0%) (100.0%) (0.0%)
## 3. c 52 ( 6.6%) I
## 4. d 104 (13.1%) 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.5%) IIIIIIIIIIIIII
##
## 7 antibio_type_b 1. c 1 ( 0.1%) 791 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 777 (98.2%) IIIIIIIIIIIIIIIIIII
##
## 8 antibio_term Mean (sd) : 1.4 (3.3) 23 distinct values : 791 0
## [integer] min < med < max: : (100.0%) (0.0%)
## 0 < 0 < 35 :
## IQR (CV) : 1 (2.4) :
## : .
##
## 9 age_num Mean (sd) : 61.1 (13) 65 distinct values : . 791 0
## [integer] min < med < max: . : : : (100.0%) (0.0%)
## 24 < 62 < 93 : : : : .
## IQR (CV) : 19 (0.2) : : : : : : :
## . : : : : : : : .
##
## 10 sex 1. M 599 (75.7%) IIIIIIIIIIIIIII 791 0
## [factor] 2. F 192 (24.3%) IIII (100.0%) (0.0%)
##
## 11 bmi_num Mean (sd) : 23.4 (4) 553 distinct values : 705 86
## [numeric] min < med < max: : (89.1%) (10.9%)
## 13.7 < 22.8 < 44.7 :
## IQR (CV) : 5 (0.2) : : :
## : : : .
##
## 12 smoking Mean (sd) : 257.5 (390) 116 distinct values : 694 97
## [integer] min < med < max: : (87.7%) (12.3%)
## 0 < 0 < 2200 :
## IQR (CV) : 410 (1.5) :
## : : . .
##
## 13 brthel_cate 1. 0 270 (41.3%) IIIIIIII 654 137
## [factor] 2. 1 215 (32.9%) IIIIII (82.7%) (17.3%)
## 3. 2 169 (25.8%) IIIII
##
## 14 child_numl Mean (sd) : 8.3 (2) 11 distinct values : . : 689 102
## [integer] min < med < max: : : : : (87.1%) (12.9%)
## 5 < 8 < 15 : : : : :
## IQR (CV) : 3 (0.2) : : : : : . .
## : : : : : : : .
##
## 15 child_score 1. a 135 (19.0%) III 711 80
## [factor] 2. b 376 (52.9%) IIIIIIIIII (89.9%) (10.1%)
## 3. c 200 (28.1%) IIIII
##
## 16 cci_num Mean (sd) : 4.5 (1.3) 11 distinct values : 791 0
## [integer] min < med < max: : (100.0%) (0.0%)
## 3 < 4 < 13 :
## IQR (CV) : 1 (0.3) : :
## : : : .
##
## 17 malignancy 1. 0 697 (88.1%) IIIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 94 (11.9%) II (100.0%) (0.0%)
##
## 18 hd 1. 0 780 (98.6%) IIIIIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 11 ( 1.4%) (100.0%) (0.0%)
##
## 19 hcc 1. 0 641 (81.0%) IIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 150 (19.0%) III (100.0%) (0.0%)
##
## 20 alocohol 1. 0 418 (52.8%) IIIIIIIIII 791 0
## [factor] 2. 1 373 (47.2%) IIIIIIIII (100.0%) (0.0%)
##
## 21 past_varix_rup 1. 0 601 (76.0%) IIIIIIIIIIIIIII 791 0
## [factor] 2. 1 190 (24.0%) IIII (100.0%) (0.0%)
##
## 22 antiplate_user_01 1. 0 780 (98.6%) IIIIIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 11 ( 1.4%) (100.0%) (0.0%)
##
## 23 anticoag_user01 1. 0 778 (98.4%) IIIIIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 13 ( 1.6%) (100.0%) (0.0%)
##
## 24 nsaids_user01 1. 0 773 (97.7%) IIIIIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 18 ( 2.3%) (100.0%) (0.0%)
##
## 25 steroid_user01 1. 0 789 (99.7%) IIIIIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 2 ( 0.3%) (100.0%) (0.0%)
##
## 26 ppi_user01 1. 0 91 (11.5%) II 791 0
## [factor] 2. 1 700 (88.5%) IIIIIIIIIIIIIIIII (100.0%) (0.0%)
##
## 27 map_day0 Mean (sd) : 2.9 (2.9) 13 distinct values : 791 0
## [integer] min < med < max: : (100.0%) (0.0%)
## 0 < 4 < 26 : :
## IQR (CV) : 4 (1) : : .
## : : : .
##
## 28 shock 1. 0 476 (62.1%) IIIIIIIIIIII 767 24
## [factor] 2. 1 291 (37.9%) IIIIIII (97.0%) (3.0%)
##
## 29 t_bil Mean (sd) : 2.1 (1.9) 254 distinct values : 766 25
## [numeric] min < med < max: : (96.8%) (3.2%)
## 0.2 < 1.4 < 13.7 :
## IQR (CV) : 1.6 (0.9) : :
## : : : .
##
## 30 ast Mean (sd) : 77.3 (90.1) 195 distinct values : 774 17
## [numeric] min < med < max: : (97.9%) (2.1%)
## 13 < 49 < 984 :
## IQR (CV) : 56 (1.2) :
## : .
##
## 31 alt Mean (sd) : 39.9 (44.5) 124 distinct values : 774 17
## [numeric] min < med < max: : (97.9%) (2.1%)
## 7 < 28 < 562 :
## IQR (CV) : 25 (1.1) :
## : .
##
## 32 wbc Mean (sd) : 8413.7 (4638.1) 345 distinct values : 777 14
## [numeric] min < med < max: : : (98.2%) (1.8%)
## 1400 < 7400 < 58400 : :
## IQR (CV) : 5020 (0.6) : :
## : : :
##
## 33 hb Mean (sd) : 8.8 (2.5) 120 distinct values . . : 777 14
## [numeric] min < med < max: : : : (98.2%) (1.8%)
## 2.1 < 8.6 < 16.5 : : : :
## IQR (CV) : 3.3 (0.3) : : : : : .
## . : : : : : : . .
##
## 34 plt Mean (sd) : 116.4 (69.3) 217 distinct values : 777 14
## [numeric] min < med < max: : (98.2%) (1.8%)
## 18 < 102 < 1073 : .
## IQR (CV) : 68 (0.6) : :
## : : .
##
## 35 alb Mean (sd) : 2.9 (0.6) 35 distinct values : 755 36
## [numeric] min < med < max: : : (95.4%) (4.6%)
## 1.1 < 2.9 < 4.9 : : :
## IQR (CV) : 0.8 (0.2) : : : :
## : : : : : .
##
## 36 eGFR Mean (sd) : 71.6 (30.7) 737 distinct values : . 776 15
## [numeric] min < med < max: : : . (98.1%) (1.9%)
## 4.3 < 69.1 < 196 : : :
## IQR (CV) : 40.7 (0.4) : : : : .
## : : : : : : . .
##
## 37 eGFR30 1. 0 717 (92.3%) IIIIIIIIIIIIIIIIII 777 14
## [factor] 2. 1 60 ( 7.7%) I (98.2%) (1.8%)
##
## 38 crp Mean (sd) : 0.7 (1.4) 360 distinct values : 754 37
## [numeric] min < med < max: : (95.3%) (4.7%)
## 0 < 0.3 < 18.6 :
## IQR (CV) : 0.7 (1.9) :
## : .
##
## 39 pt Mean (sd) : 57.1 (16.7) 384 distinct values . . : 745 46
## [numeric] min < med < max: : : : (94.2%) (5.8%)
## 13.4 < 58 < 107.9 : : : : :
## IQR (CV) : 23.8 (0.3) . : : : : : .
## . : : : : : : : . .
##
## 40 aptt Mean (sd) : 32.7 (14.3) 234 distinct values : 701 90
## [numeric] min < med < max: : (88.6%) (11.4%)
## 17.6 < 30.2 < 200 :
## IQR (CV) : 7.6 (0.4) :
## : :
##
## 41 death_day2 1. 0 743 (93.9%) IIIIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 48 ( 6.1%) I (100.0%) (0.0%)
##
## 42 ivr_day2 1. 0 787 (99.5%) IIIIIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 4 ( 0.5%) (100.0%) (0.0%)
##
## 43 ope_day2 1. 0 787 (99.5%) IIIIIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 4 ( 0.5%) (100.0%) (0.0%)
##
## 44 ivr_ope_day2 1. 0 783 (99.0%) IIIIIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 8 ( 1.0%) (100.0%) (0.0%)
##
## 45 rebleed_end_hemo 1. 0 680 (86.0%) IIIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 111 (14.0%) II (100.0%) (0.0%)
##
## 46 sepsis 1. 0 782 (98.9%) IIIIIIIIIIIIIIIIIII 791 0
## [factor] 2. 1 9 ( 1.1%) (100.0%) (0.0%)
##
## 47 outcome 1. 0 629 (79.5%) IIIIIIIIIIIIIII 791 0
## [factor] 2. 1 162 (20.5%) IIII (100.0%) (0.0%)
##
## 48 los Mean (sd) : 13 (15) 58 distinct values : 791 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//RtmpzS2jMH/file69f39ddaee.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.
素データでのtable1の作成
col_continuous <-
c("age_num","child_numl","cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "alb", "crp", "pt", "aptt")
col_factors <-
c("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")
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 233
## age_num (mean (SD)) 61.54 (12.63) 60.05 (13.75) 0.141
## 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.967
## map_day0 (mean (SD)) 2.86 (2.92) 3.12 (2.98) 0.254
## t_bil (mean (SD)) 2.01 (1.88) 2.26 (2.02) 0.107
## ast (mean (SD)) 74.97 (82.14) 82.95 (106.36) 0.259
## alt (mean (SD)) 38.83 (42.49) 42.48 (48.91) 0.295
## wbc (mean (SD)) 8194.22 (4242.15) 8929.18 (5431.81) 0.043
## hb (mean (SD)) 8.69 (2.47) 9.00 (2.40) 0.111
## plt (mean (SD)) 116.49 (61.73) 116.29 (84.69) 0.971
## alb (mean (SD)) 2.91 (0.57) 2.93 (0.64) 0.761
## crp (mean (SD)) 0.76 (1.39) 0.72 (1.57) 0.679
## pt (mean (SD)) 57.00 (16.56) 57.44 (17.09) 0.740
## aptt (mean (SD)) 31.95 (12.69) 34.40 (17.30) 0.036
## sex = F (%) 141 (25.3) 51 (21.9) 0.358
## brthel_cate (%) 0.494
## 0 186 (41.5) 84 (40.8)
## 1 152 (33.9) 63 (30.6)
## 2 110 (24.6) 59 (28.6)
## malignancy = 1 (%) 65 (11.6) 29 (12.4) 0.845
## 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.3) 0.258
## alocohol = 1 (%) 246 (44.1) 127 (54.5) 0.009
## past_varix_rup = 1 (%) 127 (22.8) 63 (27.0) 0.233
## antiplate_user_01 = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## anticoag_user01 = 1 (%) 8 ( 1.4) 5 ( 2.1) 0.681
## nsaids_user01 = 1 (%) 13 ( 2.3) 5 ( 2.1) 1.000
## steroid_user01 = 1 (%) 2 ( 0.4) 0 ( 0.0) 0.890
## ppi_user01 = 1 (%) 486 (87.1) 214 (91.8) 0.074
## shock = 1 (%) 197 (36.5) 94 (41.2) 0.255
## eGFR30 = 1 (%) 40 ( 7.3) 20 ( 8.6) 0.642
## Stratified by antibio_exposure
## SMD Missing
## n
## age_num (mean (SD)) 0.113 0.0
## child_numl (mean (SD)) 0.054 12.9
## cci_num (mean (SD)) 0.003 0.0
## map_day0 (mean (SD)) 0.089 0.0
## t_bil (mean (SD)) 0.126 3.2
## ast (mean (SD)) 0.084 2.1
## alt (mean (SD)) 0.080 2.1
## wbc (mean (SD)) 0.151 1.8
## hb (mean (SD)) 0.126 1.8
## plt (mean (SD)) 0.003 1.8
## alb (mean (SD)) 0.024 4.6
## crp (mean (SD)) 0.032 4.7
## pt (mean (SD)) 0.026 5.8
## aptt (mean (SD)) 0.162 11.4
## sex = F (%) 0.080 0.0
## brthel_cate (%) 0.100 17.3
## 0
## 1
## 2
## malignancy = 1 (%) 0.025 0.0
## child_score (%) 0.058 10.1
## a
## b
## c
## hd = 1 (%) 0.013 0.0
## hcc = 1 (%) 0.098 0.0
## alocohol = 1 (%) 0.210 0.0
## past_varix_rup = 1 (%) 0.099 0.0
## antiplate_user_01 = 1 (%) 0.013 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.155 0.0
## shock = 1 (%) 0.096 3.0
## eGFR30 = 1 (%) 0.047 1.8連続量の可視化
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 412 rows containing non-finite values (`stat_bin()`).
単変量の解析をみてみる
# 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"))
sepsis_model <- glm(sepsis ~ 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()
sepsis_or <- coef(sepsis_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...sepsis_ci <- confint(sepsis_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.3491 95% CI: 0.9302 to 1.9421"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.9852 95% CI: 0.5029 to 1.8344"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.6308 95% CI: 1.0691 to 2.4663"print("Sepsis:")## [1] "Sepsis:"print(paste("OR:", round(sepsis_or["antibio_exposure"], 4), "95% CI:", round(sepsis_ci$conf.low[2], 4), "to", round(sepsis_ci$conf.high[2], 4)))## [1] "OR: 0.2963 95% CI: 0.0159 to 1.628"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.0132 95% CI: 0.9023 to 1.1393"多重代入
imp1 <- mice(df,
m=5, #mは代入の結果として作成するデータセットの数
maxit = 5, #maxitは反復回数
method="pmm",
seed = 10
)##
## iter imp variable
## 1 1 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 1 2 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 1 3 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 1 4 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 1 5 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 2 1 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt
## 2 2 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 2 3 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp* pt* aptt*
## 2 4 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp* pt* aptt*
## 2 5 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 3 1 bmi_num smoking brthel_cate child_numl* child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 3 2 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp* pt* aptt
## 3 3 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 3 4 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 3 5 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp* pt* aptt*
## 4 1 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp* pt* aptt*
## 4 2 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp* pt* aptt
## 4 3 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 4 4 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 4 5 bmi_num smoking brthel_cate* child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 5 1 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp* pt* aptt*
## 5 2 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 5 3 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*
## 5 4 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt
## 5 5 bmi_num smoking brthel_cate child_numl child_score* shock t_bil ast alt wbc* hb* plt* alb eGFR* eGFR30* crp pt* aptt*## Warning: Number of logged events: 655主要評価解析
# 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:5) {
imputed_data <- complete(imp1, i)
# Fit GEE propensity score model
propensity_model <- geeglm(antibio_exposure ~ age_num + 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 + 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:5) {
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", "sepsis")
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!## 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!
## 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!
## 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)
)
# Reorder forest_plot_data by the desired order of outcomes
forest_plot_data$Outcome <- factor(forest_plot_data$Outcome, levels = c("outcome", "death_day2", "rebleed_end_hemo", "sepsis"))
# Display results as a table
print(forest_plot_data)## Outcome OR Lower_CI Upper_CI P_value
## 1 outcome 1.3083878 1.0361863 1.6520953 0.023901862
## 2 death_day2 0.8770159 0.5869584 1.3104115 0.521840146
## 3 rebleed_end_hemo 1.6330407 1.2463294 2.1397409 0.000374995
## 4 sepsis 0.3071999 0.1008344 0.9359087 0.037844259# 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("sepsis", "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## 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.5984
##
## [[2]]
## Area under the curve: 0.6128
##
## [[3]]
## Area under the curve: 0.6155
##
## [[4]]
## Area under the curve: 0.6411
##
## [[5]]
## Area under the curve: 0.6179# Plot the propensity score overlap for each imputed dataset
for (i in 1:5) {
ggplot(results[[1]]$PS_data[[i]], aes(x = ps, fill = as.factor(antibio_exposure))) +
geom_histogram(position = "identity", alpha = 0.5, bins = 30) +
labs(
title = paste0("Imputed dataset ", i),
x = "Propensity Score",
y = "Count",
fill = "Antibio Exposure"
) +
theme_bw() +
scale_fill_discrete(labels = c("No Exposure", "Exposure")) +
xlim(0, 1) +
guides(fill = guide_legend(reverse = TRUE)) -> p
print(p)
}## 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()`).
## 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" = "Complex outcome",
"death_day2" = "Death in hospitalization",
"rebleed_end_hemo" = "Rebleeding",
"sepsis" = "Bacteraemia"))
# Reorder forest_plot_data by the desired order of outcomes
forest_plot_data$Outcome <- factor(forest_plot_data$Outcome,
levels = c("Complex outcome",
"Death in hospitalization",
"Rebleeding",
"Bacteraemia"))
# Display results as a table
print(forest_plot_data)## Outcome OR Lower_CI Upper_CI P_value
## 1 Complex outcome 1.3083878 1.0361863 1.6520953 0.023901862
## 2 Death in hospitalization 0.8770159 0.5869584 1.3104115 0.521840146
## 3 Rebleeding 1.6330407 1.2463294 2.1397409 0.000374995
## 4 Bacteraemia 0.3071999 0.1008344 0.9359087 0.037844259# 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("Bacteraemia", "Rebleeding", "Death in hospitalization", "Complex outcome"))
# Display forest plot
forest_plot_ordered
副次評価 解析
#LOS
# 本解析: IPTW models on imputed data with Gamma distribution
gee_results <- list()
summary_tables <- list()
conf_ints <- list()
weighted_data <- list()
for (i in 1:5) {
imputed_data <- complete(imp1, i)
# Fit GEE propensity score model
propensity_model <- geeglm(antibio_exposure ~ age_num + 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 + 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 gamma regression model
stabilized_weights <- with(imputed_data, antibio_exposure/(ps) + (1 - antibio_exposure)/(1 - ps))
iptw_model <- glm(los ~ antibio_exposure, data = imputed_data, weights = stabilized_weights, family = Gamma("log"))
# Store the results
gee_results[[i]] <- iptw_model
summary_tables[[i]] <- summary(iptw_model)
conf_ints[[i]] <- confint(iptw_model, level = 0.95, method = "Wald")
# Store weighted data
imputed_data$weight <- stabilized_weights
weighted_data[[i]] <- imputed_data
}## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...# Save the results from each imputed dataset
params <- list()
for (i in 1:5) {
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 Gamma distribution)\n")##
## Pooled results (IPTW with Gamma distribution)cat("\nRate Ratio:", pooled_rr, "\n")##
## Rate Ratio: 1.020283cat("\n95% Confidence Interval: (", pooled_lower, ", ", pooled_upper, ")\n")##
## 95% Confidence Interval: ( 0.8485043 , 1.226837 )cat("\nP-value:", pooled_p, "\n")##
## P-value: 0.8309559負の二項分布で解析してみる
# 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:5) {
imputed_data <- complete(imp1, i)
# Fit GEE propensity score model
propensity_model <- geeglm(antibio_exposure ~ age_num + 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 + 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...
## 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:5) {
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.020283cat("\n95% Confidence Interval: (", pooled_lower, ", ", pooled_upper, ")\n")##
## 95% Confidence Interval: ( 0.9456284 , 1.100831 )cat("\nP-value:", pooled_p, "\n")##
## P-value: 0.6044956サマリー表の作成 補完後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 233
## age_num (mean (SD)) 61.54 (12.63) 60.05 (13.75) 0.141
## child_numl (mean (SD)) 8.12 (2.02) 8.29 (2.09) 0.281
## cci_num (mean (SD)) 4.51 (1.32) 4.52 (1.33) 0.967
## map_day0 (mean (SD)) 2.86 (2.92) 3.12 (2.98) 0.254
## t_bil (mean (SD)) 1.99 (1.86) 2.23 (2.01) 0.107
## ast (mean (SD)) 76.27 (83.60) 82.97 (106.04) 0.345
## alt (mean (SD)) 39.28 (43.00) 42.31 (48.75) 0.386
## wbc (mean (SD)) 8120.04 (4229.65) 8917.08 (5423.23) 0.027
## hb (mean (SD)) 8.82 (2.59) 9.02 (2.43) 0.302
## plt (mean (SD)) 114.65 (62.22) 116.10 (84.56) 0.789
## alb (mean (SD)) 2.93 (0.58) 2.93 (0.64) 0.976
## crp (mean (SD)) 1.19 (2.96) 0.78 (1.85) 0.050
## pt (mean (SD)) 59.23 (18.16) 58.52 (17.64) 0.615
## aptt (mean (SD)) 45.16 (43.06) 44.53 (40.28) 0.848
## sex = F (%) 141 (25.3) 51 (21.9) 0.358
## brthel_cate (%) 0.372
## 0 232 (41.6) 97 (41.6)
## 1 191 (34.2) 70 (30.0)
## 2 135 (24.2) 66 (28.3)
## malignancy = 1 (%) 65 (11.6) 29 (12.4) 0.845
## child_score (%) 0.895
## a 93 (16.7) 42 (18.0)
## b 266 (47.7) 110 (47.2)
## c 199 (35.7) 81 (34.8)
## hd = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## hcc = 1 (%) 112 (20.1) 38 (16.3) 0.258
## alocohol = 1 (%) 246 (44.1) 127 (54.5) 0.009
## past_varix_rup = 1 (%) 127 (22.8) 63 (27.0) 0.233
## antiplate_user_01 = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## anticoag_user01 = 1 (%) 8 ( 1.4) 5 ( 2.1) 0.681
## nsaids_user01 = 1 (%) 13 ( 2.3) 5 ( 2.1) 1.000
## steroid_user01 = 1 (%) 2 ( 0.4) 0 ( 0.0) 0.890
## ppi_user01 = 1 (%) 486 (87.1) 214 (91.8) 0.074
## shock = 1 (%) 204 (36.6) 95 (40.8) 0.301
## eGFR30 = 1 (%) 40 ( 7.2) 20 ( 8.6) 0.591
## Stratified by antibio_exposure
## SMD Missing
## n
## age_num (mean (SD)) 0.113 0.0
## child_numl (mean (SD)) 0.084 0.0
## cci_num (mean (SD)) 0.003 0.0
## map_day0 (mean (SD)) 0.089 0.0
## t_bil (mean (SD)) 0.124 0.0
## ast (mean (SD)) 0.070 0.0
## alt (mean (SD)) 0.066 0.0
## wbc (mean (SD)) 0.164 0.0
## hb (mean (SD)) 0.082 0.0
## plt (mean (SD)) 0.020 0.0
## alb (mean (SD)) 0.002 0.0
## crp (mean (SD)) 0.166 0.0
## pt (mean (SD)) 0.040 0.0
## aptt (mean (SD)) 0.015 0.0
## sex = F (%) 0.080 0.0
## brthel_cate (%) 0.110 0.0
## 0
## 1
## 2
## malignancy = 1 (%) 0.025 0.0
## child_score (%) 0.037 0.0
## a
## b
## c
## hd = 1 (%) 0.013 0.0
## hcc = 1 (%) 0.098 0.0
## alocohol = 1 (%) 0.210 0.0
## past_varix_rup = 1 (%) 0.099 0.0
## antiplate_user_01 = 1 (%) 0.013 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.155 0.0
## shock = 1 (%) 0.087 0.0
## eGFR30 = 1 (%) 0.053 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
## 1 age_num -0.1145
## 26 steroid_user01 0.0740
## 5 t_bil 0.0720
## 2 child_numl 0.0647
## 14 aptt -0.0533
## 12 crp -0.0329
## 4 map_day0 0.0324
## 10 plt -0.0285
## 13 pt -0.0282
## 3 cci_num 0.0276
## 9 hb 0.0275
## 8 wbc 0.0236
## 16 brthel_cate -0.0056
## 11 alb 0.0049
## 18 child_score 0.0046
## 6 ast -0.0039
## 7 alt 0.0039
## 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 shock 0.0000
## 29 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))
)
# Create Love plot
ggplot(smd_df_pre_post, 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))
# Covariate name mapping
covariate_name_mapping <- c(
"age_num" = "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",
"shock" = "Shock index > 1",
"eGFR30" = "eGFR < 30"
)
# Update Covariate names in smd_df_pre_post
smd_df_pre_post$Covariate <- covariate_name_mapping[as.character(smd_df_pre_post$Covariate)]
smd_df_pre_post$Covariate <- factor(smd_df_pre_post$Covariate, levels = covariate_name_mapping[sort(c(col_continuous, col_factors))])
# Create Love plot
ggplot(smd_df_pre_post, 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))
素データでの背景サマリーの作成
# Rename columns
col_continuous <- c("Age", "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("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", "Shock index > 1", "eGFR < 30")
# Update column names in the dataframe
colnames(df)[colnames(df) %in% c("age_num", "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("sex", "brthel_cate", "child_score", "hd", "hcc", "malignancy", "alocohol", "past_varix_rup", "antiplate_user_01", "anticoag_user01", "nsaids_user01", "steroid_user01", "ppi_user01", "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", "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 (mean (SD)) 61.54 (12.63)
## 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)
## 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 233
## Age (mean (SD)) 60.05 (13.75) 0.141
## Sex = F (%) 51 (21.9) 0.358
## Barthel index (%) 0.494
## 0 84 (40.8)
## 1 63 (30.6)
## 2 59 (28.6)
## 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.967
## Maintenance hemodialysis = 1 (%) 29 (12.4) 0.845
## Hepatic cancer = 1 (%) 3 ( 1.3) 1.000
## Malignant tomor history = 1 (%) 38 (16.3) 0.258
## Alcohol related disease = 1 (%) 127 (54.5) 0.009
## Past varix rupture history = 1 (%) 63 (27.0) 0.233
## Antiplatelet use = 1 (%) 3 ( 1.3) 1.000
## Anticoagulant use = 1 (%) 5 ( 2.1) 0.681
## NSAIDs use = 1 (%) 5 ( 2.1) 1.000
## Corticosteroid use = 1 (%) 0 ( 0.0) 0.890
## Acid blocker use = 1 (%) 214 (91.8) 0.074
## RBC transfusion (mean (SD)) 3.12 (2.98) 0.254
## Total bilirubin (mean (SD)) 2.26 (2.02) 0.107
## Aspartate aminotransferase (mean (SD)) 82.95 (106.36) 0.259
## Alanine aminotransferase (mean (SD)) 42.48 (48.91) 0.295
## Albumin (mean (SD)) 2.93 (0.64) 0.761
## White blood cell (mean (SD)) 8929.18 (5431.81) 0.043
## Hemoglobin (mean (SD)) 9.00 (2.40) 0.111
## Platelet (mean (SD)) 116.29 (84.69) 0.971
## eGFR < 30 = 1 (%) 20 ( 8.6) 0.642
## C-reactive protein (mean (SD)) 0.72 (1.57) 0.679
## Prothrombin time (mean (SD)) 57.44 (17.09) 0.740
## Activated partial thromboplastin time (mean (SD)) 34.40 (17.30) 0.036
## Shock index > 1 = 1 (%) 94 (41.2) 0.255
## antibio_exposure (mean (SD)) 1.00 (0.00) <0.001
## Stratified by antibio_exposure
## test SMD Missing
## n
## Age (mean (SD)) 0.113 0.0
## Sex = F (%) 0.080 0.0
## Barthel index (%) 0.100 17.3
## 0
## 1
## 2
## Child-Pugh classification (%) 0.058 10.1
## a
## b
## c
## Child-Pugh score (mean (SD)) 0.054 12.9
## Charlson Comorbidity Index (mean (SD)) 0.003 0.0
## Maintenance hemodialysis = 1 (%) 0.025 0.0
## Hepatic cancer = 1 (%) 0.013 0.0
## Malignant tomor history = 1 (%) 0.098 0.0
## Alcohol related disease = 1 (%) 0.210 0.0
## Past varix rupture history = 1 (%) 0.099 0.0
## Antiplatelet use = 1 (%) 0.013 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.155 0.0
## RBC transfusion (mean (SD)) 0.089 0.0
## Total bilirubin (mean (SD)) 0.126 3.2
## Aspartate aminotransferase (mean (SD)) 0.084 2.1
## Alanine aminotransferase (mean (SD)) 0.080 2.1
## Albumin (mean (SD)) 0.024 4.6
## White blood cell (mean (SD)) 0.151 1.8
## Hemoglobin (mean (SD)) 0.126 1.8
## Platelet (mean (SD)) 0.003 1.8
## eGFR < 30 = 1 (%) 0.047 1.8
## C-reactive protein (mean (SD)) 0.032 4.7
## Prothrombin time (mean (SD)) 0.026 5.8
## Activated partial thromboplastin time (mean (SD)) 0.162 11.4
## Shock index > 1 = 1 (%) 0.096 3.0
## antibio_exposure (mean (SD)) Inf 0.0library(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("Age", "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: 28 × 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.5 - 35.4)
## 3 Age (62, 52 - 70) (60, 49 - 71)
## 4 Alanine aminotransferase (27, 19 - 42) (30, 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, 32.5 - 94.5)
## 10 Barthel index 0% 0%
## # … with 18 more rows# Continuous variables for which we want median and IQR
continuous_vars <- c("Age", "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: 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.5 - 35.4)
## 3 Age (62, 52 - 70) (60, 49 - 71)
## 4 Alanine aminotransferase (27, 19 - 42) (30, 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, 32.5 - 94.5)
## 10 Barthel index 0% 0%
## # … with 19 more rows