Untitled

必要なlibraryの読み込み

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(readxl)
library(tableone)
library(lubridate)
library(skimr)
library(svglite)
library(haven)
library(summarytools)

## 
##  次のパッケージを付け加えます: 'summarytools' 
## 
##  以下のオブジェクトは 'package:tibble' からマスクされています:
## 
##     view

library(naniar)

## 
##  次のパッケージを付け加えます: 'naniar' 
## 
##  以下のオブジェクトは 'package:skimr' からマスクされています:
## 
##     n_complete

library(devtools)

##  要求されたパッケージ usethis をロード中です

library(reader)

##  要求されたパッケージ NCmisc をロード中です 
## 
##  次のパッケージを付け加えます: 'reader' 
## 
##  以下のオブジェクトは 'package:NCmisc' からマスクされています:
## 
##     cat.path, get.ext, rmv.ext

library(stringr)
library(missForest)
library(mice)

## 
##  次のパッケージを付け加えます: 'mice' 
## 
##  以下のオブジェクトは 'package:stats' からマスクされています:
## 
##     filter
## 
##  以下のオブジェクトは 'package:base' からマスクされています:
## 
##     cbind, rbind

library(geepack)
library(ggplot2)
library(cobalt)

##  cobalt (Version 4.4.1, Build Date: 2022-11-03)

library(tableone)
library(sandwich)
library(survey)

##  要求されたパッケージ grid をロード中です 
##  要求されたパッケージ Matrix をロード中です 
## 
##  次のパッケージを付け加えます: 'Matrix' 
## 
##  以下のオブジェクトは 'package:tidyr' からマスクされています:
## 
##     expand, pack, unpack
## 
##  要求されたパッケージ survival をロード中です 
## 
##  次のパッケージを付け加えます: 'survey' 
## 
##  以下のオブジェクトは 'package:graphics' からマスクされています:
## 
##     dotchart

library(WeightIt)
library(twang)

## To reproduce results from prior versions of the twang package, please see the version="legacy" option described in the documentation.

library(naniar)
pacman::p_load(norm2, miceadds, lmtest, car, ROCR, pROC,future,furrr)

データ読み込み

 #data <- read_excel("varix_icu_ok_3.xlsx")

data <- read_excel("varix_icu_ok_3_0926.xlsx")

データの整理

df <-
  data|> 　#miceパッケージが順序ありの因子型でないと読めないためこれを指定。
  mutate(
        hosp_id=as.integer(hosp_id),
        year=as.integer(year),
        pt_id=as.integer(pt_id),
        ad_num=as.integer(ad_num),
        antibio_exposure=as.integer(antibio_exposure),　#geeglmを回すにはfactorでは回らない
        antibio_type_a= factor(antibio_type_a, levels = c("a", "b","c","d","e","f","g","h","i","99")),
        antibio_type_b= factor(antibio_type_b, levels = c("c","d","f","h","i","99")),
        antibio_term= as.integer(antibio_term ),
        age_num=as.integer(age_num),
        age_cate=factor(age_cate,levels = c("0", "1", "2", "3")),
        sex= factor(sex, levels = c("M", "F")),
        smoking= as.integer(smoking),
        brthel_cate= factor(brthel_cate, levels = c("0", "1", "2")),
        child_numl= as.integer(child_numl),
        child_score=factor(child_score, levels = c("a", "b", "c")),
        jcs_all=factor(jcs_all, levels = c("0", "1", "2","3","10","20","30","100","200","300")),
        cci_num=as.integer(cci_num),
        malignancy=factor(malignancy),
        hd=factor(hd),
        hcc=factor(hcc),
        alocohol=factor(alocohol),
        past_varix_rup=factor(past_varix_rup),
        antiplate_user_01=factor(antiplate_user_01),
        anticoag_user01=factor(anticoag_user01),
        nsaids_user01=factor(nsaids_user01),
        steroid_user01=factor(steroid_user01),
        ppi_user01=factor(ppi_user01),
        beta_user01=factor(beta_user01),   
        vaso=factor(vaso),
        map_day0= as.integer(map_day0),
        shock=factor(shock),
        alb_cate= factor(alb_cate, levels = c("0", "1", "2")),
        pt_cate= factor(pt_cate, levels = c("0", "1", "2")),
        aptt_cate=factor(aptt_cate,levels = c("0", "1", "2")),
        eGFR30=factor(eGFR30),
        outcome=factor(outcome),
        death_day2=factor(death_day2),
        rebleed_end_hemo=factor(rebleed_end_hemo),
        sbp=factor(sbp),
        infection=factor(infection),
        infection_all=factor(infection_all),
        los=as.integer(los),
        cdi=factor(cdi)
         )

多重代入のロード

load("impglm.Rdata")

