s_ATE_GEE→GLM
市田 親正
2023-05-23
必要なlibraryの読み込み
library(tidyverse)
library(readxl)
library(tableone)
library(lubridate)
library(skimr)
library(svglite)
library(haven)
library(summarytools)
library(naniar)
library(devtools)
library(reader)
library(stringr)
library(missForest)
library(mice)
library(geepack)
library(ggplot2)
library(cobalt)
library(tableone)
library(sandwich)
library(survey)
library(WeightIt)
library(MatchIt)
library(twang)
library(naniar)
pacman::p_load(norm2, miceadds, lmtest, car, ROCR, pROC,future,furrr)データ読み込み
#data <- read_excel("varix_icu_ok_3.xlsx")data <- read_excel("varix_icu_ok_3_0926.xlsx")データの整理
df <-
data|> #miceパッケージが順序ありの因子型でないと読めないためこれを指定。
mutate(
hosp_id=as.integer(hosp_id),
year=as.integer(year),
pt_id=as.integer(pt_id),
ad_num=as.integer(ad_num),
antibio_exposure=as.integer(antibio_exposure), #geeglmを回すにはfactorでは回らない
antibio_type_a= factor(antibio_type_a, levels = c("a", "b","c","d","e","f","g","h","i","99")),
antibio_type_b= factor(antibio_type_b, levels = c("c","d","f","h","i","99")),
antibio_term= as.integer(antibio_term ),
age_num=as.integer(age_num),
age_cate=factor(age_cate,levels = c("0", "1", "2", "3")),
sex= factor(sex, levels = c("M", "F")),
smoking= as.integer(smoking),
brthel_cate= factor(brthel_cate, levels = c("0", "1", "2")),
child_numl= as.integer(child_numl),
child_score=factor(child_score, levels = c("a", "b", "c")),
jcs_all=factor(jcs_all, levels = c("0", "1", "2","3","10","20","30","100","200","300")),
cci_num=as.integer(cci_num),
malignancy=factor(malignancy),
hd=factor(hd),
hcc=factor(hcc),
alocohol=factor(alocohol),
past_varix_rup=factor(past_varix_rup),
antiplate_user_01=factor(antiplate_user_01),
anticoag_user01=factor(anticoag_user01),
nsaids_user01=factor(nsaids_user01),
steroid_user01=factor(steroid_user01),
ppi_user01=factor(ppi_user01),
beta_user01=factor(beta_user01),
vaso=factor(vaso),
map_day0= as.integer(map_day0),
shock=factor(shock),
alb_cate= factor(alb_cate, levels = c("0", "1", "2")),
pt_cate= factor(pt_cate, levels = c("0", "1", "2")),
aptt_cate=factor(aptt_cate,levels = c("0", "1", "2")),
eGFR30=factor(eGFR30),
outcome=factor(outcome),
death_day2=factor(death_day2),
rebleed_end_hemo=factor(rebleed_end_hemo),
sbp=factor(sbp),
infection=factor(infection),
infection_all=factor(infection_all),
los=as.integer(los),
cdi=factor(cdi)
)素データでのtable1の作成
col_continuous <-
c("child_numl","cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "crp")
col_factors <-
c("age_cate","sex", "brthel_cate", "malignancy", "child_score", "hd", "hcc", "alocohol", "past_varix_rup", "antiplate_user_01", "anticoag_user01", "nsaids_user01" ,"steroid_user01", "ppi_user01","beta_user01", "vaso", "shock", "alb_cate","pt_cate", "aptt_cate", "eGFR30")
df|> select(c(col_continuous, col_factors, "antibio_exposure"))|>
CreateTableOne(vars = c(col_continuous, col_factors), strata = "antibio_exposure", factorVars = col_factors) -> tableone## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(col_continuous)
##
## # Now:
## data %>% select(all_of(col_continuous))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(col_factors)
##
## # Now:
## data %>% select(all_of(col_factors))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.print(tableone, smd = T, missing = T, test = T, explain = T) ## Stratified by antibio_exposure
## 0 1 p test
## n 558 232
## child_numl (mean (SD)) 8.30 (1.96) 8.41 (2.03) 0.508
## cci_num (mean (SD)) 4.51 (1.32) 4.52 (1.33) 0.917
## map_day0 (mean (SD)) 2.86 (2.92) 3.13 (2.98) 0.231
## t_bil (mean (SD)) 2.01 (1.88) 2.26 (2.03) 0.098
## ast (mean (SD)) 74.97 (82.14) 83.16 (106.55) 0.248
## alt (mean (SD)) 38.83 (42.49) 42.57 (49.00) 0.285
## wbc (mean (SD)) 8194.22 (4242.15) 8945.32 (5438.02) 0.039
## hb (mean (SD)) 8.69 (2.47) 8.98 (2.39) 0.133
## plt (mean (SD)) 116.49 (61.73) 116.46 (84.84) 0.996
## crp (mean (SD)) 0.76 (1.39) 0.72 (1.57) 0.699
## age_cate (%) 0.155
## 0 322 (57.7) 143 (61.6)
## 1 144 (25.8) 46 (19.8)
## 2 75 (13.4) 39 (16.8)
## 3 17 ( 3.0) 4 ( 1.7)
## sex = F (%) 141 (25.3) 51 (22.0) 0.374
## brthel_cate (%) 0.486
## 0 186 (41.5) 83 (40.5)
## 1 152 (33.9) 63 (30.7)
## 2 110 (24.6) 59 (28.8)
## malignancy = 1 (%) 65 (11.6) 29 (12.5) 0.829
## child_score (%) 0.776
## a 93 (18.8) 42 (19.4)
## b 266 (53.7) 110 (50.9)
## c 136 (27.5) 64 (29.6)
## hd = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## hcc = 1 (%) 112 (20.1) 38 (16.4) 0.269
## alocohol = 1 (%) 246 (44.1) 127 (54.7) 0.008
## past_varix_rup = 1 (%) 127 (22.8) 63 (27.2) 0.221
## antiplate_user_01 = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## anticoag_user01 = 1 (%) 8 ( 1.4) 5 ( 2.2) 0.675
## nsaids_user01 = 1 (%) 13 ( 2.3) 5 ( 2.2) 1.000
## steroid_user01 = 1 (%) 2 ( 0.4) 0 ( 0.0) 0.892
## ppi_user01 = 1 (%) 486 (87.1) 213 (91.8) 0.077
## beta_user01 = 1 (%) 26 ( 4.7) 26 (11.2) 0.001
## vaso = 1 (%) 19 ( 3.4) 7 ( 3.0) 0.953
## shock = 1 (%) 197 (36.5) 94 (41.4) 0.236
## alb_cate (%) 0.377
## 0 72 (13.6) 37 (16.5)
## 1 256 (48.3) 97 (43.3)
## 2 202 (38.1) 90 (40.2)
## pt_cate (%) 0.977
## 0 113 (21.7) 48 (21.5)
## 1 324 (62.3) 138 (61.9)
## 2 83 (16.0) 37 (16.6)
## aptt_cate (%) 0.564
## 0 436 (89.3) 185 (87.3)
## 1 47 ( 9.6) 23 (10.8)
## 2 5 ( 1.0) 4 ( 1.9)
## eGFR30 = 1 (%) 40 ( 7.3) 19 ( 8.3) 0.769
## Stratified by antibio_exposure
## SMD Missing
## n
## child_numl (mean (SD)) 0.054 12.8
## cci_num (mean (SD)) 0.008 0.0
## map_day0 (mean (SD)) 0.093 0.0
## t_bil (mean (SD)) 0.129 3.2
## ast (mean (SD)) 0.086 2.2
## alt (mean (SD)) 0.082 2.2
## wbc (mean (SD)) 0.154 1.8
## hb (mean (SD)) 0.119 1.8
## plt (mean (SD)) <0.001 1.8
## crp (mean (SD)) 0.030 4.7
## age_cate (%) 0.183 0.0
## 0
## 1
## 2
## 3
## sex = F (%) 0.077 0.0
## brthel_cate (%) 0.101 17.3
## 0
## 1
## 2
## malignancy = 1 (%) 0.026 0.0
## child_score (%) 0.058 10.0
## a
## b
## c
## hd = 1 (%) 0.012 0.0
## hcc = 1 (%) 0.096 0.0
## alocohol = 1 (%) 0.214 0.0
## past_varix_rup = 1 (%) 0.102 0.0
## antiplate_user_01 = 1 (%) 0.012 0.0
## anticoag_user01 = 1 (%) 0.054 0.0
## nsaids_user01 = 1 (%) 0.012 0.0
## steroid_user01 = 1 (%) 0.085 0.0
## ppi_user01 = 1 (%) 0.154 0.0
## beta_user01 = 1 (%) 0.244 0.0
## vaso = 1 (%) 0.022 0.0
## shock = 1 (%) 0.100 3.0
## alb_cate (%) 0.111 4.6
## 0
## 1
## 2
## pt_cate (%) 0.017 5.9
## 0
## 1
## 2
## aptt_cate (%) 0.084 11.4
## 0
## 1
## 2
## eGFR30 = 1 (%) 0.034 1.9多重代入
#imp1 <- mice(df,
# m=100, #mは代入の結果として作成するデータセットの数
# maxit =200, #maxitは反復回数
# method="pmm",
# printFlag = FALSE,
# seed = 10
# )多重代入のロード
load("impglm.Rdata")主要評価解析
# Function to perform IPTW analysis on different outcomes
perform_iptw_analysis <- function(outcome_name) {
glm_results <- list()
summary_tables <- list()
conf_ints <- list()
weighted_data <- list()
ps_data <- list()
for (i in 1:100) {
imputed_data <- complete(impglm, i)
# Fit GLM propensity score model
propensity_model <- glm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hd + hcc + alocohol + past_varix_rup + antiplate_user_01 + anticoag_user01 + nsaids_user01 + steroid_user01 + ppi_user01 + beta_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt_cate + aptt_cate,
data = imputed_data,
family = binomial("logit"))
# 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
glm_results[[i]] <- iptw_model
summary_tables[[i]] <- summary(iptw_model)
conf_ints[[i]] <- confint(iptw_model, level = 0.95, method = "Wald")
# Store weighted data
imputed_data$weight <- stabilized_weights
weighted_data[[i]] <- imputed_data
ps_data[[i]] <- imputed_data[, c("ps", "antibio_exposure")]
}
# Save the results from each imputed dataset
params <- list()
for (i in 1:100) {
params[[i]] <- data.frame(Estimate = coef(glm_results[[i]])[2],
Std.Error = coef(summary_tables[[i]])[2, "Std. Error"],
df = Inf)
}
# Calculate the pooled estimates
Qbar <- mean(sapply(params, function(x) x$Estimate))
se_Qbar <- sqrt(sum(sapply(params, function(x) x$Std.Error^2))/length(params))
# Calculate the pooled odds ratio and confidence intervals
pooled_or <- exp(Qbar)
pooled_lower <- exp(Qbar - 1.96 * se_Qbar)
pooled_upper <- exp(Qbar + 1.96 * se_Qbar)
# Calculate the pooled p-value (two-tailed)
pooled_p <- 2 * (1 - pnorm(abs(Qbar / se_Qbar)))
return(list(OR = pooled_or, Lower_CI = pooled_lower, Upper_CI = pooled_upper, P_value = pooled_p, PS_data = ps_data))
}
# Perform IPTW analysis for each outcome
outcomes <- c("outcome", "death_day2", "rebleed_end_hemo", "sbp", "cdi")
results <- lapply(outcomes, perform_iptw_analysis)## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...# Combine results into a data frame for plotting
forest_plot_data <- data.frame(
Outcome = outcomes,
OR = sapply(results, function(x) x$OR),
Lower_CI = sapply(results, function(x) x$Lower_CI),
Upper_CI = sapply(results, function(x) x$Upper_CI),
P_value = sapply(results, function(x) x$P_value)
)
# Display results as a table
print(forest_plot_data)## Outcome OR Lower_CI Upper_CI P_value
## 1 outcome 1.2635199 0.9183459 1.738433 0.15078467
## 2 death_day2 0.9991775 0.6744053 1.480350 0.99672647
## 3 rebleed_end_hemo 1.5795179 0.9310448 2.679653 0.09006373
## 4 sbp 1.6304365 0.8337447 3.188414 0.15311315
## 5 cdi 2.8982187 0.3869646 21.706565 0.30028935# Create forest plot with the specified order and adjusted p-value position
forest_plot_ordered <- ggplot(forest_plot_data, aes(x = Outcome, y = OR, ymin = Lower_CI, ymax = Upper_CI)) +
geom_pointrange(size = 0.5) +
geom_hline(yintercept = 1, linetype = "dashed") +
scale_y_log10() +
coord_flip() +
theme_bw() +
labs(
x = "Outcome",
y = "Odds Ratio (95% Confidence Interval)"
) +
geom_text(aes(label = sprintf("p = %.3f", P_value)), vjust = -1.0, size = 3) +
scale_x_discrete(limits = c( "sbp", "rebleed_end_hemo", "death_day2","outcome"), labels = c( "SBP", "Rebleeding", "In hospital mortality", "Composite outcome"))
# Display forest plot
forest_plot_ordered
# Calculate and display c-statistics for each imputed dataset
c_statistics <- lapply(results[[1]]$PS_data, function(data) {
roc_obj <- roc(data$antibio_exposure, data$ps)
auc(roc_obj)
})## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases# Print c-statistics
print(c_statistics)## [[1]]
## Area under the curve: 0.6369
##
## [[2]]
## Area under the curve: 0.6376
##
## [[3]]
## Area under the curve: 0.6383
##
## [[4]]
## Area under the curve: 0.6387
##
## [[5]]
## Area under the curve: 0.6391
##
## [[6]]
## Area under the curve: 0.6352
##
## [[7]]
## Area under the curve: 0.6384
##
## [[8]]
## Area under the curve: 0.6397
##
## [[9]]
## Area under the curve: 0.6413
##
## [[10]]
## Area under the curve: 0.6403
##
## [[11]]
## Area under the curve: 0.6377
##
## [[12]]
## Area under the curve: 0.6421
##
## [[13]]
## Area under the curve: 0.6361
##
## [[14]]
## Area under the curve: 0.6369
##
## [[15]]
## Area under the curve: 0.6393
##
## [[16]]
## Area under the curve: 0.6368
##
## [[17]]
## Area under the curve: 0.6438
##
## [[18]]
## Area under the curve: 0.6345
##
## [[19]]
## Area under the curve: 0.6433
##
## [[20]]
## Area under the curve: 0.6408
##
## [[21]]
## Area under the curve: 0.6356
##
## [[22]]
## Area under the curve: 0.6396
##
## [[23]]
## Area under the curve: 0.6414
##
## [[24]]
## Area under the curve: 0.6397
##
## [[25]]
## Area under the curve: 0.6364
##
## [[26]]
## Area under the curve: 0.6419
##
## [[27]]
## Area under the curve: 0.6399
##
## [[28]]
## Area under the curve: 0.6414
##
## [[29]]
## Area under the curve: 0.6378
##
## [[30]]
## Area under the curve: 0.6396
##
## [[31]]
## Area under the curve: 0.6373
##
## [[32]]
## Area under the curve: 0.6406
##
## [[33]]
## Area under the curve: 0.6384
##
## [[34]]
## Area under the curve: 0.6407
##
## [[35]]
## Area under the curve: 0.6455
##
## [[36]]
## Area under the curve: 0.6424
##
## [[37]]
## Area under the curve: 0.6395
##
## [[38]]
## Area under the curve: 0.6418
##
## [[39]]
## Area under the curve: 0.6378
##
## [[40]]
## Area under the curve: 0.6394
##
## [[41]]
## Area under the curve: 0.6348
##
## [[42]]
## Area under the curve: 0.6393
##
## [[43]]
## Area under the curve: 0.6386
##
## [[44]]
## Area under the curve: 0.6396
##
## [[45]]
## Area under the curve: 0.6387
##
## [[46]]
## Area under the curve: 0.6388
##
## [[47]]
## Area under the curve: 0.6463
##
## [[48]]
## Area under the curve: 0.6396
##
## [[49]]
## Area under the curve: 0.6435
##
## [[50]]
## Area under the curve: 0.641
##
## [[51]]
## Area under the curve: 0.6381
##
## [[52]]
## Area under the curve: 0.6362
##
## [[53]]
## Area under the curve: 0.6356
##
## [[54]]
## Area under the curve: 0.6384
##
## [[55]]
## Area under the curve: 0.6393
##
## [[56]]
## Area under the curve: 0.6398
##
## [[57]]
## Area under the curve: 0.6401
##
## [[58]]
## Area under the curve: 0.6394
##
## [[59]]
## Area under the curve: 0.6493
##
## [[60]]
## Area under the curve: 0.6333
##
## [[61]]
## Area under the curve: 0.6417
##
## [[62]]
## Area under the curve: 0.6392
##
## [[63]]
## Area under the curve: 0.6432
##
## [[64]]
## Area under the curve: 0.6387
##
## [[65]]
## Area under the curve: 0.6388
##
## [[66]]
## Area under the curve: 0.635
##
## [[67]]
## Area under the curve: 0.6415
##
## [[68]]
## Area under the curve: 0.6404
##
## [[69]]
## Area under the curve: 0.6477
##
## [[70]]
## Area under the curve: 0.6381
##
## [[71]]
## Area under the curve: 0.6458
##
## [[72]]
## Area under the curve: 0.6444
##
## [[73]]
## Area under the curve: 0.6347
##
## [[74]]
## Area under the curve: 0.6392
##
## [[75]]
## Area under the curve: 0.6376
##
## [[76]]
## Area under the curve: 0.644
##
## [[77]]
## Area under the curve: 0.6424
##
## [[78]]
## Area under the curve: 0.6379
##
## [[79]]
## Area under the curve: 0.6416
##
## [[80]]
## Area under the curve: 0.6394
##
## [[81]]
## Area under the curve: 0.6367
##
## [[82]]
## Area under the curve: 0.6356
##
## [[83]]
## Area under the curve: 0.6489
##
## [[84]]
## Area under the curve: 0.6401
##
## [[85]]
## Area under the curve: 0.6452
##
## [[86]]
## Area under the curve: 0.6425
##
## [[87]]
## Area under the curve: 0.6414
##
## [[88]]
## Area under the curve: 0.6399
##
## [[89]]
## Area under the curve: 0.6382
##
## [[90]]
## Area under the curve: 0.6403
##
## [[91]]
## Area under the curve: 0.638
##
## [[92]]
## Area under the curve: 0.6424
##
## [[93]]
## Area under the curve: 0.6382
##
## [[94]]
## Area under the curve: 0.6372
##
## [[95]]
## Area under the curve: 0.6414
##
## [[96]]
## Area under the curve: 0.6449
##
## [[97]]
## Area under the curve: 0.6406
##
## [[98]]
## Area under the curve: 0.6402
##
## [[99]]
## Area under the curve: 0.641
##
## [[100]]
## Area under the curve: 0.6382# 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)
}
負の二項分布で解析してみる
# Install and load the MASS package
if (!requireNamespace("MASS", quietly = TRUE)) {
install.packages("MASS")
}
library(MASS)##
## 次のパッケージを付け加えます: 'MASS'## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## select# LOS
# 本解析: IPTW models on imputed data with Negative Binomial distribution
glm_results <- list()
summary_tables <- list()
conf_ints <- list()
weighted_data <- list()
for (i in 1:100) {
imputed_data <- complete(impglm, i)
# Fit GLM propensity score model
propensity_model <- glm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hd + hcc + alocohol + past_varix_rup + antiplate_user_01 + anticoag_user01 + nsaids_user01 + steroid_user01 + ppi_user01 + beta_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt_cate + aptt_cate,
data = imputed_data,
family = binomial("logit"))
# 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
glm_results[[i]] <- iptw_model
summary_tables[[i]] <- summary(iptw_model)
conf_ints[[i]] <- confint(iptw_model, level = 0.95, method = "Wald")
# Store weighted data
imputed_data$weight <- stabilized_weights
weighted_data[[i]] <- imputed_data
}## Waiting for profiling to be done...## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...# Save the results from each imputed dataset
params <- list()
for (i in 1:100) {
params[[i]] <- data.frame(estimate = coef(glm_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.059285cat("\n95% Confidence Interval: (", pooled_lower, ", ", pooled_upper, ")\n")##
## 95% Confidence Interval: ( 0.9797491 , 1.145279 )cat("\nP-value:", pooled_p, "\n")##
## P-value: 0.1481046補完後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)
mean_smd_df$Mean_SMD<- round(abs(mean_smd_df$Mean_SMD), 4)
mean_smd_df %>% mutate(Imputation ="After IPTW")->mean_smd_df
mean_smd_df %>% rename(SMD=Mean_SMD)->mean_smd_df
rownames(mean_smd_df) <- NULL
print(mean_smd_df)## Covariate SMD Imputation
## 1 steroid_user01 0.0712 After IPTW
## 2 child_numl 0.0413 After IPTW
## 3 cci_num 0.0392 After IPTW
## 4 crp 0.0323 After IPTW
## 5 t_bil 0.0241 After IPTW
## 6 wbc 0.0241 After IPTW
## 7 alb_cate 0.0155 After IPTW
## 8 age_cate 0.0113 After IPTW
## 9 hb 0.0077 After IPTW
## 10 map_day0 0.0068 After IPTW
## 11 aptt_cate 0.0049 After IPTW
## 12 pt_cate 0.0046 After IPTW
## 13 plt 0.0044 After IPTW
## 14 child_score 0.0029 After IPTW
## 15 alt 0.0015 After IPTW
## 16 brthel_cate 0.0015 After IPTW
## 17 ast 0.0008 After IPTW
## 18 sex 0.0000 After IPTW
## 19 malignancy 0.0000 After IPTW
## 20 hd 0.0000 After IPTW
## 21 hcc 0.0000 After IPTW
## 22 alocohol 0.0000 After IPTW
## 23 past_varix_rup 0.0000 After IPTW
## 24 antiplate_user_01 0.0000 After IPTW
## 25 anticoag_user01 0.0000 After IPTW
## 26 nsaids_user01 0.0000 After IPTW
## 27 ppi_user01 0.0000 After IPTW
## 28 beta_user01 0.0000 After IPTW
## 29 vaso 0.0000 After IPTW
## 30 shock 0.0000 After IPTW
## 31 eGFR30 0.0000 After IPTWplot作成のためのpreデータ作成 before2.csvから引っ張る
# パッケージの読み込みは不要
before_table <-read.csv("before2.csv", col.names = c("covariate", "non_user", "user", "smd", "missing"))
before_table %>% filter(!is.na(smd))-> before_table
before_table %>% mutate(Imputation ="Before IPTW")->before_table
before_table %>% dplyr::select(covariate,SMD=smd,Imputation)->before_table
# 新しい行名を指定
new_row_names <- c("child_numl", "cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "crp", "age_cate", "sex", "brthel_cate", "malignancy", "child_score", "hd", "hcc", "alocohol", "past_varix_rup", "antiplate_user_01", "anticoag_user01", "nsaids_user01", "steroid_user01", "ppi_user01", "beta_user01", "vaso", "shock","alb_cate", "pt_cate", "aptt_cate", "eGFR30")
# 行名を変更
before_table %>% mutate(covariate, new_row_names)->before_table
before_table %>% dplyr::select(Covariate=new_row_names,SMD,Imputation)->before_table
before_table## Covariate SMD Imputation
## 1 child_numl 0.054 Before IPTW
## 2 cci_num 0.008 Before IPTW
## 3 map_day0 0.093 Before IPTW
## 4 t_bil 0.129 Before IPTW
## 5 ast 0.086 Before IPTW
## 6 alt 0.082 Before IPTW
## 7 wbc 0.154 Before IPTW
## 8 hb 0.119 Before IPTW
## 9 plt 0.001 Before IPTW
## 10 crp 0.030 Before IPTW
## 11 age_cate 0.183 Before IPTW
## 12 sex 0.077 Before IPTW
## 13 brthel_cate 0.101 Before IPTW
## 14 malignancy 0.026 Before IPTW
## 15 child_score 0.058 Before IPTW
## 16 hd 0.012 Before IPTW
## 17 hcc 0.096 Before IPTW
## 18 alocohol 0.214 Before IPTW
## 19 past_varix_rup 0.102 Before IPTW
## 20 antiplate_user_01 0.012 Before IPTW
## 21 anticoag_user01 0.054 Before IPTW
## 22 nsaids_user01 0.012 Before IPTW
## 23 steroid_user01 0.085 Before IPTW
## 24 ppi_user01 0.154 Before IPTW
## 25 beta_user01 0.244 Before IPTW
## 26 vaso 0.022 Before IPTW
## 27 shock 0.100 Before IPTW
## 28 alb_cate 0.111 Before IPTW
## 29 pt_cate 0.017 Before IPTW
## 30 aptt_cate 0.084 Before IPTW
## 31 eGFR30 0.034 Before IPTWpre,postの結合
smd_df_pre_post_main <- rbind(before_table, mean_smd_df)
smd_df_pre_post_main## Covariate SMD Imputation
## 1 child_numl 0.0540 Before IPTW
## 2 cci_num 0.0080 Before IPTW
## 3 map_day0 0.0930 Before IPTW
## 4 t_bil 0.1290 Before IPTW
## 5 ast 0.0860 Before IPTW
## 6 alt 0.0820 Before IPTW
## 7 wbc 0.1540 Before IPTW
## 8 hb 0.1190 Before IPTW
## 9 plt 0.0010 Before IPTW
## 10 crp 0.0300 Before IPTW
## 11 age_cate 0.1830 Before IPTW
## 12 sex 0.0770 Before IPTW
## 13 brthel_cate 0.1010 Before IPTW
## 14 malignancy 0.0260 Before IPTW
## 15 child_score 0.0580 Before IPTW
## 16 hd 0.0120 Before IPTW
## 17 hcc 0.0960 Before IPTW
## 18 alocohol 0.2140 Before IPTW
## 19 past_varix_rup 0.1020 Before IPTW
## 20 antiplate_user_01 0.0120 Before IPTW
## 21 anticoag_user01 0.0540 Before IPTW
## 22 nsaids_user01 0.0120 Before IPTW
## 23 steroid_user01 0.0850 Before IPTW
## 24 ppi_user01 0.1540 Before IPTW
## 25 beta_user01 0.2440 Before IPTW
## 26 vaso 0.0220 Before IPTW
## 27 shock 0.1000 Before IPTW
## 28 alb_cate 0.1110 Before IPTW
## 29 pt_cate 0.0170 Before IPTW
## 30 aptt_cate 0.0840 Before IPTW
## 31 eGFR30 0.0340 Before IPTW
## 32 steroid_user01 0.0712 After IPTW
## 33 child_numl 0.0413 After IPTW
## 34 cci_num 0.0392 After IPTW
## 35 crp 0.0323 After IPTW
## 36 t_bil 0.0241 After IPTW
## 37 wbc 0.0241 After IPTW
## 38 alb_cate 0.0155 After IPTW
## 39 age_cate 0.0113 After IPTW
## 40 hb 0.0077 After IPTW
## 41 map_day0 0.0068 After IPTW
## 42 aptt_cate 0.0049 After IPTW
## 43 pt_cate 0.0046 After IPTW
## 44 plt 0.0044 After IPTW
## 45 child_score 0.0029 After IPTW
## 46 alt 0.0015 After IPTW
## 47 brthel_cate 0.0015 After IPTW
## 48 ast 0.0008 After IPTW
## 49 sex 0.0000 After IPTW
## 50 malignancy 0.0000 After IPTW
## 51 hd 0.0000 After IPTW
## 52 hcc 0.0000 After IPTW
## 53 alocohol 0.0000 After IPTW
## 54 past_varix_rup 0.0000 After IPTW
## 55 antiplate_user_01 0.0000 After IPTW
## 56 anticoag_user01 0.0000 After IPTW
## 57 nsaids_user01 0.0000 After IPTW
## 58 ppi_user01 0.0000 After IPTW
## 59 beta_user01 0.0000 After IPTW
## 60 vaso 0.0000 After IPTW
## 61 shock 0.0000 After IPTW
## 62 eGFR30 0.0000 After IPTWpre,postの結合と名称変更
smd_df_pre_post_main <- rbind(before_table, mean_smd_df)
# Create a mapping of covariate names
covariate_name_mapping <- c(
"age_cate"="Age",
"child_numl" = "Child-Pugh score",
"cci_num" = "Charlson Comorbidity Index",
"map_day0" = "RBC transfusion",
"t_bil" = "Total bilirubin",
"ast" = "Aspartate aminotransferase",
"alt" = "Alanine aminotransferase",
"wbc" = "White blood cell",
"hb" = "Hemoglobin",
"plt" = "Platelet",
"alb_cate" = "Albumin",
"crp" = "C-reactive protein",
"pt_cate" = "Prothrombin time",
"aptt_cate" = "Activated partial thromboplastin time",
"sex" = "Sex",
"brthel_cate" = "Barthel index",
"malignancy" = "Malignant tomor history",
"child_score" = "Child-Pugh classification",
"hd" = "Maintenance hemodialysis",
"hcc" = "Hepatic cancer",
"alocohol" = "Alcohol related disease",
"past_varix_rup" = "Past varix rupture history",
"antiplate_user_01" = "Antiplatelet use",
"anticoag_user01" = "Anticoagulant use",
"nsaids_user01" = "NSAIDs use",
"steroid_user01" = "Corticosteroid use",
"ppi_user01" = "Acid blocker use",
"beta_user01" = "β blocker use",
"vaso" = "Vasopressor",
"shock" = "Shock index > 1",
"eGFR30" = "eGFR < 30"
)
# Replace covariate names with the new names
smd_df_pre_post_main$Covariate <- covariate_name_mapping[smd_df_pre_post_main$Covariate]
smd_df_pre_post_main## Covariate SMD Imputation
## 1 Child-Pugh score 0.0540 Before IPTW
## 2 Charlson Comorbidity Index 0.0080 Before IPTW
## 3 RBC transfusion 0.0930 Before IPTW
## 4 Total bilirubin 0.1290 Before IPTW
## 5 Aspartate aminotransferase 0.0860 Before IPTW
## 6 Alanine aminotransferase 0.0820 Before IPTW
## 7 White blood cell 0.1540 Before IPTW
## 8 Hemoglobin 0.1190 Before IPTW
## 9 Platelet 0.0010 Before IPTW
## 10 C-reactive protein 0.0300 Before IPTW
## 11 Age 0.1830 Before IPTW
## 12 Sex 0.0770 Before IPTW
## 13 Barthel index 0.1010 Before IPTW
## 14 Malignant tomor history 0.0260 Before IPTW
## 15 Child-Pugh classification 0.0580 Before IPTW
## 16 Maintenance hemodialysis 0.0120 Before IPTW
## 17 Hepatic cancer 0.0960 Before IPTW
## 18 Alcohol related disease 0.2140 Before IPTW
## 19 Past varix rupture history 0.1020 Before IPTW
## 20 Antiplatelet use 0.0120 Before IPTW
## 21 Anticoagulant use 0.0540 Before IPTW
## 22 NSAIDs use 0.0120 Before IPTW
## 23 Corticosteroid use 0.0850 Before IPTW
## 24 Acid blocker use 0.1540 Before IPTW
## 25 β blocker use 0.2440 Before IPTW
## 26 Vasopressor 0.0220 Before IPTW
## 27 Shock index > 1 0.1000 Before IPTW
## 28 Albumin 0.1110 Before IPTW
## 29 Prothrombin time 0.0170 Before IPTW
## 30 Activated partial thromboplastin time 0.0840 Before IPTW
## 31 eGFR < 30 0.0340 Before IPTW
## 32 Corticosteroid use 0.0712 After IPTW
## 33 Child-Pugh score 0.0413 After IPTW
## 34 Charlson Comorbidity Index 0.0392 After IPTW
## 35 C-reactive protein 0.0323 After IPTW
## 36 Total bilirubin 0.0241 After IPTW
## 37 White blood cell 0.0241 After IPTW
## 38 Albumin 0.0155 After IPTW
## 39 Age 0.0113 After IPTW
## 40 Hemoglobin 0.0077 After IPTW
## 41 RBC transfusion 0.0068 After IPTW
## 42 Activated partial thromboplastin time 0.0049 After IPTW
## 43 Prothrombin time 0.0046 After IPTW
## 44 Platelet 0.0044 After IPTW
## 45 Child-Pugh classification 0.0029 After IPTW
## 46 Alanine aminotransferase 0.0015 After IPTW
## 47 Barthel index 0.0015 After IPTW
## 48 Aspartate aminotransferase 0.0008 After IPTW
## 49 Sex 0.0000 After IPTW
## 50 Malignant tomor history 0.0000 After IPTW
## 51 Maintenance hemodialysis 0.0000 After IPTW
## 52 Hepatic cancer 0.0000 After IPTW
## 53 Alcohol related disease 0.0000 After IPTW
## 54 Past varix rupture history 0.0000 After IPTW
## 55 Antiplatelet use 0.0000 After IPTW
## 56 Anticoagulant use 0.0000 After IPTW
## 57 NSAIDs use 0.0000 After IPTW
## 58 Acid blocker use 0.0000 After IPTW
## 59 β blocker use 0.0000 After IPTW
## 60 Vasopressor 0.0000 After IPTW
## 61 Shock index > 1 0.0000 After IPTW
## 62 eGFR < 30 0.0000 After IPTWlove plot作成
# Sort the dataframe by absolute SMD of Pre IPTW in descending order
smd_df_pre_post_main_sorted <- smd_df_pre_post_main %>%
filter(Imputation == "Before IPTW") %>%
arrange(SMD) %>%
mutate(Covariate = factor(Covariate, levels = Covariate))
# Merge sorted Pre IPTW data with Post IPTW data
smd_df_pre_post_main_sorted <- smd_df_pre_post_main %>%
filter(Imputation == "After IPTW") %>%
mutate(Covariate = factor(Covariate, levels = smd_df_pre_post_main_sorted$Covariate)) %>%
bind_rows(smd_df_pre_post_main_sorted)
# Create Love plot
love_plot_main <- ggplot(smd_df_pre_post_main_sorted, aes(x = Covariate, y = SMD, color = Imputation)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
coord_flip() +
theme_bw() +
labs(
x = "Covariate",
y = "Absolute Standardized Mean Difference"
) +
guides(color = guide_legend(title = NULL))
# Show the love plot
print(love_plot_main)