PSM_メイン_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_2.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),
pt_cate= factor(pt_cate, levels = c("0", "1", "2")),
aptt_cate=factor(aptt_cate,levels = c("0", "1", "2", "3")),
eGFR30=factor(eGFR30),
outcome=factor(outcome),
death_day2=factor(death_day2),
rebleed_end_hemo=factor(rebleed_end_hemo),
sbp=factor(sbp),
infection=factor(infection),
infection_all=factor(infection_all),
los=as.integer(los),
cdi=factor(cdi)
)素データでのtable1の作成
col_continuous <-
c("child_numl","cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "alb", "crp")
col_factors <-
c("age_cate","sex", "brthel_cate", "malignancy", "child_score", "hd", "hcc", "alocohol", "past_varix_rup", "antiplate_user_01", "anticoag_user01", "nsaids_user01" ,"steroid_user01", "ppi_user01","beta_user01", "vaso", "shock", "pt_cate", "aptt_cate", "eGFR30")
df|> select(c(col_continuous, col_factors, "antibio_exposure"))|>
CreateTableOne(vars = c(col_continuous, col_factors), strata = "antibio_exposure", factorVars = col_factors) -> tableone## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(col_continuous)
##
## # Now:
## data %>% select(all_of(col_continuous))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(col_factors)
##
## # Now:
## data %>% select(all_of(col_factors))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.print(tableone, smd = T, missing = T, test = T, explain = T) ## Stratified by antibio_exposure
## 0 1 p test
## n 558 232
## child_numl (mean (SD)) 8.30 (1.96) 8.41 (2.03) 0.508
## cci_num (mean (SD)) 4.51 (1.32) 4.52 (1.33) 0.917
## map_day0 (mean (SD)) 2.86 (2.92) 3.13 (2.98) 0.231
## t_bil (mean (SD)) 2.01 (1.88) 2.26 (2.03) 0.098
## ast (mean (SD)) 74.97 (82.14) 83.16 (106.55) 0.248
## alt (mean (SD)) 38.83 (42.49) 42.57 (49.00) 0.285
## wbc (mean (SD)) 8194.22 (4242.15) 8945.32 (5438.02) 0.039
## hb (mean (SD)) 8.69 (2.47) 8.98 (2.39) 0.133
## plt (mean (SD)) 116.49 (61.73) 116.46 (84.84) 0.996
## alb (mean (SD)) 2.91 (0.57) 2.92 (0.64) 0.789
## crp (mean (SD)) 0.76 (1.39) 0.72 (1.57) 0.699
## age_cate (%) 0.155
## 0 322 (57.7) 143 (61.6)
## 1 144 (25.8) 46 (19.8)
## 2 75 (13.4) 39 (16.8)
## 3 17 ( 3.0) 4 ( 1.7)
## sex = F (%) 141 (25.3) 51 (22.0) 0.374
## brthel_cate (%) 0.486
## 0 186 (41.5) 83 (40.5)
## 1 152 (33.9) 63 (30.7)
## 2 110 (24.6) 59 (28.8)
## malignancy = 1 (%) 65 (11.6) 29 (12.5) 0.829
## child_score (%) 0.776
## a 93 (18.8) 42 (19.4)
## b 266 (53.7) 110 (50.9)
## c 136 (27.5) 64 (29.6)
## hd = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## hcc = 1 (%) 112 (20.1) 38 (16.4) 0.269
## alocohol = 1 (%) 246 (44.1) 127 (54.7) 0.008
## past_varix_rup = 1 (%) 127 (22.8) 63 (27.2) 0.221
## antiplate_user_01 = 1 (%) 8 ( 1.4) 3 ( 1.3) 1.000
## anticoag_user01 = 1 (%) 8 ( 1.4) 5 ( 2.2) 0.675
## nsaids_user01 = 1 (%) 13 ( 2.3) 5 ( 2.2) 1.000
## steroid_user01 = 1 (%) 2 ( 0.4) 0 ( 0.0) 0.892
## ppi_user01 = 1 (%) 486 (87.1) 213 (91.8) 0.077
## beta_user01 = 1 (%) 26 ( 4.7) 26 (11.2) 0.001
## vaso = 1 (%) 19 ( 3.4) 7 ( 3.0) 0.953
## shock = 1 (%) 197 (36.5) 94 (41.4) 0.236
## pt_cate (%) 0.977
## 0 113 (21.7) 48 (21.5)
## 1 324 (62.3) 138 (61.9)
## 2 83 (16.0) 37 (16.6)
## aptt_cate (%) 0.757
## 0 436 (89.3) 185 (87.3)
## 1 47 ( 9.6) 23 (10.8)
## 2 1 ( 0.2) 1 ( 0.5)
## 3 4 ( 0.8) 3 ( 1.4)
## eGFR30 = 1 (%) 40 ( 7.3) 19 ( 8.3) 0.769
## Stratified by antibio_exposure
## SMD Missing
## n
## child_numl (mean (SD)) 0.054 12.8
## cci_num (mean (SD)) 0.008 0.0
## map_day0 (mean (SD)) 0.093 0.0
## t_bil (mean (SD)) 0.129 3.2
## ast (mean (SD)) 0.086 2.2
## alt (mean (SD)) 0.082 2.2
## wbc (mean (SD)) 0.154 1.8
## hb (mean (SD)) 0.119 1.8
## plt (mean (SD)) <0.001 1.8
## alb (mean (SD)) 0.021 4.6
## crp (mean (SD)) 0.030 4.7
## age_cate (%) 0.183 0.0
## 0
## 1
## 2
## 3
## sex = F (%) 0.077 0.0
## brthel_cate (%) 0.101 17.3
## 0
## 1
## 2
## malignancy = 1 (%) 0.026 0.0
## child_score (%) 0.058 10.0
## a
## b
## c
## hd = 1 (%) 0.012 0.0
## hcc = 1 (%) 0.096 0.0
## alocohol = 1 (%) 0.214 0.0
## past_varix_rup = 1 (%) 0.102 0.0
## antiplate_user_01 = 1 (%) 0.012 0.0
## anticoag_user01 = 1 (%) 0.054 0.0
## nsaids_user01 = 1 (%) 0.012 0.0
## steroid_user01 = 1 (%) 0.085 0.0
## ppi_user01 = 1 (%) 0.154 0.0
## beta_user01 = 1 (%) 0.244 0.0
## vaso = 1 (%) 0.022 0.0
## shock = 1 (%) 0.100 3.0
## pt_cate (%) 0.017 5.9
## 0
## 1
## 2
## aptt_cate (%) 0.085 11.4
## 0
## 1
## 2
## 3
## eGFR30 = 1 (%) 0.034 1.9多重代入
#imp1 <- mice(df,
# m=5, #mは代入の結果として作成するデータセットの数
# maxit = 5, #maxitは反復回数
# method="pmm",
# printFlag = FALSE,
# seed = 10
# )多重代入のロード
load("imp.Rdata")傾向スコアマッチング メインアウトカム
# Function to perform propensity score matching analysis on different outcomes
perform_psm_analysis <- function(outcome_name) {
matching_results <- list()
summary_tables <- list()
conf_ints <- list()
matched_data <- list()
ps_data <- list()
for (i in 1:100) {
imputed_data <- complete(imp1, i)
# Fit GEE propensity score model
propensity_model <- geeglm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hd + hcc + alocohol + past_varix_rup + antiplate_user_01 + anticoag_user01 + nsaids_user01 + steroid_user01 + ppi_user01 + beta_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt + aptt,
data = imputed_data,
id = hosp_id,
family = binomial("logit"),
corstr = "exchangeable")
# Calculate propensity scores
ps <- predict(propensity_model, type = "response", newdata = imputed_data)
imputed_data$ps <- ps
# Perform propensity score matching
match_it <- matchit(antibio_exposure ~ ps, data = imputed_data, method = "nearest", caliper = 0.2, replace = FALSE)
# Extract the matched data
matched_data[[i]] <- match.data(match_it)
# Fit logistic regression model on the matched data
logistic_model <- glm(as.formula(paste0(outcome_name, " ~ antibio_exposure")), data = matched_data[[i]], family = binomial("logit"))
# Store the results
matching_results[[i]] <- logistic_model
summary_tables[[i]] <- summary(logistic_model)
conf_ints[[i]] <- confint(logistic_model, level = 0.95, method = "Wald")
# Store the propensity score data
ps_data[[i]] <- matched_data[[i]][, c("ps", "antibio_exposure")]
}
# Save the results from each imputed dataset
params <- list()
for (i in 1:100) {
params[[i]] <- data.frame(Estimate = coef(matching_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 PSM analysis for each outcome
outcomes <- c("outcome", "death_day2", "rebleed_end_hemo", "sbp", "cdi")
results_psm <- lapply(outcomes, perform_psm_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_psm <- data.frame(
Outcome = outcomes,
OR = sapply(results_psm, function(x) x$OR),
Lower_CI = sapply(results_psm, function(x) x$Lower_CI),
Upper_CI = sapply(results_psm, function(x) x$Upper_CI),
P_value = sapply(results_psm, function(x) x$P_value)
)
# Reorder forest_plot_data_psm by the desired order of outcomes
forest_plot_data_psm$Outcome <- factor(forest_plot_data_psm$Outcome, levels = c("outcome", "death_day2", "rebleed_end_hemo", "sbp", "cdi"))
# Display results as a table
print(forest_plot_data_psm)## Outcome OR Lower_CI Upper_CI P_value
## 1 outcome 1.089212e+00 5.729463e-01 2.070668e+00 0.7943098
## 2 death_day2 8.593616e-01 4.025717e-01 1.834462e+00 0.6952450
## 3 rebleed_end_hemo 1.392099e+00 3.883179e-01 4.990604e+00 0.6115591
## 4 sbp 1.763475e+00 9.719782e-167 3.199498e+166 0.9976826
## 5 cdi 2.732668e+05 0.000000e+00 Inf 0.9963505# Create forest plot with the specified order and adjusted p-value position
forest_plot_ordered_psm <- ggplot(forest_plot_data_psm, 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_psm
# Save the forest plot as a PNG
ggsave(filename = "PSM_main_forest.png", plot = forest_plot_ordered_psm, width = 8, height = 6, units = "in", dpi = 300)
# Calculate and display c-statistics for each imputed dataset
c_statistics_psm <- lapply(results_psm[[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_psm)## [[1]]
## Area under the curve: 0.4975
##
## [[2]]
## Area under the curve: 0.5037
##
## [[3]]
## Area under the curve: 0.4964
##
## [[4]]
## Area under the curve: 0.5017
##
## [[5]]
## Area under the curve: 0.5023
##
## [[6]]
## Area under the curve: 0.4984
##
## [[7]]
## Area under the curve: 0.5022
##
## [[8]]
## Area under the curve: 0.4979
##
## [[9]]
## Area under the curve: 0.504
##
## [[10]]
## Area under the curve: 0.5044
##
## [[11]]
## Area under the curve: 0.5012
##
## [[12]]
## Area under the curve: 0.5033
##
## [[13]]
## Area under the curve: 0.5032
##
## [[14]]
## Area under the curve: 0.4979
##
## [[15]]
## Area under the curve: 0.5014
##
## [[16]]
## Area under the curve: 0.4983
##
## [[17]]
## Area under the curve: 0.499
##
## [[18]]
## Area under the curve: 0.5036
##
## [[19]]
## Area under the curve: 0.5023
##
## [[20]]
## Area under the curve: 0.4949
##
## [[21]]
## Area under the curve: 0.4985
##
## [[22]]
## Area under the curve: 0.5015
##
## [[23]]
## Area under the curve: 0.5032
##
## [[24]]
## Area under the curve: 0.4981
##
## [[25]]
## Area under the curve: 0.5027
##
## [[26]]
## Area under the curve: 0.5029
##
## [[27]]
## Area under the curve: 0.5047
##
## [[28]]
## Area under the curve: 0.4962
##
## [[29]]
## Area under the curve: 0.5021
##
## [[30]]
## Area under the curve: 0.5066
##
## [[31]]
## Area under the curve: 0.4983
##
## [[32]]
## Area under the curve: 0.5022
##
## [[33]]
## Area under the curve: 0.5035
##
## [[34]]
## Area under the curve: 0.4967
##
## [[35]]
## Area under the curve: 0.5025
##
## [[36]]
## Area under the curve: 0.4977
##
## [[37]]
## Area under the curve: 0.5037
##
## [[38]]
## Area under the curve: 0.5042
##
## [[39]]
## Area under the curve: 0.506
##
## [[40]]
## Area under the curve: 0.5019
##
## [[41]]
## Area under the curve: 0.5039
##
## [[42]]
## Area under the curve: 0.5057
##
## [[43]]
## Area under the curve: 0.5016
##
## [[44]]
## Area under the curve: 0.4978
##
## [[45]]
## Area under the curve: 0.4974
##
## [[46]]
## Area under the curve: 0.4975
##
## [[47]]
## Area under the curve: 0.5026
##
## [[48]]
## Area under the curve: 0.5032
##
## [[49]]
## Area under the curve: 0.5025
##
## [[50]]
## Area under the curve: 0.5019
##
## [[51]]
## Area under the curve: 0.4972
##
## [[52]]
## Area under the curve: 0.5039
##
## [[53]]
## Area under the curve: 0.503
##
## [[54]]
## Area under the curve: 0.5023
##
## [[55]]
## Area under the curve: 0.5027
##
## [[56]]
## Area under the curve: 0.5023
##
## [[57]]
## Area under the curve: 0.5017
##
## [[58]]
## Area under the curve: 0.4987
##
## [[59]]
## Area under the curve: 0.4975
##
## [[60]]
## Area under the curve: 0.5029
##
## [[61]]
## Area under the curve: 0.5033
##
## [[62]]
## Area under the curve: 0.5013
##
## [[63]]
## Area under the curve: 0.5029
##
## [[64]]
## Area under the curve: 0.4961
##
## [[65]]
## Area under the curve: 0.5035
##
## [[66]]
## Area under the curve: 0.5023
##
## [[67]]
## Area under the curve: 0.5031
##
## [[68]]
## Area under the curve: 0.4975
##
## [[69]]
## Area under the curve: 0.5024
##
## [[70]]
## Area under the curve: 0.5025
##
## [[71]]
## Area under the curve: 0.5015
##
## [[72]]
## Area under the curve: 0.5019
##
## [[73]]
## Area under the curve: 0.5037
##
## [[74]]
## Area under the curve: 0.5029
##
## [[75]]
## Area under the curve: 0.4972
##
## [[76]]
## Area under the curve: 0.5028
##
## [[77]]
## Area under the curve: 0.5035
##
## [[78]]
## Area under the curve: 0.4976
##
## [[79]]
## Area under the curve: 0.503
##
## [[80]]
## Area under the curve: 0.4985
##
## [[81]]
## Area under the curve: 0.5032
##
## [[82]]
## Area under the curve: 0.5024
##
## [[83]]
## Area under the curve: 0.5034
##
## [[84]]
## Area under the curve: 0.499
##
## [[85]]
## Area under the curve: 0.5047
##
## [[86]]
## Area under the curve: 0.5022
##
## [[87]]
## Area under the curve: 0.4957
##
## [[88]]
## Area under the curve: 0.5018
##
## [[89]]
## Area under the curve: 0.503
##
## [[90]]
## Area under the curve: 0.498
##
## [[91]]
## Area under the curve: 0.5026
##
## [[92]]
## Area under the curve: 0.4973
##
## [[93]]
## Area under the curve: 0.4971
##
## [[94]]
## Area under the curve: 0.5042
##
## [[95]]
## Area under the curve: 0.4989
##
## [[96]]
## Area under the curve: 0.5027
##
## [[97]]
## Area under the curve: 0.4975
##
## [[98]]
## Area under the curve: 0.4983
##
## [[99]]
## Area under the curve: 0.4961
##
## [[100]]
## Area under the curve: 0.5037# Plot the propensity score overlap for each imputed dataset
for (i in 1:5) {
ggplot(results_psm[[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)
}
before matchindのSMD算出 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 Matching")->before_table
before_table %>% dplyr::select(covariate,SMD=smd,Imputation)->before_table
# 新しい行名を指定
new_row_names <- c("child_numl", "cci_num", "map_day0", "t_bil", "ast", "alt", "wbc", "hb", "plt", "alb", "crp", "age_cate", "sex", "brthel_cate", "malignancy", "child_score", "hd", "hcc", "alocohol", "past_varix_rup", "antiplate_user_01", "anticoag_user01", "nsaids_user01", "steroid_user01", "ppi_user01", "beta_user01", "vaso", "shock", "pt_cate", "aptt_cate", "eGFR30")
# 行名を変更
before_table %>% mutate(covariate, new_row_names)->before_table
before_table %>% dplyr::select(Covariate=new_row_names,SMD,Imputation)->before_table
before_table## Covariate SMD Imputation
## 1 child_numl 0.054 Before Matching
## 2 cci_num 0.008 Before Matching
## 3 map_day0 0.093 Before Matching
## 4 t_bil 0.129 Before Matching
## 5 ast 0.086 Before Matching
## 6 alt 0.082 Before Matching
## 7 wbc 0.154 Before Matching
## 8 hb 0.119 Before Matching
## 9 plt 0.001 Before Matching
## 10 alb 0.021 Before Matching
## 11 crp 0.030 Before Matching
## 12 age_cate 0.183 Before Matching
## 13 sex 0.077 Before Matching
## 14 brthel_cate 0.101 Before Matching
## 15 malignancy 0.026 Before Matching
## 16 child_score 0.058 Before Matching
## 17 hd 0.012 Before Matching
## 18 hcc 0.096 Before Matching
## 19 alocohol 0.214 Before Matching
## 20 past_varix_rup 0.102 Before Matching
## 21 antiplate_user_01 0.012 Before Matching
## 22 anticoag_user01 0.054 Before Matching
## 23 nsaids_user01 0.012 Before Matching
## 24 steroid_user01 0.085 Before Matching
## 25 ppi_user01 0.154 Before Matching
## 26 beta_user01 0.244 Before Matching
## 27 vaso 0.022 Before Matching
## 28 shock 0.100 Before Matching
## 29 pt_cate 0.017 Before Matching
## 30 aptt_cate 0.085 Before Matching
## 31 eGFR30 0.034 Before MatchingafterSMD作成
library(effsize)
library(dplyr)
library(ggplot2)
# smd関数の定義
smd <- function(x, y) {
m_x <- mean(x, na.rm = TRUE)
m_y <- mean(y, na.rm = TRUE)
sd_x <- sd(x, na.rm = TRUE)
sd_y <- sd(y, na.rm = TRUE)
pooled_sd <- sqrt(((length(x) - 1) * sd_x^2 + (length(y) - 1) * sd_y^2) / (length(x) + length(y) - 2))
return(abs(m_x - m_y) / pooled_sd)
}
# 標準化平均差 (SMD) を計算する関数 (カテゴリ変数対応版)
calc_smd <- function(data_list, treat_col, continuous_covariates, factor_covariates) {
treat_data <- data_list[[1]]
control_data <- data_list[[2]]
smd_values <- c()
for (var in continuous_covariates) {
smd_values <- c(smd_values, smd(treat_data[[var]], control_data[[var]]))
}
for (var in factor_covariates) {
treat_table <- table(treat_data[[var]], useNA = "no")
control_table <- table(control_data[[var]], useNA = "no")
# Ensure both tables have the same levels
all_levels <- sort(unique(c(levels(treat_data[[var]]), levels(control_data[[var]]))))
treat_table <- treat_table[all_levels]
control_table <- control_table[all_levels]
treat_prop <- treat_table / sum(treat_table)
control_prop <- control_table / sum(control_table)
# Compute SMD for categorical variables
smd_cat <- sum((treat_prop - control_prop)^2) / 4
smd_values <- c(smd_values, smd_cat)
}
names(smd_values) <- c(continuous_covariates, factor_covariates)
return(smd_values)
}
# マッチング前後のSMDを計算する
smd_before <- calc_smd(list(df[df$antibio_exposure == 1, ], df[df$antibio_exposure == 0, ]), "antibio_exposure", col_continuous, col_factors)
# perform_psm_analysis関数
perform_psm_analysis <- function(outcome_name) {
matching_results <- list()
summary_tables <- list()
conf_ints <- list()
matched_data <- list()
ps_data <- list()
for (i in 1:100) {
imputed_data <- complete(imp1, i)
# Fit GEE propensity score model
propensity_model <- geeglm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hd + hcc + alocohol + past_varix_rup + antiplate_user_01 + anticoag_user01 + nsaids_user01 + steroid_user01 + ppi_user01+ beta_user01 + vaso + map_day0 + 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
# Perform propensity score matching
match_it <- matchit(antibio_exposure ~ ps, data = imputed_data, method = "nearest", caliper = 0.2, replace = FALSE)
# Extract the matched data
matched_data[[i]] <- match.data(match_it)
# Fit logistic regression model on the matched data
logistic_model <- glm(as.formula(paste0(outcome_name, " ~ antibio_exposure")), data = matched_data[[i]], family = binomial("logit"))
# Store the results
matching_results[[i]] <- logistic_model
summary_tables[[i]] <- summary(logistic_model)
conf_ints[[i]] <- confint(logistic_model, level = 0.95, method = "Wald")
# Store the propensity score data
ps_data[[i]] <- matched_data[[i]][, c("ps", "antibio_exposure")]
}
# Save the results from each imputed dataset
params <- list()
for (i in 1:100) {
params[[i]] <- data.frame(Estimate = coef(matching_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, matched_data = matched_data))
}
# Perform PSM analysis for each outcome
outcomes <- c("outcome", "death_day2", "rebleed_end_hemo", "sbp", "cdi")
results_psm <- lapply(outcomes, perform_psm_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...# Extract matched_data from results_psm and flatten the list
matched_data <- unlist(lapply(results_psm, function(x) x$matched_data), recursive = FALSE)
# マッチング前後のSMDを計算する
smd_before <- calc_smd(list(df[df$antibio_exposure == 1, ], df[df$antibio_exposure == 0, ]), "antibio_exposure", col_continuous, col_factors)
# 各代入データに対するマッチング後のSMDを計算し、平均を求める
smd_after_list <- lapply(matched_data, function(x) {
treat_data <- subset(x, x[["antibio_exposure"]] == 1)
control_data <- subset(x, x[["antibio_exposure"]] == 0)
calc_smd(list(treat_data, control_data), "antibio_exposure", col_continuous, col_factors)
})
smd_after <- colMeans(do.call(rbind, smd_after_list))after SMD作成続き
smd_data_after <- data.frame(
Covariate = names(smd_after),
SMD = unlist(smd_after),
Imputation = "After Matching"
)
smd_df_pre_post2 <- rbind(before_table, smd_data_after)
# Create Love plot
love_plot <- ggplot(smd_df_pre_post2, 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)
# SMDをテーブル形式で出力する
smd_table <- data.frame(
Covariate = c(names(smd_after)),
SMD_after = c(rep(NA, length(smd_before)), round(smd_after, 4))
)
# Show the SMD table
print(smd_table)## Covariate SMD_after
## 1 child_numl NA
## 2 cci_num NA
## 3 map_day0 NA
## 4 t_bil NA
## 5 ast NA
## 6 alt NA
## 7 wbc NA
## 8 hb NA
## 9 plt NA
## 10 alb NA
## 11 crp NA
## 12 age_cate NA
## 13 sex NA
## 14 brthel_cate NA
## 15 malignancy NA
## 16 child_score NA
## 17 hd NA
## 18 hcc NA
## 19 alocohol NA
## 20 past_varix_rup NA
## 21 antiplate_user_01 NA
## 22 anticoag_user01 NA
## 23 nsaids_user01 NA
## 24 steroid_user01 NA
## 25 ppi_user01 NA
## 26 beta_user01 NA
## 27 vaso NA
## 28 shock NA
## 29 pt_cate NA
## 30 aptt_cate NA
## 31 eGFR30 NA
## 32 child_numl 0.0429
## 33 cci_num 0.0420
## 34 map_day0 0.0495
## 35 t_bil 0.0416
## 36 ast 0.0360
## 37 alt 0.0398
## 38 wbc 0.0483
## 39 hb 0.0355
## 40 plt 0.0683
## 41 alb 0.0377
## 42 crp 0.0653
## 43 age_cate 0.0008
## 44 sex 0.0003
## 45 brthel_cate 0.0015
## 46 malignancy 0.0003
## 47 child_score 0.0005
## 48 hd 0.0000
## 49 hcc 0.0002
## 50 alocohol 0.0008
## 51 past_varix_rup 0.0015
## 52 antiplate_user_01 0.0000
## 53 anticoag_user01 0.0000
## 54 nsaids_user01 0.0000
## 55 steroid_user01 0.0000
## 56 ppi_user01 0.0001
## 57 beta_user01 0.0007
## 58 vaso 0.0003
## 59 shock 0.0006
## 60 pt_cate 0.0005
## 61 aptt_cate 0.0002
## 62 eGFR30 0.0001love plotの名称・順番を変更
# 名前の変更
covariate_name_mapping <- c(
"age_cate"="Age",
"child_numl" = "Child-Pugh score",
"cci_num" = "Charlson Comorbidity Index",
"map_day0" = "RBC transfusion",
"t_bil" = "Total bilirubin",
"ast" = "Aspartate aminotransferase",
"alt" = "Alanine aminotransferase",
"wbc" = "White blood cell",
"hb" = "Hemoglobin",
"plt" = "Platelet",
"alb" = "Albumin",
"crp" = "C-reactive protein",
"pt_cate" = "Prothrombin time",
"aptt_cate" = "Activated partial thromboplastin time",
"sex" = "Sex",
"brthel_cate" = "Barthel index",
"malignancy" = "Malignant tomor history",
"child_score" = "Child-Pugh classification",
"hd" = "Maintenance hemodialysis",
"hcc" = "Hepatic cancer",
"alocohol" = "Alcohol related disease",
"past_varix_rup" = "Past varix rupture history",
"antiplate_user_01" = "Antiplatelet use",
"anticoag_user01" = "Anticoagulant use",
"nsaids_user01" = "NSAIDs use",
"steroid_user01" = "Corticosteroid use",
"ppi_user01" = "Acid blocker use",
"beta_user01" = "β blocker use",
"vaso" = "Vasopressor",
"shock" = "Shock index > 1",
"eGFR30" = "eGFR < 30"
)
# 行名を変更
before_table <- before_table %>%
mutate(Covariate = recode_factor(Covariate, !!!covariate_name_mapping))
smd_data_after <- smd_data_after %>%
mutate(Covariate = recode_factor(Covariate, !!!covariate_name_mapping))
# SMDデータを結合
smd_df_pre_post3 <- rbind(before_table, smd_data_after)
smd_df_pre_post3 %>% arrange(SMD)->smd_df_pre_post3_sorted
# "Before Matching"のSMDデータを昇順に並べ替え
smd_df_pre_before_sorted <- smd_df_pre_post3 %>%
filter(Imputation == "Before Matching") %>%
arrange(SMD)
# "Before Matching"の順序を保持するための順序ベクトルを作成
before_order <- smd_df_pre_before_sorted$Covariate
# 全体のデータフレームでCovariateの順序を変更
smd_df_pre_post3_ordered <- smd_df_pre_post3 %>%
mutate(Covariate = factor(Covariate, levels = before_order))
# Create Love plot for ordered data
love_plot_ordered <- ggplot(smd_df_pre_post3_ordered, 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_ordered)
負の二項分布を用いてlos算出
library(MASS)##
## 次のパッケージを付け加えます: 'MASS'## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## selectmatched_data <- list()
for (i in 1:100) {
imputed_data <- complete(imp1, i)
# 傾向スコアモデルを推定
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 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt + aptt,
data = imputed_data,
family = binomial("logit"))
# マッチングを実行
m_out <- matchit(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,
method = "nearest",
distance = propensity_model$fitted.values)
# マッチング後のデータセットを取得し、リストに格納
matched_data[[i]] <- match.data(m_out)
}
# 負の二項分布を用いた回帰分析の結果を格納するリスト
nb_results <- list()
# 各マッチング後のデータセットに対して負の二項分布を用いた回帰分析を実行
for (i in 1:100) {
nb_model <- glm.nb(los ~ antibio_exposure, data = matched_data[[i]])
nb_results[[i]] <- nb_model
}
# 各データセットからの結果を格納するリスト
params_nb <- list()
# 各データセットからの結果を取得
for (i in 1:100) {
params_nb[[i]] <- data.frame(estimate = coef(nb_results[[i]])[2],
std.error = coef(summary(nb_results[[i]]))[2, "Std. Error"],
df = Inf) # 無限大の自由度を追加
}
# 結果を手動でプール
Qbar_nb <- mean(sapply(params_nb, function(x) x$estimate))
se_Qbar_nb <- sqrt(sum(sapply(params_nb, function(x) x$std.error^2))/length(params_nb))
# プールされたレート比および信頼区間を計算
pooled_rr_nb <- exp(Qbar_nb)
pooled_lower_nb <- exp(Qbar_nb - 1.96 * se_Qbar_nb)
pooled_upper_nb <- exp(Qbar_nb + 1.96 * se_Qbar_nb)
# プールされたz値を計算
pooled_z_nb <- Qbar_nb / se_Qbar_nb
# プールされたp値（両側）を計算
pooled_p_nb <- 2 * (1 - pnorm(abs(pooled_z_nb)))
cat("\nPooled results (PSM with Negative Binomial distribution)\n")##
## Pooled results (PSM with Negative Binomial distribution)cat("\nRate Ratio:", pooled_rr_nb, "\n")##
## Rate Ratio: 1.041677cat("\n95% Confidence Interval: (", pooled_lower_nb, ", ", pooled_upper_nb, ")\n")##
## 95% Confidence Interval: ( 0.9045576 , 1.199582 )cat("\nP-value:", pooled_p_nb, "\n")##
## P-value: 0.5706982