imputed_data <- complete(impglm, 1)

abまとめ

library(dplyr)

# `child_score` をフィルタリングして 'a' および 'b' を取得
filtered_data <- imputed_data %>% filter(child_score %in% c("a", "b"))

# Propensity score modelの適合
propensity_model <- glm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hcc + alocohol + past_varix_rup + anticoag_user01 + nsaids_user01 + ppi_user01 + beta_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt_cate + aptt_cate, 
                        data = filtered_data, 
                        family = binomial("logit"))

# Propensity scoresを計算
ps <- predict(propensity_model, type = "response", newdata = filtered_data)

# ATT-IPTWを計算
att_weights <- ifelse(filtered_data$antibio_exposure == 1, 1/ps, 1/(1-ps))
iptw_model <- glm(outcome ~ antibio_exposure, data = filtered_data, weights = att_weights, family = binomial("logit"))

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# 結果を取得
estimate <- coef(iptw_model)[2]
std_error <- coef(summary(iptw_model))[2, "Std. Error"]
conf_int <- confint(iptw_model, level = 0.95, method = "Wald")

## Waiting for profiling to be done...

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# Odds ratioとconfidence intervalsを計算
or <- exp(estimate)
lower <- exp(estimate - 1.96 * std_error)
upper <- exp(estimate + 1.96 * std_error)

# p-valueを計算
p <- 2 * (1 - pnorm(abs(estimate / std_error)))

# 結果をdata frameに保存
result_a_b <- data.frame(
  Child_Score = "a,b",
  OR = or,
  Lower_CI = lower,
  Upper_CI = upper,
  P_value = p
)

print(result_a_b)

##                  Child_Score        OR  Lower_CI Upper_CI   P_value
## antibio_exposure         a,b 0.9009794 0.5535244 1.466537 0.6748441

c

# `child_score` をフィルタリングして 'c' を取得
imputed_data_c <- imputed_data %>% filter(child_score == "c")

# Propensity score modelの適合
propensity_model_c <- glm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hcc + alocohol + past_varix_rup + anticoag_user01 + nsaids_user01 + ppi_user01 + beta_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt_cate + aptt_cate, 
                          data = imputed_data_c, 
                          family = binomial("logit"))

# Propensity scoresを計算
ps_c <- predict(propensity_model_c, type = "response", newdata = imputed_data_c)

# ATT-IPTWを計算
att_weights_c <- ifelse(imputed_data_c$antibio_exposure == 1, 1/ps_c, 1/(1-ps_c))
iptw_model_c <- glm(outcome ~ antibio_exposure, data = imputed_data_c, weights = att_weights_c, family = binomial("logit"))

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# 結果を取得
estimate_c <- coef(iptw_model_c)[2]
std_error_c <- coef(summary(iptw_model_c))[2, "Std. Error"]
conf_int_c <- confint(iptw_model_c, level = 0.95, method = "Wald")

## Waiting for profiling to be done...

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# Odds ratioとconfidence intervalsを計算
or_c <- exp(estimate_c)
lower_c <- exp(estimate_c - 1.96 * std_error_c)
upper_c <- exp(estimate_c + 1.96 * std_error_c)

# p-valueを計算
p_c <- 2 * (1 - pnorm(abs(estimate_c / std_error_c)))

# 結果をdata frameに保存
result_c <- data.frame(
  Child_Score = "c",
  OR = or_c,
  Lower_CI = lower_c,
  Upper_CI = upper_c,
  P_value = p_c
)

print(result_c)

##                  Child_Score       OR Lower_CI Upper_CI     P_value
## antibio_exposure           c 1.905181 1.199943 3.024908 0.006280322

library(dplyr)

# `child_score` をフィルタリングして 'a' を取得
data_a <- imputed_data %>% filter(child_score == "a")

# 各因子変数のレベルの数を計算
factor_vars <- sapply(data_a, is.factor)  # 因子変数を識別
level_counts <- sapply(data_a[, factor_vars], function(x) length(unique(x)))

print(level_counts)

##    antibio_type_a    antibio_type_b          age_cate               sex 
##                 8                 2                 4                 2 
##       brthel_cate       child_score           jcs_all        malignancy 
##                 3                 1                 3                 2 
##                hd               hcc          alocohol    past_varix_rup 
##                 2                 2                 2                 2 
## antiplate_user_01   anticoag_user01     nsaids_user01    steroid_user01 
##                 2                 2                 1                 1 
##        ppi_user01       beta_user01              vaso             shock 
##                 2                 2                 2                 2 
##          alb_cate            eGFR30           pt_cate         aptt_cate 
##                 3                 2                 2                 1 
##           outcome        death_day2  rebleed_end_hemo               sbp 
##                 2                 2                 2                 2 
##         infection     infection_all               cdi 
##                 2                 2                 2

