stabilized_ATE_main_pt_aptt修正
市田 親正
2023-05-13
必要なlibraryの読み込み
library(tidyverse)
library(readxl)
library(tableone)
library(lubridate)
library(skimr)
library(svglite)
library(haven)
library(summarytools)## Warning in system2("/usr/bin/otool", c("-L", shQuote(DSO)), stdout = TRUE):
## 命令 ''/usr/bin/otool' -L
## '/Library/Frameworks/R.framework/Resources/library/tcltk/libs//tcltk.so''
## の実行は状態 1 を持ちましたlibrary(naniar)
library(devtools)
library(reader)
library(stringr)
library(missForest)
library(mice)
library(geepack)
library(ggplot2)
library(cobalt)
library(tableone)
library(sandwich)
library(survey)
library(WeightIt)
library(MatchIt)
library(twang)
library(naniar)
pacman::p_load(norm2, miceadds, lmtest, car, ROCR, pROC)データ読み込み
data <- read_excel("varix_icu_ok_1.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),
vaso=factor(vaso),
map_day0= as.integer(map_day0),
shock=factor(shock),
pt_cate= factor(pt_cate, levels = c("0", "1", "2")),
aptt_cate=factor(aptt_cate,levels = c("0", "1", "2", "3")),
eGFR30=factor(eGFR30),
outcome=factor(outcome),
death_day2=factor(death_day2),
rebleed_end_hemo=factor(rebleed_end_hemo),
sbp=factor(sbp),
infection=factor(infection),
infection_all=factor(infection_all),
los=as.integer(los)
)素データでのtable1の作成
col_continuous <-
c("child_numl","cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "alb", "crp")
col_factors <-
c("age_cate","sex", "brthel_cate", "malignancy", "child_score", "hd", "hcc", "alocohol", "past_varix_rup", "antiplate_user_01", "anticoag_user01", "nsaids_user01" ,"steroid_user01", "ppi_user01","vaso", "shock", "pt_cate", "aptt_cate", "eGFR30")
df|> select(c(col_continuous, col_factors, "antibio_exposure"))|>
CreateTableOne(vars = c(col_continuous, col_factors), strata = "antibio_exposure", factorVars = col_factors) -> tableone## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(col_continuous)
##
## # Now:
## data %>% select(all_of(col_continuous))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(col_factors)
##
## # Now:
## data %>% select(all_of(col_factors))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.print(tableone, smd = T, missing = T, test = T, explain = T) ## Stratified by antibio_exposure
## 0 1 p test
## n 558 232
## child_numl (mean (SD)) 8.30 (1.96) 8.41 (2.03) 0.508
## cci_num (mean (SD)) 4.51 (1.32) 4.52 (1.33) 0.917
## map_day0 (mean (SD)) 2.86 (2.92) 3.13 (2.98) 0.231
## t_bil (mean (SD)) 2.01 (1.88) 2.26 (2.03) 0.098
## ast (mean (SD)) 74.97 (82.14) 83.16 (106.55) 0.248
## alt (mean (SD)) 38.83 (42.49) 42.57 (49.00) 0.285
## wbc (mean (SD)) 8194.22 (4242.15) 8945.32 (5438.02) 0.039
## hb (mean (SD)) 8.69 (2.47) 8.98 (2.39) 0.133
## plt (mean (SD)) 116.49 (61.73) 116.46 (84.84) 0.996
## alb (mean (SD)) 2.91 (0.57) 2.92 (0.64) 0.789
## crp (mean (SD)) 0.76 (1.39) 0.72 (1.57) 0.699
## age_cate (%) 0.235
## 0 322 (57.7) 143 (61.6)
## 1 144 (25.8) 46 (19.8)
## 2 90 (16.1) 43 (18.5)
## 3 2 ( 0.4) 0 ( 0.0)
## sex = F (%) 141 (25.3) 51 (22.0) 0.374
## brthel_cate (%) 0.486
## 0 186 (41.5) 83 (40.5)
## 1 152 (33.9) 63 (30.7)
## 2 110 (24.6) 59 (28.8)
## malignancy = 1 (%) 65 (11.6) 29 (12.5) 0.829
## child_score (%) 0.776
## a 93 (18.8) 42 (19.4)
## b 266 (53.7) 110 (50.9)
## c 136 (27.5) 64 (29.6)
## hd = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## hcc = 1 (%) 112 (20.1) 38 (16.4) 0.269
## alocohol = 1 (%) 246 (44.1) 127 (54.7) 0.008
## past_varix_rup = 1 (%) 127 (22.8) 63 (27.2) 0.221
## antiplate_user_01 = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## anticoag_user01 = 1 (%) 8 ( 1.4) 5 ( 2.2) 0.675
## nsaids_user01 = 1 (%) 13 ( 2.3) 5 ( 2.2) 1.000
## steroid_user01 = 1 (%) 2 ( 0.4) 0 ( 0.0) 0.892
## ppi_user01 = 1 (%) 486 (87.1) 213 (91.8) 0.077
## vaso = 1 (%) 19 ( 3.4) 7 ( 3.0) 0.953
## shock = 1 (%) 197 (36.5) 94 (41.4) 0.236
## pt_cate (%) 0.977
## 0 113 (21.7) 48 (21.5)
## 1 324 (62.3) 138 (61.9)
## 2 83 (16.0) 37 (16.6)
## aptt_cate (%) 0.757
## 0 436 (89.3) 185 (87.3)
## 1 47 ( 9.6) 23 (10.8)
## 2 1 ( 0.2) 1 ( 0.5)
## 3 4 ( 0.8) 3 ( 1.4)
## eGFR30 = 1 (%) 40 ( 7.3) 19 ( 8.3) 0.769
## Stratified by antibio_exposure
## SMD Missing
## n
## child_numl (mean (SD)) 0.054 12.8
## cci_num (mean (SD)) 0.008 0.0
## map_day0 (mean (SD)) 0.093 0.0
## t_bil (mean (SD)) 0.129 3.2
## ast (mean (SD)) 0.086 2.2
## alt (mean (SD)) 0.082 2.2
## wbc (mean (SD)) 0.154 1.8
## hb (mean (SD)) 0.119 1.8
## plt (mean (SD)) <0.001 1.8
## alb (mean (SD)) 0.021 4.6
## crp (mean (SD)) 0.030 4.7
## age_cate (%) 0.170 0.0
## 0
## 1
## 2
## 3
## sex = F (%) 0.077 0.0
## brthel_cate (%) 0.101 17.3
## 0
## 1
## 2
## malignancy = 1 (%) 0.026 0.0
## child_score (%) 0.058 10.0
## a
## b
## c
## hd = 1 (%) 0.012 0.0
## hcc = 1 (%) 0.096 0.0
## alocohol = 1 (%) 0.214 0.0
## past_varix_rup = 1 (%) 0.102 0.0
## antiplate_user_01 = 1 (%) 0.012 0.0
## anticoag_user01 = 1 (%) 0.054 0.0
## nsaids_user01 = 1 (%) 0.012 0.0
## steroid_user01 = 1 (%) 0.085 0.0
## ppi_user01 = 1 (%) 0.154 0.0
## vaso = 1 (%) 0.022 0.0
## shock = 1 (%) 0.100 3.0
## pt_cate (%) 0.017 5.9
## 0
## 1
## 2
## aptt_cate (%) 0.085 11.4
## 0
## 1
## 2
## 3
## eGFR30 = 1 (%) 0.034 1.9多重代入
imp1 <- mice(df,
m=100, #mは代入の結果として作成するデータセットの数
maxit =200, #maxitは反復回数
method="pmm",
printFlag = FALSE,
seed = 10
)## Warning: Number of logged events: 400000主要評価解析
# Function to perform IPTW analysis on different outcomes
perform_iptw_analysis <- function(outcome_name) {
gee_results <- list()
summary_tables <- list()
conf_ints <- list()
weighted_data <- list()
ps_data <- list()
for (i in 1:100) {
imputed_data <- complete(imp1, i)
# Fit GEE propensity score model
propensity_model <- geeglm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hd + hcc + alocohol + past_varix_rup + antiplate_user_01 + anticoag_user01 + nsaids_user01 + steroid_user01 + ppi_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt + aptt,
data = imputed_data,
id = hosp_id,
family = binomial("logit"),
corstr = "exchangeable")
# Calculate propensity scores
ps <- predict(propensity_model, type = "response", newdata = imputed_data)
imputed_data$ps <- ps
# Calculate IPTW and fit the logistic regression model
stabilized_weights <- with(imputed_data, antibio_exposure/(ps) + (1 - antibio_exposure)/(1 - ps))
iptw_model <- glm(as.formula(paste0(outcome_name, " ~ antibio_exposure")), data = imputed_data, weights = stabilized_weights, family = binomial("logit"))
# Store the results
gee_results[[i]] <- iptw_model
summary_tables[[i]] <- summary(iptw_model)
conf_ints[[i]] <- confint(iptw_model, level = 0.95, method = "Wald")
# Store weighted data
imputed_data$weight <- stabilized_weights
weighted_data[[i]] <- imputed_data
ps_data[[i]] <- imputed_data[, c("ps", "antibio_exposure")]
}
# Save the results from each imputed dataset
params <- list()
for (i in 1:100) {
params[[i]] <- data.frame(Estimate = coef(gee_results[[i]])[2],
Std.Error = coef(summary_tables[[i]])[2, "Std. Error"],
df = Inf)
}
# Calculate the pooled estimates
Qbar <- mean(sapply(params, function(x) x$Estimate))
se_Qbar <- sqrt(sum(sapply(params, function(x) x$Std.Error^2))/length(params))
# Calculate the pooled odds ratio and confidence intervals
pooled_or <- exp(Qbar)
pooled_lower <- exp(Qbar - 1.96 * se_Qbar)
pooled_upper <- exp(Qbar + 1.96 * se_Qbar)
# Calculate the pooled p-value (two-tailed)
pooled_p <- 2 * (1 - pnorm(abs(Qbar / se_Qbar)))
return(list(OR = pooled_or, Lower_CI = pooled_lower, Upper_CI = pooled_upper, P_value = pooled_p, PS_data = ps_data))
}
# Perform IPTW analysis for each outcome
outcomes <- c("outcome", "death_day2", "rebleed_end_hemo", "sbp")
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...# 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.2690801 0.8983942 1.792714 0.1763568
## 2 death_day2 0.9165571 0.6060258 1.386206 0.6797521
## 3 rebleed_end_hemo 1.7026495 0.8576406 3.380222 0.1282418
## 4 sbp 1.5962912 0.8076067 3.155181 0.1785185# 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.6269
##
## [[2]]
## Area under the curve: 0.6355
##
## [[3]]
## Area under the curve: 0.6311
##
## [[4]]
## Area under the curve: 0.6395
##
## [[5]]
## Area under the curve: 0.6352
##
## [[6]]
## Area under the curve: 0.6227
##
## [[7]]
## Area under the curve: 0.6281
##
## [[8]]
## Area under the curve: 0.6323
##
## [[9]]
## Area under the curve: 0.6254
##
## [[10]]
## Area under the curve: 0.6286
##
## [[11]]
## Area under the curve: 0.6348
##
## [[12]]
## Area under the curve: 0.6294
##
## [[13]]
## Area under the curve: 0.6333
##
## [[14]]
## Area under the curve: 0.6308
##
## [[15]]
## Area under the curve: 0.6352
##
## [[16]]
## Area under the curve: 0.6307
##
## [[17]]
## Area under the curve: 0.6315
##
## [[18]]
## Area under the curve: 0.6302
##
## [[19]]
## Area under the curve: 0.6234
##
## [[20]]
## Area under the curve: 0.6314
##
## [[21]]
## Area under the curve: 0.6338
##
## [[22]]
## Area under the curve: 0.6283
##
## [[23]]
## Area under the curve: 0.6249
##
## [[24]]
## Area under the curve: 0.6334
##
## [[25]]
## Area under the curve: 0.6317
##
## [[26]]
## Area under the curve: 0.6299
##
## [[27]]
## Area under the curve: 0.6353
##
## [[28]]
## Area under the curve: 0.6342
##
## [[29]]
## Area under the curve: 0.63
##
## [[30]]
## Area under the curve: 0.6285
##
## [[31]]
## Area under the curve: 0.6351
##
## [[32]]
## Area under the curve: 0.6263
##
## [[33]]
## Area under the curve: 0.6236
##
## [[34]]
## Area under the curve: 0.6238
##
## [[35]]
## Area under the curve: 0.6288
##
## [[36]]
## Area under the curve: 0.6314
##
## [[37]]
## Area under the curve: 0.6327
##
## [[38]]
## Area under the curve: 0.6329
##
## [[39]]
## Area under the curve: 0.6298
##
## [[40]]
## Area under the curve: 0.6273
##
## [[41]]
## Area under the curve: 0.6246
##
## [[42]]
## Area under the curve: 0.6297
##
## [[43]]
## Area under the curve: 0.6271
##
## [[44]]
## Area under the curve: 0.6298
##
## [[45]]
## Area under the curve: 0.6357
##
## [[46]]
## Area under the curve: 0.6335
##
## [[47]]
## Area under the curve: 0.6374
##
## [[48]]
## Area under the curve: 0.633
##
## [[49]]
## Area under the curve: 0.6282
##
## [[50]]
## Area under the curve: 0.6347
##
## [[51]]
## Area under the curve: 0.6267
##
## [[52]]
## Area under the curve: 0.6278
##
## [[53]]
## Area under the curve: 0.635
##
## [[54]]
## Area under the curve: 0.6239
##
## [[55]]
## Area under the curve: 0.6239
##
## [[56]]
## Area under the curve: 0.63
##
## [[57]]
## Area under the curve: 0.6289
##
## [[58]]
## Area under the curve: 0.6339
##
## [[59]]
## Area under the curve: 0.6277
##
## [[60]]
## Area under the curve: 0.6313
##
## [[61]]
## Area under the curve: 0.6359
##
## [[62]]
## Area under the curve: 0.6301
##
## [[63]]
## Area under the curve: 0.6272
##
## [[64]]
## Area under the curve: 0.6335
##
## [[65]]
## Area under the curve: 0.6281
##
## [[66]]
## Area under the curve: 0.6256
##
## [[67]]
## Area under the curve: 0.6287
##
## [[68]]
## Area under the curve: 0.63
##
## [[69]]
## Area under the curve: 0.625
##
## [[70]]
## Area under the curve: 0.6331
##
## [[71]]
## Area under the curve: 0.6281
##
## [[72]]
## Area under the curve: 0.6298
##
## [[73]]
## Area under the curve: 0.6298
##
## [[74]]
## Area under the curve: 0.628
##
## [[75]]
## Area under the curve: 0.628
##
## [[76]]
## Area under the curve: 0.6295
##
## [[77]]
## Area under the curve: 0.6243
##
## [[78]]
## Area under the curve: 0.6259
##
## [[79]]
## Area under the curve: 0.6341
##
## [[80]]
## Area under the curve: 0.626
##
## [[81]]
## Area under the curve: 0.6279
##
## [[82]]
## Area under the curve: 0.6315
##
## [[83]]
## Area under the curve: 0.6278
##
## [[84]]
## Area under the curve: 0.6285
##
## [[85]]
## Area under the curve: 0.6296
##
## [[86]]
## Area under the curve: 0.6331
##
## [[87]]
## Area under the curve: 0.6298
##
## [[88]]
## Area under the curve: 0.6295
##
## [[89]]
## Area under the curve: 0.6314
##
## [[90]]
## Area under the curve: 0.6268
##
## [[91]]
## Area under the curve: 0.6305
##
## [[92]]
## Area under the curve: 0.6281
##
## [[93]]
## Area under the curve: 0.6339
##
## [[94]]
## Area under the curve: 0.6319
##
## [[95]]
## Area under the curve: 0.6303
##
## [[96]]
## Area under the curve: 0.6272
##
## [[97]]
## Area under the curve: 0.6296
##
## [[98]]
## Area under the curve: 0.6325
##
## [[99]]
## Area under the curve: 0.6294
##
## [[100]]
## Area under the curve: 0.6246# 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
gee_results <- list()
summary_tables <- list()
conf_ints <- list()
weighted_data <- list()
for (i in 1:100) {
imputed_data <- complete(imp1, i)
# Fit GEE propensity score model
propensity_model <- geeglm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hd + hcc + alocohol + past_varix_rup + antiplate_user_01 + anticoag_user01 + nsaids_user01 + steroid_user01 + ppi_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt + aptt,
data = imputed_data,
id = hosp_id,
family = binomial("logit"),
corstr = "exchangeable")
# Calculate propensity scores
ps <- predict(propensity_model, type = "response", newdata = imputed_data)
imputed_data$ps <- ps
# Calculate IPTW and fit the negative binomial regression model
stabilized_weights <- with(imputed_data, antibio_exposure/(ps) + (1 - antibio_exposure)/(1 - ps))
iptw_model <- glm.nb(los ~ antibio_exposure, data = imputed_data, weights = stabilized_weights)
# Store the results
gee_results[[i]] <- iptw_model
summary_tables[[i]] <- summary(iptw_model)
conf_ints[[i]] <- confint(iptw_model, level = 0.95, method = "Wald")
# Store weighted data
imputed_data$weight <- stabilized_weights
weighted_data[[i]] <- imputed_data
}## Waiting for profiling to be done...## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...# Save the results from each imputed dataset
params <- list()
for (i in 1:100) {
params[[i]] <- data.frame(estimate = coef(gee_results[[i]])[2],
std.error = coef(summary_tables[[i]])[2, "Std. Error"],
df = Inf) # Add infinite degrees of freedom
}
# Manually pool the results
Qbar <- mean(sapply(params, function(x) x$estimate))
se_Qbar <- sqrt(sum(sapply(params, function(x) x$std.error^2))/length(params))
# Calculate the pooled rate ratio and confidence intervals
pooled_rr <- exp(Qbar)
pooled_lower <- exp(Qbar - 1.96 * se_Qbar)
pooled_upper <- exp(Qbar + 1.96 * se_Qbar)
# Calculate the pooled z-value
pooled_z <- Qbar / se_Qbar
# Calculate the pooled p-value (two-tailed)
pooled_p <- 2 * (1 - pnorm(abs(pooled_z)))
cat("\nPooled results (IPTW with Negative Binomial distribution)\n")##
## Pooled results (IPTW with Negative Binomial distribution)cat("\nRate Ratio:", pooled_rr, "\n")##
## Rate Ratio: 1.036824cat("\n95% Confidence Interval: (", pooled_lower, ", ", pooled_upper, ")\n")##
## 95% Confidence Interval: ( 0.9591131 , 1.120832 )cat("\nP-value:", pooled_p, "\n")##
## P-value: 0.3629483補完後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 ="Post 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.0713 Post IPTW
## 2 crp 0.0467 Post IPTW
## 3 aptt_cate 0.0464 Post IPTW
## 4 alb 0.0440 Post IPTW
## 5 child_numl 0.0435 Post IPTW
## 6 t_bil 0.0271 Post IPTW
## 7 wbc 0.0209 Post IPTW
## 8 cci_num 0.0195 Post IPTW
## 9 map_day0 0.0150 Post IPTW
## 10 hb 0.0122 Post IPTW
## 11 ast 0.0073 Post IPTW
## 12 pt_cate 0.0064 Post IPTW
## 13 plt 0.0046 Post IPTW
## 14 age_cate 0.0038 Post IPTW
## 15 brthel_cate 0.0016 Post IPTW
## 16 alt 0.0010 Post IPTW
## 17 child_score 0.0002 Post IPTW
## 18 sex 0.0000 Post IPTW
## 19 malignancy 0.0000 Post IPTW
## 20 hd 0.0000 Post IPTW
## 21 hcc 0.0000 Post IPTW
## 22 alocohol 0.0000 Post IPTW
## 23 past_varix_rup 0.0000 Post IPTW
## 24 antiplate_user_01 0.0000 Post IPTW
## 25 anticoag_user01 0.0000 Post IPTW
## 26 nsaids_user01 0.0000 Post IPTW
## 27 ppi_user01 0.0000 Post IPTW
## 28 vaso 0.0000 Post IPTW
## 29 shock 0.0000 Post IPTW
## 30 eGFR30 0.0000 Post 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 ="Pre IPTW")->before_table
before_table %>% dplyr::select(covariate,SMD=smd,Imputation)->before_table
# 新しい行名を指定
new_row_names <- c("child_numl", "cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "alb", "crp", "age_cate", "sex", "brthel_cate", "malignancy", "child_score", "hd", "hcc", "alocohol", "past_varix_rup", "antiplate_user_01", "anticoag_user01", "nsaids_user01", "steroid_user01", "ppi_user01", "vaso", "shock", "pt_cate", "aptt_cate", "eGFR30")
# 行名を変更
before_table %>% mutate(covariate, new_row_names)->before_table
before_table %>% dplyr::select(Covariate=new_row_names,SMD,Imputation)->before_table
before_table## Covariate SMD Imputation
## 1 child_numl 0.054 Pre IPTW
## 2 cci_num 0.008 Pre IPTW
## 3 map_day0 0.093 Pre IPTW
## 4 t_bil 0.129 Pre IPTW
## 5 ast 0.086 Pre IPTW
## 6 alt 0.082 Pre IPTW
## 7 wbc 0.154 Pre IPTW
## 8 hb 0.119 Pre IPTW
## 9 plt 0.001 Pre IPTW
## 10 alb 0.021 Pre IPTW
## 11 crp 0.030 Pre IPTW
## 12 age_cate 0.170 Pre IPTW
## 13 sex 0.077 Pre IPTW
## 14 brthel_cate 0.101 Pre IPTW
## 15 malignancy 0.026 Pre IPTW
## 16 child_score 0.058 Pre IPTW
## 17 hd 0.012 Pre IPTW
## 18 hcc 0.096 Pre IPTW
## 19 alocohol 0.214 Pre IPTW
## 20 past_varix_rup 0.102 Pre IPTW
## 21 antiplate_user_01 0.012 Pre IPTW
## 22 anticoag_user01 0.054 Pre IPTW
## 23 nsaids_user01 0.012 Pre IPTW
## 24 steroid_user01 0.085 Pre IPTW
## 25 ppi_user01 0.154 Pre IPTW
## 26 vaso 0.022 Pre IPTW
## 27 shock 0.100 Pre IPTW
## 28 pt_cate 0.017 Pre IPTW
## 29 aptt_cate 0.085 Pre IPTW
## 30 eGFR30 0.034 Pre 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 Pre IPTW
## 2 cci_num 0.0080 Pre IPTW
## 3 map_day0 0.0930 Pre IPTW
## 4 t_bil 0.1290 Pre IPTW
## 5 ast 0.0860 Pre IPTW
## 6 alt 0.0820 Pre IPTW
## 7 wbc 0.1540 Pre IPTW
## 8 hb 0.1190 Pre IPTW
## 9 plt 0.0010 Pre IPTW
## 10 alb 0.0210 Pre IPTW
## 11 crp 0.0300 Pre IPTW
## 12 age_cate 0.1700 Pre IPTW
## 13 sex 0.0770 Pre IPTW
## 14 brthel_cate 0.1010 Pre IPTW
## 15 malignancy 0.0260 Pre IPTW
## 16 child_score 0.0580 Pre IPTW
## 17 hd 0.0120 Pre IPTW
## 18 hcc 0.0960 Pre IPTW
## 19 alocohol 0.2140 Pre IPTW
## 20 past_varix_rup 0.1020 Pre IPTW
## 21 antiplate_user_01 0.0120 Pre IPTW
## 22 anticoag_user01 0.0540 Pre IPTW
## 23 nsaids_user01 0.0120 Pre IPTW
## 24 steroid_user01 0.0850 Pre IPTW
## 25 ppi_user01 0.1540 Pre IPTW
## 26 vaso 0.0220 Pre IPTW
## 27 shock 0.1000 Pre IPTW
## 28 pt_cate 0.0170 Pre IPTW
## 29 aptt_cate 0.0850 Pre IPTW
## 30 eGFR30 0.0340 Pre IPTW
## 31 steroid_user01 0.0713 Post IPTW
## 32 crp 0.0467 Post IPTW
## 33 aptt_cate 0.0464 Post IPTW
## 34 alb 0.0440 Post IPTW
## 35 child_numl 0.0435 Post IPTW
## 36 t_bil 0.0271 Post IPTW
## 37 wbc 0.0209 Post IPTW
## 38 cci_num 0.0195 Post IPTW
## 39 map_day0 0.0150 Post IPTW
## 40 hb 0.0122 Post IPTW
## 41 ast 0.0073 Post IPTW
## 42 pt_cate 0.0064 Post IPTW
## 43 plt 0.0046 Post IPTW
## 44 age_cate 0.0038 Post IPTW
## 45 brthel_cate 0.0016 Post IPTW
## 46 alt 0.0010 Post IPTW
## 47 child_score 0.0002 Post IPTW
## 48 sex 0.0000 Post IPTW
## 49 malignancy 0.0000 Post IPTW
## 50 hd 0.0000 Post IPTW
## 51 hcc 0.0000 Post IPTW
## 52 alocohol 0.0000 Post IPTW
## 53 past_varix_rup 0.0000 Post IPTW
## 54 antiplate_user_01 0.0000 Post IPTW
## 55 anticoag_user01 0.0000 Post IPTW
## 56 nsaids_user01 0.0000 Post IPTW
## 57 ppi_user01 0.0000 Post IPTW
## 58 vaso 0.0000 Post IPTW
## 59 shock 0.0000 Post IPTW
## 60 eGFR30 0.0000 Post IPTWpre,postの結合と名称変更
smd_df_pre_post_main <- rbind(before_table, mean_smd_df)
# Create a mapping of covariate names
covariate_name_mapping <- c(
"age_cate"="Age",
"child_numl" = "Child-Pugh score",
"cci_num" = "Charlson Comorbidity Index",
"map_day0" = "RBC transfusion",
"t_bil" = "Total bilirubin",
"ast" = "Aspartate aminotransferase",
"alt" = "Alanine aminotransferase",
"wbc" = "White blood cell",
"hb" = "Hemoglobin",
"plt" = "Platelet",
"alb" = "Albumin",
"crp" = "C-reactive protein",
"pt_cate" = "Prothrombin time",
"aptt_cate" = "Activated partial thromboplastin time",
"sex" = "Sex",
"brthel_cate" = "Barthel index",
"malignancy" = "Malignant tomor history",
"child_score" = "Child-Pugh classification",
"hd" = "Maintenance hemodialysis",
"hcc" = "Hepatic cancer",
"alocohol" = "Alcohol related disease",
"past_varix_rup" = "Past varix rupture history",
"antiplate_user_01" = "Antiplatelet use",
"anticoag_user01" = "Anticoagulant use",
"nsaids_user01" = "NSAIDs use",
"steroid_user01" = "Corticosteroid use",
"ppi_user01" = "Acid blocker use",
"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 Pre IPTW
## 2 Charlson Comorbidity Index 0.0080 Pre IPTW
## 3 RBC transfusion 0.0930 Pre IPTW
## 4 Total bilirubin 0.1290 Pre IPTW
## 5 Aspartate aminotransferase 0.0860 Pre IPTW
## 6 Alanine aminotransferase 0.0820 Pre IPTW
## 7 White blood cell 0.1540 Pre IPTW
## 8 Hemoglobin 0.1190 Pre IPTW
## 9 Platelet 0.0010 Pre IPTW
## 10 Albumin 0.0210 Pre IPTW
## 11 C-reactive protein 0.0300 Pre IPTW
## 12 Age 0.1700 Pre IPTW
## 13 Sex 0.0770 Pre IPTW
## 14 Barthel index 0.1010 Pre IPTW
## 15 Malignant tomor history 0.0260 Pre IPTW
## 16 Child-Pugh classification 0.0580 Pre IPTW
## 17 Maintenance hemodialysis 0.0120 Pre IPTW
## 18 Hepatic cancer 0.0960 Pre IPTW
## 19 Alcohol related disease 0.2140 Pre IPTW
## 20 Past varix rupture history 0.1020 Pre IPTW
## 21 Antiplatelet use 0.0120 Pre IPTW
## 22 Anticoagulant use 0.0540 Pre IPTW
## 23 NSAIDs use 0.0120 Pre IPTW
## 24 Corticosteroid use 0.0850 Pre IPTW
## 25 Acid blocker use 0.1540 Pre IPTW
## 26 Vasopressor 0.0220 Pre IPTW
## 27 Shock index > 1 0.1000 Pre IPTW
## 28 Prothrombin time 0.0170 Pre IPTW
## 29 Activated partial thromboplastin time 0.0850 Pre IPTW
## 30 eGFR < 30 0.0340 Pre IPTW
## 31 Corticosteroid use 0.0713 Post IPTW
## 32 C-reactive protein 0.0467 Post IPTW
## 33 Activated partial thromboplastin time 0.0464 Post IPTW
## 34 Albumin 0.0440 Post IPTW
## 35 Child-Pugh score 0.0435 Post IPTW
## 36 Total bilirubin 0.0271 Post IPTW
## 37 White blood cell 0.0209 Post IPTW
## 38 Charlson Comorbidity Index 0.0195 Post IPTW
## 39 RBC transfusion 0.0150 Post IPTW
## 40 Hemoglobin 0.0122 Post IPTW
## 41 Aspartate aminotransferase 0.0073 Post IPTW
## 42 Prothrombin time 0.0064 Post IPTW
## 43 Platelet 0.0046 Post IPTW
## 44 Age 0.0038 Post IPTW
## 45 Barthel index 0.0016 Post IPTW
## 46 Alanine aminotransferase 0.0010 Post IPTW
## 47 Child-Pugh classification 0.0002 Post IPTW
## 48 Sex 0.0000 Post IPTW
## 49 Malignant tomor history 0.0000 Post IPTW
## 50 Maintenance hemodialysis 0.0000 Post IPTW
## 51 Hepatic cancer 0.0000 Post IPTW
## 52 Alcohol related disease 0.0000 Post IPTW
## 53 Past varix rupture history 0.0000 Post IPTW
## 54 Antiplatelet use 0.0000 Post IPTW
## 55 Anticoagulant use 0.0000 Post IPTW
## 56 NSAIDs use 0.0000 Post IPTW
## 57 Acid blocker use 0.0000 Post IPTW
## 58 Vasopressor 0.0000 Post IPTW
## 59 Shock index > 1 0.0000 Post IPTW
## 60 eGFR < 30 0.0000 Post 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 == "Pre 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 == "Post 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)