a　ステロイド、nsaids、aptt cateは抜かないと回らない。

# `child_score` をフィルタリングして 'a' を取得
imputed_data_c <- imputed_data %>% filter(child_score == "a")

# Propensity score modelの適合
propensity_model_c <- glm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hcc + alocohol + past_varix_rup + anticoag_user01 +  ppi_user01 + beta_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt_cate , 
                          data = imputed_data_c, 
                          family = binomial("logit"))




# Propensity scoresを計算
ps_c <- predict(propensity_model_c, type = "response", newdata = imputed_data_c)

# ATT-IPTWを計算
att_weights_c <- ifelse(imputed_data_c$antibio_exposure == 1, 1/ps_c, 1/(1-ps_c))
iptw_model_c <- glm(outcome ~ antibio_exposure, data = imputed_data_c, weights = att_weights_c, family = binomial("logit"))

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# 結果を取得
estimate_c <- coef(iptw_model_c)[2]
std_error_c <- coef(summary(iptw_model_c))[2, "Std. Error"]
conf_int_c <- confint(iptw_model_c, level = 0.95, method = "Wald")

## Waiting for profiling to be done...

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# Odds ratioとconfidence intervalsを計算
or_c <- exp(estimate_c)
lower_c <- exp(estimate_c - 1.96 * std_error_c)
upper_c <- exp(estimate_c + 1.96 * std_error_c)

# p-valueを計算
p_c <- 2 * (1 - pnorm(abs(estimate_c / std_error_c)))

# 結果をdata frameに保存
result_c <- data.frame(
  Child_Score = "a",
  OR = or_c,
  Lower_CI = lower_c,
  Upper_CI = upper_c,
  P_value = p_c
)

print(result_c)

##                  Child_Score        OR  Lower_CI Upper_CI   P_value
## antibio_exposure           a 0.8665251 0.2248918 3.338786 0.8350943

b

# `child_score` をフィルタリングして 'b' を取得
imputed_data_c <- imputed_data %>% filter(child_score == "b")

# Propensity score modelの適合
propensity_model_c <- glm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hcc + alocohol + past_varix_rup + anticoag_user01 + nsaids_user01 + ppi_user01 + beta_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt_cate + aptt_cate, 
                          data = imputed_data_c, 
                          family = binomial("logit"))

# Propensity scoresを計算
ps_c <- predict(propensity_model_c, type = "response", newdata = imputed_data_c)

# ATT-IPTWを計算
att_weights_c <- ifelse(imputed_data_c$antibio_exposure == 1, 1/ps_c, 1/(1-ps_c))
iptw_model_c <- glm(outcome ~ antibio_exposure, data = imputed_data_c, weights = att_weights_c, family = binomial("logit"))

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# 結果を取得
estimate_c <- coef(iptw_model_c)[2]
std_error_c <- coef(summary(iptw_model_c))[2, "Std. Error"]
conf_int_c <- confint(iptw_model_c, level = 0.95, method = "Wald")

## Waiting for profiling to be done...

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# Odds ratioとconfidence intervalsを計算
or_c <- exp(estimate_c)
lower_c <- exp(estimate_c - 1.96 * std_error_c)
upper_c <- exp(estimate_c + 1.96 * std_error_c)

# p-valueを計算
p_c <- 2 * (1 - pnorm(abs(estimate_c / std_error_c)))

# 結果をdata frameに保存
result_c <- data.frame(
  Child_Score = "b",
  OR = or_c,
  Lower_CI = lower_c,
  Upper_CI = upper_c,
  P_value = p_c
)

print(result_c)

##                  Child_Score        OR  Lower_CI Upper_CI   P_value
## antibio_exposure           b 0.7932464 0.4564749 1.378476 0.4113442

# `child_score` をフィルタリングして 'c' を取得
imputed_data_c <- imputed_data %>% filter(child_score == "b")

# Propensity score modelの適合
propensity_model_c <- glm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hcc + alocohol + past_varix_rup + anticoag_user01 + nsaids_user01 + ppi_user01 + beta_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt_cate + aptt_cate, 
                          data = imputed_data_c, 
                          family = binomial("logit"))

# Propensity scoresを計算
ps_c <- predict(propensity_model_c, type = "response", newdata = imputed_data_c)

# ATT-IPTWを計算
att_weights_c <- ifelse(imputed_data_c$antibio_exposure == 1, 1/ps_c, 1/(1-ps_c))
iptw_model_c <- glm(outcome ~ antibio_exposure, data = imputed_data_c, weights = att_weights_c, family = binomial("logit"))

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# 結果を取得
estimate_c <- coef(iptw_model_c)[2]
std_error_c <- coef(summary(iptw_model_c))[2, "Std. Error"]
conf_int_c <- confint(iptw_model_c, level = 0.95, method = "Wald")

## Waiting for profiling to be done...

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# Odds ratioとconfidence intervalsを計算
or_c <- exp(estimate_c)
lower_c <- exp(estimate_c - 1.96 * std_error_c)
upper_c <- exp(estimate_c + 1.96 * std_error_c)

# p-valueを計算
p_c <- 2 * (1 - pnorm(abs(estimate_c / std_error_c)))

# 結果をdata frameに保存
result_c <- data.frame(
  Child_Score = "b",
  OR = or_c,
  Lower_CI = lower_c,
  Upper_CI = upper_c,
  P_value = p_c
)

print(result_c)

##                  Child_Score        OR  Lower_CI Upper_CI   P_value
## antibio_exposure           b 0.7932464 0.4564749 1.378476 0.4113442

interaction ab vs c

# 単一の完成データセットを取得
imputed_data <- complete(impglm, 1)

# `child_score` を「a, b」と「c」で2つのカテゴリに再カテゴリ化
imputed_data$child_grouped <- ifelse(imputed_data$child_score %in% c("a", "b"), "ab", "c")

# GLM propensity score modelを適合
propensity_model <- glm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hcc + alocohol + past_varix_rup + anticoag_user01 + nsaids_user01 + ppi_user01 + beta_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt_cate + aptt_cate, 
                        data = imputed_data, 
                        family = binomial("logit"))

# Propensity scoresを計算
ps <- predict(propensity_model, type = "response", newdata = imputed_data)
imputed_data$ps <- ps

# ATT-IPTWを計算
iptw <- with(imputed_data, ifelse(antibio_exposure == 1, 1/ps, 1/(1 - ps)))
imputed_data$iptw <- iptw

# IPTWを用いてlogistic regression modelを適合
iptw_model <- glm(outcome ~ antibio_exposure * child_grouped, data = imputed_data, weights = iptw, family = binomial("logit"))

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# 交互作用項のp値を取得
interaction_p <- coef(summary(iptw_model))["antibio_exposure:child_groupedc", "Pr(>|z|)"]

print(paste("Interaction P-value:", round(interaction_p, 4)))

## [1] "Interaction P-value: 0.326"

library(dplyr)

# Child-Pughスコアを因子として扱う
imputed_data$child_score <- as.factor(imputed_data$child_score)

# 交互作用項を含むpropensity score modelの適合
propensity_model_interaction <- glm(antibio_exposure ~ age_cate + sex + brthel_cate + child_numl + cci_num + malignancy + hcc + alocohol + past_varix_rup + anticoag_user01 + nsaids_user01 + ppi_user01 + beta_user01 + vaso + map_day0 + shock + t_bil + ast + alt + wbc + hb + plt + alb + eGFR30 + crp + pt_cate + aptt_cate + child_score, 
                                    data = imputed_data, 
                                    family = binomial("logit"))

# Propensity scoresを計算
ps_interaction <- predict(propensity_model_interaction, type = "response", newdata = imputed_data)

# ATT-IPTWを計算
att_weights_interaction <- ifelse(imputed_data$antibio_exposure == 1, 1/ps_interaction, 1/(1-ps_interaction))

# 交互作用項を持つモデルを適合
iptw_interaction_model <- glm(outcome ~ antibio_exposure * child_score, 
                              data = imputed_data, 
                              weights = att_weights_interaction, 
                              family = binomial("logit"))

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# 結果のサマリー
summary(iptw_interaction_model)

## 
## Call:
## glm(formula = outcome ~ antibio_exposure * child_score, family = binomial("logit"), 
##     data = imputed_data, weights = att_weights_interaction)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3141  -0.7453  -0.4362  -0.3881   6.7077  
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    -3.1749     0.4290  -7.401 1.35e-13 ***
## antibio_exposure               -0.2997     0.6684  -0.448 0.653842    
## child_scoreb                    0.4933     0.4720   1.045 0.295944    
## child_scorec                    1.7854     0.4613   3.871 0.000108 ***
## antibio_exposure:child_scoreb   0.4435     0.7206   0.615 0.538249    
## antibio_exposure:child_scorec   0.7035     0.7055   0.997 0.318680    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1096.52  on 789  degrees of freedom
## Residual deviance:  998.48  on 784  degrees of freedom
## AIC: 975.57
## 
## Number of Fisher Scoring iterations: 5

# 交互作用項のp値を取得
interaction_p_value <- coef(summary(iptw_interaction_model))["antibio_exposure:child_scorec", "Pr(>|z|)"]

print(paste("交互作用項のp値: ", interaction_p_value))

## [1] "交互作用項のp値:  0.318679895035494"